text:
using Test

## General ######################################################
sp_2_cqs = Dict{AbstractString, Symbol}("AAA" => :cqs_0,
                                 "AA" => :cqs_1,
                                 "A" => :cqs_2,
                                 "BBB" => :cqs_3,
                                 "BB" => :cqs_4,
                                 "other" => :cqs_5,
                                 "na" => :na)

## SCR Parameters ###############################################
## SCR BSCR =====================================================
ds2 =
  Dict{Symbol, Any}(:corr => [1.00  0.25  0.25  0.25  0.25;
                               0.25  1.00  0.25  0.25  0.50;
                               0.25  0.25  1.00  0.25  0.00;
                               0.25  0.25  0.25  1.00  0.00;
                               0.25  0.50  0.00  0.00  1.00],
                     :mdl => [:S2Mkt, :S2Def, :S2Life,
                              :S2Health, :S2NonLife],
                     :coc => 0.06)

## SCR Op =======================================================
ds2_op = Dict{Symbol, Float64}(:bscr => 0.3,
                               :cost => 0.25,
                               :prem => 0.04,
                               :prem_py => 1.2,
                               :tp => 0.0045)

## SCR Mkt ======================================================
## The indices 21, 31, 41, 12, 13, 14 need to be adjusted
corr_mkt_raw =
  [1.00  0.00  0.00  0.00  0.25  0.00;
   0.00  1.00  0.75  0.75  0.25  0.00;
   0.00  0.75  1.00  0.50  0.25  0.00;
   0.00  0.75  0.50  1.00  0.25  0.00;
   0.25  0.25  0.25  0.25  1.00  0.00;
   0.00  0.00  0.00  0.00  0.00  1.00]

ind_adj = Array{Vector{Int}}(undef, 0)
push!(ind_adj, [1,2], [1,3], [1,4], [2,1], [3,1], [4,1])

function corrmkt(raw::Matrix, ind_adj, updown::Symbol)
  a = (updown == :up ? 0.0 : 0.5)
  corr = deepcopy(raw)
  for 𝑝𝑎𝑖𝑟 ∈ ind_adj
    corr[𝑝𝑎𝑖𝑟[1], 𝑝𝑎𝑖𝑟[2]]= a
  end
  return corr
end

ds2_mkt =
  Dict{Symbol, Any}(:raw => corr_mkt_raw,
                    :adj => ind_adj,
                    :corr => corrmkt,
                    :mdl => [:S2MktInt, :S2MktEq, :S2MktProp,
                             :S2MktSpread, :S2MktFx, :S2MktConc])

## SCR Mkt Int --------------------------------------------------
## only first 9 deltas, SCR.5.25
δ_down =
  [-0.75, -0.65, -0.56, -0.50, -0.46, -0.42, -0.39, -0.36, -0.33]
δ_up =
  [ 0.70,  0.70,  0.64,  0.59,  0.55,  0.52,  0.49,  0.47,  0.44]

ds2_mkt_int =
  Dict{Symbol, Any}(:shock =>
                      Dict{Symbol, Any}(:spot_up => δ_up,
                                        :spot_down => δ_down),
                    :spot_up_abs_min => 0.01)

## SCR Mkt Eq ---------------------------------------------------
eq_shock_raw_1, eq_shock_raw_2 = -0.39, -0.49
eq_adj_sym = -0.075
eq_shock =
  Dict{Symbol, Float64}(:type_1 => eq_shock_raw_1 + eq_adj_sym,
                        :type_2 => eq_shock_raw_2 + eq_adj_sym)

ds2_mkt_eq = Dict{Symbol, Any}(:shock => eq_shock,
                               :corr => [1  0.75;  0.75  1])

## SCR Mkt all --------------------------------------------------
ds2_mkt_all =
  Dict{Symbol, Dict{Symbol, Any}}(:mkt => ds2_mkt,
                                  :mkt_int => ds2_mkt_int,
                                  :mkt_eq => ds2_mkt_eq)

## SCR Def ======================================================
ds2_def = Dict{Symbol, Any}(:corr => [1.00  0.75; 0.75  1.00])

## SCR Def Type 1  ----------------------------------------------
df_def_1 =
  DataFrame(quant = ["ProbGD"],
            cqs_0 =      [0.00002],
            cqs_1 =      [0.00010],
            cqs_2 =      [0.00050],
            cqs_3 =      [0.00240],
            cqs_4 =      [0.01200],
            cqs_5 =      [0.04200],
            cqs_6 =      [0.04200],
            cqs_unrat =  [0.04200])
df_def_scr =
  DataFrame(range = [:low, :medium, :high],
            multiplier = [3.0, 5.0, NaN],
            threshold_upper = [0.07, 0.20, 1.00])

ds2_def_1 = Dict{Symbol, Any}(:prob => df_def_1,
                              :scr_par => df_def_scr)

## SCR Def all  -------------------------------------------------
ds2_def_all = Dict{Symbol, Any}(:def => ds2_def,
                                :def_1 => ds2_def_1)

## SCR Life =====================================================
ds2_life =
  Dict{Symbol, Any}(
  :corr => [ 1.00  -0.25   0.25   0.00   0.25   0.00   0.25;
             -0.25   1.00   0.00   0.25   0.25   0.25   0.00;
             0.25   0.00   1.00   0.00   0.50   0.00   0.25;
             0.00   0.25   0.00   1.00   0.50   0.00   0.25;
             0.25   0.25   0.50   0.50   1.00   0.50   0.25;
             0.00   0.25   0.00   0.00   0.50   1.00   0.00;
             0.25   0.00   0.25   0.25   0.25   0.00   1.00],
   :mdl => [:S2LifeQx, :S2LifePx, :S2LifeMorb, :S2LifeSx,
            :S2LifeCost,  :S2LifeRevision,  :S2LifeCat])

@test ds2_life[:corr] == ds2_life[:corr]'

## SCR Life Mort ------------------------------------------------
ds2_life_qx = Dict{Symbol, Any}(:shock => Dict(:qx => 0.15))

## SCR Life Long ------------------------------------------------
ds2_life_px = Dict{Symbol, Any}(:shock => Dict(:px => -0.20))

## SCR Life Lapse -----------------------------------------------
ds2_life_sx =
  Dict{Symbol, Any}(:shock => Dict(:sx_up => 0.5,
                                   :sx_down => -0.5,
                                   :sx_mass_other => 0.4,
                                   :sx_mass_pension => 0.7),
                    :shock_param =>
                      Dict(:sx_down_threshold => -0.2))
## SCR Life Expense ---------------------------------------------
ds2_life_cost =
  Dict{Symbol, Any}(:shock => Dict(:cost => 0.10),
                    :shock_param => Dict(:infl => 0.01))

## SCR Life Cat -------------------------------------------------
ds2_life_cat = Dict{Symbol, Any}(:shock => Dict(:cat => 0.15))

## SCR all ======================================================
ds2_life_all = Dict{Symbol, Any}(:life => ds2_life,
                                 :life_qx => ds2_life_qx,
                                 :life_px => ds2_life_px,
                                 :life_sx => ds2_life_sx,
                                 :life_cost => ds2_life_cost,
                                 :life_cat => ds2_life_cat)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "x264"
version = v"2020.07.14"

# Collection of sources required to build x264
sources = [
    ArchiveSource("https://code.videolan.org/videolan/x264/-/archive/db0d417728460c647ed4a847222a535b00d3dbcb/x264-db0d417728460c647ed4a847222a535b00d3dbcb.tar.gz",
                  "b79b7038ce083f152fdd35da2f0d770ac9189fa2319fd09012567bc3a33737af"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/x264-*
if [[ "${target}" == x86_64* ]] || [[ "${target}" == i686* ]]; then
    apk add nasm
    export AS=nasm
else
    export AS="${CC}"
fi
./configure --prefix=$prefix --host=$target --enable-shared --enable-pic --disable-static
# Remove unsafe compilation flag
sed -i 's/ -ffast-math//g' config.mak
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("x264", :x264),
    LibraryProduct("libx264", :libx264),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

ext: jl

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# This should be run at all cli top function
function _run_cli_os()
    try
        # backup
        update_remote_backup(;force = false, loglevel = :upload)
    catch ignored; end
end


text:
module StaticTaylorSeries

export StaticTaylor

"""
Static Taylor series with variable V, length N (i.e. order N-1) and element type T.
"""
struct StaticTaylor{N,T}
    coeffs::NTuple{N,T}
end

function coeffstring(t::StaticTaylor, i, variable=:x)

    if i == 1  # order 0
        return string(t.coeffs[i])

    elseif i == 2  # order 1
        return string(t.coeffs[i], variable)

    else
        return string(t.coeffs[i], variable, "^", i-1)
    end
end

function print_taylor(io::IO, t::StaticTaylor, variable=:x)

    print(io, "(" * join([coeffstring(t, i, variable) for i in 1:length(t.coeffs)], " + ") * ")")

end

function Base.show(io::IO, t::StaticTaylor)
    print_taylor(io, t)
end

function Base.show(io::IO, t::StaticTaylor{N,T}) where {N, T<:StaticTaylor}
    print_taylor(io, t, :y)
end

#StaticTaylor(iterable...) = StaticTaylor(SVector(iterable...))

StaticTaylor{N}(v::NTuple{N,T}) where {N,T} = StaticTaylor(v)

#StaticTaylor{N}(iterable) where {N} = StaticTaylor{N}(iterable)

StaticTaylor{N}(iterable...) where {N} = StaticTaylor{N}(iterable)

# StaticTaylor(v) = StaticTaylor(v)

StaticTaylor(iterable...) = StaticTaylor{length(iterable)}(iterable)

import Base:getindex, length, eltype

getindex(s::StaticTaylor, i::Integer) = s.coeffs[i+1]

length(s::StaticTaylor{N,T}) where {N,T} = N

eltype(s::StaticTaylor{N,T}) where {N,T} = T

import Base: +, -, * ,^

function +(s::StaticTaylor{N,T}, t::StaticTaylor{N,T}) where {N,T}
    return StaticTaylor(s.coeffs .+ t.coeffs)
end

function +(s::StaticTaylor{N,T}, α::Real) where {N,T}
    return StaticTaylor{N,T}(ntuple(i -> i == 1 ? s.coeffs[1] + α : s.coeffs[i], N))
end

+(α::Real, s::StaticTaylor) = s + α

-(s::StaticTaylor) = StaticTaylor(.-(s.coeffs))

function -(s::StaticTaylor{N,T}, t::StaticTaylor{N,T}) where {N,T}
    return StaticTaylor(s.coeffs .- t.coeffs)
end

-(s::StaticTaylor, α::Real) = s + (-α)

-(α::Real, s::StaticTaylor) = -(s - a)

Base.literal_pow(::typeof(^), x::StaticTaylor, ::Val{p}) where {p} = x^p

^(x::StaticTaylor, n::Integer) = Base.power_by_squaring(x, n)

# function *(s::StaticTaylor{N,T}, t::StaticTaylor{N,T}) where {N,T}
#     v = SVector(ntuple(k->sum(s[i]*t[k-1-i] for i in 0:k-1), Val(N)))
#     return StaticTaylor(v)
# end

# The following is modified from StaticUnivariatePolynomials.jl
@generated function Base.:*(p1::StaticTaylor{N,T}, p2::StaticTaylor{N,T}, max_degree=N) where {N, T}
    exprs = Any[nothing for i in 1:N]
    for i in 0 : N-1   # order is N-1
        for j in 0 : N-1
            k = i + j + 1  # setindex does not have offset

            if k > max_degree
                continue
            end

            if exprs[k] === nothing
                exprs[k] = :(p1[$i] * p2[$j])
            else
                exprs[k] = :(muladd(p1[$i], p2[$j], $(exprs[k])))
            end
        end
    end

    # Core.println("Generated code with N=$N:")
    # Core.println(exprs)
    # Core.println()

    return quote
        Base.@_inline_meta
        StaticTaylor{N,T}(tuple($(exprs...)))
    end
end

@generated function mult(p1::StaticTaylor{N,T}, p2::StaticTaylor{N,T}, ::Val{max_degree}) where {N, T, max_degree}
    exprs = Any[nothing for i in 1:max_degree]
    for i in 0 : N-1   # order is N-1
        for j in 0 : N-1
            k = i + j + 1  # setindex does not have offset

            if k > max_degree
                continue
            end

            if exprs[k] === nothing
                exprs[k] = :(p1[$i] * p2[$j])
            else
                exprs[k] = :(muladd(p1[$i], p2[$j], $(exprs[k])))
            end
        end
    end

    # Core.println("Generated code with N=$N:")
    # Core.println(exprs)
    # Core.println()

    return quote
        Base.@_inline_meta
        StaticTaylor{max_degree,T}(tuple($(exprs...)))
    end
end
#
# @generated function *(p1::StaticTaylor{N,T}, p2::StaticTaylor{N,T}) where {N, T}
#     exprs = Any[nothing for i in 1:N]
#     for i in 0 : N-1   # order is N-1
#         for j in 0 : N-1
#             k = i + j + 1  # setindex does not have offset
#
#             if k > N
#                 continue
#             end
#
#             if exprs[k] === nothing
#                 exprs[k] = :(p1[$i] * p2[$j])
#             else
#                 exprs[k] = :(muladd(p1[$i], p2[$j], $(exprs[k])))
#             end
#         end
#     end
#
#     # Core.println("Generated code with N=$N:")
#     # Core.println(exprs)
#     # Core.println()
#
#     return quote
#         Base.@_inline_meta
#         StaticTaylor{N,T}(tuple($(exprs...)))
#     end
# end

*(s::StaticTaylor, α::Real) = StaticTaylor(α .* s.coeffs)
*(α::Real, s::StaticTaylor) = s * α

/(s::StaticTaylor, α::Real) = StaticTaylor(s.coeffs ./ α)

end # module

max_stars_repo_name: dpsanders/StaticTaylorSeries

max_issues_repo_name: dpsanders/StaticTaylorSeries

max_forks_repo_name: dpsanders/StaticTaylorSeries


text:
"""
AU - Unrolled axis
F - Factor, step size per unroll
N - How many times is it unrolled
AV - Vectorized axis
W - vector width
M - bitmask indicating whether each factor is masked
i::I - index
"""
struct Unroll{AU,F,N,AV,W,M,I}
    i::I
end
@inline Unroll{AU,F,N,AV,W,M}(i::I) where {AU,F,N,AV,W,M,I} = Unroll{AU,F,N,AV,W,M,I}(i)

const VectorIndexCore{W} = Union{Vec{W},MM{W},Unroll{<:Any,<:Any,<:Any,<:Any,W}}
const VectorIndex{W} = Union{VectorIndexCore{W},LazyMulAdd{<:Any,<:Any,<:VectorIndexCore{W}}}

# const BoolVec = Union{Mask,VecUnroll{<:Any, <:Any, Bool, <: Mask}}

const SCOPE_METADATA = """
!1 = !{!\"noaliasdomain\"}
!2 = !{!\"noaliasscope\", !1}
!3 = !{!2}
"""
const LOAD_SCOPE_FLAGS = ", !alias.scope !3";
const STORE_SCOPE_FLAGS = ", !noalias !3";

const USE_TBAA = false
# use TBAA?
# define: LOAD_SCOPE_TBAA, LOAD_SCOPE_TBAA_FLAGS, SCOPE_METADATA, STORE_TBAA, SCOPE_FLAGS, STORE_TBAA_FLAGS
let
    LOAD_TBAA = """
    !4 = !{!"jtbaa"}
    !5 = !{!6, !6, i64 0, i64 0}
    !6 = !{!"jtbaa_arraybuf", !4, i64 0}
    """;
    LOAD_TBAA_FLAGS = ", !tbaa !5";
    global const LOAD_SCOPE_TBAA = USE_TBAA ? SCOPE_METADATA * LOAD_TBAA : SCOPE_METADATA;
    global const LOAD_SCOPE_TBAA_FLAGS = USE_TBAA ? LOAD_SCOPE_FLAGS * LOAD_TBAA_FLAGS : LOAD_SCOPE_FLAGS

    global const STORE_TBAA = USE_TBAA ? """
    !4 = !{!"jtbaa", !5, i64 0}
    !5 = !{!"jtbaa"}
    !6 = !{!"jtbaa_data", !4, i64 0}
    !7 = !{!8, !8, i64 0}
    !8 = !{!"jtbaa_arraybuf", !6, i64 0}
    """ : ""
    global const STORE_TBAA_FLAGS = USE_TBAA ? ", !tbaa !7" : ""
end

"""
An omnibus offset constructor.

The general motivation for generating the memory addresses as LLVM IR rather than combining multiple lllvmcall Julia functions is
that we want to minimize the `inttoptr` and `ptrtoint` calculations as we go back and fourth. These can get in the way of some
optimizations, such as memory address calculations.
It is particulary import for `gather` and `scatter`s, as these functions take a `Vec{W,Ptr{T}}` argument to load/store a
`Vec{W,T}` to/from. If `sizeof(T) < sizeof(Int)`, converting the `<W x \$(typ)*` vectors of pointers in LLVM to integer
vectors as they're represented in Julia will likely make them too large to fit in a single register, splitting the operation
into multiple operations, forcing a corresponding split of the `Vec{W,T}` vector as well.
This would all be avoided by not promoting/widenting the `<W x \$(typ)>` into a vector of `Int`s.

For this last issue, an alternate workaround would be to wrap a `Vec` of 32-bit integers with a type that defines it as a pointer for use with
internal llvmcall functions, but I haven't really explored this optimization.
"""
function offset_ptr(
    ::Type{T}, ind_type::Symbol, indargname, ibits::Int, W::Int, X::Int, M::Int, O::Int, forgep::Bool
) where {T}
    T_sym = JULIA_TYPES[T]
    offset_ptr(T_sym, sizeof_T, ind_type, indargname, ibits, W, X, M, O, forgep)
end
function offset_ptr(
    T_sym::Symbol, ind_type::Symbol, indargname, ibits::Int, W::Int, X::Int, M::Int, O::Int, forgep::Bool
)
    sizeof_T = JULIA_TYPE_SIZE[T_sym]
    i = 0
    Morig = M
    isbit = T_sym === :Bit
    if isbit
        typ = "i1"
        vtyp = isone(W) ? typ : "<$W x i1>"
        M = max(1, M >> 3)
        O >>= 3
        @assert ((isone(X) | iszero(X)) && (ind_type !== :Vec)) "indexing BitArrays with a vector not currently supported."
    else
        typ = LLVM_TYPES_SYM[T_sym]
        vtyp = vtype(W, typ) # vtyp is dest type
    end
    instrs = String[]
    if M == 0
        ind_type = :StaticInt
    elseif ind_type === :StaticInt
        M = 0
    end
    if isone(W)
        X = 1
    end
    if iszero(M)
        tz = intlog2(sizeof_T)
        tzf = sizeof_T
        index_gep_typ = typ
    else
        tz = min(trailing_zeros(M), 3)
        tzf = 1 << tz
        index_gep_typ = ((tzf == sizeof_T) | iszero(M)) ? typ : "i$(tzf << 3)"
        M >>= tz
    end
    # after this block, we will have a index_gep_typ pointer
    if iszero(O)
        push!(instrs, "%ptr.$(i) = inttoptr $(JULIAPOINTERTYPE) %0 to $(index_gep_typ)*"); i += 1
    else # !iszero(O)
        if !iszero(O & (tzf - 1)) # then index_gep_typ works for the constant offset
            offset_gep_typ = "i8"
            offset = O
        else # then we need another intermediary
            offset_gep_typ = index_gep_typ
            offset = O >>> tz
        end
        push!(instrs, "%ptr.$(i) = inttoptr $(JULIAPOINTERTYPE) %0 to $(offset_gep_typ)*"); i += 1
        push!(instrs, "%ptr.$(i) = getelementptr inbounds $(offset_gep_typ), $(offset_gep_typ)* %ptr.$(i-1), i32 $(offset)"); i += 1
        if forgep && iszero(M) && (iszero(X) || isone(X))
            push!(instrs, "%ptr.$(i) = ptrtoint $(offset_gep_typ)* %ptr.$(i-1) to $(JULIAPOINTERTYPE)"); i += 1
            return instrs, i
        elseif offset_gep_typ != index_gep_typ
            push!(instrs, "%ptr.$(i) = bitcast $(offset_gep_typ)* %ptr.$(i-1) to $(index_gep_typ)*"); i += 1
        end
    end
    # will do final type conversion
    if ind_type === :Vec
        if isone(M)
            indname = indargname
        else
            indname = "indname"
            constmul = llvmconst(W, "i$(ibits) $M")
            push!(instrs, "%$(indname) = mul nsw <$W x i$(ibits)> %$(indargname), $(constmul)")
        end
        push!(instrs, "%ptr.$(i) = getelementptr inbounds $(index_gep_typ), $(index_gep_typ)* %ptr.$(i-1), <$W x i$(ibits)> %$(indname)"); i += 1
        if forgep
            push!(instrs, "%ptr.$(i) = ptrtoint <$W x $index_gep_typ*> %ptr.$(i-1) to <$W x $JULIAPOINTERTYPE>"); i += 1
        elseif index_gep_typ != vtyp
            push!(instrs, "%ptr.$(i) = bitcast <$W x $index_gep_typ*> %ptr.$(i-1) to <$W x $typ*>"); i += 1
        end
        return instrs, i
    end
    if ind_type === :Integer
        if isbit
            scale = 8
            if ispow2(Morig)
                if abs(Morig) ≥ 8
                    M = Morig >> 3
                    if M != 1
                        indname = "indname"
                        push!(instrs, "%$(indname) = mul nsw i$(ibits) %$(indargname), $M")
                    else
                        indname = indargname
                    end
                else
                    shifter = 3 - intlog2(Morig)
                    push!(instrs, "%shiftedind = ashr i$(ibits) %$(indargname), $shifter")
                    if Morig > 0
                        indname = "shiftedind"
                    else
                        indname = "indname"
                        push!(instrs, "%$(indname) = mul i$(ibits) %shiftedind, -1")
                    end
                end
            else
                @assert iszero(Morig) "Scale factors on bit accesses must be 0 or a power of 2."
                indname = "0"
            end
        else
            if isone(M)
                indname = indargname
            elseif iszero(M)
                indname = "0"
            else
                indname = "indname"
                push!(instrs, "%$(indname) = mul nsw i$(ibits) %$(indargname), $M")
            end
            # TODO: if X != 1 and X != 0, check if it is better to gep -> gep, or broadcast -> add -> gep
        end
        push!(instrs, "%ptr.$(i) = getelementptr inbounds $(index_gep_typ), $(index_gep_typ)* %ptr.$(i-1), i$(ibits) %$(indname)"); i += 1
    end
    # ind_type === :Integer || ind_type === :StaticInt
    if !(isone(X) | iszero(X)) # vec
        vibytes = min(4, REGISTER_SIZE ÷ W)
        vityp = "i$(8vibytes)"
        vi = join((X*w for w ∈ 0:W-1), ", $vityp ")
        if typ !== index_gep_typ
            push!(instrs, "%ptr.$(i) = bitcast $(index_gep_typ)* %ptr.$(i-1) to $(typ)*"); i += 1
        end
        push!(instrs, "%ptr.$(i) = getelementptr inbounds $(typ), $(typ)* %ptr.$(i-1), <$W x $(vityp)> <$vityp $vi>"); i += 1
        if forgep
            push!(instrs, "%ptr.$(i) = ptrtoint <$W x $typ*> %ptr.$(i-1) to <$W x $JULIAPOINTERTYPE>"); i += 1
        end
        return instrs, i
    end
    if forgep # if forgep, just return now
        push!(instrs, "%ptr.$(i) = ptrtoint $(index_gep_typ)* %ptr.$(i-1) to $JULIAPOINTERTYPE"); i += 1
    elseif index_gep_typ != vtyp
        push!(instrs, "%ptr.$(i) = bitcast $(index_gep_typ)* %ptr.$(i-1) to $(vtyp)*"); i += 1
    end
    instrs, i
end

gep_returns_vector(W::Int, X::Int, M::Int, ind_type::Symbol) = (!isone(W) && ((ind_type === :Vec) || !(isone(X) | iszero(X))))
#::Type{T}, ::Type{I}, W::Int = 1, ivec::Bool = false, constmul::Int = 1) where {T <: NativeTypes, I <: Integer}
function gep_quote(
    ::Type{T}, ind_type::Symbol, ::Type{I}, W::Int, X::Int, M::Int, O::Int, forgep::Bool
) where {T, I}
    T_sym = JULIA_TYPES[T]
    I_sym = JULIA_TYPES[I]
    gep_quote(T_sym, ind_type, I_sym, W, X, M, O, forgep)
end
function gep_quote(
    T_sym::Symbol, ind_type::Symbol, I_sym::Symbol, W::Int, X::Int, M::Int, O::Int, forgep::Bool
)
    sizeof_T = JULIA_TYPE_SIZE[T_sym]
    sizeof_I = JULIA_TYPE_SIZE[I_sym]
    if W > 1 && ind_type !== :Vec
        X, Xr = divrem(X, sizeof_T)
        @assert iszero(Xr)
    end
    if iszero(O) && (iszero(X) | isone(X)) && (iszero(M) || ind_type === :StaticInt)
        return Expr(:block, Expr(:meta, :inline), :ptr)
    end
    ibits = 8sizeof_I
    # ::Type{T}, ind_type::Symbol, indargname = '1', ibytes::Int, W::Int = 1, X::Int = 1, M::Int = 1, O::Int = 0, forgep::Bool = false
    instrs, i = offset_ptr(T_sym, ind_type, '1', ibits, W, X, M, O, true)
    ret = Expr(:curly, :Ptr, T_sym)
    lret = JULIAPOINTERTYPE
    if gep_returns_vector(W, X, M, ind_type)
        ret = Expr(:curly, :_Vec, W, ret)
        lret = "<$W x $lret>"
    end

    args = Expr(:curly, :Tuple, Expr(:curly, :Ptr, T_sym))
    largs = String[JULIAPOINTERTYPE]
    arg_syms = Union{Symbol,Expr}[:ptr]

    if !(iszero(M) || ind_type === :StaticInt)
        push!(arg_syms, Expr(:call, :data, :i))
        if ind_type === :Integer
            push!(args.args, I_sym)
            push!(largs, "i$(ibits)")
        else
            push!(args.args, Expr(:curly, :_Vec, W, I_sym))
            push!(largs, "<$W x i$(ibits)>")
        end
    end
    push!(instrs, "ret $lret %ptr.$(i-1)")
    llvmcall_expr("", join(instrs, "\n"), ret, args, lret, largs, arg_syms)
end

@generated function gep(ptr::Ptr{T}, i::I) where {I <: Integer, T <: NativeTypes}
    gep_quote(T, :Integer, I, 1, 1, 1, 0, true)
end
@generated function gep(ptr::Ptr{T}, ::StaticInt{N}) where {N, T <: NativeTypes}
    gep_quote(T, :StaticInt, Int, 1, 1, 0, N, true)
end
@generated function gep(ptr::Ptr{T}, i::LazyMulAdd{M,O,I}) where {T <: NativeTypes, I <: Integer, O, M}
    gep_quote(T, :Integer, I, 1, 1, M, O, true)
end
@generated function gep(ptr::Ptr{T}, i::Vec{W,I}) where {W, T <: NativeTypes, I <: Integer}
    gep_quote(T, :Vec, I, W, 1, 1, 0, true)
end
@generated function gep(ptr::Ptr{T}, i::LazyMulAdd{M,O,Vec{W,I}}) where {W, T <: NativeTypes, I <: Integer, M, O}
    gep_quote(T, :Vec, I, W, 1, M, O, true)
end
@inline gesp(ptr::AbstractStridedPointer, i) = similar_no_offset(ptr, gep(ptr, i))

function vload_quote(
    ::Type{T}, ::Type{I}, ind_type::Symbol, W::Int, X::Int, M::Int, O::Int, mask::Bool, align::Bool = false
) where {T <: NativeTypes, I <: Integer}
    T_sym = JULIA_TYPES[T]
    I_sym = JULIA_TYPES[I]
    vload_quote(T_sym, I_sym, ind_type, W, X, M, O, mask, align)
end
function vload_quote(
    T_sym::Symbol, I_sym::Symbol, ind_type::Symbol, W::Int, X::Int, M::Int, O::Int, mask::Bool, align::Bool = false
)
    isbit = T_sym === :Bit
    if !isbit && W > 1
        return vload_quote(T_sym, I_sym, ind_type, W, X, M, O, mask, align, :(_Vec{$W,$T_sym}))
    end
    jtyp = isbit ? (isone(W) ? :Bool : mask_type_symbol(W)) : T_sym
    vload_quote(T_sym, I_sym, ind_type, W, X, M, O, mask, align, jtyp)
end
function vload_split_quote(W::Int, sizeof_T::Int, mask::Bool, align::Bool)
    D, r1 = divrem(W * sizeof_T, REGISTER_SIZE)
    Wnew, r2 = divrem(W, D)
    @assert (iszero(r1) & iszero(r2)) "If loading more than a vector, Must load a multiple of the vector width."
    q = Expr(:block,Expr(:meta,:inline))
    # ind_type = :StaticInt, :Integer, :Vec
    push!(q.args, :(isplit = splitvectortotuple(StaticInt{$D}(), StaticInt{$Wnew}(), i)))
    mask && push!(q.args, :(msplit = splitvectortotuple(StaticInt{$D}(), StaticInt{$Wnew}(), m)))
    t = Expr(:tuple)
    alignval = Expr(:call, Expr(:curly, :Val, align))
    for d ∈ 1:D
        call = Expr(:call, :vload, :ptr)
        push!(call.args, Expr(:ref, :isplit, d))
        mask && push!(call.args, Expr(:ref, :msplit, d))
        push!(call.args, alignval)
        push!(t.args, call)
    end
    push!(q.args, :(VecUnroll($t)))
    q
end
function vload_quote(
    T_sym::Symbol, I_sym::Symbol, ind_type::Symbol, W::Int, X::Int, M::Int, O::Int, mask::Bool, align::Bool, ret
)
    sizeof_T = JULIA_TYPE_SIZE[T_sym]
    if W * sizeof_T > REGISTER_SIZE
        return vload_split_quote(W, sizeof_T, mask, align)
    end
    sizeof_I = JULIA_TYPE_SIZE[I_sym]
    ibits = 8sizeof_I
    if W > 1 && ind_type !== :Vec
        X, Xr = divrem(X, sizeof_T)
        @assert iszero(Xr)
    end
    instrs, i = offset_ptr(T_sym, ind_type, '1', ibits, W, X, M, O, false)
    grv = gep_returns_vector(W, X, M, ind_type)
    # considers booleans to only occupy 1 bit in memory, so they must be handled specially
    isbit = T_sym === :Bit
    if isbit
        # @assert !grv "gather's not are supported with `BitArray`s."
        mask = false # TODO: not this?
        # typ = "i$W"
        alignment = (align & (!grv)) ? cld(W,8) : 1
        typ = "i1"
    else
        alignment = (align & (!grv)) ? _get_alignment(W, T_sym) : _get_alignment(0, T_sym)
        typ = LLVM_TYPES_SYM[T_sym]
    end

    decl = LOAD_SCOPE_TBAA
    dynamic_index = !(iszero(M) || ind_type === :StaticInt)

    vtyp = vtype(W, typ)
    mask && truncate_mask!(instrs, '1' + dynamic_index, W, 0)
    if grv
        loadinstr = "$vtyp @llvm.masked.gather." * suffix(W, T_sym) * '.' * ptr_suffix(W, T_sym)
        decl *= "declare $loadinstr(<$W x $typ*>, i32, <$W x i1>, $vtyp)"
        m = mask ? m = "%mask.0" : llvmconst(W, "i1 1")
        passthrough = mask ? "zeroinitializer" : "undef"
        push!(instrs, "%res = call $loadinstr(<$W x $typ*> %ptr.$(i-1), i32 $alignment, <$W x i1> $m, $vtyp $passthrough)" * LOAD_SCOPE_TBAA_FLAGS)
    elseif mask
        suff = suffix(W, T_sym)
        loadinstr = "$vtyp @llvm.masked.load." * suff * ".p0" * suff
        decl *= "declare $loadinstr($vtyp*, i32, <$W x i1>, $vtyp)"
        push!(instrs, "%res = call $loadinstr($vtyp* %ptr.$(i-1), i32 $alignment, <$W x i1> %mask.0, $vtyp zeroinitializer)" * LOAD_SCOPE_TBAA_FLAGS)
    else
        push!(instrs, "%res = load $vtyp, $vtyp* %ptr.$(i-1), align $alignment" * LOAD_SCOPE_TBAA_FLAGS)
    end
    if isbit
        lret = string('i', max(8,W))
        if W > 1
            if W < 8
                push!(instrs, "%resint = bitcast <$W x i1> %res to i$(W)")
                push!(instrs, "%resfinal = zext i$(W) %resint to i8")
            else
                push!(instrs, "%resfinal = bitcast <$W x i1> %res to i$(W)")
            end
        else
            push!(instrs, "%resfinal = zext i1 %res to i8")
        end
        push!(instrs, "ret $lret %resfinal")
    else
        lret = vtyp
        push!(instrs, "ret $vtyp %res")
    end
    args = Expr(:curly, :Tuple, Expr(:curly, :Ptr, T_sym))
    largs = String[JULIAPOINTERTYPE]
    arg_syms = Union{Symbol,Expr}[:ptr]
    if dynamic_index
        push!(arg_syms, :(data(i)))
        if ind_type === :Integer
            push!(args.args, I_sym)
            push!(largs, "i$(ibits)")
        else
            push!(args.args, :(_Vec{$W,$I_sym}))
            push!(largs, "<$W x i$(ibits)>")
        end
    end
    if mask
        push!(arg_syms, :(data(m)))
        push!(args.args, mask_type(W))
        push!(largs, "i$(max(8,nextpow2(W)))")
    end
    if isbit && W > 1
        quote
            $(Expr(:meta,:inline))
            Mask{$W}($(llvmcall_expr(decl, join(instrs, "\n"), ret, args, lret, largs, arg_syms, true)))
        end
    else
        llvmcall_expr(decl, join(instrs, "\n"), ret, args, lret, largs, arg_syms)
    end
end
# vload_quote(T, ::Type{I}, ind_type::Symbol, W::Int, X, M, O, mask, align = false)
@generated function vload(ptr::Ptr{T}, ::Val{A}) where {T <: NativeTypes, A}
    vload_quote(T, Int, :StaticInt, 1, 1, 0, 0, false, A)
end
@generated function vload(ptr::Ptr{T}, i::I, ::Val{A}) where {A, T <: NativeTypes, I <: IntegerTypes}
    vload_quote(T, I, :Integer, 1, 1, 1, 0, false, A)
end
@generated function vload(ptr::Ptr{T}, ::StaticInt{N}, ::Val{A}) where {A, N, T <: NativeTypes}
    vload_quote(T, Int, :StaticInt, 1, 1, 0, N, false, A)
end
@generated function vload(ptr::Ptr{T}, i::Vec{W,I}, ::Val{A}) where {A, W, T <: NativeTypes, I <: IntegerTypes}
    vload_quote(T, I, :Vec, W, 1, 1, 0, false, A)
end
@generated function vload(ptr::Ptr{T}, i::MM{W,X,I}, ::Val{A}) where {A, W, X, T <: NativeTypes, I <: IntegerTypes}
    vload_quote(T, I, :Integer, W, X, 1, 0, false, A)
end
@generated function vload(ptr::Ptr{T}, ::MM{W,X,StaticInt{N}}, ::Val{A}) where {A, W, X, T <: NativeTypes, N}
    vload_quote(T, Int, :StaticInt, W, X, 1, N, false, A)
end
@generated function vload(ptr::Ptr{T}, i::LazyMulAdd{M,O,I}, ::Val{A}) where {A, M, O, T <: NativeTypes, I <: IntegerTypes}
    vload_quote(T, I, :Integer, 1, 1, M, O, false, A)
end
@generated function vload(ptr::Ptr{T}, i::LazyMulAdd{M,O,Vec{W,I}}, ::Val{A}) where {A, W, M, O, T <: NativeTypes, I <: IntegerTypes}
    vload_quote(T, I, :Vec, W, 1, M, O, false, A)
end
@generated function vload(ptr::Ptr{T}, i::LazyMulAdd{M,O,MM{W,X,I}}, ::Val{A}) where {A, W, M, O, X, T <: NativeTypes, I <: IntegerTypes}
    vload_quote(T, I, :Integer, W, X*M, M, O, false, A)
end
@generated function vload(ptr::Ptr{T}, i::LazyMulAdd{M,O,MM{W,X,StaticInt{I}}}, ::Val{A}) where {A, W, M, O, X, T <: NativeTypes, I}
    vload_quote(T, Int, :StaticInt, W, X*M, M, O + I*M, false, A)
end

@generated function vload(ptr::Ptr{T}, i::Vec{W,I}, m::Mask{W}, ::Val{A}) where {A, W, T <: NativeTypes, I <: IntegerTypes}
    vload_quote(T, I, :Vec, W, 1, 1, 0, true, A)
end
@generated function vload(ptr::Ptr{T}, i::MM{W,X,I}, m::Mask{W}, ::Val{A}) where {A, W, X, T <: NativeTypes, I <: IntegerTypes}
    vload_quote(T, I, :Integer, W, X, 1, 0, true, A)
end
@generated function vload(ptr::Ptr{T}, ::MM{W,X,StaticInt{N}}, m::Mask{W}, ::Val{A}) where {A, W, X, T <: NativeTypes, N}
    vload_quote(T, Int, :StaticInt, W, X, 1, N, true, A)
end
@generated function vload(ptr::Ptr{T}, i::LazyMulAdd{M,O,Vec{W,I}}, m::Mask{W}, ::Val{A}) where {A, W, M, O, T <: NativeTypes, I <: IntegerTypes}
    vload_quote(T, I, :Vec, W, 1, M, O, true, A)
end
@generated function vload(ptr::Ptr{T}, i::LazyMulAdd{M,O,MM{W,X,I}}, m::Mask{W}, ::Val{A}) where {A, W, M, O, X, T <: NativeTypes, I <: IntegerTypes}
    vload_quote(T, I, :Integer, W, X*M, M, O, true, A)
end
@generated function vload(ptr::Ptr{T}, i::LazyMulAdd{M,O,MM{W,X,StaticInt{I}}}, m::Mask{W}, ::Val{A}) where {A, W, M, O, X, T <: NativeTypes, I}
    vload_quote(T, Int, :StaticInt, W, X*M, M, O + I*M, true, A)
end

@inline function _vload_scalar(ptr::Ptr{Bit}, i::Integer, ::Val{A}) where {A}
    d = i >> 3; r = i & 7;
    u = vload(Base.unsafe_convert(Ptr{UInt8}, ptr), d, Val{A}())
    (u >> r) % Bool
end
@inline vload(ptr::Ptr{Bit}, i::IntegerTypesHW, ::Val{A}) where {A} = _vload_scalar(ptr, i, Val{A}())
# avoid ambiguities
@inline vload(ptr::Ptr{Bit}, ::StaticInt{N}, ::Val{A}) where {A,N} = _vload_scalar(ptr, StaticInt{N}(), Val{A}())

@inline vload(ptr::Union{Ptr,AbstractStridedPointer}) = vload(ptr, Val{false}())
@inline vloada(ptr::Union{Ptr,AbstractStridedPointer}) = vload(ptr, Val{true}())
@inline vload(ptr::Union{Ptr,AbstractStridedPointer}, i::Union{Number,Tuple,Unroll}) = vload(ptr, i, Val{false}())
@inline vloada(ptr::Union{Ptr,AbstractStridedPointer}, i::Union{Number,Tuple,Unroll}) = vload(ptr, i, Val{true}())
@inline vload(ptr::Union{Ptr,AbstractStridedPointer}, i::Union{Number,Tuple,Unroll}, m::Mask) = vload(ptr, i, m, Val{false}())
@inline vloada(ptr::Union{Ptr,AbstractStridedPointer}, i::Union{Number,Tuple,Unroll}, m::Mask) = vload(ptr, i, m, Val{true}())
@inline vload(ptr::Union{Ptr,AbstractStridedPointer}, i::Union{Number,Tuple,Unroll}, b::Bool) = vload(ptr, i, b, Val{false}())
@inline vloada(ptr::Union{Ptr,AbstractStridedPointer}, i::Union{Number,Tuple,Unroll}, b::Bool) = vload(ptr, i, b, Val{true}())

@inline function vload(ptr::Ptr{T}, i::Number, b::Bool, ::Val{A}) where {T,A}
    if b
        vload(ptr, i, Val{A}())
    else
        zero(T)
    end
end
@inline vwidth_from_ind(i::Tuple) = vwidth_from_ind(i, StaticInt(1), StaticInt(0))
@inline vwidth_from_ind(i::Tuple{}, ::StaticInt{W}, ::StaticInt{U}) where {W,U} = (StaticInt{W}(), StaticInt{U}())
@inline function vwidth_from_ind(i::Tuple{<:AbstractSIMDVector{W},Vararg}, ::Union{StaticInt{1},StaticInt{W}}, ::StaticInt{U}) where {W,U}
    vwidth_from_ind(Base.tail(i), StaticInt{W}(), StaticInt{W}(U))
end
@inline function vwidth_from_ind(
    i::Tuple{<:VecUnroll{U,W},Vararg}, ::Union{StaticInt{1},StaticInt{W}}, ::Union{StaticInt{0},StaticInt{U}}
) where {U,W}
    vwidth_from_ind(Base.tail(i), StaticInt{W}(), StaticInt{W}(U))
end
@inline zero_init(::Type{T}, ::StaticInt{1}, ::StaticInt{0}) where {T} = zero(T)
@inline zero_init(::Type{T}, ::StaticInt{W}, ::StaticInt{0}) where {W,T} = vzero(Val(W), T)
@inline zero_init(::Type{T}, ::StaticInt{W}, ::StaticInt{U}) where {W,U,T} = zero(VecUnroll{U,W,T,Vec{W,T}})

@inline zero_init(::Type{T}, ::Tuple{StaticInt{W},StaticInt{U}}) where {W,U,T} = zero_init(T, StaticInt{W}(), StaticInt{U}())

@inline function vload(ptr::Ptr{T}, i::Tuple, b::Bool, ::Val{A}) where {T,A}
    if b
        vload(ptr, i, Val{A}())
    else
        zero_init(T, vwidth_from_ind(i))
    end
end
@generated function zero_vecunroll(::Val{N}, ::Val{W}, ::Type{T}) where {N,W,T}
    Expr(:block, Expr(:meta, :inline), :(zero(VecUnroll{$(N-1),$W,$T,Vec{$W,$T}})))
end
@inline function vload(ptr::Ptr{T}, i::Unroll{AU,F,N,AV,W,M,I}, b::Bool, ::Val{A}) where {T,AU,F,N,AV,W,M,I,A}
    if b
        vload(ptr, i, Val{A}())
    else
        zero_vecunroll(Val{N}(), Val{W}(), T)
        # VecUnroll(ntuple(@inline(_ -> vzero(Val{W}(), T)), Val{N}()))
    end
end

function vstore_quote(
    ::Type{T}, ::Type{I}, ind_type::Symbol, W::Int, X::Int, M::Int, O::Int, mask::Bool, align::Bool, noalias::Bool, nontemporal::Bool
) where {T <: NativeTypes, I <: Integer}
    T_sym = JULIA_TYPES[T]
    I_sym = JULIA_TYPES[I]
    vstore_quote(T_sym, I_sym, ind_type, W, X, M, O, mask, align, noalias, nontemporal)
end
function vstore_quote(
    T_sym::Symbol, I_sym::Symbol, ind_type::Symbol, W::Int, X::Int, M::Int, O::Int, mask::Bool, align::Bool, noalias::Bool, nontemporal::Bool
)
    if W > 1
        vstore_quote(T_sym, I_sym, ind_type, W, X, M, O, mask, align, noalias, nontemporal, :(_Vec{$W,$T_sym}))
    else
        vstore_quote(T_sym, I_sym, ind_type, W, X, M, O, mask, align, noalias, nontemporal, T_sym)
    end
end
function vstore_quote(
    T_sym::Symbol, I_sym::Symbol, ind_type::Symbol, W::Int, X::Int, M::Int, O::Int, mask::Bool, align::Bool, noalias::Bool, nontemporal::Bool, jtyp
)
    sizeof_T = JULIA_TYPE_SIZE[T_sym]
    sizeof_I = JULIA_TYPE_SIZE[I_sym]
    ibits = 8sizeof_I
    if W > 1 && ind_type !== :Vec
        X, Xr = divrem(X, sizeof_T)
        @assert iszero(Xr)
    end
    instrs, i = offset_ptr(T_sym, ind_type, '2', ibits, W, X, M, O, false)

    grv = gep_returns_vector(W, X, M, ind_type)

    align != nontemporal # should I do this?
    alignment = (align & (!grv)) ? _get_alignment(W, T_sym) : _get_alignment(0, T_sym)

    decl = noalias ? SCOPE_METADATA * STORE_TBAA : STORE_TBAA
    metadata = noalias ? STORE_SCOPE_FLAGS * STORE_TBAA_FLAGS : STORE_TBAA_FLAGS
    dynamic_index = !(iszero(M) || ind_type === :StaticInt)

    typ = LLVM_TYPES_SYM[T_sym]
    lret = vtyp = vtype(W, typ)
    mask && truncate_mask!(instrs, '2' + dynamic_index, W, 0)
    if grv
        storeinstr = "void @llvm.masked.scatter." * suffix(W, T_sym) * '.' * ptr_suffix(W, T_sym)
        decl *= "declare $storeinstr($vtyp, <$W x $typ*>, i32, <$W x i1>)"
        m = mask ? m = "%mask.0" : llvmconst(W, "i1 1")
        push!(instrs, "call $storeinstr($vtyp %1, <$W x $typ*> %ptr.$(i-1), i32 $alignment, <$W x i1> $m)" * metadata)
        # push!(instrs, "call $storeinstr($vtyp %1, <$W x $typ*> %ptr.$(i-1), i32 $alignment, <$W x i1> $m)")
    elseif mask
        suff = suffix(W, T_sym)
        storeinstr = "void @llvm.masked.store." * suff * ".p0" * suff
        decl *= "declare $storeinstr($vtyp, $vtyp*, i32, <$W x i1>)"
        push!(instrs, "call $storeinstr($vtyp %1, $vtyp* %ptr.$(i-1), i32 $alignment, <$W x i1> %mask.0)" * metadata)
    elseif nontemporal
        push!(instrs, "store $vtyp %1, $vtyp* %ptr.$(i-1), align $alignment, !nontemporal !{i32 1}" * metadata)
    else
        push!(instrs, "store $vtyp %1, $vtyp* %ptr.$(i-1), align $alignment" * metadata)
    end
    push!(instrs, "ret void")
    ret = :Cvoid; lret = "void"
    ptrtyp = Expr(:curly, :Ptr, T_sym)
    args = if W > 1
        Expr(:curly, :Tuple, ptrtyp, Expr(:curly, :NTuple, W, Expr(:curly, :VecElement, T_sym)))
    else
        Expr(:curly, :Tuple, ptrtyp, T_sym)
    end
    largs = String[JULIAPOINTERTYPE, vtyp]
    arg_syms = Union{Symbol,Expr}[:ptr, Expr(:call, :data, :v)]
    if dynamic_index
        push!(arg_syms, :(data(i)))
        if ind_type === :Integer
            push!(args.args, I_sym)
            push!(largs, "i$(ibits)")
        else
            push!(args.args, :(_Vec{$W,$I_sym}))
            push!(largs, "<$W x i$(ibits)>")
        end
    end
    if mask
        push!(arg_syms, :(data(m)))
        push!(args.args, mask_type(W))
        push!(largs, "i$(max(8,nextpow2(W)))")
    end
    llvmcall_expr(decl, join(instrs, "\n"), ret, args, lret, largs, arg_syms)
end

@generated function vstore!(
    ptr::Ptr{T}, v::T, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T <: NativeTypes, A, S, NT}
    vstore_quote(T, Int, :StaticInt, 1, 1, 0, 0, false, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::T, i::I, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T <: NativeTypes, I <: IntegerTypes, A, S, NT}
    vstore_quote(T, I, :Integer, 1, 1, 1, 0, false, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::T, ::StaticInt{N}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T <: NativeTypes, N, A, S, NT}
    vstore_quote(T, Int, :StaticInt, 1, 1, 0, N, false, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, i::Vec{W,I}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T <: NativeTypes, I <: IntegerTypes, A, S, NT}
    vstore_quote(T, I, :Vec, W, 1, 1, 0, false, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, i::MM{W,X,I}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, X, T <: NativeTypes, I <: IntegerTypes, A, S, NT}
    vstore_quote(T, I, :Integer, W, X, 1, 0, false, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, i::I, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T <: NativeTypes, I <: IntegerTypes, A, S, NT}
    vstore_quote(T, I, :Integer, W, sizeof(T), 1, 0, false, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T <: NativeTypes, A, S, NT}
    vstore_quote(T, Int, :StaticInt, W, sizeof(T), 0, 0, false, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, ::StaticInt{N}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T <: NativeTypes, N, A, S, NT}
    vstore_quote(T, Int, :StaticInt, W, sizeof(T), 0, N, false, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, ::MM{W,X,StaticInt{N}}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, X, T <: NativeTypes, N, A, S, NT}
    vstore_quote(T, Int, :StaticInt, W, X, 0, N, false, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::T, i::LazyMulAdd{M,O,I}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T <: NativeTypes, M, O, I <: IntegerTypes, A, S, NT}
    vstore_quote(T, I, :Integer, 1, 1, M, O, false, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, i::LazyMulAdd{M,O,Vec{W,I}}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T <: NativeTypes, M, O, I <: IntegerTypes, A, S, NT}
    vstore_quote(T, I, :Vec, W, 1, M, O, false, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, i::LazyMulAdd{M,O,MM{W,X,I}}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, X, T <: NativeTypes, M, O, I <: IntegerTypes, A, S, NT}
    vstore_quote(T, I, :Integer, W, X*M, M, O, false, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, i::LazyMulAdd{M,O,MM{W,X,StaticInt{I}}}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, X, T <: NativeTypes, M, O, I, A, S, NT}
    vstore_quote(T, Int, :StaticInt, W, X*M, M, O + M*I, false, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, i::LazyMulAdd{M,O,I}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T <: NativeTypes, M, O, I <: IntegerTypes, A, S, NT}
    vstore_quote(T, I, :Integer, W, sizeof(T), M, O, false, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, i::Vec{W,I}, m::Mask{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T <: NativeTypes, I <: IntegerTypes, A, S, NT}
    vstore_quote(T, I, :Vec, W, 1, 1, 0, true, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, i::MM{W,X,I}, m::Mask{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, X, T <: NativeTypes, I <: IntegerTypes, A, S, NT}
    vstore_quote(T, I, :Integer, W, X, 1, 0, true, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, i::I, m::Mask{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T <: NativeTypes, I <: IntegerTypes, A, S, NT}
    vstore_quote(T, I, :Integer, W, sizeof(T), 1, 0, true, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, ::StaticInt{N}, m::Mask{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T <: NativeTypes, N, A, S, NT}
    vstore_quote(T, Int, :StaticInt, W, sizeof(T), 0, N, true, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, ::MM{W,X,StaticInt{N}}, m::Mask{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, X, T <: NativeTypes, N, A, S, NT}
    vstore_quote(T, Int, :StaticInt, W, X, 0, N, true, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, i::LazyMulAdd{M,O,Vec{W,I}}, m::Mask{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T <: NativeTypes, M, O, I <: IntegerTypes, A, S, NT}
    vstore_quote(T, I, :Vec, W, 1, M, O, true, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, i::LazyMulAdd{M,O,MM{W,X,I}}, m::Mask{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, X, T <: NativeTypes, M, O, I <: IntegerTypes, A, S, NT}
    vstore_quote(T, I, :Integer, W, X*M, M, O, true, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, i::LazyMulAdd{M,O,MM{W,X,StaticInt{I}}}, m::Mask{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, X, T <: NativeTypes, M, O, I, A, S, NT}
    vstore_quote(T, Int, :StaticInt, W, X*M, M, O + M*I, true, A, S, NT)
end
@generated function vstore!(
    ptr::Ptr{T}, v::Vec{W,T}, i::LazyMulAdd{M,O,I}, m::Mask{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T <: NativeTypes, M, O, I <: IntegerTypes, A, S, NT}
    vstore_quote(T, I, :Integer, W, sizeof(T), M, O, true, A, S, NT)
end

# broadcasting scalar stores
@inline function vstore!(
    ptr::Ptr{T}, v::Base.HWReal, i::VectorIndex{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T, A, S, NT}
    vstore!(ptr, convert(Vec{W,T}, v), i, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::Ptr{T}, v::Base.HWReal, i::VectorIndex{W}, m::Mask, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T, A, S, NT}
    vstore!(ptr, convert(Vec{W,T}, v), i, m, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::Ptr{T}, v::Base.HWReal, i::MM{1}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T, A, S, NT}
    vstore!(ptr, convert(T, v), i.i, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::Ptr{T}, v::Base.HWReal, i::MM{1}, m::Mask{1}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T, A, S, NT}
    Bool(m) && vstore!(ptr, convert(T, v), i.i, Val{A}(), Val{S}(), Val{NT}())
    nothing
end
@inline function vstore!(
    ptr::Ptr{T}, v::Base.HWReal, i::LazyMulAdd{M,O,MM{1,X,I}}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T, A, S, NT, M,O,X,I}
    vstore!(ptr, convert(T, v), _materialize(i).i, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::Ptr{T}, v::Base.HWReal, i::LazyMulAdd{M,O,MM{1,X,I}}, m::Mask{1}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T, A, S, NT, M,O,X,I}
    Bool(m) && vstore!(ptr, convert(T, v), _materialize(i).i, Val{A}(), Val{S}(), Val{NT}())
    nothing
end

# @inline function vstore!(ptr::Ptr{Bit}, m::Mask{W,U}, i, ::Val{A}, ::Val{S}, ::Val{NT}) where {W,U,A,S,NT}
#     @assert W == 8sizeof(U)
#     vstore!(Base.unsafe_convert(Ptr{U}, ptr), data(m), i, Val{A}(), Val{S}(), Val{NT}())
# end

# BitArray stores
@inline function vstore!(ptr::Ptr{Bit}, v::Mask{W,U}, ::Val{A}, ::Val{S}, ::Val{NT}) where {W, U, A, S, NT}
    vstore!(Base.unsafe_convert(Ptr{U}, ptr), data(v), Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::Ptr{Bit}, v::Mask{W,U}, i::VectorIndex{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, U, A, S, NT}
    vstore!(Base.unsafe_convert(Ptr{U}, ptr), data(v), data(i) >> 3, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::Ptr{Bit}, v::Mask{W,U}, i::VectorIndex{W}, m::Mask, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, U, A, S, NT}
    ishift = data(i) >> 3
    u = bitselect(data(m), vload(Base.unsafe_convert(Ptr{U}, ptr), ishift, Val{A}()), data(v))
    vstore!(Base.unsafe_convert(Ptr{U}, ptr), u, ishift, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(f::F, ptr::Ptr{Bit}, v::Mask{W,U}, ::Val{A}, ::Val{S}, ::Val{NT}) where {W, U, A, S, NT, F<:Function}
    vstore!(f, Base.unsafe_convert(Ptr{U}, ptr), data(v), Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    f::F, ptr::Ptr{Bit}, v::Mask{W,U}, i::VectorIndex{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, U, A, S, NT, F<:Function}
    vstore!(f, Base.unsafe_convert(Ptr{U}, ptr), data(v), data(i) >> 3, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    f::F, ptr::Ptr{Bit}, v::Mask{W,U}, i::VectorIndex{W}, m::Mask, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, U, A, S, NT, F<:Function}
    ishift = data(i) >> 3
    u = bitselect(data(m), vload(Base.unsafe_convert(Ptr{U}, ptr), ishift, Val{A}()), data(v))
    vstore!(f, Base.unsafe_convert(Ptr{U}, ptr), u, ishift, Val{A}(), Val{S}(), Val{NT}())
end

@inline function vstore!(ptr::Ptr{T}, v, ::Val{A}, ::Val{S}, ::Val{NT}) where {T, A, S, NT}
    vstore!(ptr, convert(T, v), Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::Ptr{T}, v, i, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T, A, S, NT}
    vstore!(ptr, convert(T, v), i, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::Ptr{T}, v, i::VectorIndex{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W,T, A, S, NT}
    vstore!(ptr, convert(Vec{W,T}, v), i, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    ptr::Ptr{T}, v, i, m::Mask{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T, A, S, NT}
    vstore!(ptr, convert(Vec{W,T}, v), i, m, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(f::F, ptr::Ptr{T}, v, ::Val{A}, ::Val{S}, ::Val{NT}) where {T, A, S, NT, F<:Function}
    vstore!(f, ptr, convert(T, v), Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    f::F, ptr::Ptr{T}, v, i, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T, A, S, NT, F<:Function}
    vstore!(f, ptr, convert(T, v), i, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    f::F, ptr::Ptr{T}, v, i::VectorIndex{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T, A, S, NT, F<:Function}
    vstore!(f, ptr, convert(Vec{W,T}, v), i, Val{A}(), Val{S}(), Val{NT}())
end
@inline function vstore!(
    f::F, ptr::Ptr{T}, v, i, m::Mask{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {W, T, A, S, NT, F<:Function}
    vstore!(f, ptr, convert(Vec{W,T}, v), i, m, Val{A}(), Val{S}(), Val{NT}())
end

for (store,align,alias,nontemporal) ∈ [
    (:vstore!,false,false,false),
    (:vstorea!,true,false,false),
    (:vstorent!,true,false,true),
    (:vnoaliasstore!,false,true,false),
    (:vnoaliasstorea!,true,true,false),
    (:vnoaliasstorent!,true,true,true)
]
    @eval begin
        @inline function $store(ptr::Union{Ptr,AbstractStridedPointer}, v::Number)
            vstore!(ptr, v, Val{$align}(), Val{$alias}(), Val{$nontemporal}())
        end
        @inline function $store(ptr::Union{Ptr,AbstractStridedPointer}, v::Number, i::Union{Number,Tuple,Unroll})
            vstore!(ptr, v, i, Val{$align}(), Val{$alias}(), Val{$nontemporal}())
        end
        @inline function $store(ptr::Union{Ptr,AbstractStridedPointer}, v::Number, i::Union{Number,Tuple,Unroll}, m::Mask)
            vstore!(ptr, v, i, m, Val{$align}(), Val{$alias}(), Val{$nontemporal}())
        end
        @inline function $store(ptr::Union{Ptr,AbstractStridedPointer}, v::Number, i::Union{Number,Tuple,Unroll}, b::Bool)
            b && vstore!(ptr, v, i, Val{$align}(), Val{$alias}(), Val{$nontemporal}())
        end

        @inline function $store(f::F, ptr::Union{Ptr,AbstractStridedPointer}, v::Number) where {F<:Function}
            vstore!(f, ptr, v, Val{$align}(), Val{$alias}(), Val{$nontemporal}())
        end
        @inline function $store(f::F, ptr::Union{Ptr,AbstractStridedPointer}, v::Number, i::Union{Number,Tuple,Unroll}) where {F<:Function}
            vstore!(f, ptr, v, i, Val{$align}(), Val{$alias}(), Val{$nontemporal}())
        end
        @inline function $store(f::F, ptr::Union{Ptr,AbstractStridedPointer}, v::Number, i::Union{Number,Tuple,Unroll}, m::Mask) where {F<:Function}
            vstore!(f, ptr, v, i, m, Val{$align}(), Val{$alias}(), Val{$nontemporal}())
        end
        @inline function $store(f::F, ptr::Union{Ptr,AbstractStridedPointer}, v::Number, i::Union{Number,Tuple,Unroll}, b::Bool) where {F<:Function}
            b && vstore!(f, ptr, v, i, Val{$align}(), Val{$alias}(), Val{$nontemporal}())
        end
    end
end

# unroll
@inline Base.Broadcast.broadcastable(u::Unroll) = (u,)

"""
Returns a vector of expressions for a set of unrolled indices.

"""
function unrolled_indicies(D::Int, AU::Int, F::Int, N::Int, AV::Int, W::Int)
    baseind = Expr(:tuple)
    for d in 1:D
        i = Expr(:call, :Zero)
        if d == AV && W > 1
            i = Expr(:call, Expr(:curly, :MM, W), i)
        end
        push!(baseind.args, i)
    end
    WF = AU == AV ? W : 1
    inds = Vector{Expr}(undef, N)
    inds[1] = baseind
    for n in 1:N-1
        ind = copy(baseind)
        i = Expr(:call, Expr(:curly, :StaticInt, n*F*WF))
        if AU == AV && W > 1
            i = Expr(:call, Expr(:curly, :MM, W), i)
        end
        ind.args[AU] = i
        inds[n+1] = ind
    end
    inds
end

function vload_unroll_quote(D::Int, AU::Int, F::Int, N::Int, AV::Int, W::Int, M::UInt, mask::Bool, align::Bool)
    t = Expr(:tuple)
    inds = unrolled_indicies(D, AU, F, N, AV, W)
    # TODO: Consider doing some alignment checks before accepting user's `align`?
    alignval = Expr(:call, Expr(:curly, :Val, align))
    for n in 1:N
        l = Expr(:call, :vload, :gptr, inds[n])
        (mask && (M % Bool)) && push!(l.args, :m)
        push!(l.args, alignval)
        M >>= 1
        push!(t.args, l)
    end
    quote
        $(Expr(:meta, :inline))
        gptr = gesp(ptr, u.i)
        VecUnroll($t)
    end
end

@generated function vload(ptr::AbstractStridedPointer{T,D}, u::Unroll{AU,F,N,AV,W,M,I}, ::Val{A}) where {A,AU,F,N,AV,W,M,I,T,D}
    vload_unroll_quote(D, AU, F, N, AV, W, M, false, A)
end
@generated function vload(ptr::AbstractStridedPointer{T,D}, u::Unroll{AU,F,N,AV,W,M,I}, m::Mask{W}, ::Val{A}) where {A,AU,F,N,AV,W,M,I,T,D}
    vload_unroll_quote(D, AU, F, N, AV, W, M, true, A)
end

function vstore_unroll_quote(
    D::Int, AU::Int, F::Int, N::Int, AV::Int, W::Int, M::UInt, mask::Bool, align::Bool, noalias::Bool, nontemporal::Bool
)
    t = Expr(:tuple)
    inds = unrolled_indicies(D, AU, F, N, AV, W)
    q = quote
        $(Expr(:meta, :inline))
        gptr = gesp(ptr, u.i)
        t = data(v)
    end
    alignval = Expr(:call, Expr(:curly, :Val, align))
    noaliasval = Expr(:call, Expr(:curly, :Val, noalias))
    nontemporalval = Expr(:call, Expr(:curly, :Val, nontemporal))
    for n in 1:N
        l = Expr(:call, :vstore!, :gptr, Expr(:ref, :t, n), inds[n])
        (mask && (M % Bool)) && push!(l.args, :m)
        push!(l.args, alignval, noaliasval, nontemporalval)
        M >>= 1
        push!(q.args, l)
    end
    q
end
@generated function vstore!(
    ptr::AbstractStridedPointer{T,D}, v::VecUnroll{Nm1,W}, u::Unroll{AU,F,N,AV,W,M,I}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {AU,F,N,AV,W,M,I,T,D,Nm1,S,A,NT}
    @assert Nm1+1 == N
    vstore_unroll_quote(D, AU, F, N, AV, W, M, false, A, S, NT)
end
@generated function vstore!(
    ptr::AbstractStridedPointer{T,D}, v::VecUnroll{Nm1,W}, u::Unroll{AU,F,N,AV,W,M,I}, m::Mask{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {AU,F,N,AV,W,M,I,T,D,Nm1,S,A,NT}
    @assert Nm1+1 == N
    vstore_unroll_quote(D, AU, F, N, AV, W, M, true, A, S, NT)
end
function vstore_unroll_i_quote(Nm1, Wsplit, W, A, S, NT, mask::Bool)
    N = Nm1 + 1
    @assert N*Wsplit == W
    q = Expr(:block, Expr(:meta, :inline), :(vt = data(v)), :(im = _materialize(i)))
    if mask
        let U = mask_type_symbol(Wsplit)
            push!(q.args, :(mt = data(vconvert(VecUnroll{$Nm1,$Wsplit,Bit,Mask{$Wsplit,$U}}, m))))
        end
    end
    j = 0
    alignval = Expr(:call, Expr(:curly, :Val, A))
    aliasval = Expr(:call, Expr(:curly, :Val, S))
    notmpval = Expr(:call, Expr(:curly, :Val, NT))
    for n ∈ 1:N
        shufflemask = Expr(:tuple)
        for w ∈ 1:Wsplit
            push!(shufflemask.args, j)
            j += 1
        end
        ex = :(vstore!(ptr, vt[$n], shufflevector(im, Val{$shufflemask}())))
        mask && push!(ex.args, Expr(:ref, :mt, n))
        push!(ex.args, alignval)
        push!(ex.args, aliasval)
        push!(ex.args, notmpval)
        push!(q.args, ex)
    end
    q
end
@generated function vstore!(
    ptr::Ptr{T}, v::VecUnroll{Nm1,Wsplit}, i::VectorIndex{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T,Nm1,Wsplit,W,A,S,NT}
    vstore_unroll_i_quote(Nm1, Wsplit, W, A, S, NT, false)
end
@generated function vstore!(
    ptr::Ptr{T}, v::VecUnroll{Nm1,Wsplit}, i::VectorIndex{W}, m::Mask{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T,Nm1,Wsplit,W,A,S,NT}
    vstore_unroll_i_quote(Nm1, Wsplit, W, A, S, NT, true)
end
function vstorebit_unroll_i_quote(Nm1, Wsplit, W, A, S, NT, mask::Bool)
    N = Nm1 + 1
    @assert N*Wsplit == W
    @assert W == 8
    # q = Expr(:block, Expr(:meta, :inline), :(vt = data(v)), :(im = _materialize(i)), :(u = 0x00))
    q = Expr(:block, Expr(:meta, :inline), :(vt = data(v)), :(u = 0x00))
    j = 0
    while true
        push!(q.args, :(u |= data($(Expr(:ref, :vt, (N-j))))))
        j += 1
        j == N && break
        push!(q.args, :(u <<= $Wsplit))
    end
    alignval = Expr(:call, Expr(:curly, :Val, A))
    aliasval = Expr(:call, Expr(:curly, :Val, S))
    notmpval = Expr(:call, Expr(:curly, :Val, NT))
    mask && push!(q.args, :(u = bitselect(data(m), vload(Base.unsafe_convert(Ptr{$(mask_type_symbol(W))}, ptr), (i.i >> 3), Val{$A}()), u)))
    call = Expr(:call, :vstore!, :(reinterpret(Ptr{UInt8}, ptr)), :u, :(i.i >> 3))
    push!(call.args, alignval, aliasval, notmpval)
    push!(q.args, call)
    q
end
@generated function vstore!(
    ptr::Ptr{Bit}, v::VecUnroll{Nm1,Wsplit,Bit,Mask{Wsplit,UInt8}}, i::MM{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {Nm1,Wsplit,W,A,S,NT}
    # 1 + 1
    vstorebit_unroll_i_quote(Nm1, Wsplit, W, A, S, NT, false)
end
@generated function vstore!(
    ptr::Ptr{Bit}, v::VecUnroll{Nm1,Wsplit,Bit,Mask{Wsplit,UInt8}}, i::MM{W}, m::Mask{W}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {Nm1,Wsplit,W,A,S,NT}
    vstorebit_unroll_i_quote(Nm1, Wsplit, W, A, S, NT, true)
end

# @inline vstore!(::typeof(identity), ptr, v, u) = vstore!(ptr, v, u)
# @inline vstore!(::typeof(identity), ptr, v, u, m) = vstore!(ptr, v, u, m)
# @inline vnoaliasstore!(::typeof(identity), ptr, v, u) = vnoaliasstore!(ptr, v, u)
# @inline vnoaliasstore!(::typeof(identity), ptr, v, u, m) = vnoaliasstore!(ptr, v, u, m)

# If `::Function` vectorization is masked, then it must not be reduced by `::Function`.
@generated function vstore!(
    ::Function, ptr::AbstractStridedPointer{T,D,C}, vu::VecUnroll{U,W}, u::Unroll{AU,F,N,AV,W,M,I}, m, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T,D,C,U,AU,F,N,W,M,I,AV,A,S,NT}
    @assert N == U + 1
    # mask means it isn't vectorized
    @assert AV > 0 "AV ≤ 0, but masking what, exactly?"
    Expr(:block, Expr(:meta, :inline), :(vstore!(ptr, vu, u, m, Val{$A}(), Val{$S}(), Val{$NT}())))
end

function transposeshuffle(split, W, offset::Bool)
    tup = Expr(:tuple)
    w = 0
    S = 1 << split
    i = offset ? S : 0
    while w < W
        for s ∈ 0:S-1
            push!(tup.args, w + s + i)
        end
        for s ∈ 0:S-1
            # push!(tup.args, w + W + s)
            push!(tup.args, w + W + s + i)
        end
        w += 2S
    end
    Expr(:call, Expr(:curly, :Val, tup))
end

function horizonal_reduce_store_expr(W, Ntotal, (C,D,AU,F), op = :+, reduct = :vsum, noalias::Bool = false)
    N = ((C == AU) && isone(F)) ? prevpow2(Ntotal) : 0
    q = Expr(:block, Expr(:meta, :inline), :(v = data(vu)))
    store = noalias ? :vnoaliasstore! : :vstore!
    @assert ispow2(W)
    if N > 1
        if N < Ntotal
            push!(q.args, :(gptr = gesp(ptr, u.i)))
            push!(q.args, :(bptr = pointer(gptr)))
        else
            push!(q.args, :(bptr = gep(ptr, u.i)))
        end
        extractblock = Expr(:block)
        vectors = [Symbol(:v_, n) for n ∈ 0:N-1]
        for n ∈ 1:N
            push!(extractblock.args, Expr(:(=), vectors[n], Expr(:ref, :v, n)))
        end
        push!(q.args, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), extractblock))
        ncomp = 0
        minWN = min(W,N)
        while ncomp < N
            Nt = minWN;
            Wt = W
            splits = 0
            while Nt > 1
                Nt >>>= 1
                shuffle0 = transposeshuffle(splits, Wt, false)
                shuffle1 = transposeshuffle(splits, Wt, true)
                splits += 1
                for nh ∈ 1:Nt
                    n1 = 2nh
                    n0 = n1 - 1
                    v0 = vectors[n0 + ncomp]; v1 = vectors[n1 + ncomp]; vh = vectors[nh + ncomp];
                    # combine n0 and n1
                    push!(q.args, Expr(
                        :(=), vh, Expr(
                            :call, op,
                            Expr(:call, :shufflevector, v0, v1, shuffle0),
                            Expr(:call, :shufflevector, v0, v1, shuffle1))
                    ))
                end
            end
            # v0 is now the only vector
            v0 = vectors[ncomp + 1]
            while Wt > minWN
                Wh = Wt >>> 1
                v0new = Symbol(v0, Wt)
                push!(q.args, Expr(
                    :(=), v0new, Expr(
                        :call, op,
                        Expr(:call, :shufflevector, v0, Expr(:call, Expr(:curly, :Val, Expr(:tuple, [w for w ∈ 0:Wh-1]...)))),
                        Expr(:call, :shufflevector, v0, Expr(:call, Expr(:curly, :Val, Expr(:tuple, [w for w ∈ Wh:Wt-1]...)))))
                )
                      )
                v0 = v0new
                Wt = Wh
            end
            push!(q.args, Expr(:call, store, :bptr, v0))
            ncomp += minWN
        end
    else
        push!(q.args, :(gptr = gesp(ptr, u.i)))
    end
    if N < Ntotal
        zeroexpr = Expr(:call, Expr(:curly, :StaticInt, 0))
        ind = Expr(:tuple); foreach(_ -> push!(ind.args, zeroexpr), 1:D)
        for n ∈ N+1:Ntotal
            (n > N+1) && (ind = copy(ind)) # copy to avoid overwriting old
            ind.args[AU] = Expr(:call, Expr(:curly, :StaticInt, F*(n-1)))
            scalar = Expr(:call, reduct, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), Expr(:ref, :v, n)))
            push!(q.args, Expr(:call, store, :gptr, scalar, ind))
        end
    end
    q
end
@generated function vstore!(
    ::G, ptr::AbstractStridedPointer{T,D,C}, vu::VecUnroll{U,W}, u::Unroll{AU,F,N,AV,W,M,I}, ::Val{A}, ::Val{S}, ::Val{NT}
) where {T,D,C,U,AU,F,N,W,M,I,G<:Function,AV,A,S,NT}
    @assert N == U + 1
    if G === typeof(identity) || AV > 0
        return Expr(:block, Expr(:meta, :inline), :(vstore!(ptr, vu, u, Val{$A}(), Val{$S}(), Val{$NT}())))
    elseif G === typeof(vsum)
        op = :+; reduct = :vsum
    elseif G === typeof(vprod)
        op = :*; reduct = :vprod
    else
        throw("Function $f not recognized.")
    end
    horizonal_reduce_store_expr(W, N, (C,D,AU,F), op, reduct, S)
end

function lazymulunroll_load_quote(M,O,N,mask,align)
    t = Expr(:tuple)
    alignval = Expr(:call, Expr(:curly, :Val, align))
    for n in 1:N+1
        ind = if (M != 1) | (O != 0)
            :(LazyMulAdd{$M,$O}(u[$n]))
        else
            Expr(:ref, :u, n)
        end
        call = if mask
            Expr(:call, :vload, :ptr, ind, Expr(:ref, :mt, n), alignval)
        else
            Expr(:call, :vload, :ptr, ind, alignval)
        end
        push!(t.args, call)
    end
    q = Expr(:block, Expr(:meta, :inline), :(u = data(um)))
    mask && push!(q.args, :(mt = data(m)))
    push!(q.args, Expr(:call, :VecUnroll, t))
    q
end
@generated function vload(ptr::Ptr{T}, um::VecUnroll{N,W,I,V}, ::Val{A}) where {T,N,W,I,V,A}
    lazymulunroll_load_quote(1,0,N,false,A)
end
@generated function vload(ptr::Ptr{T}, um::VecUnroll{N,W,I,V}, m::VecUnroll{N,W,Bit,Mask{W,U}}, ::Val{A}) where {T,N,W,I,V,A,U}
    lazymulunroll_load_quote(1,0,N,true,A)
end
@inline function vload(ptr::Ptr{T}, um::VecUnroll{N,W,I,V}, m::Mask, ::Val{A}) where {T,N,W,I,V,A}
    vload(ptr, um, VecUnroll(splitvectortotuple(StaticInt{N}() + One(), StaticInt{W}(), m)), Val{A}())
end
@generated function vload(ptr::Ptr{T}, um::LazyMulAdd{M,O,VecUnroll{N,W,I,V}}, ::Val{A}) where {T,M,O,N,W,I,V,A}
    lazymulunroll_load_quote(M,O,N,false,A)
end
@generated function vload(ptr::Ptr{T}, um::LazyMulAdd{M,O,VecUnroll{N,W,I,V}}, m::VecUnroll{N,W,Bit,Mask{W,U}}, ::Val{A}) where {T,M,O,N,W,I,V,A,U}
    lazymulunroll_load_quote(M,O,N,true,A)
end
@inline function vload(ptr::Ptr{T}, um::LazyMulAdd{M,O,VecUnroll{N,W,I,V}}, m::Mask, ::Val{A}) where {T,M,O,N,W,I,V,A}
    vload(ptr, um, VecUnroll(splitvectortotuple(StaticInt{N}() + One(), StaticInt{W}(), m)), Val{A}())
end
function lazymulunroll_store_quote(M,O,N,mask,align,noalias,nontemporal)
    q = Expr(:block, Expr(:meta, :inline), :(u = um.data.data), :(v = vm.data.data))
    alignval = Expr(:call, Expr(:curly, :Val, align))
    noaliasval = Expr(:call, Expr(:curly, :Val, noalias))
    nontemporalval = Expr(:call, Expr(:curly, :Val, nontemporal))
    for n in 1:N+1
        push!(q.args, Expr(:call, :vstore!, :ptr, Expr(:ref, :v, n), :(LazyMulAdd{$M,$O}(u[$n])), alignval, noaliasval, nontemporalval))
    end
    q
end

@generated function prefetch(ptr::Ptr{Cvoid}, ::Val{L}, ::Val{R}) where {L, R}
    @assert L ∈ (0,1,2,3)
    @assert R ∈ (0, 1)
    decl = "declare void @llvm.prefetch(i8*, i32, i32, i32)"
    instrs = """
        %addr = inttoptr $JULIAPOINTERTYPE %0 to i8*
        call void @llvm.prefetch(i8* %addr, i32 $R, i32 $L, i32 1)
        ret void
    """
    llvmcall_expr(decl, instrs, Cvoid, :(Tuple{Ptr{Cvoid}}), "void", [JULIAPOINTERTYPE], [:ptr])
end
@inline prefetch(ptr::Ptr{T}, ::Val{L}, ::Val{R}) where {T,L,R} = prefetch(Base.unsafe_convert(Ptr{Cvoid}, ptr), Val{L}(), Val{R}())

@inline function prefetch(ptr::Union{AbstractStridedPointer,Ptr}, i, ::Val{Locality}, ::Val{ReadOrWrite}) where {Locality, ReadOrWrite}
    prefetch(gep(ptr, i), Val{Locality}(), Val{ReadOrWrite}())
end
@inline prefetch(ptr::Ptr) = prefetch(ptr, Val{3}(), Val{0}())
@inline prefetch(ptr::Ptr, ::Val{L}) where {L} = prefetch(ptr, Val{L}(), Val{0}())
@inline prefetch(ptr::Ptr, i) = prefetch(ptr, i, Val{3}(), Val{0}())
@inline prefetch(ptr::Ptr, i, ::Val{L}) where {L} = prefetch(ptr, i, Val{L}(), Val{0}())

@inline prefetch0(x, i) = prefetch(gep(stridedpointer(x), (data(i),)), Val{3}(), Val{0}())
@inline prefetch0(x, I::Tuple) = prefetch(gep(stridedpointer(x), data.(I)), Val{3}(), Val{0}())
@inline prefetch0(x, i, j) = prefetch(gep(stridedpointer(x), (data(i), data(j))), Val{3}(), Val{0}())
# @inline prefetch0(x, i, j, oi, oj) = prefetch(gep(stridedpointer(x), (data(i) + data(oi) - 1, data(j) + data(oj) - 1)), Val{3}(), Val{0}())
@inline prefetch1(x, i) = prefetch(gep(stridedpointer(x), (data(i),)), Val{2}(), Val{0}())
@inline prefetch1(x, i, j) = prefetch(gep(stridedpointer(x), (data(i), data(j))), Val{2}(), Val{0}())
# @inline prefetch1(x, i, j, oi, oj) = prefetch(gep(stridedpointer(x), (data(i) + data(oi) - 1, data(j) + data(oj) - 1)), Val{2}(), Val{0}())
@inline prefetch2(x, i) = prefetch(gep(stridedpointer(x), (data(i),)), Val{1}(), Val{0}())
@inline prefetch2(x, i, j) = prefetch(gep(stridedpointer(x), (data(i), data(j))), Val{1}(), Val{0}())
# @inline prefetch2(x, i, j, oi, oj) = prefetch(gep(stridedpointer(x), (data(i) + data(oi) - 1, data(j) + data(oj) - 1)), Val{1}(), Val{0}())

@generated function lifetime_start!(ptr::Ptr{T}, ::Val{L}) where {L,T}
    decl = "declare void @llvm.lifetime.start(i64, i8* nocapture)"
    instrs = "%ptr = inttoptr $JULIAPOINTERTYPE %0 to i8*\ncall void @llvm.lifetime.start(i64 $(L*sizeof(T)), i8* %ptr)\nret void"
    llvmcall_expr(decl, instrs, :Cvoid, :(Tuple{Ptr{$T}}), "void", [JULIAPOINTERTYPE], [:ptr])
end
@generated function lifetime_end!(ptr::Ptr{T}, ::Val{L}) where {L,T}
    decl = "declare void @llvm.lifetime.end(i64, i8* nocapture)"
    instrs = "%ptr = inttoptr $JULIAPOINTERTYPE %0 to i8*\ncall void @llvm.lifetime.end(i64 $(L*sizeof(T)), i8* %ptr)\nret void"
    llvmcall_expr(decl, instrs, :Cvoid, :(Tuple{Ptr{$T}}), "void", [JULIAPOINTERTYPE], [:ptr])
end

@inline lifetime_start!(ptr::Ptr) = lifetime_start!(ptr, Val{-1}())
@inline lifetime_end!(ptr::Ptr) = lifetime_end!(ptr, Val{-1}())
# Fallback is to do nothing. Intention is (e.g.) for PaddedMatrices/StackPointers.
@inline lifetime_start!(::Any) = nothing
@inline lifetime_end!(::Any) = nothing

@generated function compressstore!(ptr::Ptr{T}, v::Vec{W,T}, mask::Mask{W,U}) where {W,T <: NativeTypes, U<:Unsigned}
    @assert 8sizeof(U) >= W
    typ = LLVM_TYPES[T]
    vtyp = "<$W x $typ>"
    mtyp_input = LLVM_TYPES[U]
    mtyp_trunc = "i$W"
    instrs = String["%ptr = inttoptr $JULIAPOINTERTYPE %1 to $typ*"]
    truncate_mask!(instrs, '2', W, 0)
    decl = "declare void @llvm.masked.compressstore.$(suffix(W,T))($vtyp, $typ*, <$W x i1>)"
    push!(instrs, "call void @llvm.masked.compressstore.$(suffix(W,T))($vtyp %0, $typ* %ptr, <$W x i1> %mask.0)\nret void")
    llvmcall_expr(decl, join(instrs,"\n"), :Cvoid, :(Tuple{NTuple{$W,VecElement{$T}}, Ptr{$T}, $U}), "void", [vtyp, JULIAPOINTERTYPE, "i$(8sizeof(U))"], [:(data(v)), :ptr, :(data(mask))])
end

@generated function expandload(ptr::Ptr{T}, mask::Mask{W,U}) where {W, T <: NativeTypes, U<:Unsigned}
    @assert 8sizeof(U) >= W
    typ = LLVM_TYPES[T]
    vtyp = "<$W x $typ>"
    vptrtyp = "<$W x $typ*>"
    mtyp_input = LLVM_TYPES[U]
    mtyp_trunc = "i$W"
    instrs = String[]
    push!(instrs, "%ptr = inttoptr $JULIAPOINTERTYPE %0 to $typ*")
    if mtyp_input == mtyp_trunc
        push!(instrs, "%mask = bitcast $mtyp_input %1 to <$W x i1>")
    else
        push!(instrs, "%masktrunc = trunc $mtyp_input %1 to $mtyp_trunc")
        push!(instrs, "%mask = bitcast $mtyp_trunc %masktrunc to <$W x i1>")
    end
    decl = "declare $vtyp @llvm.masked.expandload.$(suffix(W,T))($typ*, <$W x i1>, $vtyp)"
    push!(instrs, "%res = call $vtyp @llvm.masked.expandload.$(suffix(W,T))($typ* %ptr, <$W x i1> %mask, $vtyp zeroinitializer)\nret $vtyp %res")
    llvmcall_expr(decl, join(instrs,"\n"), :(NTuple{$W,VecElement{$T}}), :(Tuple{Ptr{$T}, $U}), vtyp, [JULIAPOINTERTYPE, "i$(8sizeof(U))"], [:ptr, :(data(mask))])
end

@inline vload(::StaticInt{N}, args...) where {N} = StaticInt{N}()
@inline stridedpointer(::StaticInt{N}) where {N} = StaticInt{N}()
@inline zero_offsets(::StaticInt{N}) where {N} = StaticInt{N}()


text:
#Reconstruction

function minmod(a,b)
    if a*b < 0.0
        return 0.0
    elseif abs(a) < abs(b)
        return a
    else
        return b
    end
end

function tvd_minmod_reconstruction(n, g, f, x, xi)
    fp = zeros(n)
    fm = zeros(n)

    for i = g:(n-g+2)
        dx_up = x[i] - x[i-1]
        dx_down = x[i+1] - x[i]
        dx_m = x[i] - xi[i]
        dx_p = xi[i+1] - x[i]
        df_up = (f[i] - f[i-1])/dx_up
        df_down = (f[i+1]-f[i])/dx_down
        delta = minmod(df_up, df_down)
        fp[i] = f[i] + delta*dx_p
        fm[i] = f[i] - delta*dx_m
    end

    return fp, fm
end

signum(x,y) = y >= 0.0 ? abs(x) : -abs(x)

function tvd_mc_reconstruction(n, g, f, x, xi)
    fp = zeros(n)
    fm = zeros(n)

    for i = g:(n-g+2)
        dx_up = x[i] - x[i-1]
        dx_down = x[i+1] - x[i]
        dx_m = x[i] -xi[i]
        dx_p = xi[i+1] - x[i]
        df_up = (f[i]-f[i-1]) / dx_up
        df_down = (f[i+1]-f[i]) / dx_down

        if df_up*df_down < 0.0
            delta = 0.0
        else
            delta = signum(min(2.0abs(df_up),
                               2.0abs(df_down),
                               0.5(abs(df_up)+abs(df_down))),
                           df_up + df_down)
        end

        fp[i] = f[i] + delta*dx_p
        fm[i] = f[i] - delta*dx_m
    end

    return fp, fm
end

#2-dim reconstruction
function reconstruct(hyd::data2d)

  #x-pencils
  for j = (hyd.g-1):(hyd.ny-hyd.g+2)
    rhop, rhom = tvd_minmod_reconstruction(hyd.nx,
                                       hyd.g,
                                       vec(hyd.rho[j,:]),
                                       hyd.x,
                                       hyd.xi)

    hyd.rhop[j,:,1] = rhop'
    hyd.rhom[j,:,1] = rhom'

    epsp, epsm = tvd_minmod_reconstruction(hyd.nx,
                                       hyd.g,
                                       vec(hyd.eps[j,:]),
                                       hyd.x,
                                       hyd.xi)

    hyd.epsp[j,:,1] = epsp'
    hyd.epsm[j,:,1] = epsm'

    velxp, velxm = tvd_minmod_reconstruction(hyd.nx,
                                         hyd.g,
                                         vec(hyd.velx[j,:]),
                                         hyd.x,
                                         hyd.xi)

    hyd.velxp[j,:,1] = velxp'
    hyd.velxm[j,:,1] = velxm'

    velyp, velym = tvd_minmod_reconstruction(hyd.nx,
                                         hyd.g,
                                         vec(hyd.vely[j,:]),
                                         hyd.x,
                                         hyd.xi)

    hyd.velyp[j,:,1] = velyp'
    hyd.velym[j,:,1] = velym'
  end

  #y-pencils
  for i = (hyd.g-1):(hyd.nx-hyd.g+2)
    hyd.rhop[:,i,2], hyd.rhom[:,i,2] = tvd_minmod_reconstruction(hyd.ny,
                                                             hyd.g,
                                                             hyd.rho[:,i],
                                                             hyd.y,
                                                             hyd.yi)

    hyd.epsp[:,i,2], hyd.epsm[:,i,2] = tvd_minmod_reconstruction(hyd.ny,
                                                             hyd.g,
                                                             hyd.eps[:,i,1],
                                                             hyd.y,
                                                             hyd.yi)

    hyd.velxp[:,i,2], hyd.velxm[:,i,2] = tvd_minmod_reconstruction(hyd.ny,
                                                               hyd.g,
                                                               hyd.velx[:,i],
                                                               hyd.y,
                                                               hyd.yi)

    hyd.velyp[:,i,2], hyd.velym[:,i,2] = tvd_minmod_reconstruction(hyd.ny,
                                                               hyd.g,
                                                               hyd.vely[:,i],
                                                               hyd.y,
                                                               hyd.yi)
  end

  hyd.pressp[:,:,1] = eos_press(hyd.rhop[:,:,1], hyd.epsp[:,:,1], gamma)
  hyd.pressm[:,:,1] = eos_press(hyd.rhom[:,:,1], hyd.epsm[:,:,1], gamma)
  hyd.pressp[:,:,2] = eos_press(hyd.rhop[:,:,2], hyd.epsp[:,:,2], gamma)
  hyd.pressm[:,:,2] = eos_press(hyd.rhom[:,:,2], hyd.epsm[:,:,2], gamma)

  hyd.qp[:,:,:,1] = prim2con(hyd.rhop[:,:,1], hyd.velxp[:,:,1], hyd.velyp[:,:,1], hyd.epsp[:,:,1])
  hyd.qm[:,:,:,1] = prim2con(hyd.rhom[:,:,1], hyd.velxm[:,:,1], hyd.velym[:,:,1], hyd.epsm[:,:,1])
  hyd.qp[:,:,:,2] = prim2con(hyd.rhop[:,:,2], hyd.velxp[:,:,2], hyd.velyp[:,:,2], hyd.epsp[:,:,2])
  hyd.qm[:,:,:,2] = prim2con(hyd.rhom[:,:,2], hyd.velxm[:,:,2], hyd.velym[:,:,2], hyd.epsm[:,:,2])

  return hyd
end

max_stars_repo_name: gzahl/hydro

max_issues_repo_name: gzahl/hydro

max_forks_repo_name: gzahl/hydro


text:
function NullSpace(B, r)
    #This function returns the transformation of basis rendering the matrix B
    #into the form Bp=[eye(r) 0]

    #B is an array with dim1=n+2 rows and at least dim2=n+1 columns
    #r is the rank of B

    #P is the transformation matrix and is such that its last dim2-r columns span the null space of B
    #PermC is the permutation of columns throughout the process

    dim1, dim2 = size(B)
    P = Matrix(1.0I, dim2, dim2)
    PermC = 1:dim2

    Bp = copy(B)

    for a = 1:r
        temp = abs.(Bp[a:dim1, a:dim2])
        d1 = dim1 - a + 1
        d2 = dim2 - a + 1

        ind = argmax(reshape(temp, d1 * d2, 1))[1]
        column = ceil(Int64, ind / d1)
        row = a - 1 + ind - (column - 1) * d1
        column = a - 1 + column
        M = Bp[row, column]
        PerR = collect(1:dim1)
        PerR[a] = row
        PerR[row] = a
        PerC = collect(1:dim2)
        PerC[a] = column
        PerC[column] = a
        Bp = Bp[PerR, :]
        Bp = Bp[:, PerC]
        P = P[:, PerC]
        Bp[:, a] = Bp[:, a]/M
        P[:, a] = P[:, a]/M
        index = complementary(a, dim2)
        for i = index
            coef = Bp[a, i]
            Bp[:, i] = Bp[:, i] - coef * Bp[:, a]
            P[:, i] = P[:, i] - coef * P[:, a]
        end
        PerC
        PermC
        PermC = PermC[PerC]
    end
    # At this stage the columns in B corresponding to the indices PermC(1:r)
    # are linearly independent
    # In addition, since all the columns not chosen to be part of the first r
    # stack are simply changed as C -> C - a C' (and not permuted in any way), P is such that
    # Null(PermC,:) = [P';eye(dim2-r)]

    Null = P[:, (r + 1):dim2]

    return Null, PermC
end


text:
abstract type EnergyBC end

using ..TurbulenceClosures
using SurfaceFluxes: get_energy_flux, surface_conditions, DGScheme

"""
    Insulating() :: EnergyBC
No energy flux across the boundary.
"""
struct Insulating <: EnergyBC end
function atmos_energy_boundary_state!(nf, bc_energy::Insulating, atmos, _...) end
function atmos_energy_normal_boundary_flux_second_order!(
    nf,
    bc_energy::Insulating,
    atmos,
    _...,
) end

"""
    PrescribedTemperature(fn) :: EnergyBC
Prescribe the temperature at the boundary by `fn`, a function with signature
`fn(state, aux, t)` returning the temperature (in K).
"""
struct PrescribedTemperature{FN} <: EnergyBC
    fn::FN
end
function atmos_energy_boundary_state!(
    nf,
    bc_energy::PrescribedTemperature,
    atmos,
    state⁺,
    args,
)
    @unpack aux⁻, state⁻, t = args
    FT = eltype(aux⁻)
    param_set = parameter_set(atmos)
    _T_0::FT = T_0(param_set)
    _cv_d::FT = cv_d(param_set)

    T = bc_energy.fn(state⁻, aux⁻, t)
    E_int⁺ = state⁺.ρ * _cv_d * (T - _T_0)
    state⁺.energy.ρe =
        E_int⁺ + state⁺.ρ * gravitational_potential(atmos.orientation, aux⁻)
end

function atmos_energy_normal_boundary_flux_second_order!(
    nf,
    bc_energy::PrescribedTemperature,
    atmos,
    fluxᵀn,
    args,
)
    @unpack state⁻, aux⁻, diffusive⁻, hyperdiff⁻, t, n⁻ = args
    tend_type = Flux{SecondOrder}()
    _args⁻ = (;
        state = state⁻,
        aux = aux⁻,
        t,
        diffusive = diffusive⁻,
        hyperdiffusive = hyperdiff⁻,
    )
    pargs = merge(_args⁻, (precomputed = precompute(atmos, _args⁻, tend_type),))
    total_flux =
        Σfluxes(Energy(), eq_tends(Energy(), atmos, tend_type), atmos, pargs)
    nd_ρh_tot = dot(n⁻, total_flux)
    # both sides involve projections of normals, so signs are consistent
    fluxᵀn.energy.ρe += nd_ρh_tot
end

"""
    PrescribedEnergyFlux(fn) :: EnergyBC
Prescribe the net inward energy flux across the boundary by `fn`, a function
with signature `fn(state, aux, t)`, returning the flux (in W/m^2).
"""
struct PrescribedEnergyFlux{FN} <: EnergyBC
    fn::FN
end
function atmos_energy_boundary_state!(
    nf,
    bc_energy::PrescribedEnergyFlux,
    atmos,
    _...,
) end
function atmos_energy_normal_boundary_flux_second_order!(
    nf,
    bc_energy::PrescribedEnergyFlux,
    atmos,
    fluxᵀn,
    args,
)
    @unpack state⁻, aux⁻, t = args

    # DG normal is defined in the outward direction
    # we want to prescribe the inward flux
    fluxᵀn.energy.ρe -= bc_energy.fn(state⁻, aux⁻, t)
end

"""
    Adiabaticθ(fn) :: EnergyBC
Prescribe the net inward potential temperature flux
across the boundary by `fn`, a function with signature
`fn(state, aux, t)`, returning the flux (in kgK/m^2).
"""
struct Adiabaticθ{FN} <: EnergyBC
    fn::FN
end
function atmos_energy_boundary_state!(nf, bc_energy::Adiabaticθ, atmos, _...) end
function atmos_energy_normal_boundary_flux_second_order!(
    nf,
    bc_energy::Adiabaticθ,
    atmos,
    fluxᵀn,
    args,
)
    @unpack state⁻, aux⁻, t = args

    # DG normal is defined in the outward direction
    # we want to prescribe the inward flux
    fluxᵀn.energy.ρθ_liq_ice -= bc_energy.fn(state⁻, aux⁻, t)
end

"""
    BulkFormulaEnergy(fn) :: EnergyBC
Calculate the net inward energy flux across the boundary. The drag
coefficient is `C_h = fn_C_h(atmos, state, aux, t, normu_int_tan)`. The surface
temperature and q_tot are `T, q_tot = fn_T_and_q_tot(atmos, state, aux, t)`.
Return the flux (in W m^-2).
"""
struct BulkFormulaEnergy{FNX, FNTM} <: EnergyBC
    fn_C_h::FNX
    fn_T_and_q_tot::FNTM
end
function atmos_energy_boundary_state!(
    nf,
    bc_energy::BulkFormulaEnergy,
    atmos,
    _...,
) end
function atmos_energy_normal_boundary_flux_second_order!(
    nf,
    bc_energy::BulkFormulaEnergy,
    atmos,
    fluxᵀn,
    args,
)
    @unpack state⁻, aux⁻, t, state_int⁻, aux_int⁻ = args

    u_int⁻ = state_int⁻.ρu / state_int⁻.ρ
    u_int⁻_tan = projection_tangential(atmos, aux_int⁻, u_int⁻)
    normu_int⁻_tan = norm(u_int⁻_tan)
    C_h = bc_energy.fn_C_h(atmos, state⁻, aux⁻, t, normu_int⁻_tan)
    T, q_tot = bc_energy.fn_T_and_q_tot(atmos, state⁻, aux⁻, t)
    param_set = parameter_set(atmos)

    # calculate MSE from the states at the boundary and at the interior point
    ts = PhaseEquil_ρTq(param_set, state⁻.ρ, T, q_tot)
    ts_int = recover_thermo_state(atmos, state_int⁻, aux_int⁻)
    e_pot = gravitational_potential(atmos.orientation, aux⁻)
    e_pot_int = gravitational_potential(atmos.orientation, aux_int⁻)
    MSE = moist_static_energy(ts, e_pot)
    MSE_int = moist_static_energy(ts_int, e_pot_int)

    # TODO: use the correct density at the surface
    ρ_avg = average_density(state⁻.ρ, state_int⁻.ρ)
    # NOTE: difference from design docs since normal points outwards
    fluxᵀn.energy.ρe -= C_h * ρ_avg * normu_int⁻_tan * (MSE - MSE_int)
end

"""
    NishizawaEnergyFlux(fn) :: EnergyBC
Calculate the net inward energy flux across the boundary following Nishizawa and Kitamura (2018).
Return the flux (in W m^-2).
"""
struct NishizawaEnergyFlux{FNTM, FNX} <: EnergyBC
    fn_z0::FNX
    fn_T_and_q_tot::FNTM
end
function atmos_energy_boundary_state!(
    nf,
    bc_energy::NishizawaEnergyFlux,
    atmos,
    _...,
) end
function atmos_energy_normal_boundary_flux_second_order!(
    nf,
    bc_energy::NishizawaEnergyFlux,
    atmos,
    fluxᵀn,
    args,
)
    @unpack state⁻, aux⁻, t, aux_int⁻, state_int⁻ = args
    param_set = parameter_set(atmos)
    FT = eltype(state⁻)
    # Interior state
    u_int⁻ = state_int⁻.ρu / state_int⁻.ρ
    u_int⁻_tan = projection_tangential(atmos, aux_int⁻, u_int⁻)
    normu_int⁻_tan = norm(u_int⁻_tan)
    q_tot_int =
        moisture_model(atmos) isa DryModel ? FT(0) :
        state_int⁻.moisture.ρq_tot / state_int⁻.ρ
    # recover thermo state
    ts_int = recover_thermo_state(atmos, state_int⁻, aux_int⁻)
    x_in =
        MArray{Tuple{3}, FT}(FT[normu_int⁻_tan, dry_pottemp(ts_int), q_tot_int])

    # Boundary state
    T, q_tot = bc_energy.fn_T_and_q_tot(atmos, state⁻, aux⁻, t)
    x_s = MArray{Tuple{3}, FT}(FT[FT(0), T, q_tot])

    ## Initial guesses for MO parameters, these should be a function of state.
    LMO_init = FT(100) # Initial value so that ξ_init<<1
    u_star_init = FT(0.1) * normu_int⁻_tan
    th_star_init = T
    qt_star_init = q_tot
    MO_param_guess = MArray{Tuple{4}, FT}(FT[
        LMO_init,
        u_star_init,
        th_star_init,
        qt_star_init,
    ])

    # Roughness and interior heights
    z_0 = bc_energy.fn_z0(atmos, state⁻, aux⁻, t, normu_int⁻_tan)
    z_in = altitude(atmos, aux_int⁻)

    θ_flux, q_tot_flux = get_energy_flux(surface_conditions(
        param_set,
        MO_param_guess,
        x_in,
        x_s,
        z_0,
        T,
        z_in,
        DGScheme(),
    ))

    # recover thermo state
    ts_surf = PhaseEquil_ρTq(param_set, state⁻.ρ, T, q_tot)
    # Add sensible heat flux
    fluxᵀn.energy.ρe -= state⁻.ρ * θ_flux * cp_m(ts_surf)
    # Add latent heat flux
    fluxᵀn.energy.ρe -= state⁻.ρ * q_tot_flux * latent_heat_vapor(param_set, T)
end


text:
@testset "Collections" begin
  @testset "DomainCollection" begin
    # TODO
  end

  @testset "DataCollection" begin
    # TODO
  end
end