text:
#
# We use these definitions because when we have other SIMD operations with masks
# LLVM optimizes the masks better.
function truncate_mask!(instrs, input, W, suffix, reverse_load::Bool = false)
  mtyp_input = "i$(max(8,nextpow2(W)))"
  mtyp_trunc = "i$(W)"
  if reverse_load
    bitreverse = "i$(W) @llvm.bitreverse.i$(W)"
    decl = "declare $bitreverse(i$(W))"
    bitrevmask = "bitrevmask.$(suffix)"
    if mtyp_input == mtyp_trunc
      str = """
          %$(bitrevmask) = call $(bitreverse)($(mtyp_trunc) %$(input))
          %mask.$(suffix) = bitcast $mtyp_input %$(bitrevmask) to <$W x i1>
      """
    else
      str = """
          %masktrunc.$(suffix) = trunc $mtyp_input %$input to $mtyp_trunc
          %$(bitrevmask) = call $(bitreverse)($(mtyp_trunc) %masktrunc.$(suffix))
          %mask.$(suffix) = bitcast $mtyp_trunc %$(bitrevmask) to <$W x i1>
      """
    end
  else
    decl = ""
    if mtyp_input == mtyp_trunc
      str = "%mask.$(suffix) = bitcast $mtyp_input %$input to <$W x i1>"
    else
      str = "%masktrunc.$(suffix) = trunc $mtyp_input %$input to $mtyp_trunc\n%mask.$(suffix) = bitcast $mtyp_trunc %masktrunc.$(suffix) to <$W x i1>"
    end
  end
  push!(instrs, str)
  decl
end
function zext_mask!(instrs, input, W, suffix)
  mtyp_input = "i$(max(8,nextpow2(W)))"
  mtyp_trunc = "i$(W)"
  str = if mtyp_input == mtyp_trunc
    "%res.$(suffix) = bitcast <$W x i1> %$input to $mtyp_input"
  else
    "%restrunc.$(suffix) = bitcast <$W x i1> %$input to $mtyp_trunc\n%res.$(suffix) = zext $mtyp_trunc %restrunc.$(suffix) to $mtyp_input"
  end
  push!(instrs, str)
end
function binary_mask_op_instrs(W, op)
  mtyp_input = "i$(max(8,nextpow2(W)))"
  instrs = String[]
  truncate_mask!(instrs, '0', W, 0)
  truncate_mask!(instrs, '1', W, 1)
  push!(instrs, "%combinedmask = $op <$W x i1> %mask.0, %mask.1")
  zext_mask!(instrs, "combinedmask", W, 1)
  push!(instrs, "ret $mtyp_input %res.1")
  join(instrs, "\n")
end
function binary_mask_op(W, U, op, evl::Symbol = Symbol(""))
  instrs = binary_mask_op_instrs(W, op)
  mask = Expr(:curly, evl === Symbol("") ? :Mask : :EVLMask, W)
  gf = GlobalRef(Core, :getfield)
  gf1 = Expr(:call, gf, :m1, 1, false)
  gf2 = Expr(:call, gf, :m2, 1, false)
  llvmc = Expr(:call, GlobalRef(Base, :llvmcall), instrs, U, :(Tuple{$U,$U}), gf1, gf2)
  call = Expr(:call, mask, llvmc)
  evl === Symbol("") ||
    push!(call.args, Expr(:call, evl, :($gf(m1, :evl)), :($gf(m2, :evl))))
  Expr(:block, Expr(:meta, :inline), call)
end

@inline data(m::Mask) = getfield(m, :u)
@inline data(m::EVLMask) = getfield(m, :u)
@inline Base.convert(::Type{Mask{W,U}}, m::EVLMask{W,U}) where {W,U} =
  Mask{W,U}(getfield(m, :u))
for (f, op, evl) ∈ [
  (:vand, "and", :min),
  (:vor, "or", :max),
  (:vxor, "xor", Symbol("")),
  (:veq, "icmp eq", Symbol("")),
  (:vne, "icmp ne", Symbol("")),
]
  @eval begin
    @generated function $f(m1::AbstractMask{W,U}, m2::AbstractMask{W,U}) where {W,U}
      binary_mask_op(
        W,
        U,
        $op,
        ((m1 <: EVLMask) && (m2 <: EVLMask)) ? $(QuoteNode(evl)) : Symbol(""),
      )
    end
  end
end
for f ∈ [:vand, :vor, :vxor] # ignore irrelevant bits, so just bitcast to `Bool`
  @eval @inline $f(a::Vec{W,Bool}, b::Vec{W,Bool}) where {W} =
    vreinterpret(Bool, $f(vreinterpret(UInt8, a), vreinterpret(UInt8, b)))
end
for f ∈ [:vne, :veq] # Here we truncate.
  @eval @inline $f(a::Vec{W,Bool}, b::Vec{W,Bool}) where {W} =
    convert(Bool, $f(convert(Bit, a), convert(Bit, b)))
end

@generated function vconvert(
  ::Type{Vec{W,I}},
  m::AbstractMask{W,U},
) where {W,I<:IntegerTypesHW,U<:Union{UInt8,UInt16,UInt32,UInt64}}
  bits = 8sizeof(I)
  instrs = String[]
  truncate_mask!(instrs, '0', W, 0)
  push!(
    instrs,
    "%res = zext <$W x i1> %mask.0 to <$W x i$(bits)>\nret <$W x i$(bits)> %res",
  )
  gf = Expr(:call, GlobalRef(Core, :getfield), :m, 1, false)
  llvmc = Expr(
    :call,
    GlobalRef(Base, :llvmcall),
    join(instrs, "\n"),
    :(_Vec{$W,$I}),
    :(Tuple{$U}),
    gf,
  )
  Expr(:block, Expr(:meta, :inline), Expr(:call, :Vec, llvmc))
end

@generated function splitint(
  i::S,
  ::Type{T},
) where {S<:Base.BitInteger,T<:Union{Bool,Base.BitInteger}}
  sizeof_S = sizeof(S)
  sizeof_T = sizeof(T)
  if sizeof_T > sizeof_S
    return :(i % T)
  elseif sizeof_T == sizeof_S
    return :i
  end
  W, r = divrem(sizeof_S, sizeof_T)
  @assert iszero(r)
  vtyp = "<$W x i$(8sizeof_T)>"
  instrs = """
      %split = bitcast i$(8sizeof_S) %0 to $vtyp
      ret $vtyp %split
  """
  quote
    $(Expr(:meta, :inline))
    Vec($LLVMCALL($instrs, _Vec{$W,$T}, Tuple{$S}, i))
  end
end
@generated function fuseint(v::Vec{W,I}) where {W,I<:Union{Bool,Base.BitInteger}}
  @assert ispow2(W)
  bytes = W * sizeof(I)
  bits = 8bytes
  @assert bytes ≤ 16
  T = (I <: Signed) ? Symbol(:Int, bits) : Symbol(:UInt, bits)
  vtyp = "<$W x i$(8sizeof(I))>"
  styp = "i$(bits)"
  instrs = """
      %fused = bitcast $vtyp %0 to $styp
      ret $styp %fused
  """
  quote
    $(Expr(:meta, :inline))
    $LLVMCALL($instrs, $T, Tuple{_Vec{$W,$I}}, data(v))
  end
end

function vadd_expr(W, U, instr)
  instrs = String[]
  truncate_mask!(instrs, '0', W, 0)
  truncate_mask!(instrs, '1', W, 1)
  push!(
    instrs,
    """%uv.0 = zext <$W x i1> %mask.0 to <$W x i8>
%uv.1 = zext <$W x i1> %mask.1 to <$W x i8>
%res = $instr <$W x i8> %uv.0, %uv.1
ret <$W x i8> %res""",
  )
  Expr(
    :block,
    Expr(:meta, :inline),
    :(Vec(
      $LLVMCALL(
        $(join(instrs, "\n")),
        _Vec{$W,UInt8},
        Tuple{$U,$U},
        getfield(m1, :u),
        getfield(m2, :u),
      ),
    )),
  )
end
@generated vadd_fast(m1::AbstractMask{W,U}, m2::AbstractMask{W,U}) where {W,U} =
  vadd_expr(W, U, "add")
@generated vsub_fast(m1::AbstractMask{W,U}, m2::AbstractMask{W,U}) where {W,U} =
  vadd_expr(W, U, "sub")
@inline vadd(m1::AbstractMask{W,U}, m2::AbstractMask{W,U}) where {W,U} = vadd_fast(m1, m2)
@inline vsub(m1::AbstractMask{W,U}, m2::AbstractMask{W,U}) where {W,U} = vsub_fast(m1, m2)

@inline Base.:(&)(m::AbstractMask{W,U}, b::Bool) where {W,U} =
  Mask{W,U}(Core.ifelse(b, getfield(m, :u), zero(getfield(m, :u))))
@inline Base.:(&)(b::Bool, m::AbstractMask{W,U}) where {W,U} =
  Mask{W,U}(Core.ifelse(b, getfield(m, :u), zero(getfield(m, :u))))

@inline Base.:(|)(m::AbstractMask{W,U}, b::Bool) where {W,U} =
  Mask{W,U}(Core.ifelse(b, getfield(max_mask(Mask{W,U}), :u), getfield(m, :u)))
@inline Base.:(|)(b::Bool, m::AbstractMask{W,U}) where {W,U} =
  Mask{W,U}(Core.ifelse(b, getfield(max_mask(Mask{W,U}), :u), getfield(m, :u)))

@inline function Base.:(&)(m::EVLMask{W,U}, b::Bool) where {W,U}
  EVLMask{W,U}(
    Core.ifelse(b, getfield(m, :u), zero(getfield(m, :u))),
    Core.ifelse(b, getfield(m, :evl), 0x00000000),
  )
end
@inline function Base.:(&)(b::Bool, m::EVLMask{W,U}) where {W,U}
  EVLMask{W,U}(
    Core.ifelse(b, getfield(m, :u), zero(getfield(m, :u))),
    Core.ifelse(b, getfield(m, :evl), 0x00000000),
  )
end

@inline function Base.:(|)(m::EVLMask{W,U}, b::Bool) where {W,U}
  EVLMask{W,U}(
    Core.ifelse(b, getfield(max_mask(Mask{W,U}), :u), getfield(m, :u)),
    Core.ifelse(b, W % UInt32, getfield(m, :evl)),
  )
end
@inline function Base.:(|)(b::Bool, m::EVLMask{W,U}) where {W,U}
  EVLMask{W,U}(
    Core.ifelse(b, getfield(max_mask(Mask{W,U}), :u), getfield(m, :u)),
    Core.ifelse(b, W % UInt32, getfield(m, :evl)),
  )
end

@inline Base.:(⊻)(m::AbstractMask{W,U}, b::Bool) where {W,U} =
  Mask{W,U}(Core.ifelse(b, ~getfield(m, :u), getfield(m, :u)))
@inline Base.:(⊻)(b::Bool, m::AbstractMask{W,U}) where {W,U} =
  Mask{W,U}(Core.ifelse(b, ~getfield(m, :u), getfield(m, :u)))

@inline vshl(m::AbstractMask{W,U}, i::IntegerTypesHW) where {W,U} =
  Mask{W,U}(shl(getfield(m, :u), i))
@inline vashr(m::AbstractMask{W,U}, i::IntegerTypesHW) where {W,U} =
  Mask{W,U}(shr(getfield(m, :u), i))
@inline vlshr(m::AbstractMask{W,U}, i::IntegerTypesHW) where {W,U} =
  Mask{W,U}(shr(getfield(m, :u), i))

@inline zero_mask(::AbstractSIMDVector{W}) where {W} = Mask(zero_mask(Val(W)))
@inline zero_mask(::VecUnroll{N,W}) where {N,W} = VecUnroll{N}(Mask(zero_mask(Val(W))))
@inline max_mask(::AbstractSIMDVector{W}) where {W} = Mask(max_mask(Val(W)))
@inline max_mask(::VecUnroll{N,W}) where {N,W} = VecUnroll{N}(Mask(max_mask(Val(W))))
@inline zero_mask(::NativeTypes) = false
@inline max_mask(::NativeTypes) = true

for (U, W) in [(UInt8, 8), (UInt16, 16), (UInt32, 32), (UInt64, 64)]
  @eval @inline vany(m::AbstractMask{$W,$U}) = getfield(m, :u) != $(zero(U))
  @eval @inline vall(m::AbstractMask{$W,$U}) = getfield(m, :u) == $(typemax(U))
end
# TODO: use vector reduction intrsincs
@inline function vany(m::AbstractMask{W}) where {W}
  mm = getfield(max_mask(Val{W}()), :u)
  mu = getfield(m, :u)
  (mu & mm) !== zero(mu)
end
@inline function vall(m::AbstractMask{W}) where {W}
  mm = getfield(max_mask(Val{W}()), :u)
  mu = getfield(m, :u)
  mm & mu === mm
end
@inline vany(b::Bool) = b
@inline vall(b::Bool) = b
@inline vsum(m::AbstractMask) = count_ones(getfield(m, :u))
@inline vprod(m::AbstractMask) = vall(m)

@generated function vnot(m::AbstractMask{W,U}) where {W,U}
  mtyp_input = "i$(8sizeof(U))"
  mtyp_trunc = "i$(W)"
  instrs = String[]
  truncate_mask!(instrs, '0', W, 0)
  mask = llvmconst(W, "i1 true")
  push!(instrs, "%resvec.0 = xor <$W x i1> %mask.0, $mask")
  zext_mask!(instrs, "resvec.0", W, 1)
  push!(instrs, "ret $mtyp_input %res.1")
  quote
    $(Expr(:meta, :inline))
    Mask{$W}($LLVMCALL($(join(instrs, "\n")), $U, Tuple{$U}, getfield(m, :u)))
  end
end
@inline vnot(x::Bool) = Base.not_int(x)
# @inline Base.:(~)(m::Mask) = !m

@inline Base.count_ones(m::AbstractMask) = count_ones(getfield(m, :u))
@inline vadd(m::AbstractMask, i::IntegerTypesHW) = i + count_ones(m)
@inline vadd(i::IntegerTypesHW, m::AbstractMask) = i + count_ones(m)

@generated function vzero(::Type{M}) where {W,M<:Mask{W}}
  Expr(
    :block,
    Expr(:meta, :inline),
    Expr(:call, Expr(:curly, :Mask, W), Expr(:call, :zero, mask_type_symbol(W))),
  )
end
@generated function vzero(::Type{M}) where {W,M<:EVLMask{W}}
  Expr(
    :block,
    Expr(:meta, :inline),
    Expr(
      :call,
      Expr(:curly, :EVLMask, W),
      Expr(:call, :zero, mask_type_symbol(W)),
      0x00000000,
    ),
  )
end
@inline vzero(::Mask{W,U}) where {W,U} = Mask{W}(zero(U))
@inline vzero(::EVLMask{W,U}) where {W,U} = EVLMask{W}(zero(U), 0x00000000)
@inline Base.zero(::Type{M}) where {W,M<:AbstractMask{W}} = vzero(M)
@inline zero_mask(::Union{Val{W},StaticInt{W}}) where {W} =
  EVLMask{W}(zero(VectorizationBase.mask_type(Val{W}())), 0x00000000)

@generated function max_mask(::Union{Val{W},StaticInt{W}}) where {W}
  U = mask_type(W)
  :(EVLMask{$W,$U}($(one(U) << W - one(U)), $(UInt32(W))))
end
@inline max_mask(::Type{T}) where {T} = max_mask(pick_vector_width(T))
@generated max_mask(::Type{Mask{W,U}}) where {W,U} =
  EVLMask{W,U}(one(U) << W - one(U), W % UInt32)

@generated function valrem(::Union{Val{W},StaticInt{W}}, l::T) where {W,T<:Union{Integer,StaticInt}}
  ex = ispow2(W) ? :(l & $(T(W - 1))) : Expr(:call, Base.urem_int, :l, T(W))
  Expr(:block, Expr(:meta, :inline), ex)
end

function bzhi_quote(b)
  T = b == 32 ? :UInt32 : :UInt64
  typ = 'i' * string(b)
  instr = "i$b @llvm.x86.bmi.bzhi.$b"
  decl = "declare $instr(i$b, i$b) nounwind readnone"
  instrs = "%res = call $instr(i$b %0, i$b %1)\n ret i$b %res"
  llvmcall_expr(decl, instrs, T, :(Tuple{$T,$T}), typ, [typ, typ], [:a, :b])
end
@generated bzhi(a::UInt32, b::UInt32) = bzhi_quote(32)
@generated bzhi(a::UInt64, b::UInt64) = bzhi_quote(64)

# @generated function _mask(::Union{Val{W},StaticInt{W}}, l::I, ::True) where {W,I<:Union{Integer,StaticInt}}
#   # if `has_opmask_registers()` then we can use bitmasks directly, so we create them via bittwiddling
#   M = mask_type(W)
#   quote # If the arch has opmask registers, we can generate a bitmask and then move it into the opmask register
#     $(Expr(:meta,:inline))
#     evl = valrem(Val{$W}(), (l % $M) - one($M))
#     EVLMask{$W,$M}($(typemax(M)) >>> ($(M(8sizeof(M))-1) - evl), evl + one(evl))
#   end
# end
@generated function _mask_bzhi(::Union{Val{W},StaticInt{W}}, l::I) where {W,I<:Union{Integer,StaticInt}}
  U = mask_type_symbol(W)
  T = W > 32 ? :UInt64 : :UInt32
  quote
    $(Expr(:meta, :inline))
    m = valrem(StaticInt{$W}(), l % $T)
    m = Core.ifelse((m % UInt8) == 0x00, $W % $T, m)
    EVLMask{$W,$U}(bzhi(-1 % $T, m) % $U, m)
  end
end
# @inline function _mask_bzhi(::Union{Val{W},StaticInt{W}}, l::I) where {W,I<:Union{Integer,StaticInt}}
#   U = mask_type(StaticInt(W))
#   # m = ((l) % UInt32) & ((W-1) % UInt32)
#   m = valrem(StaticInt{W}(), l % UInt32)
#   m = Core.ifelse((m % UInt8) == 0x00, W % UInt32, m)
#   # m = Core.ifelse(zero(m) == m, -1 % UInt32, m)
#   EVLMask{W,U}(bzhi(-1 % UInt32, m) % U, m)
# end
# @inline function _mask(::Union{Val{W},StaticInt{W}}, l::I, ::True) where {W,I<:Union{Integer,StaticInt}}
#   U = mask_type(StaticInt(W))
#   m = ((l-one(l)) % UInt32) & ((W-1) % UInt32)
#   m += one(m)
#   EVLMask{W,U}(bzhi(-1 % UInt32, m) % U, m)
# end
# @generated function _mask(::Union{Val{W},StaticInt{W}}, l::I, ::True) where {W,I<:Union{Integer,StaticInt}}
#   M = mask_type_symbol(W)
#   quote
#     $(Expr(:meta,:inline))
#     evl = valrem(Val{$W}(), vsub_nw((l % $M), one($M)))
#     EVLMask{$W}(data(evl ≥ MM{$W}(0)), vadd_nw(evl, one(evl)))
#   end
# end

function mask_shift_quote(W::Int, bmi::Bool)
  if (((Sys.ARCH === :x86_64) || (Sys.ARCH === :i686))) && bmi
    W ≤ 64 && return Expr(:block, Expr(:meta, :inline), :(_mask_bzhi(StaticInt{$W}(), l)))
  end
  MT = mask_type(W)
  quote # If the arch has opmask registers, we can generate a bitmask and then move it into the opmask register
    $(Expr(:meta, :inline))
    evl = valrem(Val{$W}(), (l % $MT) - one($MT))
    EVLMask{$W,$MT}($(typemax(MT)) >>> ($(MT(8sizeof(MT)) - 1) - evl), evl + one(evl))
  end
end
@generated _mask_shift(::StaticInt{W}, l, ::True) where {W} = mask_shift_quote(W, true)
@generated _mask_shift(::StaticInt{W}, l, ::False) where {W} = mask_shift_quote(W, false)
@static if Base.libllvm_version ≥ v"12"
  function active_lane_mask_quote(W::Int)
    quote
      $(Expr(:meta, :inline))
      upper = (l % UInt32) & $((UInt32(W - 1)))
      upper = Core.ifelse(upper == 0x00000000, $(W % UInt32), upper)
      mask(Val{$W}(), 0x00000000, upper)
    end
  end
else
  function active_lane_mask_quote(W::Int)
    quote
      $(Expr(:meta, :inline))
      mask(Val{$W}(), 0x00000000, vsub_nw(l % UInt32, 0x00000001) & $(UInt32(W - 1)))
    end
  end
end
function mask_cmp_quote(W::Int, RS::Int, bmi::Bool)
  M = mask_type_symbol(W)
  bytes = min(RS ÷ W, 8)
  bytes < 4 && return mask_shift_quote(W, bmi)
  T = integer_of_bytes_symbol(bytes, true)
  quote
    $(Expr(:meta, :inline))
    evl = valrem(Val{$W}(), vsub_nw((l % $T), one($T)))
    EVLMask{$W}(data(evl ≥ MM{$W}(zero($T))), vadd_nw(evl, one(evl)))
  end
end
@generated _mask_cmp(
  ::Union{Val{W},StaticInt{W}},
  l::I,
  ::StaticInt{RS},
  ::True,
) where {W,RS,I<:Union{Integer,StaticInt}} = mask_cmp_quote(W, RS, true)
@generated _mask_cmp(
  ::Union{Val{W},StaticInt{W}},
  l::I,
  ::StaticInt{RS},
  ::False,
) where {W,RS,I<:Union{Integer,StaticInt}} = mask_cmp_quote(W, RS, false)
@generated _mask(::Union{Val{W},StaticInt{W}}, l::I, ::True) where {W,I<:Union{Integer,StaticInt}} =
  mask_shift_quote(W, true)
@generated function _mask(::Union{Val{W},StaticInt{W}}, l::I, ::False) where {W,I<:Union{Integer,StaticInt}}
  # Otherwise, it's probably more efficient to use a comparison, as this will probably create some type that can be used directly for masked moves/blends/etc
  if W > 16
    Expr(
      :block,
      Expr(:meta, :inline),
      :(_mask_shift(StaticInt{$W}(), l, has_feature(Val(:x86_64_bmi)))),
    )
    # mask_shift_quote(W)
    # elseif (Base.libllvm_version ≥ v"11") && ispow2(W)
    # elseif ((Sys.ARCH ≢ :x86_64) && (Sys.ARCH ≢ :i686)) && (Base.libllvm_version ≥ v"11") && ispow2(W)
  elseif false
    # cmpval = Base.libllvm_version ≥ v"12" ? -one(I) : zero(I)
    active_lane_mask_quote(W)
  else
    Expr(
      :block,
      Expr(:meta, :inline),
      :(_mask_cmp(
        Val{$W}(),
        l,
        simd_integer_register_size(),
        has_feature(Val(:x86_64_bmi)),
      )),
    )
  end
end
# This `mask` method returns a constant, independent of `has_opmask_registers()`; that only effects method of calculating
# the constant. So it'd be safe to bake in a value.
@inline mask(::Union{Val{W},StaticInt{W}}, L) where {W} =
  _mask(StaticInt(W), L, has_feature(Val(:x86_64_avx512f)) & ge_one_fma(cpu_name()))
@inline mask(::Union{Val{W},StaticInt{W}}, ::StaticInt{L}) where {W,L} =
  _mask(StaticInt(W), L, has_feature(Val(:x86_64_avx512f)) & ge_one_fma(cpu_name()))
@inline mask(::Type{T}, l::Union{Integer,StaticInt}) where {T} =
  _mask(pick_vector_width(T), l, has_feature(Val(:x86_64_avx512f)) & ge_one_fma(cpu_name()))

# @generated function masktable(::Union{Val{W},StaticInt{W}}, rem::Union{Integer,StaticInt}) where {W}
#     masks = Expr(:tuple)
#     for w ∈ 0:W-1
#         push!(masks.args, data(mask(Val(W), w == 0 ? W : w)))
#     end
#     Expr(
#         :block,
#         Expr(:meta,:inline),
#         Expr(:call, Expr(:curly, :Mask, W), Expr(
#             :macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)),
#             Expr(:call, :getindex, masks, Expr(:call, :+, 1, Expr(:call, :valrem, Expr(:call, Expr(:curly, W)), :rem)))
#         ))
#     )
# end

@inline tomask(m::Unsigned) = Mask{sizeof(m)}(m)
@inline tomask(m::Mask) = m
@generated function tomask(v::Vec{W,Bool}) where {W}
  usize = W > 8 ? nextpow2(W) : 8
  utyp = "i$(usize)"
  U = mask_type_symbol(W)
  instrs = String[]
  push!(instrs, "%bitvec = trunc <$W x i8> %0 to <$W x i1>")
  zext_mask!(instrs, "bitvec", W, 0)
  push!(instrs, "ret i$(usize) %res.0")
  quote
    $(Expr(:meta, :inline))
    Mask{$W}($LLVMCALL($(join(instrs, "\n")), $U, Tuple{_Vec{$W,Bool}}, data(v)))
  end
end
@inline tomask(v::AbstractSIMDVector{W,Bool}) where {W} =
  tomask(vconvert(Vec{W,Bool}, data(v)))
# @inline tounsigned(m::Mask) = getfield(m, :u)
# @inline tounsigned(m::Vec{W,Bool}) where {W} = getfield(tomask(m), :u)
@inline tounsigned(v) = getfield(tomask(v), :u)

@generated function vrem(m::Mask{W,U}, ::Type{I}) where {W,U,I<:IntegerTypesHW}
  bits = 8sizeof(I)
  instrs = String[]
  truncate_mask!(instrs, '0', W, 0)
  push!(
    instrs,
    "%res = zext <$W x i1> %mask.0 to <$W x i$(bits)>\nret <$W x i$(bits)> %res",
  )
  quote
    $(Expr(:meta, :inline))
    Vec($LLVMCALL($(join(instrs, "\n")), _Vec{$W,$I}, Tuple{$U}, data(m)))
  end
end
Vec(m::Mask{W}) where {W} = m % int_type(Val{W}())

# @inline getindexzerobased(m::Mask, i) = (getfield(m, :u) >>> i) % Bool
# @inline function extractelement(m::Mask{W}, i::Union{Integer,StaticInt}) where {W}
#     @boundscheck i > W && throw(BoundsError(m, i))
#     getindexzerobased(m, i)
# end
@generated function extractelement(v::Mask{W,U}, i::I) where {W,U,I}
  instrs = String[]
  truncate_mask!(instrs, '0', W, 0)
  push!(instrs, "%res1 = extractelement <$W x i1> %mask.0, i$(8sizeof(I)) %1")
  push!(instrs, "%res8 = zext i1 %res1 to i8\nret i8 %res8")
  instrs_string = join(instrs, "\n")
  call = :($LLVMCALL($instrs_string, Bool, Tuple{$U,$I}, data(v), i))
  Expr(:block, Expr(:meta, :inline), call)
end
@generated function insertelement(
  v::Mask{W,U},
  x::T,
  i::I,
) where {W,T,U,I<:Union{Bool,IntegerTypesHW}}
  mtyp_input = "i$(max(8,nextpow2(W)))"
  instrs = String["%bit = trunc i$(8sizeof(T)) %1 to i1"]
  truncate_mask!(instrs, '0', W, 0)
  push!(instrs, "%bitvec = insertelement <$W x i1> %mask.0, i1 %bit, i$(8sizeof(I)) %2")
  zext_mask!(instrs, "bitvec", W, 1)
  push!(instrs, "ret $(mtyp_input) %res.1")
  instrs_string = join(instrs, "\n")
  call = :(Mask{$W}($LLVMCALL($instrs_string, $U, Tuple{$U,$T,$I}, data(v), x, i)))
  Expr(:block, Expr(:meta, :inline), call)
end

# @generated function Base.isodd(i::MM{W,1}) where {W}
#     U = mask_type(W)
#     evenfirst = 0xaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa % U
#     # Expr(:block, Expr(:meta, :inline), :(isodd(getfield(i, :i)) ? Mask{$W}($oddfirst) : Mask{$W}($evenfirst)))
#     Expr(:block, Expr(:meta, :inline), :(Mask{$W}($evenfirst >> (getfield(i, :i) & 0x03))))
# end
# @generated function Base.iseven(i::MM{W,1}) where {W}
#     U = mask_type(W)
#     oddfirst = 0x55555555555555555555555555555555 % U
#     # evenfirst = oddfirst << 1
#     # Expr(:block, Expr(:meta, :inline), :(isodd(getfield(i, :i)) ? Mask{$W}($evenfirst) : Mask{$W}($oddfirst)))
#     Expr(:block, Expr(:meta, :inline), :(Mask{$W}($oddfirst >> (getfield(i, :i) & 0x03))))
# end
@inline Base.isodd(i::MM{W,1}) where {W} = Mask{W}(
  (0xaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa % mask_type(Val{W}())) >>> (getfield(i, :i) & 0x03),
)
@inline Base.iseven(i::MM{W,1}) where {W} = Mask{W}(
  (0x55555555555555555555555555555555 % mask_type(Val{W}())) >>> (getfield(i, :i) & 0x03),
)

function cmp_quote(W, cond, vtyp, T1, T2 = T1)
  instrs = String["%m = $cond $vtyp %0, %1"]
  zext_mask!(instrs, 'm', W, '0')
  push!(instrs, "ret i$(max(8,nextpow2(W))) %res.0")
  U = mask_type_symbol(W)
  quote
    $(Expr(:meta, :inline))
    Mask{$W}(
      $LLVMCALL(
        $(join(instrs, "\n")),
        $U,
        Tuple{_Vec{$W,$T1},_Vec{$W,$T2}},
        data(v1),
        data(v2),
      ),
    )
  end
end
function icmp_quote(W, cond, bytes, T1, T2 = T1)
  vtyp = vtype(W, "i$(8bytes)")
  cmp_quote(W, "icmp " * cond, vtyp, T1, T2)
end
function fcmp_quote(W, cond, T)
  vtyp = vtype(W, T === Float32 ? "float" : "double")
  cmp_quote(W, "fcmp nsz arcp contract reassoc " * cond, vtyp, T)
end
# @generated function compare(::Val{cond}, v1::Vec{W,I}, v2::Vec{W,I}) where {cond, W, I}
# cmp_quote(W, cond, sizeof(I), I)
# end
# for (f,cond) ∈ [(:(==), :eq), (:(!=), :ne), (:(>), :ugt), (:(≥), :uge), (:(<), :ult), (:(≤), :ule)]
for (f, cond) ∈ [(:veq, "eq"), (:vne, "ne")]
  @eval @generated function $f(
    v1::Vec{W,T1},
    v2::Vec{W,T2},
  ) where {W,T1<:IntegerTypesHW,T2<:IntegerTypesHW}
    if sizeof(T1) != sizeof(T2)
      return Expr(
        :block,
        Expr(:meta, :inline),
        :((v3, v4) = promote(v1, v2)),
        Expr(:call, $f, :v3, :v4),
      )
    end
    icmp_quote(W, $cond, sizeof(T1), T1, T2)
  end
end
for (f, cond) ∈ [(:vgt, "ugt"), (:vge, "uge"), (:vlt, "ult"), (:vle, "ule")]
  @eval @generated function $f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Unsigned}
    icmp_quote(W, $cond, sizeof(T), T)
  end
end
for (f, cond) ∈ [(:vgt, "sgt"), (:vge, "sge"), (:vlt, "slt"), (:vle, "sle")]
  @eval @generated function $f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Signed}
    icmp_quote(W, $cond, sizeof(T), T)
  end
end

# for (f,cond) ∈ [(:veq, "oeq"), (:vgt, "ogt"), (:vge, "oge"), (:vlt, "olt"), (:vle, "ole"), (:vne, "one")]
for (f, cond) ∈ [
  (:veq, "oeq"),
  (:vgt, "ogt"),
  (:vge, "oge"),
  (:vlt, "olt"),
  (:vle, "ole"),
  (:vne, "une"),
]
  # for (f,cond) ∈ [(:veq, "ueq"), (:vgt, "ugt"), (:vge, "uge"), (:vlt, "ult"), (:vle, "ule"), (:vne, "une")]
  @eval @generated function $f(
    v1::Vec{W,T},
    v2::Vec{W,T},
  ) where {W,T<:Union{Float32,Float64}}
    fcmp_quote(W, $cond, T)
  end
end

@inline function vgt(
  v1::AbstractSIMDVector{W,S},
  v2::AbstractSIMDVector{W,U},
) where {W,S<:SignedHW,U<:UnsignedHW}
  (v1 > zero(S)) & (vconvert(U, v1) > v2)
end
@inline function vgt(
  v1::AbstractSIMDVector{W,U},
  v2::AbstractSIMDVector{W,S},
) where {W,S<:SignedHW,U<:UnsignedHW}
  (v2 < zero(S)) | (vconvert(S, v1) > v2)
end

@inline function vge(
  v1::AbstractSIMDVector{W,S},
  v2::AbstractSIMDVector{W,U},
) where {W,S<:SignedHW,U<:UnsignedHW}
  (v1 ≥ zero(S)) & (vconvert(U, v1) ≥ v2)
end
@inline function vge(
  v1::AbstractSIMDVector{W,U},
  v2::AbstractSIMDVector{W,S},
) where {W,S<:SignedHW,U<:UnsignedHW}
  (v2 < zero(S)) | (vconvert(S, v1) ≥ v2)
end

@inline vlt(
  v1::AbstractSIMDVector{W,S},
  v2::AbstractSIMDVector{W,U},
) where {W,S<:SignedHW,U<:UnsignedHW} = vgt(v2, v1)
@inline vlt(
  v1::AbstractSIMDVector{W,U},
  v2::AbstractSIMDVector{W,S},
) where {W,S<:SignedHW,U<:UnsignedHW} = vgt(v2, v1)
@inline vle(
  v1::AbstractSIMDVector{W,S},
  v2::AbstractSIMDVector{W,U},
) where {W,S<:SignedHW,U<:UnsignedHW} = vge(v2, v1)
@inline vle(
  v1::AbstractSIMDVector{W,U},
  v2::AbstractSIMDVector{W,S},
) where {W,S<:SignedHW,U<:UnsignedHW} = vge(v2, v1)
for (op, f) ∈ [(:vgt, :(>)), (:vge, :(≥)), (:vlt, :(<)), (:vle, :(≤))]
  @eval begin
    @inline function $op(
      v1::V1,
      v2::V2,
    ) where {
      V1<:Union{IntegerTypesHW,AbstractSIMDVector{<:Any,<:IntegerTypesHW}},
      V2<:Union{IntegerTypesHW,AbstractSIMDVector{<:Any,<:IntegerTypesHW}},
    }
      V3 = promote_type(V1, V2)
      $op(itosize(v1, V3), itosize(v2, V3))
    end
    @inline function $op(v1, v2)
      v3, v4 = promote(v1, v2)
      $op(v3, v4)
    end
    @inline $op(s1::IntegerTypesHW, s2::IntegerTypesHW) = $f(s1, s2)
    @inline $op(s1::Union{Float32,Float64}, s2::Union{Float32,Float64}) = $f(s1, s2)
  end
end
for (op, f) ∈ [(:veq, :(==)), (:vne, :(≠))]
  @eval begin
    @inline $op(a, b) = ((c, d) = promote(a, b); $op(c, d))
    @inline $op(s1::NativeTypes, s2::NativeTypes) = $f(s1, s2)
  end
end

@generated function vifelse(m::AbstractMask{W,U}, v1::Vec{W,T}, v2::Vec{W,T}) where {W,U,T}
  typ = LLVM_TYPES[T]
  vtyp = vtype(W, typ)
  selty = vtype(W, "i1")
  f = "select"
  if Base.libllvm_version ≥ v"9" && ((T === Float32) || (T === Float64))
    f *= " nsz arcp contract reassoc"
  end
  instrs = String[]
  truncate_mask!(instrs, '0', W, 0)
  push!(instrs, "%res = $f $selty %mask.0, $vtyp %1, $vtyp %2\nret $vtyp %res")
  quote
    $(Expr(:meta, :inline))
    Vec(
      $LLVMCALL(
        $(join(instrs, "\n")),
        _Vec{$W,$T},
        Tuple{$U,_Vec{$W,$T},_Vec{$W,$T}},
        data(m),
        data(v1),
        data(v2),
      ),
    )
  end
end

@inline vifelse(m::Vec{W,Bool}, s1::T, s2::T) where {W,T<:NativeTypes} =
  vifelse(m, Vec{W,T}(s1), Vec{W,T}(s2))
@inline vifelse(m::AbstractMask{W}, s1::T, s2::T) where {W,T<:NativeTypes} =
  vifelse(m, Vec{W,T}(s1), Vec{W,T}(s2))
@inline vifelse(m::AbstractMask{W,U}, s1, s2) where {W,U} =
  ((x1, x2) = promote(s1, s2); vifelse(m, x1, x2))
@inline vifelse(m::AbstractMask{W}, v1::VecUnroll{N,W}, v2::VecUnroll{N,W}) where {N,W} =
  VecUnroll(fmap(vifelse, m, getfield(v1, :data), getfield(v2, :data)))

@inline Base.Bool(m::AbstractMask{1,UInt8}) = (getfield(m, :u) & 0x01) === 0x01
@inline vconvert(::Type{Bool}, m::AbstractMask{1,UInt8}) = (getfield(m, :u) & 0x01) === 0x01
@inline vifelse(m::AbstractMask{1}, s1::T, s2::T) where {T<:NativeTypes} =
  Base.ifelse(Bool(m), s1, s2)
@inline vifelse(
  f::F,
  m::AbstractSIMD{W,B},
  a::Vararg{NativeTypesV,K},
) where {F<:Function,K,W,B<:Union{Bool,Bit}} = vifelse(m, f(a...), a[K])
@inline vifelse(f::F, m::Bool, a::Vararg{NativeTypesV,K}) where {F<:Function,K} =
  ifelse(m, f(a...), a[K])

@inline vconvert(::Type{EVLMask{W,U}}, b::Bool) where {W,U} = b & max_mask(StaticInt{W}())

@inline vifelse(m::AbstractMask{W}, a::AbstractMask{W}, b::AbstractMask{W}) where {W} =
  bitselect(m, b, a)

@inline Base.isnan(v::AbstractSIMD) = v != v
@inline Base.isfinite(x::AbstractSIMD) = iszero(x - x)

@inline Base.flipsign(x::AbstractSIMD, y::AbstractSIMD) = vifelse(y > zero(y), x, -x)
for T ∈ [:Float32, :Float64]
  @eval begin
    @inline Base.flipsign(x::AbstractSIMD, y::$T) = vifelse(y > zero(y), x, -x)
    @inline Base.flipsign(x::$T, y::AbstractSIMD) = vifelse(y > zero(y), x, -x)
  end
end
@inline Base.flipsign(x::AbstractSIMD, y::Real) = ifelse(y > zero(y), x, -x)
@inline Base.flipsign(x::Real, y::AbstractSIMD) = ifelse(y > zero(y), x, -x)
@inline Base.flipsign(x::Signed, y::AbstractSIMD) = ifelse(y > zero(y), x, -x)
@inline Base.isodd(x::AbstractSIMD{W,T}) where {W,T<:Union{Integer,StaticInt}} = (x & one(T)) != zero(T)
@inline Base.iseven(x::AbstractSIMD{W,T}) where {W,T<:Union{Integer,StaticInt}} = (x & one(T)) == zero(T)

@generated function vifelse(m::Vec{W,Bool}, v1::Vec{W,T}, v2::Vec{W,T}) where {W,T}
  typ = LLVM_TYPES[T]
  vtyp = vtype(W, typ)
  selty = vtype(W, "i1")
  f = "select"
  if Base.libllvm_version ≥ v"9" && ((T === Float32) || (T === Float64))
    f *= " nsz arcp contract reassoc"
  end
  instrs = String["%mask.0 = trunc <$W x i8> %0 to <$W x i1>"]
  push!(instrs, "%res = $f $selty %mask.0, $vtyp %1, $vtyp %2\nret $vtyp %res")
  quote
    $(Expr(:meta, :inline))
    Vec(
      $LLVMCALL(
        $(join(instrs, "\n")),
        _Vec{$W,$T},
        Tuple{_Vec{$W,Bool},_Vec{$W,$T},_Vec{$W,$T}},
        data(m),
        data(v1),
        data(v2),
      ),
    )
  end
end
@inline vifelse(b::Bool, w, x) = ((y, z) = promote(w, x); vifelse(b, y, z))
@inline vifelse(b::Bool, w::T, x::T) where {T<:Union{NativeTypes,AbstractSIMDVector}} =
  Core.ifelse(b, w, x)
@inline vifelse(b::Bool, w::T, x::T) where {T<:VecUnroll} =
  VecUnroll(fmap(Core.ifelse, b, getfield(w, :data), getfield(x, :data)))

@generated function vifelse(
  m::AbstractMask{W},
  vu1::VecUnroll{Nm1,Wsplit},
  vu2::VecUnroll{Nm1,Wsplit},
) where {W,Wsplit,Nm1}
  N = Nm1 + 1
  @assert N * Wsplit == W
  U = mask_type_symbol(Wsplit)
  quote
    $(Expr(:meta, :inline))
    vifelse(vconvert(VecUnroll{$Nm1,$Wsplit,Bit,Mask{$Wsplit,$U}}, m), vu1, vu2)
  end
end

@inline vmul(v::AbstractSIMDVector, m::AbstractMask) = vifelse(m, v, zero(v))
@inline vmul(m::AbstractMask, v::AbstractSIMDVector) = vifelse(m, v, zero(v))
@inline vmul(m1::AbstractMask, m2::AbstractMask) = m1 & m2
@inline vmul(v::AbstractSIMDVector, b::Bool) = b ? v : zero(v)
@inline vmul(b::Bool, v::AbstractSIMDVector) = b ? v : zero(v)
@inline vmul(v::VecUnroll{N,W,T}, b::Bool) where {N,W,T} = b ? v : zero(v)
@inline vmul(b::Bool, v::VecUnroll{N,W,T}) where {N,W,T} = b ? v : zero(v)

@static if Base.libllvm_version ≥ v"11"
  """
    mask(::Union{StaticInt{W},Val{W}}, base, N)
    mask(base::MM{W}, N)

  The two arg (`base`, `N`) method takes a base (current index) and last index of a loop.
  Idiomatic use for three-arg version may look like

  ```julia
  using VectorizationBase
  sp = stridedpointer(x);
  for i ∈ 1:8:N
      m = mask(Val(8), (MM{8}(i),), N) # if using an integer base, also needs a `Val` or `StaticInt` to indicate size.
      v = vload(sp, (MM{8}(i),), m)
      # do something with `v`
  end
  ```
  or, a full runnable example:
  ```julia
  using VectorizationBase, SLEEFPirates
  x = randn(117); y = similar(x);
  function vexp!(y, x)
      W = VectorizationBase.pick_vector_width(eltype(x));
      L = length(y);
      spx = stridedpointer(x); spy = stridedpointer(y);
      i = MM(W, 1); # use an `MM` index.
      while (m = mask(i,L); m !== VectorizationBase.zero_mask(W))
          yᵢ = exp(vload(spx, (i,), m))
          vstore!(spy, yᵢ, (i,), m)
          i += W
      end
  end

  vexp!(y, x)
  @assert y ≈ exp.(x)

  # A sum optimized for short vectors (e.g., 10-20 elements)
  function simd_sum(x)
      W = VectorizationBase.pick_vector_width(eltype(x));
      L = length(x);
      spx = stridedpointer(x);
      i = MM(W, 1); # use an `MM` index.
      s = VectorizationBase.vzero(W, eltype(x))
      while (m = mask(i,L); m !== VectorizationBase.zero_mask(W))
          s += vload(spx, (i,), m)
          i += W
      end
      VectorizationBase.vsum(s)
  end
  # or
  function simd_sum(x)
      W = VectorizationBase.pick_vector_width(eltype(x));
      L = length(x);
      spx = stridedpointer(x);
      i = MM(W, 1); # use an `MM` index.
      s = VectorizationBase.vzero(W, eltype(x))
      cond = true
      m = mask(i,L)
      while cond
          s += vload(spx, (i,), m)
          i += W
          m = mask(i,L)
          cond = m !== VectorizationBase.zero_mask(W)
      end
      VectorizationBase.vsum(s)
  end
  ```

  ```julia
  julia> VectorizationBase.mask(Val(8), 1, 6) # starting with `i = 1`, if vector is of length 6, 6 lanes are on
  Mask{8,Bool}<1, 1, 1, 1, 1, 1, 0, 0>

  julia> VectorizationBase.mask(Val(8), 81, 93) # if `i = 81` and the vector is of length 93, we want all lanes on.
  Mask{8,Bool}<1, 1, 1, 1, 1, 1, 1, 1>

  julia> VectorizationBase.mask(Val(8), 89, 93) # But after `i += 8`, we're at `i = 89`, and now want just 5 lanes on.
  Mask{8,Bool}<1, 1, 1, 1, 1, 0, 0, 0>
  ```
  """
  @generated function mask(
    ::Union{Val{W},StaticInt{W}},
    base::T,
    N::T,
  ) where {W,T<:IntegerTypesHW}
    # declare <8 x i1> @llvm.get.active.lane.mask.v8i1.i64(i64 %base, i64 %n)
    bits = 8sizeof(T)
    typ = "i$(bits)"
    decl = "declare <$W x i1> @llvm.get.active.lane.mask.v$(W)i1.$(typ)($(typ), $(typ))"
    instrs = [
      "%m = call <$W x i1> @llvm.get.active.lane.mask.v$(W)i1.$(typ)($(typ) %0, $(typ) %1)",
    ]
    zext_mask!(instrs, 'm', W, 0)
    push!(instrs, "ret i$(max(nextpow2(W),8)) %res.0")
    args = [:base, :N]
    call = llvmcall_expr(
      decl,
      join(instrs, "\n"),
      mask_type_symbol(W),
      :(Tuple{$T,$T}),
      "i$(max(nextpow2(W),8))",
      [typ, typ],
      args,
      true,
    )
    Expr(
      :block,
      Expr(:meta, :inline),
      :(EVLMask{$W}($call, ((N % UInt32) - (base % UInt32)) + 0x00000001)),
    )
  end
  @inline mask(i::MM{W}, N::T) where {W,T<:IntegerTypesHW} =
    mask(Val{W}(), getfield(i, :i), N)
end

@generated function Base.vcat(m1::AbstractMask{_W1}, m2::AbstractMask{_W2}) where {_W1,_W2}
  if _W1 == _W2
    W = _W1
  else
    W = _W1 + _W2
    U = integer_of_bytes_symbol(W, true)
    return Expr(
      :block,
      Expr(:meta, :inline),
      :(Mask{$W}(((data(m1) % $U) << $_W2) | (data(m2) % $U))),
    )
  end
  mtyp_input = "i$(max(8,nextpow2(W)))"
  instrs = String[]
  truncate_mask!(instrs, '0', W, 0)
  truncate_mask!(instrs, '1', W, 1)

  W2 = W + W
  shuffmask = Vector{String}(undef, W2)
  for w ∈ eachindex(shuffmask)
    shuffmask[w] = string(w - 1)
  end
  mask = '<' * join(map(x -> string("i32 ", x), shuffmask), ", ") * '>'

  push!(
    instrs,
    "%combinedmask = shufflevector <$W x i1> %mask.0, <$W x i1> %mask.1, <$(W2) x i32> $mask",
  )

  mtyp_output = "i$(max(8,nextpow2(W2)))"
  zext_mask!(instrs, "combinedmask", W2, 1)
  push!(instrs, "ret $mtyp_output %res.1")
  instrj = join(instrs, "\n")
  U = mask_type_symbol(W)
  U2 = mask_type_symbol(W2)
  Expr(
    :block,
    Expr(:meta, :inline),
    :(Mask{$W2}($LLVMCALL($instrj, $U2, Tuple{$U,$U}, getfield(m1, :u), getfield(m2, :u)))),
  )
end
# @inline function Base.vcat(m1::AbstractMask{W}, m2::AbstractMask{W}) where {W}
#     U = mask_type(Val(W))
#     u1 = data(m1) % U
#     u2 = data(m2) % U
#     (u1 << W) | u2
# end

@inline ifelse(b::Bool, m1::Mask{W}, m2::Mask{W}) where {W} =
  Mask{W}(Core.ifelse(b, getfield(m1, :u), getfield(m2, :u)))
@inline ifelse(b::Bool, m1::Mask{W}, m2::EVLMask{W}) where {W} =
  Mask{W}(Core.ifelse(b, getfield(m1, :u), getfield(m2, :u)))
@inline ifelse(b::Bool, m1::EVLMask{W}, m2::Mask{W}) where {W} =
  Mask{W}(Core.ifelse(b, getfield(m1, :u), getfield(m2, :u)))
@inline ifelse(b::Bool, m1::EVLMask{W}, m2::EVLMask{W}) where {W} = EVLMask{W}(
  Core.ifelse(b, getfield(m1, :u), getfield(m2, :u)),
  Core.ifelse(b, getfield(m1, :evl), getfield(m2, :evl)),
)

@inline vconvert(::Type{<:AbstractMask{W}}, b::Bool) where {W} =
  b ? max_mask(Val(W)) : zero_mask(Val(W))
@inline vconvert(::Type{Mask{W}}, b::Bool) where {W} =
  b ? max_mask(Val(W)) : zero_mask(Val(W))
@inline vconvert(::Type{EVLMask{W}}, b::Bool) where {W} =
  b ? max_mask(Val(W)) : zero_mask(Val(W))

@inline Base.max(x::AbstractMask, y::AbstractMask) = x | y
@inline Base.min(x::AbstractMask, y::AbstractMask) = x & y
@inline Base.FastMath.max_fast(x::AbstractMask, y::AbstractMask) = x | y
@inline Base.FastMath.min_fast(x::AbstractMask, y::AbstractMask) = x & y


text:
export PermMatrix, pmrand

"""
    PermMatrix{Tv, Ti}(perm::AbstractVector{Ti}, vals::AbstractVector{Tv}) where {Tv, Ti<:Integer}
    PermMatrix(perm::Vector{Ti}, vals::Vector{Tv}) where {Tv, Ti}
    PermMatrix(ds::AbstractMatrix)

PermMatrix represents a special kind linear operator: Permute and Multiply, which means `M * v = v[perm] * val`
Optimizations are used to make it much faster than `SparseMatrixCSC`.

* `perm` is the permutation order,
* `vals` is the multiplication factor.

[Generalized Permutation Matrix](https://en.wikipedia.org/wiki/Generalized_permutation_matrix)
"""
struct PermMatrix{Tv,Ti<:Integer,Vv<:AbstractVector{Tv},Vi<:AbstractVector{Ti}} <:
       AbstractMatrix{Tv}
    perm::Vi   # new orders
    vals::Vv   # multiplied values.

    function PermMatrix{Tv,Ti,Vv,Vi}(
        perm::Vi,
        vals::Vv,
    ) where {Tv,Ti<:Integer,Vv<:AbstractVector{Tv},Vi<:AbstractVector{Ti}}
        if length(perm) != length(vals)
            throw(DimensionMismatch("permutation ($(length(perm))) and multiply ($(length(vals))) length mismatch."))
        end
        new{Tv,Ti,Vv,Vi}(perm, vals)
    end
end

function PermMatrix{Tv,Ti}(perm, vals) where {Tv,Ti<:Integer}
    PermMatrix{Tv,Ti,Vector{Tv},Vector{Ti}}(Vector{Ti}(perm), Vector{Tv}(vals))
end

function PermMatrix(
    perm::Vi,
    vals::Vv,
) where {Tv,Ti<:Integer,Vv<:AbstractVector{Tv},Vi<:AbstractVector{Ti}}
    PermMatrix{Tv,Ti,Vv,Vi}(perm, vals)
end

################# Array Functions ##################

size(M::PermMatrix) = (length(M.perm), length(M.perm))
function size(A::PermMatrix, d::Integer)
    if d < 1
        throw(ArgumentError("dimension must be ≥ 1, got $d"))
    elseif d <= 2
        return length(A.perm)
    else
        return 1
    end
end
getindex(M::PermMatrix, i::Integer, j::Integer) = M.perm[i] == j ? M.vals[i] : 0
function Base.setindex!(M::PermMatrix, val, i::Integer, j::Integer)
    if M.perm[i] == j
        @inbounds M.vals[i] = val
    else
        throw(BoundsError(M, (i, j)))
    end
end

copyto!(A::PermMatrix, B::PermMatrix) =
    (copyto!(A.perm, B.perm); copyto!(A.vals, B.vals); A)

"""
    pmrand(T::Type, n::Int) -> PermMatrix

Return random PermMatrix.
"""
function pmrand end

pmrand(::Type{T}, n::Int) where {T} = PermMatrix(randperm(n), randn(T, n))
pmrand(n::Int) = pmrand(Float64, n)

similar(x::PermMatrix{Tv,Ti}) where {Tv,Ti} =
    PermMatrix{Tv,Ti}(copy(x.perm), similar(x.vals))
similar(x::PermMatrix{Tv,Ti}, ::Type{T}) where {Tv,Ti,T} =
    PermMatrix{T,Ti}(copy(x.perm), similar(x.vals, T))

# TODO: rewrite this
# function show(io::IO, M::PermMatrix)
#     println("PermMatrix")
#     for item in zip(M.perm, M.vals)
#         i, p = item
#         println("- ($i) * $p")
#     end
# end

######### sparse array interfaces  #########
nnz(M::PermMatrix) = length(M.vals)
nonzeros(M::PermMatrix) = M.vals
findnz(M::PermMatrix) = (collect(1:size(M, 1)), M.perm, M.vals)
dropzeros!(M::PermMatrix; trim::Bool = false) = M
isdense(::PermMatrix) = false

Base.zero(pm::PermMatrix) = PermMatrix(pm.perm, zero(pm.vals))

ext: jl


text:
Base.show(io::IO, v::LocalVar) = begin
    print(io, "LocalVar($(v.sym), is_mutable=$(v.is_mutable.x), is_shared=$(v.is_shared.x))")
end

PrettyPrint.is_atom_t(::Type{LocalVar}) = true

Base.show(io::IO, scope :: Scope) = pprint(io, scope)


text:
abstract type AbstractFEASolver end

abstract type AbstractDisplacementSolver <: AbstractFEASolver end
GPUUtils.whichdevice(s::AbstractDisplacementSolver) = whichdevice(s.u)

@params mutable struct DirectDisplacementSolver{T, dim, TP<:AbstractPenalty{T}} <: AbstractDisplacementSolver
    problem::StiffnessTopOptProblem{dim, T}
    globalinfo::GlobalFEAInfo{T}
    elementinfo::ElementFEAInfo{dim, T}
    u::AbstractVector{T}
    lhs::AbstractVector{T}
    rhs::AbstractVector{T}
    vars::AbstractVector{T}
    penalty::TP
    prev_penalty::TP
    xmin::T
end
function DirectDisplacementSolver(sp::StiffnessTopOptProblem{dim, T};
    xmin=T(1)/1000,
    penalty=PowerPenalty{T}(1),
    prev_penalty=copy(penalty),
    quad_order=default_quad_order(sp)) where {dim, T}

    elementinfo = ElementFEAInfo(sp, quad_order, Val{:Static})
    globalinfo = GlobalFEAInfo(sp)
    u = zeros(T, ndofs(sp.ch.dh))
    lhs = similar(u)
    rhs = similar(u)
    vars = fill(one(T), getncells(sp.ch.dh.grid) - sum(sp.black) - sum(sp.white))
    varind = sp.varind

    prev_penalty = @set prev_penalty.p = T(NaN)
    return DirectDisplacementSolver(sp, globalinfo, elementinfo, u, lhs, rhs, vars, penalty, prev_penalty, xmin)
end
function (s::DirectDisplacementSolver{T})(::Type{Val{safe}} = Val{false}, ::Type{newT} = T; assemble_f = true) where {T, safe, newT}
    rhs = assemble_f ? s.globalinfo.f : s.rhs
    lhs = assemble_f ? s.u : s.lhs
    globalinfo = GlobalFEAInfo(s.globalinfo.K, rhs)
    assemble!(globalinfo, s.problem, s.elementinfo, s.vars, getpenalty(s), s.xmin, assemble_f = assemble_f)
    K = globalinfo.K
    if safe
        m = meandiag(K)
        for i in 1:size(K,1)
            if K[i,i] ≈ zero(T)
                K[i,i] = m
            end
        end
    end
    try
        if T === newT
            lhs .= cholesky(Symmetric(K)) \ rhs
        else
            lhs .= cholesky(Symmetric(newT.(K))) \ newT.(rhs)
        end
    catch err
        lhs .= T(NaN)
    end

    nothing
end


text:
using Test

using Distributions, LinearAlgebra, Statistics
using EnKF

@testset "EnsembleState" begin

# Test size and initialize
States = [randn(10) for i = 1:4]

EnS = EnsembleState(States)

@test EnS.S == States

@test size(EnS) == (4, (10,))

@test mean(EnS) - (EnS.S[1]+EnS.S[2]+EnS.S[3]+EnS.S[4])/4 ==zeros(10)

EnS = EnsembleState(4, zeros(5))

@test size(EnS) == (4, (5,))

@test EnS.S[1] == zeros(5)

# Test mean and deviation
A = MvNormal(ones(5), I)
out = initialize(10, A)
out̄ = mean(out)

out̃1 = out.S[1] .- out̄
out̃2 = out.S[2] .- out̄

fluc = EnsembleState(10,zeros(5))
deviation(fluc, out)

@test fluc.S[1] == out̃1
@test fluc.S[2] == out̃2

# Test addition and substraction
ENS1 = initialize(10,4)

ENS2 = initialize(10,4)

ENS3 = EnsembleState(10,zeros(4))

ENS3 = ENS1 + ENS2
@test ENS1.S + ENS2.S == ENS3.S

ENS4 = EnsembleState(10,zeros(4))

ENS4 = ENS1- ENS2
@test ENS1.S - ENS2.S == ENS4.S
end

max_stars_repo_name: mleprovost/EnKF

max_issues_repo_name: mleprovost/EnKF

max_forks_repo_name: mleprovost/EnKF


text:
type NetworkStruct
    status::Cint
    ptr::Ptr{Void}
end


text:
# an AbstractIndex is a thing that can be used to look up ordered things by name, but that
# will also accept a position or set of positions or range or other things and pass them
# through cleanly.
abstract type AbstractIndex end

function Base.summary(idx::AbstractIndex)
    l = length(idx)
    "data frame with $l column$(l == 1 ? "" : "s")"
end
Base.summary(io::IO, idx::AbstractIndex) = print(io, summary(idx))

const ColumnIndex = Union{Signed, Unsigned, Symbol}

struct Index <: AbstractIndex   # an OrderedDict would be nice here...
    lookup::Dict{Symbol, Int}      # name => names array position
    names::Vector{Symbol}
end

function Index(names::AbstractVector{Symbol}; makeunique::Bool=false)
    u = make_unique(names, makeunique=makeunique)
    lookup = Dict{Symbol, Int}(zip(u, 1:length(u)))
    Index(lookup, u)
end
Index() = Index(Dict{Symbol, Int}(), Symbol[])
Base.length(x::Index) = length(x.names)
Base.names(x::Index) = copy(x.names)
_names(x::Index) = x.names
Base.copy(x::Index) = Index(copy(x.lookup), copy(x.names))
Base.isequal(x::AbstractIndex, y::AbstractIndex) = _names(x) == _names(y) # it is enough to check names
Base.:(==)(x::AbstractIndex, y::AbstractIndex) = isequal(x, y)

function rename!(x::Index, nms::AbstractVector{Symbol}; makeunique::Bool=false)
    if !makeunique
        if length(unique(nms)) != length(nms)
            dup = unique(nms[nonunique(DataFrame(nms=nms))])
            dupstr = join(string.(':', dup), ", ", " and ")
            msg = "Duplicate variable names: $dupstr. Pass makeunique=true" *
                  " to make them unique using a suffix automatically."
            throw(ArgumentError(msg))
        end
    end
    if length(nms) != length(x)
        throw(DimensionMismatch("Length of nms doesn't match length of x."))
    end
    make_unique!(x.names, nms, makeunique=makeunique)
    empty!(x.lookup)
    for (i, n) in enumerate(x.names)
        x.lookup[n] = i
    end
    return x
end

function rename!(x::Index, nms::AbstractVector{Pair{Symbol, Symbol}})
    xbackup = copy(x)
    processedfrom = Set{Symbol}()
    processedto = Set{Symbol}()
    toholder = Dict{Symbol,Int}()
    for (from, to) in nms
        if from ∈ processedfrom
            copy!(x.lookup, xbackup.lookup)
            x.names .= xbackup.names
            throw(ArgumentError("Tried renaming :$from multiple times."))
        end
        if to ∈ processedto
            copy!(x.lookup, xbackup.lookup)
            x.names .= xbackup.names
            throw(ArgumentError("Tried renaming to :$to multiple times."))
        end
        push!(processedfrom, from)
        push!(processedto, to)
        from == to && continue # No change, nothing to do
        if !haskey(xbackup, from)
            copy!(x.lookup, xbackup.lookup)
            x.names .= xbackup.names
            throw(ArgumentError("Tried renaming :$from to :$to, when :$from " *
                                "does not exist in the Index."))
        end
        if haskey(x, to)
            toholder[to] = x.lookup[to]
        end
        col = haskey(toholder, from) ? pop!(toholder, from) : pop!(x.lookup, from)
        x.lookup[to] = col
        x.names[col] = to
    end
    if !isempty(toholder)
        copy!(x.lookup, xbackup.lookup)
        x.names .= xbackup.names
        throw(ArgumentError("Tried renaming to :$(first(keys(toholder))), " *
                            "when it already exists in the Index."))
    end
    return x
end

rename!(f::Function, x::Index) = rename!(x, [(n=>Symbol(f(n))) for n in x.names])

Base.haskey(x::Index, key::Symbol) = haskey(x.lookup, key)
Base.haskey(x::Index, key::Integer) = 1 <= key <= length(x.names)
Base.haskey(x::Index, key::Bool) =
    throw(ArgumentError("invalid key: $key of type Bool"))
Base.keys(x::Index) = names(x)

# TODO: If this should stay 'unsafe', perhaps make unexported
function Base.push!(x::Index, nm::Symbol)
    x.lookup[nm] = length(x) + 1
    push!(x.names, nm)
    return x
end

function Base.merge!(x::Index, y::AbstractIndex; makeunique::Bool=false)
    adds = add_names(x, y, makeunique=makeunique)
    i = length(x)
    for add in adds
        i += 1
        x.lookup[add] = i
    end
    append!(x.names, adds)
    return x
end

Base.merge(x::Index, y::AbstractIndex; makeunique::Bool=false) =
    merge!(copy(x), y, makeunique=makeunique)

function Base.delete!(x::Index, idx::Integer)
    # reset the lookup's beyond the deleted item
    for i in (idx + 1):length(x.names)
        x.lookup[x.names[i]] = i - 1
    end
    delete!(x.lookup, x.names[idx])
    deleteat!(x.names, idx)
    return x
end

function Base.delete!(x::Index, nm::Symbol)
    if !haskey(x.lookup, nm)
        return x
    end
    idx = x.lookup[nm]
    return delete!(x, idx)
end

function Base.empty!(x::Index)
    empty!(x.lookup)
    empty!(x.names)
    x
end

function Base.insert!(x::Index, idx::Integer, nm::Symbol)
    if !(1 <= idx <= length(x.names)+1)
        throw(BoundsError(x, idx))
     end
    for i = idx:length(x.names)
        x.lookup[x.names[i]] = i + 1
    end
    x.lookup[nm] = idx
    insert!(x.names, idx, nm)
    x
end

@inline Base.getindex(x::AbstractIndex, idx::Bool) = throw(ArgumentError("invalid index: $idx of type Bool"))

@inline function Base.getindex(x::AbstractIndex, idx::Integer)
    if !(1 <= idx <= length(x))
        throw(BoundsError(x, idx))
    end
    Int(idx)
end

@inline function Base.getindex(x::AbstractIndex, idx::AbstractVector{Int})
    isempty(idx) && return idx
    minidx, maxidx = extrema(idx)
    if minidx < 1
        throw(BoundsError(x, idx))
    end
    if maxidx > length(x)
        throw(BoundsError(x, idx))
    end
    allunique(idx) || throw(ArgumentError("Elements of $idx must be unique"))
    idx
end

@inline function Base.getindex(x::AbstractIndex, idx::AbstractRange{Int})
    isempty(idx) && return idx
    minidx, maxidx = extrema(idx)
    if minidx < 1
        throw(BoundsError(x, idx))
    end
    if maxidx > length(x)
        throw(BoundsError(x, idx))
    end
    allunique(idx) || throw(ArgumentError("Elements of $idx must be unique"))
    idx
end

@inline Base.getindex(x::AbstractIndex, idx::AbstractRange{<:Integer}) =
    getindex(x, collect(Int, idx))
@inline Base.getindex(x::AbstractIndex, ::Colon) = Base.OneTo(length(x))
@inline Base.getindex(x::AbstractIndex, notidx::Not) =
    setdiff(1:length(x), getindex(x, notidx.skip))
@inline Base.getindex(x::AbstractIndex, idx::Between) = x[idx.first]:x[idx.last]
@inline Base.getindex(x::AbstractIndex, idx::All) =
    isempty(idx.cols) ? (1:length(x)) : union(getindex.(Ref(x), idx.cols)...)

@inline function Base.getindex(x::AbstractIndex, idx::AbstractVector{<:Integer})
    if any(v -> v isa Bool, idx)
        throw(ArgumentError("Bool values except for AbstractVector{Bool} are not allowed for column indexing"))
    end
    getindex(x, Vector{Int}(idx))
end

@inline Base.getindex(x::AbstractIndex, idx::AbstractRange{Bool}) = getindex(x, collect(idx))

@inline function Base.getindex(x::AbstractIndex, idx::AbstractVector{Bool})
    length(x) == length(idx) || throw(BoundsError(x, idx))
    findall(idx)
end

# catch all method handling cases when type of idx is not narrowest possible, Any in particular
@inline function Base.getindex(x::AbstractIndex, idxs::AbstractVector)
    isempty(idxs) && return Int[] # special case of empty idxs
    if idxs[1] isa Real
        if !all(v -> v isa Integer && !(v isa Bool), idxs)
            throw(ArgumentError("Only Integer values allowed when indexing by vector of numbers"))
        end
        return getindex(x, convert(Vector{Int}, idxs))
    end
    idxs[1] isa Symbol && return getindex(x, convert(Vector{Symbol}, idxs))
    throw(ArgumentError("idxs[1] has type $(typeof(idxs[1])); "*
                        "Only Integer or Symbol values allowed when indexing by vector"))
end

@inline function Base.getindex(x::AbstractIndex, rx::Regex)
    getindex(x, filter(name -> occursin(rx, String(name)), _names(x)))
end

# Fuzzy matching rules:
# 1. ignore case
# 2. maximum Levenshtein distance is 2
# 3. always show matches with 0 difference (wrong case)
# 4. on top of 3. do not show more than 8 matches in total
# Returns candidates ordered by (distance, name) pair
function fuzzymatch(l::Dict{Symbol, Int}, idx::Symbol)
        idxs = uppercase(string(idx))
        dist = [(REPL.levenshtein(uppercase(string(x)), idxs), x) for x in keys(l)]
        sort!(dist)
        c = [count(x -> x[1] <= i, dist) for i in 0:2]
        maxd = max(0, searchsortedlast(c, 8) - 1)
        [s for (d, s) in dist if d <= maxd]
end

@inline function lookupname(l::Dict{Symbol, Int}, idx::Symbol)
    i = get(l, idx, nothing)
    if i === nothing
        candidates = fuzzymatch(l, idx)
        if isempty(candidates)
            throw(ArgumentError("column name :$idx not found in the data frame"))
        end
        candidatesstr = join(string.(':', candidates), ", ", " and ")
        throw(ArgumentError("column name :$idx not found in the data frame; " *
                            "existing most similar names are: $candidatesstr"))
    end
    i
end

@inline Base.getindex(x::Index, idx::Symbol) = lookupname(x.lookup, idx)
@inline function Base.getindex(x::Index, idx::AbstractVector{Symbol})
    allunique(idx) || throw(ArgumentError("Elements of $idx must be unique"))
    [lookupname(x.lookup, i) for i in idx]
end

# Helpers

function add_names(ind::Index, add_ind::AbstractIndex; makeunique::Bool=false)
    u = names(add_ind)

    seen = Set(_names(ind))
    dups = Int[]

    for i in 1:length(u)
        name = u[i]
        in(name, seen) ? push!(dups, i) : push!(seen, name)
    end
    if length(dups) > 0
        if !makeunique
            dupstr = join(string.(':', unique(u[dups])), ", ", " and ")
            msg = "Duplicate variable names: $dupstr. Pass makeunique=true" *
                  " to make them unique using a suffix automatically."
            throw(ArgumentError(msg))
        end
    end
    for i in dups
        nm = u[i]
        k = 1
        while true
            newnm = Symbol("$(nm)_$k")
            if !in(newnm, seen)
                u[i] = newnm
                push!(seen, newnm)
                break
            end
            k += 1
        end
    end

    return u
end

@inline parentcols(ind::Index) = Base.OneTo(length(ind))
@inline parentcols(ind::Index, cols) = ind[cols]

### SubIndex of Index. Used by SubDataFrame, DataFrameRow, and DataFrameRows

struct SubIndex{I<:AbstractIndex,S<:AbstractVector{Int},T<:AbstractVector{Int}} <: AbstractIndex
    parent::I
    cols::S # columns from idx selected in SubIndex
    remap::T # reverse mapping from cols to their position in the SubIndex
end

SubIndex(parent::AbstractIndex, ::Colon) = parent

Base.copy(x::SubIndex) = Index(_names(x))

@inline parentcols(ind::SubIndex) = ind.cols

Base.@propagate_inbounds parentcols(ind::SubIndex, idx::Union{Integer,AbstractVector{<:Integer}}) =
    ind.cols[idx]
Base.@propagate_inbounds parentcols(ind::SubIndex, idx::Bool) =
    throw(ArgumentError("column indexing with Bool is not allowed"))

Base.@propagate_inbounds function parentcols(ind::SubIndex, idx::Symbol)
    parentcol = ind.parent[idx]
    @boundscheck begin
        remap = ind.remap
        remap[parentcol] == 0 && throw(ArgumentError("$idx not found"))
    end
    return parentcol
end

Base.@propagate_inbounds parentcols(ind::SubIndex, idx::AbstractVector{Symbol}) =
    [parentcols(ind, i) for i in idx]

Base.@propagate_inbounds parentcols(ind::SubIndex, idx::Regex) =
    [parentcols(ind, i) for i in _names(ind) if occursin(idx, String(i))]

Base.@propagate_inbounds parentcols(ind::SubIndex, ::Colon) = ind.cols

Base.@propagate_inbounds parentcols(ind::SubIndex, idx::Not) = parentcols(ind, ind[idx])

function SubIndex(parent::AbstractIndex, cols::AbstractUnitRange{Int})
    l = last(cols)
    f = first(cols)
    if !checkindex(Bool, Base.OneTo(length(parent)), cols)
        throw(BoundsError(parent, cols))
    end
    remap = (1:l) .- f .+ 1
    SubIndex(parent, cols, remap)
end

function SubIndex(parent::AbstractIndex, cols::AbstractVector{Int})
    ncols = length(parent)
    remap = zeros(Int, ncols)
    for (i, col) in enumerate(cols)
        if !(1 <= col <= ncols)
            throw(BoundsError(parent, cols))
        end
        if remap[col] != 0
            throw(ArgumentError("duplicate selected column detected"))
        end
        remap[col] = i
    end
    SubIndex(parent, cols, remap)
end

@inline SubIndex(parent::AbstractIndex, cols::ColumnIndex) =
    throw(ArgumentError("cols argument must be a vector (got $cols)"))

Base.@propagate_inbounds SubIndex(parent::AbstractIndex, cols) =
    SubIndex(parent, parent[cols])

Base.length(x::SubIndex) = length(x.cols)
Base.names(x::SubIndex) = copy(_names(x))
_names(x::SubIndex) = view(_names(x.parent), x.cols)

function Base.haskey(x::SubIndex, key::Symbol)
    haskey(x.parent, key) || return false
    pos = x.parent[key]
    remap = x.remap
    checkbounds(Bool, remap, pos) || return false
    remap[pos] > 0
end

Base.haskey(x::SubIndex, key::Integer) = 1 <= key <= length(x)
Base.haskey(x::SubIndex, key::Bool) =
    throw(ArgumentError("invalid key: $key of type Bool"))
Base.keys(x::SubIndex) = names(x)

function Base.getindex(x::SubIndex, idx::Symbol)
    remap = x.remap
    remap[x.parent[idx]]
end

function Base.getindex(x::SubIndex, idx::AbstractVector{Symbol})
    allunique(idx) || throw(ArgumentError("Elements of $idx must be unique"))
    [x[i] for i in idx]
end


text:
module InputTexts

Base.Experimental.@optlevel 1

using Redux
using CImGui
import ..ReduxImGui: get_label, get_string

# actions
abstract type AbstractInputTextAction <: AbstractSyncAction end

"""
    Rename(label, new_label)
Change label to `new_label`.
Note that renaming may also change the identifier, please refer to `help?> CImGui.PushID`.
"""
struct Rename <: AbstractInputTextAction
    label::String
    new_label::String
end

"""
    ChangeSize(label, size)
Change buffer size to `size`.
"""
struct ChangeSize <: AbstractInputTextAction
    label::String
    size::Int
end

# state
struct State <: AbstractImmutableState
    label::String
    size::Int
    buffer::String
    flags::CImGui.ImGuiInputTextFlags_
end
State(
    label::AbstractString;
    size = 280,
    buffer = String(fill(UInt8(0), size)),
    flags = CImGui.ImGuiInputTextFlags_None,
) = State(label, size, buffer, flags)

# reducers
input_text(state::AbstractState, action::AbstractAction) = state
input_text(state::Vector{<:AbstractState}, action::AbstractAction) = state
input_text(state::Dict{String,<:AbstractState}, action::AbstractAction) = state

function input_text(state::Dict{String,State}, action::AbstractInputTextAction)
    s = Dict{String,State}()
    for (k,v) in state
        s[k] = get_label(v) == action.label ? input_text(v, action) : v
    end
    return s
end

input_text(s::State, a::Rename) = State(a.new_label, s.size, s.buffer, s.flags)
function input_text(s::State, a::ChangeSize)
    new_buf = String(fill(UInt8(0), a.size - s.size))
    return State(s.label, a.size, string(s.buffer, new_buf), s.flags)
end

# helper
"""
    InputText(store::AbstractStore, get_state=Redux.get_state) -> Bool
Return `true` when triggered.
"""
function InputText(store::AbstractStore, get_state=Redux.get_state)
    s = get_state(store)
    is_triggered = CImGui.InputText(get_label(s), s.buffer, s.size, s.flags)
    return is_triggered
end

get_label(s::State) = s.label
get_size(s::State) = s.size
get_string(s::State) = split(s.buffer, '\0') |> first

end # module


text:
# Script to create the test wavefunction
# Dynamics not required: imprint an exact GPE vortex core onto a Thomas-Fermi state

using SpecialFunctions, Plots, Plots.Measures, JLD2, LaTeXStrings, QuadGK, VortexDistributions
using QuantumFluidSpectra
import VortexDistributions: Λ, ψa
gr()

function new_plot(;size=(400,200))
    fs,tfs,lw = 12,8,1.5
    p = plot(legend=:topright,lw=lw,size=size,frame=:box,
    foreground_color_legend = nothing,
    xtickfontsize=tfs,ytickfontsize=tfs,xguidefontsize=fs,
    tickfont=font(fs, "Times"),
    yguidefontsize=fs,legendfontsize=tfs)
    return p
end

## square domain
L = 22
N = 512
μ = 30.0
g = 0.1
w = 1

## useful methods
V(x,y) = 0.5*w^2*(x^2 + y^2)
ψ0(x,y,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y)/μ,0.0)+im*0.0)
healinglength(x,y,μ,g) = 1/sqrt(g*abs2(ψ0(x,y,μ,g)))
R(w) = sqrt(2*μ/w^2)
Rtf = R(1)
ξ = healinglength(0,0.,μ,g)
x = LinRange(-L/2,L/2,N); y = x
X = (x,y)
n0 = abs2.(ψ0(0,0,μ,g))
ek_unit = pi*n0*Rtf^3/(2μ)

##
ψ = ψ0.(x,y',μ,g)
# heatmap(x,y,abs.(ψ),aspect_ratio=1.0)

## Imprint charge 1 vortex at the origin
pv = PointVortex(0,0.,1)
# v1 = ScalarVortex(ξ,pv)

# check ansatz core agrees exactly with TFv
v1 = ScalarVortex(Ansatz(ψa,ξ,Λ),pv)
psi = Torus(copy(ψ),x,y)
vortex!(psi,v1)
gr()
heatmap(x,y,abs.(psi.ψ),aspect_ratio=1.0)

## k grid for spectra, and reference values
kx,ky = kvecs((L,L),(N,N))
K = (kx,ky)
kR = 2*pi/Rtf
kxi = 2*pi/ξ

kmin = 0.1kR
kmax = 1.2kxi
Np = 1000
k = LinRange(kmin,kmax,Np)

εki = incompressible_spectrum(k,Psi(psi.ψ,X,K))
# calc spec and save
# ψv = psi.ψ
# psiv = Psi(ψv,X,K)
# εki = incompressible_spectrum(k,psiv)
# @save joinpath(@__DIR__,"test_psi.jld2") ψv x y εki

#@load joinpath(@__DIR__,"test_psi.jld2") ψv x y εki

## Fig 3 (a) ui power spectrum plot
pgfplotsx()
ep3a = new_plot()

# Analytic form homog. [PRX]
f(x) = x*(besseli(0,x)*besselk(1,x)-besseli(1,x)*besselk(0,x))
FΛ(x) = f(x/2/Λ)^2/x
plot!(k*ξ,FΛ.(k*ξ),line=(1,:blue,0.8),label=false)

# analytic spectra: vortex in trap
Tint(x,a,b)= x*sqrt((1-x^2)/(x^2+b^2))*besselj1(a*x)
Tv(a,b) = quadgk(x->Tint(x,a,b),0.,1.)[1]
Mv(x,y) = x*abs2(Tv(x,y))
ekin_v_a(k) = Mv(k*Rtf,ξ/Λ/Rtf)*ek_unit
plot!(k*ξ,ekin_v_a.(k)*ξ/μ,line=(4,:pink),label=false)

plot!(k*ξ,εki*ξ/μ,scale=:log10,line=(1,:red),legend=:topright,label=false,minorticks=true,grid=false)
xlims!(extrema(k*ξ)[1],ξ*kxi*1.15)
ylims!(1e-3,1e2)
xlabel!(L"k \xi")
ylabel!(L"\varepsilon^i_k(k)\xi/\mu",grid=false)
vline!([1.0],label=false,line=(1,:gray))
vline!([2π],label=false,line=(1,:gray))
vline!([kR*ξ],label=false,line=(1,:gray))
ya = 6e-2
annotate!(0.115,ya,text(L"\frac{2\pi}{R}",10))
annotate!(.83,ya,text(L"\frac{1}{\xi}",10))
annotate!(4.9,ya,text(L"\frac{2\pi}{\xi}",10))
annotate!(5,35,text(L"(a)",10))
# joinpath(@__DIR__,"test_data/central_vortex_εik.pdf") |>  savefig
# gr()
plot!()

##
# εki
# ekin_v_a.(k)
# maximum relative error as function of grid size
# @256 => 1.65979
# @512 => 0.07953
# @1024 => 0.01276
# @2048 => 0.00337
maximum(abs.((εki .- ekin_v_a.(k))./ekin_v_a.(k)))

errk = abs.((εki .- ekin_v_a.(k))./ekin_v_a.(k))
plot(k*ξ, errk,yscale=:log10)

## Two point velocity correlation [check units!]
# pgfplotsx()
ep3b = new_plot()
r = LinRange(0,2*Rtf,1000)
gi = gv(r,k,ekin_v_a.(k))
plot!(r/ξ,gi,label=false,grid=false)
vline!([Rtf/ξ],line=(1,:gray),label=false)
ylabel!(L"g_k^i(r)")
plot!(yticks=0:0.5:1.0)
xlabel!(L"r/\xi")
xlims!(extrema(r/ξ)...)
annotate!(45,.75,text(L"R",10))
annotate!(82,0.92,text(L"(b)",10))

## combined plot
l = @layout [a; b]
pc = plot(ep3a,ep3b,layout=l,size=(430,350))
# "2d_trapvtf_combined.pdf" |> savefig
plot!()


text:
using Base: Int32
struct Diffresgrids
    R1::Matrix{Float64}
    R2::Matrix{Float64}
    R3::Matrix{Float64}
    I1::Matrix{Float64}
    I2::Matrix{Float64}
    I3::Matrix{Float64}
    Nx::Int32
    Ny::Int32
    Nx1::Int32
    Ny1::Int32
    Nx2::Int32
    Ny2::Int32
    Nx3::Int32
    Ny3::Int32

end

function multigrid_grids(crs_grd::Grid, walls::Vector{Wall}, deltas::Vector{Float64}, corners::Matrix{Float64}, m::Vector{Int})
    # fine_grid = Grid(crs_grd, Rectangle, thickness of step (one element vector), x/y position of bounding box (larger) and top right corner, 2, 3+)

    # defining the new grids by a change of 2 in resolution
    fine_grid = Grid(crs_grd, walls, deltas, corners, 2, m)
    finer_grid = Grid(fine_grid, walls, deltas, corners, 2, m)
    finest_grid = Grid(finer_grid, walls, deltas, corners, 2, m)

    # defining projections at each of the resolutions
    P1 = crs_grd.projection
    P1inv = crs_grd.inverse_projection

    P2 = fine_grid.projection
    P2inv = fine_grid.inverse_projection

    P3 = finer_grid.projection
    P3inv = finer_grid.inverse_projection

    P4 = finest_grid.projection
    P4inv = finest_grid.inverse_projection

    # obtaining restriction/interpolation matrices in Cartesian coordinates at all resolutions
    # need to resolve difference in resolution between given grids and the resolutions
    # calculated in the transformation_matrices function
    Nx = crs_grd.Nx
    Ny = crs_grd.Ny
    R1, I1, Nxnew, Nynew = transformation_matrices(Nx, Ny)

    Nx1 = fine_grid.Nx
    Ny1 = fine_grid.Ny
    R2, I2, Nx1new, Ny1new = transformation_matrices(Nx1, Ny1)

    Nx2 = finer_grid.Nx
    Ny2 = finer_grid.Ny
    R3, I3, Nx2new, Ny2new = transformation_matrices(Nx2, Ny2)

    # calculating generalized restriction/interpolation matrices
    R1p = P2 * R1 * P1inv
    R2p = P3 * R2 * P2inv
    R3p = P4 * R3 * P3inv

    I1p = P1 * I1 * P2inv
    I2p = P2 * I2 * P3inv
    I3p = P3 * I3 * P4inv

    return diffresgrids(R1p, R2p, R3p, I1p, I2p, I3p, Nx, Ny, Nx1, Ny1, Nx2, Ny2, Nx3, Ny3)
end


text:
abstract type EnergyComponent{D} end

"""
    get_specific_energies(energy_component::EnergyComponent) -> Vector{Float}

Return the specific energies of the particles in the energy component.
"""
function get_specific_energies end

"""
    get_specific_energies(
        energy_component::EnergyComponent,
        (i, j)::Tuple{Unsigned,Unsigned},
    ) -> (Float, Float)

Return the specific energies of the particles with the given indices.
"""
function get_specific_energies(
    energy_component::EnergyComponent,
    (i, j)::Tuple{Unsigned,Unsigned},
)
    specific_energies = get_specific_energies(energy_component)
    specific_energies[i], specific_energies[j]
end

"""
    get_specific_energy(energy_component::EnergyComponent, i::Unsigned) -> Float

Return the specific energy of the particle with the given index.
"""
function get_specific_energy(energy_component::EnergyComponent, i::Unsigned)
    specific_energies = get_specific_energies(energy_component)
    specific_energies[i]
end

"""
    initderivatives(energy_component::EnergyComponent{ED}) -> ED

Creates initial derivatives for the energy component.
"""
initderivatives(::EnergyComponent{Nothing}) = nothing

"""
    initderivatives!(
        derivatives::ED,
        ::EnergyComponent{ED},
    )

Resets the given derivatives for the energy component.
"""
function initderivatives!(::Nothing, ::EnergyComponent{Nothing}) end

"""
    updatederivatives!(
        energy_derivatives::ED,
        (i, j)::Tuple{Unsigned,Unsigned},
        vᵣ::Number,
        (Δaᵢ_pressure, Δaⱼ_pressure, Δa_diffusion)::Tuple{Number,Number,Number},
        energy_component::EnergyComponent{ED},
    )

Update energy component derivatives for particles `i` and `j`.
"""
function updatederivatives!(
    ::Nothing,
    ::Tuple{Unsigned,Unsigned},
    ::Number,
    ::Tuple{Number,Number,Number},
    ::EnergyComponent{Nothing},
) end

"""
    updatederivative!(
        energy_derivatives::ED,
        i::Unsigned,
        vᵣ::Number,
        (Δaᵢ_pressure, Δaⱼ_pressure, Δa_diffusion)::Tuple{Number,Number,Number},
        energy_component::EnergyComponent{ED},
    )

Update energy component derivative for particle `i`.
"""
function updatederivative!(
    ::Nothing,
    ::Unsigned,
    ::Number,
    ::Tuple{Number,Number,Number},
    ::EnergyComponent{Nothing},
) end

"""
    evolvevariables!(
        energy_component::EnergyComponent{ED},
        energy_derivatives::ED,
        Δt::Number,
    )

Evolve energy component with the given derivatives and time step.
"""
function evolvevariables!(::EnergyComponent{Nothing}, ::Nothing, ::Number) end

include("energy/standard.jl")
include("energy/isothermal.jl")

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
module COMPUTE

#Find your own path where you save the julia core file
include("/Users/yuguangchen/Projects/Tex/new structure/simulation/non-zero phi3/s matrix generating/all in one/core.jl")

using PyCall
#If you use Python 3, replace all "cPickle" to "pickle"
@pyimport pickle

#define global variables
N = 200
x = linspace(0,2*pi,N)
y = linspace(0,2*pi,N)

#dump the new data to filename
function mypickle(filename, obj)
    out = open(filename,"w")
    pickle.dump(obj, out, protocol=pickle.HIGHEST_PROTOCOL)
    close(out)
end

#read saved s matrix from the filename. You have to use myunpickle(\path\filename.extension) to properly call it.
function myunpickle(filename)
    r = nothing
    @pywith pybuiltin("open")(filename,"rb") as f begin
        r = pickle.load(f)
    end
    return r
end

#find your own file of "S.pickle" file, which specifies the s matrix
#Be aware of the naming system. The example uses particular the symmtric one near level crossing, and candidate number 1.
s = myunpickle("/Users/yuguangchen/Projects/Tex/new structure/simulation/non-zero phi3/s matrix generating/all in one/S matrix/S_weyl_1.pickle")

#energy independent or energy slowly varying S matrix

#Summation of small scale of the invariant
#Be aware of the naming system when you save the array.
function summation(phi3::Float64)
    InvarSum = CORE.InvarSum1Grid(x,y,phi3,s::Array{Complex{Float64},2})
    mypickle("/Users/yuguangchen/Projects/Tex/new structure/simulation/non-zero phi3/s matrix generating/all in one/Summation/Summation_S_weyl_1.pickle", InvarSum)
end

#Integrated summation over 2 phases
function chernsum(phi3::Float64)
    Chern = CORE.ChernSum(phi3,s::Array{Complex{Float64},2})
    mypickle("/Users/yuguangchen/Projects/Tex/new structure/simulation/non-zero phi3/s matrix generating/all in one/ChernSum/Chernsum_S_sym_weyl_1.pickle", Chern)
end

#Non-topologogical contribution -- the large scale conductant, integrated over energy
function landauer(phi3::Float64)
    conductance = CORE.avrgInvarGrid(x,y,phi3,s::Array{Complex{Float64},2})
    mypickle("/Users/yuguangchen/Projects/Tex/new structure/simulation/non-zero phi3/s matrix generating/all in one/Landauer/Landauer_S_weyl_1.pickle", conductance)
end

#Landuaer conductance integrated with phases.
function chernlandauer(phi3::Float64)
    Chern = CORE.ChernL(phi3,s::Array{Complex{Float64},2})
    mypickle("/Users/yuguangchen/Projects/Tex/new structure/simulation/non-zero phi3/s matrix generating/all in one/ChernLandauer/Chernlandauer_S_sym_weyl_1.pickle", Chern)
end

#Berry curvature of discrete energy levels only
function berry(phi3::Float64)
    berry = CORE.BerryBranchGrid(x,y,phi3,s::Array{Complex{Float64},2})
    mypickle("/Users/yuguangchen/Projects/Tex/new structure/simulation/non-zero phi3/s matrix generating/all in one/Berry/Berry_S_weyl_1.pickle", berry)
end

#curvature of discrete energy -- berry -- integrated over 2 phases
function chernberry(phi3::Float64)
    Chern = CORE.ChernBerry(x,y,phi3,s::Array{Complex{Float64},2})
    mypickle("/Users/yuguangchen/Projects/Tex/new structure/simulation/non-zero phi3/s matrix generating/all in one/ChernBerry/Chernberry_S_sym_weyl_1.pickle", Chern)
end

#Lowest ABS
function lowestabs(phi3::Float64)
    abs = CORE.LowestABSGrid(x,y,phi3,s::Array{Complex{Float64},2})
    mypickle("/Users/yuguangchen/Projects/Tex/new structure/simulation/non-zero phi3/s matrix generating/all in one/LowestABS/Lowestabs_S_weyl_1.pickle", abs)
end

#HighestABS
function highestabs(phi3::Float64)
    abs = CORE.HighestABSGrid(x,y,phi3,s::Array{Complex{Float64},2})
    mypickle("/Users/yuguangchen/Projects/Tex/new structure/simulation/non-zero phi3/s matrix generating/all in one/HighestABS/Highestabs_S_weyl_1.pickle", abs)
end

#energy dependent S matrix

#invariant part of a small-scale energy dependent S matriz (scale<1)
function invariantsmall(phi3::Float64,scale::Float64)
    curvature = CORE.IntInvarESGrid(x,y,phi3,scale,s::Array{Complex{Float64},2})
    mypickle("/Users/yuguangchen/Projects/Tex/new structure/simulation/non-zero phi3/s matrix generating/all in one/S_e small_curv/Invariantsmall_S_sym_weyl_1.pickle", curvature)
end

#invariant part of a large-scale energy dependent S matriz (scale>1)
function invariantlarge(phi3::Float64,scale::Float64)
    curvature = CORE.IntInvarEBGrid(x,y,phi3,scale,s::Array{Complex{Float64},2})
    mypickle("/Users/yuguangchen/Projects/Tex/new structure/simulation/non-zero phi3/s matrix generating/all in one/S_e_large_curv/Invariantlarge_S_weyl_1.pickle", curvature)
end

#Energy dependent curvature integrated over phases -- Chern number
function cherns_e(phi3::Float64, scale::Float64)
    Chern = CORE.ChernE(phi3, scale, s::Array{Complex{Float64},2})
    mypickle("/Users/yuguangchen/Projects/Tex/new structure/simulation/non-zero phi3/s matrix generating/all in one/ChernS_e/Cherns_e_S_weyl_1.pickle", Chern)
end

end

max_stars_repo_path: all in one/Compute.jl

max_issues_repo_path: all in one/Compute.jl

max_forks_repo_path: all in one/Compute.jl


text:
using Logging, Random
using JuMP, Gurobi
using DecisionProgramming

Random.seed!(13)

const N = 4
const L = [1]
const R_k = [k + 1 for k in 1:N]
const A_k = [(N + 1) + k for k in 1:N]
const F = [2*N + 2]
const T = [2*N + 3]
const L_states = ["high", "low"]
const R_k_states = ["high", "low"]
const A_k_states = ["yes", "no"]
const F_states = ["failure", "success"]
const c_k = rand(N)
const b = 0.03
fortification(k, a) = [c_k[k], 0][a]

@info("Creating the influence diagram.")
S = States([
    (length(L_states), L),
    (length(R_k_states), R_k),
    (length(A_k_states), A_k),
    (length(F_states), F)
])
C = Vector{ChanceNode}()
D = Vector{DecisionNode}()
V = Vector{ValueNode}()
X = Vector{Probabilities}()
Y = Vector{Consequences}()

for j in L
    I_j = Vector{Node}()
    X_j = zeros(S[I_j]..., S[j])
    X_j[1] = rand()
    X_j[2] = 1.0 - X_j[1]
    push!(C, ChanceNode(j, I_j))
    push!(X, Probabilities(j, X_j))
end

for j in R_k
    I_j = L
    x, y = rand(2)
    X_j = zeros(S[I_j]..., S[j])
    X_j[1, 1] = max(x, 1-x)
    X_j[1, 2] = 1.0 - X_j[1, 1]
    X_j[2, 2] = max(y, 1-y)
    X_j[2, 1] = 1.0 - X_j[2, 2]
    push!(C, ChanceNode(j, I_j))
    push!(X, Probabilities(j, X_j))
end

for (i, j) in zip(R_k, A_k)
    I_j = [i]
    push!(D, DecisionNode(j, I_j))
end

for j in F
    I_j = L ∪ A_k
    x, y = rand(2)
    X_j = zeros(S[I_j]..., S[j])
    for s in paths(S[A_k])
        d = exp(b * sum(fortification(k, a) for (k, a) in enumerate(s)))
        X_j[1, s..., 1] = max(x, 1-x) / d
        X_j[1, s..., 2] = 1.0 - X_j[1, s..., 1]
        X_j[2, s..., 1] = min(y, 1-y) / d
        X_j[2, s..., 2] = 1.0 - X_j[2, s..., 1]
    end
    push!(C, ChanceNode(j, I_j))
    push!(X, Probabilities(j, X_j))
end

for j in T
    I_j = A_k ∪ F
    Y_j = zeros(S[I_j]...)
    for s in paths(S[A_k])
        cost = sum(-fortification(k, a) for (k, a) in enumerate(s))
        Y_j[s..., 1] = cost + 0
        Y_j[s..., 2] = cost + 100
    end
    push!(V, ValueNode(j, I_j))
    push!(Y, Consequences(j, Y_j))
end

validate_influence_diagram(S, C, D, V)
sort!.((C, D, V, X, Y), by = x -> x.j)

P = DefaultPathProbability(C, X)
U = DefaultPathUtility(V, Y)

@info("Creating the decision model.")
U⁺ = PositivePathUtility(S, U)
model = Model()
z = DecisionVariables(model, S, D)
x_s = PathCompatibilityVariables(model, z, S, P, probability_cut = false)
EV = expected_value(model, x_s, U⁺, P)
@objective(model, Max, EV)

@info("Starting the optimization process.")
optimizer = optimizer_with_attributes(
    () -> Gurobi.Optimizer(Gurobi.Env()),
    "IntFeasTol"      => 1e-9,
)
set_optimizer(model, optimizer)
optimize!(model)

@info("Extracting results.")
Z = DecisionStrategy(z)

@info("Printing decision strategy:")
print_decision_strategy(S, Z)

@info("Printing state probabilities:")
sprobs = StateProbabilities(S, P, Z)
print_state_probabilities(sprobs, L)
print_state_probabilities(sprobs, R_k)
print_state_probabilities(sprobs, A_k)
print_state_probabilities(sprobs, F)

@info("Computing utility distribution.")
udist = UtilityDistribution(S, P, U, Z)

@info("Printing utility distribution.")
print_utility_distribution(udist)

@info("Printing statistics")
print_statistics(udist)


text:
# ===================================================
# The following schemes are described in the paper
#   Weak exponential schemes for stochastic
#   differential equations with additive noise
#   CM Mora - IMA journal of numerical analysis, 2005
# ===================================================

_drift_jacobian(model::AbstractSDE{D,M}, t0, x0) where {D,M} = ForwardDiff.jacobian(x->drift(model, t0, x), x0)
_drift_jacobian(model::AbstractSDE{1,M}, t0, x0) where {M} = ForwardDiff.derivative(x->drift(model, t0, x), x0)

function step(model::StateIndependentDiffusion, scheme::EulerExponential1, t0, x0, Δw)
  Jμ = _drift_jacobian(model, t0, x0)
  μ = drift(model, t0, x0)
  σ = diffusion(model, t0, x0)
  x = exp(Jμ * scheme.Δt) * (x0 + (μ - Jμ * x0) * scheme.Δt + σ * Δw)
  x
end

function step(model::StateIndependentDiffusion, scheme::EulerExponential2, t0, x0, Δw)
  Jμ = _drift_jacobian(model, t0, x0)
  μ = drift(model, t0, x0)
  σ = diffusion(model, t0, x0)
  x = (I - Jμ * scheme.Δt) \ (x0 + (μ - Jμ * x0) * scheme.Δt + σ * Δw)
  x
end

function step(model::StateIndependentDiffusion, scheme::EulerExponential3, t0, x0, Δw)
  Jμ = _drift_jacobian(model, t0, x0)
  expmJμ = exp(0.5 * Jμ * scheme.Δt)
  μ = drift(model, t0, x0)
  σ = diffusion(model, t0, x0)
  x = expmJμ * (expmJμ * x0 + (μ - Jμ * x0) * scheme.Δt + σ * Δw)
  x
end

function _normal_transition_params(model, scheme::EulerExponential3, t0, x0, x1)
  Jμ = _drift_jacobian(model, t0, x0)
  expmJμ = exp(0.5 * Jμ * scheme.Δt)
  μ = expmJμ * (expmJμ * x0 + (drift(model, t0, x0) - Jμ * x0) * scheme.Δt)
  σ = expmJμ * diffusion(model, t0, x0)
  Σ = scheme.Δt * σ * σ'
  μ, Σ
end


text:
function bisect_left(A, p, r, v)
	i = p
	if p < r
	   q = floor(Int, (p + r) / 2)
	   if v <= A[q]
		   i = bisect_left(A, p, q, v)
	   else
		   i = bisect_left(A, q + 1, r, v)
	   end
	end
	return i
end

function bisect_right(A, p, r, v)
	i = p
	if p < r
	   q = floor(Int, (p + r) / 2)
	   if v < A[q]
		   i = bisect_right(A, p, q, v)
	   else
		   i = bisect_right(A, q + 1, r, v)
	   end
	end
	return i
end

function find_median(a, low, upp)
	i_low = bisect_left(a, 1, length(a) + 1, low)
	i_upp = bisect_right(a, 1, length(a) + 1, upp)

	mid = floor(Int, (i_upp + i_low) / 2)
	if (i_upp - i_low) % 2 == 0
		return (a[mid] + a[mid - 1]) / 2
	else
		return a[mid]
	end
end


text:
function Draw_P_vs_samples(samples,P_samples,file_name::String)
    @show n_chains = length(samples)
    @show n_samples = length(samples[1])
    @show n_features = length(samples[1][1])

    sample_mat = Array{Float64,2}(undef,n_chains*n_samples,n_features)
    P_mat =  Array{Float64,1}(undef,n_chains*n_samples)

    for i = 1:n_chains
        for j = 1:n_samples
            sample_mat[j+(i-1)*n_samples,:] = samples[i][j]
            P_mat[j+(i-1)*n_samples] = P_samples[i][j]
        end
    end

    plot(;layout=n_features,legends=false)
    for i = 1:n_features
        plot!(sample_mat[:,i],P_mat,subplot=i)
    end

    savefig("./plots/"*file_name*"_P_vs_Sample.svg")
end

function Draw_trace_plot(samples,file_name::String)
    @show n_chains = length(samples)
    @show n_samples = length(samples[1])
    @show n_parameters = length(samples[1][1])

    #plot(layout = n_parameters,legends=false)
    for i = 1:n_parameters
        plot(;legends=false)
        for j = 1:n_chains
            x_val = 1:n_samples
            y_val = [samples[j][k][i] for k = 1:n_samples]

            plot!(x_val,y_val)
        end
        savefig("./plots/trace_plot/"*file_name*"_$(i).svg")
    end
end

function Draw_histogram(samples,file_name::String)
    @show n_chains = length(samples)
    @show n_samples = length(samples[1])
    @show n_features = length(samples[1][1])

    sample_mat = Array{Float64,2}(undef,n_chains*n_samples,n_features)

    for i = 1:n_chains
        for j = 1:n_samples
            sample_mat[j+(i-1)*n_samples,:] = samples[i][j]
        end
    end

    histogram(sample_mat,layout=n_features,legend=false)
    savefig("./plots/"*file_name*".svg")
end

function Find_max_sample(samples,P_samples)
    @show n_chains = length(samples)
    @show n_samples = length(samples[1])
    @show n_features = length(samples[1][1])

    max_P_vec = Array{Float64,1}(undef,n_chains)
    max_ind_vec = Array{Int,1}(undef,n_chains)

    for i = 1:n_chains
        P_vex = [P_samples[i][j] for j = 1:n_samples]
        max_P_vec[i] = maximum(P_vex)
        max_ind_vec[i] = argmax(P_vex)
    end

    max_chains = argmax(max_P_vec)
    max_sample = max_ind_vec[max_chains]
    max_features = copy(samples[max_chains][max_sample])

    return max_features
end


text:
#funtio

function updateC!(C,c_1,c_μ,c_c,p_c,w,hsig,aux)
    C=(1-c_1-c_μ)*C+
            c_1*(p_c*p_c'+(1-hsig)*c_c*(2-c_c)*C)+
            c_μ*aux*Diagonal(w)*aux'
end

function new_cma_es(f, x0 ; tol=1e-5,maxiter=1e3)
      N=length(x0)
      λ=4+floor(3*log(N))
      m=x0
      σ=0.3
      μ=floor(λ/2)
      μ_Int=Integer(μ)
      w=(log(μ+0.5).-log.(1:μ))
      w=w/sum(w)
      μ_w=(sum(w)^2)/sum(w.^2)
      c_c=(4+μ_w/N)/(N+4+2*μ_w/N)
      c_σ=(μ_w+2)/(N+μ_w+5)
      c_1=2/((N+1.3)^2+μ_w)
      c_μ=min(1-c_1,
              2*(μ_w-2+1/μ_w)/((N+2)^2+μ_w) )
      damp_σ=1+c_σ+ 2*max(0.0,
                  sqrt( (μ_w-1)/(N+1) )-1    )
      p_c=zeros(N,1)
      p_σ=p_c
      B=I(N)
      D=ones(N,1)
      C=(B*Diagonal(D.^2)*B')
      invsqrtC=B*Diagonal( D.^-1 )*B'
      eigeneval=0
      χ_N=sqrt(N)*(1-1/(4*N) +1/(21*N^2) )
      count=0
      normdiff=Inf
      N=length(x0)
      fx=zeros(Integer(λ))

      C=(Diagonal(ones(N)))
      X=(zeros(Integer(λ),N))

      @inbounds for count in 1:maxiter

          @inbounds @simd for i in 1:Integer(λ)
            X[i,:]=(σ.*B*(D.*randn(N,1))+m)'
            fx[i]=f(X[i,:])
            count=count+1
          end

          temp=@view (X[sortperm(fx),:])[1:μ_Int,:]
          dm=(temp'-repeat(m,1,μ_Int))*w
          m=m+dm
          p_σ=(1-c_σ)*p_σ+
                  ((sqrt(c_σ*(2-c_σ)*μ_w)*invsqrtC)*dm)/σ

          hsig=norm(p_σ)/sqrt(1-(1-c_σ)^(2*count/λ))χ_N<1.4+2/(N+1)

          p_c=(1-c_c)*p_c+(hsig*sqrt(c_c*(2-c_c)*μ)/σ) *dm

          aux=(1/σ).*(temp'-repeat(m,1,μ_Int) )

          updateC!(C,c_1,c_μ,c_c,p_c,w,hsig,aux)

          σ=σ*exp( (c_σ/damp_σ)*(norm(p_σ)/χ_N-1) )

          if count-eigeneval>λ/(c_1+c_μ)/N/10
            eigeneval=count
            C=(C+C')./2
            eig=eigen(C)
            invsqrtC=eig.vectors'*Diagonal( eig.values.^-0.5 )*eig.vectors
          end
          count=count+1
          normdiff=norm(dm)
          if normdiff< tol #|| maximum(D)>1e7*minimum(D)
                  break
          end
          end
          return(f(X[1,:]),X[1,:])

end