text:
module TriangularSolve

using VectorizationBase, LinearAlgebra #LoopVectorization
using VectorizationBase: vfnmadd_fast, AbstractStridedPointer, AbstractMask, stridedpointer_preserve, zero_offsets, gesp, StridedPointer
using CloseOpenIntervals: CloseOpen, SafeCloseOpen
using Static
using IfElse: ifelse
using LoopVectorization
using Polyester

@generated function solve_AU(A::VecUnroll{Nm1}, spu::AbstractStridedPointer, noff, ::Val{UNIT}) where {Nm1, UNIT}
  A_n_expr = UNIT ? :nothing : :(A_n = Base.FastMath.div_fast(A_n, U_n_n))
  N = Nm1 + 1
  quote
    $(Expr(:meta,:inline))
    Ad = VectorizationBase.data(A)
    Base.Cartesian.@nexprs $N n -> begin
      A_n = Ad[n]
      Base.Cartesian.@nexprs $(UNIT ? :(n-1) : :n) m -> begin
        U_m_n = vload(spu, (noff+(m-1),noff+(n-1)))
      end
    end
    Base.Cartesian.@nexprs $N n -> begin
      Base.Cartesian.@nexprs n-1 k -> begin
        A_n = Base.FastMath.sub_fast(A_n, Base.FastMath.mul_fast(A_k, U_k_n))
      end
      $A_n_expr
    end
    VecUnroll(Base.Cartesian.@ntuple $N A)
  end
end

# @generated function nmuladd(A::VecUnroll{Nm1},B::AbstractStridedPointer,C::VecUnroll{Nm1}) where {Nm1}
#   N = Nm1 + 1
#   quote
#     $(Expr(:meta,:inline))
#     Ad = VectorizationBase.data(A);
#     Cd = VectorizationBase.data(C);
#     bp = stridedpointer(B)
#     Base.Cartesian.@nexprs $N n -> C_n = Cd[n]
#     Base.Cartesian.@nexprs $N k -> begin
#       A_k = Ad[k]
#       Base.Cartesian.@nexprs $N n -> begin
#         C_n = Base.FastMath.sub_fast(C_n, Base.FastMath.mul_fast(A_k, vload(B, (k-1,n-1))))
#       end
#     end
#     VecUnroll(Base.Cartesian.@ntuple $N C)
#   end
# end

# @inline function solve_Wx3W(A11::V, A12::V, A13::V, U::AbstractMatrix, ::StaticInt{W}) where {V<:VecUnroll,W}
#   WS = StaticInt{W}()

#   U11 = view(U,StaticInt(1):WS,StaticInt(1):WS)
#   A11 = solve_AU(A11, U11)

#   U12 = view(U,StaticInt(1):WS, StaticInt(1)+WS:WS*StaticInt(2))
#   A12 = nmuladd(A11, U12, A12)
#   U22 = view(U,StaticInt(1)+WS:WS*StaticInt(2),StaticInt(1)+WS:WS*StaticInt(2))
#   A12 = solve_AU(A12, U22)

#   U13 = view(U,StaticInt(1):WS, StaticInt(1)+WS*StaticInt(2):WS*StaticInt(3))
#   A13 = nmuladd(A11, U13, A13)
#   U23 = view(U,StaticInt(1)+WS:WS*StaticInt(2),StaticInt(1)+WS*StaticInt(2):WS*StaticInt(3))
#   A13 = nmuladd(A12, U23, A13)
#   U33 = view(U,StaticInt(1)+WS*StaticInt(2):WS*StaticInt(3),StaticInt(1)+WS*StaticInt(2):WS*StaticInt(3))
#   A13 = solve_AU(A13, U33)

#   return A11, A12, A13
# end

# @inline function solve_Wx3W!(ap::AbstractStridedPointer{T}, bp::AbstractStridedPointer{T}, U, rowoffset, coloffset) where {T}
#   WS = VectorizationBase.pick_vector_width(T)
#   W = Int(WS)
#   A11 = vload(bp, Unroll{2,1,W,1,W,0x0000000000000000,1}((rowoffset,coloffset)))
#   A12 = vload(bp, Unroll{2,1,W,1,W,0x0000000000000000,1}((rowoffset,coloffset+WS)))
#   A13 = vload(bp, Unroll{2,1,W,1,W,0x0000000000000000,1}((rowoffset,coloffset+WS+WS)))

#   A11, A12, A13 = solve_Wx3W(A11, A12, A13, U, WS)

#   vstore!(ap, A11, Unroll{2,1,W,1,W,0x0000000000000000,1}((rowoffset,coloffset)))
#   vstore!(ap, A12, Unroll{2,1,W,1,W,0x0000000000000000,1}((rowoffset,coloffset+WS)))
#   vstore!(ap, A13, Unroll{2,1,W,1,W,0x0000000000000000,1}((rowoffset,coloffset+WS+WS)))
# end
# @inline function solve_Wx3W!(ap::AbstractStridedPointer{T}, bp::AbstractStridedPointer{T}, U, rowoffset, coloffset, m::VectorizationBase.AbstractMask) where {T}
#   WS = VectorizationBase.pick_vector_width(T)
#   W = Int(WS)
#   A11 = vload(bp, Unroll{2,1,W,1,W,0xffffffffffffffff,1}((rowoffset,coloffset)), m)
#   A12 = vload(bp, Unroll{2,1,W,1,W,0xffffffffffffffff,1}((rowoffset,coloffset+WS)), m)
#   A13 = vload(bp, Unroll{2,1,W,1,W,0xffffffffffffffff,1}((rowoffset,coloffset+WS+WS)), m)

#   A11, A12, A13 = solve_Wx3W(A11, A12, A13, U, WS)

#   vstore!(ap, A11, Unroll{2,1,W,1,W,0xffffffffffffffff,1}((rowoffset,coloffset)), m)
#   vstore!(ap, A12, Unroll{2,1,W,1,W,0xffffffffffffffff,1}((rowoffset,coloffset+WS)), m)
#   vstore!(ap, A13, Unroll{2,1,W,1,W,0xffffffffffffffff,1}((rowoffset,coloffset+WS+WS)), m)
# end

# solve_3Wx3W!(A,B,U::UpperTriangular) = solve_3Wx3W!(A,B,parent(U))
# function solve_3Wx3W!(A::AbstractMatrix{T},B,U) where {T}
#   W = VectorizationBase.pick_vector_width(T)
#   ap = stridedpointer(A);
#   bp = stridedpointer(B);
#   solve_Wx3W!(ap, bp, U, StaticInt(1),         StaticInt(1))
#   solve_Wx3W!(ap, bp, U, StaticInt(1) + W,     StaticInt(1))
#   solve_Wx3W!(ap, bp, U, StaticInt(1) + W + W, StaticInt(1))
# end

@inline maybestore!(p, v, i) = vstore!(p, v, i)
@inline maybestore!(::Nothing, v, i) = nothing

@inline maybestore!(p, v, i, m) = vstore!(p, v, i, m)
@inline maybestore!(::Nothing, v, i, m) = nothing

@inline function store_small_kern!(spa, sp, v, spu, i, n, mask, ::Val{true})
  vstore!(spa, v, i, mask)
  vstore!(sp, v, i, mask)
end
@inline store_small_kern!(spa, ::Nothing, v, spu, i, n, mask, ::Val{true}) = vstore!(spa, v, i, mask)

@inline function store_small_kern!(spa, sp, v, spu, i, n, mask, ::Val{false})
  x = v / vload(spu, (n,n))
  vstore!(spa, x, i, mask)
  vstore!(sp, x, i, mask)
end
@inline store_small_kern!(spa, ::Nothing, v, spu, i, n, mask, ::Val{false}) = vstore!(spa, v / vload(spu, (n,n)), i, mask)

@inline function store_small_kern!(spa, sp, v, spu, i, n, ::Val{true})
  vstore!(spa, v, i)
  vstore!(sp, v, i)
end
@inline store_small_kern!(spa, ::Nothing, v, spu, i, n, ::Val{true}) = vstore!(spa, v, i)

@inline function store_small_kern!(spa, sp, v, spu, i, n, ::Val{false})
  x = v / vload(spu, (n,n))
  vstore!(spa, x, i)
  vstore!(sp, x, i)
end
@inline store_small_kern!(spa, ::Nothing, v, spu, i, n, ::Val{false}) = vstore!(spa, v / vload(spu, (n,n)), i)

function BdivU_small_kern!(spa::AbstractStridedPointer{T}, sp, spb::AbstractStridedPointer{T}, spu::AbstractStridedPointer{T}, N, mask, ::Val{UNIT}) where {T,UNIT}
  W = VectorizationBase.pick_vector_width(T)
  for n ∈ CloseOpen(N)
    Amn = vload(spb, (MM(W, StaticInt(0)),n), mask)
    for k ∈ SafeCloseOpen(n)
      Amn = vfnmadd_fast(vload(spa, (MM(W, StaticInt(0)),k), mask), vload(spu, (k,n)), Amn)
    end
    store_small_kern!(spa, sp, Amn, spu, (MM(W, StaticInt(0)),n), n, mask, Val{UNIT}())
  end
end
function BdivU_small_kern_u!(spa::AbstractStridedPointer{T}, sp, spb::AbstractStridedPointer{T}, spu::AbstractStridedPointer{T}, N, ::StaticInt{U}, ::Val{UNIT}) where {T,U,UNIT}
  W = Int(VectorizationBase.pick_vector_width(T))
  for n ∈ CloseOpen(N)
    Amn = vload(spb, Unroll{1,W,U,1,W,0x0000000000000000,1}((StaticInt(0),n)))
    for k ∈ SafeCloseOpen(n)
      Amk = vload(spa, Unroll{1,W,U,1,W,0x0000000000000000,1}((StaticInt(0),k)))
      Amn = vfnmadd_fast(Amk, vload(spu, (k,n)), Amn)
    end
    store_small_kern!(spa, sp, Amn, spu, Unroll{1,W,U,1,W,0x0000000000000000,1}((StaticInt(0),n)), n, Val{UNIT}())
  end
end
# function BdivU_small!(A::AbstractMatrix{T}, B::AbstractMatrix{T}, U::AbstractMatrix{T}) where {T}
#   W = VectorizationBase.pick_vector_width(T)
#   M, N = size(A)
#   m = 0
#   spa = stridedpointer(A)
#   spb = stridedpointer(B)
#   spu = stridedpointer(U)
#   while m < M
#     ml = m+1
#     mu = m+W
#     maskiter = mu > M
#     mask = maskiter ? VectorizationBase.mask(W, M) : VectorizationBase.max_mask(W)
#     for n ∈ 1:N
#       Amn = vload(spb, (MM(W, ml),n), mask)
#       for k ∈ 1:n-1
#         Amn = vfnmadd_fast(vload(spa, (MM(W, ml),k), mask), vload(spu, (k,n)), Amn)
#       end
#       vstore!(spa, Amn / vload(spu, (n,n)), (MM(W, ml),n), mask)
#     end
#     m = mu
#   end
#   # @inbounds @fastmath for m ∈ 1:M
#   #   for n ∈ 1:N
#   #     Amn = B[m,n]
#   #     for k ∈ 1:n-1
#   #       Amn -= A[m,k]*U[k,n]
#   #     end
#   #     A[m,n] = Amn / U[n,n]
#   #   end
#   # end
# end
# function nmuladd!(C,A,B,D)
#   @turbo for n ∈ axes(C,2), m ∈ axes(C,1)
#     Cmn = D[m,n]
#     for k ∈ axes(B,1)
#       Cmn -= A[m,k]*B[k,n]
#     end
#     C[m,n] = Cmn
#   end
# end

@generated function rdiv_solve_W_u!(spc, spb, spa, spu, n, ::StaticInt{W}, ::StaticInt{U}, ::Val{UNIT}) where {W, U, UNIT}
  quote
    # $(Expr(:meta,:inline))
    # here, we just want to load the vectors
    C11 = VectorizationBase.data(vload(spa, Unroll{2,1,$W,1,$W,0x0000000000000000,1}(Unroll{1,$W,$U,1,$W,0x0000000000000000,1}((StaticInt(0),n)))))
    Base.Cartesian.@nexprs $W c -> C11_c = C11[c]
    for nk ∈ SafeCloseOpen(n) # nmuladd
      A11 = vload(spc, Unroll{1,$W,$U,1,$W,0x0000000000000000,1}((StaticInt(0),nk)))
      Base.Cartesian.@nexprs $W c -> C11_c = vfnmadd_fast(A11, vload(spu, (nk,n+(c-1))), C11_c)
    end
    C11vu = solve_AU(VecUnroll((Base.Cartesian.@ntuple $W C11)), spu, n, Val{$UNIT}())
    i = Unroll{2,1,$W,1,$W,0x0000000000000000,1}(Unroll{1,$W,$U,1,$W,0x0000000000000000,1}((StaticInt(0),n)))
    vstore!(spc, C11vu, i)
    maybestore!(spb, C11vu, i)
  end
end
@generated function rdiv_solve_W!(spc, spb, spa, spu, n, storec::B, mask::AbstractMask{W}, ::Val{UNIT}) where {W, UNIT, B}
  storecexpr = if (B <: Bool)
    :(storec && vstore!(spc, C11, i, mask))
  else
    :(vstore!(spc, C11, i, mask))
  end
  quote
    # $(Expr(:meta,:inline))
    # here, we just want to load the vectors
    C11 = VectorizationBase.data(vload(spa, Unroll{2,1,$W,1,$W,0xffffffffffffffff,1}((StaticInt(0),n)), mask))
    Base.Cartesian.@nexprs $W c -> C11_c = C11[c]
    for nk ∈ SafeCloseOpen(n) # nmuladd
      A11 = vload(spc, (MM{$W}(StaticInt(0)),nk), mask)
      Base.Cartesian.@nexprs $W c -> C11_c =  vfnmadd_fast(A11, vload(spu, (nk,n+(c-1))), C11_c)
    end
    C11 = VecUnroll((Base.Cartesian.@ntuple $W C11))
    C11 = solve_AU(C11, spu, n, Val{$UNIT}())
    i = Unroll{2,1,$W,1,$W,0xffffffffffffffff,1}((StaticInt(0),n))
    $storecexpr
    maybestore!(spb, C11, i, mask)
  end
end

function rdiv_U!(spc::AbstractStridedPointer{T}, spa, spu, M, N, ::StaticInt{1}, ::Val{UNIT}) where {T,UNIT}
  WS = pick_vector_width(T)
  W = Int(WS)
  UF = unroll_factor(WS)
  WU = UF*WS
  MU = UF > 1 ? M : 0
  Nd, Nr = VectorizationBase.vdivrem(N, WS)
  m = 0
  while m < MU - WU + 1
    n = Nr
    if n > 0
      BdivU_small_kern_u!(spc, nothing, spa, spu, n, UF, Val(UNIT))
    end
    for i ∈ 1:Nd
      rdiv_solve_W_u!(spc, nothing, spa, spu, n, WS, UF, Val(UNIT))
      n += W
    end
    m += WU
    spa = gesp(spa, (WU,StaticInt(0)))
    spc = gesp(spc, (WU,StaticInt(0)))
  end
  finalmask = VectorizationBase.mask(WS, M)
  while m < M
    ubm = m+W
    nomaskiter = ubm < M
    mask = nomaskiter ? VectorizationBase.max_mask(WS) : finalmask
    n = Nr
    if n > 0
      BdivU_small_kern!(spc, nothing, spa, spu, n, mask, Val(UNIT))
    end
    for _ ∈ 1:Nd
      rdiv_solve_W!(spc, nothing, spa, spu, n, nothing, mask, Val(UNIT))
      n += W
    end
    spa = gesp(spa, (WS,StaticInt(0)))
    spc = gesp(spc, (WS,StaticInt(0)))
    m = ubm
  end
  nothing
end

const LDIVBUFFERS = Vector{UInt8}[]
function lubuffer(::Val{T}, ::StaticInt{UF}, N) where {T, UF}
  buff = LDIVBUFFERS[Threads.threadid()]
  RSUF = StaticInt{UF}()*VectorizationBase.register_size()
  L = RSUF*N
  L > length(buff) && resize!(buff, L)
  ptr = Base.unsafe_convert(Ptr{T}, buff)
  StridedPointer{T,2,1,0,(1,2)}(ptr, (VectorizationBase.static_sizeof(T), RSUF), (StaticInt(0),StaticInt(0)))
end
function rdiv!(A::AbstractMatrix{T}, U::UpperTriangular{T}, ::Val{THREAD} = Val(true)) where {T<:Union{Float32,Float64},THREAD}
  div_dispatch!(A, A, parent(U), Val(false), Val(THREAD))
  return A
end
function rdiv!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, U::UpperTriangular{T}, ::Val{THREAD} = Val(true)) where {T<:Union{Float32,Float64},THREAD}
  div_dispatch!(C, A, parent(U), Val(false), Val(THREAD))
  return C
end
function rdiv!(A::AbstractMatrix{T}, U::UnitUpperTriangular{T}, ::Val{THREAD} = Val(true)) where {T<:Union{Float32,Float64},THREAD}
  div_dispatch!(A, A, parent(U), Val(true), Val(THREAD))
  return A
end
function rdiv!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, U::UnitUpperTriangular{T}, ::Val{THREAD} = Val(true)) where {T<:Union{Float32,Float64},THREAD}
  div_dispatch!(C, A, parent(U), Val(true), Val(THREAD))
  return C
end
function ldiv!(U::LowerTriangular{T}, A::AbstractMatrix{T}, ::Val{THREAD} = Val(true)) where {T<:Union{Float32,Float64},THREAD}
  div_dispatch!(transpose(A), transpose(A), transpose(parent(U)), Val(false), Val(THREAD))
  return A
end
function ldiv!(C::AbstractMatrix{T}, U::LowerTriangular{T}, A::AbstractMatrix{T}, ::Val{THREAD} = Val(true)) where {T<:Union{Float32,Float64},THREAD}
  div_dispatch!(transpose(C), transpose(A), transpose(parent(U)), Val(false), Val(THREAD))
  return C
end
function ldiv!(U::UnitLowerTriangular{T}, A::AbstractMatrix{T}, ::Val{THREAD} = Val(true)) where {T<:Union{Float32,Float64},THREAD}
  div_dispatch!(transpose(A), transpose(A), transpose(parent(U)), Val(true), Val(THREAD))
  return A
end
function ldiv!(C::AbstractMatrix{T}, U::UnitLowerTriangular{T}, A::AbstractMatrix{T}, ::Val{THREAD} = Val(true)) where {T<:Union{Float32,Float64},THREAD}
  div_dispatch!(transpose(C), transpose(A), transpose(parent(U)), Val(true), Val(THREAD))
  return C
end

ldiv!(A, B, ::Val = Val(true)) = LinearAlgebra.ldiv!(A, B)
ldiv!(Y, A, B, ::Val = Val(true)) = LinearAlgebra.ldiv!(Y, A, B)
rdiv!(A, B, ::Val = Val(true)) = LinearAlgebra.rdiv!(A, B)

function block_size(::Val{T}) where {T}
  elements_l2 = (VectorizationBase.cache_size(StaticInt(2))*StaticInt(19)) ÷ (VectorizationBase.static_sizeof(T)*StaticInt(60))
  Static.floortostaticint(sqrt(elements_l2))
end

function nmuladd!(C, A, U, M, K, N)
  @turbo for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
    Cmn = A[m,n]
    for k ∈ CloseOpen(K)
      Cmn -= C[m,k]*U[k,n]
    end
    C[m,K+n] = Cmn
  end
end

function rdiv_block_N!(
  spc::AbstractStridedPointer{T}, spa, spu, M, N, ::Val{UNIT}, ::StaticInt{X}, Bsize = nothing
) where {T,UNIT,X}
  spa_rdiv = spa
  spc_base = spc
  n = 0
  W = VectorizationBase.pick_vector_width(T)
  B_normalized = Bsize === nothing ? VectorizationBase.vcld(N, VectorizationBase.vcld(N, block_size(Val(T)))*W)*W : Bsize
  repeat = N > B_normalized
  N_temp = Core.ifelse(repeat, B_normalized, N)
  while true
    # println("Solve with N_temp = $N_temp and n = $n")
    rdiv_U!(spc, spa_rdiv, gesp(spu, (n,StaticInt{0}())), M, N_temp, StaticInt{X}(), Val(UNIT))
    repeat || break

    spa = gesp(spa, (StaticInt(0), B_normalized))
    spc = gesp(spc, (StaticInt(0), B_normalized))
    spu = gesp(spu, (StaticInt(0), B_normalized))
    nnext = n + B_normalized
    # N_temp =
    n += B_normalized
    repeat = n + B_normalized < N
    N_temp = repeat ? N_temp : N - n
    # N_temp = min(n + B_normalized, N) - n
    # println("nmuladd with N_temp = $N_temp and n = $n")
    # mul!(
    #   copyto!(view(C, :, n+1:n+N_temp), view(A, :, n+1:n+N_temp)),
    #   view(C, :, 1:n),
    #   view(U, 1:n, n+1:n+N_temp),
    #   -1.0, 1.0
    # )
    nmuladd!(spc_base, spa, spu, M, n, N_temp)
    spa_rdiv = spc
  end
end
function rdiv_block_MandN!(
  spc::AbstractStridedPointer{T}, spa, spu, M, N, ::Val{UNIT}, ::StaticInt{X}
) where {T,UNIT,X}
  B = block_size(Val(T))
  W = VectorizationBase.pick_vector_width(T)
  B_normalized = VectorizationBase.vcld(N, VectorizationBase.vcld(N, B)*W)*W
  WUF = W*unroll_factor(W)
  B_m = VectorizationBase.vcld(M, VectorizationBase.vcld(M, B)*WUF)*WUF
  m = 0
  while m < M
    mu = m + B_m
    Mtemp = min(M, mu) - m
    rdiv_block_N!(
      spc, spa, spu, Mtemp, N, Val(UNIT), StaticInt{X}(),
      VectorizationBase.vcld(N, VectorizationBase.vcld(N, B)*W)*W
    )
    spa = gesp(spa, (B_m, StaticInt{0}()))
    spc = gesp(spc, (B_m, StaticInt{0}()))
    m = mu
  end
end
function _nthreads()
  nc = VectorizationBase.num_cores()
  nt = VectorizationBase.num_threads()
  ifelse(Static.lt(nc,nt),nc,nt)
end
function m_thread_block_size(M, N, ::Val{T}) where {T}
  W = VectorizationBase.pick_vector_width(T)
  WUF = W * unroll_factor(W)
  nb = clamp(VectorizationBase.vdiv(M * N, StaticInt{256}() * W), 1, Int(_nthreads()))
  min(M, VectorizationBase.vcld(M, nb*W)*W)
end

function multithread_rdiv!(
  spc::AbstractStridedPointer{T}, spa, spu, M, N, mtb, ::Val{UNIT}, ::StaticInt{X}
) where {X,T,UNIT}
  (Md, Mr) = VectorizationBase.vdivrem(M, mtb)
  Nblock = Md + (Mr ≠ 0)
  Mrem = Core.ifelse(Mr ≠ 0, Mr, mtb)
  # @show mtb, Nblock, Mrem, Md, Mr
  # return
  let (Md, Mr) = VectorizationBase.vdivrem(M, mtb), Nblock = Md + (Mr ≠ 0), Mrem = Core.ifelse(Mr ≠ 0, Mr, mtb)
    @batch for block in CloseOpen(Nblock)
    # let block = 0
      Mtemp = Core.ifelse(block == Nblock-1, Mrem, mtb)
      rdiv_block_MandN!(
      # rdiv_block_N!(
        gesp(spc, (mtb*block, StaticInt{0}())),
        gesp(spa, (mtb*block, StaticInt{0}())),
        spu, Mtemp, N, Val{UNIT}(), StaticInt{X}()
        # spu, M, N, Val{UNIT}(), StaticInt{X}()
      )
    end
  end
  # nlaunch = Md - (Mr == 0)
  # threads, torelease = Polyester.request_threads(Base.Threads.threadid(), nlaunch)
  # nthread = length(threads)
  # if (nthread % Int32) ≤ zero(Int32)
  #   return rdiv_block_MandN!(spc, spa, spu, M, N, Val(UNIT), StaticInt{X}())
  # end
  # nbatch = nthread + one(nthread)

end

function div_dispatch!(C::AbstractMatrix{T}, A, U, ::Val{UNIT}, ::Val{THREAD}) where {UNIT,T,THREAD}
  M, N = size(A)
  ((N == 0) | (M == 0)) && return C
  _spa, spap = stridedpointer_preserve(A)
  _spc, spcp = stridedpointer_preserve(C)
  _spu, spup = stridedpointer_preserve(U)
  spa = zero_offsets(_spa)
  spc = zero_offsets(_spc)
  spu = zero_offsets(_spu)
  GC.@preserve spap spcp spup begin
    mtb = m_thread_block_size(M, N, Val(T))
    (THREAD & (M > mtb)) && return multithread_rdiv!(spc, spa, spu, M, N, mtb, Val(UNIT), VectorizationBase.contiguous_axis(A))
    if VectorizationBase.num_threads() == 1
      N > block_size(Val(T)) && return rdiv_block_MandN!(spc, spa, spu, M, N, Val(UNIT), VectorizationBase.contiguous_axis(A))
    end
    rdiv_U!(spc, spa, spu, M, N, VectorizationBase.contiguous_axis(A), Val(UNIT))
  end
end

# We're using `W x W` blocks, consuming `W` registers
# For each block we need to load 1 more value, plus another register is used for `B`. So:
# remaining_registers == register_count() - num_blocks * (W + 1) - 1
# 0 < register_count() - num_blocks * (W + 1) - 1
# num_blocks < (register_count() - 1) / (W + 1)
# num_blocks = (register_count() - 1) ÷ (W + 1)
function unroll_factor(::StaticInt{W}) where {W}
  num_blocks = (VectorizationBase.register_count() - StaticInt{1}()) ÷ (StaticInt{W}() + StaticInt{1}())
  ifelse(Static.lt(num_blocks, StaticInt{1}()), StaticInt{1}(), num_blocks)
end

function rdiv_U!(spc::AbstractStridedPointer{T}, spa, spu, M, N, ::StaticInt, ::Val{UNIT}) where {T,UNIT}
  WS = pick_vector_width(T)
  W = Int(WS)
  UF = unroll_factor(WS)
  WU = UF*WS
  MU = UF > 1 ? M : 0
  Nd, Nr = VectorizationBase.vdivrem(N, WS)
  spb = lubuffer(Val(T), UF, N)
  m = 0
  while m < MU - WU + 1
    n = Nr
    if n > 0
      BdivU_small_kern_u!(spb, spc, spa, spu, n, UF, Val(UNIT))
    end
    for i ∈ 1:Nd
      rdiv_solve_W_u!(spb, spc, spa, spu, n, WS, UF, Val(UNIT))
      n += W
    end
    m += WU
    spa = gesp(spa, (WU,StaticInt(0)))
    spc = gesp(spc, (WU,StaticInt(0)))
  end
  finalmask = VectorizationBase.mask(WS, M)
  while m < M
    ubm = m+W
    nomaskiter = ubm < M
    mask = nomaskiter ? VectorizationBase.max_mask(WS) : finalmask
    n = Nr
    if n > 0
      BdivU_small_kern!(spb, spc, spa, spu, n, mask, Val(UNIT))
    end
    for i ∈ 1:Nd
      # @show C, n
      rdiv_solve_W!(spb, spc, spa, spu, n, i ≠ Nd, mask, Val(UNIT))
      n += W
    end
    spa = gesp(spa, (WS,StaticInt(0)))
    spc = gesp(spc, (WS,StaticInt(0)))
    m = ubm
  end
  nothing
end

function __init__()
  nthread = Threads.nthreads()
  resize!(LDIVBUFFERS, nthread)
  for i ∈ 1:nthread
    LDIVBUFFERS[i] = Vector{UInt8}(undef, 3VectorizationBase.register_size()*128)
  end
end

end


text:
module DICE

#Don't polute the namespace. We want to use `solve` ourselves.
import JuMP;
import JuMP: @variable, @constraint, @NLconstraint, @objective, @NLparameter;
import JuMP: getvalue, setvalue, getdual, setupperbound, setlowerbound;
#Ipopt is distributed under the EPL.
#We don't package it here though, and assume you have this set up on your system already.
using Ipopt;

include("Abstractions.jl")

#DICE Versions
abstract type Version end

#Configuration options & parameters, model settings and output
include("BaseTypes.jl")

#Scenarios
abstract type Scenario end
include("Scenarios.jl")

struct DICENarrative
    constants::Options
    parameters::Parameters
    model::JuMP.Model
    scenario::Scenario
    version::Version
    variables::Variables
    equations::Equations
    results::Results
end

function Base.show(io::IO, ::MIME"text/plain", dice::DICENarrative)
    println(io, "$(dice.scenario) scenario using $(dice.version).")
    show(io, dice.model)
end

function solve end
function options end

export solve, options

# Include all version implementations
include("2013R.jl")
include("2016R.jl")
include("CJL.jl")

end

ext: jl


text:
using Distributed

struct DistributedSum{T}
  data::Vector{T}
end

# Construct from a distributed sum from a
# function that returns the elements of the sum.
function DistributedSum(f::Function, n::Integer)
  # Assign the matrices to the first `n` workers.
  return DistributedSum([@spawnat(workers()[ni], f(ni)) for ni in 1:n])
end

# Functions needed for @distributed loop and reduction
Base.length(A::DistributedSum) = length(A.data)
Base.firstindex(A::DistributedSum) = firstindex(A.data)
Base.lastindex(A::DistributedSum) = lastindex(A.data)
Base.getindex(A::DistributedSum, args...) = getindex(A.data, args...)
Base.iterate(A::DistributedSum, args...) = iterate(A.data, args...)

# Apply the sum
function (A::DistributedSum)(v)
  return @distributed (+) for An in A
    fetch(An)(v)
  end
end

max_stars_repo_name: mtfishman/ITensorParallel


text:
##########################################################################
# The TriangularMacroNetwork is a macrocell network with three sites,
# placed on the vertices of an equilateral triangle. The simulation
# parameters are based on 3GPP Case 1 (TR 25.814 and TR 36.814) macrocell
# simulation environment. Distance-dependent pathloss and log-normal
# shadow fading is used. The shadow fading is correlated (with correlation
# coefficient 0.5) from each MS to the three BSs, but independent
# between MSs.
#
# The typical transmit power is 46 dBm over 10 MHz (600 subcarriers).
# This leads to 10*log10(10^(46/10)/600) = 18.2 dBm for one 15 kHz subcarrier.

##########################################################################
# Six sector antenna params from 3gpp TR 25.996
immutable SixSector3gppAntennaParams <: AntennaParams
    bore_sight_angle::Float64
    angle_with_3dB_gain::Float64
    min_antenna_gain_dB::Float64
end

get_antenna_gain(antenna_params::SixSector3gppAntennaParams, angle) =
    10^(-min(12*(angle/antenna_params.angle_with_3dB_gain)^2, antenna_params.min_antenna_gain_dB)/10)

get_angle{AntennaParams_t <: SixSector3gppAntennaParams}(MS, BS::PhysicalBS{AntennaParams_t}) =
    get_angle(MS.position, BS.position) - BS.antenna_params.bore_sight_angle

##########################################################################
# Network definition
type TriangularMacroNetwork{MS_t <: PhysicalMS, BS_t <: PhysicalBS, System_t <: System, PropagationEnvironment_t <: PropagationEnvironment} <: PhysicalNetwork
    MSs::Vector{MS_t}
    BSs::Vector{BS_t}

    system::System_t
    num_MSs_per_cell::Int
    propagation_environment::PropagationEnvironment_t
    inter_site_distance::Float64
    guard_distance::Float64
    aux_network_params::AuxNetworkParams

    assignment::Assignment
end

# Convenience constructor without network params and assignments
TriangularMacroNetwork(MSs, BSs, system, num_MSs_per_cell, propagation_environment, inter_site_distance, guard_distance) =
    TriangularMacroNetwork(MSs, BSs, system, num_MSs_per_cell, propagation_environment, inter_site_distance, guard_distance, AuxNetworkParams(), Assignment())

Base.show(io::IO, x::TriangularMacroNetwork) =
    print(io, "TriangularMacro(I = $(length(x.BSs)), Kc = $(x.num_MSs_per_cell), ISD = $(x.inter_site_distance), GD = $(x.guard_distance))")
Base.showcompact(io::IO, x::TriangularMacroNetwork) =
    print(io, "TriangularMacro($(length(x.BSs)), $(x.num_MSs_per_cell), $(x.inter_site_distance), $(x.guard_distance))")

# The default parameter values are taken from 3GPP Case 1
# (TR 25.814 and TR 36.814).
function setup_triangularmacro_network(
    num_MSs_per_cell, num_MS_antennas, num_BS_antennas;
    system = SinglecarrierSystem(2e9, 15e3),
    propagation_environment = SimpleLargescaleFadingEnvironment(37.6, 15.3, 20, 8),
    inter_site_distance = 500.,
    guard_distance = 35.,
    transmit_power = 10^(18.2/10), transmit_powers = transmit_power*ones(Float64, 3),
    BS_antenna_gain_params =
        [SixSector3gppAntennaParams(-90/180*pi, 35/180*pi, 23),
         SixSector3gppAntennaParams( 30/180*pi, 35/180*pi, 23),
         SixSector3gppAntennaParams(150/180*pi, 35/180*pi, 23)],
    user_priority = 1., user_priorities = user_priority*ones(Float64, 3*num_MSs_per_cell),
    num_streams = 1, num_streamss = num_streams*ones(Int, 3*num_MSs_per_cell),
    MS_antenna_gain_dB = 0., MS_antenna_gains_dB = MS_antenna_gain_dB*ones(Float64, 3*num_MSs_per_cell),
    receiver_noise_figure = 9., receiver_noise_figures = receiver_noise_figure*ones(Float64, 3*num_MSs_per_cell))

    isa(num_MS_antennas, Vector) || (num_MS_antennas = num_MS_antennas*ones(Int, 3*num_MSs_per_cell))
    isa(num_BS_antennas, Vector) || (num_BS_antennas = num_BS_antennas*ones(Int, 3))

    BSs = [
        PhysicalBS(num_BS_antennas[1],
            Position(0, sqrt((inter_site_distance/2)^2 + (inter_site_distance/(2*sqrt(3)))^2)),
            transmit_powers[1], BS_antenna_gain_params[1]),
        PhysicalBS(num_BS_antennas[2],
            Position(-inter_site_distance/2, -inter_site_distance/(2*sqrt(3))),
            transmit_powers[2], BS_antenna_gain_params[2]),
        PhysicalBS(num_BS_antennas[3],
            Position(+inter_site_distance/2, -inter_site_distance/(2*sqrt(3))),
            transmit_powers[3], BS_antenna_gain_params[3])
    ]
    MSs = [ PhysicalMS(num_MS_antennas[k], Position(0, 0), Velocity(0, 0), user_priorities[k], num_streamss[k], MS_antenna_gains_dB[k], receiver_noise_figures[k], SimpleLargescaleFadingEnvironmentState(zeros(Float64, 3), falses(3))) for k = 1:3*num_MSs_per_cell ]

    TriangularMacroNetwork(MSs, BSs, system, num_MSs_per_cell,
        propagation_environment, inter_site_distance, guard_distance)
end

##########################################################################
# Cell assignment
function IDCellAssignment!(channel, network::TriangularMacroNetwork)
    Kc = network.num_MSs_per_cell; I = get_num_BSs(network)
    cell_assignment = Array(Int, I*Kc)

    for i = 1:I
        cell_assignment[(i-1)*Kc+1:i*Kc] = i
    end

    network.assignment = Assignment(cell_assignment, I)

    return AssignmentResults()
end

# Do this for now, even though the shadowing might actually give a certain MS
# better connection to a BS in another cell.
LargeScaleFadingCellAssignment!(channel, network::TriangularMacroNetwork) =
    IDCellAssignment!(channel, network)

##########################################################################
# Simulation functions
function draw_user_drop!{MS_t <: PhysicalMS, BS_t <: PhysicalBS, System_t <: System}(network::TriangularMacroNetwork{MS_t, BS_t, System_t, SimpleLargescaleFadingEnvironment})
    I = get_num_BSs(network); K = get_num_MSs(network)

    # Shadow fading covariance (correlation coefficient 0.5 between cells)
    Cov_shadow_sqrtm = chol(network.propagation_environment.shadow_sigma_dB^2*(ones(3,3) + eye(3))/2)

    for k = 1:K
        # Generate user position within standard triangle [0,30] degrees, to
        # the right of the base. Apply guard distance to x coordinate.
        xtri = sqrt(((network.inter_site_distance/2)^2 - network.guard_distance^2)*rand() + network.guard_distance^2)
        ytri = (xtri/sqrt(3))*rand()

        # Flip it over with probability 0.5
        if rand() < 0.5
            pos = [xtri;ytri]
        else
            pos = [xtri;-ytri]

            theta = 60*pi/180
            pos = [cos(theta) -sin(theta);sin(theta) cos(theta)]*pos
        end

        # Generate rotation angle
        i = div(k - 1, network.num_MSs_per_cell) + 1 # serving BS id
        if i == 1
            theta = 240*pi/180
        elseif i == 2
            theta = 0.
        elseif i == 3
            theta = 120*pi/180
        end

        # Get final user position
        position = [cos(theta) -sin(theta);sin(theta) cos(theta)]*pos
        network.MSs[k].position = (Position(position[1], position[2])
                                   + network.BSs[i].position)

        # Shadow fading
        network.MSs[k].propagation_environment_state =
            SimpleLargescaleFadingEnvironmentState(Cov_shadow_sqrtm*randn(I), falses(I)) # LoS not used in this model
    end
end

function draw_channel{MS_t <: PhysicalMS, BS_t <: PhysicalBS}(network::TriangularMacroNetwork{MS_t, BS_t, SinglecarrierSystem, SimpleLargescaleFadingEnvironment})
    K = get_num_MSs(network); I = get_num_BSs(network)
    Ns = get_num_MS_antennas(network); Ms = get_num_BS_antennas(network)

    coefs = Array(Matrix{Complex128}, K, I)
    large_scale_fading_factor = Array(Float64, K, I)

    distances = get_distances(network)
    angles = get_angles(network)

    for k = 1:K
        for i = 1:I
            # Small scale fading
            coefs[k,i] = (1/sqrt(2))*(randn(Ns[k], Ms[i])
                                 + im*randn(Ns[k], Ms[i]))

            # Pathloss
            pathloss_factor = sqrt(10^(-(network.propagation_environment.pathloss_beta + network.propagation_environment.pathloss_alpha*log10(distances[k,i]))/10))

            # Shadow fading
            shadow_factor = sqrt(10^(network.MSs[k].propagation_environment_state.shadow_realizations_dB[i]/10))

            # Penetration loss
            penetration_loss_factor = sqrt(10^(-network.propagation_environment.penetration_loss_dB/10))

            # BS antenna gain
            bs_antenna_gain = sqrt(get_antenna_gain(network.BSs[i].antenna_params, angles[k,i]))

            # MS antenna gain
            ms_antenna_gain = sqrt(10^(network.MSs[k].antenna_gain_dB/10))

            # Apply scale factors
            large_scale_fading_factor[k,i] = pathloss_factor*shadow_factor*penetration_loss_factor*bs_antenna_gain*ms_antenna_gain
            coefs[k,i] *= large_scale_fading_factor[k,i]
        end
    end

    return SinglecarrierChannel(coefs, Ns, Ms, K, I, large_scale_fading_factor)
end

##########################################################################
# Visualization functions
function plot_network_layout(network::TriangularMacroNetwork)
    I = get_num_BSs(network); K = get_num_MSs(network)

    fig = PyPlot.figure()
    ax = fig[:add_subplot](1, 1, 1)

    # Equilateral triangle, with cell edges
    cm = mean([ [network.BSs[i].position.x, network.BSs[i].position.y] for i = 1:I ])
    for i = 1:I
        pos1 = network.BSs[i].position
        pos2 = network.BSs[mod(i, I) + 1].position
        ax[:plot]([pos1.x, pos2.x], [pos1.y, pos2.y]; color="b", linestyle="-")
        ax[:plot]([mean([pos1.x, pos2.x]), cm[1]], [mean([pos1.y, pos2.y]), cm[2]]; color="b", linestyle="--")
    end

    # BSs
    for i = 1:I
        pos = network.BSs[i].position
        ax[:plot](pos.x, pos.y; marker="x", color="b", markersize=10)
    end

    # MSs
    for k = 1:K
        pos = network.MSs[k].position
        ax[:plot](pos.x, pos.y; marker="o", color="r", markersize=10)
    end

    display(fig)
end


text:
module ITensorsVisualization

using GLMakie
using ITensors
using LinearAlgebra
using NetworkLayout
using SparseArrays

export @visualize

include("utils.jl")
include("visualize_tensornetwork.jl")
include("visualize_macro.jl")

end

max_stars_repo_name: mtfishman/ITensorsVisualization

max_issues_repo_name: mtfishman/ITensorsVisualization

max_forks_repo_name: mtfishman/ITensorsVisualization


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -33.59531116254298 <= q <= 31.312008728451346)

                     add_NL_constraint(m, :(softplus(0.7464435003747418 + -0.7676428720221384*softplus(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + -0.5937285886485864*softplus(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + 0.70936536775501*softplus(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.9647291374973435*softplus(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.527810271934432*softplus(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + softplus(-0.35636357579153266 + 0.8665784119884141*softplus(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + 0.03175667478129274*softplus(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + -0.7799384758999897*softplus(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.4104070349373172*softplus(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + -0.7416668360909755*softplus(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + softplus(0.5038932713674513 + 0.7223377417739489*softplus(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + -0.17110213406275054*softplus(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + 0.7726394930743403*softplus(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.6118790711423934*softplus(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.16250007841990932*softplus(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + softplus(-0.7190092552996492 + 0.6952008118484776*softplus(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + 0.975536624348428*softplus(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + 0.03367235210941821*softplus(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.8059867473566515*softplus(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.5874603423029647*softplus(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + softplus(-0.5715134646750415 + 0.6325466436773857*softplus(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + -0.2633742432221431*softplus(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + -0.3786600783719858*softplus(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.5923927477083897*softplus(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.9114166052727355*softplus(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Distributions, Plots; pyplot()

function anovaFStat(allData)
    xBarArray = mean.(allData)
    nArray = length.(allData)
    xBarTotal = mean(vcat(allData...))
    L = length(nArray)

    ssBetween = sum( [nArray[i]*(xBarArray[i] - xBarTotal)^2 for i in 1:L] )
    ssWithin = sum([sum([(ob - xBarArray[i])^2 for ob in allData[i]])
				for i in 1:L])
    return (ssBetween/(L-1))/(ssWithin/(sum(nArray)-L))
end

case1 = [13.4, 13.4, 13.4, 13.4, 13.4]
case2 = [12.7, 11.8, 13.4, 12.7, 12.9]
stdDevs = [2, 2, 2, 2, 2]
numObs = [24, 15, 13, 23, 9]
L = length(case1)

N = 10^5

mcFstatsH0 = Array{Float64}(undef, N)
for i in 1:N
    mcFstatsH0[i] = anovaFStat([ rand(Normal(case1[j],stdDevs[j]),numObs[j])
		for j in 1:L ])
end

mcFstatsH1 = Array{Float64}(undef, N)
for i in 1:N
    mcFstatsH1[i] = anovaFStat([ rand(Normal(case2[j],stdDevs[j]),numObs[j])
		for j in 1:L ])
end

stephist(mcFstatsH0, bins=100,
	c=:blue, normed=true, label="Equal group means case")
stephist!(mcFstatsH1, bins=100,
	c=:red, normed=true, label="Unequal group means case")

dfBetween = L - 1
dfError = sum(numObs) - 1
xGrid = 0:0.01:10
plot!(xGrid, pdf.(FDist(dfBetween, dfError),xGrid),
            c=:black, label="F-statistic analytic")
critVal = quantile(FDist(dfBetween, dfError),0.95)
plot!([critVal, critVal],[0,0.8],
	c=:black, ls=:dash, label="Critical value boundary",
	xlims=(0,10), ylims=(0,0.8), xlabel="F-value", ylabel="Density")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
export dropunused, dropunused!
import DataFrames
import SummarizedExperiments

"""
    dropunused!(x; samples = true, experiments = true, colnames = true, mapping = true)

Drop unused samples, experiments and/or column names from the `MultiAssayExperiment` `x`.
A reference to the modified `x` is returned.

If `samples = true`, `sampledata(x)` is filtered to only retain samples that are present in the sample mapping.

If `experiments = true`, `experiments(x)` is filtered to only retain experiments that are present in the sample mapping.

If `colnames = true`, each entry of `experiments(x)` is filtered to only retain column names that are present in the sample mapping for that experiment.

If `mapping = true`, the sample mapping is filtered to remove rows that contain samples, experiments or column names that do not exist in `x`.

# Examples
```jldoctest
julia> using MultiAssayExperiments

julia> x = exampleobject();

julia> filtersamplemap!(x; experiments = "bar"); # Only keeping experiment 'bar'

julia> dropunused!(x) # We see that 'foo' is dropped
MultiAssayExperiment object
  experiments(1): bar
  sampledata(2): name disease
  metadata(1): version
```
"""
function dropunused!(x::MultiAssayExperiment; samples = true, experiments = true, colnames = true, mapping = true)
    sm = samplemap(x)

    # Filtering the sample mapping first.
    if mapping
        sampset = Set(sm[!,"sample"])
        allcolnames = Dict{AbstractString,Set{<:String}}()
        for (k, v) in x.experiments
            allcolnames[k] = Set(SummarizedExperiments.coldata(v)[!,"name"])
        end

        filterfun = function(row)
            return (
                row.sample in sampset &&
                haskey(allcolnames, row.experiment) &&
                row.colname in allcolnames[row.experiment]
            )
        end
        sm = DataFrames.filter(filterfun, sm)
        x.samplemap = sm
    end

    # Applying filters on the samples.
    if samples
        sampset = Set(sm[!,"sample"])
        x.sampledata = DataFrames.filter(row -> row.name in sampset, sampledata(x))
    end

    # Applying filters on the experiments.
    ecopy = nothing
    if experiments
        expset = Set(sm[!,"experiment"])
        ecopy = copy(x.experiments) # making a copy so that we can use this function on a shallow copy.
        for e in keys(x.experiments)
            if !(e in expset)
                delete!(ecopy, e)
            end
        end
        x.experiments = ecopy
    end

    # Applying filters on the column names.
    if colnames
        valid_colnames = Dict{AbstractString,Set{<:AbstractString}}()
        for row in eachrow(sm)
            if haskey(valid_colnames, row.experiment)
                push!(valid_colnames[row.experiment], row.colname)
            else
                valid_colnames[row.experiment] = Set([row.colname])
            end
        end

        # If we didn't already make a copy, we do so. Again, this is
        # just to avoid modifying x.experiments in place in case
        # 'x' was a shallow copy of some not-to-be-modified MAE.
        if ecopy == nothing
            x.experiments = copy(x.experiments)
        end

        for (k, v) in x.experiments
            names = SummarizedExperiments.coldata(v)[!,"name"]
            valid = valid_colnames[k]

            keep = Vector{Bool}(undef, length(names))
            counter = 0
            for i in 1:length(names)
                keep[i] = names[i] in valid
                counter += Integer(keep[i])
            end

            # Avoid creating a copy unless it is strictly necessary.
            if counter != length(keep)
                x.experiments[k] = v[:,keep]
            end
        end
    end

    return x
end

"""
    dropunused(x; kwargs...)

Return a new `MultiAssayExperiment` where unused samples, experiments or column names are removed.
This makes a copy of `x` and passes it (and any keyword arguments in `kwargs`) to [`dropunused!`](@ref);
see the latter function for more details.

# Examples
```jldoctest
julia> using MultiAssayExperiments

julia> x = exampleobject();

julia> y = filtersamplemap(x; experiments = "bar"); # Only keeping experiment 'bar'

julia> dropunused(y) # We see that 'foo' is dropped
MultiAssayExperiment object
  experiments(1): bar
  sampledata(2): name disease
  metadata(1): version
```
"""
function dropunused(x::MultiAssayExperiment; kwargs...)
    y = copy(x)
    dropunused!(y; kwargs...)
    return y
end


text:
_stresstest_cpu_size_default(dtype) = floor(Int, sqrt(L2_cachesize() / sizeof(dtype)))

"""
    stresstest_cpu(core_or_cores)

Run a CPU stress test (matrix multiplication) on one or multiple CPU cores, as specified by the positional argument.
If no argument is provided (only) the currently active CPU core will be used.

**Keyword arguments:**

* `duration`: stress test will take about the given time in seconds.
* `dtype` (default: `Float64`): element type of the matrices
* `size` (default: `floor(Int, sqrt(L2_cachesize() / sizeof(dtype)))`): matrices of size `(size, size)` will be used
* `verbose` (default: `true`): toggle printing of information
* `parallel` (default: `true`): If `true`, will (try to) run each CPU core test on a different Julia thread. Make sure to have enough Julia threads.
* `threads` (default: `nothing`): If `parallel == true`, this argument may be used to specify the Julia threads to use.
"""
function stresstest_cpu(
    cores;
    duration=2,
    dtype=Float64,
    size=_stresstest_cpu_size_default(dtype),
    verbose=true,
    kwargs...,
)
    if !all(c -> 0 ≤ c < Sys.CPU_THREADS, cores)
        throw(ArgumentError("CPU core IDs must all be ≥ 0 and ≤ Sys.CPU_THREADS."))
    end
    if BLAS.get_num_threads() != 1 && verbose
        @info("Setting BLAS.set_num_threads(1)!")
        BLAS.set_num_threads(1)
    end
    verbose && @info("Will run for about $(duration) seconds on each CPU core!")
    Δt = @elapsed _run_stresstests_cpu(cores; verbose, duration, size, kwargs...)
    verbose && @info("Took $(round(Δt; digits=2)) seconds to run the tests.")
    return nothing
end

function _stresstest_cpu_default_threads(ncores)
    nthreads = Threads.nthreads()
    if nthreads >= ncores + 1
        # if there are enough "worker" threads
        # use them and skip the main thread (1).
        return collect(2:(ncores + 1))
    elseif nthreads == ncores
        return collect(1:ncores)
    elseif nthreads < ncores
        # not enough threads
        @warn("There are too few Julia threads to run all tests in parallel!")
        return [mod1(i, nthreads) for i in 1:ncores] # round-robin distribution
    end
end

function _run_stresstests_cpu(
    cores;
    parallel=true,
    threads=_stresstest_cpu_default_threads(length(cores)),
    verbose=true,
    kwargs...,
)
    if parallel == true
        if length(threads) != length(cores)
            throw(ArgumentError("length(threads) != length(cores)"))
        end
        @sync for (i, core) in enumerate(cores)
            if verbose
                @info("Julia thread $(threads[i]) runs test on CPU core $(core).")
            end
            @tspawnat threads[i] begin
                core_before = getcpuid()
                pinthread(core)
                stresstest_cpu_kernel(; verbose, kwargs...)
                pinthread(core_before)
            end
        end
    else
        for core in cores
            core_before = getcpuid()
            pinthread(core)
            stresstest_cpu_kernel(; verbose, kwargs...)
            pinthread(core_before)
        end
    end
    return nothing
end
stresstest_cpu(core::Integer=getcpuid(); kwargs...) = stresstest_cpu([core]; kwargs...)

function stresstest_cpu_kernel(;
    duration, dtype=Float64, size=_stresstest_cpu_size_default(dtype), verbose=true
)
    A = rand(dtype, size, size)
    B = rand(dtype, size, size)
    C = zeros(dtype, size, size)

    t_start = time()
    while (time() - t_start) < duration
        mul!(C, A, B)
    end
    return nothing
end


text:
export Task_TestTask

mutable struct Task_TestTask <: Task
end

function init_task(myt::Task_TestTask,rhd::RHD2000,han,fpga)
    setTtlOut(rhd,ones(Int32,16))
    nothing
end

function do_task(myt::Task_TestTask,rhd::RHD2000,myread,han,fpga)

    for i=1:SAMPLES_PER_DATA_BLOCK
        if rhd.adc[i,1]>10000
            setTtlOut(rhd,zeros(Int32,16))
            break
        end
    end

    nothing
end

function save_task(myt::Task_TestTask,rhd::RHD2000)

    f=open("adc1.bin","a+")
    for i=1:SAMPLES_PER_DATA_BLOCK
        write(f,rhd.adc[i,1])
    end
    close(f)
    nothing
end


text:
using LinearAlgebra

using GLMakie; GLMakie.activate!()
using AbstractPlotting

using Delaunay
using DEC2D

# convert vector of vectors into a matrix
vv2m(vv) = reduce(vcat, transpose.(vv))

# Simples case
begin
    vertices = vcat([[0, 0]], [[cos(2π * i / 6), sin(2π * i / 6)] for i = 1:6])
    edges = [
        [1, 2],
        [1, 3],
        [1, 4],
        [1, 5],
        [1, 6],
        [1, 7],
        [2, 3],
        [3, 4],
        [4, 5],
        [5, 6],
        [6, 7],
        [7, 2],
    ]
    triangles = [[1, 2, 3], [1, 3, 4], [1, 4, 5], [1, 5, 6], [1, 6, 7], [1, 7, 2]]

    # return the mesh and the indices of the interior vertices
    vint, mesh = DECMesh(vertices, edges, triangles)
    vext = setdiff(1:length(vertices), vint)

    # d₀ = boundary1(mesh)
    # d₁ = boundary2d(mesh)
    # *₁ = hodgestar1(mesh)
    # *₂ = hodgestar2(mesh)
    # Δ = *₂ * d₁ * *₁ * d₀
    Δ = laplacian(mesh) # calculate the laplacian operator

    b = [-1 - 10 * sum(Δ[i, vext]) for i = 1:length(vint)]

    # numerical solution
    zf = Δ[vint, vint] \ b

    x = vv2m(vertices)[:, 1]
    y = vv2m(vertices)[:, 2]
    zn = similar(vertices, Float64)
    zn[vint] = zf
    zn[vext] .= 10

    # exact solution
    f(x, y) = (1 - x^2 - y^2) / 4 + 10
    z = [f(v...) for v in vertices]

    println(
        "La norma de la diferencia entre la solución numérica y la exacta es ",
        norm(zn - z),
    )

    #   ____  _       _   _   _
    #  |  _ \| | ___ | |_| |_(_)_ __   __ _
    #  | |_) | |/ _ \| __| __| | '_ \ / _` |
    #  |  __/| | (_) | |_| |_| | | | | (_| |
    #  |_|   |_|\___/ \__|\__|_|_| |_|\__, |
    #                                 |___/

    # plot numerical solution
    scene = Scene()
    mesh!(
        scene,
        hcat(x, y, zn),
        vv2m(mesh.tri.triangles),
        color = zn,
        scale_plot = false,
    )
    xlims!(scene, -1, 1)
    ylims!(scene, -1, 1)
    zlims!(scene, 10, 10.25)
    wireframe!(scene[end][1], color = (:black, 0.6), linewidth = 2)
    display(scene)
    sleep(5.0)
end

# Caso medio
begin
    # we get the vertices of a more fine mesh
    vertices = [
        [cos(2π * i / (50 * r + 1)), sin(2π * i / (50 * r + 1))] .* r for r = 0:0.1:1
        for i = 1:round(Int, 50 * r + 1, RoundDown)
    ]

    # get thetriangulation using Delaunay algorithm
    meshd = delaunay(vv2m(vertices))
    triangles = [meshd.simplices[i, :] for i = 1:size(meshd.simplices, 1)]
    edges = reduce(
        union,
        [[sort([t[1], t[2]]), sort([t[2], t[3]]), sort([t[3], t[1]])] for t in triangles],
    )

    vint, mesh = DECMesh(vertices, edges, triangles)
    vext = setdiff(1:length(vertices), vint)

    Δ = laplacian(mesh)

    b = [-1 - 10 * sum(Δ[i, vext]) for i = 1:length(vint)]
    zf = Δ[vint, vint] \ b

    x = vv2m(vertices)[:, 1]
    y = vv2m(vertices)[:, 2]
    zn = similar(vertices, Float64)
    zn[vint] = zf
    zn[vext] .= 10

    f(x, y) = (1 - x^2 - y^2) / 4 + 10
    z = [f(v...) for v in vertices]

    println(
        "La norma de la diferencia entre la solución numérica y la exacta es ",
        norm(zn - z),
    )

    scene = Scene()
    mesh!(
        scene,
        hcat(x, y, zn),
        vv2m(mesh.tri.triangles),
        color = zn,
        scale_plot = false,
    )
    xlims!(scene, -1, 1)
    ylims!(scene, -1, 1)
    zlims!(scene, 10, 10.25)
    wireframe!(scene[end][1], color = (:black, 0.6), linewidth = 2)
    display(scene)
    sleep(5.0)
end

# Caso denso
begin
    # we get the vertices of a more fine mesh
    vertices = [
        [cos(2π * i / (200 * r + 1)), sin(2π * i / (200 * r + 1))] .* r for r = 0:0.025:1 for i = 1:round(Int, 200 * r + 1, RoundDown)
    ]

    # get thetriangulation using Delaunay algorithm
    meshd = delaunay(vv2m(vertices))
    triangles = [meshd.simplices[i, :] for i = 1:size(meshd.simplices, 1)]
    edges = reduce(
        union,
        [[sort([t[1], t[2]]), sort([t[2], t[3]]), sort([t[3], t[1]])] for t in triangles],
    )

    vint, mesh = DECMesh(vertices, edges, triangles)
    vext = setdiff(1:length(vertices), vint)

    Δ = laplacian(mesh)

    b = [-1 - 10 * sum(Δ[i, vext]) for i = 1:length(vint)]
    zf = Δ[vint, vint] \ b

    x = vv2m(vertices)[:, 1]
    y = vv2m(vertices)[:, 2]
    zn = similar(vertices, Float64)
    zn[vint] = zf
    zn[vext] .= 10

    f(x, y) = (1 - x^2 - y^2) / 4 + 10
    z = [f(v...) for v in vertices]

    println(
        "La norma de la diferencia entre la solución numérica y la exacta es ",
        norm(zn - z),
    )

    scene = Scene()
    mesh!(
        scene,
        hcat(x, y, zn),
        vv2m(mesh.tri.triangles),
        color = zn,
        scale_plot = false,
    )
    xlims!(scene, -1, 1)
    ylims!(scene, -1, 1)
    zlims!(scene, 10, 10.25)
    wireframe!(scene[end][1], color = (:black, 0.6), linewidth = 2)
    display(scene)
end


text:
module FHLampWire

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/WireLamps" WireLamps(x::Integer, y::Integer, above::Bool=false, colors::String="Red,Yellow,Blue,Green,Orange", wireColor::String="595866",
                                                 lightCount::Integer=3, lightAlpha::Number=1.0, lightStartFade::Integer=8, lightEndFade::Integer=16)

const placements = Ahorn.PlacementDict(
    "Wire Lamps (Frost Helper)" => Ahorn.EntityPlacement(
        WireLamps,
        "line"
    )
)

Ahorn.nodeLimits(entity::WireLamps) = 1, 1

function Ahorn.selection(entity::WireLamps)
    nodes = get(entity.data, "nodes", ())
    x, y = Ahorn.position(entity)

    res = Ahorn.Rectangle[Ahorn.Rectangle(x - 4, y - 4, 8, 8)]

    for node in nodes
        nx, ny = node

        push!(res, Ahorn.Rectangle(nx - 4, ny - 4, 8, 8))
    end

    return res
end

wireColorHex = "595866"
wireColor = (89, 88, 102, 1) ./ (255, 255, 255, 1)
wireColorSelected = (Ahorn.colors.selection_selected_fc)

function renderWire(ctx::Ahorn.Cairo.CairoContext, entity::WireLamps, color::Ahorn.colorTupleType=wireColor)
    x, y = Ahorn.position(entity)

    start = (x, y)
    stop = get(entity.data, "nodes", [start])[1]
    control = (start .+ stop) ./ 2 .+ (0, 24)

    curve = Ahorn.SimpleCurve(start, stop, control)
    Ahorn.drawSimpleCurve(ctx, curve, color, thickness=1)
end

Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::WireLamps) = renderWire(ctx, entity, wireColorSelected)

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::WireLamps, room::Maple.Room)
    # Make sure Alpha is 1
    rawColor = Ahorn.argb32ToRGBATuple(parse(Int, wireColorHex, base=16))[1:3] ./ 255
    color = (rawColor..., 1.0)

    renderWire(ctx, entity, color)
end

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
include("module.jl")
using ApolloniusProblems

let test = [Circle(0.0, 0.0, 1.0), Circle(4.0, 0.0, 1.0), Circle(2.0, 4.0, 2.0)]
    println("The defining circles are: \n - ", join(test, "\n - "))
    println("The internal circle is:\n\t", ApolloniusProblems.solve(test))
    println("The external circle is:\n\t", ApolloniusProblems.solve(test, 1:3))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module SpikeTrains

export SpikeTrain, draw_uncorrelated_spikes, draw_correlated_spikes, draw_coincident_spikes, length, iterate, convert, vcat, merge, make_exponentialShift, correlation_code, coincidence_code, plot_spike_trains

using Distributions, Plots#, PlotRecipes

struct SpikeTrain
    times::Vector{Float64}
    SpikeTrain(times=Float64[]) = new(sort(times))
end

Base.length(s::SpikeTrain) = length(s.times)
Base.iterate(s::SpikeTrain, args...) = iterate(s.times, args...)
Base.convert(SpikeTrain, v) = SpikeTrain(v)
Base.vcat(s::SpikeTrain...) = SpikeTrain(vcat(getfield.(s,:times)...))
Base.merge(s::Array{SpikeTrain}; dims=Base.OneTo(ndims(s))) = dropdims(mapslices(x->vcat(x...), s, dims=dims), dims=dims)

function has_min_distance(times, min_distance)
    keep = ones(Bool, length(times))
    last_time = -Inf
    for (i,t) ∈ enumerate(times)
        if t < last_time
            keep[i] = false
        else
            last_time = t+min_distance
        end
    end
    return keep
end

ensure_min_distance(times, min_distance) = times[has_min_distance(times, min_distance)]

"""
    draw_uncorrelated_spikes(trange, rates; sorted=true)

Draw uncorrelated spike trains in the interval `trange` with the respective `rates`.
"""
function draw_uncorrelated_spikes(trange, rates)
    duration = trange[2]-trange[1]
    source = Vector{SpikeTrain}(undef, length(rates))
    # draw uncorrelated spike trains with high rate
    for (i,r) ∈ enumerate(rates)
        num_spikes = rand(Poisson(duration*r))
        spikes = rand(num_spikes).*duration.+trange[1]
        sort!(spikes)

        source[i] = SpikeTrain(spikes)
    end

    return source
end

make_exponentialShift(τ) = t -> t+rand(Exponential(τ))

"""
    draw_correlated_spikes(trange, rates, c, [shift]; min_master_distance=0)

Draws spike trains with given `rates` within the given time interval `trange` and
pairwise correlation coefficients c. A function to `shift` the drawn spikes can be
included. A minimum distance between master spikes can be required. The algorithm implemented
is the offline version of the mixture process (4.6.1) by Brette 2008 [1].

[1](http://romainbrette.fr/WordPress3/wp-content/uploads/2014/06/Brette2008NC.pdf)

E[wt] = E[n_good]*min_master_distance
P(n_good|E[wt]) = Poisson(n_good|rate*(1-E[n_good]*min_master_distance))
E[n_good] = rate*1-rate*E[n_good]*min_master_distance
E[n_good] = rate/(1+rate*min_master_distance) = 1/(1/rate+min_master_distance)

rate_desired =!= E[n_good] = 1/(1/rate+min_master_distance)
=> rate := rate_desired/(1-rate_desired*min_master_distance)
"""
function draw_correlated_spikes(trange, rates, c, shift = make_exponentialShift(1.0); min_master_distance=-0)
    # correct the rates to account for the expected loss due to pruning overlapping spikes
    frac_overlap = sum(rates)

    # draw the uncorrelated source spiketrains
    source = if c≈0.0
        @assert all(0 .<= rates*min_master_distance .< 1) "Cannot draw $(maximum(rates)) master spikes per unit-time with $(min_master_distance) minimum distance (rates must sum to < $(1/min_master_distance))."

        # correct sampling rate for expected loss of master spikes due to enforced minimum distance
        rates .= rates./(1 .- rates.*min_master_distance)
        draw_uncorrelated_spikes(trange, rates)
    else
        # determine master spike sampling rate
        v = 1/c * mean(rates.^2)/sum(rates)
        @assert 0<=length(rates)*v*min_master_distance<1 "Cannot draw ~$(length(rates)*v) master spikes per unit-time with $(min_master_distance) minimum distance (rates must sum to < $(1/min_master_distance))."

        # correct sampling rate for expected loss of master spikes due to enforced minimum distance
        v /= (1-length(rates)*v*min_master_distance)

        # determine subsampling probability
        p = (c/mean(rates.^2)) .* (rates * rates')
        @assert all(0 .<= p .<= 1) "P not all valid probabilities."

        draw_uncorrelated_spikes(trange, fill(v, length(rates)))
    end

    # ensure minimum distance between all source spike trains
    if min_master_distance > 0

        n = sum(length, source)
        times = Vector{Float64}(undef, n)
        sources = Vector{Int}(undef, n)
        i = 1
        for (s,spikes) ∈ enumerate(source)
            inew = i+length(spikes)
            times[i:inew-1] = spikes.times
            sources[i:inew-1] .= s
            i=inew
        end

        # sort times and corresponding sources
        idx=sortperm(times)
        times=times[idx]
        sources=sources[idx]

        # determine which spikes to keep ...
        keep=has_min_distance(times, min_master_distance)

        # ... and keep only those spikes for the respective source
        for (s,spikes) ∈ enumerate(source)
            source[s] = SpikeTrain(times[(sources.==s) .& keep])
        end
    end

    if c≈0.0
        return source
    else
        target = Vector{SpikeTrain}(undef, length(rates))

        # draw correlated spike trains
        for i ∈ eachindex(target)
            t = target[i] = SpikeTrain()
            for (k,s) ∈ enumerate(source)
                num_spikes = rand(Binomial(length(s.times), p[i,k]))
                append!(t.times, sample(s.times, num_spikes; replace=false))
            end

            # shift each value randomly
            t.times .= shift.(t.times)

            sort!(t.times)
        end
        return target
    end
end

"""
    draw_coincident_spikes(trange, num_spikes, p=ones(Float64,1), shift = make_exponentialShift(1.0); sorted=true)

Draws `num_spikes` Poisson spikes in the interval `trange`, each of which appears
in each spike-train `i` with the respective probability `p[i]`.
"""
function draw_coincident_spikes(trange, num_spikes, p=ones(Float64,1), shift = make_exponentialShift(1.0); sorted=true)
    times = rand(num_spikes).*(trange[2]-trange[1]).+trange[1]

    if sorted
        sort!(times)
    end

    target = Vector{SpikeTrain}(undef, length(p))
    for i ∈ eachindex(target)
        n = rand(Binomial(num_spikes, p[i]))
        target[i] = SpikeTrain(sample(times, n; replace=false))
    end

    return target
end

"""
    correlation_code(trange, stimulus, rates, correlations; kwargs...)

Each slice `stimulus[:,...,:,i]` contains an array indicating the class that the
spike-train with corresponding index represents in the time interval `trange[i:i+1]`.
Each class `c` results in spike-trains with rate `rates[c]` and mutual pairwise
correlation `correlations[c]`.
"""
function correlation_code(trange, stimulus, rates, correlations; kwargs...)
    @assert length(trange)-1 == size(stimulus)[end] "trange must include n+1 steps for stimuli with last dimension n"
    spiketrains = Array{SpikeTrain}(undef, size(stimulus))
    classes = unique(stimulus)

    for c ∈ classes
        for i ∈ Base.OneTo(length(trange)-1)
            idx = findall(s->c==s, view(stimulus, fill(Colon(), ndims(stimulus)-1)...,i))
            spiketrains[idx,i] = draw_correlated_spikes(trange[i:i+1], fill(rates[c], length(idx)), correlations[c]; kwargs...)
        end
    end

    return spiketrains
end

"""
    coincidence_code(trange, stimulus, background_rate, p; kwargs...)

Each slice `stimulus[:,...,:,i]` contains a boolean array corresponding to
whether or not a spike train with corresponding index participates or not
in the coincident spiking during the time interval `trange[i:i+1]`.
A single spike is drawn within each interval, that appears with probability `p`
in each of the spike  trains for which the corresponding entry in `stimulus` is `true`.
In addition to the signal spike for each interval, uncorrelated poisson spike-trains,
drawn with a given `background_rate`, are added as background noise.
"""
function coincidence_code(trange, stimulus, background_rate, p; kwargs...)
    @assert length(trange)-1 == size(stimulus)[end] "trange must include n+1 steps for stimuli with last dimension n"
    spiketrains = reshape(draw_uncorrelated_spikes((trange[1],trange[end]), fill(background_rate, length(stimulus))), size(stimulus))

    for i ∈ Base.OneTo(length(trange)-1)
        idx = findall(view(stimulus, fill(Colon(), ndims(stimulus)-1)...,i))
        nd = ndims(idx)+1
        spiketrains[idx,i] = merge(cat(spiketrains[idx,i], draw_coincident_spikes(trange[i:i+1], 1, fill(p, length(idx)); kwargs...), dims=nd),dims=nd)
    end

    return spiketrains
end

"""
    plot_spike_trains(spiketrains::Array{SpikeTrain}, colors=fill(:auto, length(spiketrains)), plt=plot(); kwargs...)

Draw an array of `spiketrains` with given `colors` within a single given or new plot `plt`.
"""
function plot_spike_trains(spiketrains::Array{SpikeTrain}, colors=fill(:auto, length(spiketrains)), plt=plot(); kwargs...)
    for (i,(spiketrain, color)) ∈ enumerate(zip(spiketrains, colors))
        plot!(plt, ([spiketrain.times spiketrain.times fill(NaN, length(spiketrain.times))])'[:], i .- 0.5 .+ ([zeros(length(spiketrain.times)) ones(length(spiketrain.times)) fill(NaN, length(spiketrain.times))])'[:], lc=color; kwargs...)
    end
    return plt
end

end


text:
@testset "rotation generator" begin
    all_types = (RotMatrixGenerator{3}, RotationVecGenerator,
                 RotMatrixGenerator{2}, Angle2dGenerator)
    types_2d = (RotMatrixGenerator{2}, Angle2dGenerator)
    types_3d = (RotMatrixGenerator{3}, RotationVecGenerator)

    @testset "constructor-2d" begin
        m = rand(2,2)
        s1 = RotMatrixGenerator{2}(m - m')
        s2 = RotMatrixGenerator{2, BigFloat}(m - m')
        s3 = RotMatrixGenerator{2}(0)
        s4 = RotMatrixGenerator{2, BigFloat}(0)
        s5 = RotMatrixGenerator(0)
        s6 = RotMatrixGenerator(BigFloat(0))
        s7 = Angle2dGenerator(0)
        s8 = Angle2dGenerator(BigFloat(0))
        s9 = Angle2dGenerator{BigFloat}(0)
        @test s1 isa RotMatrixGenerator{2, Float64}
        @test s2 isa RotMatrixGenerator{2, BigFloat}
        @test s3 isa RotMatrixGenerator{2, Int}
        @test s4 isa RotMatrixGenerator{2, BigFloat}
        @test s5 isa RotMatrixGenerator{2, Int}
        @test s6 isa RotMatrixGenerator{2, BigFloat}
        @test s7 isa Angle2dGenerator{Int}
        @test s8 isa Angle2dGenerator{BigFloat}
        @test s9 isa Angle2dGenerator{BigFloat}
    end

    @testset "constructor-3d" begin
        m = rand(3,3)
        s1 = RotMatrixGenerator{3}(m - m')
        s2 = RotMatrixGenerator{3, BigFloat}(m - m')
        @test s1 isa RotMatrixGenerator{3, Float64}
        @test s2 isa RotMatrixGenerator{3, BigFloat}
    end

    @testset "zero" begin
        for T in all_types
            r = zero(T)
            @test r isa RotationGenerator
            @test r === zero(r)
            @test zeros(T,2,3) == [zero(T) for i in 1:2, j in 1:3]
            @test zeros(T{BigFloat},2,3) == [zero(T{BigFloat}) for i in 1:2, j in 1:3]
            @test zeros(T,(2,3)) == [zero(T) for i in 1:2, j in 1:3]
            @test zeros(T{BigFloat},(2,3)) == [zero(T{BigFloat}) for i in 1:2, j in 1:3]
        end
    end

    @testset "one" begin
        for T in all_types
            r = one(T)
            @test r isa SMatrix
            @test r == one(r)
            @test ones(T,2,3) == [one(T) for i in 1:2, j in 1:3]
            @test ones(T{BigFloat},2,3) == [one(T{BigFloat}) for i in 1:2, j in 1:3]
            @test ones(T,(2,3)) == [one(T) for i in 1:2, j in 1:3]
            @test ones(T{BigFloat},(2,3)) == [one(T{BigFloat}) for i in 1:2, j in 1:3]
        end

        for T in types_2d
            r = one(T{BigFloat})
            @test r isa SMatrix{2,2,BigFloat}
        end

        for T in types_3d
            r = one(T{BigFloat})
            @test r isa SMatrix{3,3,BigFloat}
        end

        @test one(RotationGenerator{2}) isa SMatrix{2, 2, Float64}
        @test one(RotationGenerator{2,BigFloat}) isa SMatrix{2, 2, BigFloat}
        @test one(RotationGenerator{3}) isa SMatrix{3, 3, Float64}
        @test one(RotationGenerator{3,BigFloat}) isa SMatrix{3, 3, BigFloat}
    end

    @testset "real" begin
        @testset "$(R)" for R in all_types
            @test real(R) === R
            @test real(one(R)) === one(R)
        end
        @test real(RotationGenerator) === RotationGenerator
        @test real(RotMatrixGenerator) === RotMatrixGenerator
    end

    @testset "minus" begin
        for T in all_types
            # TODO: These should be replaced with `r = rand(T)`
            if T in types_2d
                r = T(Angle2dGenerator(1.2))
            elseif T in types_3d
                r = T(RotationVecGenerator(1.2,-0.8,0.1))
            end
            @test r isa T
            @test -r isa T
            @test r' isa T
            @test transpose(r) isa T

            @test -r == -SMatrix(r)
            @test r' == SMatrix(r)'
            @test transpose(r) == transpose(SMatrix(r))

            @test r-r == r+(-r) == zero(r)
            @test -r == r' == transpose(r)
        end
    end

    @testset "multiply" begin
        for T in all_types
            # TODO: These should be replaced with `r = rand(T)`
            if T in types_2d
                r = T(Angle2dGenerator(1.2))
            elseif T in types_3d
                r = T(RotationVecGenerator(1.2,-0.8,0.1))
            end
            a = 4.2
            @test r isa T
            @test r*1 isa T
            @test 1*r isa T
            @test a*r isa T

            @test r*1 == 1*r == r
            @test r*2 == 2*r == r+r
            @test r*a == a*r == a*SMatrix(r)
        end
    end

    @testset "division" begin
        for T in all_types
            # TODO: These should be replaced with `r = rand(T)`
            if T in types_2d
                r = T(Angle2dGenerator(1.2))
            elseif T in types_3d
                r = T(RotationVecGenerator(1.2,-0.8,0.1))
            end
            a = 4.2
            @test r isa T
            @test r/1 isa T
            @test r/a isa T

            @test r/1 == r
            @test r/2 + r/2 ≈ r
            @test r/a == SMatrix(r)/a
        end
    end

    @testset "matrix multiplication" begin
        for T in all_types
            # TODO: These should be replaced with `r = rand(T)`
            if T in types_2d
                r = T(Angle2dGenerator(1.2))
            elseif T in types_3d
                r = T(RotationVecGenerator(1.2,-0.8,0.1))
            end
            @test r isa T
            @test r*r isa SMatrix
            @test r/r isa SMatrix
            @test r^2 isa SMatrix
        end
    end

    @testset "error case" begin
        for T in types_2d
            @test_throws BoundsError zero(Angle2dGenerator)[5]
            @test_throws BoundsError zero(Angle2dGenerator)[2,3]
            @test_throws BoundsError zero(Angle2dGenerator)[3,1]
        end

        for T in types_3d
            @test_throws BoundsError zero(Angle2dGenerator)[10]
            @test_throws BoundsError zero(Angle2dGenerator)[2,4]
            @test_throws BoundsError zero(Angle2dGenerator)[4,1]
        end

        @test_throws ErrorException zero(RotationGenerator)
        @test_throws ErrorException one(RotationGenerator)
        @test_throws ErrorException zero(RotMatrixGenerator)
        @test_throws ErrorException one(RotMatrixGenerator)

        @test_throws DimensionMismatch Angle2dGenerator(1) + RotationVecGenerator(2,3,4)
    end

    @testset "params" begin
        @test Rotations.params(Angle2dGenerator(1)) == [1]
        @test Rotations.params(RotationVecGenerator(2,3,4)) == [2,3,4]
    end

    @testset "type promotion" begin
        for (T, N) in ((Angle2dGenerator, 2), (RotationVecGenerator, 3))
            R = RotMatrixGenerator

            @test zero(T)           + zero(R{N})        isa R{N, Float64}
            @test zero(R{N})        + zero(T)           isa R{N, Float64}
            @test zero(T{Int})      + zero(R{N,Int})    isa R{N, Int}
            @test zero(R{N,Int})    + zero(T{Int})      isa R{N, Int}
            @test zero(T{Int})      + zero(T{BigFloat}) isa T{BigFloat}
            @test zero(T{BigFloat}) + zero(T{Int})      isa T{BigFloat}

            @test zero(T)           - zero(R{N})        isa R{N, Float64}
            @test zero(R{N})        - zero(T)           isa R{N, Float64}
            @test zero(T{Int})      - zero(R{N,Int})    isa R{N, Int}
            @test zero(R{N,Int})    - zero(T{Int})      isa R{N, Int}
            @test zero(T{Int})      - zero(T{BigFloat}) isa T{BigFloat}
            @test zero(T{BigFloat}) - zero(T{Int})      isa T{BigFloat}

            @test 42  * zero(T)      isa T{Float64}
            @test 42  * zero(T{Int}) isa T{Int}
            @test 1.2 * zero(T{Int}) isa T{Float64}
        end
    end

    #########################################################################
    # Check that isrotationgenerator works
    #########################################################################
    @testset "Testing isrotationgenerator" begin
        # Rotation generator around x-axis
        a=[0.0 0.0 0.0
           0.0 0.0 -1.0
           0.0 1.0 0.0]
        @test isrotationgenerator(a)

        # Scaling in x-axis is not rotation generator
        a=[4.0 0.0 0.0
           0.0 0.0 0.0
           0.0 0.0 0.0]
        @test !isrotationgenerator(a)

        # Non-square matrix is not rotation generator
        @test !isrotationgenerator(zeros(2,3))
        @test !isrotationgenerator(@SMatrix zeros(2,3))

        # isrotationgenerator should work for integer
        @test isrotationgenerator([0 1 0; -1 0 0; 0 0 0])

        # zero matrix is rotation generator
        @test isrotationgenerator(zeros(1,1))
        @test isrotationgenerator(zeros(3,3))
        @test isrotationgenerator(zeros(4,4))
        @test isrotationgenerator(zero(SMatrix{1,1}))
        @test isrotationgenerator(zero(SMatrix{3,3}))
        @test isrotationgenerator(zero(SMatrix{4,4}))
        @test isrotationgenerator(zero(RotMatrix{1}))
        @test isrotationgenerator(zero(RotMatrix{3}))
        @test isrotationgenerator(zero(RotMatrix{4}))

        # Rotation matrix can be AbstractMatrix{<:Complex}
        @test isrotationgenerator(zeros(3,3) .+ 0im)
        @test isrotationgenerator(zero(SMatrix{4,4}) .+ 0im)

        # Including NaNs are not rotaion
        @test !isrotationgenerator(RotationVecGenerator(0, 0, NaN))
        @test !isrotationgenerator(Angle2d(NaN))

        # Complex unitary matrix is not rotation
        M = randn(Complex{Float64},3,3)
        @test !isrotationgenerator(M+M')
        @test !isrotationgenerator(M+transpose(M))
        @test !isrotationgenerator(real(M+M'))
        @test !isrotationgenerator(real(M+transpose(M)))
    end

    @testset "Testing show" begin
        io = IOBuffer()
        r = zero(RotMatrixGenerator{2})
        show(io, MIME("text/plain"), r)
        str = String(take!(io))
        @test startswith(str, "2×2 RotMatrixGenerator2{Float64}")

        rvec = RotationVecGenerator(1.0, 2.0, 3.0)
        show(io, MIME("text/plain"), rvec)
        str = String(take!(io))
        @test startswith(str, "3×3 RotationVecGenerator{Float64}") && occursin("(1.0, 2.0, 3.0):", str)
    end
end


text:
(nworkers() < ncpu) && (addprocs(ncpu - nprocs() + 1))
require("CUDArt")
@everywhere CUDArt.device((myid()-1) % CUDArt.devcount())
require("CUBLAS")
require("KUnet")
require("KUparser")


text:
#=
This file calculates the raw video sizes for several screen resolutions.
The output is included in PupperSim.jl as a comment block for the loadmodel function.

> julia calculate_raw_videos_sizes.jl
=#

function calculate_raw_videos_sizes(max_frames)
    whparams(w, h) = "$(lpad(w, 4)), $(lpad(h, 4))" # helper function
    wharray = [1920 1080; 1600 900; 1200 900; 1024 768; 800 600; 512 384; 400 300]   # array of resolution options

    for i = 1:size(wharray)[1]          # 1:number of rows
        wh = wharray[i, :]
        aspect_ratio = wh[1,1]//wh[2,1] # rational representation of w/h (w//h)
        @assert aspect_ratio == 16//9 || aspect_ratio == 4//3  "unrecognized aspect ratio: $aspect_ratio"
        pixels = *(wh...)               # number of pixels per image
        bytes = pixels * 3              # 3 bytes per pixel
        mb_per_frame = round(bytes/1048576, digits=2)
        raw_video_size = round((max_frames * mb_per_frame)/1024, digits=1)
        println("# At resolution ($(whparams(wh...))): $mb_per_frame MB/frame, total raw video size: $(lpad(raw_video_size, 4)) GB")
    end
end

# These constants copied from PupperSim.jl
const max_video_duration = 60       # max video duration in seconds
const video_frames_per_second = 30  # determined by GLFW.GetPrimaryMonitor refresh rate
const max_video_frames = video_frames_per_second * max_video_duration

calculate_raw_videos_sizes(max_video_frames)


text:
function arclength(φ₀)
# A terms
    Acoef =  [  -3/2*[ 1.0  0.0  -1/8      0.0      -1/64]; # A₁
               15/16*[ 0.0  1.0   0.0     -1/4        0.0]; # A₂
              -35/48*[ 0.0  0.0   1.0      0.0      -5/16]; # A₃
                     [ 0.0  0.0   0.0  315/512        0.0]; # A₄
                     [ 0.0  0.0   0.0      0.0  -693/1280]; # A₅
             ]
    Aj = Acoef * [n^k for k=1:5]

    ΣA = sum([Aj[j]*sin(2j*φ₀) for j=1:5])
    S̅φ₀ = (m₀*a)/(1+n) * (A₀*φ₀ + ΣA)

    return S̅φ₀
end


text:
import WebIO: camel2kebab

# Get median elements of ranges, used for initialising sliders.
# Differs from median(r) in that it always returns an element of the range
medianidx(r) = (1+length(r)) ÷ 2
medianelement(r::AbstractArray) = r[medianidx(r)]
medianval(r::Associative) = medianelement(collect(values(r)))
medianelement(r::Associative) = medianval(r)

_values(r::AbstractArray) = r
_values(r::Associative) = values(r)

_keys(r::AbstractArray) = 1:length(r)
_keys(r::Associative) = keys(r)

inverse_dict(d::Associative) = Dict(zip(values(d), keys(d)))

const Propkey = Union{Symbol, String}
const PropDict = Dict{Propkey, Any}

function slap_design!(w::Scope, args)
    for arg in args
        import!(w, arg)
    end
    w
end

slap_design!(w::Scope, args::AbstractString...) = slap_design!(w::Scope, args)

slap_design!(w::Scope, args::WidgetTheme = gettheme()) =
    slap_design!(w::Scope, libraries(args))

slap_design!(n::Node, args...) = slap_design!(Scope()(n), args...)

slap_design!(w::Widget, args...) = (slap_design!(scope(w), args...); w)

isijulia() = isdefined(Main, :IJulia) && Main.IJulia.inited


text:
using PyPlot

struct Vertex
    neighbors::Vector{Int}       # Indices of neighbors of this Vertex
    coordinates::Vector{Float64} # 2D coordinates of this Vertex - only for plotting
    Vertex(neighbors; coordinates=[0,0]) = new(neighbors, coordinates)
end

function Base.show(io::IO, v::Vertex)
    print(io, "Neighbors = ", v.neighbors)
end

struct Graph
    vertices::Vector{Vertex}
end

function Base.show(io::IO, g::Graph)
    for i = 1:length(g.vertices)
        println(io, "Vertex $i, ", g.vertices[i])
    end
end

function PyPlot.plot(g::Graph; scale=1.0)
    fig, ax = subplots()
    ax.set_aspect("equal")

    xmin = minimum(v.coordinates[1] for v in g.vertices)
    xmax = maximum(v.coordinates[1] for v in g.vertices)
    ymin = minimum(v.coordinates[2] for v in g.vertices)
    ymax = maximum(v.coordinates[2] for v in g.vertices)
    sz = max(xmax-xmin, ymax-ymin)
    cr = scale*0.05sz
    hw = cr/2
    axis([xmin-2cr,xmax+2cr,ymin-2cr,ymax+2cr])
    axis("off")

    for i in 1:length(g.vertices)
        c = g.vertices[i].coordinates
        ax.add_artist(matplotlib.patches.Circle(c, cr, facecolor="none", edgecolor="k"))
        ax.text(c[1], c[2], string(i),
                horizontalalignment="center", verticalalignment="center", fontsize=round(Int, 14*scale))
        for nb in g.vertices[i].neighbors
            cnb = g.vertices[nb].coordinates
            dc = cnb .- c
            L = sqrt(sum(dc.^2))
            c1 = c .+ cr/L * dc
            c2 = cnb .- cr/L * dc
            arrow(c1[1], c1[2], c2[1]-c1[1], c2[2]-c1[2],
                  head_width=hw, length_includes_head=true, facecolor="k")
        end
    end
end

function circle_graph(nv=8)
    g = Graph([])
    for i = 1:nv
        th = 2π*i/nv
        v = Vertex([mod(i,nv)+1, mod(i-2,nv)+1], coordinates=[cos(th), sin(th)])
        push!(g.vertices, v)
    end
    return g
end

g = circle_graph(10)

# Add edges by growing the neighbors array
push!(g.vertices[1].neighbors, 6);       # Add edge 1->6
append!(g.vertices[4].neighbors, [2,7]); # Add 2 new edges

# Change edges by modifying the neighbors array
g.vertices[3].neighbors[1] = 1           # Change edge, 3->1 instead of 3->4
pop!(g.vertices[5].neighbors)            # Remove edge

# Add a new vertex by growing the vertices array, and if needed
# updating some of the existing neighbors array
# Note that the new vertex should be numbered last, otherwise
# all node numbers must be reordered
newv = Vertex([1,4,9], coordinates=[-.1,.3])
push!(g.vertices, newv)
push!(g.vertices[2].neighbors, 11)

# Deleting a vertex is difficult using adjancency lists
# (unless it is the last vertex), since it requires renumbering
# of the vertices. But as a compromise, we can remove all edges
# to/from a vertex:
resize!(g.vertices[8].neighbors,0)
for v in g.vertices
    filter!(i -> i != 8, v.neighbors)
end


text:
const Deckinfo = Array{NamedTuple{(:name, :amount, :prints),Tuple{String,Int64,Array{NamedTuple{(:index, :set, :rarity),Tuple{Int64,Symbol,Int64}},1}}},1}

"""
    SimParameters

Parameters to adjust certain aspects of the simulation:
- `pitytimer_distribution`: Collection of samples determining the distribution of parameters
   used in the pack opening and pitytimers.  See `?PityTimerSamples`
- `wc_timer_track`: Vector representing the WC rarities associated with the timer tracks.  Each pack
   opened advances the position and then resets at the end.  0 is associated with no wildcard, 2 to 4
   correspond to uncommon through mythic.
- `icrs_per_pack`: A `Cycle` of vectors of 4-dimensional categorical distributions representing
   probabilistic ICR rarities associated with each pack, e.g. to simulate doing constructed events
   in between buying packs.  Each pack opened advances to the next vector of distributions, and then an
   ICR rarity is drawn from each distribution in the vector.
- `pack_vault_pct`: Vault percentage associated with opening a pack.  Defaults to 0% as of M19.
- `duplicate_vpct`: Vault percentage associated with an excess copy of a card of each rarity
- `noopen_cards`: Cards that don't show up in packs on account of being promos or banned
- `starter_cards`: Cards from NPE decks assumed to be in everyone's collection relatively quickly
- `bonus_packs`: One-time packs from, e.g. NPE, that don't count against the total packs
"""
@with_kw mutable struct SimParameters
    pitytimer_distribution::PityTimerSamples = pity_timer_samps
    wc_timer_track::Vector{Int} = [0,0,2,0,0,3,0,0,2,0,0,3,0,0,2,0,0,4,0,0,2,0,0,3,0,0,2,0,0,3]
    icrs_per_pack::Cycle{Vector{Categorical{Float64}}} = icrgen_qc(0.5,0,0)
    pack_vault_pct::Float64 = 0.0
    duplicate_vpct::NTuple{4,Float64} = (1/9, 1/3, 5/9, 1+1/9)
    noopen_cards::Vector{Int} = map(x -> first(x.prints).index, deckinfo(deckreader_mtga_format(joinpath(@__DIR__, "..", "data", "noopen_cards.txt"))))
    starter_cards::Vector{Int} = tenstarters
    fixed_starter::Bool = true
    bonus_packs::Dict{Symbol,Int} = Dict( :GRN => 4, :M19 => 5 )
    nextset::Function = next_set!
    wc_upgrade_rate::Vector{Float64} = [Inf,Inf,Inf,Inf]
    wc_upgrade_threshold::Vector{Int} = [20,20,100,0]
    welcome_bundle::Bool = false
    kaladesh_grant::Bool = false
    prevent_duplicates::Bool = false
    starting_wc_count::Vector{Int} = [8,4,2,1]
    welcome_bundle_set::Symbol = :M19
    fixed_pack_track::Cycle{Union{Missing, Symbol}} = Cycle([missing,:GRN,missing,:GRN,missing, :GRN, missing, missing, missing, missing])
    openable_schedule::Dict{Int,Vector{Symbol}} = Dict{Int,Vector{Symbol}}()
    track_collection_progress::Bool = false
    max_track_progress_packs::Int = 500
end

"""
    AccountState

Represents non-card aspects of an account: Vault percentage, wildcard counts,
bonus packs that don't count against output statistics (e.g. from NPE), pitytimer
for pack-based wildcards, and the position on the wildcard timer tracks
"""
@with_kw mutable struct AccountState
    vault_pct::Float64 = 0.0
    wc_count::Vector{Int} = [0,0,0,0]
    bonus_packs::Dict{Symbol,Int} = Dict( (s, 0) for s in sets)
    pitytimer::Vector{Int} = [0,0,0,0]
    wc_timer::Cycle{Int} = Cycle(SimParameters().wc_timer_track, 0)
    fixed_pack_track::Cycle{Union{Missing, Symbol}} = Cycle(SimParameters().fixed_pack_track.data, 0)
end

function AccountState(parameters::SimParameters)
    account = AccountState()
    account.wc_count = deepcopy(parameters.starting_wc_count)

    for set in keys(parameters.bonus_packs)
        account.bonus_packs[set] += parameters.bonus_packs[set]
    end

    account.wc_timer = Cycle(parameters.wc_timer_track, 0)
    account.fixed_pack_track = Cycle(parameters.fixed_pack_track.data,0)

    account
end

"""
    SimOutput

Summarizes simulation output in terms of:
- `packs_opened`: (Non-bonus) packs opened of each set
- `wcs_from_vault`: How many wildcards came from Vault
- `wcs_spent_on_decks`: How many wildcards were spent on decks of each rarity
- `card_sources`: An estimate of card sources for each deck. Array dimensions:
    1. Simulation rep
    2. Deck order
    3. Card rarity (1 = common, ..., 4 = mythic)
    4. Various statistics, corresponding to this 9 'columns' in this dimension:
        1. Total cards of this rarity in deck (non-stochastic)
        2. Amount needed of this rarity before starting this deck
        3. Amount needed of this rarity after finishing this deck (should be 0)
        4. Cards of this rarity in this deck crafted by wildcards during its construction
        5. Cards of this rarity in this deck opened from packs during its construction
        6. Cards of this rarity in this deck opened from ICRs during its construction
        7. As in 4, except during all decks so far
        8. As in 5, except during all decks so far
        9. As in 6, except during all decks so far
- `total_pack_wcs`: Running total across an entire iteration of wildcards for each rarity
- `total_packs`: Running total across an entire iteration of packs, free or not, for each rarity
"""
@with_kw mutable struct SimOutput
    packs_opened::Dict{Symbol,Array{Int,2}}
    bonus_packs_opened::Dict{Symbol,Array{Int,2}}
    wcs_from_vault::Array{Int,2}
    wcs_spent_on_decks::Array{Int,3}
    card_sources::Array{Float64, 4}
    total_pack_wcs::Array{Int,2}
    total_packs::Vector{Int}
    total_earned_packs::Vector{Int}
    ending_collection::Array{Int,2}
    collection_progress::Array{Union{Missing,Int},3}
end

function SimOutput(reps::Int, decks::Int, sets::Vector{Symbol}, collection::Vector{Int})::SimOutput
    SimOutput(
        packs_opened = Dict( (s, zeros(Int, (reps, decks))) for s in sets),
        bonus_packs_opened = Dict( (s, zeros(Int, (reps, decks))) for s in sets),
        wcs_from_vault = zeros(Int, (reps, 4)),
        wcs_spent_on_decks = zeros(Int, (reps, decks, 4)),
        card_sources = zeros(Float64, (reps, decks, 4, 9)),
        total_pack_wcs = zeros(Int, (reps, 4)),
        total_packs = zeros(Int, reps),
        total_earned_packs = zeros(Int, reps),
        ending_collection = zeros(Int, (reps, length(collection))),
        collection_progress = Array{Union{Missing,Int}}(missing, 0,0,0)
    )
end

function copy_into_output!(total::SimOutput, partial::SimOutput; startindex::Int)
    endindex = startindex + size(partial.card_sources,1) - 1
    range = startindex:endindex

    total.wcs_from_vault[range,:] = partial.wcs_from_vault
    total.wcs_spent_on_decks[range,:,:] = partial.wcs_spent_on_decks
    total.card_sources[range,:,:,:] = partial.card_sources
    total.total_pack_wcs[range,:] = partial.total_pack_wcs
    total.total_packs[range] = partial.total_packs
    total.total_earned_packs[range] = partial.total_earned_packs
    total.ending_collection[range,:] = partial.ending_collection

    if length(total.collection_progress) > 0 && length(partial.collection_progress) > 0
        total.collection_progress[range,:,:] = partial.collection_progress
    end

    for set in keys(total.packs_opened)
        total.packs_opened[set][range,:] = partial.packs_opened[set]
    end

    for set in keys(total.bonus_packs_opened)
        total.bonus_packs_opened[set][range,:] = partial.bonus_packs_opened[set]
    end

    endindex + 1
end

max_stars_repo_name: asieev/ArenaSim

max_issues_repo_name: asieev/ArenaSim

max_forks_repo_name: asieev/ArenaSim


text:
# add more julia processes
nprocs() < 4 ? addprocs(4-nprocs()) : nothing

using IncrementalInference
using HDF5, JLD, Gadfly, Colors, Cairo
using KernelDensityEstimate, Distributions
using Caesar, RoME
using RoMEPlotting
using DocStringExtensions

## load all the model data
# d = odometry information
# f = laser scanner detections
# MM = multi-modal individual id references
# MMr = reworked to map to only one previous feature
# examplefolder, datafolder

include(joinpath(Pkg.dir("Caesar"),"examples","wheeled","loadVicPrkData.jl"))

# include(joinpath(ENV["HOME"],"Documents","blandauthlocal.jl"))
# backend_config, user_config = standardcloudgraphsetup(addrdict=addrdict)
# # Start new session
# addrdict["sessionId"] = "VICPRK_VID"
# addrdict["robotId"] = "Ute"
# fg = initfg(sessionname=user_config["sessionId"], robotname=user_config["robotId"], cloudgraph=backend_config)
# deleteServerSession!(fg.cg, user_config["sessionId"])

include(joinpath(Pkg.dir("Caesar"),"examples","wheeled","VicPrk_FrontEndUtils.jl"))

## Example Construction and pieces

Podo=diagm([0.5;0.5;0.005])
lsrNoise=diagm([0.1;1.0])

# start factor graph
N=100
fg = initfg()
lmoccur = Dict{Symbol, Int}()

# init pose
prevn = initFactorGraph!(fg, init=d[1][1:3])[1]

measurements = f[1]
findAddNewLandmarks!(fg, lmoccur, measurements)
# there might not be any landmarks yet
findAddLandmFactorsByPose!(fg, :x0, f[1], N=N)

# pbr = Pose2Point2BearingRange(Normal(br[1], cov[1,1]), Normal(br[2],  cov[2,2]))  #{Normal, Normal}
# f = addFactor!(fg, [vps;vlm], pbr, ready=ready, autoinit=true ) #[vps;vlm],

drive!(fg, 2, d, f, lmoccur, Podo, N=N)

ls(fg)

writeGraphPdf(fg)

tree = wipeBuildNewTree!(fg, drawpdf=true)
inferOverTree!(fg, tree,N=N)

@async run(`evince bt.pdf`)

PL = Dict{Int, Any}()

PL[2] = drawPosesLandms(fg)

PL[3] = donextframe!(fg, 3, d, f, lmoccur, Podo, N=N)

PL[4] = donextframe!(fg, 4, d, f, lmoccur, Podo, N=N)

PL[5] = donextframe!(fg, 5, d, f, lmoccur, Podo, N=N)

PL[6] = donextframe!(fg, 6, d, f, lmoccur, Podo, N=N)

PL[7] = donextframe!(fg, 7, d, f, lmoccur, Podo, N=N)

PL[8] = donextframe!(fg, 8, d, f, lmoccur, Podo, N=N)

PL[9] = donextframe!(fg, 9, d, f, lmoccur, Podo, N=N)

PL[10] = donextframe!(fg, 10, d, f, lmoccur, Podo, N=N)

PL[11] = donextframe!(fg, 11, d, f, lmoccur, Podo, N=N)

PL[12] = donextframe!(fg, 12, d, f, lmoccur, Podo, N=N)

PL[13] = donextframe!(fg, 13, d, f, lmoccur, Podo, N=N)

PL[14] = donextframe!(fg, 14, d, f, lmoccur, Podo, N=N)

PL[15] = donextframe!(fg, 15, d, f, lmoccur, Podo, N=N)

PL[16] = donextframe!(fg, 16, d, f, lmoccur, Podo, N=N)

PL[17] = donextframe!(fg, 17, d, f, lmoccur, Podo, N=N)

PL[18] = donextframe!(fg, 18, d, f, lmoccur, Podo, N=N)

PL[19] = donextframe!(fg, 19, d, f, lmoccur, Podo, N=N)

PL[20] = donextframe!(fg, 20, d, f, lmoccur, Podo, N=N)

PL[21] = donextframe!(fg, 21, d, f, lmoccur, Podo, N=N)

PL[22] = donextframe!(fg, 22, d, f, lmoccur, Podo, N=N)

PL[23] = donextframe!(fg, 23, d, f, lmoccur, Podo, N=N)

# importall RoMEPlotting
# drawSubmaps(fg, [1; ], spread=5, m2hist=true)

# get almost smallest neighboring landmarks to last pose
nhbls = ls2(fg, ls(fg)[1][end])
nhbls = nhbls[symbolmatch.(nhbls, 'l')]
lmn = landmarknumber.(nhbls)
minwindow = round(Int, quantile(lmn,0.01, sorted=false))

# largest landmark number
lmmax = landmarknumber(ls(fg)[2][end])

drawSubmaps(fg, [1  minwindow-1; minwindow lmmax], m2hist=true)

# on robot front-end acquistion

# Figure export folder
currdirtime = now()
imgdir = joinpath(ENV["HOME"], "Pictures", "vicprkimgs", "$(currdirtime)")
mkdir(imgdir)

lcmode=:mmodal

#build :# poses for the factorGraph
for idx=9:30
  prev, X, nextn = getLastPose2D(fg)
  vp, fp = addOdoFG!(fg, nextn, d[idx][1:3], Podo, N=N)
  # add landmarks
  addLandmarksFactorGraph!(fg, f, idx, prevn, nextn, lcmode=lcmode, lsrNoise=lsrNoise, N=N, MM=MM)
  prevn = nextn
  if (idx%10==0)
     # Solve
     tree = wipeBuildNewTree!(fg, drawpdf=false);
    @time inferOverTree!(fg,tree, N=100);
  end
  pl=drawPosesLandms(fg,window=(nextn,25))
  Gadfly.draw(PNG(joinpath(imgdir,"$(nextn).png"),20cm,20cm),pl)
end

0

# run(`ffmpeg -y -i x%d.png -threads 4 -vcodec libx264 -s 1920x1080 -b:v 2M -filter:v "setpts=10*PTS" /home/dehann/Videos/vicprk_test.mp4`)

# pl=drawPosesLandms(fg);
# Gadfly.draw(PDF("/tmp/before.pdf",20cm,20cm),pl)

# batchSolve(fg)
# on ssh terminal, run slamindb(iterations=1)

# fetch a local copy

fg = initfg(sessionname=user_config["sessionId"], cloudgraph=backend_config) #, robotname=user_config["robotId"]
fullLocalGraphCopy!(fg)
# pl1=drawPosesLandms(fg)

pl1=drawPosesLandms(fg,window=(:x9,25))

Gadfly.draw(PDF("/tmp/after.pdf",20cm,20cm),pl1)

@async run(`evince /tmp/after.pdf`)

subLocalGraphCopy!(fg, string.([:x1,:x2,:x3,:x4,:x5]), neighbors=1)

X,L = ls(fg)
d=0
t=[]
z=zeros(length(L)+1,length(L)+1)

j=1;
for l1 in L
  i=1
  for l2 in L
    b=getKDEMax(getVertKDE(fg, l1))
    c=evalLikelihood(fg, l2 , b)
    if ( c>0.002 && l1!=l2 )
      d+=1
      #Replace this with A delete and replace
      #what should i do now
      addSoftEqualityPoint2D(fg, l1, l2 )

      push!(t,(l1,l2))
      z[i,j]=1
    elseif (l1==l2)
      z[i,j]=2
    end
    i+=1
  end
  j+=1
end
#run Server Code
fg = initfg(sessionname=user_config["sessionId"], cloudgraph=backend_config)
fullLocalGraphCopy!(fg)
pl2=drawPosesLandms(fg)
draw(PDF("daniel.pdf",20cm,20cm),pl2)

# Remove the new session from the server
deleteServerSession!(fg.cg, user_config["sessionId"])

#


text:
@testset "88.merge-sorted-array.jl" begin
    nums1 = [1, 2, 3, 0, 0, 0]
    m = 3
    nums2 = [2, 5, 6]
    n = 3

    merge_sorted_array(nums1, m, nums2, n)
    @test nums1 == [1, 2, 2, 3, 5, 6]
end


text:
#quality of life constants
const SIZE = NTuple{N,Int} where N
const LENGTH = SIZE{1}

#prototype abstract, stores the storage type in S
abstract type AbstractPrototype{S} end

#prototype stores the type to dispatch and the data necesary to create a new instance of other element types.
struct Prototype{S,T}
    type::S
    data::T #a tuple to store associated data
end

prototype(x::T) where T<:Array = Prototype(T,(eltype(x),size(x)))
prototype(d::T) where T<:Dict = Prototype(T,(valtype(d),d.slots, d.keys, d.ndel, d.count, d.age, d.idxfloor, d.maxprobe))
prototype(sv::T) where T<:SparseArrays.SparseVector = Prototype(T,(eltype(sv),sv.n,sv.nzind))
prototype(sm::T) where T<:SparseArrays.SparseMatrixCSC = Prototype(T,(eltype(sm),size(sm),sm.rowval,sm.colptr))

function make_cache!(f!,p::P,elem_type) where {P <: Prototype}
    return make_cache!(p.type,f!,p,elem_type)
end

function make_cache!(::Type{T1},f!,p::P,elem_type::Type{T2}) where {T1 <: Array, P <: Prototype, T2}
    s = p.data[2]
    n = length(s)
    return Array{T2,n}(undef,s)
end

#the approach here is more like Dictionaries.jl
function make_cache!(::Type{T1},f!,p::P,elem_type::Type{T2}) where {T1 <: Dict, P <: Prototype, T2}
    t, slots, keys, ndel, count, age, idxfloor, maxprobe = p.data
    n = length(slots)
    K = eltype(keys)
    vals = Vector{T2}(undef,n)
    return Dict{K,T2}(slots,keys,vals,ndel,count,age,idxfloor,maxprobe)
end

function make_cache!(::Type{T1},f!,p::P,elem_type::Type{T2}) where {T1 <: SparseArrays.SparseVector, P <: Prototype, T2,K,V}
    t,n,nzind = p.data
    nzval = Vector{T2}(undef,n)
    return SparseArrays.SparseVector{T2,eltype(nzind)}(n,nzind,nzval)
end

function make_cache!(::Type{T1},f!,p::P,elem_type::Type{T2}) where {T1 <: SparseArrays.SparseMatrixCSC, P <: Prototype, T2,K,V}
    t,s,rowval,colptr = p.data
    m,n = s
    l = length(rowval)
    t,n,nzind = p.data
    nzval = Vector{T2}(undef,l)
    return SparseArrays.SparseVector{T2,eltype(nzind)}(m,n,colptr,nzind,nzval)
end


text:
function part1(filename::AbstractString)
    sum((scoreline(line) for line in eachline(filename)), init = 0)
end

scoreline(line::AbstractString) = scoreline(lastunmatched(line))
scoreline(_) = 0
scoreline(c::Char) = if c in ')'
    3
elseif c == ']'
    57
elseif c == '}'
    1197
elseif c == '>'
    25137
else
    0
end

partnerof(c::Char) = ")(][}{><"[findfirst(c, "()[]{}<>")]

"""
Returns one of:
- `c::Char` an unmatched right bracket
- `cs::Vector{Char}` the unmatched left brackets
- `nothing` brackets are matched
"""
function lastunmatched(line::AbstractString)
    stack = Char[]
    for c in line
        if c in "{[(<"
            push!(stack, c)
        elseif c in "}])>"
            if pop!(stack) != partnerof(c)
                return c
            end
        end
    end
    if !isempty(stack)
        stack
    end
end

if basename(pwd()) == "aoc"
    cd("2021/10")
end
@assert part1("example.txt") == 26397
part1("input.txt")

function part2(filename::AbstractString)
    scores = autocompletescores.(
        filter!(
            Base.Fix2(isa, Vector{Char}),
            lastunmatched.(readlines(filename))
        )
    )
    sort!(scores)[1 + length(scores) ÷ 2]
end

function autocompletescores(chars)
    mapfoldr(c -> findfirst(c, "([{<"), (a, b) -> a + 5b, chars, init=0)
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
for b = 1:nbins
    img = b*ones(10,10)
    integral = integralimage(img, nbins, minval = 0, maxval = 255)
    nrow, ncol = size(img)
    for r = 1:nrow
        for c = 1:ncol
            @test integral[r,c,b] == r*c
        end
    end
end


text:
#! /usr/bin/julia

# Rosetta Code, Find limit of recursion

function divedivedive()
    global depth
    depth += 1
    divedivedive()
end

function divedivedive(d::Int)
    try
        divedivedive(d+1)
    catch
        return d
    end
end

depth = divedivedive(0)
println("A clean dive reaches a depth of ", depth, ".")

depth = 0
try
    divedivedive()
end
println("A dirty dive reaches a depth of ", depth, ".")

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code