text:
@testset "orc" begin

@testset "Undefined Symbol" begin
    tm  = JITTargetMachine()
    OrcJIT(tm) do orc
       Context() do ctx
            mod = LLVM.Module("jit"; ctx)
            T_Int32 = LLVM.Int32Type(ctx)
            ft = LLVM.FunctionType(T_Int32, [T_Int32, T_Int32])
            fn = LLVM.Function(mod, "mysum", ft)
            linkage!(fn, LLVM.API.LLVMExternalLinkage)

            fname = mangle(orc, "wrapper")
            wrapper = LLVM.Function(mod, fname, ft)
            # generate IR
            Builder(ctx) do builder
                entry = BasicBlock(wrapper, "entry"; ctx)
                position!(builder, entry)

                tmp = call!(builder, fn, [parameters(wrapper)...])
                ret!(builder, tmp)
            end

            triple!(mod, triple(tm))
            ModulePassManager() do pm
                add_library_info!(pm, triple(mod))
                add_transform_info!(pm, tm)
                run!(pm, mod)
            end
            verify(mod)

            orc_mod = compile!(orc, mod)
            @test_throws ErrorException address(orc, fname) == C_NULL

            delete!(orc, orc_mod)
        end
    end
end

@testset "Custom Resolver" begin
    tm  = JITTargetMachine()
    OrcJIT(tm) do orc
       Context() do ctx
            known_functions = Dict{String, OrcTargetAddress}()
            fnames = Dict{String, Int}()
            function lookup(name, ctx)
                name = unsafe_string(name)
                try
                    if !haskey(fnames, name)
                        fnames[name] = 0
                    end
                    fnames[name] += 1

                    return known_functions[name].ptr
                catch ex
                    @error "Exception during lookup" name exception=(ex, catch_backtrace())
                    error("OrcJIT: Could not find symbol")
                end
            end

            mod = LLVM.Module("jit"; ctx)
            T_Int32 = LLVM.Int32Type(ctx)
            ft = LLVM.FunctionType(T_Int32, [T_Int32, T_Int32])
            fn = LLVM.Function(mod, "mysum", ft)
            linkage!(fn, LLVM.API.LLVMExternalLinkage)

            fname = mangle(orc, "wrapper")
            wrapper = LLVM.Function(mod, fname, ft)
            # generate IR
            Builder(ctx) do builder
                entry = BasicBlock(wrapper, "entry"; ctx)
                position!(builder, entry)

                tmp = call!(builder, fn, [parameters(wrapper)...])
                ret!(builder, tmp)
            end

            triple!(mod, triple(tm))
            ModulePassManager() do pm
                add_library_info!(pm, triple(mod))
                add_transform_info!(pm, tm)
                run!(pm, mod)
            end
            verify(mod)

            mysum_name = mangle(orc, "mysum")
            known_functions[mysum_name] = OrcTargetAddress(@cfunction(+, Int32, (Int32, Int32)))

            f_lookup = @cfunction($lookup, UInt64, (Cstring, Ptr{Cvoid}))
            GC.@preserve f_lookup begin
                orc_mod = compile!(orc, mod, f_lookup, C_NULL, lazy=true) # will capture f_lookup

                addr = address(orc, fname)
                @test errormsg(orc) == ""
                addr2 = addressin(orc, orc_mod, fname)

                @test addr == addr2
                @test addr.ptr != 0
                @test !haskey(fnames, mysum_name)

                r = ccall(pointer(addr), Int32, (Int32, Int32), 1, 2) # uses f_lookup
                @test r == 3
            end

            @test haskey(fnames, mysum_name)
            @test fnames[mysum_name] == 1

            empty!(fnames)
            delete!(orc, orc_mod)
        end
    end
end

@testset "Default Resolver + Stub" begin
    tm  = JITTargetMachine()
    OrcJIT(tm) do orc
            Context() do ctx
            mod = LLVM.Module("jit"; ctx)
            T_Int32 = LLVM.Int32Type(ctx)
            ft = LLVM.FunctionType(T_Int32, [T_Int32, T_Int32])
            fn = LLVM.Function(mod, "mysum", ft)
            linkage!(fn, LLVM.API.LLVMExternalLinkage)

            fname = mangle(orc, "wrapper")
            wrapper = LLVM.Function(mod, fname, ft)
            # generate IR
            Builder(ctx) do builder
                entry = BasicBlock(wrapper, "entry"; ctx)
                position!(builder, entry)

                tmp = call!(builder, fn, [parameters(wrapper)...])
                ret!(builder, tmp)
            end

            triple!(mod, triple(tm))
            ModulePassManager() do pm
                add_library_info!(pm, triple(mod))
                add_transform_info!(pm, tm)
                run!(pm, mod)
            end
            verify(mod)

            create_stub!(orc, mangle(orc, "mysum"), OrcTargetAddress(@cfunction(+, Int32, (Int32, Int32))))

            orc_mod = compile!(orc, mod)

            addr = address(orc, fname)
            @test errormsg(orc) == ""

            r = ccall(pointer(addr), Int32, (Int32, Int32), 1, 2)
            @test r == 3

            delete!(orc, orc_mod)
        end
    end
end

@testset "Default Resolver + Global Symbol" begin
    tm  = JITTargetMachine()
    OrcJIT(tm) do orc
        Context() do ctx
            mod = LLVM.Module("jit"; ctx)
            T_Int32 = LLVM.Int32Type(ctx)
            ft = LLVM.FunctionType(T_Int32, [T_Int32, T_Int32])
            mysum = mangle(orc, "mysum")
            fn = LLVM.Function(mod, mysum, ft)
            linkage!(fn, LLVM.API.LLVMExternalLinkage)

            fname = mangle(orc, "wrapper")
            wrapper = LLVM.Function(mod, fname, ft)
            # generate IR
            Builder(ctx) do builder
                entry = BasicBlock(wrapper, "entry"; ctx)
                position!(builder, entry)

                tmp = call!(builder, fn, [parameters(wrapper)...])
                ret!(builder, tmp)
            end

            triple!(mod, triple(tm))
            ModulePassManager() do pm
                add_library_info!(pm, triple(mod))
                add_transform_info!(pm, tm)
                run!(pm, mod)
            end
            verify(mod)

            # Should do pretty much the same as `@ccallable`
            LLVM.add_symbol(mysum, @cfunction(+, Int32, (Int32, Int32)))
            ptr = LLVM.find_symbol(mysum)
            @test ptr !== C_NULL
            @test ccall(ptr, Int32, (Int32, Int32), 1, 2) == 3

            orc_mod = compile!(orc, mod, lazy=true)

            addr = address(orc, fname)
            @test errormsg(orc) == ""

            r = ccall(pointer(addr), Int32, (Int32, Int32), 1, 2)
            @test r == 3

            delete!(orc, orc_mod)
        end
    end
end

@testset "Loading ObjectFile" begin
    tm  = JITTargetMachine()
    OrcJIT(tm) do orc
        Context() do ctx
            sym = mangle(orc, "SomeFunction")

            mod = LLVM.Module("jit"; ctx)
            ft = LLVM.FunctionType(LLVM.VoidType(ctx))
            fn = LLVM.Function(mod, sym, ft)

            Builder(ctx) do builder
                entry = BasicBlock(fn, "entry"; ctx)
                position!(builder, entry)
                ret!(builder)
            end
            verify(mod)

            obj = emit(tm, mod, LLVM.API.LLVMObjectFile)

            orc_m = add!(orc, MemoryBuffer(obj))
            addr = address(orc, sym)

            @test addr.ptr != 0
            delete!(orc, orc_m)
        end
    end
end

@testset "Stubs" begin
    tm  = JITTargetMachine()
    OrcJIT(tm) do orc
        Context() do ctx
            toggle = Ref{Bool}(false)
            on()  = (toggle[] = true; nothing)
            off() = (toggle[] = false; nothing)

            # Note that `CFunction` objects can be GC'd (???)
            # and we capture them below.
            func_on = @cfunction($on, Cvoid, ())
            GC.@preserve func_on begin
                ptr = Base.unsafe_convert(Ptr{Cvoid}, func_on)

                create_stub!(orc, "mystub", OrcTargetAddress(ptr))
                addr = address(orc, "mystub")

                @test addr.ptr != 0
                @test toggle[] == false

                ccall(pointer(addr), Cvoid, ())
                @test toggle[] == true
            end

            func_off = @cfunction($off, Cvoid, ())
            GC.@preserve func_off begin
                ptr = Base.unsafe_convert(Ptr{Cvoid}, func_off)

                set_stub!(orc, "mystub", OrcTargetAddress(ptr))

                @test addr == address(orc, "mystub")
                @test toggle[] == true

                ccall(pointer(addr), Cvoid, ())
                @test toggle[] == false
            end
        end
    end
end

@testset "callback!" begin
    tm  = JITTargetMachine()
    OrcJIT(tm) do orc
        triggered = Ref{Bool}(false)

        # Setup the lazy callback for creating a module
        function callback(orc_ref::LLVM.API.LLVMOrcJITStackRef, callback_ctx::Ptr{Cvoid})
            orc = OrcJIT(orc_ref)
            sym = mangle(orc, "SomeFunction")

            # 1. IRGen & Optimize the module
            orc_mod = Context() do ctx
                mod = LLVM.Module("jit"; ctx)
                ft = LLVM.FunctionType(LLVM.VoidType(ctx))
                fn = LLVM.Function(mod, sym, ft)

                Builder(ctx) do builder
                    entry = BasicBlock(fn, "entry"; ctx)
                    position!(builder, entry)
                    ret!(builder)
                end
                verify(mod)

                triple!(mod, triple(tm))
                ModulePassManager() do pm
                    add_library_info!(pm, triple(mod))
                    add_transform_info!(pm, tm)
                    run!(pm, mod)
                end
                verify(mod)

                # 2. Add the IR module to the JIT
                return compile!(orc, mod)
            end

            # 3. Obtain address of compiled module
            addr = addressin(orc, orc_mod, sym)

            # 4. Update the stub pointer to point to the recently compiled module
            set_stub!(orc, "lazystub", addr)

            # 5. Return the address of tie implementation, since we are going to call it now
            triggered[] = true
            return addr.ptr
        end
        c_callback = @cfunction($callback, UInt64, (LLVM.API.LLVMOrcJITStackRef, Ptr{Cvoid}))

        GC.@preserve c_callback begin
            initial_addr = callback!(orc, c_callback, C_NULL)
            create_stub!(orc, "lazystub", initial_addr)
            addr = address(orc, "lazystub")

            ccall(pointer(addr), Cvoid, ()) # Triggers compilation
        end
        @test triggered[]
    end
end

end


text:
#= This is The same implementation of schelling.jl, except that we run replicate simulations in parallel.

To run simulations in parallel, you need to define all types and functions on all the workers, i.e. processing cores. Use the `@everywhere` for that.
=#

using Distributed
addprocs(3)  # Add the number of parallel processing units
@everywhere using Agents

# Create agent, model, and grid types
@everywhere mutable struct SchellingAgent{T<:Integer} <: AbstractAgent # An agent
# object should always be a subtype of AbstractAgent
  id::T
  pos::Tuple{T, T}
  mood::Bool # true is happy and false is unhappy
  group::T
end

@everywhere mutable struct SchellingModel{T<:Integer, Y<:AbstractArray, Z<:AbstractSpace} <: AbstractModel  # A model
	# object should always be a subtype of AbstractModel
	space::Z  # A space object, which is a field
	# of the model object is always subtype of AbstractSpace
	agents::Y  # a list of agents
	scheduler::Function
	min_to_be_happy::T  # minimum number of neighbors
	#to be of the same kind so that they are happy
end

@everywhere mutable struct MyGrid{T<:Integer, Y<:AbstractArray} <: AbstractSpace
  dimensions::Tuple{T, T}
  space::SimpleGraph
  agent_positions::Y  # an array of arrays for each grid node
end

# instantiate the model
@everywhere function instantiate_model(;numagents=320, griddims=(20, 20), min_to_be_happy=3)
  agent_positions = [Int64[] for i in 1:gridsize(griddims)]
  mygrid = MyGrid(griddims, grid(griddims, false, true), agent_positions)  # create a 2D grid with nodes have a max of 8 neighbors
  model = SchellingModel(mygrid, SchellingAgent[], random_activation, min_to_be_happy)

  agents = vcat([SchellingAgent(Int(i), (1,1), false, 0) for i in 1:(numagents/2)], [SchellingAgent(Int(i), (1,1), false, 1) for i in (numagents/2)+1:numagents])
  for agent in agents
    add_agent_single!(agent, model)
  end
  return model
end

@everywhere function agent_step!(agent, model)
  if agent.mood == true
    return
  end
  while agent.mood == false
    neighbor_cells = node_neighbors(agent, model)
    same = 0
    for nn in neighbor_cells
      nsid = get_node_contents(nn, model)
      if length(nsid) == 0
        continue
      else
        nsid = nsid[1]
      end
      ns = model.agents[nsid].group
      if ns == agent.group
        same += 1
      end
    end
    if same >= model.min_to_be_happy
      agent.mood = true
    else
      move_agent_single!(agent, model)
    end
  end
end

@everywhere agent_properties = [:pos, :mood, :group]
@everywhere steps_to_collect_data = collect(1:2)
model = instantiate_model(numagents=370, griddims=(20,20), min_to_be_happy=3);
nsteps = 3
nreplicates = 10
all_data = batchrunner_parallel(agent_step!, model, nsteps, agent_properties, steps_to_collect_data, nreplicates);

ext: jl


text:
module RessourceManager

export PATHS

export SetWorkingDir
export AddPath
export GetPath
export RessourcePath
export CurrentDay

PATHS = Dict{Any,Any}(:ROOT => "")

SetWorkingDir() = (PATHS[:ROOT] = string(dirname(Base.source_path()),"/"))
RessourcePath(path::AbstractString) = (!in(':',path) ? path = string(PATHS[:ROOT], path, path[length(path)] != '/' ? "/" : "") : path)
AddPath(id::Any, path::AbstractString) = PATHS[id] = RessourcePath(path)
GetPath(key::Any) = PATHS[haskey(PATHS, key) ? key : :ROOT]
CurrentDay() = replace(replace(string(Dates.now()),":","-"),"T","_")

SetWorkingDir()
end # RessourceManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
#!/usr/bin/env julia
include("se_dseason.jl")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{fontspec}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\setmainfont{Asana-Math}")
pgfsave("../../images/se_dseason_81_19_kv.svg", plot_seasontypes("women", 1981, 2019))
pgfsave("../../images/se_dseason_81_19_m.svg", plot_seasontypes("men", 1981, 2019))

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
module MaxHelpingHandMultiRoomStrawberrySeed

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/MultiRoomStrawberrySeed" MultiRoomStrawberrySeed(x::Integer, y::Integer,
    strawberryName::String="multi_room_strawberry", sprite::String="MaxHelpingHand/miniberry/miniberry", ghostSprite::String="MaxHelpingHand/miniberry/ghostminiberry", index::Int=-1)

const bundledSprites = String["strawberry/seed", "MaxHelpingHand/miniberry/miniberry"]
const bundledGhostSprites = String["ghostberry/seed", "MaxHelpingHand/miniberry/ghostminiberry"]

const placements = Ahorn.PlacementDict(
    "Multi-Room Strawberry Seed (max480's Helping Hand)" => Ahorn.EntityPlacement(
        MultiRoomStrawberrySeed
    )
)

Ahorn.editingOptions(entity::MultiRoomStrawberrySeed) = Dict{String,Any}(
    "sprite" => bundledSprites,
    "ghostSprite" => bundledGhostSprites
)

function Ahorn.selection(entity::MultiRoomStrawberrySeed)
    x, y = Ahorn.position(entity)
    sprite = "collectables/" * get(entity.data, "sprite", "strawberry/seed") * "00"

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::MultiRoomStrawberrySeed, room::Maple.Room)
    sprite = "collectables/" * get(entity.data, "sprite", "strawberry/seed") * "00"

    Ahorn.drawSprite(ctx, sprite, 0, 0)
end

end

max_stars_repo_name: swoolcock/MaxHelpingHand

max_issues_repo_name: swoolcock/MaxHelpingHand

max_forks_repo_name: swoolcock/MaxHelpingHand


text:
# QuEST_jl/src/base/calculations.jl
#

function calcDensityInnerProduct(ϱ1 ::QuEST_Types.Qureg, ϱ2 ::QuEST_Types.Qureg) ::Qreal

    w = ccall(:calcDensityInnerProduct,
              Qreal,
              (QuEST_Types.Qureg, QuEST_Types.Qureg),
              ϱ1,
              ϱ2)

    return Qreal(w)
end

function calcExpecDiagonalOp(qureg      ::QuEST_Types.Qureg, op         ::QuEST_Types.DiagonalOp)  ::QuEST_Types.Complex

    ret= ccall(:calcExpecDiagonalOp,
               QuEST_Types.Complex,
               (QuEST_Types.Qureg, QuEST_Types.DiagonalOp),
               qureg,
               op)

    return ret
end

function calcExpecPauliHamil(qureg           ::QuEST_Types.Qureg,
                             hamil           ::QuEST_Types.PauliHamil,
                             workspace       ::QuEST_Types.Qureg)       Qreal

    ret = ccall(:calcExpecPauliHamil, Qreal, (QuEST_Types.Qureg, QuEST_Types.PauliHamil, QuEST_Types.Qureg), qureg, hamil, workspace)

    return ret

end

function calcExpecPauliProd(qureg        ::QuEST_Types.Qureg,
                            targetQubits ::Vector{QubitIdx},
                            pauliCodes   ::Vector{QuEST_Types.pauliOpType},
                            workspace    ::QuEST_Types.Qureg)          ::Qreal

    @assert length(targetQubits) == length(pauliCodes)
    #@assert all( σ -> 0 ≤ σ ≤ 3,   pauliCodes )

    expval = ccall(:calcExpecPauliProd,
                   Qreal,
                   (QuEST_Types.Qureg, Ptr{QubitIdx}, Ptr{QuEST_Types.pauliOpType}, Cint,                 QuEST_Types.Qureg),
                   qureg,              targetQubits,  pauliCodes,                   length(targetQubits), workspace)

    return expval
end

function calcExpecPauliSum(qureg         ::QuEST_Types.Qureg,
                           allPauliCodes ::Vector{QuEST_Types.pauliOpType},
                           termCoeffs    ::Vector{Qreal},
                           workspace     ::QuEST_Types.Qureg)          ::Float64

    @assert length(allPauliCodes) ==  length(termCoeffs) * getNumQubits(qureg)
    #@assert all( σ -> 0 ≤ σ ≤ 3,  allPauliCodes )

    ex = ccall(:calcExpecPauliSum,
               Qreal,
               (QuEST_Types.Qureg, Ptr{QuEST_Types.pauliOpType}, Ptr{Qreal}, Cint, QuEST_Types.Qureg),
               qureg,
               allPauliCodes,
               termCoeffs,
               Cint(length(termCoeffs)),
               workspace)

    return ex

end

function calcFidelity(qureg ::QuEST_Types.Qureg, pureState ::QuEST_Types.Qureg) ::Qreal

    fi = ccall(:calcFidelity, Qreal, (QuEST_Types.Qureg, QuEST_Types.Qureg), qureg, pureState)
    return fi

end

function calcHilbertSchmidtDistance(a ::QuEST_Types.Qureg, b ::QuEST_Types.Qureg) ::Qreal

    hsd = ccall(:calcHilbertSchmidtDistance, Qreal, (QuEST_Types.Qureg, QuEST_Types.Qureg), a, b)
    return hsd

end

function calcInnerProduct(bra ::QuEST_Types.Qureg, ket ::QuEST_Types.Qureg) ::Complex{Qreal}

    w = ccall(:calcInnerProduct, QuEST_Types.Complex, (QuEST_Types.Qureg, QuEST_Types.Qureg), bra, ket)
    return Complex{Qreal}(w.real,w.imag)

end

function calcProbOfOutcome(qureg        ::QuEST_Types.Qureg,
                           measureQubit ::Integer,
                           outcome      ::Integer)   ::Qreal

    p = ccall(:calcProbOfOutcome,
              Qreal,
              (QuEST_Types.Qureg, QubitIdx,      Cint),
              qureg,              measureQubit,  outcome)
    return p
end

function calcPurity(qureg ::QuEST_Types.Qureg) ::Qreal

    pu = ccall(:calcPurity, Qreal, (QuEST_Types.Qureg,), qureg)
    return pu

end

function calcTotalProb(qureg ::QuEST_Types.Qureg) ::Qreal

    one = ccall(:calcTotalProb, Qreal, (QuEST_Types.Qureg,), qureg)
    return one

end

function getAmp(qureg ::QuEST_Types.Qureg,  idx ::Integer) ::Complex{Qreal}

    α = ccall(:getAmp, QuEST_Types.Complex,
              (QuEST_Types.Qureg, Clonglong),
              qureg,              idx)
    return Complex{Qreal}(α.real,α.imag)

end

function getDensityAmp(qureg ::QuEST_Types.Qureg, row ::Integer, col ::Integer) ::Complex{Qreal}

    α = ccall(:getDensityAmp,
              QuEST_Types.Complex,
              (QuEST_Types.Qureg, Clonglong, Clonglong),
              qureg,
              Clonglong(row),
              Clonglong(col))

    return Complex{Qreal}(α.real, α.imag)

end

function getImagAmp(qureg      ::QuEST_Types.Qureg,
                    index      ::Integer)      ::Qreal

    ret = ccall(:getImagAmp, Qreal, (QuEST_Types.Qureg, Clonglong), qureg, Clonglong(index))
    return ret

end

function getNumAmps(qureg ::QuEST_Types.Qureg) ::Clonglong

    return ccall(:getNumAmps, Clonglong, (QuEST_Types.Qureg,), qureg)

end

function getNumQubits(qureg ::QuEST_Types.Qureg) ::QubitIdx

    return ccall(:getNumQubits, Cint, (QuEST_Types.Qureg,), qureg)

end

function getProbAmp(qureg ::QuEST_Types.Qureg,
                    idx   ::Integer) :: Qreal

    p = ccall(:getProbAmp, Qreal, (QuEST_Types.Qureg, Clonglong), qureg, Clonglong(idx))
    return p

end

function getRealAmp(qureg ::QuEST_Types.Qureg,
                    idx   ::Integer) :: Qreal

    p = ccall(:getRealAmp, Qreal, (QuEST_Types.Qureg, Clonglong), qureg, Clonglong(idx))
    return p

end

#EOF


text:
# winner from benchmarks/which_is_the_fastest_UInt_histogram.jl
"""
    uint_hist(bits, [RADIX_SIZE = 16, RADIX_MASK = 0xffff])

Computes a histogram (counts) for the vector RADIX_SIZE bits at a time. E.g. if eltype(bits) is UInt64 and RADIX_SIZE is 16
then 4 histograms are created for each of the 16 bit chunks.
"""
function uint_hist(bits::AbstractVector{T}, RADIX_SIZE = 16, RADIX_MASK = 0xffff) where T
    iter = ceil(Integer, sizeof(T)*8/RADIX_SIZE)
    hist = zeros(UInt32, 2^RADIX_SIZE, iter)

    for j = 1:length(bits)
        for i = 0:iter-1
            @inbounds hist[1+Int((uint_mapping(bits[j]) >> (i * RADIX_SIZE)) & RADIX_MASK), i+1] += 1
        end
    end
    hist
end


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

# using Markdown
# using InteractiveUtils

# ╔═╡ a1e38635-d1cb-4184-8e77-dab98a057e4c
#=╠═╡ notebook_exclusive
begin
	import Pkg
	Pkg.activate("..")
end
  ╠═╡ notebook_exclusive =#

# ╔═╡ 24ee155b-af1a-4915-a23b-3ff16e7ba844
begin
	using Revise
	using PlutoDevMacros
end

# ╔═╡ fa2b4748-4ec4-475a-a255-40b16cee11be
#=╠═╡ notebook_exclusive
using BenchmarkTools
  ╠═╡ notebook_exclusive =#

# ╔═╡ c5a1f2f0-c2b8-4ff7-b33c-467495b2317c
@plutoinclude "test2.jl" "all"

# ╔═╡ 7f9edbbe-2f8d-4b94-8199-6debf0143cca
#=╠═╡ notebook_exclusive
struct TestStruct3 end
  ╠═╡ notebook_exclusive =#

# ╔═╡ 25c3a29c-3a29-4de6-b66b-aaeb7c788d29
#=╠═╡ notebook_exclusive
asd(::TestStruct3) = "TESTSTRUCT3"
  ╠═╡ notebook_exclusive =#

# ╔═╡ 6630b18d-b98b-49b8-a416-e793a1d3b474
#=╠═╡ notebook_exclusive
asd(1)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 0542d7f4-4104-4884-b79e-d94d3fd3041c
#=╠═╡ notebook_exclusive
asd(1.0)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 25924828-2ec2-4df0-8afd-3bcb6e3aad28
#=╠═╡ notebook_exclusive
asd(TestStruct1())
  ╠═╡ notebook_exclusive =#

# ╔═╡ c0448179-1dac-4c61-9fb3-5b2c14f44f30
#=╠═╡ notebook_exclusive
asd(TestStruct2())
  ╠═╡ notebook_exclusive =#

# ╔═╡ 0d7ab7ce-1ce5-4249-b8de-8aeed5bec07e
#=╠═╡ notebook_exclusive
asd(TestStruct3())
  ╠═╡ notebook_exclusive =#

# ╔═╡ 720caf23-72de-4161-a7b0-0a10d690f0c3
#=╠═╡ notebook_exclusive
asd(TestStruct2())
  ╠═╡ notebook_exclusive =#

# ╔═╡ 26b39429-5766-45a5-a173-ad632c31718c
#=╠═╡ notebook_exclusive
@benchmark asd(TestStruct1())
  ╠═╡ notebook_exclusive =#

# ╔═╡ Cell order:
# ╠═a1e38635-d1cb-4184-8e77-dab98a057e4c
# ╠═24ee155b-af1a-4915-a23b-3ff16e7ba844
# ╠═fa2b4748-4ec4-475a-a255-40b16cee11be
# ╠═c5a1f2f0-c2b8-4ff7-b33c-467495b2317c
# ╠═6630b18d-b98b-49b8-a416-e793a1d3b474
# ╠═0542d7f4-4104-4884-b79e-d94d3fd3041c
# ╠═25924828-2ec2-4df0-8afd-3bcb6e3aad28
# ╠═c0448179-1dac-4c61-9fb3-5b2c14f44f30
# ╠═7f9edbbe-2f8d-4b94-8199-6debf0143cca
# ╠═25c3a29c-3a29-4de6-b66b-aaeb7c788d29
# ╠═0d7ab7ce-1ce5-4249-b8de-8aeed5bec07e
# ╠═720caf23-72de-4161-a7b0-0a10d690f0c3
# ╠═26b39429-5766-45a5-a173-ad632c31718c


text:
"""
    mpi_one_sided(data, comm = mpi_comm(), root = mpi_root; capacity)

Declare `data` as mpi-distributed and set communication strategy to one-sided with remote
memory access (RMA). `capacity` sets the capacity of the RMA windows.

Sets up the [`MPIData`](@ref) structure with [`MPIOneSided`](@ref) strategy.
"""
function mpi_one_sided(data, comm = mpi_comm(), root = mpi_root; capacity)
    MPI.Initialized() || error("MPI needs to be initialised first.")
    np = MPI.Comm_size(comm)
    id = MPI.Comm_rank(comm)
    P = eltype(storage(data))

    # compute the required capacity for the communication buffer as a
    # fraction of the capacity of `data`
    # id == root && println("on rank $id, capacity = ",cap)
    s = MPIOneSided(np, id, comm, P, Int32(capacity))
    return MPIData(data, comm, root, s)
end

"""
    MPIOneSided(nprocs, myrank, comm, ::Type{T}, capacity)

Communication buffer for use with MPI one-sided communication (remote memory access). Up to
`capacity` elements of type `T` can be exchanged between MPI ranks via [`put`](@ref). It is
important that `isbitstype(T) == true`. Objects of type `MPIOneSided` have to be
freed manually with a (blocking) call to [`free()`](@ref).
"""
mutable struct MPIOneSided{T}  <: DistributeStrategy
    mpiid::Int # unique ID for MPI-distributed objects
    np::Int32 # number of MPI processes
    id::Int32 # number of this rank
    comm::MPI.Comm # MPI communicator
    DT_b::MPI.Datatype # datatype for data
    b_win::MPI.Win # window for data buffer
    DT_l::MPI.Datatype # datatype for lengths of buffered data
    l_win::MPI.Win # window for length of buffered data
    capacity::Int32 # capacity of buffer
    b_vec::Vector{T} # local array for MPI window b_win
    l_vec::Vector{Int32} # local array for MPI window l_win

    function MPIOneSided(nprocs, myrank, comm, ::Type{T}, capacity) where T
        @assert isbitstype(T) "Buffer type for MPIOneSided needs to be isbitstype; found $T"
        b_vec = Vector{T}(undef, capacity)
        l_vec = zeros(Int32,1)
        DT_b = MPI.Datatype(T)
        b_win = MPI.Win_create(b_vec, comm) # separate windows for buffer and length
        DT_l = MPI.Datatype(Int32)
        l_win = MPI.Win_create(l_vec, comm)
        mpiid = next_mpiID()
        obj = new{T}(mpiid,nprocs,myrank,comm, DT_b,b_win, DT_l,l_win,capacity,b_vec,l_vec)
        mpi_registry[mpiid] = Ref(obj) # register the object to avoid
        # arbitrary garbage collection
        # ccall(:jl_, Cvoid, (Any,), "installing finalizer on MPIOneSided")
        finalizer(myclose, obj) # install finalizer
        return obj
    end
end
Base.eltype(::MPIOneSided{T}) where T = T

function Base.show(io::IO, s::MPIOneSided{T}) where T
    println(io, "MPIOneSided{$T}")
    println(io, "  mpiid: ",s.mpiid)
    println(io, "  np: ",s.np)
    println(io, "  id: ",s.id)
    println(io, "  comm: ",s.comm)
    println(io, "  MPI.Win and b_vec::Vector{$T} with capacity ",s.capacity)
    println(io, "  MPI.Win for number of elements in buffer")
end

function myclose(obj::MPIOneSided)
    MPI.free(obj.b_win)
    MPI.free(obj.l_win)
    return nothing
end

"""
    free(obj::MPIOneSided)
De-reference the object, call finalizer and the garbage collector immediately.
This is a syncronizing MPI call. Make sure that the object is not used later.
MPI syncronizing.
"""
function free(obj::MPIOneSided)
    global mpi_registry
    MPI.Barrier(obj.comm)
    delete!(mpi_registry, obj.mpiid)
    finalize(obj)
    GC.gc()
end
function free(d::MPIData{D,S}) where {D,S<:MPIOneSided}
    free(d.s)
    finalize(d)
    GC.gc()
end
function free(d::MPIData) # default
    finalize(d)
    GC.gc()
end

function fence(s::MPIOneSided, assert = 0)
    MPI.Win_fence(assert, s.b_win)
    MPI.Win_fence(assert, s.l_win) # joint fences
end

"""
    put(buf::Vector{T}, [len,] targetrank, s::MPIOneSided{T})

Deposit a vector `buf` into the MPI window `s` on rank `targetrank`. If
`len` is given, only the first `len` elements are transmitted.
"""
@inline function put(buf::Vector{T}, len, targetrank, s::MPIOneSided{T}) where T
    @boundscheck len ≤ length(buf) && len ≤ s.capacity ||
        error("Not enough space left in buffer")
    b_buffer = MPI.Buffer(buf, len, s.DT_b)
    MPI.Put(b_buffer, targetrank, s.b_win)
    l_buffer = MPI.Buffer([len,], 1, s.DT_l)
    MPI.Put(l_buffer, targetrank, s.l_win)
end
@inline function put(buf::Vector{T}, targetrank, s::MPIOneSided{T}) where T
    len = length(buf)
    put(buf, len, targetrank, s)
end

function sbuffer!(s::MPIOneSided) # unsafe version - avoids allocations
    # println("from sbuffer $(s.b_vec), $(s.l_vec), $(s.l_vec[1])")
    # return s.b_vec[1:s.l_vec[1]]
    return view(s.b_vec, 1 : s.l_vec[1])
end
function sbuffer(s::MPIOneSided) # safe version for reading shared buffer - returns local array
    MPI.Win_lock(MPI.LOCK_SHARED, s.id, 0, s.b_win)
    MPI.Win_lock(MPI.LOCK_SHARED, s.id, 0, s.l_win)
    res = s.b_vec[1 : s.l_vec[1]]
    MPI.Win_unlock(s.id,s.l_win)
    MPI.Win_unlock(s.id,s.b_win)
    return res
end

function mpi_combine_walkers!(target, source, s::MPIOneSided{P}) where P
    # now target is just a local data structure, e.g. DVec
    # allocate local buffers for sorting
    bufs = [Vector{P}(undef,length(source)) for i in 1:(s.np-1)] # type-stable
    lens = zeros(Int,(s.np-1))
    # sort source into send buffers
    @inbounds for (key,val) in pairs(source)
        tr = targetrank(key, s.np)
        if tr < s.id
            lens[tr+1] +=1
            bufs[tr+1][lens[tr+1]] = key => val
        elseif tr > s.id
            lens[tr] +=1
            bufs[tr][lens[tr]] = key => val
        else # tr == s.id
            deposit!(target, key, val, nothing)
        end
    end
    # send data with RMA communications to higher rank by `offset`
    for offset in 1 : s.np-1
        # println("$(s.id) before first fence")
        fence(s)
        # println("jumped first fence on $(s.id)")

        tr = (s.id + offset) % s.np # compute target rank for sending
        bufindex = tr > s.id ? tr : tr + 1 # and buffer index
        # println("Sending $(lens[bufindex]) pairs from $(s.id) to $tr")

        put(bufs[bufindex], lens[bufindex], tr, s)
        # println("sending done on $(s.id)")

        fence(s)
        # println("jumped the fence on $(s.id)")
        for (key, val) in sbuffer(s)
            deposit!(target, key, val, nothing)
        end
    end
    fence(s)
    return target
end


text:
using ParquetFiles
using Parquet
using IteratorInterfaceExtensions
using TableTraits
using Test

@testset "ParquetFiles" begin

    parquet_pkg_dir = joinpath(dirname(pathof(Parquet)), "..")

    include(joinpath(parquet_pkg_dir, "test", "get_parcompat.jl"))

    test_filename = joinpath(parquet_pkg_dir, "test", "parquet-compatibility", "parquet-testdata", "impala", "1.1.1-NONE", "nation.impala.parquet")

    pqf = load(test_filename)

    @test IteratorInterfaceExtensions.isiterable(pqf) == true
    @test TableTraits.isiterabletable(pqf) == true

    it = IteratorInterfaceExtensions.getiterator(pqf)

    ar = collect(it)

    @test length(ar) == 25
    @test ar[1] == (n_nationkey = 0, n_name = "ALGERIA", n_regionkey = 0, n_comment = " haggle. carefully final deposits detect slyly agai")

    @test sprint((stream, data) -> show(stream, "text/html", data), pqf)[1:100] == "

n_nationkey	n_name	n_regionkey	n_comment
 show(stream, "application/vnd.dataresource+json", data), pqf)[1:100] == "{\"schema\":{\"fields\":[{\"name\":\"n_nationkey\",\"type\":\"integer\"},{\"name\":\"n_name\",\"type\":\"string\"},{\"nam"
    @test sprint(show, pqf)[1:100] == "25x4 Parquet file\nn_nationkey │ n_name    │ n_regionkey\n────────────┼─"
    @test showable("text/html", pqf) == true
    @test showable("application/vnd.dataresource+json", pqf) == true

end


text:
@doc raw"""
    galerkin_projection(
        f, B::AbstractBSplineBasis,
        [deriv = Derivative(0)], [VectorType = Vector{Float64}],
    )

Perform Galerkin projection of a function `f` onto the given basis.

By default, returns a vector with values

```math
φ_i = ⟨ b_i, f ⟩
= ∫_a^b b_i(x) \, f(x) \, \mathrm{d}x,
```

where ``a`` and ``b`` are the boundaries of the B-spline basis
``\{ b_i \}_{i = 1}^N``.

The integrations are performed using Gauss--Legendre quadrature.
The number of quadrature nodes is chosen so that the result is exact when ``f``
is a polynomial of degree ``k - 1`` (or, more generally, a spline belonging to
the space spanned by the basis `B`).
Here ``k`` is the order of the B-spline basis.
In the more general case, this function returns a quadrature approximation of
the projection.

See also [`galerkin_projection!`](@ref) for the in-place operation, and
[`galerkin_matrix`](@ref) for more details.
"""
function galerkin_projection(
        f, B::AbstractBSplineBasis,
        deriv = Derivative(0),
        ::Type{V} = Vector{Float64},
    ) where {V <: AbstractVector}
    N = length(B)
    φ = V(undef, N)
    galerkin_projection!(f, φ, B, deriv)
end

"""
    galerkin_projection!(
        f, φ::AbstractVector, B::AbstractBSplineBasis, [deriv = Derivative(0)],
    )

Compute Galerkin projection ``φ_i = ⟨ b_i, f ⟩``.

See [`galerkin_projection`](@ref) for details.
"""
function galerkin_projection!(
        f, φ::AbstractVector, B::AbstractBSplineBasis, deriv = Derivative(0),
    )
    N = length(B)
    if length(φ) != N
        throw(DimensionMismatch("incorrect length of output vector φ"))
    end
    Base.require_one_based_indexing(φ)

    k = order(B)
    ts = knots(B)

    # Quadrature information (nodes, weights).
    quadx, quadw = _quadrature_prod(Val(2k - 2))
    @assert length(quadx) == k  # we need k quadrature points per knot segment

    fill!(φ, 0)

    nlast = last(eachindex(ts))

    # We loop over all knot segments Ω[n] = (ts[n], ts[n + 1]).
    # For all B-splines with support in this segment, we integrate the product
    # B[i] * f over this segment, adding the result to φ[i].
    @inbounds for n in eachindex(ts)
        n == nlast && break
        tn, tn1 = ts[n], ts[n + 1]
        tn1 == tn && continue  # interval of length = 0

        metric = QuadratureMetric(tn, tn1)

        # Unnormalise quadrature nodes, such that xs ∈ [tn, tn1]
        xs = metric .* quadx
        # @assert all(x -> tn ≤ x ≤ tn1, xs)

        is = nonzero_in_segment(B, n)

        # This is a property of B-spline bases, which should be preserved by
        # derived (recombined) bases.
        @assert 0 < length(is) ≤ k

        # Evaluate all required basis functions on quadrature nodes.
        bis = eval_basis_functions(B, is, xs, deriv)

        fs = ntuple(i -> f(xs[i]), Val(length(xs)))  # this works fine!

        # For some reason, the alternatives below allocate and give bad
        # performance... (on Julia 1.7-beta2)
        # fs = map(f, xs)
        # fs = f.(xs)

        for (ni, i) in enumerate(is)
            bs = bis[ni]
            φ[i] += metric.α * ((bs .* fs) ⋅ quadw)
        end
    end

    φ
end


text:
import PredictMD

if PredictMD.is_travis_ci() && !haskey(ENV, "PREDICTMD_TEST_GROUP")
    ENV["PREDICTMD_TEST_GROUP"] = lowercase(
        strip(
            get(
                ENV,
                "GROUP",
                "",
                )
            )
        )
end

const _test_group_environment_variable = lowercase(
    strip(
        get(
            ENV,
            "PREDICTMD_TEST_GROUP",
            ""
            )
        )
    )

if length(_test_group_environment_variable) == 0
    const _test_group_value = "default"
else
    const _test_group_value = _test_group_environment_variable
end

if haskey(TEST_GROUP_STRING_TO_INSTANCE, _test_group_value)
    const TEST_GROUP = TEST_GROUP_STRING_TO_INSTANCE[_test_group_value]
    @debug(
        string(
            "PREDICTMD_TEST_GROUP: \"",
            TEST_GROUP_INSTANCE_TO_STRING[TEST_GROUP],
            "\"",
            )
        )
else
    const _valid_test_group_values = string(
        "\"",
        join(
            sort(collect(keys(TEST_GROUP_STRING_TO_INSTANCE))),
            "\", \"",
            ),
        "\"",
        )
    error(
        string(
            "\"",
            _test_group_value,
            "\" is not a valid value for PREDICTMD_TEST_GROUP. ",
            "Valid values are: ",
            _valid_test_group_values,
            ".",
            )
        )
end


text:
function clean(config::Config)
    function rm_deadlinks(dir)
        for f in readdir(dir)
            file_path = abspath(joinpath(dir, f))
            src = file_path

            while islink(src)
                src = joinpath(dirname(src), readlink(src))
            end

            src = abspath(src)
            if !ispath(src)
                rm(file_path)
            end
        end
    end
    rm_deadlinks(config.dir.store)
    rm_deadlinks(config.dir.bin)
end

function Base.rm(config::Config; name::AbstractString="", dir::AbstractString="")
    if name != "" && dir == ""
        # If we find the name in the bin folder then we should just delete the julia symlink
        if name in readdir(abspath(config.dir.bin)) && name != "playground"
            rm(joinpath(abspath(config.dir.bin), name))
            return true
        # Otherwise the name should be in
        elseif name in readdir(abspath(config.dir.store))
            dir = get_playground_dir(config, "", name)

            # The dir returned could be a link
            # so we attempt to read that link.
            if islink(dir)
                try
                    dir = readlink(dir)
                catch
                    # If it fails just assume that we have a dead link
                    # and run clean_link and return
                    rm(dir)
                    return true
                end
            end
        else
            error("Unknown name $name")
        end
    elseif dir == "" && name == ""
        error("No julia-version, playground name or directory provided.")
    end

    # By this point dir should be valid or the function should have already exited.
    # DeclarativePackages creates a read-only directory so in case we run into that
    # during deletion we recursively chmod the path with write permissions.
    run(`chmod -R +w $(abspath(dir))`)
    warn("Recusively deleting $(abspath(dir))...")
    rm(abspath(dir), recursive=true)

    # Just to be safe run clean_links
    clean(config)
    return true
end


text:
"""
    DQMC(m::M; kwargs...) where M<:Model

Create a determinant quantum Monte Carlo simulation for model `m` with
keyword parameters `kwargs`.

### Keyword Arguments:
- `seed`: The random seed used by the simulation.
- `checkerboard=false`: If true, the simulation uses a generic checkerboard
decomposition.
- `thermalization_measurements::Dict{Symbol, AbstractMeasurement}`: A collection
of measurements run during the thermalization stage. By default, none are used.
- `measurements::Dict{Symbol, AbstractMeasurement}`: A collection of measurements
run during the measurement stage. Calls `default_measurements` if not specified.
- `recorder = ConfigRecorder`: Type of recorder used for saving configurations
generated during the simulation. Used (by default) when `replay!`ing simulations.
 (`Discarder` or `ConfigRecorder`)
- `recording_rate = measure_rate`: Rate at which configurations are recorded.
- `thermalization = 100`: Number of thermalization sweeps
- `sweeps`: Number of measurement sweeps
- `all_checks = true`: Check for Propagation instabilities and sign problems.
- `safe_mult = 10`: Number of "safe" matrix multiplications. Every `safe_mult`
multiplications, a UDT decomposition is used to stabilize the product.
- `delta_tau = 0.1`: Time discretization of the path integral
- `beta::Float64`: Inverse temperature used in the simulation
- `slices::Int = beta / delta_tau`: Number of imaginary time slice in the
simulation
- `measure_rate = 10`: Number of sweeps discarded between every measurement.
- `last_sweep = 0`: Sets the index of the last finished sweep. The simulation
will start with sweep `last_sweep + 1`.
"""
function DQMC(model::M;
        seed::Int=-1,
        checkerboard::Bool=false,
        thermalization_measurements = Dict{Symbol, AbstractMeasurement}(),
        measurements = Dict{Symbol, AbstractMeasurement}(),
        last_sweep = 0,
        measure_rate = 10,
        recording_rate = measure_rate,
        recorder = ConfigRecorder(DQMC, M, recording_rate),
        scheduler = SimpleScheduler(LocalSweep()),
        kwargs...
    ) where M<:Model
    # default params
    # paramskwargs = filter(kw->kw[1] in fieldnames(DQMCParameters), kwargs)
    parameters = DQMCParameters(measure_rate = measure_rate; kwargs...)

    HET = hoppingeltype(DQMC, model)
    GET = greenseltype(DQMC, model)
    HMT = hopping_matrix_type(DQMC, model)
    GMT = greens_matrix_type(DQMC, model)
    IMT = interaction_matrix_type(DQMC, model)
    stack = DQMCStack{GET, HET, GMT, HMT, IMT}()
    ut_stack = UnequalTimeStack{GET, GMT}()

    seed == -1 || Random.seed!(seed)
    conf = rand(DQMC, model, parameters.slices)
    analysis = DQMCAnalysis()
    CB = checkerboard ? CheckerboardTrue : CheckerboardFalse

    mc = DQMC(
        CB,
        model, conf, deepcopy(conf), last_sweep,
        stack, ut_stack, scheduler,
        parameters, analysis,
        recorder, thermalization_measurements, measurements
    )

    init!(mc)
    return mc
end

"""
    DQMC(m::M, params::Dict)
    DQMC(m::M, params::NamedTuple)

Create a determinant quantum Monte Carlo simulation for model `m` with
(keyword) parameters as specified in the dictionary/named tuple `params`.
"""
DQMC(m::Model, params::Dict{Symbol}) = DQMC(m; params...)
DQMC(m::Model, params::NamedTuple) = DQMC(m; params...)

# convenience
@inline beta(mc::DQMC) = mc.parameters.beta
@inline nslices(mc::DQMC) = mc.parameters.slices
@inline model(mc::DQMC) = mc.model
@inline conf(mc::DQMC) = mc.conf
@inline current_slice(mc::DQMC) = mc.stack.current_slice
@inline last_sweep(mc::DQMC) = mc.last_sweep
@inline configurations(mc::DQMC) = mc.recorder
@inline lattice(mc::DQMC) = lattice(mc.model)
@inline parameters(mc::DQMC) = merge(parameters(mc.parameters), parameters(mc.model))
@inline parameters(p::DQMCParameters) = (
    beta = p.beta, delta_tau = p.delta_tau, thermalization = p.thermalization, sweeps = p.sweeps
)

# cosmetics
import Base.summary
import Base.show
Base.summary(mc::DQMC) = "DQMC simulation of $(summary(mc.model))"
function Base.show(io::IO, mc::DQMC)
    print(io, "Determinant quantum Monte Carlo simulation\n")
    print(io, "Model: ", mc.model, "\n")
    print(io, "Beta: ", mc.parameters.beta, " (T ≈ $(round(1/mc.parameters.beta, sigdigits=3)))\n")
    N_th_meas = length(mc.thermalization_measurements)
    N_me_meas = length(mc.measurements)
    print(io, "Measurements: ", N_th_meas + N_me_meas, " ($N_th_meas + $N_me_meas)")
end
Base.show(io::IO, m::MIME"text/plain", mc::DQMC) = print(io, mc)

function init!(mc::DQMC)
    init_hopping_matrices(mc, mc.model)
    initialize_stack(mc, mc.stack)
    init_scheduler!(mc, mc.scheduler)
    # generate_communication_functions(mc.conf)
    nothing
end
@deprecate resume_init!(mc::DQMC) init!(mc) false

"""
    run!(mc::DQMC[; kwargs...])

Runs the given Monte Carlo simulation `mc`. Returns true if the run finished and
false if it cancelled early to generate a resumable save-file.

### Keyword Arguments:
- `verbose = true`: If true, print progress messaged to stdout.
- `thermalization`: Number of thermalization sweeps. Uses the value passed to
`DQMC` by default.
- `sweeps`: Number of measurement sweeps. Uses the value passed to `DQMC` by
default.
- `safe_every::TimePeriod`: Set the interval for regularly scheduled saves.
- `safe_before::Date`: If this date is passed, `run!` will generate a resumable
save file and exit
- `grace_period = Minute(5)`: Buffer between the current time and `safe_before`.
The time required to generate a save file should be included here.
- `resumable_filename`: Name of the resumable save file. The default is based on
`safe_before`.
- `overwrite = false`: If set to true a file with the same name as
`resumable_filename` will be overwritten. (This will create a temporary backup)

See also: [`resume!`](@ref)
"""
@bm function run!(
        mc::DQMC;
        verbose::Bool = true,
        ignore = tuple(),
        sweeps::Int = mc.parameters.sweeps,
        thermalization = mc.parameters.thermalization,
        safe_before::TimeType = now() + Year(100),
        safe_every::TimePeriod = Hour(10000),
        grace_period::TimePeriod = Minute(5),
        resumable_filename::String = "resumable_$(Dates.format(safe_before, "d_u_yyyy-HH_MM")).jld2",
        overwrite = false
    )

    # Update number of sweeps
    if (mc.parameters.thermalization != thermalization) || (mc.parameters.sweeps != sweeps)
        verbose && println("Rebuilding DQMCParameters with new number of sweeps.")
        mc.parameters = DQMCParameters(mc.parameters, thermalization = thermalization, sweeps = sweeps)
    end
    total_sweeps = sweeps + thermalization

    # Generate measurement groups
    th_groups = generate_groups(
        mc, mc.model,
        [mc.measurements[k] for k in keys(mc.thermalization_measurements) if !(k in ignore)]
    )
    groups = generate_groups(
        mc, mc.model,
        [mc.measurements[k] for k in keys(mc.measurements) if !(k in ignore)]
    )

    start_time = now()
    last_checkpoint = now()
    max_sweep_duration = 0.0
    verbose && println("Started: ", Dates.format(start_time, "d.u yyyy HH:MM"))

    # fresh stack
    verbose && println("Preparing Green's function stack")

    init!(mc)
    reverse_build_stack(mc, mc.stack)
    propagate(mc)

    _time = time() # for step estimations
    t0 = time() # for analysis.runtime, may need to reset
    verbose && println("\n\nThermalization stage - ", thermalization)

    next_print = (div(mc.last_sweep, mc.parameters.print_rate) + 1) * mc.parameters.print_rate
    while mc.last_sweep < total_sweeps
        verbose && (mc.last_sweep == thermalization + 1) && println("\n\nMeasurement stage - ", sweeps)

        # Perform whatever update is scheduled next
        update(mc.scheduler, mc, mc.model)

        # Trigger measurements
        if mc.last_sweep ≤ thermalization
            if iszero(mc.last_sweep % mc.parameters.measure_rate)
                for (requirement, group) in th_groups
                    apply!(requirement, group, mc, mc.model, mc.last_sweep)
                end
            end
            if mc.last_sweep == thermalization
                mc.analysis.th_runtime += time() - t0
                t0 = time()
            end
        else
            push!(mc.recorder, mc, mc.model, mc.last_sweep)
            if iszero(mc.last_sweep % mc.parameters.measure_rate)
                for (requirement, group) in groups
                    apply!(requirement, group, mc, mc.model, mc.last_sweep)
                end
            end
        end

        # Show sweep statistics - i.e. time/sweep, acceptance rates
        if mc.last_sweep >= next_print
            next_print += mc.parameters.print_rate
            sweep_dur = (time() - _time)/mc.parameters.print_rate
            max_sweep_duration = max(max_sweep_duration, sweep_dur)
            if verbose
                println("\t", mc.last_sweep)
                @printf("\t\tsweep dur: %.3fs\n", sweep_dur)
                show_statistics(stdout, mc.scheduler, "\t\t")
            end

            flush(stdout)
            _time = time()
        end

        # Trigger checkpoint safe/early exit save
        if safe_before - now() < Millisecond(grace_period) +
                Millisecond(round(Int, 2e3max_sweep_duration))

            println("Early save initiated for sweep #$(mc.last_sweep).\n")
            verbose && println("Current time: ", Dates.format(now(), "d.u yyyy HH:MM"))
            verbose && println("Target time:  ", Dates.format(safe_before, "d.u yyyy HH:MM"))
            save(resumable_filename, mc, overwrite = overwrite, rename = false)
            verbose && println("\nEarly save finished")
            disconnect(connected_ids)
            return false
        elseif (now() - last_checkpoint) > safe_every
            verbose && println("Performing scheduled save.")
            last_checkpoint = now()
            save(resumable_filename, mc, overwrite = overwrite, rename = false)
        end
    end

    if mc.last_sweep ≤ thermalization
        mc.analysis.th_runtime += time() - t0
    else
        mc.analysis.me_runtime += time() - t0
    end

    disconnect(connected_ids)

    # Print (numerical) error information
    if verbose
        if length(mc.analysis.imaginary_probability) > 0
            s = mc.analysis.imaginary_probability
            println("\nImaginary Probability Errors: ($(s.count))")
            @printf("\tmax  = %0.3e\n", max(s))
            @printf("\tmean = %0.3e\n", mean(s))
            @printf("\tmin  = %0.3e\n\n", min(s))
        end
        if length(mc.analysis.negative_probability) > 0
            s = mc.analysis.negative_probability
            println("\nNegative Probability Errors: ($(s.count))")
            @printf("\tmax  = %0.3e\n", max(s))
            @printf("\tmean = %0.3e\n", mean(s))
            @printf("\tmin  = %0.3e\n\n", min(s))
        end
        if length(mc.analysis.propagation_error) > 0
            s = mc.analysis.propagation_error
            println("\nPropagation Errors: ($(s.count))")
            @printf("\tmax  = %0.3e\n", max(s))
            @printf("\tmean = %0.3e\n", mean(s))
            @printf("\tmin  = %0.3e\n\n", min(s))
        end
    end

    # Total timings
    end_time = now()
    if verbose
        println("\nEnded: ", Dates.format(end_time, "d.u yyyy HH:MM"))
        @printf("Duration: %.2f minutes", (end_time - start_time).value/1000. /60.)
        println()
    end

    return true
end

"""
    replay(mc::DQMC[; configurations::Iterable = mc.recorder; kwargs...])

Replays previously generated configurations and measures observables along the
way.

### Keyword Arguments (both):
- `verbose = true`: If true, print progress messaged to stdout.
- `safe_every::TimePeriod`: Set the interval for regularly scheduled saves.
- `safe_before::Date`: If this date is passed, `replay!` will generate a
resumable save file and exit
- `grace_period = Minute(5)`: Buffer between the current time and `safe_before`.
The time required to generate a save file should be included here.
- `filename`: Name of the save file. The default is based on `safe_before`.
- `start=1`: The first sweep in the simulation. This will be changed when using
`resume!(save_file)`.
- `ignore`: A collection of measurement keys to ignore. Defaults to the key of
the configuration measurement.
- `measure_rate = 1`: Rate at which measurements are taken. Note that this is
based on the recorded configurations, not actual sweeps.
"""
function replay!(
        mc::DQMC, configurations = mc.recorder;
        ignore = tuple(),
        verbose::Bool = true,
        safe_before::TimeType = now() + Year(100),
        safe_every::TimePeriod = Hour(10000),
        grace_period::TimePeriod = Minute(5),
        resumable_filename::String = "resumable_$(Dates.format(safe_before, "d_u_yyyy-HH_MM")).jld",
        overwrite = false,
        measure_rate = 1
    )
    start_time = now()
    last_checkpoint = now()
    max_sweep_duration = 0.0
    verbose && println("Started: ", Dates.format(start_time, "d.u yyyy HH:MM"))

    # Check for measurements
    isempty(mc.thermalization_measurements) && @debug(
        "There is no thermalization process in a replayed simulation."
    )
    isempty(mc.measurements) && @warn(
        "There are no measurements set up for this simulation!"
    )
    # Generate measurement groups
    groups = generate_groups(
        mc, mc.model,
        [mc.measurements[k] for k in keys(mc.measurements) if !(k in ignore)]
    )

    if measure_rate != mc.parameters.measure_rate
        mc.parameters = DQMCParameters(mc.parameters, measure_rate = measure_rate)
    end

    verbose && println("Preparing Green's function stack")
    init!(mc)
    build_stack(mc, mc.stack)
    propagate(mc)
    mc.stack.current_slice = 1
    mc.conf = rand(DQMC, mc.model, nslices(mc))

    _time = time()
    verbose && println("\n\nReplaying measurement stage - ", length(configurations))
    prepare!(mc.measurements, mc, mc.model)
    for i in mc.last_sweep+1:mc.parameters.measure_rate:length(configurations)
        copyto!(mc.conf, decompress(mc, mc.model, configurations[i]))
        calculate_greens(mc, 0) # outputs to mc.stack.greens
        for (requirement, group) in groups
            apply!(requirement, group, mc, mc.model, i)
        end
        mc.last_sweep = i

        if mod(i, mc.parameters.print_rate) == 0
            sweep_dur = (time() - _time)/mc.parameters.print_rate
            max_sweep_duration = max(max_sweep_duration, sweep_dur)
            if verbose
                println("\t", i)
                @printf("\t\tsweep dur: %.3fs\n", sweep_dur)
            end
            flush(stdout)
            _time = time()
        end

        if safe_before - now() < Millisecond(grace_period) +
                Millisecond(round(Int, 2e3max_sweep_duration))

            println("Early save initiated for sweep #$i.\n")
            verbose && println("Current time: ", Dates.format(now(), "d.u yyyy HH:MM"))
            verbose && println("Target time:  ", Dates.format(safe_before, "d.u yyyy HH:MM"))
            save(resumable_filename, mc, overwrite = overwrite, rename = false)
            verbose && println("\nEarly save finished")

            return false
        elseif (now() - last_checkpoint) > safe_every
            verbose && println("Performing scheduled save.")
            last_checkpoint = now()
            save(resumable_filename, mc, overwrite = overwrite, rename = false)
        end
    end
    finish!(mc.measurements, mc, mc.model)

    end_time = now()
    verbose && println("Ended: ", Dates.format(end_time, "d.u yyyy HH:MM"))
    verbose && @printf("Duration: %.2f minutes", (end_time - start_time).value/1000. /60.)

    return true
end


text:
printstyled("  processing of NodalData\n", color=:light_green)

fstr = path*"/SampleFiles/Nodal/Node1_UTC_20200307_170738.006.tdms"
S1 = read_nodal("silixa", fstr)

# these should all work
for f in (:convert_seis!, :demean!, :detrend!, :sync!, :taper!, :ungap!, :unscale!)
  printstyled(string("    ", f, "\n"), color=:light_green)
  getfield(SeisIO, f)(S1)
end

@test_throws ErrorException merge!(S1)

# test resampling
printstyled(string("    ", "resample!", "\n"), color=:light_green)
f0 = 500.
S2 = resample(S1,f0)
resample!(S1,f0)
@test S2.data == S1.data

# test filtering
printstyled(string("    ", "filtfilt!", "\n"), color=:light_green)
S2 = filtfilt(S1,rt="Bandpass",fl=100.,fh=200.)
filtfilt!(S1,rt="Bandpass",fl=100.,fh=200.)
@test S2.data == S1.data


text:
using JSON
using BSON

function pad_sequences(l, maxlen=500)
    if length(l) <= maxlen
        res = zeros(maxlen - length(l))
        for ele in l
            push!(res, ele)
        end
        return res
    else
        length(1) > maxlen || throw("String length exceeds maximum length")
        end
end

function read_weights(filename=sentiment_weights)
    return BSON.load(filename)
end

function read_word_ids(filename=sentiment_words)
    return JSON.parse(String(read(open(filename, "r"))))
end

function embedding(embedding_matrix, x)
    temp = embedding_matrix[:, Int64(x[1])+1]
    for i=2:length(x)
        temp = hcat(temp, embedding_matrix[:, Int64(x[i])+1])
    end
    return reshape(temp, reverse(size(temp)))
end

function flatten(x)
    l = prod(size(x))
    x = permutedims(x, reverse(range(1, length=ndims(x))))
    return reshape(x, (l, 1))
end

function get_sentiment(handle_unknown, ip::Array{T, 1}, weight, rwi) where T <: AbstractString
    model = (x,) -> begin
    a_1 = embedding(weight[:embedding_1]["embedding_1"]["embeddings:0"], x)
    a_2 = flatten(a_1)
    a_3 = Flux.Dense(weight[:dense_1]["dense_1"]["kernel:0"], weight[:dense_1]["dense_1"]["bias:0"], Flux.relu)(a_2)
    a_4 = Flux.Dense(weight[:dense_2]["dense_2"]["kernel:0"], weight[:dense_2]["dense_2"]["bias:0"], Flux.sigmoid)(a_3)
    return a_4
    end
    res = Array{Int, 1}()
    for ele in ip
	if ele in keys(rwi) && rwi[ele] <= size(weight[:embedding_1]["embedding_1"]["embeddings:0"])[2]   # there are only 5000 unique embeddings
            push!(res, rwi[ele])
	else
	    for words in handle_unknown(ele)
		if words in keys(rwi) && rwi[words] <= size(weight[:embedding_1]["embedding_1"]["embeddings:0"])[2]
		    push!(res, rwi[words])
		end
	    end

	end
    end
    return model(pad_sequences(res))[1]
end

struct SentimentModel
    weight
    words

    function SentimentModel()
        new(read_weights(), read_word_ids())
    end
end

struct SentimentAnalyzer
    model::SentimentModel

    SentimentAnalyzer() = new(SentimentModel())
end

function Base.show(io::IO, s::SentimentAnalyzer)
    print(io, "Sentiment Analysis Model Trained on IMDB with a $(length(s.model.words)) word corpus")
end

function(m::SentimentModel)(handle_unknown, text::Array{T, 1}) where T <: AbstractString
    return get_sentiment(handle_unknown, text, m.weight, m.words)
end

"""
    model = SentimentAnalyzer()
    model(doc)
    model(doc, handle_unknown)

Predict sentiment of the input doc in range 0 to 1, 0 being least sentiment score and 1 being the highest.

# Arguments
-  doc              = Input Document for calculating document (`AbstractDocument` type)
-  handle_unknown   = A function for handling unknown words. Should return an array (default x->tuple())
"""
function(m::SentimentAnalyzer)(d::AbstractDocument, handle_unknown = x->tuple())
    m.model(handle_unknown, tokens(d))
end


text:
function ismonotonic(A::AA{T,1}) where T<:AF
    return (all(diff(A) .>= zero(T)) | all(diff(A) .<= zero(T)))
end

function strip_columns(A::AA{T,2}, cols::AA{Bool,1}) where T<:AF
    @assert length(cols) == size(A,2)
    return A[:, .!cols]
end

# credit: traktofon @ https://discourse.julialang.org/t/findnearest-function/4143/4
function searchsortednearest(a::AbstractVector{T}, x::T) where T
    idx = searchsortedfirst(a,x)
    if (idx==1); return idx; end
    if (idx>length(a)); return length(a); end
    if (a[idx]==x); return idx; end
    if (abs(a[idx]-x) < abs(a[idx-1]-x))
        return idx
    else
        return idx-1
    end
end

function searchsortednearest(x::T, a::AbstractVector{T}) where T
    return searchsortednearest(a, x)
end

"""
    calc_rms(A)

Calculate the RMS of given data about the observed mean.

# Arguments
- `A::AbstractArray{Float64,1}`: 1D array containing data for RMS calculation
"""
function calc_rms(A::AA{T,1}) where T<:AF
    return sqrt(sum((A.-mean(A)).^2)/length(A))
end

function strip_nans_by_column(A::AA{T,1}) where T<:AF
    nans = isnan.(A)
    return A[.!vec(nans)]
end

function strip_nans_by_column(A::AA{T,2}) where T<:AF
    nancols = all(isnan.(A), dims=1)
    return A[:,.!vec(nancols)]
end

function strip_nans_by_row(A::AA{T,2}) where T<:AF
    nanrows = all(isnan.(A), dims=2)
    return A[.!vec(nanrows),:]
end

function findnearest(A::AA{T,1}, t) where T<:Real
    return findmin(abs.(A.-t))[2]
end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
module TestImb

using Lale
using DataFrames
using PyCall
using Test

const Hyperopt = laleoperator("Hyperopt","lale")

function testpipe(pipe,train_X,train_y,test_X,test_y)
   optimizer = Hyperopt(pipe, max_evals = 10, scoring="roc_auc")
   trained_optimizer = fit(optimizer,train_X, train_y)
   predictions = predict(trained_optimizer,test_X)
   lalesummary(trained_optimizer)
   @test score(:accuracy,test_y,predictions) > 0.50
end

function imb_test()
   fetch = pyimport("lale.datasets.openml").fetch
   (train_X, train_y), (test_X, test_y) = fetch("breast-cancer", "classification", preprocess=true)
   train_X = DataFrame(train_X,:auto)
   test_X = DataFrame(test_X,:auto)

   MinMaxScaler = laleoperator("MinMaxScaler")
   PCA = laleoperator("PCA")
   NoOp = laleoperator("NoOp","lale")
   LR= laleoperator("LogisticRegression")
   RF= laleoperator("RandomForestClassifier")
   SMOTE = laleoperator("SMOTE","imblearn")
   ADASYN = laleoperator("ADASYN","imblearn")
   CondensedNearestNeighbour = laleoperator("CondensedNearestNeighbour","imblearn")
   SMOTEENN = laleoperator("SMOTEENN","imblearn")
   accuracy_score=pyimport("sklearn.metrics").accuracy_score

   basepipeline =  MinMaxScaler() >> PCA() >> RF()
   testpipe(basepipeline,train_X,train_y,test_X,test_y)

   # SMOTE
   correctedpipeline =  SMOTE(basepipeline)
   testpipe(correctedpipeline,train_X,train_y,test_X,test_y)

   # CondensedNearestNeighbour
   correctedpipeline =  CondensedNearestNeighbour(basepipeline)
   testpipe(correctedpipeline,train_X,train_y,test_X,test_y)

   # SMOTEENN
   correctedpipeline =  SMOTEENN(basepipeline)
   testpipe(correctedpipeline,train_X,train_y,test_X,test_y)

   # ADASYN
   correctedpipeline =  ADASYN(basepipeline)
   testpipe(correctedpipeline,train_X,train_y,test_X,test_y)

end
@testset "Imbalanced Dataset" begin
   imb_test()
end

end


text:
# ---
# title: 397. Integer Replacement
# id: problem397
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, Bit Manipulation
# link: <https://leetcode.com/problems/integer-replacement/description/>
# hidden: true
# ---
#
# Given a positive integer `n`, you can apply one of the following operations:
#
#   1. If `n` is even, replace `n` with `n / 2`.
#   2. If `n` is odd, replace `n` with either `n + 1` or `n - 1`.
#
# Return _the minimum number of operations needed for`n` to become `1`_.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 8
#     Output: 3
#     Explanation: 8 -> 4 -> 2 -> 1
#
#
# **Example 2:**
#
#
#
#     Input: n = 7
#     Output: 4
#     Explanation: 7 -> 8 -> 4 -> 2 -> 1
#     or 7 -> 6 -> 3 -> 2 -> 1
#
#
# **Example 3:**
#
#
#
#     Input: n = 4
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 231 - 1`
#
#
## @lc code=start
using LeetCode

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


text:
function unfold_particle_filter(num_particles::Int, locations::Array{Coordinate,1}, movement_times::Array{Float64,1}, num_samples::Int)
    #init_obs = Gen.choicemap((:chain => 1 => :next_location, locations[1]))
    #init_obs[:chain => 1 => :find_best => :movement_time] = movement_times[1] #and movement times
    init_obs = Gen.choicemap()
    state = Gen.initialize_particle_filter(unfold_model, (0,), init_obs, num_particles)
    for i = 1:num_particles
        println(state.log_weights[i])
    end

    for t=1:length(locations)
        println(t)
        ess = effective_sample_size(normalize_weights(state.log_weights)[2])
        println("ess ", ess)
        maybe_resample!(state, ess_threshold=num_particles)#used to be /2. now always resampling becasue I want to get rid of -Inf before they become NANs
        obs = Gen.choicemap((:chain => t => :next_location, locations[t])) #put the location
        obs[:chain => t => :find_best => :movement_time] = movement_times[t] #and movement times
        Gen.particle_filter_step!(state, (t,), (UnknownChange(),), obs)
        for i = 1:num_particles
            println(state.log_weights[i])
            #println(state.traces[i][:chain => t => :find_best => :distracted])
            println(state.traces[i][:chain => t])
        end
    end

    # return a sample of traces from the weighted collection:
    return Gen.sample_unweighted_traces(state, num_samples)
end

function effective_sample_size(log_normalized_weights::Vector{Float64})
    log_ess = -logsumexp(2. * log_normalized_weights)
    return exp(log_ess)
end

function normalize_weights(log_weights::Vector{Float64})
    log_total_weight = logsumexp(log_weights)
    log_normalized_weights = log_weights .- log_total_weight
    return (log_total_weight, log_normalized_weights)
end

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
using AxisAlgorithms
using Test, LinearAlgebra, SparseArrays, Random

include("tridiag.jl")
include("matmul.jl")
include("woodbury.jl")

max_stars_repo_name: timholy/FilterMD

max_issues_repo_name: timholy/FilterMD

max_forks_repo_name: timholy/FilterMD


text:
include("../examples/lo1.jl")
include("../examples/cqo1.jl")
include("../examples/qo1.jl")
include("../examples/qcqo1.jl")
include("../examples/milo1.jl")
include("../examples/sdo1.jl")
include("../examples/eo1.jl")
include("../examples/nlo1.jl")
include("../examples/solvebasis.jl")
include("../examples/portfolio.jl")
include("../examples/production.jl")
include("../examples/callback.jl")

push!(ARGS, "$(dirname(dirname(@__FILE__)))/examples/feasrepair.lp")
include("../examples/feasrepairex1.jl")


text:
using CRNT
using Documenter

makedocs(;
    modules=[CRNT],
    authors="Laura Brustenga i Moncusí <brust@math.ku.dk> and contributors",
    repo="https://github.com/LauraBMo/CRNT.jl/blob/{commit}{path}#L{line}",
    sitename="CRNT.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://LauraBMo.github.io/CRNT.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/LauraBMo/CRNT.jl",
)


text:
# # Support Vector Machine
#

using Distributions
using KernelFunctions
using LIBSVM
using LinearAlgebra
using Plots
using Random

## Set plotting theme
theme(:wong)

## Set seed
Random.seed!(1234);

# Number of samples:
N = 100;

# Select randomly between two classes:
y_train = rand([-1, 1], N);

# Random attributes for both classes:
X = Matrix{Float64}(undef, 2, N)
rand!(MvNormal(randn(2), I), view(X, :, y_train .== 1))
rand!(MvNormal(randn(2), I), view(X, :, y_train .== -1));
x_train = ColVecs(X);

# Create a 2D grid:
test_range = range(floor(Int, minimum(X)), ceil(Int, maximum(X)); length=100)
x_test = ColVecs(mapreduce(collect, hcat, Iterators.product(test_range, test_range)));

# Create kernel function:
k = SqExponentialKernel() ∘ ScaleTransform(2.0)

# [LIBSVM](https://github.com/JuliaML/LIBSVM.jl) can make use of a pre-computed kernel matrix.
# KernelFunctions.jl can be used to produce that.
# Precomputed matrix for training (corresponds to linear kernel)
model = svmtrain(kernelmatrix(k, x_train), y_train; kernel=LIBSVM.Kernel.Precomputed)

# Precomputed matrix for prediction
y_pr, _ = svmpredict(model, kernelmatrix(k, x_train, x_test));

# Compute prediction on a grid:
contourf(test_range, test_range, y_pr)
scatter!(X[1, :], X[2, :]; color=y_train, lab="data", widen=false)


text:
function correctness_check(
            ring, exps, coeffs_zz, gbexps, gbcoeffs_qq,
            moduli, randomized, rng)
    if randomized_correctness_check(ring, exps, coeffs_zz, gbexps, gbcoeffs_qq, moduli, rng)
        @info "Randomized check passed!"
        # if guaranteed checks are required
        if !randomized
            if guaranteed_correctness_check(ring, exps, coeffs_zz, gbexps, gbcoeffs_qq, moduli, rng)
                @info "Guaranteed checks pass!"
                return true
            end
            @info "Guaranteed checks didn't pass."
            return false
        else
            return true
        end
    end
    @info "Randomized check failed."
    return false
end

function randomized_correctness_check(
            ring, exps, coeffs_zz, gbexps,
            gbcoeffs_qq, moduli, rng)

    goodprime = nextgoodprime(coeffs_zz, moduli, 2^30)
    # @warn "goodprime" coeffs_zz moduli goodprime gbcoeffs_qq
    ring_ff, coeffs_ff = reduce_modulo(coeffs_zz, ring, goodprime)

    gbcoeffs_zz = scale_denominators!(gbcoeffs_qq)

    # @warn "gbzz" gbcoeffs_zz

    ring_ff, gbcoeffs_ff = reduce_modulo(gbcoeffs_zz, ring, goodprime)

    # check that initial ideal contains in the computed groebner basis modulo goodprime
    expsreduced, _ = normal_form_f4(ring_ff, gbexps, gbcoeffs_ff, exps, coeffs_ff, rng)
    for es in expsreduced
        if !isempty(es)
            return false
        end
    end

    # check that the basis is groebner basis modulo goodprime
    if !isgroebner_f4(ring_ff, gbexps, gbcoeffs_ff, rng)
        return false
    end

    return true
end

function guaranteed_correctness_check(
            ring, exps, coeffs_zz, gbexps,
            gbcoeffs_qq, moduli, rng)

end


text:
# SPDX-License-Identifier: X11
# 2020-11-21 (Timed)
# Simple Math (300pt)

const MODULUS = 998244353

function main()
  A,B,C = parse.(Int, (split ∘ readline)())
  println(mod(sumto(A) * sumto(B) * sumto(C), MODULUS))
end

sumto(n::Integer)::BigInt = n * (n + 1) ÷ 2

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
module CSTParser
global debug = true

using Tokenize
import Base: length, first, last, getindex, setindex!
import Tokenize.Tokens
import Tokenize.Tokens: RawToken, AbstractToken, iskeyword, isliteral, isoperator, untokenize
import Tokenize.Lexers: Lexer, peekchar, iswhitespace

export ParseState, parse_expression

include("lexer.jl")
include("spec.jl")
include("utils.jl")
include("recovery.jl")
include("components/internals.jl")
include("components/keywords.jl")
include("components/lists.jl")
include("components/operators.jl")
include("components/strings.jl")
include("location.jl")
include("conversion.jl")
include("display.jl")
include("interface.jl")

"""
    parse_expression(ps)

Parses an expression until `closer(ps) == true`. Expects to enter the
`ParseState` the token before the the beginning of the expression and ends
on the last token.

Acceptable starting tokens are:
+ A keyword
+ An opening parentheses or brace.
+ An operator.
+ An instance (e.g. identifier, number, etc.)
+ An `@`.

"""
@addctx :expr function parse_expression(ps::ParseState)
    if ps.nt.kind == Tokens.COMMA
        push!(ps.errors, Error((ps.nt.startbyte:ps.nws.endbyte) .+ 1, "Expression began with a comma."))
        ret = ErrorToken(PUNCTUATION(next(ps)))
    elseif ps.nt.kind ∈ term_c && ps.nt.kind != Tokens.END
        push!(ps.errors, Error((ps.nt.startbyte:ps.nws.endbyte) .+ 1, "Expression began with a terminal token: $(ps.nt.kind)."))
        ret = ErrorToken(INSTANCE(next(ps)))
    else
        next(ps)
        if iskeyword(ps.t.kind) && ps.t.kind != Tokens.DO
            ret = parse_kw(ps)
        elseif ps.t.kind == Tokens.LPAREN
            ret = parse_paren(ps)
        elseif ps.t.kind == Tokens.LSQUARE
            ret = @default ps parse_array(ps)
        elseif ps.t.kind == Tokens.LBRACE
            ret = @default ps @closebrace ps parse_braces(ps)
        elseif isinstance(ps.t) || isoperator(ps.t)
            if ps.t.kind == Tokens.WHERE
                ret = IDENTIFIER(ps)
            else
                ret = INSTANCE(ps)
            end
            if is_colon(ret) && ps.nt.kind != Tokens.COMMA
                ret = parse_unary(ps, ret)
            end
        elseif ps.t.kind == Tokens.AT_SIGN
            ret = parse_macrocall(ps)
        else
            ret = ErrorToken(INSTANCE(ps))
            push!(ps.errors, Error((ps.nt.startbyte:ps.nws.endbyte) .+ 1, "Expression began with a : $(ps.nt.kind)."))
        end

        while !closer(ps)
            ret = parse_compound(ps, ret)
        end
    end
    return ret
end

function parse_compound(ps::ParseState, @nospecialize ret)
    if ps.nt.kind == Tokens.FOR
        ret = parse_generator(ps, ret)
    elseif ps.nt.kind == Tokens.DO
        ret = @default ps @closer ps block parse_do(ps, ret)
    elseif isajuxtaposition(ps, ret)
        op = OPERATOR(0, 0, Tokens.STAR, false)
        ret = parse_operator(ps, ret, op)
    elseif (ret isa EXPR{x_Str} ||  ret isa EXPR{x_Cmd}) && ps.nt.kind == Tokens.IDENTIFIER
        arg = IDENTIFIER(next(ps))
        push!(ret, LITERAL(arg.fullspan, arg.span, val(ps.t, ps), Tokens.STRING))
    elseif (ret isa IDENTIFIER || (ret isa BinarySyntaxOpCall && is_dot(ret.op))) && (ps.nt.kind == Tokens.STRING || ps.nt.kind == Tokens.TRIPLE_STRING || ps.nt.kind == Tokens.CMD)
        next(ps)
        arg = parse_string_or_cmd(ps, ret)
        head = arg.kind == Tokens.CMD ? x_Cmd : x_Str
        ret = EXPR{head}(Any[ret, arg])
    elseif ps.nt.kind == Tokens.LPAREN
        no_ws = !isemptyws(ps.ws)
        err_rng = ps.t.endbyte + 2:ps.nt.startbyte
        ret = @closeparen ps parse_call(ps, ret)
        if no_ws && !(ret isa UnaryOpCall || ret isa UnarySyntaxOpCall)
            push!(ps.errors, Error(err_rng, "White space in function call."))
            ret = ErrorToken(ret)
        end
    elseif ps.nt.kind == Tokens.LBRACE
        if isemptyws(ps.ws)
            ret = @default ps @nocloser ps inwhere @closebrace ps parse_curly(ps, ret)
        else
            push!(ps.errors, Error(ps.t.endbyte + 2:ps.nt.startbyte , "White space in brace call."))
            ret = ErrorToken(@default ps @nocloser ps inwhere @closebrace ps parse_curly(ps, ret))
        end
    elseif ps.nt.kind == Tokens.LSQUARE && isemptyws(ps.ws) && !(ret isa OPERATOR)
        ret = @default ps @nocloser ps block parse_ref(ps, ret)
    elseif ps.nt.kind == Tokens.COMMA
        ret = parse_tuple(ps, ret)
    elseif isunaryop(ret) && ps.nt.kind != Tokens.EQ
        ret = parse_unary(ps, ret)
    elseif isoperator(ps.nt)
        op = OPERATOR(next(ps))
        ret = parse_operator(ps, ret, op)
    elseif ret isa UnarySyntaxOpCall && is_prime(ret.arg2)
        # prime operator followed by an identifier has an implicit multiplication
        nextarg = @precedence ps 11 parse_expression(ps)
        ret = BinaryOpCall(ret, OPERATOR(0, 0, Tokens.STAR,false), nextarg)
################################################################################
# Everything below here is an error
################################################################################
    elseif ps.nt.kind in (Tokens.RPAREN, Tokens.RSQUARE, Tokens.RBRACE)
        push!(ps.errors, Error((ps.t.startbyte:ps.nt.endbyte) .+ 1 , "Disallowed compound expression."))
        ret = EXPR{ErrorToken}([ret, ErrorToken(PUNCTUATION(next(ps)))])
    else
        push!(ps.errors, Error((ps.t.startbyte:ps.nt.endbyte) .+ 1 , "Disallowed compound expression."))
        nextarg = parse_expression(ps)
        ret = EXPR{ErrorToken}([ret, nextarg])
    end
    return ret
end

"""
    parse_paren(ps, ret)

Parses an expression starting with a `(`.
"""
@addctx :paren function parse_paren(ps::ParseState)
    args = Any[PUNCTUATION(ps)]
    @closeparen ps @default ps @nocloser ps inwhere parse_comma_sep(ps, args, false, true, true)

    if length(args) == 2 && ((ps.ws.kind != SemiColonWS || (length(args) == 2 && args[2] isa EXPR{Block})) && !(args[2] isa EXPR{Parameters}))
        accept_rparen(ps, args)
        ret = EXPR{InvisBrackets}(args)
    else
        accept_rparen(ps, args)
        ret = EXPR{TupleH}(args)
    end
    return ret
end

"""
    parse(str, cont = false)

Parses the passed string. If `cont` is true then will continue parsing until the end of the string returning the resulting expressions in a TOPLEVEL block.
"""
function parse(str::String, cont = false)
    ps = ParseState(str)
    x, ps = parse(ps, cont)
    return x
end

function parse_doc(ps::ParseState)
    if (ps.nt.kind == Tokens.STRING || ps.nt.kind == Tokens.TRIPLE_STRING) && !isemptyws(ps.nws)
        doc = LITERAL(next(ps))
        if (ps.nt.kind == Tokens.ENDMARKER || ps.nt.kind == Tokens.END)
            return doc
        elseif isbinaryop(ps.nt) && !closer(ps)
            ret = parse_compound(ps, doc)
            return ret
        end

        ret = parse_expression(ps)
        ret = EXPR{MacroCall}(Any[GlobalRefDOC, doc, ret])
    elseif ps.nt.kind == Tokens.IDENTIFIER && val(ps.nt, ps) == "doc" && (ps.nnt.kind == Tokens.STRING || ps.nnt.kind == Tokens.TRIPLE_STRING)
        doc = IDENTIFIER(next(ps))
        next(ps)
        arg = parse_string_or_cmd(ps, doc)
        doc = EXPR{x_Str}(Any[doc, arg])
        ret = parse_expression(ps)
        ret = EXPR{MacroCall}(Any[GlobalRefDOC, doc, ret])
    else
        ret = parse_expression(ps)
    end
    return ret
end

function parse(ps::ParseState, cont = false)
    if ps.l.io.size == 0
        return (cont ? EXPR{FileH}(Any[]) : nothing), ps
    end
    last_line = 0
    curr_line = 0

    if cont
        top = EXPR{FileH}(Any[])
        if ps.nt.kind == Tokens.WHITESPACE || ps.nt.kind == Tokens.COMMENT
            next(ps)
            push!(top, LITERAL(ps.nt.startbyte, ps.nt.startbyte, "", Tokens.NOTHING))
        end

        while !ps.done && !ps.errored
            curr_line = ps.nt.startpos[1]
            ret = parse_doc(ps)

            # join semicolon sep items
            if curr_line == last_line && last(top.args) isa EXPR{TopLevel}
                push!(last(top.args), ret)
            elseif ps.ws.kind == SemiColonWS
                push!(top, EXPR{TopLevel}(Any[ret]))
            else
                push!(top, ret)
            end
            last_line = curr_line
        end
    else
        if ps.nt.kind == Tokens.WHITESPACE || ps.nt.kind == Tokens.COMMENT
            next(ps)
            top = LITERAL(ps.nt.startbyte, ps.nt.startbyte, "", Tokens.NOTHING)
        else
            top = parse_doc(ps)
            last_line = ps.nt.startpos[1]
            if ps.ws.kind == SemiColonWS
                top = EXPR{TopLevel}(Any[top])
                while ps.ws.kind == SemiColonWS && ps.nt.startpos[1] == last_line && ps.nt.kind != Tokens.ENDMARKER
                    ret = parse_doc(ps)
                    push!(top, ret)
                    last_line = ps.nt.startpos[1]
                end
            end
        end
    end

    return top, ps
end

function parse_file(path::String)
    x = parse(read(path, String), true)
    File([], [], path, x, [])
end

function parse_directory(path::String, proj = Project(path, []))
    for f in readdir(path)
        if isfile(joinpath(path, f)) && endswith(f, ".jl")
            try
                push!(proj.files, parse_file(joinpath(path, f)))
            catch
                println("$f failed to parse")
            end
        elseif isdir(joinpath(path, f))
            parse_directory(joinpath(path, f), proj)
        end
    end
    proj
end
end


text:
using Halide
using Test

@testset "Halide.jl" begin
    # Write your tests here.
end


text:
#==============================================================================

		B_Switch Code

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

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

    Now we are going to construct a matrix summarizing the evolution of V_z

    This comes from the following discretized Bellman equation:

    ρv_ij = u(c_ij) + v_k(zF(k_{i}) -δk-c) + v_z⋅μ(z)
                                   + 1/2(v_{zz})σ^2(z)

                                   or

    ρv_ij = u(c_ij) + v_k(zF(k_{i}) -δk-c) + ((v_{i,j+1}-v_{i,j})/Δz)μ(z)
                            + 1/2((v_{i,j+1}-2v_{i,j}+v_{ij-1})/Δz^2)σ^2(z)

    Assume forward difference because of boundary conditions

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

function B_switch(μ, Σ_sq, dz, H)
	 dz_sq = dz^2
	 yy = (-Σ_sq/dz_sq .- μ/dz)
	 χ = Σ_sq/(2*dz_sq)
	 ζ = μ/dz .+ Σ_sq/(2*dz_sq)

	 # Define the Diagonals of this matrix
	 updiag_z = zeros(H,1)
	 	for j = 1:J
			updiag_z =[updiag_z; repeat([ζ[j]], H, 1)]
		end
	 updiag_z =(updiag_z[:])

	 centerdiag_z=repeat([χ[1]+yy[1]],H,1)
		for j = 2:J-1
			centerdiag_z = [centerdiag_z; repeat([yy[j]], H, 1)]
		end
	 centerdiag_z=[centerdiag_z; repeat([yy[J]+ζ[J]], H, 1)]
	  centerdiag_z = centerdiag_z[:]

	lowdiag_z = repeat([χ[2]], H, 1)
		for j=3:J
			lowdiag_z = [lowdiag_z; repeat([χ[j]],H,1)]
		end
	lowdiag_z=lowdiag_z[:]

	# spdiags in Matlab allows for automatic trimming/adding of zeros
	    # spdiagm does not do this

	B_switch = sparse(Diagonal(centerdiag_z))+ [zeros(H,H*J);  sparse(Diagonal(lowdiag_z)) zeros(H*(J-1), H)]+ sparse(Diagonal(updiag_z))[(H+1):end,1:(H*J)] # trim off rows of zeros

return B_switch

end


text:
using Distributions, Random, JuMP, Ipopt
using LinearAlgebra
using Suppressor

struct ExperimentData
    y::Array{Float64}
    x
    ξ
end

function frankel_generator(β, ξ)
    n, d = size(ξ)
    β = β[2:(d+1)]
    γ = rand(Uniform(0, 8), n)
    m = 1
    z = rand(Normal(), (n, d))
    e = rand(Normal(), n)
    x = z .+ m.*γ.*(ξ.+β')
    y = z*ones(d) .+ e
    return ExperimentData(y, x, ξ)
end

function mock_survey_coder(β, ξ, recode, n)
    m = length(recode)
    γ = rand(Uniform(0, 8), n)
    z = rand(1:m, n)
    x = Int.(floor.(z .+ γ.*(ξ.+β)))
    x[x.>m] .= m
    x[x .< m] .=1
    x = recode[x]
    z = recode[z]
    return z, x
end

function survey_generator(β, ξ)
    n, d = size(ξ)
    β = β[2:(d+1)]
    z1, x1 = mock_survey_coder(β[1], ξ[:,1], [21.0, 29.0, 39.0, 49.0, 59.0, 69.0, 79.0, 85.0], n)
    z2, x2 = mock_survey_coder(β[2], ξ[:,2], [10, 12, 13, 14, 16, 18, 20], n)
    z3, x3 = mock_survey_coder(β[3], ξ[:,3], [0, 1], n)
    #y = floor.((10.0 .+ [z1 z2 z3]*[0.1, 0.12, 10.0])/10.0)*10
    y = 10.0 .+ [z1 z2 z3]*[0.1, 0.12, 10.0]
    return ExperimentData(y, [x1 x2 x3], ξ)
end

function update_β_ols(data::ExperimentData, β, α=0)
    β_n = linReg(data.y, data.x)
    return β_n
end

function mse(data, β)
    return mean((data.y .- [ones(length(data.y)) data.x]*β).^2)
end

#20000, 100, 4
function update_β_robust(data::ExperimentData, β, α=[0.5, 1/80000.0, 1/150.0, 1/3.0])
    n, d = size(data.ξ)
    yhat = β[1] .+ sum((data.ξ .+ β[2:(d+1)]').*data.x, dims=2)

    ehat = data.y .- yhat
    dβ = [-2/n*sum((data.y .- yhat)), linReg(ehat.^2, data.ξ)[2:(d+1)]...]
    println(dβ)
    β_n = β .- α.*dβ
    β_n[1] = mean(data.y)  - mean(sum(β_n[2:(d+1)]'.*data.x, dims=2))
    return β_n
end

function run_simulation(generator, updater, β0::Array{Float64, 1}, n::Int, steps::Int)
    d = length(β0) - 1
    β_history = zeros(Float64, (steps, d+1))
    mse_history = zeros(Float64, steps)
    β_history[1, :] = β0
    for s in 2:steps
        η = β_history[s-1, 2:(d+1)]
        ξ = rand([-1,1], (n, d)).*η'
        data = generator(β_history[s-1, :], ξ)
        mse_history[s-1] = mse(data, β_history[s-1, :])
        β_history[s, :] = updater(data, β_history[s-1, :])
    end

    mse_history[steps] = mse_history[steps-1]
    return β_history, mse_history
end

function fk_solution()
    γ = rand([0, 0.5, 8], n)
    z = rand(Normal(), n)
    e = rand(Normal(0, 0.5), n)
    y = z .+ e
    model =  Model(Ipopt.Optimizer)
    @variable(model, b1)
    @variable(model, b0)
    @NLobjective(model, Min, 1/n*sum((y[i] - b1*(z[i]+γ[i]*b1)-b0)^2 for i in 1:n))
    @suppress_out begin
        optimize!(model)
    end
    return [value(b0) value(b1)]
end

function linReg(y, x)
    x = [ones(length(y)) x ]
    β = inv(x'*x)*x'*y
    return β
end

max_stars_repo_name: evanmunro/personalized-policy

max_issues_repo_name: evanmunro/personalized-policy

max_forks_repo_name: evanmunro/personalized-policy


text:
import FoldsChainRules
import FoldsChainRulesBenchmarks
SUITE = FoldsChainRulesBenchmarks.setup()


text:
# Adapted from https://ccrma.stanford.edu/~jos/filters/Partial_Fraction_Expansion_residued_m.html

using DSP # requires DSP module for filt, equivalent Julia function to filter

function residued(b, a, toler)
if nargin<3, toler=0.001; end
NUM = vec(b);
DEN = vec(a);
nb = length(NUM);
na = length(DEN);
f = [];
if na<=nb
  f = filt(NUM,DEN,[1,zeros(nb-na)]);
  NUM = NUM - conv(DEN,f);
  NUM = NUM(nb-na+2:end);
end
[r,p,f2,e] = residuez(NUM,DEN,toler);

return [r, p, f, e]