text:
using Sundials, Test, SparseArrays, DiffEqOperators

# Test for Jacobian usage
function Lotka(du, u, p, t)
    du[1] = u[1] - u[1] * u[2] # REPL[7], line 3:
    du[2] = -3 * u[2] + 1 * u[1] * u[2]
    nothing
end

jac_called = false
function Lotka_jac(J, u, p, t)
    global jac_called
    jac_called = true
    J[1, 1] = 1.0 - u[2]
    J[1, 2] = -u[1]
    J[2, 1] = 1 * u[2]
    J[2, 2] = -3 + u[1]
    nothing
end

Lotka_f = ODEFunction(Lotka, jac = Lotka_jac)
prob = ODEProblem(Lotka_f, ones(2), (0.0, 10.0))
good_sol = solve(prob, CVODE_BDF())
@test jac_called == true

Lotka_f = ODEFunction(
    Lotka,
    jac = Lotka_jac,
    jac_prototype = sparse([1, 2, 1, 2], [1, 1, 2, 2], zeros(4)),
)

prob = ODEProblem(Lotka_f, ones(2), (0.0, 10.0))
jac_called = false
sol9 = solve(prob, CVODE_BDF(linear_solver = :KLU))
@test jac_called == true
@test Array(sol9) ≈ Array(good_sol)

Lotka_fj = ODEFunction(Lotka, jac_prototype = JacVecOperator{Float64}(Lotka, ones(2)))

prob = ODEProblem(Lotka_fj, ones(2), (0.0, 10.0))
sol9 = solve(prob, CVODE_BDF(linear_solver = :GMRES))

function f2!(res, du, u, p, t)
    res[1] = 1.01du[1]
    return
end

jac_called = false
function f2_jac!(out, du, u, p, gamma, t)
    global jac_called
    jac_called = true
    out[1] = 1.01
end

f2_f = DAEFunction(f2!, jac = f2_jac!)

u0 = [0.0]
tspan = (0.0, 10.0)
du0 = [0.0]
dae_prob = DAEProblem(f2_f, du0, u0, tspan, differential_vars = [true])
good_sol = solve(dae_prob, IDA())
@test jac_called == true

function testjac(res, du, u, p, t)
    res[1] = du[1] - 1.5 * u[1] + 1.0 * u[1] * u[2]
    res[2] = du[2] + 3 * u[2] - u[1] * u[2]
end

jac_called = false
function testjac_jac(J, du, u, p, gamma, t)
    global jac_called
    jac_called = true
    J[1, 1] = gamma - 1.5 + 1.0 * u[2]
    J[1, 2] = 1.0 * u[1]
    J[2, 1] = -1 * u[2]
    J[2, 2] = gamma + 3 - u[1]
    nothing
end

testjac_f = DAEFunction(testjac, jac = testjac_jac)

prob3 = DAEProblem(
    testjac_f,
    [0.5, -2.0],
    ones(2),
    (0.0, 10.0),
    differential_vars = [true, true],
)
sol3 = solve(prob3, IDA())

@test jac_called == true
jac_called = false
prob4 = DAEProblem(testjac, [0.5, -2.0], ones(2), (0.0, 10.0))
sol4 = solve(prob4, IDA())

@test jac_called == false

println("Jacobian vs no Jacobian difference:")
println(maximum(sol3 - sol4))
@test maximum(sol3 - sol4) < 1e-6


text:
function _replace(s, pairs::Pair...)
    for p in pairs
        s = replace(s, p)
    end
    return s
end

# Convert flisp definitions and comments to psuedo-Julia to reflect the
# structure of the existing flisp parser.
#
# Surrounded with all this compiler technology, but still resorting to a pile
# of regexs? 😂😱
function juliafy_flisp(fl_input, jl_output)
    prev_newline = false
    had_comment = false
    for line in readlines(fl_input)
        if occursin(r"^\(define *\(", line)
            had_comment && println(jl_output, "#")
            println(jl_output, "# flisp: $line")
            m = match(r"\(define *\(([-a-zA-Z?_=0-9*><:!]+) *([^)]*)", replace(line, "-"=>"_"))
            isnothing(m) && @error "no match for line" line
            funcname = m[1]
            funcname = _replace(funcname,
                r"(.*)\?"=>s"is_\1",
                "=" => "equals",
                "*" => "_star",
                ">" => "_gt",
                "<" => "_lt",
                ":" => "_",
            )
            funcargs = _replace(m[2],
                r" *\("                 => ";",
                r" +"                   => ", ",
                "."                     => "_",
                r"([-a-zA-Z?_=]+)\?"    => s"is_\1",
                r", *#t"                => "=true",
                r", *#f"                => "=false",
                ";"                     => "; ",
            )
            if startswith(funcname, "parse_")
                funcargs = replace(funcargs, r"^ *s\b"=>"ps::ParseState")
            end
            text = """
            function $funcname($funcargs)
                TODO("$funcname unimplemented")
            end
            """
            ex = Meta.parse(text, raise=false)
            if Meta.isexpr(ex, :error)
                @warn "Generated bad code"  message=ex.args[1] code=Text(text)
            end
            print(jl_output, text)
            prev_newline = false
            had_comment = false
        elseif occursin(r"^;", line)
            println(jl_output, replace(line, r"^;+" => "#"))
            prev_newline = false
            had_comment = true
        elseif line == ""
            if !prev_newline
                println(jl_output)
            end
            prev_newline = true
            had_comment = false
        end
    end
end

open("/home/chris/dev/julia/src/julia-parser.scm", "r") do fl_input
    open(joinpath(@__DIR__, "julia_parser_scm.jl"), "w") do jl_output
        juliafy_flisp(fl_input, jl_output)
    end
end

ext: jl


text:
import NonUniformRandomVariateGeneration.sampleMultinomial!

@inline function _iotaParallel!(array::Vector{Int64}, N::Int64, nthreads::Int64,
  Nperthread::Int64)
  Threads.@threads for i in 1:nthreads
    start::Int64 = Nperthread * (i - 1) + 1
    finish::Int64 = start + Nperthread - 1
    for i = start:finish
      @inbounds array[i] = i
    end
  end
end

@inline function _cpeHelper(smcio::SMCIO)
  Threads.@threads for i = 1:smcio.nthreads
    ip::_SMCInternalParallel = smcio.internal.parallel
    oldEves::Vector{Int64} = smcio.internal.oldEves
    @inbounds localAs::SubVectorI64 = ip.localAs[i]
    @inbounds localEves::SubVectorI64 = ip.localEves[i]
    _setEves!(localEves, oldEves, localAs)
  end
end

@inline function _copyParticlesEves!(smcio::SMCIO{Particle}) where Particle
  Threads.@threads for i = 1:smcio.nthreads
    ip::_SMCInternalParallel = smcio.internal.parallel
    zetas::Vector{Particle} = smcio.zetas
    @inbounds localZetaAncs::SubVector{Particle} = ip.localZetaAncs[i]
    @inbounds localAs::SubVectorI64 = ip.localAs[i]
    @inbounds localOldEves::SubVectorI64 = ip.localOldEves[i]
    @inbounds localEves::SubVectorI64 = ip.localEves[i]
    _copyParticles!(localZetaAncs, zetas, localAs)
    localOldEves .= localEves
  end
  _cpeHelper(smcio)
end

@inline function _MutateWeightHelper1(p::Int64, smcio::SMCIO)
  ip::_SMCInternalParallel= smcio.internal.parallel
  maxlw::Float64 = maximum(ip.maxlws)
  smcio.internal.maxlw = maxlw

  for i = 1:smcio.nthreads
    @inbounds tmp::Float64 = exp(ip.maxlws[i] - maxlw)
    @inbounds ip.partialSums[i] *= tmp
    @inbounds ip.partialSumSqs[i] *= tmp * tmp
  end

  smcio.internal.sws = sum(ip.partialSums)
  smcio.internal.mws = smcio.internal.sws / smcio.N
  @inbounds smcio.esses[p] = smcio.internal.mws * smcio.internal.mws *
    smcio.N / sum(ip.partialSumSqs)
end

@inline function _MutateWeightHelper2(smcio::SMCIO)
  Threads.@threads for i = 1:smcio.nthreads
    ip::_SMCInternalParallel = smcio.internal.parallel
    @inbounds G::Float64 = exp(ip.maxlws[i] - smcio.internal.maxlw) /
      smcio.internal.mws
    @inbounds ip.localWs[i] .*= G
  end
end

@inline function _MutateWeight!(p::Int64, model::SMCModel,
  smcio::SMCIO{Particle, ParticleScratch},
  ref::Union{Nothing, Particle} = nothing) where {Particle, ParticleScratch}
  Threads.@threads for i = 1:smcio.nthreads
    ip::_SMCInternalParallel{Particle, ParticleScratch} = smcio.internal.parallel
    @inbounds localZetas::SubVector{Particle} = ip.localZetas[i]
    @inbounds localZetaAncs::SubVector{Particle} = ip.localZetaAncs[i]
    @inbounds localLogWs::SubVectorF64 = ip.localLogWs[i]
    @inbounds localWs::SubVectorF64 = ip.localWs[i]
    @inbounds pScratch::ParticleScratch = ip.particleScratches[i]
    engine::RNG = getRNG()

    if i == 1 && ref != nothing
      _mutateParticles!(localZetas, engine, p, model.M!, localZetaAncs,
        pScratch, ref)
    else
      _mutateParticles!(localZetas, engine, p, model.M!, localZetaAncs,
        pScratch)
    end
    @inbounds noresample::Bool = p == 1 || smcio.resample[p-1]
    _logWeightParticles!(localLogWs, p, model.lG, localZetas, pScratch,
      noresample, localWs)

    localmaxlw::Float64 = maximum(localLogWs)
    @inbounds ip.maxlws[i] = localmaxlw
    localWs .= exp.(localLogWs .- localmaxlw)
    @inbounds ip.partialSums[i] = sum(localWs)
    @inbounds ip.partialSumSqs[i] = mapreduce(x -> x * x, +, localWs)
  end

  _MutateWeightHelper1(p, smcio)
  _MutateWeightHelper2(smcio)
end

@inline function _resampleHelper(smcio::SMCIO, conditional::Bool)
  Threads.@threads for i = 1:smcio.nthreads
    ip::_SMCInternalParallel = smcio.internal.parallel
    @inbounds if ip.Ns[i] > 0
      @inbounds localScratch1 = ip.localScratch1s[i]
      @inbounds localWs = ip.localWs[i]
      @inbounds localN = ip.Ns[i]
      engine::RNG = getRNG()
      start::Int64 = ip.Nperthread * (i - 1)
      @inbounds offset::Int64 = i == 1 ? 0 : ip.NsPartial[i-1]
      if i == 1 && conditional
        @inbounds smcio.internal.as[1] = 1
        if localN > 1
          sampleSerial!(smcio.internal.as, engine, localWs, localN - 1,
            smcio.internal.scratch2, localScratch1, offset + 1, start)
        end
      else
        sampleSerial!(smcio.internal.as, engine, localWs, localN,
          smcio.internal.scratch2, localScratch1, offset, start)
      end
    end
  end
end

@inline function _resampleParallel!(smcio::SMCIO, p::Int64, conditional::Bool)
  if p != smcio.n-1 && smcio.esses[p] > smcio.essThreshold
    smcio.resample[p] = false
    _iotaParallel!(smcio.internal.as, smcio.N, smcio.nthreads,
      smcio.internal.parallel.Nperthread)
    return
  end
  @inbounds smcio.resample[p] = true
  smcio.internal.nresamples += 1

  ip::_SMCInternalParallel = smcio.internal.parallel
  Ns::Vector{Int64} = ip.Ns
  NsPartial::Vector{Int64} = ip.NsPartial
  partialSums::Vector{Float64} = ip.partialSums
  partialSums ./= smcio.internal.sws
  if conditional
    sampleMultinomial!(smcio.N-1, partialSums, Ns, getRNG())
    @inbounds Ns[1] += 1
  else
    sampleMultinomial!(smcio.N, partialSums, Ns, getRNG())
  end

  cumsum!(NsPartial,Ns)
  _resampleHelper(smcio, conditional)
end

@inline function _Vhat1ParallelHelper(smcio::SMCIO)
  N::Int64 = smcio.N
  vs::Vector{Float64} = smcio.internal.scratch1

  v::Float64 = 0.0
  for t = 1:smcio.nthreads
    @inbounds v += vs[t]
  end
  q::Int64 = 1 + smcio.internal.nresamples
  return 1 - (1-v)*(N/(N-1))^q
end

@inline function _Vhat1Parallel(smcio::SMCIO)
  # N::Int64 = smcio.N
  # vs::Vector{Float64} = smcio.internal.scratch1
  Threads.@threads for t = 1:smcio.nthreads
    N::Int64 = smcio.N
    vs::Vector{Float64} = smcio.internal.scratch1
    Nperthread::Int64 = smcio.internal.parallel.Nperthread
    eves::Vector{Int64} = smcio.eves
    ws::Vector{Float64} = smcio.ws
    i::Int64 = searchsortedfirst(eves, Nperthread * (t - 1) + 1)
    finish::Int64 = searchsortedlast(eves, Nperthread * t)
    vLocal::Float64 = 0.0
    while i <= finish
      @inbounds e = eves[i]
      tmp::Float64 = 0.0
      @inbounds while (i <= N && eves[i] == e)
        @inbounds tmp += ws[i]
        i += 1
      end
      tmp /= N
      vLocal += tmp*tmp
    end
    vs[t] = vLocal
  end
  return _Vhat1ParallelHelper(smcio)
end

@inline function _intermediateOutputParallel!(smcio::SMCIO{Particle},
  p::Int64) where Particle
  Threads.@threads for i = 1:smcio.nthreads
    ip::_SMCInternalParallel{Particle} = smcio.internal.parallel
    @inbounds _copyParticles!(ip.localAllZetas[p][i], ip.localZetas[i])
    @inbounds ip.localAllWs[p][i] .= ip.localWs[i]
    @inbounds ip.localAllEves[p][i] .= ip.localEves[i]
    @inbounds p < smcio.n && (ip.localAllAs[p][i] .= ip.localAs[i])
  end
end

# Parallel implementation of SMC
@inline function _smcParallel!(model::SMCModel, smcio::SMCIO{Particle},
  ref::Union{Nothing, Vector{Particle}} = nothing,
  refout::Union{Nothing, Vector{Particle}} = nothing) where Particle
  smcio.internal.nresamples = 0
  lZ::Float64 = 0.0
  _iotaParallel!(smcio.eves, smcio.N, smcio.nthreads,
    smcio.internal.parallel.Nperthread)

  for p = 1:smcio.n
    p > 1 && _copyParticlesEves!(smcio)
    if ref == nothing
      _MutateWeight!(p, model, smcio)
    else
      _MutateWeight!(p, model, smcio, ref[p])
    end

    lZ += smcio.internal.maxlw + log(smcio.internal.mws)
    @inbounds smcio.logZhats[p] = lZ

    @inbounds smcio.Vhat1s[p] = _Vhat1Parallel(smcio)

    p < smcio.n && _resampleParallel!(smcio, p, ref != nothing)
    smcio.fullOutput && _intermediateOutputParallel!(smcio, p)
  end

  refout != nothing && SequentialMonteCarlo.pickParticle!(refout, smcio)
  return
end


text:
using MLFlowClient
using Test
using UUIDs
using Dates

function mlflow_server_is_running(mlf::MLFlow)
    try
        response = MLFlowClient.mlfget(mlf, "experiments/list")
        return isa(response, Dict)
    catch e
        return false
    end
end

# creates an instance of mlf
# skips test if mlflow is not available on default location, http://localhost:5000
macro ensuremlf()
    e = quote
        mlf = MLFlow()
        mlflow_server_is_running(mlf) || return nothing
    end
    eval(e)
end

max_stars_repo_name: JuliaAI/MLFlowClient

max_issues_repo_name: JuliaAI/MLFlowClient

max_forks_repo_name: JuliaAI/MLFlowClient


text:
"""
     orderSOPJobsPrecedenceToEdge(model::SequencingModel)

Find a variable ordering for the model based on the edge weights and precedence graph
Return the variable ordering

# Arguments
- `model::SequencingModel`: the sequencing problem to be evaluated
"""
function orderSOPPrecedenceToEdge(model::SequencingModel)
    objType = valtype(model.objective.f.components)
    sortedKeys = orderEdgeWeight(model)
    #put jobs into the order by precedence and then average edge weight
    finalOrder = Vector{objType}()
    followerNumbers = getPrecedenceNumbers(model.precedenceP2F)
    maxLength = maximum(values(followerNumbers))
    while maxLength>0
        for i in 1:length(sortedKeys)
            if sortedKeys[i] in keys(followerNumbers) && followerNumbers[sortedKeys[i]]==maxLength
                push!(finalOrder, sortedKeys[i])
            end
        end
        maxLength = maxLength - 1
    end
    for i in 1:length(sortedKeys)
        if !(sortedKeys[i] in finalOrder)
            push!(finalOrder, sortedKeys[i])
        end
    end
    return finalOrder
end


text:
########################################
## File Name: cost_models.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/13
## Description: Cost Models for SACBP
########################################

using StatsFuns

abstract type CostModel end

struct CostPosRangeLocalization <: CostModel
end
function instCost(model::CostPosRangeLocalization,s::AugState{PhysPos{U},VecBelMvNormal{T}},
                  u::PosControl,Cu::Matrix{Float64}) where T <: Real where U <: Real
    return 1/2*vec(u)'*Cu*vec(u)
end
function instCost_grad_p(model::CostPosRangeLocalization,s::AugState{PhysPos{U},VecBelMvNormal{T}}) where T <: Real where U <: Real
    return zeros(length(vec(s.p)))
end
function instCost_grad_b(model::CostPosRangeLocalization,s::AugState{PhysPos{U},VecBelMvNormal{T}}) where T <: Real where U <: Real
    gVecs = Vector{Vector{T}}(undef, length(s.b));
    for ii = 1:length(s.b)
        gVecs[ii] = zeros(length(s.b[ii].params))
    end
    return gVecs
end
function termCost(model::CostPosRangeLocalization,
                  s::AugState{PhysPos{U},VecBelMvNormal{T}}) where T <: Real where U <: Real
    DArray_end = Distributions.MvNormal.(s.b);
    entropies = Distributions.entropy.(DArray_end);
    weighted = exp.(entropies); # Check coefficient!
    return sum(weighted);
    #return sum(Distributions.entropy.(DArray_end))
end
function termCost_grad_p(model::CostPosRangeLocalization,s::AugState{PhysPos{U},VecBelMvNormal{T}}) where T <: Real where U <: Real
    return zeros(length(vec(s.p)))
end
function termCost_grad_b(model::CostPosRangeLocalization,s::AugState{PhysPos{U},VecBelMvNormal{T}}) where T <: Real where U <: Real
    gVecs = Vector{Vector{T}}(undef, length(s.b));
    entropies = Distributions.entropy.(Distributions.MvNormal.(s.b));
    for ii = 1:length(s.b)
        gVecs[ii] = [zeros(s.b[ii].dim);exp(entropies[ii])*vec(0.5*inv(Matrix(Distributions.MvNormal(s.b[ii]).Σ)))]  # Check coefficient!
        #gVecs[ii] = [zeros(s.b[ii].dim);vec(0.5*inv(full(Distributions.MvNormal(s.b[ii]).Σ)))];
    end
    return gVecs
end

struct CostManipulate2D <: CostModel
end
function instCost(model::CostManipulate2D,b::BelMvNormal,u::MControl2D,Cs::Matrix{Float64},Cu::Matrix{Float64})
    μ,Σ = b.params[1:11],reshape(b.params[12:end],11,11); # \theta is between 0 and 2*pi;
    μ_target = [0.,0.,pi,0.,0.,0.,0.,0.,0.,0.,0.];    # Changing target \theta to pi.
    return 1/2*tr(Cs*Σ) + 1/2*(μ - μ_target)'*Cs*(μ - μ_target) + 1/2*vec(u)'*Cu*vec(u);
end
function instCost_grad(model::CostManipulate2D,b::BelMvNormal,Cs::Matrix{Float64})
    μ = b.params[1:11];
    μ_target = [0.,0.,pi,0.,0.,0.,0.,0.,0.,0.,0.];    # Changing target \theta to pi.
    return [Cs*(μ - μ_target);1/2*vec(Cs)]
end
function termCost(model::CostManipulate2D,b::BelMvNormal,Cs::Matrix{Float64})
    μ,Σ = b.params[1:11],reshape(b.params[12:end],11,11);
    μ_target = [0.,0.,pi,0.,0.,0.,0.,0.,0.,0.,0.];    # Changing target \theta to pi.
    return 1/2*tr(Cs*Σ) + 1/2*(μ - μ_target)'*Cs*(μ - μ_target);
end
termCost_grad(model::CostManipulate2D,b::BelMvNormal,Cs::Matrix{Float64}) = instCost_grad(model,b,Cs);


text:
function relocation_table!(mod)
    i64 = LLVM.IntType(64; ctx=LLVM.context(mod))
    d = IdDict{Any, Tuple{String, LLVM.GlobalVariable}}()

    for func ∈ LLVM.functions(mod), bb ∈ LLVM.blocks(func), inst ∈ LLVM.instructions(bb)
        if isa(inst, LLVM.LoadInst) && occursin("inttoptr", string(inst))
            get_pointers!(d, mod, inst)
        elseif isa(inst, LLVM.StoreInst) && occursin("inttoptr", string(inst))
            @debug "Relocating StoreInst" inst
            get_pointers!(d, mod, inst)
        elseif inst isa LLVM.RetInst && occursin("inttoptr", string(inst))
            @debug "Relocating RetInst" inst LLVM.operands(inst)
            get_pointers!(d, mod, inst)
        elseif isa(inst, LLVM.BitCastInst) && occursin("inttoptr", string(inst))
            @debug "Relocating BitCastInst" inst LLVM.operands(inst)
            get_pointers!(d, mod, inst)
        elseif isa(inst, LLVM.CallInst)
            @debug "Relocating CallInst" inst LLVM.operands(inst)
            dest = LLVM.called_value(inst)
            if occursin("inttoptr", string(dest)) && length(LLVM.operands(dest)) > 0
                @debug "Relocating CallInst inttoptr" dest LLVM.operands(dest) LLVM.operands(inst)
                ptr_arg = first(LLVM.operands(dest))
                ptr_val = convert(Int, ptr_arg)
                ptr = Ptr{Cvoid}(ptr_val)

                frames = ccall(:jl_lookup_code_address, Any, (Ptr{Cvoid}, Cint,), ptr, 0)
                if length(frames) >= 1
                    fn, file, line, linfo, fromC, inlined = last(frames)
                    fn = string(fn)
                    if ptr == cglobal(:jl_alloc_array_1d)
                        fn = "jl_alloc_array_1d"
                    end
                    if ptr == cglobal(:jl_alloc_array_2d)
                        fn = "jl_alloc_array_2d"
                    end
                    if ptr == cglobal(:jl_alloc_array_3d)
                        fn = "jl_alloc_array_3d"
                    end
                    if ptr == cglobal(:jl_new_array)
                        fn = "jl_new_array"
                    end
                    if ptr == cglobal(:jl_array_copy)
                        fn = "jl_array_copy"
                    end
                    if ptr == cglobal(:jl_alloc_string)
                        fn = "jl_alloc_string"
                    end
                    if ptr == cglobal(:jl_in_threaded_region)
                        fn = "jl_in_threaded_region"
                    end
                    if ptr == cglobal(:jl_enter_threaded_region)
                        fn = "jl_enter_threaded_region"
                    end
                    if ptr == cglobal(:jl_exit_threaded_region)
                        fn = "jl_exit_threaded_region"
                    end
                    if ptr == cglobal(:jl_set_task_tid)
                        fn = "jl_set_task_tid"
                    end
                    if ptr == cglobal(:jl_new_task)
                        fn = "jl_new_task"
                    end
                    if ptr == cglobal(:malloc)
                        fn = "malloc"
                    end
                    if ptr == cglobal(:memmove)
                        fn = "memmove"
                    end
                    if ptr == cglobal(:jl_array_grow_beg)
                        fn = "jl_array_grow_beg"
                    end
                    if ptr == cglobal(:jl_array_grow_end)
                        fn = "jl_array_grow_end"
                    end
                    if ptr == cglobal(:jl_array_grow_at)
                        fn = "jl_array_grow_at"
                    end
                    if ptr == cglobal(:jl_array_del_beg)
                        fn = "jl_array_del_beg"
                    end
                    if ptr == cglobal(:jl_array_del_end)
                        fn = "jl_array_del_end"
                    end
                    if ptr == cglobal(:jl_array_del_at)
                        fn = "jl_array_del_at"
                    end
                    if ptr == cglobal(:jl_array_ptr)
                        fn = "jl_array_ptr"
                    end
                    if ptr == cglobal(:jl_value_ptr)
                        fn = "jl_value_ptr"
                    end
                    if ptr == cglobal(:jl_get_ptls_states)
                        fn = "jl_get_ptls_states"
                    end
                    if ptr == cglobal(:jl_gc_add_finalizer_th)
                        fn = "jl_gc_add_finalizer_th"
                    end
                    if ptr == cglobal(:jl_symbol_n)
                        fn = "jl_symbol_n"
                    end
                end

                if length(fn) > 1 && fromC
                    mod = LLVM.parent(LLVM.parent(LLVM.parent(inst)))
                    lfn = LLVM.API.LLVMGetNamedFunction(mod, fn)

                    if lfn == C_NULL
                        lfn = LLVM.API.LLVMAddFunction(mod, fn, LLVM.API.LLVMGetCalledFunctionType(inst))
                    else
                        lfn = LLVM.API.LLVMConstBitCast(lfn, LLVM.PointerType(LLVM.FunctionType(LLVM.API.LLVMGetCalledFunctionType(inst))))
                    end
                    LLVM.API.LLVMSetOperand(inst, LLVM.API.LLVMGetNumOperands(inst)-1, lfn)
                end
            end
            get_pointers!(d, mod, inst)
        end
    end
    IdDict{Any, String}(val => name for (val, (name, _)) ∈ d)
end

function get_pointers!(d, mod, inst)
    jl_t = (LLVM.StructType(LLVM.LLVMType[]; ctx=LLVM.context(mod)))
    for (i, arg) ∈ enumerate(LLVM.operands(inst))
        if occursin("inttoptr", string(arg)) && arg isa LLVM.ConstantExpr
            op1 = LLVM.Value(LLVM.API.LLVMGetOperand(arg, 0))
            if op1 isa LLVM.ConstantExpr
                op1 = LLVM.Value(LLVM.API.LLVMGetOperand(op1, 0))
            end
            ptr = Ptr{Cvoid}(convert(Int, op1))
            frames = ccall(:jl_lookup_code_address, Any, (Ptr{Cvoid}, Cint,), ptr, 0)
            if length(frames) >= 1
                fn, file, line, linfo, fromC, inlined = last(frames)
                if (isempty(String(fn)) && isempty(String(file))) || fn == :jl_system_image_data
                    val = unsafe_pointer_to_objref(ptr)
                    if val ∈ keys(d)
                        _, gv = d[val]
                        LLVM.API.LLVMSetOperand(inst, i-1, gv)
                    else
                        gv_name = GPUCompiler.safe_name(String(gensym(repr(Core.Typeof(val)))))
                        gv = LLVM.GlobalVariable(mod, llvmeltype(arg), gv_name, LLVM.addrspace(llvmtype(arg)))

                        LLVM.extinit!(gv, true)
                        LLVM.API.LLVMSetOperand(inst, i-1, gv)

                        d[val] = (gv_name, gv)
                    end
                else
                    @warn "Found data we don't know how to relocate." frames
                end
            end
        end
    end
end

llvmeltype(x::LLVM.Value) = eltype(LLVM.llvmtype(x))

function pointer_patching_diff(mod::LLVM.Module, path1=tempname(), path2=tempname(); show_reloc_table=false)
    s1 = string(mod)
    write(path1, s1)

    d = StaticCompiler.relocation_table!(mod)
    if show_reloc_table
        @show d
    end

    s2 = string(mod)
    write(path2, s2)

    try
        # this always ends in an error for me for some reason
        run(`diff $path1 $path2`)
    catch e;
        nothing
    end
end


text:
println("Loading Libraries...")
include("post_process_defs_pyplot.jl")

# RESULTS_DIR = "results/vb-sim-paper/"
RESULTS_DIR = ARGS[1]

# Collecting paths to output
paths_to_output = String[]
for (root, dirs, files) in walkdir(RESULTS_DIR)
  if "output.bson" in files
    append!(paths_to_output, ["$root/output.bson"])
  end
end

# Post processing
println("Doing post processing...")
for path_to_output in paths_to_output
  println("Processing: $(path_to_output)")
  post_process(path_to_output)
end


text:
module CollisionDetection

using Compat

export Octree
export boxes, fitsinbox, boudingbox, boxesoverlap
export searchtree

include("octree.jl")
include("searcheq.jl")

end # module

max_stars_repo_name: shakisparki/CollisionDetection

max_issues_repo_name: shakisparki/CollisionDetection

max_forks_repo_name: shakisparki/CollisionDetection


text:
const A = 1.0
const B = 1.5
const B1 = B + 1

@taylorize function brusselator!(du, u, p, t)
    x, y = u
    x² = x * x
    aux = x² * y
    du[1] = A + aux - B1*x
    du[2] = B*x - aux
    return du
end

U0(r) = Singleton([1.0, 1.0]) ⊕ BallInf(zeros(2), r)

bruss(r) = @ivp(u' = brusselator!(u), u(0) ∈ U0(r), dim: 2)

sol_1 = solve(bruss(0.1), T=30.0, alg=TMJets20(orderT=6, orderQ=2))


text:
### Domain [mm]
x_max = 1

### Models ###
# The considered PSF for the considered imaging tecnique
	# sigma also defines the resolution of the imaging system.
	sigma = 0.04
	filter = (x, y) -> exp(-(x^2 + y^2)/2/sigma^2)
	filter_dx = (x, y) -> -x/sigma^2*filter(x, y)
	filter_dy = (x, y) -> -y/sigma^2*filter(x, y)
# The number of time samples for reconstruction: 2K+1
	K = 2
# the maximum allowed speed for any considered particle
	v_max = 15 #[mm/s]
        particle_velocity = 1 #[mm/s]
# The sampling rate
	tau = 1/500

### Noise ###
sigma_noise = 0.01

### Medium ###
# total number of the experimets measurements (gives the length of the experience)
	n_im = 500
# Vessels separation
	dx = 0.02
# microbubbles distance to booundary
	bndry_sep = 0.05

### Inverse problem ###
# Considered norm to minimize
	lp_norm = 1
# Threshold on when to consider that there was a difference on the number of particles
# at two consecutive time steps.
	jump_threshold = 0.1

### Sequence ###
# Probability of a particle to activate or deactivate
	activation_probability=0.02

# Building the vessel network
x0 = [0.1,0.1]
# arrays of t_max, relative direction, curve_type, parameters
curve1 = [0.3, -pi/6,1, "straight", [] ]
curve2 = [0.4, 0,1,"ellipse", [1,2,0.5]]
curve3 = [0.55, 0,1, "ellipse", [-1,2,1]]
curve4 = [0.45, pi/2,1, "ellipse", [-1,2,1]]
curve5 = [0.35, 0,1, "ellipse", [1,2,1]]

Vessel = Tree(curve1)      #curve 1
addchild(Vessel,1, curve2) #curve 2
addchild(Vessel,1, curve3) #curve 3
addchild(Vessel,2, curve4) #curve 4
addchild(Vessel,2, curve5) #curve 5

# The tree is transform to obtain the curves precisely,
# i.e. not defined relative to the neighbouring ones.
Absolute_Vessel(Vessel, 1, x0, 0)


text:
using .BifurcationsBase: BifurcationProblem, BifurcationSolver
using .Codim1: AbstractSpecialPoint
using .Continuations: ContinuationOptions
import .Continuations: init

init(prob::BifurcationProblem; kwargs...) =
    BifurcationSolver(prob, ContinuationOptions(; kwargs...))

init(point::AbstractSpecialPoint, args...; kwargs...) =
    init(BifurcationProblem(point, args...); kwargs...)


text:
using Dates; now()
using Revise, MathematicalSystems, Reachability, LinearAlgebra, HybridSystems
using LazySets
using Reachability: solve
using Plots
using TaylorIntegration

@taylorize function bball_up!(t, x, dx)
    dx[1] = x[2]
    dx[2] = -9.8 - 0.1*(x[1])^2
    return dx
end

@taylorize function bball_down!(t, x, dx)
    dx[1] = x[2]
    dx[2] = -9.8 + 0.1*(x[1])^2
    return dx
end

function bouncing_ball()

    automaton = LightAutomaton(2) # two modes

    inv_up = HPolyhedron([HalfSpace([-1.0, 0.0], 0.0),  # x >= 0
                          HalfSpace([0.0, -1.0], 0.0)]) # v >= 0

    inv_down = HPolyhedron([HalfSpace([-1.0, 0.0], 0.0),  # x >= 0
                         HalfSpace([ 0.0, 1.0], 0.0)])  # v <= 0

    m1 = ConstrainedBlackBoxContinuousSystem(bball_up!, 2, inv_up)

    m2 = ConstrainedBlackBoxContinuousSystem(bball_down!, 2, inv_down)

    modes = [m1,m2]  #modes

    add_transition!(automaton, 2, 1, 1)    #alpha transition
    add_transition!(automaton, 1, 2, 2)    #beta transition

    guard_alpha = HPolyhedron([HalfSpace([1.0, 0.0], 0.0),   #x>=0
                              HalfSpace([-1.0, 0.0], 0.0)])  #x<=0
    guard_beta =  HPolyhedron([HalfSpace([0.0, 1.0], 0.0),   #v<=0
                               HalfSpace([0.0,-1.0], 0.0)])  #v>=0

    A = [1.0 0.0; 0.0 -0.8]
    t1 = ConstrainedLinearMap(A, guard_alpha)
    t2 = ConstrainedLinearMap(A, guard_beta)

    #resetmaps
    resetmaps = [t1,t2]

    # switching
    switching = AutonomousSwitching()
    switchings = fill(switching, 2)

    ℋ = HybridSystem(automaton, modes, resetmaps, switchings)

    #initial condition in mode one
    X0= Hyperrectangle(low=[4.9, -0.2], high=[5.1, 0.0])

    initial_condition = [(2, X0)]  # initial condition in the "down" mode

    system = InitialValueProblem(ℋ, initial_condition)

    options = Options(:mode=>"reach", :T=>6.0, :plot_vars=>[1, 2], :project_reachset=>false)

    return (system, options)
end

problem, options = bouncing_ball();

options = Options(:mode=>"reach", :T=>3.0, :plot_vars=>[1, 2], :project_reachset=>false, :verbosity => "info")

@time sol_TMJets = solve(problem, options, TMJets(:orderT=>5, :orderQ=>2, :abs_tol=>1e-10),LazyDiscretePost(:check_invariant_intersection=>true));

plot(sol_TMJets, use_subindices=false, aspectratio=1, alpha=.5, color=:lightblue)

max_stars_repo_path: models/bouncing ball with air friction/Hybrid_model_bouncing_ball.jl

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
function generate_namedtuple(::Type{NamedTuple{names,types}}, q) where {names,types}
    if @generated
        vals = Tuple(:(getvalue(q, $i, $(fieldtype(types, i)))) for i = 1:fieldcount(types))
        return :(NamedTuple{names,types}(($(vals...),)))
    else
        return NamedTuple{names,types}(Tuple(getvalue(q, i, fieldtype(types, i)) for i = 1:fieldcount(types)))
    end
end

# This section of the code hijacks SQLite internals to make queryverse-compatible iterators
# TODO: submit as a PR to SQLite
struct SQLiteCursor{Row}
    statement::Stmt
    status::RefValue{Cint}
    cursor_row::RefValue{Int}
end

function eltype(::SQLiteCursor{Row}) where {Row}
    Row
end
function IteratorSize(::Type{<:SQLiteCursor})
    SizeUnknown()
end

function isdone(cursor::SQLiteCursor)
    status = cursor.status[]
    if status == SQLITE_DONE
        true
    elseif status == SQLITE_ROW
        false
    elseif sqliteerror(cursor.statement.db)
        false
    else
        error("Unknown SQLite cursor status")
    end
end

function getvalue(cursor::SQLiteCursor, column_number::Int, ::Type{Value}) where {Value}
    handle = cursor.statement.handle
    column_type = sqlite3_column_type(handle, column_number)
    if column_type == SQLITE_NULL
        Value()
    else
        julia_type = juliatype(column_type) # native SQLite Int, Float, and Text types
        sqlitevalue(
            if julia_type === Any
            if !isbitstype(Value)
                Value
            else
                julia_type
            end
        else
            julia_type
        end, handle, column_number)
    end
end

function iterate(cursor::SQLiteCursor{Row}) where {Row}
    if isdone(cursor)
        nothing
    else
        named_tuple = generate_namedtuple(Row, cursor)
        cursor.cursor_row[] = 1
        named_tuple, 1
    end
end

function iterate(cursor::SQLiteCursor{Row}, state) where {Row}
    if state != cursor.cursor_row[]
        error("State does not match SQLiteCursor cursor_row")
    else
        cursor.status[] = sqlite3_step(cursor.statement.handle)
        if isdone(cursor)
            nothing
        else
            named_tuple = generate_namedtuple(Row, cursor)
            cursor.cursor_row[] = state + 1
            named_tuple, state + 1
        end
    end
end

function isiterable(::SourceCode)
    true
end
function isiterabletable(::SourceCode)
    true
end
function collect(source::SourceCode)
    collect(getiterator(source))
end

function second((value_1, value_2))
    value_2
end

function name_and_type(handle, column_number, nullable=true, strict_types=true)
    Symbol(unsafe_string(sqlite3_column_name(handle, column_number))),
    if strict_types
        julia_type = juliatype(handle, column_number)
        if nullable
            DataValue{julia_type}
        else
            julia_type
        end
    else
        Any
    end
end

function getiterator(source_code::SourceCode)
    statement = Stmt(
        source_code.source,
        string(finalize(translate(source_code.code)))
    )
    # bind!(statement, values)
    status = execute(statement)
    handle = statement.handle
    schema = ntuple(
        let handle = handle
        column_number -> name_and_type(handle, column_number)
    end,
        sqlite3_column_count(handle)
    )
    SQLiteCursor{NamedTuple{Tuple(map(first, schema)),Tuple{map(second, schema)...}}}(statement, Ref(status), Ref(0))
end

# Use default show methods from the queryverse
function show(stream::IO, source::SourceCode)
    printtable(stream, getiterator(source), "SQLite query result")
end

function showable(::MIME"text/html", source::SourceCode)
    true
end
function show(stream::IO, ::MIME"text/html", source::SourceCode)
    printHTMLtable(stream, getiterator(source))
end

function showable(::MIME"application/vnd.dataresource+json", source::SourceCode)
    true
end
function show(stream::IO, ::MIME"application/vnd.dataresource+json", source::SourceCode)
    printdataresource(stream, getiterator(source))
end


text:
using SparseArrays: spzeros
import JuMP
const MOI = JuMP.MOI

"""
    NonNegKSparseLMO{T}(K::Int, right_hand_side::T)
LMO for the non-negative K-sparse polytope:
```
C = B_1(τK) ∩ B_∞(τ) ∩ R^n_+
```
with `τ` the `right_hand_side` parameter.
The LMO results in a vector with the K smallest negative values
of direction, taking values `τ`.
"""
struct NonNegKSparseLMO{T} <: FrankWolfe.LinearMinimizationOracle
    K::Int
    right_hand_side::T
end

function FrankWolfe.compute_extreme_point(lmo::NonNegKSparseLMO{T}, direction; kwargs...) where {T}
    K = min(lmo.K, length(direction))
    K_indices = sortperm(direction[1:K])
    K_values = direction[K_indices]
    for idx in K+1:length(direction)
        new_val = direction[idx]
        # new smaller value: shift everything right
        if new_val < K_values[1]
            K_values[2:end] .= K_values[1:end-1]
            K_indices[2:end] .= K_indices[1:end-1]
            K_indices[1] = idx
            K_values[1] = new_val
            # new value in the interior
        elseif new_val < K_values[K]
            # NOTE: not out of bound since unreachable with K=1
            j = K - 1
            while new_val < K_values[j]
                j -= 1
            end
            K_values[j+1:end] .= K_values[j:end-1]
            K_indices[j+1:end] .= K_indices[j:end-1]
            K_values[j+1] = new_val
            K_indices[j+1] = idx
        end
    end
    v = spzeros(T, length(direction))
    for (idx, val) in zip(K_indices, K_values)
        if val < 0
            v[idx] = lmo.right_hand_side
        end
    end
    return v
end

function FrankWolfe.convert_mathopt(
    lmo::NonNegKSparseLMO{T},
    optimizer::OT;
    dimension::Integer,
    kwargs...,
) where {T,OT}
    τ = lmo.right_hand_side
    n = dimension
    K = min(lmo.K, n)
    MOI.empty!(optimizer)
    (x, _) = MOI.add_constrained_variables(optimizer, fill(MOI.Interval(0.0, τ), n))
    MOI.add_constraint(optimizer, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(n), x), 0.0), MOI.LessThan(τ * K))
    return FrankWolfe.MathOptLMO(optimizer)
end

max_stars_repo_name: jmaces/fw-rde

max_issues_repo_name: ZIB-IOL/fw-rde

max_forks_repo_name: ZIB-IOL/fw-rde


text:
## karhunen_loeve.jl : approximate Gaussian random field generator using a Karhunen-Lo\`eve decomposition

"""
    KarhunenLoeve{n} <: GaussianRandomFieldGenerator

A [`GaussiandRandomFieldGenerator`](@ref) using a Karhunen-Lo\u00e8ve (KL) expansion with `n` terms.

# Examples
```jldoctest
julia> m = Matern(0.1,1.0)
Matérn (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> c = CovarianceFunction(2,m)
2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> pts1 = 0:0.02:1; pts2 = 0:0.02:1
0.0:0.02:1.0

julia> grf = GaussianRandomField(c,KarhunenLoeve(300),pts1,pts2)
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51x51 structured grid, using a KL expansion with 300 terms

julia> plot_eigenvalues(grf) # plot the eigenvalue decay
[...]

julia> plot_eigenfunction(grf,4) # plots the fourth eigenfunction
[...]

```
The more terms are retained in the expansion, the better the approximation will be.
```jldoctest
julia> nterms = [1 2 5 10 20 50 100 200 500 1000]
1×10 Array{Int64,2}:
1  2  5  10  20  50  100  200  500  1000

julia> for n in nterms
grf = GaussianRandomField(c,KarhunenLoeve(n),pts1,pts2)
@show rel_error(grf)
end
rel_error(grf) = 0.7499982529722711
rel_error(grf) = 0.49999825591379987
rel_error(grf) = 0.4425751861338164
rel_error(grf) = 0.35789475278408045
rel_error(grf) = 0.16805079842673853
rel_error(grf) = 0.11187098338277579
rel_error(grf) = 0.05130466343704787
rel_error(grf) = 0.017343327476498027
rel_error(grf) = 0.0034278579378175245
rel_error(grf) = 0.0007216777927243623

```
The KL expansion is computed using the Nystrom method. Optional argument are the quadrature rule `quad` and the number of quadrature points in each direction `nq`. Possible values for `quad` are `GaussLegendre()`, `Midpoint()`, `Simpson()`, `Trapezoidal` and `EOLE()` (default). The total number of quadrature points must be larger than or equal to the requested number of terms. Default is `ceil(n^(1/d))`. This value is best left untouched.

```jldoctest
julia> grf = GaussianRandomField(c,KarhunenLoeve(300),pts1,pts2,quad=GaussLegendre())
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51x51 structured grid, using a KL expansion with 300 terms

julia> grf = GaussianRandomField(c,KarhunenLoeve(300),pts1,pts2,quad=Simpson(),nq=20)
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51x51 structured grid, using a KL expansion with 300 terms

```
See also: [`Cholesky`](@ref), [`Spectral`](@ref), [`CirculantEmbedding`](@ref)
"""
struct KarhunenLoeve{n} <: GaussianRandomFieldGenerator
    function KarhunenLoeve{n}() where n
        n > 0 || throw(DomainError(n, "in KarhunenLoeve{n}, number of terms n must be positive"))
        new{n}()
    end
end

KarhunenLoeve(n::Integer) = KarhunenLoeve{n}()

function _GaussianRandomField(mean, cov::CovarianceFunction{d}, method::KarhunenLoeve{n},
                              pts...;
                              nq::Integer=ceil(typeof(n), n^(1/d)),
                              quad::QuadratureRule=EOLE(),
                              eigensolver::AbstractEigenSolver=EigsSolver()) where {d,n}
    # check if number of terms and number of quadrature points are compatible
    nq = nq > 0 ? ntuple(i -> nq, d) : length.(pts)
    prod_nq = prod(nq)
    prod_nq < n && throw(ArgumentError("too many terms requested, increase nq or lower n"))
    # adjustment for ARPACK error when looking for ALL eigenvalues
    eigensolver isa EigsSolver && prod_nq == n && (nq = nq .+ 1)

    # determine bounding box when irregular mesh
    is_irr = pts[1] isa AbstractMatrix
    a = is_irr ? Tuple(minimum(pts[1],dims=2)) : minimum.(pts)
    b = is_irr ? Tuple(maximum(pts[1],dims=2)) : maximum.(pts)

    # compute quadrature nodes and weights
    struc = get_nodes_and_weights.(nq, a, b, Ref(quad))
    nodes = first.(struc)
    weights = last.(struc)

    # eigenvalue problem
    C = apply(cov,nodes,nodes)
    W = d == 1 ? Diagonal(weights...) : Diagonal(kron(weights...))
    W .= sqrt.(W)
    B = Symmetric(W * C * W) # should be symmetric and positive semi-definite
    isposdef(B) || @warn "equivalent eigenvalue problem is not SPD, results may be wrong or inaccurate"

    # solve
    eigenval, eigenfunc = compute(B, n, eigensolver)

    # compute eigenfunctions in nodes
    K = apply(cov, pts, nodes)
    Λ = Diagonal(1 ./ eigenval)
    eigenfunc = K * W * eigenfunc * Λ

    m = findfirst(x -> x < 0, eigenval)
    if m != nothing
        m -= 1
        @warn begin
            "$(length(eigenval) - m) negative eigenvalues ≥ $(eigenval[end]) detected, "
            "Gaussian random field will be approximated (ignoring all negative eigenvalues)"
        end
    else
        m = n
    end

    resize!(eigenval, m)
    eigenfunc = eigenfunc[:, 1:m]

    # store eigenvalues and eigenfunctions
    eigenval .= sqrt.(eigenval) # note: store sqrt of eigenval for more efficient sampling
    data = SpectralData(eigenval, eigenfunc)

    GaussianRandomField{KarhunenLoeve{length(eigenval)},typeof(cov),typeof(pts),typeof(mean),typeof(data)}(mean,cov,pts,data)
end

# returns the required dimension of the random points
randdim(grf::GaussianRandomField{KarhunenLoeve{n}}) where n = n

# relative error in the KL approximation
function rel_error(grf::GaussianRandomField{<:KarhunenLoeve})
    Leb = prod(grf.pts) do point
        a, b = extrema(point)
        b - a
    end
    1 - sum(abs2, grf.data.eigenval) / (std(grf.cov)^2 * Leb)
end

# sample function for both Spectral() and KarhunenLoeve(n) type
function _sample(grf::GaussianRandomField{<:Union{Spectral,KarhunenLoeve}}, xi)
    grf.mean + std(grf.cov) * reshape((grf.data.eigenfunc * Diagonal(grf.data.eigenval)) * xi, size(grf.mean))
end

Base.show(io::IO, ::KarhunenLoeve{n}) where n = print(io, "KL expansion with $(n) terms")


text:
#!/usr/bin/env julia

using Luxor, Test

using Random
Random.seed!(42)

function general_tests()
    pt1 = Point(rand() * 4, rand() * 4)

    # arithmetic with tuples
    @test pt1 + (5, 8) == Point(pt1.x + 5, pt1.y + 8)
    @test pt1 * (5, 8) == Point(pt1.x * 5, pt1.y * 8)
    @test pt1 - (5, 8) == Point(pt1.x - 5, pt1.y - 8)
    @test pt1 / (5, 8) == Point(pt1.x / 5, pt1.y / 8)

    @test -pt1 == Point(-pt1.x, -pt1.y)

    # is point/4 inside a box
    # a: we now have to wrap arguments with Ref() to ensure they  broadcast as scalar
    @test isinside(Ref(pt1) ./ 4, box(O, 10, 10, vertices=true))

    # b: but we defined Point as broadcastable, so we shouldn't need that
    @test isinside(pt1 ./ 4, box(O, 10, 10, vertices=true))

    # is point not in every corner of box
    @test all(Ref(Point(1, 1)) .< box(O, 10, 10, vertices=true)) == false
    @test any(Ref(Point(0, 0)) .< [Point(1, 1), Point(1, 2), Point(2, 3)]) == true
    # is point outside every corner of box
    @test all(Ref(Point(10, 10))  .> box(O, 10, 10, vertices=true)) == true

    # ? do these need Ref()?
    @test .>=(Point(5, 5),  Point(5, 5))
    @test .>=(Point(5, 5),  Point(5, 5))

    @test perpendicular(Point(10, 5)) == Point(5.0, -10.0)
    @test perpendicular(Point(-10, -5)) == Point(-5.0,10.0)

    @test crossproduct(Point(2, 3), Point(3, 2)) == -5.0
    @test crossproduct(Point(-20, 30), Point(60, 20)) < 2000

    # shared end points
    # intersection of (A == C) || (B == C) || (A == D) || (B == D)
    pt1 = Point(5, 5)
    pt2 = Point(6, 5)
    @test pt1 + pt2 == Point(11.0,10.0)
    @test intersectionlines(pt1, pt1, pt1, pt1)[1] == false
    @test intersectionlines(pt1, pt2, pt1, pt2, crossingonly=true)[1] == false
    @test intersectionlines(pt1, pt2, pt2, pt1)[1] == false
    pt3 = Point(6, 6)
    pt4 = Point(5, 6)
    @test intersectionlines(pt1, pt2, pt3, pt4)[1] == false

    @test intersectionlines(pt1, pt2, pt3, pt2) == (true,Luxor.Point(6.0,5.0))

    # not crossing
    pt1 = Point(5, 5)
    pt2 = Point(5, 5)
    pt3 = Point(5, 5)
    pt4 = Point(5, 5)
    @test intersectionlines(shuffle([pt1, pt2, pt3, pt4])...)[1] == false

    # parallel
    pt1 = Point(5, 5)
    pt2 = Point(5, 6)
    pt3 = Point(6, 5)
    pt4 = Point(6, 6)
    @test intersectionlines(pt1, pt2, pt3, pt4)[1] == false

    @test pt1 * pt2 == Luxor.Point(25.0, 30.0)
    @test pt1 ^ 2 == Luxor.Point(25.0, 25.0)
    @test pt1 ^ 3 == Luxor.Point(125.0,125.0)

    @test pt1 ^ 1.5 < pt1 ^ 1.6
    @test pt1 ^ 1.6 > pt1 ^ 1.5

    # test between interpolation
    @test isequal(midpoint(pt1, pt3), between(pt1, pt3, 0.5))
    @test isequal(between(pt1, pt3, 0.0), pt1)
    @test isequal(between(pt1, pt3, 1.0), pt3)

    # test uniqueness and those pesky -0.0s
    @test length(unique([pt1, pt2, pt3, pt4])) == 4
    @test length(unique([pt1, pt2, pt3, pt3])) == 3
    @test length(unique([pt1, pt3, pt3, pt3])) == 2
    @test length(unique([pt3, pt3, pt3, pt3])) == 1
    @test length(unique([])) == 0
    @test length(unique([Point(0.0, 0.0), Point(0.0, 0.0)])) == 1
    @test length(unique([Point(0.0, -0.0), Point(0.0, 0.0)])) == 1
    @test length(unique([Point(0.0, -0.0), Point(-0.0, 0.0)])) == 1
    @test length(unique([Point(0.0, -0.0), Point(-0.0, -0.0)])) == 1
    @test length(unique([Point(0.0, 0.0), Point(0, 0)])) == 1
    @test length(unique([Point(0, 0), Point(-0, -0)])) == 1

    # test drop perpendicular
    p1 = Point(0, 0)
    p2 = Point(0, 50)
    @test perpendicular(p1, p2, 10) == Luxor.Point(-10.0, 0.0)
    @test perpendicular(p1, p2, 20) == Luxor.Point(-20.0, 0.0)

    # test pointinverse
    p1 = Point(10, 10)
    cpt = O
    rad = 50
    flag, antipoint = pointinverse(p1, cpt, rad)
    @test flag == true
    @test isapprox(antipoint, Point(125, 125))

    p2 = Point(20, 0.0)
    flag, antipoint = pointinverse(p2, cpt, rad)
    @test flag == true
    @test isapprox(antipoint.y, 0.0)

    # point shouldn't be centerpoint
    @test_throws ErrorException pointinverse(cpt, cpt, rad)
end

function point_arithmetic_test(fname, npoints=20)
    Drawing(1200, 1200, fname)
    background("white")
    sethue("thistle")
    box(BoundingBox(), :stroke)
    scale(0.5, 0.5)
    setline(2.5)
    setopacity(0.5)
    fontsize(8)
    origin()
    randompoints = randompointarray(Point(-600, -600), Point(600, 600), npoints)

    # +
    pl1 = map(pt -> pt + 2, randompoints)
    pl1a = map(pt -> 2 + pt, randompoints)

    # -
    pl2 = map(pt -> pt - 2, randompoints)
    pl2a = map(pt -> 2 - pt, randompoints)

    # *
    pl3 = map(pt -> pt * rand(), randompoints)
    pl3a = map(pt -> rand() * pt, randompoints)

    # /
    pl4 = map(pt -> pt / rand(), randompoints)

    # .*
    pl5 = 1.012 .* randompoints
    pl5a = randompoints .* 1.012
    pl5b = [1.03, 0.97, -1.05] .* randompoints[1:3]

    # ./
    pl6 = randompoints[1:3] ./ [1.03, 0.97, -1.05]
    pl6a =  [1.03, 0.97, -1.05] .* randompoints[1:3]

    # ^
    pl7 = map(pt -> ^(pt/2, 2), randompoints)

    # ^ FAILS
    # pl8 = map(pt -> ^(pt, 1.2), randompoints)

    # issue #20615 now fixed in v0.6
    sethue("red")
    for p in zip(pl1, pl1a, pl2, pl2a, pl3, pl3a, pl4, pl5, pl5a, pl5b, pl6, pl6a, pl7)
       map(pt -> circle(pt, 3, :fill), p)
       poly(collect(p), :stroke)
    end

    # comparisons
    testfunctions = [isequal, isless, <, >, ==]
    sethue("green")
    for f in testfunctions
        for i in 1:length(randompoints)-1
            if f(randompoints[i], randompoints[i+1])
                circle(randompoints[i], 18, :fill)
                circle(randompoints[i+1], 7, :fill)
                text(string("$f"), randompoints[i] - 12)
            end
        end
    end

    testfunctions = [<, >, >=, <=]
    sethue("orange")
    for f in testfunctions
        for i in 1:length(randompoints)-1
            if f.(randompoints[1:npoints], randompoints[npoints:-1:1]) != true
                ellipse(randompoints[i], 19, 32, :fill)
                text(string("v6"), randompoints[i] - 6)
            end
        end
    end

    sethue("purple")
    for i in 1:length(randompoints)-1
        v = cmp(randompoints[i], randompoints[i+1])
        text(string(round(v, digits=1)), randompoints[i] + 15)
    end

    sethue("cyan")
    for i in 1:length(randompoints)-1
        n = distance(randompoints[i], randompoints[i+1])
        text(string(round(n, digits=1)), randompoints[i])
    end

    sethue("magenta")
    map(pt -> circle(pt, 6, :stroke), [midpoint(randompoints), randompoints[1], randompoints[2]])

    if all(randompoints .== randompoints)
        text("the points compare elementwise")
    else
        error("elementwise comparison failed")
    end

    if any(randompoints .!= randompoints)
        error("elementwise comparison failed")
    else
        text("the points really do compare elementwise", O + 15)
    end

    @test finish() == true

    println("...finished test: output in $(fname)")
end

general_tests()
point_arithmetic_test("point-arithmetic.pdf", 100)


text:
using CJKFrequencies

using Test

@testset "all tests" begin
    include("charfreq_tests.jl")
    include("freq_datasets_tests.jl")
    include("lexicon_tests.jl")
end


text:
using StatsBase.sample
using TermWin

include("utils.jl")
include("board.jl")
include("player.jl")
include("draw.jl")
include("game.jl")

max_stars_repo_name: RatanRSur/robots-from-scratch

max_issues_repo_name: RatanRSur/robots-from-scratch

max_forks_repo_name: RatanRSur/robots-from-scratch


text:
using KitBase, FluxReconstruction, OrdinaryDiffEq

begin
    x0 = -1
    x1 = 1
    ncell = 100
    nface = ncell + 1
    dx = (x1 - x0) / ncell
    deg = 2 # polynomial degree
    nsp = deg + 1
    cfl = 0.1
    dt = cfl * dx
    t = 0.0
    a = 1.0
end

pspace = FRPSpace1D(x0, x1, ncell, deg)

begin
    xFace = collect(x0:dx:x1)
    xGauss = legendre_point(deg)
    xsp = global_sp(xFace, xGauss)
    ll = lagrange_point(xGauss, -1.0)
    lr = lagrange_point(xGauss, 1.0)
    lpdm = ∂lagrange(xGauss)
end

u = zeros(ncell, nsp)
for i = 1:ncell, ppp1 = 1:nsp
    u[i, ppp1] = exp(-20.0 * xsp[i, ppp1]^2)
end
uold = deepcopy(u)
f = zeros(ncell, nsp)
u_face = zeros(ncell, 2)
f_face = zeros(ncell, 2)
au = zeros(nface);
f_interaction = zeros(nface)

e2f = zeros(Int, ncell, 2)
for i = 1:ncell
    if i == 1
        e2f[i, 2] = nface
        e2f[i, 1] = i + 1
    elseif i == ncell
        e2f[i, 2] = i
        e2f[i, 1] = 1
    else
        e2f[i, 2] = i
        e2f[i, 1] = i + 1
    end
end
f2e = zeros(Int, nface, 2)
for i = 1:nface
    if i == 1
        f2e[i, 1] = i
        f2e[i, 2] = ncell
    elseif i == nface
        f2e[i, 1] = 1
        f2e[i, 2] = i - 1
    else
        f2e[i, 1] = i
        f2e[i, 2] = i - 1
    end
end

function mol!(du, u, p, t) # method of lines
    xFace, e2f, f2e, a, deg, ll, lr, lpdm = p

    ncell = size(u, 1)
    nsp = size(u, 2)

    f = zeros(ncell, nsp)
    for i = 1:ncell, j = 1:nsp
        J = (xFace[i+1] - xFace[i]) / 2.0
        f[i, j] = advection_flux(u[i, j], a) / J
    end

    u_face = zeros(ncell, nsp)
    f_face = zeros(ncell, nsp)

    u_face[:, 1] .= u * lr
    f_face[:, 1] .= f * lr
    u_face[:, 2] .= u * ll
    f_face[:, 2] .= f * ll
    #=for i = 1:ncell, j = 1:nsp
        # right face of element i
        u_face[i, 1] += u[i, j] * lr[j]
        f_face[i, 1] += f[i, j] * lr[j]

        # left face of element i
        u_face[i, 2] += u[i, j] * ll[j]
        f_face[i, 2] += f[i, j] * ll[j]
    end=#

    au = zeros(nface)
    for i = 1:nface
        au[i] =
            (f_face[f2e[i, 1], 2] - f_face[f2e[i, 2], 1]) /
            (u_face[f2e[i, 1], 2] - u_face[f2e[i, 2], 1] + 1e-6)
    end

    f_interaction = zeros(nface)
    for i = 1:nface
        f_interaction[i] = (
            0.5 * (f_face[f2e[i, 2], 1] + f_face[f2e[i, 1], 2]) -
            0.5 * abs(au[i]) * (u_face[f2e[i, 1], 2] - u_face[f2e[i, 2], 1])
        )
    end

    dgl, dgr = ∂radau(deg, xGauss)

    rhs1 = zeros(ncell, nsp)
    for i = 1:ncell, ppp1 = 1:nsp, k = 1:nsp
        rhs1[i, ppp1] += f[i, k] * lpdm[ppp1, k]
    end

    for i = 1:ncell, ppp1 = 1:nsp
        du[i, ppp1] =
            -(
                rhs1[i, ppp1] +
                (f_interaction[e2f[i, 2]] - f_face[i, 2]) * dgl[ppp1] +
                (f_interaction[e2f[i, 1]] - f_face[i, 1]) * dgr[ppp1]
            )
    end
end

tspan = (0.0, 2.0)
p = (xFace, e2f, f2e, a, deg, ll, lr, lpdm)
prob = ODEProblem(mol!, u, tspan, p)
sol = solve(prob, Tsit5(), reltol = 1e-8, abstol = 1e-8, progress = true)

using Plots
plot(xsp[:, 2], sol.u[end][:, 2], label = "t=2")
scatter!(xsp[:, 2], exp.(-20 .* xsp[:, 2] .^ 2), label = "t=0")


text:
using HorizonSideRobots

function corner(r::Robot)
    s=0
    w=0
    while !isborder(r,Sud)
        move!(r,Sud)
        s+=1
    end
    while !isborder(r,West)
        move!(r,West)
        w+=1
    end
    diag(r)
    move!(r,West)
    putmarker!(r)
    paint(r)
    while !isborder(r,Sud)
        move!(r,Sud)
    end
    while !isborder(r,West)
        move!(r,West)
    end
    for a in 1:w
        move!(r,Ost)
    end
    for a in 1:s
        move!(r,Nord)
    end
end

function paint(r::Robot)
    while !isborder(r,Sud) || !isborder(r,West)
        while !ismarker(r)
            putmarker!(r)
            move!(r,Ost)
        end
        move!(r,Sud)
        while !isborder(r,West)
            move!(r,West)
        end
    end
    while !isborder(r,Ost)
        putmarker!(r)
        move!(r,Ost)
    end
end

function diag(r)
    while !isborder(r,Ost)
        move!(r,Ost)
    end
    putmarker!(r)
    while !isborder(r,Nord) && !isborder(r,West)
        move!(r,Nord)
        move!(r,West)
        putmarker!(r)
    end
end

corner(r)

max_stars_repo_name: Cheburkova/mirea-progs

max_issues_repo_name: Cheburkova/mirea-progs

max_forks_repo_name: Cheburkova/mirea-progs


text:
using SpawnerRecruitModels
using Test

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


text:
@testset "param" begin

# Need to test symmetry of param3dmat for different boundary conditions.
@testset "symmetry of create_paramop, uniform grid" begin
    # Create a grid.
    isbloch = [true, true, true]
    # isbloch = [true, false, false]
    lprim = ([-1.5, -0.5, 0.5, 1.5], [-1.5, -0.5, 0.5, 1.5], [-1.5, -0.5, 0.5, 1.5])
    g3 = Grid(lprim, isbloch)
    N = g3.N

    # Create materials.
    εvac = 1.0
    vac = Material("Vacuum", ε=εvac)

    εdiel = 2.0
    diel = Material("Dielectric", ε=εdiel)

    # Create objects.
    dom_vac = Object(Cuboid(g3.bounds), vac)
    obj_diel = Object(Cuboid([0,0,0], [1,1,1]), diel)
    # obj_diel = Object(Ball([0,0,0], 1), diel)

    # Add objects.
    oind2shp = Shape3[]
    oind2εind = ParamInd[]
    oind2μind = ParamInd[]
    εind2ε = S²ComplexF3[]
    μind2μ = S²ComplexF3[]

    add_obj!(oind2shp, (oind2εind,oind2μind), (εind2ε,μind2μ), dom_vac, obj_diel)

    # Construct arguments and call assign_param!.
    N = g3.N
    ε3d = create_param_array(N)
    εxx_oind3d = create_oind_array(N)
    εyy_oind3d = create_oind_array(N)
    εzz_oind3d = create_oind_array(N)
    εoo_oind3d = create_oind_array(N)

    μ3d = create_param_array(N)
    μxx_oind3d = create_oind_array(N)
    μyy_oind3d = create_oind_array(N)
    μzz_oind3d = create_oind_array(N)
    μoo_oind3d = create_oind_array(N)

    boundft = SVec(EE,EE,EE)
    assign_param!(ε3d, (μxx_oind3d,μyy_oind3d,μzz_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.ghosted.τl, g3.isbloch)
    assign_param!(ε3d, tuple(μoo_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.ghosted.τl, g3.isbloch)
    assign_param!(μ3d, (εxx_oind3d,εyy_oind3d,εzz_oind3d), oind2shp, oind2μind, μind2μ, ft2gt.(HH,boundft), g3.ghosted.τl, g3.isbloch)
    assign_param!(μ3d, tuple(εoo_oind3d), oind2shp, oind2μind, μind2μ, ft2gt.(HH,boundft), g3.ghosted.τl, g3.isbloch)

    # Perform smoothing.
    smooth_param!(ε3d, (εxx_oind3d,εyy_oind3d,εzz_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.l, g3.ghosted.l, g3.σ, g3.ghosted.∆τ)
    smooth_param!(ε3d, tuple(εoo_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.l, g3.ghosted.l, g3.σ, g3.ghosted.∆τ)

    ∆l = g3.∆l[nDL]
    ∆l′ = g3.∆l[nPR]
    ∆l′⁻¹ = map(x->inv.(x), ∆l′)

    Mε = create_paramop(ε3d, boundft.!=EE, ∆l, ∆l′⁻¹, g3.isbloch, order_cmpfirst=true)

    # @test issymmetric(Mε)
    @test Mε ≈ transpose(Mε)

    # Code I used to reveal the matrix structure while debugging:
    # (i,j) = (nX,nY); full(real.(Mε3))[27(i-1)+1:27i, 27(j-1)+1:27j]
end  # @testset "create_paramop"

@testset "symmetry of create_paramop, nonuniform grid" begin
    # Create a grid.
    isbloch = [true, true, true]
    # isbloch = [true, false, false]
    lprim = ([-1.5, -0.4, 0.3, 1.5], [-1.5, -0.4, 0.3, 1.5], [-1.5, -0.4, 0.3, 1.5])
    g3 = Grid(lprim, isbloch)
    N = g3.N

    # Create materials.
    εvac = 1.0
    vac = Material("Vacuum", ε=εvac)

    εdiel = 2.0
    diel = Material("Dielectric", ε=εdiel)

    # Create objects.
    dom_vac = Object(Cuboid(g3.bounds), vac)
    obj_diel = Object(Cuboid([0,0,0], [1,1,1]), diel)
    # obj_diel = Object(Ball([0,0,0], 1), diel)

    # Add objects.
    oind2shp = Shape3[]
    oind2εind = ParamInd[]
    oind2μind = ParamInd[]
    εind2ε = S²ComplexF3[]
    μind2μ = S²ComplexF3[]

    add_obj!(oind2shp, (oind2εind,oind2μind), (εind2ε,μind2μ), dom_vac, obj_diel)

    # Construct arguments and call assign_param!.
    N = g3.N
    ε3d = create_param_array(N)
    εxx_oind3d = create_oind_array(N)
    εyy_oind3d = create_oind_array(N)
    εzz_oind3d = create_oind_array(N)
    εoo_oind3d = create_oind_array(N)

    μ3d = create_param_array(N)
    μxx_oind3d = create_oind_array(N)
    μyy_oind3d = create_oind_array(N)
    μzz_oind3d = create_oind_array(N)
    μoo_oind3d = create_oind_array(N)

    boundft = SVec(EE,EE,EE)
    assign_param!(ε3d, (μxx_oind3d,μyy_oind3d,μzz_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.ghosted.τl, g3.isbloch)
    assign_param!(ε3d, tuple(μoo_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.ghosted.τl, g3.isbloch)
    assign_param!(μ3d, (εxx_oind3d,εyy_oind3d,εzz_oind3d), oind2shp, oind2μind, μind2μ, ft2gt.(HH,boundft), g3.ghosted.τl, g3.isbloch)
    assign_param!(μ3d, tuple(εoo_oind3d), oind2shp, oind2μind, μind2μ, ft2gt.(HH,boundft), g3.ghosted.τl, g3.isbloch)

    # Perform smoothing.
    smooth_param!(ε3d, (εxx_oind3d,εyy_oind3d,εzz_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.l, g3.ghosted.l, g3.σ, g3.ghosted.∆τ)
    smooth_param!(ε3d, tuple(εoo_oind3d), oind2shp, oind2εind, εind2ε, ft2gt.(EE,boundft), g3.l, g3.ghosted.l, g3.σ, g3.ghosted.∆τ)

    ∆l = g3.∆l[nDL]
    ∆l′ = g3.∆l[nPR]
    ∆l′⁻¹ = map(x->inv.(x), ∆l′)

    Mε = create_paramop(ε3d, boundft.!=EE, ∆l, ∆l′⁻¹, g3.isbloch, order_cmpfirst=true)

    # Ml
    ∆lprim = g3.∆l[nPR]
    ∆ldual = g3.∆l[nDL]

    Nx, Ny, Nz = N
    nX, nY, nZ = 1, 2, 3
    ∆X = repeat(∆ldual[nX], outer=(1,Ny,Nz))
    ∆Y = repeat(reshape(∆ldual[nY], (1,Ny,1)), outer=(Nx,1,Nz))
    ∆Z = repeat(reshape(∆ldual[nZ], (1,1,Nz)), outer=(Nx,Ny,1))

    ∆L = [reshape(∆X,1,:); reshape(∆Y,1,:); reshape(∆Z,1,:)]
    ∆l = ∆L[:]
    Ml⁻¹ = sparse(Diagonal(1 ./ ∆l))

    # Ma
    ∆X = repeat(∆lprim[nX], outer=(1,Ny,Nz))
    ∆Y = repeat(reshape(∆lprim[nY], (1,Ny,1)), outer=(Nx,1,Nz))
    ∆Z = repeat(reshape(∆lprim[nZ], (1,1,Nz)), outer=(Nx,Ny,1))

    ∆YZ = ∆Y .* ∆Z
    ∆ZX = ∆Z .* ∆X
    ∆XY = ∆X .* ∆Y

    ∆A = [reshape(∆YZ,1,:); reshape(∆ZX,1,:); reshape(∆XY,1,:)]
    ∆a = ∆A[:]
    Ma = sparse(Diagonal(∆a))

    Mε_sym = Ma * Mε * Ml⁻¹

    @test Mε_sym ≈ transpose(Mε_sym)

    # Code I used to reveal the matrix structure while debugging:
    # (i,j) = (nX,nY); full(real.(Mε3))[27(i-1)+1:27i, 27(j-1)+1:27j]
end  # @testset "create_paramop"

# Need to test if PML still works after the field-averaging operator applied to the material
# parameters, because effectively it makes the material inhomogeneous along the direction
# normal to the PML interface.

end  # @testset "param"


text:
using StatisticalRethinkingBase
using Test


text:
"""
this kernel will prepare da
"""
module PrepareArrtoBool
using CUDA, Logging,..CUDAGpuUtils, Logging,StaticArrays, ..IterationUtils,..BitWiseUtils, ..ReductionUtils, ..CUDAAtomicUtils,..MetaDataUtils,..HFUtils
export @planeIter,getLargeForBoolKernel,getSmallForBoolKernel,@getBoolCubeKernel,@localAllocations,@uploadLocalfpFNCounters,@uploadMinMaxesToShmem,@uploadDataToMetaData,@finalGlobalSet

"""
allocates in the local, register in shared memory
"""
macro localAllocations()

    return esc(quote
    anyPositive = false # true If any bit will bge positive in this array - we are not afraid of data race as we can set it multiple time to true
    #creates shared memory and initializes it to 0
    # shmemSum = @cuStaticSharedMem(Float32,(32,2))
    locFps= UInt32(0)
    locFns= UInt32(0)
    offsetIter= UInt8(0)
    #storing data about block in a forrmat where each Int32 number is representing a part of data block with constant x and y and varia ble z position
    # shmemblockData = @cuDynamicSharedMem(Float32,dataBdim[1], dataBdim[2])
    shmemblockData = @cuDynamicSharedMem(UInt32,(dataBdim[1], dataBdim[2]))

    ######## needed for establishing min and max values of blocks that are intresting us
     minX =@cuStaticSharedMem(Float32, 1)
     maxX= @cuStaticSharedMem(Float32, 1)
     minY = @cuStaticSharedMem(Float32, 1)
     maxY= @cuStaticSharedMem(Float32, 1)
     minZ = @cuStaticSharedMem(Float32, 1)
     maxZ= @cuStaticSharedMem(Float32, 1)

     #resetting
     minX[1]= Float32(1110.0)
     maxX[1]= Float32(0.0)
     minY[1]= Float32(1110.0)
     maxY[1]= Float32(0.0)
     minZ[1]= Float32(1110.0)
     maxZ[1]= Float32(0.0)
     #in shared memory

#####needed for fp fn sums
     #1 - false negative; 2- false positive
     #locArr= (UInt32(0.0), UInt32(0.0))# for global fp fn sums
     #locArrB= (Float32(0.0), Float32(0.0))# for local fp fn sums
    #  1)   Left FP
    #  2)   Left FN
    #  3)   Right FP
    #  4)   Right FN
    #  5)   Posterior FP
    #  6)   Posterior FN
    #  7)   Anterior FP
    #  8)   Anterior FN
    #  9)   Top FP
    #  10)   Top FN
    #  11)   Bottom FP
    #  12)   Bottom FN
    #13)   main part FP
    #14)   main Part FN
    localQuesValues= @cuStaticSharedMem(UInt32, 14)

    #making sure they are initialized all to zeros

     sync_threads()
end)
end

"""
invoked on each lane and on the basis of its position will update the number of fp or fn in given queue
"""
macro uploadLocalfpFNCounters()
   return esc(quote
   atomicallyAddToSpot(localQuesValues,getIndexOfQueue(xpos,ypos,zpos,dataBdim,boolSegm),1)
    end)
end

"""
invoked after we gone through data block and now we save data into shared memory
"""
macro uploadMinMaxesToShmem()
    return  esc(quote

    # if( anyPositive)
    #     @ifXY 1 1  CUDA.@cuprint "xMeta+1 $(xMeta+1) yMeta+1 $(yMeta+1) zMeta+1 $(zMeta+1) anyPositive $(anyPositive) \n"
    #     @ifXY 1 1  CUDA.@cuprint "xMeta+1 $(xMeta+1) yMeta+1 $(yMeta+1) zMeta+1 $(zMeta+1) anyPositive $(anyPositive) \n"
    # end
    # @ifXY 1 1 if(anyPositive) CUDA.@cuprint "aaaaaa minX[1] $(minX[1]) xMeta+1 $(xMeta+1)  " end

        @ifXY 1 1 if(anyPositive) minX[1]= min(minX[1],xMeta+1) end
        @ifXY 1 2 if(anyPositive) maxX[1]= max(maxX[1],xMeta+1) end
        @ifXY 2 1 if(anyPositive) minY[1]= min(minY[1],yMeta+1) end
        @ifXY 2 2 if(anyPositive) maxY[1]= max(maxY[1],yMeta+1) end
        @ifXY 3 1 if(anyPositive) minZ[1]= min(minZ[1],zMeta+1) end
        @ifXY 3 2 if(anyPositive) maxZ[1]= max(maxZ[1],zMeta+1) end
    end)

end

"""
invoked after we gone through data block and now we save data into appropriate spots in metadata of this metadata block
"""
macro uploadDataToMetaData()
    esc(quote
    #now we should also add the total value by adding all fp or fn values required
    @ifY 1 if(threadIdxX()<15 && anyPositive)
        @setMeta(getBeginingOfFpFNcounts()+ threadIdxX(),localQuesValues[threadIdxX()])

    end
    @ifXY 15 1 if(anyPositive)
        # metaData[xMeta+1,yMeta+1,zMeta+1,(getBeginingOfFpFNcounts()+ 15)]= (localQuesValues[1]+localQuesValues[3]+localQuesValues[5]+localQuesValues[7] +localQuesValues[9]+localQuesValues[11]+localQuesValues[13] )
        @setMeta(getBeginingOfFpFNcounts()+ 15, (localQuesValues[1]+localQuesValues[3]+localQuesValues[5]+localQuesValues[7] +localQuesValues[9]+localQuesValues[11]+localQuesValues[13] )  )
        end
    @ifXY 16 1 if(anyPositive)
        # metaData[xMeta+1,yMeta+1,zMeta+1,(getBeginingOfFpFNcounts()+ 16)]= (localQuesValues[2]+localQuesValues[4]+localQuesValues[6]+localQuesValues[8]+localQuesValues[10]+localQuesValues[12]+localQuesValues[14])
        @setMeta(getBeginingOfFpFNcounts()+ 16, (localQuesValues[2]+localQuesValues[4]+localQuesValues[6]+localQuesValues[8]+localQuesValues[10]+localQuesValues[12]+localQuesValues[14]) )
        end

    end)

end#uploadDataToMetaData

"""
invoked after all of the data was scanned so after we will do atomics between blocks we will know
    the minimal and maximal in each dimensions
"""
macro  finalGlobalSet()
    esc(quote
        offsetIter=1
        @redWitAct(offsetIter,shmemblockData,  locFns,+,     locFps,+   )
        @addAtomic(shmemblockData,fn,fp)
        # if(minX[1]<100)
        #     # CUDA.@cuprint "aaaaaaaa  minxRes[1] $(minxRes[1])  minX $(minX[1])\n"
        # end
        @ifXY 1 1 atomicMinSet(minxRes,minX[1])
        @ifXY 1 2 atomicMaxSet(maxxRes,maxX[1])

        @ifXY 2 1 atomicMinSet(minyRes,minY[1])
        @ifXY 2 2 atomicMaxSet(maxyRes,maxY[1])

        @ifXY 3 1 atomicMinSet(minzRes,minZ[1])
        @ifXY 3 2 atomicMaxSet(maxzRes,maxZ[1])
    end)
end

"""
we need to give back number of false positive and false negatives and min,max x,y,x of block containing all data
adapted from https://github.com/JuliaGPU/CUDA.jl/blob/afe81794038dddbda49639c8c26469496543d831/src/mapreduce.jl
goldGPU - array holding 3 dimensional data of gold standard bollean array
segmGPU - array with 3 dimensional the data we want to compare with gold standard
reducedGold - the smallest boolean block (3 dim array) that contains all positive entris from both masks
reducedSegm - the smallest boolean block (3 dim array) that contains all positive entris from both masks
numberToLooFor - number we will analyze whether is the same between two sets
cuda arrays holding just single value wit atomically reduced result
,fn,fp
,minxRes,maxxRes
,minyRes,maxyRes
,minZres,maxZres
dataBdim - dimensions of data block
metaDataDims - dimensions of the metadata
loopXMeta,loopYZMeta- indicates how many times we need to iterate over the metadata
inBlockLoopX,inBlockLoopY,inBlockLoopZ - indicates how many times we need to iterate over the data block using our size of thread block
                                          basically data block size will be established by the thread block size of main kernel
"""
# function getBoolCubeKernel(goldGPU
#         ,segmGPU
#         ,numberToLooFor::T
#         ,reducedGoldA
#         ,reducedSegmA
#         ,reducedGoldB
#         ,reducedSegmB
#         ,fn
#         ,fp
#         ,minxRes
#         ,maxxRes
#         ,minyRes
#         ,maxyRes
#         ,minzRes
#         ,maxzRes
#         ,dataBdim
#         ,metaData
#         ,metaDataDims
#         ,mainArrDims
#         ,loopXMeta,loopYZMeta
#         ,inBlockLoopX,inBlockLoopY,inBlockLoopZ
# ) where T

macro getBoolCubeKernel()
 return esc(quote
    @localAllocations()
    #we need nested x,y,z iterations so we will iterate over the matadata and on its basis over the  data in the main arrays
    #first loop over the metadata

    HFUtils.@iter3dOuter(metaDataDims,loopMeta,metaDataLength,
         begin
        # inner loop is over the data indicated by metadata
        @iterDataBlockZdeepest(mainArrDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,xMeta,yMeta,zMeta
                         ,begin
                                boolGold=goldGPU[x,y,z]==numberToLooFor
                                boolSegm=segmGPU[x,y,z]==numberToLooFor

                                #we set all bits so we do not need to reset
                                @setBitTo(offsetIter,zpos,boolGold)

                                # @setBitTo((shmemblockData[xpos,ypos]),zpos,boolGold)
                                #we need to also collect data about how many fp and fn we have in main part and borders
                                #important in case of corners we will first analyze z and y dims and z dim on last resort only !

                                # #in case some is positive we can go futher with looking for max,min in dims and add to the new reduced boolean arrays waht we are intrested in
                                if(boolGold  || boolSegm)
                                        anyPositive=true
                                        if((boolGold  ⊻ boolSegm))
                                            @uploadLocalfpFNCounters()
                                            locFps+=boolSegm
                                            locFns+=boolGold
                                        end# if (boolGold  ⊻ boolSegm)
                                    #now

                                end#if boolGold  || boolSegm
                            end#ex
                          #additional after X
                            ,begin
                        # CUDA.@cuprint "xpos $(xpos) ypos $(ypos) loopXdim $(inBlockLoopX) xdim $(xdim) \n "
                        #here we iterated over all z dimension so shmemblockData[xpos,ypos] is ready to be uploaded to global memory
                        if(offsetIter>0)
                           @inbounds reducedGoldA[x,y,(zMeta+1)]=offsetIter
                        end
                        end)

                  # #now we are just after we iterated over a single data block  we need to we save data about border data blocks
                  anyPositive = sync_threads_or(anyPositive)

                 @uploadMinMaxesToShmem()

                 offsetIter=0
                 #in order to reduce used shared memory we are setting values of output array separately
                 @iterDataBlockZdeepest(mainArrDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,xMeta,yMeta,zMeta
                 ,begin
                    #if((segmGPU[x,y,z]==numberToLooFor)) @setBitTo1(offsetIter,zpos) end
                      @setBitTo(offsetIter,zpos,(segmGPU[x,y,z]==numberToLooFor))
                       #we set all bits so we do not need to reset
                    #@setBitTo(shmemblockData[xpos,ypos],zpos,(segmGPU[x,y,z]==numberToLooFor))
                end,begin
                #here we iterated over all z dimension so shmemblockData[xpos,ypos] is ready to be uploaded to global memory
                if(offsetIter>0)
                    @inbounds reducedSegmA[x,y,(zMeta+1)]=offsetIter
                end
                end)
                sync_threads()

                    #we want to invoke this only once per data block
                    #save the data about number of fp and fn of this block and accumulate also this sum for global sum
                    #doing all on first warp
                    @uploadDataToMetaData()
                    # if(anyPositive)
                    #     @ifXY 4 1  CUDA.@cuprint "xMeta $(xMeta+1)  ,yMeta $(yMeta+1),zMeta $(zMeta+1) \n"
                    # end
                    #invoked after we gone through data block and now we save data into shared memory

                    sync_threads()
                    #resetting
                    anyPositive= false  #reset

                    @ifY 2 if(threadIdxX()<15)
                        localQuesValues[threadIdxX()]=0
                    end
               sync_threads()

            end) #outer loop
    #             #consider ceating tuple structure where we will have  number of outer tuples the same as z dim then inner tuples the same as y dim and most inner tuples will have only the entries that are fp or fn - this would make us forced to put results always in correct spots

     @finalGlobalSet()

   return
end)#quote
   end

"""
creates small memory footprint GPU variables for getBoolCubeKernel
  return  minX,maxX,minY,maxY,minZ,maxZ,fn,fp
"""
function getSmallForBoolKernel()
    return (CuArray([Float32(1110.0) ])
    , CuArray([Float32(0.0)])
    , CuArray([Float32(1110.0)])
    , CuArray([Float32(0.0)    ])
    , CuArray([Float32(1110.0)])
    , CuArray([Float32(0.0) ])
    ,CuArray([UInt32(0)])
    ,CuArray([UInt32(0)]))
end

"""
creates large memory footprint GPU variables for getBoolCubeKernel
    return reducedGoldA,reducedSegmA,reducedGoldB,reducedSegmB
"""
function getLargeForBoolKernel(mainArrDims,dataBdim)
    #this is in order to be sure that array is divisible by data block so we reduce necessity of boundary checks
    xDim= cld(mainArrDims[1],dataBdim[1])*dataBdim[1]
    yDim = cld(mainArrDims[2],dataBdim[2])*dataBdim[2]
    zDim = cld(mainArrDims[3],dataBdim[3])
    newDims = (xDim,yDim,zDim)
return (
    CUDA.zeros(UInt32,(newDims)),CUDA.zeros(UInt32,(newDims))
    )

end

# """
# iterating over shmemblockData
# """
# macro planeIter(loopXinPlane,loopYinPlane,maxXdim, maxYdim,ex)
#     mainExp = generalizedItermultiDim(
#     arrDims=:()
#     ,loopXdim=loopXinPlane
#     ,loopYdim=loopYinPlane
#     ,yCheck = :(y <=$maxYdim)
#     ,xCheck = :(x <=$maxXdim )
#     ,is3d = false
#     , ex = ex)
#       return esc(:( $mainExp))
# end

end#module


text:
export EncodingOp, lrEncodingOp, fourierEncodingOp, encodingOps_simple,
       encodingOps_parallel, encodingOp_multiEcho, encodingOp_multiEcho_parallel

"""
    encodingOps_simple(acqData::AcquisitionData, shape::NTuple{D,Int64}
                              ; kargs...) where D

generates an Array of LinearOperators which describe signal encoding of the individual
contrasts in an MRI acquisition (for a given slice).

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `shape::NTuple{D,Int64}`              - size of image to be encoded/reconstructed
"""
function encodingOps_simple(acqData::AcquisitionData{T}, shape::NTuple{D,Int64}; kargs...) where {T,D}
  numContr = numContrasts(acqData)
  tr = [trajectory(acqData,i) for i=1:numContr]
  idx = acqData.subsampleIndices
  return [fourierEncodingOp(shape, tr[i], "fast",
                 subsampleIdx=idx[i]; kargs...) for i=1:numContr]
end

"""
    encodingOps_parallel(acqData::AcquisitionData{T}, shape::NTuple{D,Int64}
                              , senseMaps::Array{Complex{T}}
                              ; kargs...) where {T,D}

generates an Array of LinearOperators which describe signal encoding of the individual
contrasts in an MRI acquisition. The different coils are taken into account
in terms of their sensitivities

# Arguments
* `acqData::AcquisitionData{T}`            - AcquisitionData object
* `shape::NTuple{D,Int64}`              - size of image to be encoded/reconstructed
* `senseMaps::Array{Complex{T}}`        - coil sensitivities
"""
function encodingOps_parallel(acqData::AcquisitionData{T}, shape::NTuple{D,Int64}
                                , senseMaps::Array{Complex{T},4}
                                ; slice=1, kargs...) where {T,D}

  smaps = ( D==2 ? senseMaps[:,:,slice,:] : senseMaps )

  numContr, numChan = numContrasts(acqData), numChannels(acqData)
  # fourier operators
  ft = encodingOps_simple(acqData, shape; slice=slice, kargs...)
  S = SensitivityOp(reshape(smaps,:,numChan),1)
  Op = [ diagOp(ft[i], numChan) ∘ S for i=1:numContr]

  return Op
end

"""
    encodingOp_multiEcho(acqData::AcquisitionData, shape::NTuple{D,Int64}; kargs...) where D

generates a LinearOperator which describe combined signal encoding of all
the contrasts in an MRI acquisition (for a given slice).

# Arguments
* `acqData::AcquisitionData`            - AcquisitionData object
* `shape::NTuple{D,Int64}`              - size of image to be encoded/reconstructed
"""
function encodingOp_multiEcho(acqData::AcquisitionData, shape::NTuple{D,Int64}
                                ; kargs...) where D
  # fourier operators
  ft = encodingOps_simple(acqData, shape; kargs...)
  return diagOp(ft...)
end

"""
    encodingOp_multiEcho_parallel(acqData::AcquisitionData{T}, shape::NTuple{D,Int64}
                                          , senseMaps::Array{Complex{T}}
                                          ; kargs...) where {T,D}

generates a LinearOperator which describe combined signal encoding of all
the contrasts in an MRI acquisition (for a given slice). The different coils are taken into account
in terms of their sensitivities

# Arguments
* `acqData::AcquisitionData{T}`            - AcquisitionData object
* `shape::NTuple{2,Int64}`              - size of image to be encoded/reconstructed
* `senseMaps::Array{Complex{T}}`        - coil sensitivities
"""
function encodingOp_multiEcho_parallel(acqData::AcquisitionData{T}, shape::NTuple{D,Int64}
                                          , senseMaps::Array{Complex{T}}
                                          ; slice::Int64=1, kargs...) where {T,D}

  smaps = ( D==2 ? senseMaps[:,:,slice,:] : senseMaps )

  numChan = numChannels(acqData)
  # fourier operators
  ft = encodingOps_simple(acqData, shape; kargs...)
  S = SensitivityOp(reshape(smaps,:,numChan),numContrasts(acqData))
  return diagOp(repeat(ft, outer=numChan)...) ∘ S
end

###################################
# Encoding with low rank projection
###################################
function lrEncodingOp(acqData::AcquisitionData, shape, params::Dict; numContr::Int64=1, parallel::Bool=false)

  numChan = numChannels(acqData)
  # low rank operator
  N=prod(shape)
  subspace = get( params, :phi, Matrix{Float64}(I,N,N) )
  K = size(subspace,2)
  Φ = MapSliceOp(subspace[:,:,1],2,(N, K, numChan), (N, numContr, numChan))

  # Fourier Operator
  tr = trajectory(acqData,1)
  ft = fourierEncodingOp2d(shape, tr, "fast"; params...)

  # coil sensitivities in case of SENSE-reconstruction
  if parallel
    S = SensitivityOp(params[:senseMaps][:,:],K)
    E = diagOp( [ft for i=1:numChan*K]... )*S
  else
    E = diagOp( [ft for i=1:K]... )
  end

  #  TODO: sampling operator in case of undersampled cartesian acquisitions
  # vectorize sampling pattern
  if get(params,:sampling, nothing) == "binary"
    M = SamplingOp(Array{Bool}(acqData.subsampleIndices))
  else
    # sampling idx for all contrasts but one coil
    subIdx = hcat(acqData.subsampleIndices...)
    M = SamplingOp( hcat([subIdx for c=1:numChan]...), (N, numContr, numChan))
  end

  return M ∘ (Φ ∘ E)

end

"""
return Fourier encoding operator (either Explicit or NFFT)
  * `opname` : "explicit" or "fast"
  * `slice` : slice to which the operator will be applied
  echoImage : calculate signal evolution relative to the echo time
"""
function fourierEncodingOp(shape::NTuple{D,Int64}, tr::Trajectory{T}, opName::String;
          subsampleIdx::Vector{Int64}=Int64[], slice::Int64=1, correctionMap::Array{Complex{T}}=Complex{T}[],
          echoImage::Bool=true,  kargs...) where {T,D}

  # extract proper portion of correctionMap
  if !isempty(correctionMap)
    cmap = ( D==2 ? correctionMap[:,:,slice] : correctionMap )
  end
  # Fourier transformations
  if opName=="explicit"
    @debug "ExplicitOp"
    ftOp = ExplicitOp(shape, tr, cmap, echoImage=echoImage)
  elseif opName=="fast"
    @debug "NFFT-based Op"
    if !isempty(correctionMap) && correctionMap!=zeros(Complex{T},size(correctionMap))
      ftOp = FieldmapNFFTOp(shape, tr, cmap, echoImage=echoImage; kargs...)
    elseif isCartesian(tr)
      @debug "FFTOp"
      ftOp = FFTOp(Complex{T}, shape; unitary=false)
    else
      ftOp = NFFTOp(shape, tr; kargs...)
    end
  else
    @error "opName $(opName) is not known"
  end

  # subsampling
  if !isempty(subsampleIdx) && (subsampleIdx != collect(1:size(tr,2))) && isCartesian(tr)
    if D==2
      S = SamplingOp(subsampleIdx,(tr.numSamplingPerProfile,tr.numProfiles))
    else
      S = SamplingOp(subsampleIdx,(tr.numSamplingPerProfile,tr.numProfiles,tr.numSlices))
    end
    return S ∘ ftOp
  else
    return ftOp
  end
end


text:
using MsgPackRpcServer
using Base.Test

module RemoteFunctions
  function hello()
    "hi from remote server"
  end

  function reverse(bool)
    !bool
  end

  function add(x, y)
    x+y
  end

  function array_map(a,b,c...)
    a=map([a,b,c...]) do x
      x*x
    end
    (a...)
  end

  function double_hash_count(h)
    h["count"] = h["count"]*h["count"]
    h
  end
end
@async MsgPackRpcServer.run(8080, RemoteFunctions)

sleep(1) # temporarily wait here until MsgPackRpcServer boot finished
run(`ruby test_client.rb`)

max_stars_repo_name: remore/MsgPackRpcServer

max_issues_repo_name: remore/MsgPackRpcServer

max_forks_repo_name: remore/MsgPackRpcServer


text:
## convenient way to get library information
using BinaryProvider

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    println("Deps path: $depsjl_path")
    error("VideoIO not installed properly, run Pkg.build(\"VideoIO\"), restart Julia and try again")
end
include(depsjl_path)

include("version.jl")

av_load_path = joinpath(dirname(@__FILE__), "ffmpeg")


text:
using StaticJulia
using Base.Test

fcall(a, b) = a + b
function test(a, b)
    c = fcall(a, b)
    c == b
end

mod = llvm_module("Test2")

f = llvm_compile(test, (Int, Int), mod)

@code_llvm test(5, 2)

@test f(5, 2) == test(5, 2)
@test f(0, 5) == test(0, 5)

using BenchmarkTools
b1 = @benchmark f(5, 2)
b2 = @benchmark test(5, 2)
# BenchmarkTools is really cool! :P
x = judge(minimum(b1), minimum(b2))
# We can see, our own compiled function has a 600% regression in minimum runtime.
# This is because we haven't run any optimization passes yet!
# Take a look at both IR's to figure out what got optimized!


text:
function get_random_string()
    return string(utc_to_string(now_localzone()), "-", rand(UInt32))::String
end


text:
"""
    wins_and_losses(create_warrior::Type{<:Warrior}, create_opponent::Type{<:Warrior})

Compute number of times `warror` wins against `opponent`. Why is this useful?
Because there is a lot of randomness involved one type of warrior may not always
win over another in a single match. Hence we run numerous battles because.

Whatever warrior has the best odds of winning will have more wins when running
numerous battles.
"""
function wins_and_losses(createwarrior::Type{<:Warrior}, createopponent::Type{<:Warrior})
   wins = 0
   losses = 0

   for _ in 1:150
       warrior = createwarrior("warrior", 30)
       opponent = createopponent("opponent", 30)

       while true
          attack!(warrior, opponent)

          if opponent.health <= 0
              wins += 1
              break
          end

          attack!(opponent, warrior)
          if warrior.health <= 0
              losses += 1
              break
          end
       end

   end

   wins, losses
end

max_stars_repo_name: ordovician/BattleSim

max_issues_repo_name: ordovician/BattleSim

max_forks_repo_name: ordovician/BattleSim


text:
using SparseArrays

"""
    PartitionMask{A}(A_1::A, A_2::A, A_3::A) where {A}

This is used to partition and recombine a vector into 3 disjoint "subvectors".

Implements
- `partition(m::PartitionMask, x)`: partitions `x` into 3 disjoint "subvectors"
- `combine(m::PartitionMask, x_1, x_2, x_3)`: combines 3 disjoint vectors into a single one

Note that `PartitionMask` is _not_ a `Bijector`. It is indeed a bijection, but
does not follow the `Bijector` interface.

Its main use is in `Coupling` where we want to partition the input into 3 parts,
one part to transform, one part to map into the parameter-space of the transform applied
to the first part, and the last part of the vector is not used for anything.

# Examples
```julia-repl
julia> using Bijectors: PartitionMask, partition, combine

julia> m = PartitionMask(3, [1], [2]) # <= assumes input-length 3
PartitionMask{Bool,SparseArrays.SparseMatrixCSC{Bool,Int64}}(
  [1, 1]  =  true,
  [2, 1]  =  true,
  [3, 1]  =  true)

julia> # Partition into 3 parts; the last part is inferred to be indices `[3, ]` from
       # the fact that `[1]` and `[2]` does not make up all indices in `1:3`.
       x1, x2, x3 = partition(m, [1., 2., 3.])
([1.0], [2.0], [3.0])

julia> # Recombines the partitions into a vector
       combine(m, x1, x2, x3)
3-element Array{Float64,1}:
 1.0
 2.0
 3.0
```
Note that the underlying `SparseMatrix` is using `Bool` as the element type. We can also
specify this to be some other type using the `sp_type` keyword:
```julia-repl
julia> m = PartitionMask{Float32}(3, [1], [2])
PartitionMask{Float32,SparseArrays.SparseMatrixCSC{Float32,Int64}}(
  [1, 1]  =  1.0,
  [2, 1]  =  1.0,
  [3, 1]  =  1.0)
```
"""
struct PartitionMask{T, A}
    A_1::A
    A_2::A
    A_3::A

    # Only make it possible to construct using matrices
    PartitionMask(A_1::A, A_2::A, A_3::A) where {T<:Real, A <: AbstractMatrix{T}} = new{T, A}(A_1, A_2, A_3)
end

PartitionMask(args...) = PartitionMask{Bool}(args...)

function PartitionMask{T}(
    n::Int,
    indices_1::AbstractVector{Int},
    indices_2::AbstractVector{Int},
    indices_3::AbstractVector{Int}
) where {T<:Real}
    A_1 = sparse(indices_1, 1:length(indices_1), one(T), n, length(indices_1))
    A_2 = sparse(indices_2, 1:length(indices_2), one(T), n, length(indices_2))
    A_3 = sparse(indices_3, 1:length(indices_3), one(T), n, length(indices_3))

    return PartitionMask(A_1, A_2, A_3)
end

PartitionMask{T}(
    n::Int,
    indices_1::AbstractVector{Int},
    indices_2::AbstractVector{Int};
) where {T} = PartitionMask{T}(n, indices_1, indices_2, nothing)

PartitionMask{T}(
    n::Int,
    indices_1::AbstractVector{Int},
    indices_2::AbstractVector{Int},
    indices_3::Nothing,
) where {T} = PartitionMask{T}(n, indices_1, indices_2, setdiff(1:n, indices_1, indices_2))

PartitionMask{T}(
    n::Int,
    indices_1::AbstractVector{Int},
    indices_2::Nothing,
    indices_3::AbstractVector{Int},
) where {T} = PartitionMask{T}(n, indices_1, setdiff(1:n, indices_1, indices_3), indices_3)

"""
    PartitionMask(n::Int, indices)

Assumes you want to _split_ the vector, where `indices` refer to the
parts of the vector you want to apply the bijector to.
"""
function PartitionMask{T}(n::Int, indices) where {T}
    indices_2 = setdiff(1:n, indices)

    # sparse arrays <3
    A_1 = sparse(indices, 1:length(indices), one(T), n, length(indices))
    A_2 = sparse(indices_2, 1:length(indices_2), one(T), n, length(indices_2))

    return PartitionMask(A_1, A_2, spzeros(T, n, 0))
end
function PartitionMask{T}(x::AbstractVector, indices) where {T}
    return PartitionMask{T}(length(x), indices)
end

"""
    combine(m::PartitionMask, x_1, x_2, x_3)

Combines `x_1`, `x_2`, and `x_3` into a single vector.
"""
@inline combine(m::PartitionMask, x_1, x_2, x_3) = m.A_1 * x_1 .+ m.A_2 * x_2 .+ m.A_3 * x_3

"""
    partition(m::PartitionMask, x)

Partitions `x` into 3 disjoint subvectors.
"""
@inline partition(m::PartitionMask, x) = (transpose(m.A_1) * x, transpose(m.A_2) * x, transpose(m.A_3) * x)

# Coupling

"""
    Coupling{F, M}(θ::F, mask::M)

Implements a coupling-layer as defined in [1].

# Examples
```julia-repl
julia> m = PartitionMask(3, [1], [2]) # <= going to use x[2] to parameterize transform of x[1]
PartitionMask{SparseArrays.SparseMatrixCSC{Float64,Int64}}(
  [1, 1]  =  1.0,
  [2, 1]  =  1.0,
  [3, 1]  =  1.0)

julia> cl = Coupling(θ -> Shift(θ[1]), m) # <= will do `y[1:1] = x[1:1] + x[2:2]`;

julia> x = [1., 2., 3.];

julia> cl(x)
3-element Array{Float64,1}:
 3.0
 2.0
 3.0

julia> inv(cl)(cl(x))
3-element Array{Float64,1}:
 1.0
 2.0
 3.0

julia> coupling(cl) # get the `Bijector` map `θ -> b(⋅, θ)`
Shift

julia> couple(cl, x) # get the `Bijector` resulting from `x`
Shift{Array{Float64,1},1}([2.0])
```

# References
[1] Kobyzev, I., Prince, S., & Brubaker, M. A., Normalizing flows: introduction and ideas, CoRR, (),  (2019).
"""
struct Coupling{F, M} <: Bijector{1} where {F, M <: PartitionMask}
    θ::F
    mask::M
end

function Coupling(θ, n::Int)
    idx = div(n, 2)
    return Coupling(θ, PartitionMask(n, 1:idx))
end

function Coupling(cl::Coupling, mask::PartitionMask)
    return Coupling(cl.θ, mask)
end

"Returns the constructor of the coupling law."
coupling(cl::Coupling) = cl.θ

"Returns the coupling law constructed from `x`."
function couple(cl::Coupling, x::AbstractVector)
    # partition vector using `cl.mask::PartitionMask`
    x_1, x_2, x_3 = partition(cl.mask, x)

    # construct bijector `B` using θ(x₂)
    b = cl.θ(x_2)

    return b
end

function (cl::Coupling)(x::AbstractVector)
    # partition vector using `cl.mask::PartitionMask`
    x_1, x_2, x_3 = partition(cl.mask, x)

    # construct bijector `B` using θ(x₂)
    b = cl.θ(x_2)

    # recombine the vector again using the `PartitionMask`
    return combine(cl.mask, b(x_1), x_2, x_3)
end
(cl::Coupling)(x::AbstractMatrix) = eachcolmaphcat(cl, x)

function (icl::Inverse{<:Coupling})(y::AbstractVector)
    cl = icl.orig

    y_1, y_2, y_3 = partition(cl.mask, y)

    b = cl.θ(y_2)
    ib = inv(b)

    return combine(cl.mask, ib(y_1), y_2, y_3)
end
(icl::Inverse{<:Coupling})(y::AbstractMatrix) = eachcolmaphcat(icl, y)

function logabsdetjac(cl::Coupling, x::AbstractVector)
    x_1, x_2, x_3 = partition(cl.mask, x)
    b = cl.θ(x_2)

    # `B` might be 0-dim in which case it will treat `x_1` as a batch
    # therefore we sum to ensure such a thing does not happen
    return sum(logabsdetjac(b, x_1))
end

function logabsdetjac(cl::Coupling, x::AbstractMatrix)
    return [logabsdetjac(cl, view(x, :, i)) for i in axes(x, 2)]
end


text:
let
    #= line 1 =#
    var"|a|b_op_Getter"::Function = identity::Function
    var"|a|c_f"::Any = (x->begin
                    #= none:1 =#
                    #= none:4 =#
                    x
                end)::Any
    var"|a|d_x"::Any = 1::Any
    var"|a|e_g"::Function = function (var"|a|f_x",)
                #= /home/tainez/github/MLFS/test/a.jl:20 =#
                function (var"|a|g_y",)
                        #= /home/tainez/github/MLFS/test/a.jl:20 =#
                        (var"|a|b_op_Getter"::Function)((var"|a|c_f"::Any)(var"|a|g_y"::Any)::Any)::Any
                    end::Function
            end::Function
    var"|a|"::Tuple = ()
    a::Tuple = (("f32", "f16", "i64", "i32", "op_Getter", "str", "f64", "bool", "i16", "x", "module", "field", "char", "fieldnames", "g", "i8", "f", "Type"), (Type, Type, Type, Type, var"|a|b_op_Getter", Type, Type, Type, Type, var"|a|d_x", Type, Type, Type, Type, var"|a|e_g", Type, var"|a|c_f", Type))
    0
end