text:
@inline cld_fast(x, y) = cld(x, y)
@inline function cld_fast(x::I, y) where {I <: Integer}
    # ux = unsigned(x); uy = unsigned(y)
    # d = Base.udiv_int(ux, uy)
    d = div_fast(x, y)
    (d + (d * unsigned(y) != unsigned(x))) % I
    # ifelse(d * uy == ux, d, d + one(d)) % I
end
cld_fast(::StaticInt{N}, ::StaticInt{M}) where {N,M}= (StaticInt{N}() + StaticInt{M}() + One()) ÷ StaticInt{M}()
@inline function divrem_fast(x::I, y) where {I <: Integer}
    ux = unsigned(x); uy = unsigned(y)
    d = Base.udiv_int(ux, uy)
    r = ux - d * uy
    # r = vsub(ux, vmul(d, uy))
    d % I, r % I
end
@inline divrem_fast(::StaticInt{x}, y::I) where {x, I <: Integer} = divrem_fast(x % I, y)
@inline div_fast(x::I, y::Integer) where {I <: Integer} = Base.udiv_int(unsigned(x), unsigned(y)) % I
@inline div_fast(::StaticInt{x}, y::I) where {x, I <: Integer} = Base.udiv_int(unsigned(x), unsigned(y)) % I
# @inline div_fast(x::I, ::StaticInt{x}) where {x, I <: Integer} = Base.udiv_int(x % UInt32, y % UInt32) % I
divrem_fast(::StaticInt{N}, ::StaticInt{M}) where {N,M}= divrem(StaticInt{N}(), StaticInt{M}())
div_fast(::StaticInt{N}, ::StaticInt{M}) where {N,M}= StaticInt{N}() ÷ StaticInt{M}()
@generated function div_fast(x::I, ::StaticInt{M}) where {I<:Integer,M}
    if VectorizationBase.ispow2(M)
        lm = VectorizationBase.intlog2(M)
        Expr(:block, Expr(:meta,:inline), :(x >>> $lm))
    else
        Expr(:block, Expr(:meta,:inline), :(div_fast(x, $(I(M)))))
    end
end

function Base.copyto!(B::AbstractStrideArray{<:Any,<:Any,<:Any,N}, A::AbstractStrideArray{<:Any,<:Any,<:Any,N}) where {N}
    @avx for I ∈ eachindex(A, B)
        B[I] = A[I]
    end
    B
end

@generated _max(::StaticInt{N}, ::StaticInt{M}) where {N,M} = :(StaticInt{$(max(N,M))}())
const MᵣW_mul_factor = VectorizationBase.REGISTER_SIZE === 64 ? StaticInt{4}() : StaticInt{9}()

function effective_l2(::Type{T}) where {T}
    if CACHE_INCLUSIVITY[2]
        something(core_cache_size(T, Val(2)), StaticInt{262144}() ÷ static_sizeof(T)) - (something(core_cache_size(T, Val(1)), StaticInt{32768}() ÷ static_sizeof(T)) * StaticInt{CACHE_COUNT[1]}()) ÷ StaticInt{CACHE_COUNT[2]}()
    else
        something(core_cache_size(T, Val(2)), StaticInt{262144}() ÷ static_sizeof(T))
    end
end
function effective_l3(::Type{T}) where {T}
    if CACHE_INCLUSIVITY[3]
        something(cache_size(T, Val(3)), StaticInt{3145728}() ÷ static_sizeof(T)) - (effective_l2(T) * StaticInt{CACHE_COUNT[2]}()) ÷ StaticInt{CACHE_COUNT[3]}()
    else
        something(cache_size(T, Val(3)), StaticInt{3145728}() ÷ static_sizeof(T))
    end
end
effective_cache(::Type{T}, ::Val{1}) where {T} = something(core_cache_size(T, Val(1)), StaticInt{32768}() ÷ static_sizeof(T))
effective_cache(::Type{T}, ::Val{2}) where {T} = effective_l2(T)
effective_cache(::Type{T}, ::Val{3}) where {T} = effective_l3(T)
function matmul_params(::Type{T}) where {T}
    Mᵣ = StaticInt{mᵣ}()
    Nᵣ = StaticInt{nᵣ}()
    W = VectorizationBase.pick_vector_width_val(T)
    MᵣW = Mᵣ * W
    mc = MᵣW * MᵣW_mul_factor # TODO: make this smarter/less heuristic
    # L₁ = something(core_cache_size(T, Val(1)), StaticInt{32768}() ÷ static_sizeof(T))
    L₂ = effective_l2(T)
    # _kc = ((StaticInt{795}() * L₂) ÷ StaticInt{1024}() - StaticInt{30420}()) ÷ mc
    _kc = ((StaticInt{795}() * L₂) ÷ StaticInt{1024}() - StaticInt{15250}()) ÷ mc
    # _kc = ((StaticInt{3}() * L₂) ÷ (StaticInt{4}() * mc)) - StaticInt{130}()
    kc = _max(_kc, StaticInt{120}())
    L₃ = effective_l3(T)
    _nc = ((((StaticInt{132}() * L₃) ÷ StaticInt{125}()) - StaticInt{256651}()) ÷ (kc * Nᵣ)) * Nᵣ
    # _nc = ((StaticInt{3}() * L₃) ÷ (StaticInt{4}() * kc * Nᵣ)) * Nᵣ
    nc = _max(_nc, StaticInt{400}())
    mc, kc, nc
end

@generated function dense_dims_subset(::DenseDims{D}, ::StrideRank{R}) where {D,R}
    t = Expr(:tuple)
    for n in eachindex(R)
        push!(t.args, D[n] & (R[n] == 1))
    end
    Expr(:call, Expr(:curly, :DenseDims, t))
end

"""
Only packs `A`. Primitively does column-major packing: it packs blocks of `A` into a column-major temporary.
"""
function jmulpackAonly!(
    C::AbstractStridedPointer{T}, A::AbstractStridedPointer, B::AbstractStridedPointer,
    α, β, ::StaticInt{mc}, ::StaticInt{kc}, ::StaticInt{nc}, (M,K,N)
) where {T, mc, kc, nc}
    W = VectorizationBase.pick_vector_width_val(T)
    mᵣW = StaticInt{mᵣ}() * W

    Miter = cld_fast(M, StaticInt{mc}())
    Mbsize, Mrem, Mremfinal, Mblocks = split_m(M, Miter, W)
    Mblock_Mrem, Mblock_ = promote(Mbsize + W, Mbsize)

    Kiter = cld_fast(K, StaticInt{kc}())
    Kblock, Krem = divrem_fast(K, Kiter)
    Kblock_Krem, Kblock_ = promote(Kblock + One(), Kblock)

    Aptr = zero_offsets(A); Bptr = zero_offsets(B); Cptr = zero_offsets(C);
    # @show Miter, Kiter, M, K, N
    for ko ∈ CloseOpen(Kiter)
        ksize = ifelse(ko < Krem, Kblock_Krem, Kblock_)
        let Aptr = Aptr, Cptr = Cptr
            for mo in CloseOpen(Miter)
                msize = ifelse(mo == (Miter - 1), Mremfinal, ifelse(mo < Mrem, Mblock_Mrem, Mblock_))
                if ko == 0
                    packaloopmul!(Cptr, Aptr, Bptr, α, β, (msize, ksize, N))
                else
                    packaloopmul!(Cptr, Aptr, Bptr, α, One(), (msize, ksize, N))
                end
                Aptr = gesp(Aptr, (msize, Zero()))
                Cptr = gesp(Cptr, (msize, Zero()))
            end
        end
        Aptr = gesp(Aptr, (Zero(), ksize))
        Bptr = gesp(Bptr, (ksize, Zero()))
    end
    # end # GC.@preserve
    nothing
end
"""
Packs both arrays `A` and `B`.
Primitely packs both `A` and `B` into column major temporaries.

Column-major `B` is preferred over row-major, because without packing the stride across `k` iterations of `B` becomes excessive, and without `nᵣ` being a multiple of the cacheline size, we would fail to make use of 100% of the loaded cachelines.
Unfortunately, using column-major `B` does mean that we are starved on integer registers within the macrokernel.

Once `LoopVectorization` adds a few features to make it easy to abstract away tile-major memory layouts, we will switch to those, probably improving performance for larger matrices.
"""
function jmulpackAB!(
    C::AbstractStridedPointer{T}, A::AbstractStridedPointer, B::AbstractStridedPointer,
    α, β, ::StaticInt{mc}, ::StaticInt{kc}, ::StaticInt{nc}, (M, K, N)# = matmul_sizes(C, A, B)
) where {T, mc, kc, nc}
    W = VectorizationBase.pick_vector_width_val(T)
    mᵣW = StaticInt{mᵣ}() * W

    num_n_iter = cld_fast(N, StaticInt{nc}())
    _nreps_per_iter = div_fast(N, num_n_iter) + StaticInt{nᵣ}() - One()
    nreps_per_iter = div_fast(_nreps_per_iter, StaticInt{nᵣ}()) * StaticInt{nᵣ}()
    Niter = num_n_iter - One()
    Nrem = N - Niter * nreps_per_iter

    num_m_iter = cld_fast(M, StaticInt{mc}())
    _mreps_per_iter = div_fast(M, num_m_iter) + mᵣW - One()
    mcrepetitions = div_fast(_mreps_per_iter, mᵣW)
    mreps_per_iter = mᵣW * mcrepetitions
    Miter = num_m_iter - One()
    Mrem = M - Miter * mreps_per_iter
    if (Mrem < mᵣW) & (Miter > 0)
        Miter -= one(Miter)
        Mrem += mreps_per_iter
    end

    num_k_iter = cld_fast(K, StaticInt{kc}())
    kreps_per_iter = ((div_fast(K, num_k_iter)) + 3) & -4
    Krem = K - (num_k_iter - One()) * kreps_per_iter
    Kiter = num_k_iter - One()

    Aptr = zero_offsets(A)
    Bptr = zero_offsets(B)
    Cptr = zero_offsets(C)
    # L3ptr = Base.unsafe_convert(Ptr{Tb}, pointer(BCACHE) + (Threads.threadid()-1)*BSIZE*8)
    # bc = _use_bcache()
    # L3ptr = Base.unsafe_convert(Ptr{T}, BCACHE)
    bcache = _use_bcache()
    L3ptr = Base.unsafe_convert(Ptr{T}, bcache)
    noffset = 0
    _Nrem, _nreps_per_iter = promote(Nrem, nreps_per_iter)
    _Mrem, _mreps_per_iter = promote(Mrem, mreps_per_iter)
    _Krem, _kreps_per_iter = promote(Krem, kreps_per_iter)
    # Cb = preserve_buffer(C); Ab = preserve_buffer(A); Bb = preserve_buffer(B);
    # GC.@preserve Cb Ab Bb BCACHE begin
    # @show T Niter, _Nrem, _nreps_per_iter Kiter, _Krem, _kreps_per_iter
    # return _free_bcache!(bcache)
    GC.@preserve BCACHE begin
        for no in 0:Niter
            # Krem
            # pack kc x nc block of B
            nsize = ifelse(no == Niter, _Nrem, _nreps_per_iter)
            # koffset = 0
            let Aptr = Aptr, Bptr = Bptr
                for ko ∈ 0:Kiter
                    ksize = ifelse(ko == 0, _Krem, _kreps_per_iter)
                    # _β = ifelse(ko == 0, convert(T, β), one(T))
                    # Bsubset2 = PtrArray(gesp(Bptr, (koffset, noffset)), (ksize, nsize), dense_dims_subset(dense_dims(B), stride_rank(B)))
                    # Bsubset2 = PtrArray(Bptr, (ksize, nsize), dense_dims_subset(dense_dims(B), stride_rank(B)))
                    Bsubset2 = PtrArray(Bptr, (ksize, nsize), none_dense(Val{2}()))
                    Bpacked2 = ptrarray0(L3ptr, (ksize, nsize))
                    # @show offsets(Bpacked2), offsets(Bsubset2)
                    copyto!(Bpacked2, Bsubset2)
                    let Aptr = Aptr, Cptr = Cptr, Bptr = zstridedpointer(Bpacked2)
                        # moffset = 0
                        for mo in 0:Miter
                            msize = ifelse(mo == Miter, _Mrem, _mreps_per_iter)
                            # pack mreps_per_iter x kreps_per_iter block of A
                            # Asubset = PtrArray(gesp(Aptr, (moffset, koffset)), (msize, ksize), dense_dims_subset(dense_dims(A), stride_rank(A)))

                            # Cpmat = PtrArray(gesp(Cptr, (moffset, noffset)), (msize, nsize), dense_dims_subset(dense_dims(C), stride_rank(C)))
                            # packaloopmul!(Cpmat, Asubset, Bpacked2, α, _β, (zrange(msize), zrange(ksize), zrange(nsize)))
                            if ko == 0
                                packaloopmul!(Cptr, Aptr, Bptr, α, β, (msize, ksize, nsize))
                            else
                                packaloopmul!(Cptr, Aptr, Bptr, α, One(), (msize, ksize, nsize))
                            end
                            # moffset += mreps_per_iter
                            Aptr = gesp(Aptr, (msize, Zero()))
                            Cptr = gesp(Cptr, (msize, Zero()))
                        end
                    end
                    # koffset += ksize
                    Aptr = gesp(Aptr, (Zero(), ksize))
                    Bptr = gesp(Bptr, (ksize, Zero()))
                end
            end
            # noffset += nreps_per_iter
            Bptr = gesp(Bptr, (Zero(), nsize))
            Cptr = gesp(Cptr, (Zero(), nsize))
        end
    end # GC.@preserve
    _free_bcache!(bcache)
    nothing
end

@inline contiguousstride1(A) = ArrayInterface.contiguous_axis(A) === ArrayInterface.Contiguous{1}()
@inline contiguousstride1(A::AbstractStridedPointer{T,N,1}) where {T,N} = true
@inline firstbytestride(A::AbstractStridedPointer) = VectorizationBase.bytestrides(A)[One()]
# @inline firstbytestride(A::AbstractStrideArray) = bytestride(A, One())
# @inline firstbytestride(A::PermutedDimsArray) = LinearAlgebra.stride1(A)
# @inline firstbytestride(A::Adjoint{<:Any,<:AbstractMatrix}) = stride(parent(A), 2)
# @inline firstbytestride(A::Transpose{<:Any,<:AbstractMatrix}) = stride(parent(A), 2)
# @inline firstbytestride(::Any) = typemax(Int)

@inline function vectormultiple(bytex, ::Type{Tc}, ::Type{Ta}) where {Tc,Ta}
    Wc = VectorizationBase.pick_vector_width_val(Tc) * static_sizeof(Ta) - One()
    iszero(bytex & (VectorizationBase.REGISTER_SIZE - 1))
end
@inline function dontpack(pA::AbstractStridedPointer{Ta}, M, K, ::StaticInt{mc}, ::StaticInt{kc}, ::Type{Tc}) where {mc, kc, Tc, Ta}
    (contiguousstride1(pA) &&
         ((((VectorizationBase.AVX512F ? 9 : 13) * VectorizationBase.pick_vector_width(Tc)) ≥ M) ||
          (vectormultiple(bytestride(pA, StaticInt{2}()), Tc, Ta) && ((M * K) ≤ (mc * kc)) && iszero(reinterpret(Int, pointer(pA)) & (VectorizationBase.REGISTER_SIZE - 1)))))
end

@inline function matmul_serial(A::AbstractMatrix, B::AbstractMatrix)
    m = size(A, StaticInt{1}())
    p = size(B, StaticInt{2}())
    C = StrideArray{promote_type(eltype(A),eltype(B))}(undef, (m,p))
    matmul_serial!(C, A, B)
    return C
end

@inline function matmul_serial!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
    maybeinline(C, A) && return inlineloopmul!(C, A, B, One(), Zero())
    matmul_serial!(C, A, B, One(), Zero(), ArrayInterface.is_column_major(C))
end
@inline function matmul_serial!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α)
    maybeinline(C, A) && return inlineloopmul!(C, A, B, α, Zero())
    matmul_serial!(C, A, B, α, Zero(), ArrayInterface.is_column_major(C))
end
@inline function matmul_serial!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α, β)
    maybeinline(C, A) && return inlineloopmul!(C, A, B, α, β)
    matmul_serial!(C, A, B, α, β, ArrayInterface.is_column_major(C))
end
@inline matmul_serial!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α, β, ::Val{false}) = (matmul_serial!(C', B', A', α, β, nothing); return C)
@inline matmul_serial!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α, β, ::Val{true}) = matmul_serial!(C, A, B, α, β, nothing)

"""
  matmul_serial!(C, A, B[, α = 1, β = 0])

Calculates `C = α * (A * B) + β * C` in place.

A single threaded matrix-matrix-multiply implementation.
Supports dynamically and statically sized arrays.

Organizationally, `matmul_serial!` checks the arrays properties to try and dispatch to an appropriate implementation.
If the arrays are small and statically sized, it will dispatch to an inlined multiply.

Otherwise, based on the array's size, whether they are transposed, and whether the columns are already aligned, it decides to not pack at all, to pack only `A`, or to pack both arrays `A` and `B`.
"""
@inline function matmul_serial!(
    C::AbstractMatrix{T}, A::AbstractMatrix, B::AbstractMatrix, α, β, MKN::Union{Nothing,Tuple{Vararg{Integer,3}}}
) where {T}
    pA = zstridedpointer(A); pB = zstridedpointer(B); pC = zstridedpointer(C);
    Cb = preserve_buffer(C); Ab = preserve_buffer(A); Bb = preserve_buffer(B);
    (M,K,N) = MKN === nothing ? matmul_sizes(C, A, B) : MKN
    Mc, Kc, Nc = matmul_params(T)
    GC.@preserve Cb Ab Bb begin
        if VectorizationBase.CACHE_SIZE[2] === nothing ||  (nᵣ ≥ N) || dontpack(pA, M, K, Mc, Kc, T)
            loopmul!(pC, pA, pB, α, β, (CloseOpen(M),CloseOpen(K),CloseOpen(N)))
            return C
        else
            jmul_singlethread_pack!(pC, pA, pB, α, β, Mc, Kc, Nc, (M,K,N))
            return C
        end
    end
    return C
end # function

function jmul_singlethread_pack!(pC, pA, pB, α, β, ::StaticInt{Mc}, ::StaticInt{Kc}, ::StaticInt{Nc}, (M,K,N)) where {Mc,Kc,Nc}
    if VectorizationBase.CACHE_SIZE[3] === nothing || (contiguousstride1(pB) ? (Kc * Nc ≥ K * N) : (firstbytestride(pB) ≤ 1600))
        # println("Pack A mul")
        jmulpackAonly!(pC, pA, pB, α, β, StaticInt{Mc}(), StaticInt{Kc}(), StaticInt{Nc}(), (M,K,N))
    else
        # println("Pack A and B mul")
        jmulpackAB!(pC, pA, pB, α, β, StaticInt{Mc}(), StaticInt{Kc}(), StaticInt{Nc}(), (M,K,N))
    end
    nothing
end

@inline function matmul!(
    C::AbstractStrideArray{S,D,T,2,2},
    A::AbstractMatrix,
    B::AbstractMatrix
) where {S, D, T}
    matmul!(C', B', A')
    C
end

# struct ThreadRun
#     id::UInt32
#     nthread::UInt32
# end
# ThreadRun(i::Int, n::Int) = ThreadRun(i % UInt32, n % UInt32)

struct LoopMulClosure{P,TC,TA,TB,Α,Β,Md,Kd,Nd}
    C::TC
    A::TA
    B::TB
    α::Α # \Alpha
    β::Β # \Beta
    M::Md
    K::Kd
    N::Nd
end
function LoopMulClosure{P}(
    C::TC, A::TA, B::TB, α::Α, β::Β, M::Md, K::Kd, N::Nd
) where {P,TC<:AbstractStridedPointer,TA<:AbstractStridedPointer,TB<:AbstractStridedPointer,Α,Β,Md,Kd,Nd}
    LoopMulClosure{P,TC,TA,TB,Α,Β,Md,Kd,Nd}(C, A, B, α, β, M, K, N)
end
# function LoopMulClosure{true}(
#     C::TC, A::TA, B::TB, α::Α, β::Β, Maxis::M, Kaxis::K, Naxis::N
# ) where {TC,TA,TB,Α,Β,M,K,N}
#     LoopMulClosure{true,TC,TA,TB,Α,Β,M,K,N}(C, A, B, α, β, Maxis, Kaxis, Naxis)
# end
function LoopMulClosure{P}(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α, β, M, K, N) where {P} # if not packing, discard `PtrArray` wrapper
    LoopMulClosure{P}(zstridedpointer(C), zstridedpointer(A), zstridedpointer(B), α, β, M, K, N)
end
# function LoopMulClosure{true}(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α, β, M, K, N)
#     LoopMulClosure{true}(C, A, B, α, β, M, K, N)
# end
(m::LoopMulClosure{false})() = loopmul!(m.C, m.A, m.B, m.α, m.β, (CloseOpen(m.M), CloseOpen(m.K), CloseOpen(m.N)))
# (m::LoopMulClosure{true})() = jmulpackAonly!(m.C, m.A, m.B, m.α, m.β, (m.Maxis, m.Kaxis, m.Naxis))
function (m::LoopMulClosure{true,TC})() where {T,TC <: AbstractStridedPointer{T}}
    Mc,Kc,Nc = matmul_params(T)
    jmulpackAonly!(m.C, m.A, m.B, m.α, m.β, Mc, Kc, Nc, (m.M, m.K, m.N))
end

# struct PackAClosure{TC,TA,TB,Α,Β,Md,Kd,Nd,T}
#     # lmc::LoopMulClosure{TC,TA,TB,Α,Β,M,K,N}
#     C::TC
#     A::TA
#     B::TB
#     α::Α # \Alpha
#     β::Β # \Beta
#     M::Md
#     K::Kd
#     N::Nd
#     tasks::T
# end
# # function PackAClosure(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, α, β, M, K, N, tasks)
# #     PackAClosure(stridedpointer(C), stridedpointer(A), stridedpointer(B), α, β, M, K, N, tasks)
# # end
# function (m::PackAClosure{TC})() where {T,TC<:AbstractStridedPointer{T}}
#     Mc,Kc,Nc = matmul_params(T)
#     jmultpackAonly!(m.C, m.A, m.B, m.α, m.β, Mc, Kc, Nc, m.M, m.K, m.N, m.tasks)
# end

"""
    matmul(A, B)

Multiply matrices `A` and `B`.
"""
@inline function matmul(A, B)
    m = size(A, StaticInt{1}())
    p = size(B, StaticInt{2}())
    C = StrideArray{promote_type(eltype(A),eltype(B))}(undef, (m,p))
    matmul!(C, A, B)
    return C
end

"""
    matmul!(C, A, B[, α, β, max_threads])

Calculates `C = α * A * B + β * C` in place, overwriting the contents of `A`.
It may use up to `max_threads` threads. It will not use threads when nested in other threaded code.
"""
@inline function matmul!(C, A, B)
    maybeinline(C, A) && return inlineloopmul!(C, A, B, One(), Zero())
    matmul!(C, A, B, One(), Zero(), nothing, ArrayInterface.is_column_major(C))
end
@inline function matmul!(C, A, B, α)
    maybeinline(C, A) && return inlineloopmul!(C, A, B, α, Zero())
    matmul!(C, A, B, α, Zero(), nothing, ArrayInterface.is_column_major(C))
end
@inline function matmul!(C, A, B, α, β)
    maybeinline(C, A) && return inlineloopmul!(C, A, B, α, β)
    matmul!(C, A, B, α, β, nothing, ArrayInterface.is_column_major(C))
end
@inline function matmul!(C, A, B, α, β, nthread)
    maybeinline(C, A) && return inlineloopmul!(C, A, B, α, β)
    matmul!(C, A, B, α, β, nthread, ArrayInterface.is_column_major(C))
end
@inline matmul!(C::AbstractMatrix, A, B, α, β, nthread, ::Val{false}) = (matmul!(C', B', A', α, β, nthread, nothing); return C)
@inline matmul!(C::AbstractMatrix, A, B, α, β, nthread, ::Val{true}) = matmul!(C, A, B, α, β, nthread, nothing)

@inline function dontpack(pA::AbstractStridedPointer{Ta}, M, K, ::StaticInt{mc}, ::StaticInt{kc}, ::Type{Tc}, nspawn) where {mc, kc, Tc, Ta}
    # MᵣW = VectorizationBase.pick_vector_width_val(Tc) * StaticInt{mᵣ}()
    # TODO: perhaps consider K vs kc by themselves?
    (contiguousstride1(pA) && ((M * K) ≤ (mc * kc) * nspawn >>> 1))
end

@inline function matmul!(C::AbstractMatrix{T}, A, B, α, β, nthread, matmuldims) where {T}#::Union{Nothing,Tuple{Vararg{Integer,3}}}) where {T}
    # Don't allow nested `jmult!`, and if we were padded `nthread` of 1, don't thread

    M, K, N = matmuldims === nothing ? matmul_sizes(C, A, B) : matmuldims
    W = VectorizationBase.pick_vector_width_val(T)
    # Assume 3 μs / spawn cost
    # 15W approximate GFLOPS target
    # 22500 = 3e-6μ/spawn * 15 / 2e-9
    MKN = M*K*N
    Mc, Kc, Nc = matmul_params(T)

    pA = zstridedpointer(A); pB = zstridedpointer(B); pC = zstridedpointer(C);
    Cb = preserve_buffer(C); Ab = preserve_buffer(A); Bb = preserve_buffer(B);
    N_too_small = nᵣ ≥ N
    GC.@preserve Cb Ab Bb begin
        # if N_too_small | (MKN < (StaticInt{5451776}() ÷ static_sizeof(T)))
        if N_too_small | (MKN < (StaticInt{110592}() * W))
            # Only start threading beyond about 88x88 * 88x88 for `Float64`, or 111x111 * 111x111 for `Float32`
            # This is a small-matrix fast path
            if N_too_small | dontpack(pA, M, K, Mc, Kc, T)
                loopmul!(pC, pA, pB, α, β, (CloseOpen(M),CloseOpen(K),CloseOpen(N)))
                return C
            else
                jmulpackAonly!(pC, pA, pB, α, β, Mc, Kc, Nc, (M,K,N))
                return C
            end
        end
        # Not taking the fast path
        # But maybe we don't want to thread anyway
        nt = _nthreads()
        _nthread = nthread === nothing ? nt : min(_nthreads(), nthread)
        # Maybe this is nested, or we have ≤ 1 threads
        if (!iszero(ccall(:jl_in_threaded_region, Cint, ())) || (_nthread ≤ 1))
            jmul_singlethread_pack!(pC, pA, pB, α, β, Mc, Kc, Nc, (M,K,N))
            return C
        end
        # We are threading, but how many threads?
        # nspawn = min(_nthread, div_fast(M + N, StaticInt{4}() * W))
        # nspawn = clamp(div_fast(M + N, StaticInt{4}() * W), 1, NUM_CORES)
        nspawn = clamp(div_fast((M + N)*StaticInt{11}(), StaticInt{64}() * W), 1, NUM_CORES)
        # L = StaticInt{22500}() * W
        # nspawn = min(_nthread, cld_fast(MKN, L))
        _matmul!(pC, pA, pB, α, β, nspawn, (M,K,N), Mc, Kc, Nc)
    end
    return C
end

@inline function gcd_fast(a::T, b::T) where {T<:Base.BitInteger}
    za = trailing_zeros(a)
    zb = trailing_zeros(b)
    k = min(za, zb)
    u = unsigned(abs(a >> za))
    v = unsigned(abs(b >> zb))
    while u != v
        if u > v
            u, v = v, u
        end
        v -= u
        v >>= trailing_zeros(v)
    end
    r = u << k
    r % T
end
# @inline function divide_blocks(M, N, ::StaticInt{Mb}, ::StaticInt{Nb}, nspawn) where {Mb,Nb}
#     Mfull, Mrem = divrem_fast(M, Mb)
#     Mtotal = Mfull + (Mrem > 0)

#     Nfull, Nrem = divrem_fast(N, Nb)
#     Ntotal = Nfull + (Nrem > 0)

#     divide_blocks(Mtotal, Ntotal, nspawn)
# end

@inline function find_first_acceptable(M, ::StaticInt{W}) where {W}
    Mr = StaticInt{mᵣ}() * StaticInt{W}()
    for (miter,niter) ∈ CORE_FACTORS
        if miter * ((MᵣW_mul_factor - One()) * Mr) ≤ M + StaticInt{2}() * StaticInt{W}()
        # if miter * ((MᵣW_mul_factor) * Mr) ≤ M
            return miter, niter
        end
    end
    last(CORE_FACTORS)
end

# function divide_blocks(M, Mblocks, W)
#     d, r = divrem_fast(cld_fast(M,W), Mblocks)
#     # if 2r > Mblocks
#     Mblock_m_1 = Mblocks - 1
#     rlarge = 2r ≥ Mblocks - 1
#     bsize = W * (d + rlarge)
#     bsize, M - bsize*Mblock_m_1, rlarge
#     # bsize_high = (d + (r>0))*W
#     # (bsize_low, M - bsize_low*(Mblocks-1)), (bsize_high, M - bsize_high*(Mblocks-1)), r
# end

@inline function divide_blocks(M, Ntotal, _nspawn, ::StaticInt{W}) where {W}
    _nspawn == NUM_CORES && return find_first_acceptable(M, StaticInt{W}())

    Miter = clamp(div_fast(M, StaticInt{W}()*StaticInt{mᵣ}() * MᵣW_mul_factor), 1, _nspawn)
    # Miter = clamp(cld(M, Mc), 1, _nspawn)
    nspawn = div_fast(_nspawn, Miter)
    Niter = if (nspawn ≤ 1) & (Miter < _nspawn)
        # rebalance Miter
        Miter = cld_fast(_nspawn, cld_fast(_nspawn, Miter))
        nspawn = div_fast(_nspawn, Miter)
    end
    Miter, cld_fast(Ntotal, max(2, cld_fast(Ntotal, nspawn)))
    # # Niter = cld_fast(Ntotal, max(2, cld_fast(Ntotal, nspawn)))
    # # Niter = cld_fast(Ntotal, cld_fast(Ntotal, _nspawn))
    # return Miter, Niter
    # nspawn = div_fast(_nspawn, Niter)
    # Miter = cld_fast(Mtotal, cld_fast(Mtotal, nspawn))
    # return Miter, Niter
    # Miter = gcd_fast(_nspawn, Mtotal)
    # nspawn = div_fast(_nspawn, Miter)
    # Niter = cld_fast(Ntotal, cld_fast(Ntotal, nspawn))
    # # Niter = gcd_fast(nspawn, Ntotal)
    # Niter = cld_fast(Ntotal, cld_fast(Ntotal, nspawn))
    # return Miter, Niter
end

@inline function split_m(M, _Mblocks, ::StaticInt{W}) where {W}
    Miters = cld_fast(M, StaticInt{W}())
    Mblocks = min(_Mblocks, Miters)
    Miter_per_block, Mrem = divrem_fast(Miters, Mblocks)
    Mbsize = Miter_per_block * StaticInt{W}()
    Mremfinal = M - Mbsize*(Mblocks-1) - Mrem * W
    Mbsize, Mrem, Mremfinal, Mblocks
end

function jmultsplitn!(C::AbstractStridedPointer{T}, A, B, α, β, ::StaticInt{Mc}, nspawn, (M,K,N), ::Val{PACK}) where {T, Mc, PACK}
    Mᵣ = StaticInt{mᵣ}(); Nᵣ = StaticInt{nᵣ}();
    W = VectorizationBase.pick_vector_width_val(T)
    MᵣW = Mᵣ*W

    _Mblocks, _Nblocks = divide_blocks(M, cld_fast(N, Nᵣ), nspawn, W)
    Mbsize, Mrem, Mremfinal, Mblocks = split_m(M, _Mblocks, W)
    Nblocks = min(N, _Nblocks)
    Nbsize, Nrem = divrem(N, Nblocks)

    _nspawn = Mblocks * Nblocks
    _nspawn_m1 = _nspawn - One()
    Mbsize_Mrem, Mbsize_ = promote(Mbsize + W, Mbsize)
    Nbsize_Nrem, Nbsize_ = promote(Nbsize + One(), Nbsize)

    tnum = 0
    # @show Mblocks, Nblocks, (Mbsize,Mrem), (Nbsize,Nrem)
    let _A = A, _B = B, _C = C
        for n ∈ Base.OneTo(Nblocks)
            nsize = ifelse(n ≤ Nrem, Nbsize_Nrem, Nbsize_)
            let _A = _A, _C = _C
                for m ∈ Base.OneTo(Mblocks)
                    msize = ifelse(m == Mblocks, Mremfinal, ifelse(m ≤ Mrem, Mbsize_Mrem, Mbsize_))
                    if tnum == _nspawn_m1
                        wait(runfunc(LoopMulClosure{PACK}(_C, _A, _B, α, β, msize, K, nsize), _nspawn))
                        waitonmultasks(CloseOpen(_nspawn))
                        return
                    end
                    runfunc!(LoopMulClosure{PACK}(_C, _A, _B, α, β, msize, K, nsize), (tnum += 1))
                    _A = gesp(_A, (msize, Zero()))
                    _C = gesp(_C, (msize, Zero()))
                end
            end
            _B = gesp(_B, (Zero(), nsize))
            _C = gesp(_C, (Zero(), nsize))
        end
    end
end
# function jmultsplitn!(C::AbstractStridedPointer{T}, A, B, α, β, nspawn, (M,K,N), ::StaticInt{Mc}, ::StaticInt{Kc}, ::StaticInt{Nc}) where {T, Mc, Kc, Nc}
#     Mᵣ = StaticInt{mᵣ}(); Nᵣ = StaticInt{nᵣ}();
#     W = VectorizationBase.pick_vector_width_val(T)
#     MᵣW = Mᵣ*W

#     # _Niter = cld_fast(nspawn, N)
#     # _Nper = max(cld_fast(N, _Niter), 2nᵣ)
#     # __Niter, _Nrem = divrem_fast(N, _Nper)
#     # __Niter += (_Nrem > 0)

#     Nchunks = cld_fast(N, StaticInt{2}() * Nᵣ)
#     nchunks_per_spawn = cld_fast(Nchunks, nspawn)

#     Nper = nchunks_per_spawn * (StaticInt{2}() * Nᵣ)
#     Niter, Nrem = divrem_fast(N, Nper)
#     if iszero(Nrem)
#         Nrem = Nper
#     else
#         Niter += 1
#     end

#     # Nper = cld_fast(N, cld_fast(N, __Niter))
#     # # Nrem =

#     # _Mtotal = cld_fast(M, MᵣW)
#     # # split_n = nspawn > _Mtotal
#     # _Niter = cld_fast(nspawn, _Mtotal)
#     # _Nper = cld(N, _Niter)
#     # @show Niter, nspawn, N, Nper, Nrem
#     nspawn_per = div_fast(nspawn, Niter)
#     spawn_start = 0#spawn_per_mrep
#     # error("")
#     # Nper = _Nper
#     # Nrem = N - (Nper * (_Niter - 1))
#     for i ∈ Base.OneTo(Niter)
#         Nsize = ifelse(i == Niter, Nrem, Nper)
#         next_start = ifelse(i == Niter, nspawn, spawn_start + nspawn_per)
#         task_view = CloseOpen(spawn_start, next_start)
#         jmultpackAonly!(C, A, B, α, β, StaticInt{Mc}(), StaticInt{Kc}(), StaticInt{Nc}(), M, K, Nsize, task_view)
#         spawn_start = next_start
#         B = gesp(B, (Zero(), Nper))
#         C = gesp(C, (Zero(), Nper))
#     end
#     # task_view = CloseOpen(spawn_start, nspawn)
#     # jmultpackAonly!(C, A, B, α, β, StaticInt{Mc}(), StaticInt{Kc}(), StaticInt{Nc}(), M, K, Nrem, task_view)
#     nothing
# end

function _matmul!(
    C::AbstractStridedPointer{T}, A::AbstractStridedPointer, B::AbstractStridedPointer,
    α, β, nspawn, (M,K,N), ::StaticInt{Mc}, ::StaticInt{Kc}, ::StaticInt{Nc}
) where {T,Mc,Kc,Nc}
    Mᵣ = StaticInt{mᵣ}(); Nᵣ = StaticInt{nᵣ}();
    W = VectorizationBase.pick_vector_width_val(T)
    MᵣW = Mᵣ*W
    # nkern = cld_fast(M * N,  MᵣW * Nᵣ)
    # @show nspawn
    # @show VectorizationBase.CACHE_SIZE[2] === nothing ||  (nᵣ ≥ N) ||
    #     dontpack(A, M, K, StaticInt{Mc}(), StaticInt{Kc}(), T) || # do not pack A
    #     (N * M ≤ (Nc * MᵣW) * nspawn) # do we just want to split aggressively?
    # @show (contiguousstride1(B) && (Kc * Nc ≥ K * N)) | (firstbytestride(B) > 1600) # do not pack B
    # @show M + MᵣW ≤ MᵣW * nspawn # do we want to split `N`
    # Approach:
    # Check if we don't want to pack A,
    #    if not, aggressively subdivide
    # if so, check if we don't want to pack B
    #    if not, check if we want to thread `N` loop anyway
    #       if so, divide `M` first, then use ratio of desired divisions / divisions along `M` to calc divisions along `N`
    #       if not, only thread along `M`. These don't need syncing, as we're not packing `B`
    #    if so, `jmultpackAB!`

    # MᵣW * (MᵣW_mul_factor - One()) # gives a smaller Mc, then
    # if 2M/nspawn is less than it, we don't don't `A`
    if VectorizationBase.CACHE_SIZE[2] === nothing ||  # do not pack A
        dontpack(A, M, K, StaticInt{Mc}(), StaticInt{Kc}(), T, nspawn) || (nᵣ ≥ N) || (W ≥ M) # do we just want to split aggressively?
        jmultsplitn!(C, A, B, α, β, StaticInt{Mc}(), nspawn, (M,K,N), Val{false}())
    elseif contiguousstride1(B) ? ((Kc * Nc - StaticInt{200_000}()) ≥ K * N) : (firstbytestride(B) ≤ 1600)
        # if M + MᵣW ≤ MᵣW * nspawn # do we want to split `N`
        # if 2nᵣ * nspawn ≤ N
            # jmultsplitn!(C, A, B, α, β, nspawn, (M,K,N), Val{true}())
        # else
        jmultsplitn!(C, A, B, α, β, StaticInt{Mc}(), nspawn, (M,K,N), Val{true}())
            # jmultpackAonly!(C, A, B, α, β, StaticInt{Mc}(), StaticInt{Kc}(), StaticInt{Nc}(), M, K, N, CloseOpen(0, nspawn))
        # end
    else
        jmultpackAB!(C, A, B, α, β, StaticInt{Mc}(), StaticInt{Kc}(), StaticInt{Nc}(), M, K, N, CloseOpen(nspawn), Val{CACHE_COUNT[3]}())
    end
end

# If tasks is [0,1,2,3] (e.g., `CloseOpen(0,4)`), it will wait on `MULTASKS[i]` for `i = [1,2,3]`.
function waitonmultasks(tasks::CloseOpen)
    tid = Int(first(tasks))
    while (tid += 1) < tasks.upper
        @inbounds wait(MULTASKS[tid])
    end
end

# _front(r::AbstractUnitRange) = CloseOpen(first(r), last(r))
# """
# Runs on threadids `tasks .+ 1`
# """
# function jmultpackAonly!(
#     C::AbstractStridedPointer{T}, A::AbstractStridedPointer, B::AbstractStridedPointer,
#     α, β, ::StaticInt{Mc}, ::StaticInt{Kc}, ::StaticInt{Nc}, M, K, N, tasks::CloseOpen
# ) where {T,Mc,Kc,Nc}
#     # e.g., tasks = CloseOpen(0,4) = [0,1,2,3] runs on [1,2,3,4]
#     to_spawn = length(tasks)

#     W = VectorizationBase.pick_vector_width_val(T)
#     _Mblock = cld_fast(M, to_spawn)
#     Mblock = cld_fast(_Mblock, W) * W
#     _Miter = M ÷ Mblock
#     _Mrem = M - Mblock * _Miter
#     if iszero(_Mrem)
#         Mrem = Mblock
#         Miter = _Miter - One()
#     else
#         Mrem = _Mrem
#         Miter = _Miter
#     end
#     tid = ft = first(tasks)
#     # @show Miter, Mrem, Mblock, to_spawn
#     # error("")
#     for _ ∈ Base.OneTo(Miter)
#         runfunc!(LoopMulClosure{true}(C, A, B, α, β, Mblock, K, N), (tid += 1))
#         A = gesp(A, (Mblock, Zero()))
#         C = gesp(C, (Mblock, Zero()))
#     end
#     wait(runfunc(LoopMulClosure{true}(C, A, B, α, β, Mrem, K, N), (tid += 1)))
#     # jmulpackAonly!(C, A, B, α, β, StaticInt{Mc}(), StaticInt{Kc}(), StaticInt{Nc}(), (Mrem, K, N))
#     waitonmultasks(CloseOpen(ft, tid))
# end

# open_upper(r::CloseOpen) = r.upper
function jmultpackAB!(
    C::AbstractStridedPointer{T}, A::AbstractStridedPointer, B::AbstractStridedPointer,
    α, β, ::StaticInt{Mc}, ::StaticInt{Kc}, ::StaticInt{Nc}, M, K, N, tasks::CloseOpen, ::Val{CC}
) where {T,Mc,Kc,Nc,CC}
    @assert CC > 1

    W = VectorizationBase.pick_vector_width_val(T)
    mᵣW = StaticInt{mᵣ}() * W

    to_spawn = length(tasks)
    # to_spawn = to_spawn + 1
    # nspawn = length(task_view) + 1
    # Nsplits = cld_fast(cld_fast(M, mᵣW), nspawn)
    Nsplits = min(cld_fast(M, mᵣW * to_spawn), CC)
    Nsplits > 1 || return jmultpackAB!(C, A, B, α, β, StaticInt{Mc}(), StaticInt{Kc}(), StaticInt{Nc}(), M, K, N, tasks, Val{1}())

    _Nsize, Nrem = divrem_fast(N, Nsplits)
    Nsize_Nrem, Nsize_ = promote(Nsize + One(), Nsize)

    Nspawn_per = cld_fast(Nsplits, to_spawn)
    task_start = 0
    for i ∈ CloseOpen(Nsplits)
        task_next = min(task_start + Nspawn_per, tasks.upper)
        _taskview = CloseOpen(task_start, task_next)
        task_start = task_next
        nsize = ifelse(i < Nrem, Nsize_Nrem, Nsize_)
        jmultpackAB!(C, A, B, α, β, StaticInt{Mc}(), StaticInt{Kc}(), StaticInt{Nc}(), M, K, nsize, _taskview, Val{1}())
        B = gesp(C, (Zero(), nsize))
        C = gesp(C, (Zero(), nsize))
    end
end
function jmultpackAB!(
    C::AbstractStridedPointer{T}, A::AbstractStridedPointer, B::AbstractStridedPointer,
    α, β, ::StaticInt{Mc}, ::StaticInt{Kc}, ::StaticInt{Nc}, M, K, N, tasks::CloseOpen, ::Val{1}
) where {T,Mc,Kc,Nc}
    W = VectorizationBase.pick_vector_width_val(T)
    mᵣW = StaticInt{mᵣ}() * W

    to_spawn = length(tasks)
    atomicsync = Ref{NTuple{9,UInt}}()
    p = Base.unsafe_convert(Ptr{UInt}, atomicsync)
    _atomic_min!(p, zero(UInt)); _atomic_min!(p + 8sizeof(UInt), zero(UInt))
    # unsafe_store!(p, zero(UInt), 1); unsafe_store!(p, zero(UInt), 9)
    Mbsize, Mrem, Mremfinal, _to_spawn = split_m(M, to_spawn, W) # M is guaranteed to be > W because of `W ≥ M` condition for `jmultsplitn!`...
    Mblock_Mrem, Mblock_ = promote(Mbsize + W, Mbsize)
    u_to_spawn = _to_spawn % UInt
    ft = first(tasks)
    tid = Int(ft)
    # @show Mbsize, Mrem, Mremfinal, _to_spawn
    # return

    bc = _use_bcache()
    bc_ptr = Base.unsafe_convert(typeof(pointer(C)), pointer(bc))
    # Mblock = cld_fast(cld_fast(M, to_spawn), MᵣW) * Mᵣ*W
    GC.@preserve atomicsync begin
        for i ∈ CloseOpen(One(), _to_spawn) # ...thus the fact that `CloseOpen()` iterates at least once is okay.
            Mblock = ifelse(i ≤ Mrem, Mblock_Mrem, Mblock_)
            runfunc!(SyncClosure{Mc,Kc,Nc}(C, A, B, α, β, Mblock, K, N, p, bc_ptr, i % UInt, u_to_spawn), (tid += 1))
            A = gesp(A, (Mblock, Zero()))
            C = gesp(C, (Mblock, Zero()))
        end
        wait(runfunc(SyncClosure{Mc,Kc,Nc}(C, A, B, α, β, Mremfinal, K, N, p, bc_ptr, zero(UInt), u_to_spawn), (tid += 1)))
        waitonmultasks(CloseOpen(ft, tid))
    end
    _free_bcache!(bc)
    return
end
function jmultpackAB!(
    C::AbstractStridedPointer{T}, A::AbstractStridedPointer, B::AbstractStridedPointer,
    α, β, ::Val{Mc}, ::Val{Kc}, ::Val{Nc}, M, K, N, tasks::CloseOpen, ::Val{1}
) where {T,Mc,Kc,Nc}
    W = VectorizationBase.pick_vector_width_val(T)
    mᵣW = StaticInt{mᵣ}() * W

    to_spawn = length(tasks)
    atomicsync = Ref{NTuple{9,UInt}}()
    p = Base.unsafe_convert(Ptr{UInt}, atomicsync)
    _atomic_min!(p, zero(UInt)); _atomic_min!(p + 8sizeof(UInt), zero(UInt))
    # unsafe_store!(p, zero(UInt), 1); unsafe_store!(p, zero(UInt), 9)
    Mbsize, Mrem, Mremfinal, _to_spawn = split_m(M, to_spawn, W) # M is guaranteed to be > W because of `W ≥ M` condition for `jmultsplitn!`...
    Mblock_Mrem, Mblock_ = promote(Mbsize + W, Mbsize)
    Mblock_Mrem = ifelse(Mrem > 0, Mblock_Mrem, Mblock_)
    u_to_spawn = _to_spawn % UInt
    tid = Int(first(tasks))
    # @show Mbsize, Mrem, Mremfinal, _to_spawn
    # return
    bc = _use_bcache()
    bc_ptr = Base.unsafe_convert(typeof(pointer(C)), pointer(bc))
    # Mblock = cld_fast(cld_fast(M, to_spawn), MᵣW) * Mᵣ*W
    GC.@preserve atomicsync begin
        for i ∈ CloseOpen(One(), _to_spawn) # ...thus the fact that `CloseOpen()` iterates at least once is okay.
            runfunc!(SyncClosure{Mc,Kc,Nc}(C, A, B, α, β, ((i ≤ Mrem) % Int,(Mblock_,Mblock_Mrem,Mremfinal)), K, N, p, bc_ptr, i % UInt, u_to_spawn), (tid += 1))
            Mblock = ifelse(i ≤ Mrem, Mblock_Mrem, Mblock_)
            A = gesp(A, (Mblock, Zero()))
            C = gesp(C, (Mblock, Zero()))
        end
        wait(runfunc(SyncClosure{Mc,Kc,Nc}(C, A, B, α, β, (2, (Mblock_,Mblock_Mrem,Mremfinal)), K, N, p, bc_ptr, zero(UInt), u_to_spawn), (tid += 1)))
        waitonmultasks(CloseOpen(first(tasks), tid))
    end
    _free_bcache!(bc)
    return
end

# struct PackABClosure{Mc,Kc,Nc,TC,TA,TB,Α,Β,Md,Kd,Nd,V}
#     C::TC
#     A::TA
#     B::TB
#     α::Α
#     β::Β
#     M::Md
#     K::Kd
#     N::Nd
#     tv::V
# end
# function PackABClosure{Mc,Kc,Nc}(
#     C::TC, A::TA, B::TB, α::Α, β::Β, ::StaticInt{Mc}, ::StaticInt{Kc}, ::StaticInt{Nc}, M::Md, K::Kd, N::Nd, _taskview::V
# ) where {Mc,Kc,Nc,TC,TA,TB,Α,Β,Md,Kd,Nd,V}
#     PackABClosure{Mc,Kc,Nc,TC,TA,TB,Α,Β,Md,Kd,Nd,V}(C, A, B, α, β, M, K, N, _taskview)
# end
# function (m::PackABClosure{Mc,Kc,Nc})() where {Mc,Kc,Nc}
#     jmultpackAB!(m.C, m.A, m.B, m.α, m.β, StaticInt{Mc}(), StaticInt{Kc}(), StaticInt{Nc}(), m.M, m.K, m.N, m.tv, Val{1}())
# end

struct SyncClosure{Mc,Kc,Nc,TC,TA,TB,Α,Β,Md,Kd,Nd,CA}
    C::TC
    A::TA
    B::TB
    α::Α
    β::Β
    M::Md
    K::Kd
    N::Nd
    p::Ptr{UInt}
    bc::Ptr{CA}
    id::UInt
    last_id::UInt
end

function SyncClosure{Mc,Kc,Nc}(
    C::TC, A::TA, B::TB, α::Α, β::Β, M::Md, K::Kd, N::Nd, p::Ptr{UInt}, bc::Ptr{CA}, id::UInt, last_id::UInt
) where {Mc,Kc,Nc,TC,TA,TB,Α,Β,Md,Kd,Nd,CA}
    SyncClosure{Mc,Kc,Nc,TC,TA,TB,Α,Β,Md,Kd,Nd,CA}(C, A, B, α, β, M, K, N, p, bc, id, last_id)
end

function (sc::SyncClosure{Mc,Kc,Nc,TC,TA,TB,Α,Β,Md})() where {Mc,Kc,Nc,TC,TA,TB,Α,Β,Md<:Tuple{Int,NTuple{3,Int}}}
    sync_mul!(sc.C, sc.A, sc.B, sc.α, sc.β, Val{Mc}(), Val{Kc}(), Val{Nc}(), sc.M, sc.K, sc.N, sc.p, sc.bc, sc.id, sc.last_id)
end
function (sc::SyncClosure{Mc,Kc,Nc})() where {Mc,Kc,Nc}
    sync_mul!(sc.C, sc.A, sc.B, sc.α, sc.β, StaticInt{Mc}(), StaticInt{Kc}(), StaticInt{Nc}(), sc.M, sc.K, sc.N, sc.p, sc.bc, sc.id, sc.last_id)
end

function pick_block(M, ::Val{N}, ::Val{R}, ::Type{T}) where {N, R, T}
    Base.fptosi(Int, (effective_cache(T, Val{N}()) * R) / M)
end
# function block_sizes(M, K, N, ::Val{Mc}, ::Val{Kc}, ::Val{Nc}) where {Mc,Kc,Nc}

# end

function sync_mul!(
    C::AbstractStridedPointer{T}, A::AbstractStridedPointer, B::AbstractStridedPointer,
    α, β, ::Val{Mc}, ::Val{Kc}, ::Val{Nc}, (Mid,Ms), K, N, atomicp::Ptr{UInt}, bc::Ptr, id::UInt, total_ids::UInt
) where {T, Mc, Kc, Nc}

    W = VectorizationBase.pick_vector_width_val(T)
    Base.Cartesian.@nexprs 3 i -> begin
        _M_i = Ms[i]
        _Miter_i = cld_fast(_M_i, StaticInt{Mc}())
        _Mbsize_i, _Mrem_i, _Mremfinal_i, _to_spawn_i = split_m(_M_i, _Miter_i, W) # M is guaranteed to be > W because of `W ≥ M` condition for `jmultsplitn!`...
        _Mblock_Mrem_i, _Mblock__i = promote(_Mbsize_i + W, _Mbsize_i)
        _Mblock_Mrem_i = ifelse(_Mrem_i > 0, _Mblock_Mrem_i, _Mblock__i)
    end
    biggest_mb = Base.Cartesian.@ncall 3 max i -> _Mblock_Mrem_i
    Base.Cartesian.@nif 3 i -> ((i - 1) == Mid) i -> begin
        Miter = _Miter_i
        Mrem = _Mrem_i
        Mremfinal = _Mremfinal_i
        Mblock_Mrem = _Mblock_Mrem_i
        Mblock_ = _Mblock__i
    end i -> begin
        Miter = _Miter_i
        Mrem = _Mrem_i
        Mremfinal = _Mremfinal_i
        Mblock_Mrem = _Mblock_Mrem_i
        Mblock_ = _Mblock__i
    end
    # M = @inbounds(Ms[Mid+1])
    # everyone will use the largest mb to pick values for `kc` and `nc`
    kc = pick_block(biggest_mb, Val{2}(), Val{Kc}(), T)
    Kiter = cld_fast(K, kc)
    Kblockdiv, Krem = divrem_fast(K, Kiter)
    Kblock_Krem, Kblock_ = promote(Kblockdiv + One(), Kblockdiv)

    nc = pick_block(kc, Val{3}(), Val{Nc}(), T)
    Niter = cld_fast(N, nc)
    Nblockdiv, Nrem = divrem_fast(N, Niter)
    Nblock_Nrem, Nblock_ = promote(Nblockdiv + One(), Nblockdiv)

    # @show biggest_mb, (Mblock_Mrem, Mblock_), (Miter, Mrem), kc, (Kblock_Krem, Kblock_), (Kiter, Krem), nc, (Nblock_Nrem, Nblock_), (Niter, Nrem)

    last_id = total_ids - one(UInt)
    atomics = atomicp + 8sizeof(UInt)
    sync_iters = zero(UInt64)

    Npackb_r_div, Npackb_r_rem = divrem_fast(Nblock_Nrem, total_ids)
    Npackb_r_block_rem, Npackb_r_block_ = promote(Npackb_r_div + One(), Npackb_r_div)

    Npackb___div, Npackb___rem = divrem_fast(Nblock_, total_ids)
    Npackb___block_rem, Npackb___block_ = promote(Npackb___div + One(), Npackb___div)

    pack_r_offset = Npackb_r_div * id + min(id, Npackb_r_rem)
    pack_r_view = CloseOpen(pack_r_offset, pack_r_offset + ifelse(id < Npackb_r_rem, Npackb_r_block_rem, Npackb_r_block_))
    pack___offset = Npackb___div * id + min(id, Npackb___rem)
    pack___view = CloseOpen(pack___offset, pack___offset + ifelse(id < Npackb___rem, Npackb___block_rem, Npackb___block_))

    GC.@preserve BCACHE begin
        for no in CloseOpen(Niter)
            # Krem
            # pack kc x nc block of B
            nfull = no < Nrem
            nsize = ifelse(nfull, Nblock_Nrem, Nblock_)
            pack_view = nfull ? pack_r_view : pack___view
            let A = A, B = B#, C = C
                for ko ∈ CloseOpen(Kiter)
                    ksize = ifelse(ko < Krem, Kblock_Krem, Kblock_)

                    Bsubset2 = PtrArray(B, (ksize, nsize), none_dense(Val{2}()))
                    Bpacked2 = ptrarray0(bc, (ksize, nsize))
                    copyto!(zview(Bpacked2, :, pack_view), zview(Bsubset2, :, pack_view))
                    # synchronize before starting the multiplication, to ensure `B` is packed
                    sync_iters += total_ids
                    _mv = _atomic_add!(atomicp, one(UInt))
                    while _mv < sync_iters
                        pause()
                        _mv = _atomic_max!(atomicp, zero(UInt))
                    end
                    # multiply
                    let A = A, B = zstridedpointer(Bpacked2), C = C
                        for mo in CloseOpen(Miter)
                            # @show mo, ko, no
                            msize = ifelse(mo == (Miter-1), Mremfinal, ifelse(mo < Mrem, Mblock_Mrem, Mblock_))
                            if ko == 0
                                packaloopmul!(C, A, B, α,     β, (msize, ksize, nsize))
                            else
                                packaloopmul!(C, A, B, α, One(), (msize, ksize, nsize))
                            end
                            A = gesp(A, (msize, Zero()))
                            C = gesp(C, (msize, Zero()))
                        end
                    end
                    A = gesp(A, (Zero(), ksize))
                    B = gesp(B, (ksize, Zero()))
                    # synchronize on completion of using `Bpacked`
                    _mv = _atomic_add!(atomics, one(UInt))
                    while _mv < sync_iters
                        pause()
                        _mv = _atomic_max!(atomics, zero(UInt))
                    end
                end
            end
            B = gesp(B, (Zero(), nsize))
            C = gesp(C, (Zero(), nsize))
        end
    end # GC.@preserve
    nothing
end
function sync_mul!(
    C::AbstractStridedPointer{T}, A::AbstractStridedPointer, B::AbstractStridedPointer,
    α, β, ::StaticInt{Mc}, ::StaticInt{Kc}, ::StaticInt{Nc}, M, K, N, atomicp::Ptr{UInt}, bc::Ptr, id::UInt, total_ids::UInt
) where {T, Mc, Kc, Nc}

    Niter = cld_fast(N, StaticInt{Nc}())
    Nblockdiv, Nrem = divrem_fast(N, Niter)
    Nblock_Nrem, Nblock_ = promote(Nblockdiv + One(), Nblockdiv)

    Kiter = cld_fast(K, StaticInt{Kc}())
    Kblockdiv, Krem = divrem_fast(K, Kiter)
    Kblock_Krem, Kblock_ = promote(Kblockdiv + One(), Kblockdiv)

    W = VectorizationBase.pick_vector_width_val(T)
    Miter = cld_fast(M, StaticInt{Mc}())
    Mbsize, Mrem, Mremfinal, _to_spawn = split_m(M, Miter, W) # M is guaranteed to be > W because of `W ≥ M` condition for `jmultsplitn!`...
    Mblock_Mrem, Mblock_ = promote(Mbsize + W, Mbsize)

    last_id = total_ids - one(UInt)
    atomics = atomicp + 8sizeof(UInt)
    sync_iters = zero(UInt64)

    Npackb_r_div, Npackb_r_rem = divrem_fast(Nblock_Nrem, total_ids)
    Npackb_r_block_rem, Npackb_r_block_ = promote(Npackb_r_div + One(), Npackb_r_div)

    Npackb___div, Npackb___rem = divrem_fast(Nblock_, total_ids)
    Npackb___block_rem, Npackb___block_ = promote(Npackb___div + One(), Npackb___div)

    pack_r_offset = Npackb_r_div * id + min(id, Npackb_r_rem)
    pack_r_view = CloseOpen(pack_r_offset, pack_r_offset + ifelse(id < Npackb_r_rem, Npackb_r_block_rem, Npackb_r_block_))
    pack___offset = Npackb___div * id + min(id, Npackb___rem)
    pack___view = CloseOpen(pack___offset, pack___offset + ifelse(id < Npackb___rem, Npackb___block_rem, Npackb___block_))

    GC.@preserve BCACHE begin
        for no in CloseOpen(Niter)
            # Krem
            # pack kc x nc block of B
            nfull = no < Nrem
            nsize = ifelse(nfull, Nblock_Nrem, Nblock_)
            pack_view = nfull ? pack_r_view : pack___view
            let A = A, B = B#, C = C
                for ko ∈ CloseOpen(Kiter)
                    ksize = ifelse(ko < Krem, Kblock_Krem, Kblock_)

                    Bsubset2 = PtrArray(B, (ksize, nsize), none_dense(Val{2}()))
                    Bpacked2 = ptrarray0(bc, (ksize, nsize))
                    copyto!(zview(Bpacked2, :, pack_view), zview(Bsubset2, :, pack_view))
                    # synchronize before starting the multiplication, to ensure `B` is packed
                    sync_iters += total_ids
                    _mv = _atomic_add!(atomicp, one(UInt))
                    while _mv < sync_iters
                        pause()
                        _mv = _atomic_max!(atomicp, zero(UInt))
                    end
                    # multiply
                    let A = A, B = zstridedpointer(Bpacked2), C = C
                        for mo in CloseOpen(Miter)
                            # @show mo, ko, no
                            msize = ifelse(mo == (Miter-1), Mremfinal, ifelse(mo < Mrem, Mblock_Mrem, Mblock_))
                            if ko == 0
                                packaloopmul!(C, A, B, α,     β, (msize, ksize, nsize))
                            else
                                packaloopmul!(C, A, B, α, One(), (msize, ksize, nsize))
                            end
                            A = gesp(A, (msize, Zero()))
                            C = gesp(C, (msize, Zero()))
                        end
                    end
                    A = gesp(A, (Zero(), ksize))
                    B = gesp(B, (ksize, Zero()))
                    # synchronize on completion of using `Bpacked`
                    _mv = _atomic_add!(atomics, one(UInt))
                    while _mv < sync_iters
                        pause()
                        _mv = _atomic_max!(atomics, zero(UInt))
                    end
                end
            end
            B = gesp(B, (Zero(), nsize))
            C = gesp(C, (Zero(), nsize))
        end
    end # GC.@preserve
    nothing
end

maybeinline(::Any, ::Any) = false
function maybeinline(C::AbstractStrideMatrix{<:Any,<:Any,T}, ::AbstractStrideMatrix{<:Any,<:Any,<:Any,1}) where {T}
    M, N = size(C)
    if M isa StaticInt && N isa StaticInt
        static_sizeof(T) * M * N < StaticInt{176}() * StaticInt(mᵣ) * StaticInt{nᵣ}()
    else
        false
    end
end
function maybeinline(C::AbstractStrideMatrix{<:Any,<:Any,T}, ::AbstractStrideMatrix) where {T}
    M, N = size(C)
    if M isa StaticInt && N isa StaticInt
        M * static_sizeof(T) ≤ StaticInt{2}() * StaticInt{VectorizationBase.REGISTER_SIZE}()
    else
        false
    end
end

@inline LinearAlgebra.mul!(C::AbstractStrideMatrix, A::StridedMatrix, B::StridedMatrix) = matmul!(C, A, B)
@inline LinearAlgebra.mul!(C::StridedMatrix, A::AbstractStrideMatrix, B::StridedMatrix) = matmul!(C, A, B)
@inline LinearAlgebra.mul!(C::StridedMatrix, A::StridedMatrix, B::AbstractStrideMatrix) = matmul!(C, A, B)
@inline LinearAlgebra.mul!(C::AbstractStrideMatrix, A::AbstractStrideMatrix, B::StridedMatrix) = matmul!(C, A, B)
@inline LinearAlgebra.mul!(C::AbstractStrideMatrix, A::StridedMatrix, B::AbstractStrideMatrix) = matmul!(C, A, B)
@inline LinearAlgebra.mul!(C::StridedMatrix, A::AbstractStrideMatrix, B::AbstractStrideMatrix) = matmul!(C, A, B)
@inline LinearAlgebra.mul!(C::AbstractStrideMatrix, A::AbstractStrideMatrix, B::AbstractStrideMatrix) = matmul!(C, A, B)

# @inline function Base.:*(
#     sp::StackPointer,
#     A::AbstractStrideMatrix{<:Any,<:Any,T},
#     B::AbstractStrideMatrix{<:Any,<:Any,T}
# ) where {T}
#     sp, D = PtrArray{T}(sp, (maybestaticsize(A, StaticInt{1}()),maybestaticsize(B, StaticInt{2}())))
#     sp, mul!(D, A, B)
# end
@inline function Base.:*(
    A::AbstractStrideMatrix{<:Any,<:Any,TA},
    B::AbstractStrideMatrix{<:Any,<:Any,TB}
) where {TA,TB}
    TC = promote_type(TA,TB)
    C = StrideArray{TC}(undef, (size(A, StaticInt{1}()),size(B, StaticInt{2}())))
    matmul!(C, A, B)
    C
end

@inline extract_λ(a) = a
@inline extract_λ(a::UniformScaling) = a.λ
@inline function Base.:*(A::AbstractStrideArray{S,D,T}, bλ::Union{Tb,UniformScaling{Tb}}) where {S,D,T<:VectorizationBase.NativeTypes,Tb <: Real}
    mv = similar(A)
    b = T(extract_λ(bλ))
    @avx for i ∈ eachindex(A)
        mv[i] = A[i] * b
    end
    mv
end

@inline function LinearAlgebra.mul!(
    C::AbstractStrideMatrix{<:Any,<:Any,T},
    A::LinearAlgebra.Diagonal{T,<:AbstractStrideVector{<:Any,T}},
    B::StridedMatrix{T}
) where {T}
    M, K, N = matmul_axes(C, A, B)
    MandK = ArrayInterface._pick_range(M, K)
    vA = parent(A)
    @avx for n ∈ N, m ∈ 1:MandK
        C[m,n] = vA[m] * B[m,n]
    end
    C
end
@inline function Base.:*(
    A::LinearAlgebra.Diagonal{T,<:AbstractVector{T}},
    B::AbstractStrideMatrix{<:Any,<:Any,T}
) where {T}
    mul!(similar(B), A, B)
end


text:
function getBarsFrame(leddata, audioSamp, rawspec)
    freqs = Array(domain(rawspec))
    samp = Array(audioSamp)
    spec = zeros(100)
    spec,maxspec = binFFT(rawspec, 10)
    for i in eachindex(spec)
        if isnan(spec[i])
            spec[i] = 0
        end
    end
    crossover = floor(Int, length(spec)/1.5)
    bottomEnd = mean((spec[1:crossover]))/maxspec
    topEnd = mean(spec[crossover+1:end])/maxspec
    for i in eachindex(leddata.channels)
        chan = leddata.channels[i]
        ref_len = length(leddata.channels[i])
        #print(i)
        if i % 4 == 1
            lows = !isnan(bottomEnd*ref_len) ? round(Int,bottomEnd*ref_len) : 0
            chan[1:lows] = colorant"blue"
            chan[lows+1:end] = colorant"black"
            #println(lows)
        elseif i % 4 == 2
            highs = !isnan(bottomEnd*ref_len) ? round(Int,topEnd*ref_len) : 0
            chan[1:highs] = colorant"red"
            chan[highs+1:end] = colorant"black"
            #println(highs)
        elseif i % 4 == 3
            lows = !isnan(bottomEnd*ref_len) ? round(Int,bottomEnd*ref_len) : 0
            chan[1:end-lows-1] = colorant"black"
            chan[end-lows:end] = colorant"blue"
            #println(lows)
        elseif i % 4 == 0
            highs = !isnan(bottomEnd*ref_len) ? round(Int,topEnd*ref_len) : 0
            chan[1:end-highs-1] = colorant"black"
            chan[end-highs:end] = colorant"red"
            #println(highs)
        else
            error("Weird modulo arithmetic failed.  Probably should've thrown an eror before this")
        end
        update!(chan)
    end
end
function rainbowBarsFrame(leddata, audioSamp, rawspec)
    colors = [convert(RGB{FixedPointNumbers.Normed{UInt8,8}}, HSL(x, 1, 0.5)) for x in linspace(0,360,length(leddata.channels)+1)[1:end-1]]
    return nBarsFrame(leddata, audioSamp, rawspec, colors)
end

function nBarsFrame(leddata, audioSamp, rawspec)
    colors = [convert(RGB{FixedPointNumbers.Normed{UInt8,8}}, HSL(240, x, 0.5)) for x in linspace(.25,1,length(leddata.channels))]
    return nBarsFrame(leddata, audioSamp, rawspec, colors)
end

function nBarsFrame(leddata, audioSamp, rawspec, colors)
    spec,maxspec = binFFT(rawspec[1:floor(Int, min(length(rawspec),length(rawspec)/48*length(leddata.channels)))], length(leddata.channels))
    default_color = colorant"black"
    normspec = zeros(length(spec))
    for i in eachindex(spec)
        ref_len = length(leddata.channels[i])
        tmp = spec[i]/fft_scale[]*ref_len
        if !isnan(tmp)
            normspec[i]=abs(tmp)
        end
    end
    min_cross = 0
    for i in eachindex(leddata.channels)
        chan = leddata.channels[i]
        ref_len = length(leddata.channels[i]);
        crossover = floor(Int, normspec[i])
        color_range = [i<=crossover for i in 1:ref_len]
        chan[1:crossover] = colors[i];
        chan[crossover+1:ref_len] = default_color;
        update!(chan)
    end
end

ext: jl


text:
using Plasmo
using Ipopt
using PlasmoTest

G = OptiGraph()
O = Ipopt.Optimizer

n1 = @optinode(G)
n2 = @optinode(G)
n3 = @optinode(G)
n4 = @optinode(G)

#Node 1 Model
@variable(n1,0 <= x1 <= 2)
@variable(n1,2 <= y1 <= 4)
@constraint(n1,x1 + y1 >= 5)
@objective(n1, Min, y1)

#Node 2 Model
@variable(n2,x2)
@NLnodeconstraint(n2,ref,exp(x2) >= 2)
@variable(n2,y2 >= 0)
@constraint(n2,x2 + y2 >= 3)
@objective(n2, Min, x2)

#Node 3 Model
@variable(n3,0 <= x3 <= 2)
@variable(n3,2 <= y3 <= 4)
@variable(n3,4 <= z3 <= 6)
@constraint(n3,9 <= x3 + y3 + z3 <= 11)

#Node 4 Model
@variable(n4,x4)
@NLnodeconstraint(n4,ref,x4^2 <= 16)
@variable(n4,y4 >= 0)
@variable(n4,z4 <= 4)
@constraint(n4,(x4 + y4)*z4 >= 12)

#Add linkconstraint
@linkconstraint(G,n1[:x1]==n2[:y2])
@linkconstraint(G,n2[:x2]==n3[:z3])
@linkconstraint(G,n2[:y2]==n4[:y4])

#Graph objective
@objective(G, Min, n1[:y1] + n2[:x2] + n3[:z3] + n4[:x4])

optimize!(G,O)

println("x1 = ",digit5(value(n1,n1[:x1])))
println("y1 = ",digit5(value(n1,n1[:y1])))
println("x2 = ",digit5(value(n2,n2[:x2])))
println("y2 = ",digit5(value(n2,n2[:y2])))
println("x3 = ",digit5(value(n3,n3[:x3])))
println("y3 = ",digit5(value(n3,n3[:y3])))
println("z3 = ",digit5(value(n3,n3[:z3])))
println("x4 = ",digit5(value(n4,n4[:x4])))
println("y4 = ",digit5(value(n4,n4[:y4])))
println("z4 = ",digit5(value(n4,n4[:z4])))
println("Objective = ", digit5(objective_value(G)))

#plt_G = plt_G = Plots.plot(G,node_labels = true, markersize=60, labelsize=30, linewidth=4, layout_options = Dict(:tol => 0.01, :iterations => 2))
#savefig("002.png")

max_stars_repo_name: BoyangSheng/Skill-studing

max_issues_repo_name: BoyangSheng/Skill-studying

max_forks_repo_name: BoyangSheng/Skill-studying


text:
using SpatialFields
using Base.Test
import MultiPoly

function test_linear_fit()
    x, y = MultiPoly.generators(MultiPoly.MPoly{Float64}, :x, :y)

    true_polynomial = 3 + 4x - 7y

    points = SVector{2, Float64}[[5; -2], [-1; 10], [-1; 11]]
    data = Float64[MultiPoly.evaluate(true_polynomial, p...) for p in points]

    fit_polynomial = SpatialFields.linear_fit(points, data)

    for p in points
        @test isapprox(MultiPoly.evaluate(fit_polynomial, p...), MultiPoly.evaluate(true_polynomial, p...))
    end
    @test isapprox(fit_polynomial[0,0], true_polynomial[0,0])
    @test isapprox(fit_polynomial[1,0], true_polynomial[1,0])
    @test isapprox(fit_polynomial[0,1], true_polynomial[0,1])
end

test_linear_fit()

function test_vector_linear_fit()
    x, y = MultiPoly.generators(MultiPoly.MPoly{Float64}, :x, :y)

    true_polynomials = [3 + 4x - 7y; -1 - 8x + 2y]
    points = SVector{2, Float64}[[5; -2], [-1;10], [-1; 11]]
    data = SVector{2, Float64}[[MultiPoly.evaluate(poly, p...) for poly in true_polynomials] for p in points]

    fit_polynomials = SpatialFields.linear_fit(points, data)

    for i = 1:length(true_polynomials)
        for p in points
            @test isapprox(MultiPoly.evaluate(fit_polynomials[i], p...), MultiPoly.evaluate(true_polynomials[i], p...))
        end
        @test isapprox(fit_polynomials[i][0,0], true_polynomials[i][0,0])
        @test isapprox(fit_polynomials[i][1,0], true_polynomials[i][1,0])
        @test isapprox(fit_polynomials[i][0,1], true_polynomials[i][0,1])
    end
end

test_vector_linear_fit()


text:
module KhepriTikZ
using KhepriBase

# functions that need specialization
include(khepribase_interface_file())

include("TikZ.jl")

function __init__()
  add_current_backend(tikz)
end

end


text:
const NOTEBOOKS = [("Week 1: Introduction to Julia", "introduction.jl"),
                   ("Week 2: Supervised Learning", "supervised_learning.jl"),
#                    ("Week 3: Generalized Linear Regression", "generalized_linear_regression.jl"),
#                    ("Week 4: Flexibility and Bias-Variance-Decomposition", "flexibility.jl"),
#                    ("Week 4: Model Assessment", "model_evaluation.jl"),
#                    ("Week 5: Regularization", "regularization.jl"),
#                    ("Week 5: Transformations of Input or Output", "transformations.jl"),
#                    ("Week 6: Gradient Descent", "gradient_descent.jl"),
#                    ("Week 7: Multilayer Perceptrons", "mlp.jl"),
#                    ("Week 8: Convolutional Neural Networks", "convnets.jl"),
#                    ("Week 8: Recurrent Neural Networks", "rnn.jl"),
#                    ("Week 8: Tree-Based Methods", "trees.jl"),
#                    ("Week 9: Clustering", "clustering.jl"),
#                    ("Week 9 & 10: Principal Component Analysis", "pca.jl"),
#                    ("Week 11 & 12: Reinforcement Learning", "rl.jl")
                  ]

function _linkname(path, nb, basedir)
    if haskey(ENV, "html_export") && ENV["html_export"] == "true"
        joinpath(basedir, "$(splitext(nb)[1]).html")
    else
        "open?path=" * joinpath(path, nb)
    end
end
function list_notebooks(file, basedir = "")
    sp = splitpath(file)
    path = joinpath(@__DIR__, "..", "notebooks")
    filename = split(sp[end], "#")[1]
    list = join(["1. " * (nb == filename ?
                            name * " (this notebook)" :
                            "[$name](" * _linkname(path, nb, basedir) * ")")
                 for (name, nb) in NOTEBOOKS], "\n")
Markdown.parse("""# Course Overview

               $list
               """)
end

function footer()
    html"""

 This page is part of an introductory machine learning course taught by Johanni Brea.
The course is strongly inspired by "An Introduction to Statistical Learning".

    """
end

max_stars_repo_name: ablaom/MLCourse

max_issues_repo_name: ablaom/MLCourse

max_forks_repo_name: ablaom/MLCourse


text:
function step_radiation!(
    mb :: ModelBlock,
    Δt :: Float64,
)
    fi = mb.fi
    co = mb.co

    fi._X_[:, 1] .+= Δt * ( co.mtx[:T_swflxConv_sT] * view(fi.SWFLX, :) + co.mtx[:T_nswflxConv_sT] * view(fi.NSWFLX, :)) / ρc_sw

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using HDF5, ProgressMeter, Interpolations
# pyplot()
gr()
default(show=false)

const nscale = 1
# const sz2 = round(Int, sz/nscale)
const sz3 = round(Int, sz/2)

function get_parameters(track::Track, Δx::Float64)
    x = last.(track.coordinates)
    y = first.(track.coordinates)
    n = length(track.lengths)
    m = minimum(minimum(i) for i in track.intensities)
    intensities = fill(m, n, n)
    for i in 1:n, j in 1:i
        intensities[j,i] = track.intensities[i][j]
    end
    lengths = track.lengths*Δx
    times = track.times
    return (x, y, n, intensities, lengths, times, m)
end

function saveit(calibstages::Vector{CalibStage}, tss::Vector{Vector{Track}}, output::IO)
    pm = Progress(sum(length(st.timelapse)*length(ts) for (st, ts) in zip(calibstages, tss)), desc = "Saving the results", output = output)
    path = joinpath(calibstages[1].home, calibstages[1].base)
    for (st, ts) in zip(calibstages, tss)
        saveit(st, ts, path, pm)
    end
end

function filterdone!(rs::Vector{Track}, path::String)
    # isdir(path) || return rs
    fs = readdir(path)
    ffs = joinpath.(path, fs)
    summary = any(zip(ffs, fs)) do fff
        ff, f = fff
        isfile(ff) && f == "roots.png"
    end
    summary || return rs
    filter!(rs) do r
        !any(zip(ffs, fs)) do fff
            ff, f = fff
            isdir(ff) && f == "root $(r.id)" && info("Results for root #$(r.id) already exist (delete folder $ff if you want to recalculate)…") == nothing
        end
    end
end

function saveit(st::CalibStage, rs::Vector{Track}, path::String, pm::Progress)
    path = joinpath(path, "stage $(st.id)")
    mkpath(path)
    info("Stage $(st.id)")
    filterdone!(rs, path)
    isempty(rs) && return nothing
    imgs = [Float64.(load(tl.path)) for tl in st.timelapse]
    # imgs = [Float64.(imresize(load(tl.path), (sz2, sz2))) for tl in st.timelapse]
    # mM = mean(quantile(vec(img), [.1, .995]) for img in imgs)
    # for img in imgs
        # img .= imadjustintensity(img, mM)
    # end
    formatlabel(x::T) where T <: Real = round(st.Δx*x*nscale, 1)
    saveoverview(deepcopy(imgs[end]), formatlabel, rs, path)
    for r in rs
        nodes = detect_nodes(r)
        # println(nodes)
        x, y, n, intensities, lengths, times, Imin = get_parameters(r, st.Δx)
        _path = joinpath(path, "root $(r.id)")
        mkpath(_path)
        save2csv(_path, x, y, intensities, lengths, times, nodes)
        save2gif(_path, nodes, x/nscale, y/nscale, n, intensities, lengths, times, formatlabel, imgs, st.Δt, Imin, pm)
    end
end

function saveoverview(img::Matrix{Float64}, formatlabel::Function, rs::Vector{Track}, path::String)
    img .= adjust_gamma(img, 2.4)
    mM = quantile(vec(img), [.5, .995])
    img .= imadjustintensity(img, mM)
    heatmap(flipdim(img, 1), aspect_ratio = 1, yformatter = formatlabel, xformatter = formatlabel, xlabel = "X (mm)", ylabel = "Y (mm)", color=:inferno, yflip=true, colorbar=false, legend=false)#, size=(sz,sz), dpi=50/3)
    for r in rs
        x = last.(r.coordinates)
        y = first.(r.coordinates)
        plot!(x, y, color=r.color, linestyle = :dot, linewidth = 5, linewidth = intensity_radius, annotations=(x[1],y[1],text(string(r.id), :left, r.color, 25)))
    end
    png(joinpath(path, "roots"))
end

function save2csv(path::String, x, y, intensities, lengths, times, nodes)
    z = [[nothing; lengths] [RowVector(times); intensities]]
    writecsv(joinpath(path, "intensities.csv"), z)
    writecsv(joinpath(path, "coordinates.csv"), zip(x, y))
    writecsv(joinpath(path, "nodes.csv"), [[lengths[node.length_i], times[node.time_i], node.Δ, node.p] for node in nodes])
end

#=function save2gif(path::String, x, y, n::Int, intensities, lengths, times, imgs, Δt::Float64, pm::Progress)
    _imgs = deepcopy(imgs)
    p = [(round(Int, xi), round(Int, yi)) for (xi, yi) in zip(x, y)]
    for i in 1:n
        draw!(_imgs[i], Path(p[1:i]), RGB{N0f16}(1,0,0))
        next!(pm)
    end
    save(joinpath(path, "root.gif"), cat(3, _imgs...), fps = round(Int, 1/(5/180000*60*60*Δt)))
    Imax = quantile(vec(intensities), 0.98)
    intensities .= min.(intensities, Imax)
    heatmap(times, lengths, intensities, xlabel = "Time (hrs)", ylabel = "Root length (mm)", yflip = true, colorbar = false)
    png(joinpath(path, "intensities.png"))
end=#

function save2gif(path::String, nodes, x, y, n::Int, intensities, lengths, times, formatlabel::Function, imgs, Δt::Float64, Imin::Float64, pm::Progress)
    Imax = quantile(vec(intensities), 0.95)
    intensities .= min.(intensities, Imax)
    xlim = round.(Int, extrema(x)) .+ [-10, 10]
    clamp!(xlim, 1, sz)
    ylim = round.(Int, extrema(y)) .+ [-10, 10]
    clamp!(ylim, 1, sz)
    imgs2 = [adjust_gamma(img[ylim[1]:ylim[2], xlim[1]:xlim[2]], 2.4) for img in imgs]
    mM = mean(quantile(vec(img), [.5, .99]) for img in imgs2)
    for img in imgs2
        img .= imadjustintensity(img, mM)
    end
    itp = interpolate((times, lengths), log.(intensities), Gridded(Linear()))
    xl = linspace(extrema(times)..., n)
    yl = linspace(extrema(lengths)..., n)
    zl = itp[xl, yl]
    Ibounds = extrema(zl)
    anim = Animation()
    for i in 1:n
        # h1 = heatmap(flipdim(imgs2[i],1), aspect_ratio = 1, yformatter = formatlabel, xlabel = "X (mm)", ylabel = "Y (mm)", color=:inferno, yflip=false, colorbar=false, xticks=[])
        h1 = heatmap(flipdim(imgs2[i], 1), aspect_ratio = 1, yformatter = formatlabel, xlabel = "X (mm)", ylabel = "Y (mm)", color=:inferno, yflip=true, colorbar=false, xticks=[])
        # h1 = heatmap(imadjustintensity(adjust_gamma(imgs[i], 2.4), mM), aspect_ratio = 1, color=:inferno, yflip=false, colorbar=false)
        foreach(nodes) do node
            Plots.annotate!(h1, x[node.length_i] - xlim[1], y[node.length_i] - ylim[1], text(string(round(node.Δ, 2), "-", round(node.p, 2)), RGBA(0,1,0, Int(i > node.time_i)), 8))
            # Plots.annotate!(h1, x[node.length_i] - xlim[1], y[node.length_i] - ylim[1], node.label(i))
            # Plots.scatter!(h1, [x[node.length_i] - xlim[1]], [y[node.length_i] - ylim[1]])
            # Plots.annotate!(h1, x[node.length_i], y[node.length_i], node.label(i))
            # Plots.scatter!(h1, [x[node.length_i]], [y[node.length_i]])
        end
        # plot!(x[1:i]-xlim[1], y[1:i]-ylim[1], color = :white)
        h3 = plot(xl, zl[i, :], fill = Ibounds[1], fillcolor = :blue, linecolor = :blue, ylim = Ibounds, xlim=extrema(xl), yticks = nothing, ylabel = "Intensity")
        h2 = plot([Ibounds[1]; zl[:,i]; Ibounds[1]], [yl[1]; yl; yl[end]], fill = Ibounds[1], fillcolor = :green, linecolor = :green, xlim = Ibounds, ylim=extrema(yl), xticks = nothing,  yflip = true, xlabel = "Intensity")
        h4 = heatmap(xl, yl, flipdim(zl,1), xlabel = "Time (hrs)", ylabel = "Root length (mm)", yflip = true, colorbar = false, color=:inferno)
        plot!(h4, xl[[1, end]], [yl[i], yl[i]], color = :blue)
        plot!(h4, [xl[i], xl[i]], yl[[1, end]], color = :green)
        foreach(nodes) do node
            # Plots.annotate!(h4, times[node.time_i], lengths[node.length_i], node.label(i))
        end
        plot(h3, h1, h4, h2, legend = false, size=(sz3, sz3), dpi=50/3)#, background_color=:black)
        Plots.frame(anim)
        next!(pm)
    end
    # mp4(anim, joinpath(path, "summary.mp4"), fps = round(Int, 24/10Δt))
    # for when the latest Plots get tagged:
    mp4(anim, joinpath(path, "summary.mp4"), fps = round(Int, 24/10Δt), show_msg=false)
end


text:
function tournament_selection(P, a = rand(1:length(P)))
    # chose two different solutions at random
    b = rand(1:length(P))
    while a == b
        b = rand(1:length(P))
    end

    # perform selection
    P[a].rank < P[b].rank || (P[a].rank == P[b].rank && P[a].crowding > P[b].crowding ) ? P[a] : P[b]
end

function gen_β(β, η, D, R)
    α = 2.0 .- β .^ (-  η - 1.0 )
    mask = R .<= 1.0 ./ α
    s = 1.0 / (η + 1.0)
    βq = [ mask[i] ?  (R[i] * α[i])^s : (1.0 / (2.0 - R[i]*α[i]))^s for i in 1:D]
    βq
end

function SBX_crossover(vector1, vector2, bounds, η=15, p_variable = 0.9)
    xu = view(bounds, 2,:)
    xl = view(bounds, 1,:)
    D = length(vector1)

    do_crossover = ones(Bool, D)
    do_crossover[rand(D) .> p_variable] .= false
    do_crossover[ abs.( vector2 - vector1 ) .<= eps() ] .= false

    y1 = min.( vector1, vector2 )
    y2 = max.( vector1, vector2 )
    Δ = max.(eps(), y2 - y1)

    R = rand(D)

    β = @. 1.0 + (2.0 * (y1 - xl) / Δ)
    βq = gen_β(β, η, D, R)
    c1 = @. 0.5*(y1 + y2 -  βq*Δ)

    β = @. 1.0 + (2.0 * (y1 - xl) / Δ)
    βq = gen_β(β, η, D, R)
    c2 = @. 0.5*(y1 + y2 +  βq*Δ)

    # swap
    mask = rand(Bool, D)
    cc = copy(c1)
    c1[mask] = c2[mask]
    c2[mask] = cc[mask]

    cc1 = copy(vector1)
    cc1[do_crossover] = c1[do_crossover]
    cc2 = copy(vector2)
    cc2[do_crossover] = c2[do_crossover]

    reset_to_violated_bounds!(cc1, bounds)
    reset_to_violated_bounds!(cc2, bounds)

    return cc1, cc2
end

function polynomial_mutation!(vector, bounds, η=20, prob = 1 / length(vector))
    do_mutation = rand(length(vector)) .< prob

    xu = view(bounds, 2,do_mutation)
    xl = view(bounds, 1,do_mutation)
    x = view(vector, do_mutation)

    δ1 = (x - xl) ./ (xu - xl)
    δ2 = (xu - x) ./ (xu - xl)

    D = length(xu)
    R = rand(D)
    mask = R .< 0.5
    s = η+1.0
    mut_pow = 1.0 / (η + 1.0)
    δq = [ mask[i] ?
            ^(2.0R[i] + (1. - 2.0R[i]) * ^(1.0 - δ1[i], s), mut_pow) - 1.0 :
            1.0 - (2.0 * (1.0 - R[i]) + 2.0 * (R[i] - 0.5) * ^(1.0 - δ2[i], s))^mut_pow
            for i in 1:D
        ]

    vector[do_mutation] = x + δq .* ( xu - xl)
    # correct using reset to bound
    #
    vector

end


text:
module LinearOptimizationSolver

using LinearAlgebra
using MKL

include("types.jl")
include("utils.jl")
include("log.jl")
include("interior_points.jl")
include("simplex.jl")

export create, solve, crossover

end # module


text:
@testset "base" begin
    @testset "Trig" begin
        @testset "Basics" for x = (Float64(π)-0.01, Complex(π, π/2))
            test_scalar(sin, x)
            test_scalar(cos, x)
            test_scalar(tan, x)
            test_scalar(sec, x)
            test_scalar(csc, x)
            test_scalar(cot, x)
            test_scalar(sinpi, x)
            test_scalar(cospi, x)
        end
        @testset "Hyperbolic" for x = (Float64(π)-0.01, Complex(π-0.01, π/2))
            test_scalar(sinh, x)
            test_scalar(cosh, x)
            test_scalar(tanh, x)
            test_scalar(sech, x)
            test_scalar(csch, x)
            test_scalar(coth, x)
        end
        @testset "Degrees" begin
            x = 45.0
            test_scalar(sind, x)
            test_scalar(cosd, x)
            test_scalar(tand, x)
            test_scalar(secd, x)
            test_scalar(cscd, x)
            test_scalar(cotd, x)
        end
        @testset "Inverses" for x = (0.5, Complex(0.5, 0.25))
            test_scalar(asin, x)
            test_scalar(acos, x)
            test_scalar(atan, x)
            test_scalar(asec, 1/x)
            test_scalar(acsc, 1/x)
            test_scalar(acot, 1/x)
        end
        @testset "Inverse hyperbolic" for x = (0.5, Complex(0.5, 0.25))
            test_scalar(asinh, x)
            test_scalar(acosh, x + 1)  # +1 accounts for domain
            test_scalar(atanh, x)
            test_scalar(asech, x)
            test_scalar(acsch, x)
            test_scalar(acoth, x + 1)
        end
        @testset "Inverse degrees" for x = (0.5, Complex(0.5, 0.25))
            test_scalar(asind, x)
            test_scalar(acosd, x)
            test_scalar(atand, x)
            test_scalar(asecd, 1/x)
            test_scalar(acscd, 1/x)
            test_scalar(acotd, 1/x)
        end
        @testset "Multivariate" begin
            @testset "atan2" begin
                # https://en.wikipedia.org/wiki/Atan2
                x, y = rand(2)
                ratan = atan(x, y)
                u = x^2 + y^2
                datan = y/u - 2x/u

                r, pushforward = frule(atan, x, y)
                @test r === ratan
                @test pushforward(NamedTuple(), 1, 2) === datan

                r, pullback = rrule(atan, x, y)
                @test r === ratan
                dself, df1, df2 = pullback(1)
                @test dself == NO_FIELDS
                @test df1 + 2df2 === datan
            end

            @testset "sincos" begin
                x, Δx, x̄ = randn(3)
                Δz = (randn(), randn())

                frule_test(sincos, (x, Δx))
                rrule_test(sincos, Δz, (x, x̄))
            end
        end
    end  # Trig

    @testset "math" begin
        for x in (-0.1, 6.4, 1.0+0.5im, -10.0+0im)
            test_scalar(deg2rad, x)
            test_scalar(rad2deg, x)

            test_scalar(inv, x)

            test_scalar(exp, x)
            test_scalar(exp2, x)
            test_scalar(exp10, x)

            x isa Real && test_scalar(cbrt, x)
            if (x isa Real && x >= 0) || x isa Complex
                # this check is needed because these have discontinuities between
                # `-10 + im*eps()` and `-10 - im*eps()`
                should_test_wirtinger = imag(x) != 0 && real(x) < 0
                test_scalar(sqrt, x; test_wirtinger=should_test_wirtinger)
                test_scalar(log, x; test_wirtinger=should_test_wirtinger)
                test_scalar(log2, x; test_wirtinger=should_test_wirtinger)
                test_scalar(log10, x; test_wirtinger=should_test_wirtinger)
                test_scalar(log1p, x; test_wirtinger=should_test_wirtinger)
            end
        end
    end

    @testset "Unary complex functions" begin
        for x in (-4.1, 6.4, 1.0+0.5im, -10.0+1.5im)
            test_scalar(real, x)
            test_scalar(imag, x)

            test_scalar(abs, x)
            test_scalar(hypot, x)

            test_scalar(angle, x)
            test_scalar(abs2, x)
            test_scalar(conj, x)
        end
    end

    @testset "*(x, y)" begin
        x, y = rand(3, 2), rand(2, 5)
        z, pullback = rrule(*, x, y)

        @test z == x * y

        z̄ = rand(3, 5)
        (ds, dx, dy) = pullback(z̄)

        @test ds === NO_FIELDS

        @test extern(dx) == extern(accumulate(zeros(3, 2), dx))
        @test extern(dy) == extern(accumulate(zeros(2, 5), dy))

        test_accumulation(rand(3, 2), dx)
        test_accumulation(rand(2, 5), dy)
    end

    @testset "binary trig ($f)" for f in (hypot, atan)
        rng = MersenneTwister(123456)
        x, Δx, x̄ = 10randn(rng, 3)
        y, Δy, ȳ = randn(rng, 3)
        Δz = randn(rng)

        frule_test(f, (x, Δx), (y, Δy))
        rrule_test(f, Δz, (x, x̄), (y, ȳ))
    end

    @testset "identity" begin
        rng = MersenneTwister(1)
        rrule_test(identity, randn(rng), (randn(rng), randn(rng)))
        rrule_test(identity, randn(rng, 4), (randn(rng, 4), randn(rng, 4)))
    end

    @testset "Constants" for x in (-0.1, 6.4, 1.0+0.5im, -10.0+0im, 0+200im)
        test_scalar(one, x)
        test_scalar(zero, x)
    end

    @testset "sign" begin
        @testset "at points" for x in (-1.1, -1.1, 0.5, 100)
            test_scalar(sign, x)
        end

        @testset "Zero over the point discontinuity" begin
            # Can't do finite differencing because we are lying
            # following the subgradient convention.

            _, pb = rrule(sign, 0.0)
            _, x̄ = pb(10.5)
            @test extern(x̄) == 0

            _, pf = frule(sign, 0.0)
            ẏ = pf(NamedTuple(), 10.5)
            @test extern(ẏ) == 0
        end
    end
end


text:
"""
    wh_learn(X, Y)

Widrow-Hoff Learning.

# Obligatory Arguments
- `test_mode::Symbol`: which test mode, currently supports :train_only, :pre_split, :careful_split and :random_split.

# Optional Arguments
- `eta::Float64=0.1`: the learning rate
- `n_epochs::Int64=1`: the number of epochs to be trained
- `weights::Matrix=nothing`: the initial weights
- `learn_seq::Vector=nothing`: the learning sequence
- `save_history::Bool=false`: if true, a partical training history will be saved
- `history_cols::Vector=nothing`: the list of column indices you want to saved in history, e.g. `[1,32,42]` or `[2]`
- `history_rows::Vector=nothing`: the list of row indices you want to saved in history, e.g. `[1,32,42]` or `[2]`
- `verbose::Bool = false`: if true, more information will be printed out
"""
function wh_learn(
    X,
    Y;
    eta = 0.01,
    n_epochs = 1,
    weights = nothing,
    learn_seq = nothing,
    save_history = false,
    history_cols = nothing,
    history_rows = nothing,
    verbose = false,
    )

    X = Array(X)
    Y = Array(Y)

    if size(X, 1) != size(Y, 1)
        throw(ArgumentError("X($(size(inputs,1))) and Y($(size(outputs,1))) length doesn't match"))
    end

    if isnothing(weights)
        W = zeros(Float64, (size(X, 2), size(Y, 2)))
    else
        W = weights
    end

    # construct learn_seq if nothing
    if isnothing(learn_seq)
        learn_seq = 1:size(X, 1)
    end

    if save_history
        history = zeros(Float64, n_epochs, length(history_rows), length(history_cols))
    end

    inputT = Matrix{Float64}(undef, (size(X, 2), 1))
    pred = Matrix{Float64}(undef, (1, size(Y, 2)))
    deltaW = Matrix{Float64}(undef, (size(X, 2), size(Y, 2)))
    if verbose
        if isnothing(learn_seq)
            pb = Progress(size(X, 1) * n_epochs)
        else
            pb = Progress(size(learn_seq, 1) * n_epochs)
        end
    end
    for j = 1:n_epochs # 100 epochs
        for i in learn_seq # for each events
            # pred = X[i:i, :]*W
            mul!(pred, view(X, i:i, :), W)
            # obsv = Y[i:i, :]-pred
            broadcast!(-, pred, view(Y, i:i, :), pred)
            # inputT = X[i:i, :]'
            transpose!(inputT, view(X, i:i, :))
            # update = inputT*obsv
            mul!(deltaW, inputT, pred)
            # deltaW = eta*update
            rmul!(deltaW, eta)
            # W += deltaW
            broadcast!(+, W, W, deltaW)

            verbose && ProgressMeter.next!(pb)
        end
        # push history
        if save_history
            history[j,:,:] = copy(W[history_rows, history_cols])
        end
    end
    if save_history
        return W, history
    end
    W
end

"""
    make_learn_seq(freq)

Make Widrow-Hoff learning sequence.
"""
function make_learn_seq(freq; random_seed = 314)
    if isnothing(freq)
        return nothing
    end

    learn_seq = [repeat([i], n) for (i, n) in enumerate(freq)]
    learn_seq = collect(Iterators.flatten(learn_seq))
    rng = MersenneTwister(random_seed)
    shuffle(rng, learn_seq)
end


text:
module TestFunctions

include("preamble.jl")

function makedata(V::Type, M::Integer, N::Integer)
    dims = testdims(M + N)
    inner_size, outersize = tuple_split(dims, Val(M))

    nested = Array{Array{V,M},N}(undef, outersize...)
    for I in eachindex(nested)
        nested[I] = rand!(Array{V,M}(undef, inner_size...))
    end

    return (
        nested = nested,
        flat = reshape(reduce(hcat, nested), dims),
        dims = dims,
        inner_size = inner_size,
        outersize = outersize,
        inner_axes = axes(first(nested)),
        outeraxes = axes(nested),
    )
end

@testset "M = $M, N = $N, V = $V" for M = 0:2, N = 0:2, V in (Float64,)
    @testset "inner_*" begin
        data = makedata(V, M, N)
        A = data.nested

        @test inner_eltype(A) === V
        @test_inferred inner_eltype(A)
        @test inner_eltype(typeof(A)) === V
        @test_inferred inner_eltype(typeof(A))

        @test inner_ndims(A) == M
        @test_inferred inner_ndims(A)
        @test inner_ndims(typeof(A)) == M
        @test_inferred inner_ndims(typeof(A))

        @test inner_axes(A) == data.inner_axes
        @test_inferred inner_axes(A)

        @test inner_size(A) == data.inner_size
        @test_inferred inner_size(A)

        @test inner_length(A) == prod(data.inner_size)
        @test_inferred inner_length(A)
    end
end

@testset "inner_* errors" begin
    A = [zeros(2), zeros(3)]
    @test_throws ArgumentError inner_axes(A)
    @test_throws ArgumentError inner_size(A)
    @test_throws ArgumentError inner_length(A)
end

end # module


text:
using Pkg
Pkg.activate(joinpath("test_env"))
Pkg.add("GitHub")
Pkg.instantiate()

using GitHub

ORG, REPO, PR = ENV["org"], ENV["repo"], ENV["pullrequest"]
TEST_RESULTS_FILE = "$(ORG)_$(REPO)_$(PR).txt"

# Need to add GITHUB_AUTH to your .bashrc
myauth = GitHub.authenticate(ENV["GITHUB_AUTH"])

function create_gist(authentication)

    file_content = ""
    file = open(TEST_RESULTS_FILE, "r")
    for line in readlines(file)
        file_content *= line*'\n'
    end
    close(file)

    file_dict = Dict(TEST_RESULTS_FILE => Dict("content" => file_content))
    gist = Dict{String,Any}("description" => "Test results",
                             "public" => true,
                             "files" => file_dict)

    posted_gist = GitHub.create_gist(params = gist, auth = authentication)

    return posted_gist
end

function post_gist_url_to_pr(comment::String; kwargs...)
    api = GitHub.DEFAULT_API
    repo = get_repo(api, ORG, REPO; kwargs...)
    pull_request = get_pull_request(api, ORG, repo, parse(Int, PR); kwargs...)
    GitHub.create_comment(api, repo, pull_request, comment; kwargs...)
end

function get_repo(api::GitHub.GitHubWebAPI, org::String, repo_name::String; kwargs...)
    my_params = Dict(:visibility => "all")
    # return GitHub.repo(api, repo; params = my_params, kwargs...)
    return Repo(GitHub.gh_get_json(api, "/repos/$(org)/$(repo_name)"; params = my_params, kwargs...))
end

function get_pull_request(api::GitHub.GitHubWebAPI, org::String, repo::Repo, pullrequest_id; kwargs...)
    my_params = Dict(:sort => "popularity",
                    :direction => "desc")
    pull_request = PullRequest(GitHub.gh_get_json(api, "/repos/$(org)/$(repo.name)/pulls/$(pullrequest_id)"; params=my_params, kwargs...))

    return pull_request
end

function get_comment_from_test_results()
    open(TEST_RESULTS_FILE, "r") do file
        text_to_match = r"tests passed"
        for line in readlines(file)
            if occursin(text_to_match, line)
                return "$(strip(line)): "
            end
        end
        return "Tests failed: "
    end
end

comment = get_comment_from_test_results()
post_gist_url_to_pr("$comment $(create_gist(myauth).html_url)"; auth = myauth)


text:
#=

Here is the process for updating Khepri

=#

# The Khepri Julia project itself is stored at

julia_khepri = joinpath(dirname(dirname(abspath(@__FILE__))), "Khepri")

upgrade_autocad(; phase="Debug") =
    let # 1. The dlls are updated in VisualStudio after compilation of the plugin, and they are stored in the folder.
        dlls = ["KhepriBase.dll", "KhepriAutoCAD.dll"]
        # 2. Depending on whether we are in Debug mode or Release mode,
        development_phase = phase # "Release"
        # 3. the dlls are located in a folder
        dlls_folder = joinpath("bin", "x64", development_phase)
        # 4. contained inside the Plugins folder, which has a specific location regarding this file itself
        plugin_folder = joinpath(dirname(dirname(dirname(abspath(@__FILE__)))), "Plugins", "KhepriAutoCAD", "KhepriAutoCAD")
        # 5. Besides the dlls, we also need the bundle folder
        bundle_name = "Khepri.bundle"
        bundle_dll_folder = joinpath(bundle_name, "Contents")
        bundle_xml = joinpath(bundle_name, "PackageContents.xml")
        # 6. which is contained in the Plugins folder
        bundle_path = joinpath(plugin_folder, bundle_name)
        # 7. The bundle needs to be copied to the current folder
        local_bundle_path = joinpath(julia_khepri, "Plugins", "AutoCAD", bundle_name)
        # 8. but, before, we remove any previously existing bundle
        mkpath(dirname(local_bundle_path))
        if isdir(local_bundle_path)
            rm(local_bundle_path, force=true, recursive=true)
        end
        # 9. Now we do the copy
        cp(bundle_path, local_bundle_path)
        # 10. and we copy the dlls to the local bundle Contents folder
        local_bundle_contents_path = joinpath(local_bundle_path, "Contents")
        for dll in dlls
            src = joinpath(plugin_folder, dlls_folder, dll)
            dst = joinpath(local_bundle_contents_path, dll)
            rm(dst, force=true)
            cp(src, dst)
        end
    end

export upgrade_autocad

#=
After updating something in the plugin, just call this function:

upgrade_autocad()

=#

#=
Pkg3 is not yet ready. We need to manually generate REQUIRE
=#

import Pkg
const PT = Pkg.Types

generate_require(dir=julia_khepri) =
    begin
        Pkg.activate(dir)             # current directory as the project
        ctx = PT.Context()
        pkg = ctx.env.pkg
        if pkg ≡ nothing
            @error "Not in a package, I won't generate REQUIRE."
        else
            @info "found package" pkg = pkg
        end
        deps = PT.get_deps(ctx)
        non_std_deps = sort(collect(setdiff(keys(deps), values(ctx.stdlibs))))
        open(joinpath(julia_khepri, "REQUIRE"), "w") do io
            for d in non_std_deps
                println(io, d)
                @info "listing $d"
            end
        end
    end

#=
generate_require()
=#

export generate_require

max_stars_repo_name: aptmcl/Khepri

max_issues_repo_name: aptmcl/Khepri

max_forks_repo_name: aptmcl/Khepri


text:
@testset "names" begin
    model = lp1_test()
    MOI.set(model, MOI.VariableName(), VI(1), "x1")
    MOI.set(model, MOI.VariableName(), VI(2), "x2")
    @test model.var_to_name[VI(1)] == "x1"
    @test model.var_to_name[VI(2)] == "x2"
    MOI.set(
        model,
        MOI.ConstraintName(),
        CI{SAF{Float64},MOI.LessThan{Float64}}(1),
        "lessthan",
    )
    @test model.con_to_name[CI{SAF{Float64},MOI.LessThan{Float64}}(1)] ==
          "lessthan"

    # Dualize without names
    dual_problem = dualize(model)
    dual_model = dual_problem.dual_model
    primal_dual_map = dual_problem.primal_dual_map
    # Query variable names
    vi_1 =
        primal_dual_map.primal_con_dual_var[CI{VI,MOI.GreaterThan{Float64}}(1)][1]
    vi_2 = primal_dual_map.primal_con_dual_var[CI{
        SAF{Float64},
        MOI.LessThan{Float64},
    }(
        1,
    )][1]
    @test MOI.get(dual_model, MOI.VariableName(), vi_1) == ""
    @test MOI.get(dual_model, MOI.VariableName(), vi_2) == ""
    # Query constraint names
    ci_1 = primal_dual_map.primal_var_dual_con[VI(1)]
    ci_2 = primal_dual_map.primal_var_dual_con[VI(2)]
    @test MOI.get(dual_model, MOI.ConstraintName(), ci_1) == ""
    @test MOI.get(dual_model, MOI.ConstraintName(), ci_2) == ""

    dual_problem = dualize(
        model;
        dual_names = Dualization.DualNames("dualvar_", "dualcon_"),
    )
    dual_model = dual_problem.dual_model
    primal_dual_map = dual_problem.primal_dual_map
    # Query variable names
    vi_1 =
        primal_dual_map.primal_con_dual_var[CI{VI,MOI.GreaterThan{Float64}}(1)][1]
    vi_2 = primal_dual_map.primal_con_dual_var[CI{
        SAF{Float64},
        MOI.LessThan{Float64},
    }(
        1,
    )][1]
    @test MOI.get(dual_model, MOI.VariableName(), vi_2) == "dualvar_lessthan_1"
    # Query constraint names
    ci_1 = primal_dual_map.primal_var_dual_con[VI(1)]
    ci_2 = primal_dual_map.primal_var_dual_con[VI(2)]
    @test MOI.get(dual_model, MOI.ConstraintName(), ci_1) == "dualcon_x1"
    @test MOI.get(dual_model, MOI.ConstraintName(), ci_2) == "dualcon_x2"
end


text:
using ..Tagging, ..Traits, ..Var, ..Space
# @inline hasintervene(ω) = hastag(ω, Val{:intervene})
@inline tagintervene(ω, intervention) = tag(ω, (intervene = intervention,), mergetags)
@inline (x::Intervened)(ω) = x.x(tagintervene(ω, x.i))
@inline (x::Intervened{X, <:HigherIntervention})(ω) where X =
  x.x(tagintervene(ω, x.i(ω)))

@inline function passintervene(traits, i::Intervention{X, V}, x::X, ω) where {X, V}
  # If the variable ` x` to be applied to ω is the variable to be replaced
  # then replace it, and apply its replacement to ω instead.
  # Othewise proceed as normally
  if i.x == x
    i.v(ω)
  else
    ctxapply(traits, x, ω)
  end
end

function passintervene(traits,
                       i::Union{MultiIntervention{Tuple{Intervention{X1, V1},
                                                        Intervention{X2, V2}}},
                                MultiIntervention{Tuple{Intervention{X1, V1},
                                                        Intervention{X2, V2},
                                                        Intervention{X3, V3}}},
                                MultiIntervention{Tuple{Intervention{X1, V1},
                                                        Intervention{X2, V2},
                                                        Intervention{X3, V3},
                                                        Intervention{X4, V4}}},
                                MultiIntervention{Tuple{Intervention{X1, V1},
                                                        Intervention{X2, V2},
                                                        Intervention{X3, V3},
                                                        Intervention{X4, V4},
                                                        Intervention{X5, V5}}}},
                       x::Union{X1, X2, X3, X4, X5},
                       ω) where {X1, X2, X3, X4, X5, V1, V2, V3, V4, V5}
  if x == i.is[1].x
    i.is[1].v(ω)
  elseif x == i.is[2].x
    i.is[2].v(ω)
  elseif length(i.is) >= 3 && x == i.is[3].x
    i.is[3].v(ω)
  elseif length(i.is) >= 4 && x == i.is[4].x
    i.is[4].v(ω)
  elseif length(i.is) >= 5 && x == i.is[5].x
    i.is[5].v(ω)
  else
    ctxapply(traits, x, ω)
  end
end

function passintervene(traits,
                       i::MultiIntervention{XS},
                       x,
                       ω) where XS
  index = 1;
  while (index <= length(i.is) && x != i.is[index].x)
    index += 1
  end
  if (index <= length(i.is))
    i.is[index].v(ω)
  else
    ctxapply(traits, x, ω)
  end
end

# We only consider intervention if the intervention types match
@inline passintervene(traits, i::AbstractIntervention, x, ω) =
  ctxapply(traits, x, ω)

(f::Vari)(traits::trait(Intervene), ω::AbstractΩ) =
  passintervene(traits, ω.tags.intervene, f, ω)