text:
@inline function Base.getindex(S::AbstractLowerTriangularMatrix{P,T,L}, i, j) where {P,T,L}
    @boundscheck i > P && ThrowBoundsError("i == $i > $P")
    j > i && return zero(T)
    @inbounds S.data[lt_sub2ind(P, i, j)]
end

@inline function Base.getindex(S::AbstractMutableLowerTriangularMatrix{P,T,L}, i, j) where {P,T,L}
    @boundscheck i > P && ThrowBoundsError("i == $i > $P")
    j > i && return zero(T)
    VectorizationBase.load(pointer(S) + (lt_sub2ind(P, i, j) - 1)*sizeof(T))
end

@inline function Base.getindex(S::AbstractUpperTriangularMatrix{P,T,L}, i, j) where {P,T,L}
    i > j && return zero(T)
    @boundscheck j > P && ThrowBoundsError("j == $j > $P.")
    @inbounds S.data[ut_sub2ind(P, i, j)]
end
@inline function Base.getindex(S::AbstractMutableUpperTriangularMatrix{P,T,L}, i, j) where {P,T,L}
    @boundscheck i > P && ThrowBoundsError("i == $i > $P")
    i > j && return zero(T)
    VectorizationBase.load(pointer(S) + (ut_sub2ind(P, i, j) - 1)*sizeof(T))
end
@inline function Base.getindex(S::LinearAlgebra.Adjoint{Union{},<:AbstractUpperTriangularMatrix{P,Vec{W,T},L}}, i::Int, j::Int) where {P,T,L,W}
    j > i && return zero(T)
    ind = ut_sub2ind(P, j, i)
    @boundscheck i > L && ThrowBoundsError("ind == $ind > $L.")
    @inbounds S.parent.data[ind]
end
@inline function Base.getindex(S::AbstractLowerTriangularMatrix{P,NTuple{W,Core.VecElement{T}},L}, i, j) where {P,W,T,L}
    @boundscheck i > P && ThrowBoundsError("i == $i > $P")
    j > i && return SIMDPirates.vbroadcast(Vec{W,T}, zero(T))
    @inbounds S.data[lt_sub2ind(P, i, j)]
end

@inline function Base.getindex(S::AbstractUpperTriangularMatrix{P,NTuple{W,Core.VecElement{T}},L}, i, j) where {P,W,T,L}
    @boundscheck j > P && ThrowBoundsError("j == $j > $P.")
    i > j && return SIMDPirates.vbroadcast(Vec{W,T}, zero(T))
    @inbounds S.data[ut_sub2ind(P, i, j)]
end

@generated function LinearAlgebra.det(A::AbstractTriangularMatrix{P,T,L}) where {P,T,L}
    quote
        out = one(T)
        @vvectorize for i ∈ 1:$P
            out *= A[i]
        end
        out
    end
end
"""
logdet(A) will be slower than log(det(A)),
but should be more numerically accurate.
det(A) is at high risk of over/underflow
for large matrices.
"""
@generated function LinearAlgebra.logdet(A::AbstractTriangularMatrix{P,T,L}) where {P,L,T}
    quote
        $(Expr(:meta,:inline))
        out = zero(T)
        @vvectorize $T for i ∈ 1:$P
            out += log(A[i])
        end
        out
    end
end

@generated function ∂logdet(A::AbstractTriangularMatrix{P,T,L}) where {P,T,L}
    quote
        $(Expr(:meta,:inline))
        out = zero(T)
        ∂out = PaddedMatrices.FixedSizeVector{$P,$T}(undef)
        @vvectorize $T for i ∈ 1:$P
            out += log(A[i])
            ∂out[i] = one($T) / A[i]
        end
        out, ∂out
    end
end

@inline nlogdet(A) = Base.FastMath.sub_fast(logdet(A))

#
# function lower_chol_small(P,T,L = calculate_L(P, T))
#     W, Wshift = VectorizationBase.pick_vector_width_shift(P, T)
#     V = Vec{W,T}
#     MT = VectorizationBase.mask_type(W)
#     q = quote
#         vΣ = VectorizationBase.vectorizable(Σ)
#         @inbounds @fastmath diag_1 = sqrt( vΣ[1] )
#         c_1 = SIMDPirates.vfdiv(vload($V, vΣ + $W), vbroadcast($V, diag_1))
#         diag² = vload($V, vΣ)
#         diag² = vifelse( one($MT), diag², SIMDPirates.vfnmadd(c_1, c_1, diag²) )
#         @inbounds @fastmath diag_2 = sqrt(diag²[2].value)
#     end
#     mask64 = 1
#     for p ∈ 2:P-1
#         c_p = Symbol(:c_, p)
#         push!(q.args, :(
#             $c_p = vload($V, vΣ + $W*$p)
#         ))
#         for j ∈ 1:p-1
#             i = p - j
#             c_prev = Symbol(:c_, j)
#             push!(q.args, :(
#                 $c_p = SIMDPirates.vfnmadd($c_prev, vbroadcast($V, $c_prev[$p]), $c_p)
#             ))
#         end
#         mask64 += 2^(p-1)
#         push!(q.args, quote
#             $c_p = SIMDPirates.vfdiv($c_p, vbroadcast($V, $(Symbol(:diag_, p))))
#             diag² = vifelse( $(MT(mask64)), diag², SIMDPirates.vfnmadd($c_p, $c_p, diag²) )
#             @inbounds @fastmath $(Symbol(:diag_, p+1)) = sqrt(diag²[$(p+1)].value)
#         end)
#     end
#
#     push!(q.args, quote
#         LowerTriangularMatrix{$P,$T,$L}(
#             $(diagonal_lowertri_output_tuple_expr(P, W, :diag_, :c_))
#         )
#     end)
#     q
# end
#
# @generated function lower_cholesky(Σ::SymmetricMatrixL{P,T,L}) where {P,T,L}
#     W, Wshift = VectorizationBase.pick_vector_width_shift(P, T)
#     Wm1 = W - 1
#     PpWm1 = P + Wm1
#     num_diag_loads = PpWm1 >>> Wshift
#     num_diag_loads == 1 && return lower_chol_small(P, T, L)
#     rem = P & Wm1
#     # num_diag_loads == 1 && return lower_chol_small(P, T, L)
#     num_full_diag_loads = P >>> Wshift
#     # Do we start with the full number of loads, or not?
#     # initial_diag_ind = lower_triangle_sub2ind(Val(P), T, 1, 1)
#     V = Vec{W,T}
#     MT = mask_type(W)
#     q = quote
#         vΣ = VectorizationBase.vectorizable(Σ)
#         @inbounds @fastmath diag_1 = sqrt( vΣ[1] )
#         @fastmath diag⁻¹_1 = 1 / diag_1
#         vdiag⁻¹_1 = vbroadcast($V, diag⁻¹_1)
#         $([
#             :($(Symbol(:diag²_, n)) = vload($V, vΣ + $(n*W)))
#             for n ∈ 0:num_full_diag_loads-1
#         ]...)
#     end
#     if num_full_diag_loads < num_diag_loads
#         push!(q.args, :($(Symbol(:diag²_, num_diag_loads)) = vload($V, vΣ + $(num_full_diag_loads*W), $(MT(2^(P & Wm1)-1)))))
#     end
#     num_col_loads = num_diag_loads
#     col_rem = W - (P & Wm1)
#     for col ∈ 1:P-1
#         col_rem += 1
#         if col_rem == W
#             col_rem = 0
#             num_col_loads -= 1
#         end
#         for load ∈ 1:num_col_loads
#             push!(q.args, :(
#                 $(Symbol(:c_, load, :_, 1)) = SIMDPirates.evmul(vdiag⁻¹_1, vload(vΣ + $(num_diag_loads*W + (load-1)*W)))
#             ))
#         end
#     end
#
# end

function load_packed_L_quote!(qa, P, symbol_name, extract_from)
    for p ∈ 1:P
        push!(qa, :($(PaddedMatrices.sym(symbol_name, p, p)) = $extract_from[$p] ))
    end
    ind = P
    for pc ∈ 1:P
        for pr ∈ pc+1:P
            ind += 1
            push!(qa, :($(PaddedMatrices.sym(symbol_name, pr, pc)) = $extract_from[$ind] ))
            # push!(qa, :($(PaddedMatrices.sym(symbol_name, pr, pc)) = $extract_from[$(lt_sub2ind(P,pc,pr))] ))
        end
    end
    qa
end
function store_packed_L_quote!(qa, P, symbol_name, T, L)
    outtup = Expr(:tuple,)
    for p ∈ 1:P
        push!(outtup.args, PaddedMatrices.sym(symbol_name, p, p) )
    end
    for pc ∈ 1:P
        for pr ∈ pc+1:P
            push!(outtup.args, PaddedMatrices.sym(symbol_name, pr, pc))
        end
    end
    for p ∈ binomial2(P+1)+1:L
        push!(outtup.args, zero(T))
    end
    # push!(qa, :(LowerTriangularMatrix{$P,$T,$L}($(
    #     Expr(:tuple, outtup...)
    # ))))
    :(LowerTriangularMatrix{$P,$T,$L}($outtup))
end
function store_packed_Ut_quote!(qa, P, symbol_name, T, L)
    outtup = Expr(:tuple,)
    for p ∈ 1:P
        push!(outtup.args, PaddedMatrices.sym(symbol_name, p, p) )
    end
    for pc ∈ 1:P
        for pr ∈ 1:pc-1
            push!(outtup.args, PaddedMatrices.sym(symbol_name, pc, pr))
        end
    end
    for p ∈ binomial2(P+1)+1:L
        push!(outtup.args, zero(T))
    end
    # push!(qa, :(UpperTriangularMatrix{$P,$T,$L}($(
    #     Expr(:tuple, outtup...)
    # ))))
    :(UpperTriangularMatrix{$P,$T,$L}($outtup))
end
function store_packed_U_quote!(qa, P, symbol_name, T, L)
    outtup = Expr(:tuple,)
    for p ∈ 1:P
        push!(outtup.args, PaddedMatrices.sym(symbol_name, p, p) )
    end
    for pc ∈ 1:P
        for pr ∈ 1:pc-1
            push!(outtup.args, PaddedMatrices.sym(symbol_name, pr, pc))
        end
    end
    for p ∈ binomial2(P+1)+1:L
        push!(outtup.args, zero(T))
    end
    # push!(qa, :(UpperTriangularMatrix{$P,$T,$L}($(
    #     Expr(:tuple, outtup...)
    # ))))
    :(UpperTriangularMatrix{$P,$T,$L}($outtup))
end

import PaddedMatrices: sym
∂sym(A, i, j, B, k, l) = Symbol(:∂, sym(A, i, j), :_∂, sym(B, k, l))
∂sym(A, B) = Symbol(:∂, A, :_∂, B)

function ∂inv_L_core_quote!(qa, P, output = :L, input = :S, ::Type{T} = Float64) where {T}
    # ∂tup = Expr(:tuple,)
    # The lazy way to write this function is to just use a dict
    # to find out what I have to do.
    ∂relations = Dict{Symbol,Set{Tuple{Symbol,Symbol}}}()
    defined_syms = Set{Symbol}()
    for c ∈ 1:P
        inputsym = sym(input, c, c)
        outputsym = sym(output, c, c)
        ∂outinsym = ∂sym(outputsym, inputsym)
        push!(qa, :( $outputsym = $(one(T)) / $inputsym ))
        push!(qa, :( $∂outinsym = - $outputsym*$outputsym ))
        # push!(∂tup.args, ∂outinsym)
        push!(defined_syms, ∂outinsym)
        push!(get!(() -> Set{Tuple{Symbol,Symbol}}(), ∂relations, outputsym), (∂outinsym, inputsym) )
    end
    for c ∈ 1:P
        for r ∈ c+1:P
            outputsym_r_c = sym(output, c, r)
            inputsym_r_c = sym(input, r, c)
            outputsym_c_c = sym(output, c, c)
            outputsym_r_r = sym(output, r, r)

            push!(qa, :( $outputsym_r_c = $inputsym_r_c * $outputsym_c_c ))
            push!(qa, :( $(∂sym(outputsym_r_c, inputsym_r_c)) = -$outputsym_r_r * $outputsym_c_c ))
            # expand partial to outputs with respect to inputs
            push!(qa, :( $(∂sym(outputsym_r_c, sym(input, c, c))) = -$outputsym_r_r * $inputsym_r_c * $(∂sym(outputsym_c_c, sym(input, c, c)) )))

            # if ∂sym(output, r, c, input, r, c) ∉ defined_syms
                # push!(∂tup.args, ∂sym(output, r, c, input, r, c))
                push!(defined_syms, ∂sym(output, c, r, input, r, c))
            # else
            #     throw("oops, $(∂sym(output, r, c, input, r, c)) ∈ $defined_syms")
            # end
            outrc_vec = get!(() -> Set{Tuple{Symbol,Symbol}}(), ∂relations, outputsym_r_c)
            push!(outrc_vec, (∂sym(output, c, r, input, r, c), sym(input, r, c)) )
            # if ∂sym(output, r, c, input, c, c) ∉ defined_syms
                # push!(∂tup.args, ∂sym(output, r, c, input, c, c))
                push!(defined_syms, ∂sym(output, c, r, input, c, c))
            # else
            #     throw("oops, $(∂sym(output, r, c, input, c, c)) ∈ $defined_syms")
            # end
            push!(qa, :( $(∂sym(output, c, r, input, r, r)) =  - $(sym(output, c, c)) * $(sym(input, r, c)) * $(∂sym(output, r, r, input, r, r)) ) )
            push!(outrc_vec, (∂sym(output, c, r, input, r, r), sym(input, r, r)) )
            push!(outrc_vec, (∂sym(output, c, r, input, c, c), sym(input, c, c)) )
            for cr ∈ c+1:r-1
                push!(qa, :( $outputsym_r_c += $(sym(input, r, cr)) * $(sym(output, c, cr)) ))
                push!(qa, :( $(∂sym(output, c, r, input, r, cr)) = -$(sym(output, r, r)) * $(sym(output, c, cr)) ))
                # push!(∂tup.args, ∂sym(output, r, c, input, r, cr))
                push!(outrc_vec, (∂sym(output, c, r, input, r, cr), sym(input, r, cr)) )
                # expand partial to outputs with respect to inputs
                ∂outoutsym = ∂sym(output, c, r, output, c, cr)
                push!(qa, :( $∂outoutsym = -$(sym(output, r, r)) * $(sym(input, r, cr)) ))
                outsym_cr_c = sym(output, c, cr)

                # push!(qa, :( $(∂sym(output, r, c, output, cr, c)) = -$(sym(output, r, r)) * $(sym(input, r, cr)) ))

                for (∂outinsym, inputsym) ∈ ∂relations[outsym_cr_c]
                    defsym = ∂sym(outputsym_r_c, inputsym)
                    push!(outrc_vec, (defsym, inputsym) )
                    if defsym ∈ defined_syms
                        push!(qa, :( $defsym += $∂outoutsym * $∂outinsym ) )
                    else
                        push!(qa, :( $defsym = $∂outoutsym * $∂outinsym ) )
                        # push!(∂tup.args, defsym)
                        push!(defined_syms, defsym)
                    end
                end
            end
            # expand partial to outputs with respect to inputs
            ∂outoutsym_r_r = ∂sym(output, c, r, output, r, r)
            ∂outoutsym_r_c = ∂sym(output, c, r, output, c, r)
            push!(qa, :( $∂outoutsym_r_r =  - $(sym(output, c, r)) ) )
            push!(qa, :( $∂outoutsym_r_c =  - $(sym(output, c, r)) ) )
            push!(qa, :( $outputsym_r_c *=  -$(sym(output, r, r)) ) )

            for (∂outinsym, inputsym) ∈ ∂relations[sym(output, r, r)]
                defsym = ∂sym(outputsym_r_c, inputsym)
                if defsym ∈ defined_syms
                    push!(qa, :( $defsym += $∂outoutsym_r_r * $∂outinsym ) )
                else
                    push!(qa, :( $defsym = $∂outoutsym_r_r * $∂outinsym ) )
                    # push!(∂tup.args, defsym)
                    push!(defined_syms, defsym)
                end
            end
        end
    end
    # ∂tup
end

@generated function lower_chol(Σ::SymmetricMatrixL{P,T,L}) where {P,T,L}
    q = quote end
    qa = q.args
    load_packed_L_quote!(qa, P, :Σ, :Σ)
    PaddedMatrices.chol_L_core_quote!(qa, P, :Σ, T)
    lq = store_packed_L_quote!(qa, P, :Σ, T, L)
    quote
        $(Expr(:meta,:inline))
        @fastmath @inbounds begin
            # begin
            $q
            $lq
        end
    end
end
@generated function lower_chol(Σ::AbstractFixedSizeMatrix{P,P,T,R}) where {P,T,R}
    q = quote end
    qa = q.args
    PaddedMatrices.load_L_quote!(qa, P, R, :Σ, :Σ)
    PaddedMatrices.chol_L_core_quote!(qa, P, :Σ, T)
    L = binomial2(P+1)
    Wm1 = VectorizationBase.pick_vector_width(L, T)-1
    L = (L + Wm1) & ~Wm1
    lq = store_packed_L_quote!(qa, P, :Σ, T, L)
    quote
        $(Expr(:meta,:inline))
        @fastmath @inbounds begin
            # begin
            $q
            $lq
        end
    end
end

@generated function PaddedMatrices.invchol(Σ::SymmetricMatrixL{P,T,L}) where {P,T,L}
    q = quote end
    qa = q.args
    load_packed_L_quote!(qa, P, :Σ, :Σ)
    PaddedMatrices.invchol_L_core_quote!(qa, P, :Lt, :Σ, T)
    uq = store_packed_Ut_quote!(qa, P, :Lt, T, L)
    quote
        $(Expr(:meta,:inline))
        @fastmath @inbounds begin
            # begin
            $q
            $uq
        end
    end
end
@generated function Base.inv(Lt::AbstractLowerTriangularMatrix{P,T,L}) where {P,T,L}
    q = quote end
    qa = q.args
    load_packed_L_quote!(qa, P, :Lt, :Lt)
    PaddedMatrices.inv_L_core_quote!(qa, P, :U, :Lt, T)
    uq = store_packed_Ut_quote!(qa, P, :U, T, L)
    quote
        $(Expr(:meta,:inline))
        @fastmath @inbounds begin
            # begin
            $q
            $uq
        end
    end
end
@generated function inv′(Lt::AbstractLowerTriangularMatrix{P,T,L}) where {P,T,L}
    q = quote end
    qa = q.args
    load_packed_L_quote!(qa, P, :Lt, :Lt)
    PaddedMatrices.inv_L_core_quote!(qa, P, :U, :Lt, T)
    uq = store_packed_Ut_quote!(qa, P, :U, T, L)
    quote
        $(Expr(:meta,:inline))
        @fastmath @inbounds begin
            # begin
            $q
            $uq
        end
    end
end

@generated function LinearAlgebra.mul!(
    out::AbstractMutableLowerTriangularMatrix{M,T},
    D::LinearAlgebra.Diagonal{T,<:AbstractFixedSizeVector{M,T,P}},
    L::AbstractLowerTriangularMatrix{M,T}
) where {M,T,P}
    W, Wshift = VectorizationBase.pick_vector_width_shift(M, T)
    Wm1 = W - 1
    V = Vec{W,T}
    q = quote end
    # Now, what remains of L is a (M-1) x (M-1) lower triangle
    # we will do W of these at a time.
    reps = M >>> Wshift
    rem = M & Wm1
    # handle rem first
    if rem > 0
        # If rem is smaller than half of W, we may use a smaller vector size here
        Wrem = VectorizationBase.pick_vector_width(rem, T)
        full_mask = (one(UInt) << Wrem) - one(UInt)
        rem_mask_type = VectorizationBase.mask_type(rem)
        # could be typemax(UInt) ???
        # but "unsafe_trunc" says "arbitrary value" is returned if this is greater
        # so seems like this is safer.
        miss = Wrem - rem
        base_ind = - miss
        triangle_ind = base_ind
        increment = M - 1
        initial_mask = Base.unsafe_trunc(rem_mask_type, ((one(UInt) << miss) - one(UInt)) ⊻ full_mask )
        push!(q.args, quote
            vd = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vD + $base_ind, $initial_mask )
            vl = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL + $base_ind, $initial_mask )
            SIMDPirates.vstore!(vout + $base_ind, SIMDPirates.vmul(vd, vl), $initial_mask)
        end)
        triangle_ind += increment
        increment -= 1
        for r ∈ 1:rem-1
            vl = Symbol(:vl_,r)
            mask = Base.unsafe_trunc(rem_mask_type, ((one(UInt) << (miss+r)) - one(UInt)) ⊻ full_mask )
            push!(q.args, quote
                $vl = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL + $triangle_ind, $mask )
                SIMDPirates.vstore!(
                    vout + $triangle_ind,
                    SIMDPirates.vmul(vd, $vl), $mask
                )
            end)
            triangle_ind += increment
            increment -= 1
        end
        base_ind = rem
    else
        base_ind = 0
    end
    # then do reps of W
    full_mask = (one(UInt) << W) - one(UInt)
    if reps > 0
        rem_quote = quote end
        mask_type = VectorizationBase.mask_type(W)
        for w ∈ 1:Wm1
            vl = Symbol(:vld1_,w)
            mask = Base.unsafe_trunc(mask_type, ((one(UInt) << w) - one(UInt)) ⊻ full_mask )
            push!(rem_quote.args, quote
                $vl = SIMDPirates.vload(SIMDPirates.Vec{$W,$T}, vL + triangle_ind, $mask)
                SIMDPirates.vstore!(vout + triangle_ind, SIMDPirates.vmul(vd, $vl), $mask)
                triangle_ind += increment
                increment -= 1
            end)
        end
        push!(q.args, quote
            for rep ∈ 0:$(reps-1)
                col1ind = $base_ind + $W * rep
                # load diagonals
                vd = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vD + col1ind )
                vl = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL + col1ind )
                SIMDPirates.vstore!(vout + col1ind, SIMDPirates.vmul(vd, vl))
                triangle_ind = col1ind + $(M-1)
                increment = $(M - 2)
                for r ∈ 0:($(rem-1) + $W*rep)
                    vli = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL + triangle_ind )
                    SIMDPirates.vstore!(vout + triangle_ind, SIMDPirates.vmul(vd, vli))
                    triangle_ind += increment
                    increment -= 1
                end
                $rem_quote
            end
        end)
    end
    quote
        $(Expr(:meta,:inline))
        d = D.diag
        vD = VectorizationBase.vectorizable(d)
        vL = VectorizationBase.vectorizable(L)
        vout = VectorizationBase.vectorizable(out)
        GC.@preserve d L out begin
            $q
        end
        out
    end
end
function Base.:*(
    D::LinearAlgebra.Diagonal{T,<:AbstractFixedSizeVector{M,T,P}},
    L::AbstractLowerTriangularMatrix{M,T}
) where {M,T,P}
    out = MutableLowerTriangularMatrix{M,T}(undef)
    mul!(out, D, L)
end
function Base.:*(
    sp::StackPointer,
    D::LinearAlgebra.Diagonal{T,<:AbstractFixedSizeVector{M,T,P}},
    L::AbstractLowerTriangularMatrix{M,T}
) where {M,T,P}
    sp, out = PtrLowerTriangularMatrix{M,T}(sp)
    sp, mul!(out, D, L)
end

@generated function Base.muladd(
    D::LinearAlgebra.Diagonal{T,<:AbstractFixedSizeVector{M,T,P}},
    L::AbstractLowerTriangularMatrix{M,T},
    A::AbstractLowerTriangularMatrix{M,T}
) where {M,T,P}
    W, Wshift = VectorizationBase.pick_vector_width_shift(M, T)
    Wm1 = W - 1
    V = Vec{W,T}
    q = quote end
    # Now, what remains of L is a (M-1) x (M-1) lower triangle
    # we will do W of these at a time.
    reps = M >>> Wshift
    rem = M & Wm1
    # handle rem first
    if rem > 0
        # If rem is smaller than half of W, we may use a smaller vector size here
        Wrem = VectorizationBase.pick_vector_width(rem, T)
        full_mask = (one(UInt) << Wrem) - one(UInt)
        rem_mask_type = VectorizationBase.mask_type(rem)
        # could be typemax(UInt) ???
        # but "unsafe_trunc" says "arbitrary value" is returned if this is greater
        # so seems like this is safer.

        miss = Wrem - rem
        base_ind = - miss
        triangle_ind = base_ind
        increment = M - 1

        initial_mask = Base.unsafe_trunc(rem_mask_type, ((one(UInt) << miss) - one(UInt)) ⊻ full_mask )
        push!(q.args, quote
            vd = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vD + $base_ind, $initial_mask )
            vl = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL + $base_ind, $initial_mask )
            va = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vA + $base_ind, $initial_mask )
            SIMDPirates.vstore!(vout + $base_ind, SIMDPirates.vmuladd(vd, vl, va), $initial_mask)
        end)
        triangle_ind += increment
        increment -= 1
        for r ∈ 1:rem-1
            vl = Symbol(:vl_,r)
            va = Symbol(:va_,r)
            mask = Base.unsafe_trunc(rem_mask_type, ((one(UInt) << (miss+r)) - one(UInt)) ⊻ full_mask )
            push!(q.args, quote
                $vl = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL + $triangle_ind, $mask )
                $va = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vA + $triangle_ind, $mask )
                SIMDPirates.vstore!(
                    vout + $triangle_ind,
                    SIMDPirates.vmuladd(vd, $vl, $va), $mask
                )
            end)
            triangle_ind += increment
            increment -= 1
        end
        base_ind = rem
    else
        base_ind = 0
    end
    # then do reps of W
    full_mask = (one(UInt) << W) - one(UInt)
    if reps > 0
        rem_quote = quote end
        mask_type = VectorizationBase.mask_type(W)
        for w ∈ 1:Wm1
            vl = Symbol(:vl_,w)
            va = Symbol(:va_,w)
            mask = Base.unsafe_trunc(mask_type, ((one(UInt) << w) - one(UInt)) ⊻ full_mask )
            push!(rem_quote.args, quote
                $vl = SIMDPirates.vload(SIMDPirates.Vec{$W,$T}, vL + triangle_ind, $mask)
                $va = SIMDPirates.vload(SIMDPirates.Vec{$W,$T}, vA + triangle_ind, $mask)
                SIMDPirates.vstore!(vout + triangle_ind, SIMDPirates.vmuladd(vd, $vl, $va), $mask)
                triangle_ind += increment
                increment -= 1
            end)
        end
        push!(q.args, quote
            for rep ∈ 0:$(reps-1)
                col1ind = $base_ind + $W * rep
                # load diagonals
                vd = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vD + col1ind )
                vl = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL + col1ind )
                va = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vA + col1ind )
                SIMDPirates.vstore!(vout + col1ind, SIMDPirates.vmuladd(vd, vl, va))
                triangle_ind = col1ind + $(M-1)
                increment = $(M - 2)
                for r ∈ 0:($(rem-1) + $W*rep)
                    vli = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL + triangle_ind )
                    vai = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vA + triangle_ind )
                    SIMDPirates.vstore!(vout + triangle_ind, SIMDPirates.vmuladd(vd, vli, vai))
                    triangle_ind += increment
                    increment -= 1
                end
                $rem_quote
            end
        end)
    end
    quote
        $(Expr(:meta,:inline))
        d = D.diag
        out = MutableLowerTriangularMatrix{$M,$T}(undef)
        vD = VectorizationBase.vectorizable(d)
        vL = VectorizationBase.vectorizable(L)
        vA = VectorizationBase.vectorizable(A)
        vout = VectorizationBase.vectorizable(out)
        GC.@preserve d L vA out begin
            $q
        end
        out
    end
end

@generated function row_sum_prod(
    L1::AbstractLowerTriangularMatrix{M,T},
    L2::AbstractLowerTriangularMatrix{M,T}
) where {M,T}
    W, Wshift = VectorizationBase.pick_vector_width_shift(M, T)
    Wm1 = W - 1
    V = Vec{W,T}
    q = quote end
    # Now, what remains of L is a (M-1) x (M-1) lower triangle
    # we will do W of these at a time.
    reps = M >>> Wshift
    rem = M & Wm1
    # handle rem first
    if rem > 0
        # If rem is smaller than half of W, we may use a smaller vector size here
        Wrem = VectorizationBase.pick_vector_width(rem, T)
        full_mask = (one(UInt) << Wrem) - one(UInt)
        rem_mask_type = VectorizationBase.mask_type(rem)
        # could be typemax(UInt) ???
        # but "unsafe_trunc" says "arbitrary value" is returned if this is greater
        # so seems like this is safer.
        miss = Wrem - rem
        base_ind = - miss
        triangle_ind = base_ind
        increment = M - 1
        initial_mask = Base.unsafe_trunc(rem_mask_type, ((one(UInt) << miss) - one(UInt)) ⊻ full_mask )
        push!(q.args, quote
            vld1 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL1 + $base_ind, $initial_mask )
            vld2 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL2 + $base_ind, $initial_mask )
            vcumulative = SIMDPirates.vmul(vld1, vld2)
        end)
        triangle_ind += increment
        increment -= 1
        for r ∈ 1:rem-1
            vld1 = Symbol(:vld1_,r)
            vld2 = Symbol(:vld2_,r)
            mask = Base.unsafe_trunc(rem_mask_type, ((one(UInt) << (miss+r)) - one(UInt)) ⊻ full_mask )
            push!(q.args, quote
                $vld1 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL1 + $triangle_ind, $mask )
                $vld2 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL2 + $triangle_ind, $mask )
                vcumulative = SIMDPirates.vmuladd($vld1, $vld2, vcumulative)
            end)
            triangle_ind += increment
            increment -= 1
        end
        push!(q.args, :(SIMDPirates.vstore!(vout + $base_ind, vcumulative, $initial_mask)))
        base_ind = rem
    else
        base_ind = 0
    end
    # then do reps of W
    full_mask = (one(UInt) << W) - one(UInt)
    if reps > 0
        rem_quote = quote end
        mask_type = VectorizationBase.mask_type(W)
        for w ∈ 1:Wm1
            vld1 = Symbol(:vld1_,w)
            vld2 = Symbol(:vld2_,w)
            mask = Base.unsafe_trunc(mask_type, ((one(UInt) << w) - one(UInt)) ⊻ full_mask )
            push!(rem_quote.args, quote
                $vld1 = SIMDPirates.vload(SIMDPirates.Vec{$W,$T}, vL1 + triangle_ind, $mask)
                $vld2 = SIMDPirates.vload(SIMDPirates.Vec{$W,$T}, vL2 + triangle_ind, $mask)
                vcumulative = SIMDPirates.vmuladd($vld1, $vld2, vcumulative)
                triangle_ind += increment
                increment -= 1
            end)
        end
        push!(q.args, quote
            for rep ∈ 0:$(reps-1)
                col1ind = $base_ind + $W * rep
                # load diagonals
                vld1 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL1 + col1ind )
                vld2 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL2 + col1ind )
                vcumulative = SIMDPirates.vmul(vld1, vld2)

                triangle_ind = col1ind + $(M-1)
                increment = $(M - 2)

                for r ∈ 0:($(rem-1) + $W*rep)
                    vld1 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL1 + triangle_ind )
                    vld2 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL2 + triangle_ind )
                    vcumulative = SIMDPirates.vmuladd(vld1, vld2, vcumulative)
                    triangle_ind += increment
                    increment -= 1
                end
                $rem_quote
                SIMDPirates.vstore!(vout + col1ind, vcumulative)
            end
        end)
    end
    quote
        $(Expr(:meta,:inline))
        out = FixedSizeVector{$M,$T}(undef)
        vL1 = VectorizationBase.vectorizable(L1)
        vL2 = VectorizationBase.vectorizable(L2)
        vout = VectorizationBase.vectorizable(out)
        GC.@preserve L1 L2 out begin
            $q
        end
        out
    end
end

@generated function row_sum_prod_add(
    L1::AbstractLowerTriangularMatrix{M,T},
    L2::AbstractLowerTriangularMatrix{M,T},
    v::AbstractFixedSizeVector{M,T}
) where {M,T}
    W, Wshift = VectorizationBase.pick_vector_width_shift(M, T)
    Wm1 = W - 1
    V = Vec{W,T}
    q = quote end
    # Now, what remains of L is a (M-1) x (M-1) lower triangle
    # we will do W of these at a time.
    reps = M >>> Wshift
    rem = M & Wm1
    # handle rem first
    if rem > 0
        # If rem is smaller than half of W, we may use a smaller vector size here
        Wrem = VectorizationBase.pick_vector_width(rem, T)
        full_mask = (one(UInt) << Wrem) - one(UInt)
        rem_mask_type = VectorizationBase.mask_type(rem)
        # could be typemax(UInt) ???
        # but "unsafe_trunc" says "arbitrary value" is returned if this is greater
        # so seems like this is safer.
        miss = Wrem - rem
        base_ind = - miss
        triangle_ind = base_ind
        increment = M - 1
        initial_mask = Base.unsafe_trunc(rem_mask_type, ((one(UInt) << miss) - one(UInt)) ⊻ full_mask )
        push!(q.args, quote
            vld1 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL1 + $base_ind, $initial_mask )
            vld2 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL2 + $base_ind, $initial_mask )
            vcumulative = SIMDPirates.vmuladd(vld1, vld2, SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vv + $base_ind, $initial_mask ))
        end)
        triangle_ind += increment
        increment -= 1
        for r ∈ 1:rem-1
            vld1 = Symbol(:vld1_,r)
            vld2 = Symbol(:vld2_,r)
            mask = Base.unsafe_trunc(rem_mask_type, (one(UInt) << (miss+r) - one(UInt)) ⊻ full_mask )
            push!(q.args, quote
                $vld1 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL1 + $triangle_ind, $mask )
                $vld2 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL2 + $triangle_ind, $mask )
                vcumulative = SIMDPirates.vmuladd($vld1, $vld2, vcumulative)
            end)
            triangle_ind += increment
            increment -= 1
        end
        push!(q.args, :(SIMDPirates.vstore!(vout + $base_ind, vcumulative, $initial_mask)))
        base_ind = rem
    else
        base_ind = 0
    end
    # then do reps of W
    full_mask = (one(UInt) << W) - one(UInt)
    if reps > 0
        rem_quote = quote end
        mask_type = VectorizationBase.mask_type(W)
        for w ∈ 1:Wm1
            vld1 = Symbol(:vld1_,w)
            vld2 = Symbol(:vld2_,w)
            mask = Base.unsafe_trunc(mask_type, ((one(UInt) << w) - one(UInt)) ⊻ full_mask )
            push!(rem_quote.args, quote
                $vld1 = SIMDPirates.vload(SIMDPirates.Vec{$W,$T}, vL1 + triangle_ind, $mask)
                $vld2 = SIMDPirates.vload(SIMDPirates.Vec{$W,$T}, vL2 + triangle_ind, $mask)
                vcumulative = SIMDPirates.vmuladd($vld1, $vld2, vcumulative)
                triangle_ind += increment
                increment -= 1
            end)
        end
        push!(q.args, quote
            for rep ∈ 0:$(reps-1)
                col1ind = $base_ind + $W * rep
                # load diagonals
                vld1 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL1 + col1ind )
                vld2 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL2 + col1ind )
                vcumulative = SIMDPirates.vmuladd(vld1, vld2, SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vv + col1ind))
                triangle_ind = col1ind + $(M-1)
                increment = $(M - 2)
                for r ∈ 0:($(rem-1) + $W*rep)
                    vld1 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL1 + triangle_ind )
                    vld2 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL2 + triangle_ind )
                    vcumulative = SIMDPirates.vmuladd(vld1, vld2, vcumulative)
                    triangle_ind += increment
                    increment -= 1
                end
                $rem_quote
                SIMDPirates.vstore!(vout + col1ind, vcumulative)
            end
        end)
    end
    quote
        $(Expr(:meta,:inline))
        out = FixedSizeVector{$M,$T}(undef)
        vv = VectorizationBase.vectorizable(v)
        vL1 = VectorizationBase.vectorizable(L1)
        vL2 = VectorizationBase.vectorizable(L2)
        vout = VectorizationBase.vectorizable(out)
        GC.@preserve v L1 L2 out begin
            $q
        end
        out
    end
end

function muladd_quote(M,T,P,N,Asym = :A)
    vAsym = Symbol(:v, Asym)
    W, Wshift = VectorizationBase.pick_vector_width_shift(M, T)
    Wm1 = W - 1
    V = Vec{W,T}
    q = quote end
    # Now, what remains of L is a (M-1) x (M-1) lower triangle
    # we will do W of these at a time.
    reps = M >>> Wshift
    rem = M & Wm1

    # handle rem first
    if rem > 0
        # If rem is smaller than half of W, we may use a smaller vector size here
        Wrem = VectorizationBase.pick_vector_width(rem, T)
        full_mask = (one(UInt) << Wrem) - one(UInt)
        rem_mask_type = VectorizationBase.mask_type(rem)
        # could be typemax(UInt) ???
        # but "unsafe_trunc" says "arbitrary value" is returned if this is greater
        # so seems like this is safer.

        miss = Wrem - rem
        base_ind = - miss
        triangle_ind = base_ind
        increment = M - 1

        initial_mask = Base.unsafe_trunc(rem_mask_type, ((one(UInt) << miss) - one(UInt)) ⊻ full_mask )
        push!(q.args, quote
            vd = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vD + $base_ind, $initial_mask )
            vl = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL + $base_ind, $initial_mask )
            va = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, $vAsym + $base_ind, $initial_mask )
            SIMDPirates.vstore!(vout + $base_ind, SIMDPirates.vmuladd(vd, vl, va), $initial_mask)
        end)
        triangle_ind += increment
        increment -= 1
        for r ∈ 1:rem-1
            vl = Symbol(:vl_,r)
            va = Symbol(:va_,r)
            mask = Base.unsafe_trunc(rem_mask_type, ((one(UInt) << (miss+r)) - one(UInt)) ⊻ full_mask )
            push!(q.args, quote
                $vl = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL + $triangle_ind, $mask )
                $va = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, $vAsym + $triangle_ind, $mask )
                SIMDPirates.vstore!(
                    vout + $triangle_ind,
                    SIMDPirates.vmuladd(vd, $vl, $va), $mask
                )
            end)
            triangle_ind += increment
            increment -= 1
        end
        base_ind = rem
    else
        base_ind = 0
    end

    # then do reps of W
    full_mask = (one(UInt) << W) - one(UInt)
    if reps > 0
        rem_quote = quote end
        mask_type = VectorizationBase.mask_type(W)
        for w ∈ 1:Wm1
            vl = Symbol(:vl_,w)
            va = Symbol(:va_,w)
            mask = Base.unsafe_trunc(mask_type, ((one(UInt) << w) - one(UInt)) ⊻ full_mask )
            push!(rem_quote.args, quote
                $vl = SIMDPirates.vload(SIMDPirates.Vec{$W,$T}, vL + triangle_ind, $mask)
                $va = SIMDPirates.vload(SIMDPirates.Vec{$W,$T}, $vAsym + triangle_ind, $mask)
                SIMDPirates.vstore!(vout + triangle_ind, SIMDPirates.vmuladd(vd, $vl, $va), $mask)
                triangle_ind += increment
                increment -= 1
            end)
        end
        push!(q.args, quote
            for rep ∈ 0:$(reps-1)
                col1ind = $base_ind + $W * rep
                # load diagonals
                vd = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vD + col1ind )
                vl = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL + col1ind )
                va = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, $vAsym + col1ind )
                SIMDPirates.vstore!(vout + col1ind, SIMDPirates.vmuladd(vd, vl, va))
                triangle_ind = col1ind + $(M-1)
                increment = $(M - 2)
                for r ∈ 0:($(rem-1) + $W*rep)
                    vli = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL + triangle_ind )
                    vai = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, $vAsym + triangle_ind )
                    SIMDPirates.vstore!(vout + triangle_ind, SIMDPirates.vmuladd(vd, vli, vai))
                    triangle_ind += increment
                    increment -= 1
                end
                $rem_quote
            end
        end)
    end
    q = quote
#        $(Expr(:meta,:inline))
        d = D.diag
        vD = VectorizationBase.vectorizable(d)
        vL = VectorizationBase.vectorizable(L)
        vout = VectorizationBase.vectorizable(out)
        GC.@preserve d L out begin
            $q
        end
        out
    end
    if Asym != :out
        pushfirst!(q.args, :($vAsym = VectorizationBase.vectorizable($Asym)))
    end
    q
end

@generated function muladd!(
    out::AbstractLowerTriangularMatrix{M,T,N},
    D::LinearAlgebra.Diagonal{T,<:AbstractFixedSizeVector{M,T,P}},
    L::AbstractLowerTriangularMatrix{M,T,N},
    A::AbstractLowerTriangularMatrix{M,T,N}
) where {M,T,P,N}
    muladd_quote(M,T,P,N, :A)
end
@generated function muladd!(
    out::AbstractLowerTriangularMatrix{M,T,N},
    D::LinearAlgebra.Diagonal{T,<:AbstractFixedSizeVector{M,T,P}},
    L::AbstractLowerTriangularMatrix{M,T,N}
) where {M,T,P,N}
    muladd_quote(M, T, P, N, :out)
end
function Base.muladd(
    sp::StackPointer,
    D::LinearAlgebra.Diagonal{T,<:AbstractFixedSizeVector{M,T}},
    L::AbstractLowerTriangularMatrix{M,T},
    A::AbstractLowerTriangularMatrix{M,T}
) where {M,T}
    (sp,out) = PtrLowerTriangularMatrix{M,T}(sp)
    muladd!(out, D, L, A)
    sp, out
end

@generated function row_sum_prod!(
    out::PaddedMatrices.AbstractMutableFixedSizeVector{M,T},
    L1::AbstractLowerTriangularMatrix{M,T},
    L2::AbstractLowerTriangularMatrix{M,T}
) where {M,T}
    W, Wshift = VectorizationBase.pick_vector_width_shift(M, T)
    Wm1 = W - 1
    V = Vec{W,T}
    q = quote end
    # Now, what remains of L is a (M-1) x (M-1) lower triangle
    # we will do W of these at a time.
    reps = M >>> Wshift
    rem = M & Wm1
    # handle rem first
    if rem > 0
        # If rem is smaller than half of W, we may use a smaller vector size here
        Wrem = VectorizationBase.pick_vector_width(rem, T)
        full_mask = (one(UInt) << Wrem) - one(UInt)
        rem_mask_type = VectorizationBase.mask_type(rem)
        # could be typemax(UInt) ???
        # but "unsafe_trunc" says "arbitrary value" is returned if this is greater
        # so seems like this is safer.

        miss = Wrem - rem
        base_ind = - miss
        triangle_ind = base_ind
        increment = M - 1

        initial_mask = Base.unsafe_trunc(rem_mask_type, ((one(UInt) << miss) - one(UInt)) ⊻ full_mask )
        push!(q.args, quote
            vld1 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL1 + $base_ind, $initial_mask )
            vld2 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL2 + $base_ind, $initial_mask )
            vcumulative = SIMDPirates.vmul(vld1, vld2)
        end)
        triangle_ind += increment
        increment -= 1
        for r ∈ 1:rem-1
            vld1 = Symbol(:vld1_,r)
            vld2 = Symbol(:vld2_,r)
            mask = Base.unsafe_trunc(rem_mask_type, ((one(UInt) << (miss+r)) - one(UInt)) ⊻ full_mask )
            push!(q.args, quote
                $vld1 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL1 + $triangle_ind, $mask )
                $vld2 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL2 + $triangle_ind, $mask )
                vcumulative = SIMDPirates.vmuladd($vld1, $vld2, vcumulative)
            end)
            triangle_ind += increment
            increment -= 1
        end
        push!(q.args, :(SIMDPirates.vstore!(vout + $base_ind, vcumulative, $initial_mask)))
        base_ind = rem
    else
        base_ind = 0
    end
    # then do reps of W
    full_mask = (one(UInt) << W) - one(UInt)
    if reps > 0
        rem_quote = quote end
        mask_type = VectorizationBase.mask_type(W)
        for w ∈ 1:Wm1
            vld1 = Symbol(:vld1_,w)
            vld2 = Symbol(:vld2_,w)
            mask = Base.unsafe_trunc(mask_type, ((one(UInt) << w) - one(UInt)) ⊻ full_mask )
            push!(rem_quote.args, quote
                $vld1 = SIMDPirates.vload(SIMDPirates.Vec{$W,$T}, vL1 + triangle_ind, $mask)
                $vld2 = SIMDPirates.vload(SIMDPirates.Vec{$W,$T}, vL2 + triangle_ind, $mask)
                vcumulative = SIMDPirates.vmuladd($vld1, $vld2, vcumulative)
                triangle_ind += increment
                increment -= 1
            end)
        end
        push!(q.args, quote
            for rep ∈ 0:$(reps-1)
                col1ind = $base_ind + $W * rep
              # load diagonals
                vld1 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL1 + col1ind )
                vld2 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL2 + col1ind )
                vcumulative = SIMDPirates.vmul(vld1, vld2)

                triangle_ind = col1ind + $(M-1)
                increment = $(M - 2)

                for r ∈ 0:($(rem-1) + $W*rep)
                    vld1 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL1 + triangle_ind )
                    vld2 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL2 + triangle_ind )
                    vcumulative = SIMDPirates.vmuladd(vld1, vld2, vcumulative)
                    triangle_ind += increment
                    increment -= 1
                end
                $rem_quote
                SIMDPirates.vstore!(vout + col1ind, vcumulative)
            end
        end)
    end
    quote
        # $(Expr(:meta,:inline))
        vL1 = VectorizationBase.vectorizable(L1)
        vL2 = VectorizationBase.vectorizable(L2)
        vout = VectorizationBase.vectorizable(out)
        GC.@preserve L1 L2 begin
            $q
        end
        out
    end
end
function row_sum_prod(
    sp::StackPointer,
    L1::AbstractLowerTriangularMatrix{M,T},
    L2::AbstractLowerTriangularMatrix{M,T}
) where {M,T}
    (sp, out) = PtrVector{M,T}(sp)
    row_sum_prod!(out, L1, L2)
    sp, out
end

function row_sum_add_quote(M, T, vaddsym = :v)
    vvaddsym = Symbol(:v, vaddsym)
    W, Wshift = VectorizationBase.pick_vector_width_shift(M, T)
    Wm1 = W - 1
    V = Vec{W,T}
    q = quote end
    # Now, what remains of L is a (M-1) x (M-1) lower triangle
    # we will do W of these at a time.
    reps = M >>> Wshift
    rem = M & Wm1
    # handle rem first # why? This misaligns us?
    if rem > 0
        # If rem is smaller than half of W, we may use a smaller vector size here
        Wrem = VectorizationBase.pick_vector_width(rem, T)
        full_mask = (one(UInt) << Wrem) - one(UInt)
        rem_mask_type = VectorizationBase.mask_type(rem)
        # could be typemax(UInt) ???
        # but "unsafe_trunc" says "arbitrary value" is returned if this is greater
        # so seems like this is safer.

        miss = Wrem - rem
        base_ind = - miss
        triangle_ind = base_ind
        increment = M - 1

        initial_mask = Base.unsafe_trunc(rem_mask_type, ((one(UInt) << miss) - one(UInt)) ⊻ full_mask )
        push!(q.args, quote
            vld1 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL1 + $base_ind, $initial_mask )
            vld2 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL2 + $base_ind, $initial_mask )
            vcumulative = SIMDPirates.vmuladd(vld1, vld2, SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, $vvaddsym + $base_ind, $initial_mask ))
        end)
        triangle_ind += increment
        increment -= 1
        for r ∈ 1:rem-1
            vld1 = Symbol(:vld1_,r)
            vld2 = Symbol(:vld2_,r)
            mask = Base.unsafe_trunc(rem_mask_type, ((one(UInt) << (miss+r)) - one(UInt)) ⊻ full_mask )
            push!(q.args, quote
                $vld1 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL1 + $triangle_ind, $mask )
                $vld2 = SIMDPirates.vload(SIMDPirates.Vec{$Wrem, $T}, vL2 + $triangle_ind, $mask )
                vcumulative = SIMDPirates.vmuladd($vld1, $vld2, vcumulative)
            end)
            triangle_ind += increment
            increment -= 1
        end
        push!(q.args, :(SIMDPirates.vstore!(vout + $base_ind, vcumulative, $initial_mask)))
        base_ind = rem
    else
        base_ind = 0
    end
    # then do reps of W
    full_mask = (one(UInt) << W) - one(UInt)
    if reps > 0
        rem_quote = quote end
        mask_type = VectorizationBase.mask_type(W)
        for w ∈ 1:Wm1
            vld1 = Symbol(:vld1_,w)
            vld2 = Symbol(:vld2_,w)
            mask = Base.unsafe_trunc(mask_type, ((one(UInt) << w) - one(UInt)) ⊻ full_mask )
            q_rem = quote
                $vld1 = SIMDPirates.vload(SIMDPirates.Vec{$W,$T}, vL1 + triangle_ind, $mask)
                $vld2 = SIMDPirates.vload(SIMDPirates.Vec{$W,$T}, vL2 + triangle_ind, $mask)
                vcumulative = SIMDPirates.vmuladd($vld1, $vld2, vcumulative)
                triangle_ind += increment
                increment -= 1
            end
            push!(rem_quote.args, q_rem)
        end
        q_rep = quote
            for rep ∈ 0:$(reps-1)
                col1ind = $base_ind + $W * rep
                # load diagonals
                vld1 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL1 + col1ind )
                vld2 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL2 + col1ind )
                vcumulative = SIMDPirates.vmuladd(vld1, vld2, SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, $vvaddsym + col1ind))
                triangle_ind = col1ind + $(M-1)
                increment = $(M - 2)
                for r ∈ 0:($(rem-1) + $W*rep)
                    vld1 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL1 + triangle_ind )
                    vld2 = SIMDPirates.vload(SIMDPirates.Vec{$W, $T}, vL2 + triangle_ind )
                    vcumulative = SIMDPirates.vmuladd(vld1, vld2, vcumulative)
                    triangle_ind += increment
                    increment -= 1
                end
                $rem_quote
                SIMDPirates.vstore!(vout + col1ind, vcumulative)
            end
        end
        push!(q.args, q_rep)
    end
    q = quote
        # $(Expr(:meta,:inline))
        vL1 = VectorizationBase.vectorizable(L1)
        vL2 = VectorizationBase.vectorizable(L2)
        vout = VectorizationBase.vectorizable(out)
        GC.@preserve L1 L2 out begin
            $q
        end
        out
    end
    if vaddsym != :out
        pushfirst!(q.args, :($vvaddsym = VectorizationBase.vectorizable($vaddsym)))
    end
    q # PaddedMatrices.simplify_expr(q)
end
@generated function row_sum_prod_add!(
    out::PaddedMatrices.AbstractMutableFixedSizeVector{M,T},
    L1::AbstractLowerTriangularMatrix{M,T},
    L2::AbstractLowerTriangularMatrix{M,T},
    v::AbstractFixedSizeVector{M,T}
) where {M,T}
    row_sum_add_quote(M, T, :v)
end
@generated function row_sum_prod_add!(
    out::PaddedMatrices.AbstractMutableFixedSizeVector{M,T},
    L1::AbstractLowerTriangularMatrix{M,T},
    L2::AbstractLowerTriangularMatrix{M,T}
) where {M,T}
    row_sum_add_quote(M, T, :out)
end
function row_sum_prod_add(
    sp::StackPointer,
    L1::AbstractLowerTriangularMatrix{M,T},
    L2::AbstractLowerTriangularMatrix{M,T},
    v::AbstractFixedSizeVector{M,T}
) where {M,T}
    (sp,out) = PtrVector{M,T}(sp)
    row_sum_prod_add!(out, L1, L2, v1)
    sp, out
end

function unrolled_update(
    P::Int, W::Int, Wshift::Int, T::DataType;
    xsym::Symbol = :ptrx, xusym::Symbol = :ptrx, Lsym::Symbol = :ptrLdiag, Lusym::Symbol = :ptrLdiag, Ltrisym::Symbol = :ptrLtri, Lutrisym::Symbol = :ptrLtri,
    rem_loop::Bool = true#, xscalar::Bool = false, track_L::Bool = false, track_x::Bool = false
)
    size_T = sizeof(T)
    V = Vec{W,T}
    q = quote end
    for c ∈ (W-P):W-1
        p = W - c
        d_c = Symbol(:d_,c)
        x_c = Symbol(:x_,c)
        c_c = Symbol(:c_,c)
        s_c = Symbol(:s_,c)
        l_c_c = Symbol(:L_,c,:_,c)
        invl_c_c = Symbol(:invL_,c,:_,c)
        push!(q.args, Expr(:(=), x_c, :(VectorizationBase.load($xsym + $c*$size_T))))
        push!(q.args, Expr(:(=), l_c_c, :(VectorizationBase.load($Lsym + $c*$size_T))))
        push!(q.args, Expr(:(=), invl_c_c, :(Base.FastMath.div_fast(one($T), $l_c_c))))
        push!(q.args, Expr(:(=), d_c, macroexpand(Base, :(@fastmath sqrt($l_c_c*$l_c_c + $x_c*$x_c)))))
        push!(q.args, Expr(:(=), c_c, :(SIMDPirates.vbroadcast($V,Base.FastMath.mul_fast($d_c, $invl_c_c)))))
        push!(q.args, Expr(:(=), s_c, :(SIMDPirates.vbroadcast($V,Base.FastMath.mul_fast($x_c, $invl_c_c)))))
        push!(q.args, Expr(:call, :(VectorizationBase.store!), :($Lusym + $c*$size_T), d_c))
        # if track_x
        #     if xscalar
        #         ∂d∂x = Symbol(:∂d_,c,:∂x)
        #         push!(q.args, Expr(:(=), ∂d∂x, :(Base.FastMath.div_fast($x_c, $d_c))))
        #         ∂s∂x = Symbol(:∂s_,c,:∂x)
        #         ∂c∂x = Symbol(:∂c_,c,:∂x)
        #         push!(q.args, Expr(:(=), ∂c∂x, :(Base.FastMath.mul_fast($∂d∂x, $invl_c_c))))
        #         push!(q.args, Expr(:(=), ∂s∂x, invl_c_c))
        #         push!(q.args, :(VectorizationBase.store!($(Symbol(:∂,xsym)) + $c*$size_T), VectorizationBase.load($(Symbol(:∂,xsym,:∂a)) + $c*$size_T) * ∂d∂x))
        #     else
        #         throw("xscalar == $xscalar while track_x == $track_x is not yet supported.")
        #     end
        # end
        # if track_L
        #     ∂d∂Lcc = Symbol(:∂d∂L_,p,:_,p)
        #     push!(q.args, Expr(:(=), ∂d∂Lcc, :(Base.FastMath.div_fast($l_c_c, $d_c))))
        #     ∂invL∂L = Symbol(:∂invL_,p,:_,p,:∂L_,p,:_,p)
        #     push!(q.args, Expr(:(=), ∂invL∂L, :(Base.FastMath.mul_fast(-one($T), $invl_c_c, $invl_c_c))))
        #     ∂s∂Lcc = Symbol(:∂s_,c,:∂L_,c,:_,c)
        #     ∂c∂Lcc = Symbol(:∂c_,c,:∂L_,c,:_,c)
        #     push!(q.args, Expr(:(=), ∂s∂Lcc, :(Base.FastMath.mul_fast($x_c, $∂invL∂L))))
        #     push!(q.args, Expr(:(=), ∂c∂Lcc, :(Base.FastMath.add_fast(Base.FastMath.mul_fast($d_c, $∂invL∂L), Base.FastMath.mul_fast($∂d∂Lcc,$invl_c_c)))))
        # end
        if p > 1
            rem = p - 1
            mask = ~ VectorizationBase.mask(T, c+1)
            vL = Symbol(:vL_,c+1,:_,c)
            vLu = Symbol(:vLu_,c+1,:_,c)
            vx = Symbol(:vL_,c+1)
            push!(q.args, Expr(:(=), vL, :(SIMDPirates.vload($V,$Ltrisym + trioffset*$size_T, $mask))))
            push!(q.args, Expr(:(=), vx, :(SIMDPirates.vload($V,$xsym, $mask))))
            push!(q.args, Expr(:(=), vLu, :(SIMDPirates.vfdiv(SIMDPirates.vadd($vL, SIMDPirates.vmul($s_c, $vx)), $c_c))))
            push!(q.args, :(SIMDPirates.vstore!($Lutrisym + trioffset*$size_T, $vLu,$mask)))
            push!(q.args, :(SIMDPirates.vstore!($xusym, SIMDPirates.vsub(SIMDPirates.vmul($c_c,$vx), SIMDPirates.vmul($s_c, $vLu)), $mask)))
        end
        if rem_loop
            rem_quote = quote
                for rep ∈ 1:repetitions
                    vL_j_p = SIMDPirates.vload($V, $Ltrisym + $size_T*trioffset + $size_T*$W*rep)
                    vx_j = SIMDPirates.vload($V, $xsym + $size_T*$W*rep)
                    vLu_j_p = SIMDPirates.vfdiv(SIMDPirates.vadd(vL_j_p, SIMDPirates.vmul($s_c,vx_j)), $c_c)
                    SIMDPirates.vstore!($Lutrisym + $size_T*trioffset + $size_T*$W*rep, vLu_j_p)
                    SIMDPirates.vstore!($xsym + $size_T*$W*rep, SIMDPirates.vsub(SIMDPirates.vmul($c_c,vx_j), SIMDPirates.vmul($s_c,vLu_j_p)))
                end
            end
            push!(q.args, rem_quote)
        end
        if p > 1
            push!(q.args, :(trioffset += triincrement))
        else
            push!(q.args, :(trioffset += triincrement+$W))
        end
        push!(q.args, :(triincrement -= 1))
    end
    push!(q.args, :(repetitions -= 1; $xsym += $W * $size_T; $Lsym += $W*$size_T;))
    Lusym == Lsym || push!(q.args, :($Lusym += $W*$size_T))
    xusym == xsym || push!(q.args, :($xusym += $W*$size_T))
    q
end

function rank_one_updated_lower_triangle_quote(P, T; xsym = :x, xusym = :x, Lsym = :L, Lusym = :L)
    W, Wshift = VectorizationBase.pick_vector_width_shift(P-1,T) # Even with avx512, we'd want W=4 for P=5, or W=2 for P=3,etc.
    Wm1 = W - 1
    V = Vec{W,T}
    reps = P >>> Wshift
    rem = P & Wm1
    size_T = sizeof(T)
    xsymp = Symbol(:ptr, xsym)
    xusymp = Symbol(:ptr, xusym)
    Lsymp = Symbol(:ptr, Lsym,:diag)
    Lusymp = Symbol(:ptr, Lusym,:diag)
    Ltrisymp = Symbol(:ptr, Lsym,:tri)
    Lutrisymp = Symbol(:ptr, Lusym,:tri)
    # We loop over batches of 8
    initial_offset = rem == 0 ? 0 : W - rem
    q = quote
        triincrement = $(P - 2)
        trioffset = 0
        repetitions = $(reps + (rem > 0) - 1)
        $xsymp = pointer($xsym) - $size_T * $initial_offset
        $Lsymp = pointer($Lsym) - $size_T * $initial_offset
        $Ltrisymp = $Lsymp + $size_T * $(P-1)
    end
    if xsym != xusym
        push!(q.args, :($xusymp = pointer($xusym) - $size_T * $initial_offset))
    end
    if Lsym != Lusym
        push!(q.args, :($Lusymp = pointer($Lusym) - $size_T * $initial_offset; $Lutrisymp = $Lusymp + $size_T * $(P-1)))
    end
    if rem > 0
        push!(q.args, unrolled_update(rem, W, Wshift, T, xsym = xsymp, xusym = xusymp, Lsym = Lsymp, Lusym = Lusymp, Ltrisym = Ltrisymp, Lutrisym = Lutrisymp, rem_loop = reps > 0))
    end
    if reps == 1
        push!(q.args, unrolled_update(W, W, Wshift, T, xsym = xsymp, xusym = xusymp, Lsym = Lsymp, Lusym = Lusymp, Ltrisym = Ltrisymp, Lutrisym = Lutrisymp, rem_loop = false))
    elseif reps > 1
        loop_quote = quote
            for _ ∈ 1:$reps
                $(unrolled_update(W, W, Wshift, T, xsym = xsymp, xusym = xusymp, Lsym = Lsymp, Lusym = Lusymp, Ltrisym = Ltrisymp, Lutrisym = Lutrisymp, rem_loop = true))
            end
        end
        push!(q.args, loop_quote)
    end
    push!(q.args, nothing)
    q
end

@generated function rank_update!(
    Lu::AbstractMutableLowerTriangularMatrix{P,T},
    L::AbstractMutableLowerTriangularMatrix{P,T},
    a::Union{<:PaddedMatrices.AbstractMutableFixedSizeVector{P,T},T}
) where {T,P}
    quote
        x = FixedSizeVector{$P,$T}(undef)
        x .= a
        $(rank_one_updated_lower_triangle_quote(P,T,Lsym = :L, Lusym = :Lu))
        Lu
    end
end
@generated function rank_update!(
    sptr::StackPointer,
    Lu::AbstractMutableLowerTriangularMatrix{P,T},
    L::AbstractMutableLowerTriangularMatrix{P,T},
    a::Union{<:PaddedMatrices.AbstractMutableFixedSizeVector{P,T},T}
) where {T,P}
    quote
        x = PtrVector{$P,$T}(pointer(sptr,$T))
        x .= a
        $(rank_one_updated_lower_triangle_quote(P,T,Lsym = :L, Lusym = :Lu))
        Lu
    end
end
function rank_update(L::AbstractLowerTriangularMatrix{P,T}, a) where {P,T}
    Lu = MutableLowerTriangularMatrix{P,T}(undef)
    rank_update!(Lu, L, a)
end
@generated function rank_update(sptr::StackPointer, L::AbstractLowerTriangularMatrix{P,T,N}, a) where {P,T,N}
    N2 = PaddedMatrices.calc_padding(N)
    ptr_offset = VectorizationBase.align(N2*sizeof(T))
    quote
        Lu = PtrLowerTriangularMatrix{$P,$T,$N2}(pointer(sptr,$T))
#        rank_update!(sptr + $(VectorizationBase.align(N*sizeof(T))), Lu, L, a)
        sptr + $ptr_offset, rank_update!(Lu, L, a)
    end
end

function phi_at_b_quote(P, T, load_A = true, store_C = true, halve_diagonal = true)
    W = VectorizationBase.pick_vector_width(P,T)
    if P > 4W && store_C && load_A
        return quote
            increment = $P
            @inbounds @fastmath for c in 1:$P
                for r in 1:c-1
                    C[r,c] = zero($T)
                end
                C_c_c = A[c] * B[c]
                for d in 1:$P-c
                    C_c_c += A[increment + d] * B[increment + d]
                end
                C[c,c] = $(halve_diagonal ? :($(T(0.5))*C_c_c) : :C_c_c)
                incr_r = increment + $P - c
                for r in 1+c:$P
                    C_r_c = A[r] * B[increment+r-c]
                    # @show C_r_c
                    incr_r -= r
                    for i in r+1:$P
                        # @show r,c, (incr_r+i), (increment+i-c)
                        C_r_c += A[incr_r + i] * B[increment+i-c]
                    end
                    incr_r += $P
                    C[r,c] = C_r_c
                end
                increment += $P - c
            end
            C
        end
    end
    q = quote end
    increment = P
    for c in 1:P
        if store_C
            for r in 1:c-1
                push!(q.args, :($(Symbol(:C_,r,:_,c)) = zero($T)))
            end
        end
        C_c_c = Symbol(:C_,c,:_,c)
        load_A ? push!(q.args, :($C_c_c = A[$c] * B[$c])) : push!(q.args, :($C_c_c = $(Symbol(:A_,c,:_,c)) * B[$c]))
        if load_A
            for d in 1:P-c
                push!(q.args, :($C_c_c += A[$(increment + d)] * B[$(increment + d)]))
            end
        else
            for d in 1:P-c
                push!(q.args, :($C_c_c += $(Symbol(:A_,c+d,:_,c)) * B[$(increment + d)]))
            end
        end
        halve_diagonal && push!(q.args, :($C_c_c *= $(T(0.5))))
        incr_r = increment + P - c
        for r in 1+c:P
            if load_A
                push!(q.args, :($(Symbol(:C_,r,:_,c)) = A[$r] * B[$(increment+r-c)]))
                incr_r -= r
                for i in r+1:P
                    push!(q.args, :($(Symbol(:C_,r,:_,c)) += A[$(incr_r + i)] * B[$(increment+i-c)]))
                end
                incr_r += P
            else
                push!(q.args, :($(Symbol(:C_,r,:_,c)) = $(Symbol(:A_,r,:_,r)) * B[$(increment+r-c)]))
                for i in r+1:P
                    push!(q.args, :($(Symbol(:C_,r,:_,c)) += $(Symbol(:A_,i,:_,r)) * B[$(increment+i-c)]))
                end
            end
        end
        increment += P - c
        if store_C
            for r in 1:P
                push!(q.args, :(C[$r,$c] = $(Symbol(:C_,r,:_,c))))
            end
        end
    end
    q
end

@generated function PhiAtB!(
    C::PaddedMatrices.AbstractMutableFixedSizeMatrix{P,P,T},
    A::AbstractMutableLowerTriangularMatrix{P,T}, # original chol
    B::AbstractMutableLowerTriangularMatrix{P,T} # Adjoint
) where {P,T}
    q = phi_at_b_quote(P, T)
    quote
        @fastmath @inbounds begin
            $q
        end
        C
    end
end

function reverse_chol_grad_expr(P, T; store_S::Bool = true, halve_diagonal = true)
    q = quote end
    inc = P
    for c in 1:P
        push!(q.args, Expr(:(=), Symbol(:A_,c,:_,c), :(A[$c])))
        push!(q.args, Expr(:(=), Symbol(:invA_,c,:_,c), :(Base.FastMath.div_fast(one($T),$(Symbol(:A_,c,:_,c))))))
        for r in c+1:P
            inc += 1
            push!(q.args, Expr(:(=), Symbol(:A_,r,:_,c), :(A[$inc])))
        end
    end
    push!(q.args, phi_at_b_quote(P, T, false, false, false)) # We treat as a reflected Symmetric matrix.
    # Now, we must calc A' \ C / A = (C / A)' / A
    for _c_ in 0:P-1
        c = P - _c_
        for r in 1:P
            X_r_c = Symbol(:X_,c,:_,r) # we transpose X; hence Symbol(:X_,c,:_,r) instead of _r_c
            clt, rlt = minmax(r, c)
            C_r_c = Symbol(:C_,rlt,:_,clt)
            push!(q.args, :($X_r_c = $C_r_c))
            for j in c+1:P
                clt, rlt = minmax(r, j)
                X_r_j = Symbol(:X_,j,:_,r)
                A_j_c = Symbol(:A_,j,:_,c)
                push!(q.args, :($X_r_c -= $X_r_j * $A_j_c ))
            end
            push!(q.args, :($X_r_c *= $(Symbol(:invA_,c,:_,c))))
        end
    end
    # Second division.
    for _c_ in 0:P-1
        c = P - _c_
        for r in 1:c
            X_r_c = Symbol(:X_,r,:_,c)
            S_r_c = Symbol(:S_,c,:_,r)
            push!(q.args, :($S_r_c = $X_r_c))
            for j in c+1:P
                S_r_j = Symbol(:S_,j,:_,r)
                A_j_c = Symbol(:A_,j,:_,c)
                push!(q.args, :($S_r_c -= $S_r_j * $A_j_c ))
            end
            push!(q.args, :($S_r_c *= $(Symbol(:invA_,c,:_,c))))
        end
    end
    if store_S
        inc = P
        for c in 1:P
            if halve_diagonal
                push!(q.args, Expr(:(=), :(S[$c]), Expr(:call, :(*), T(0.5), Symbol(:S_,c,:_,c))))
            else
                push!(q.args, Expr(:(=), :(S[$c]), Symbol(:S_,c,:_,c)))
            end
            for r in c+1:P
                inc += 1
                push!(q.args, Expr(:(=), :(S[$inc]), Symbol(:S_,r,:_,c)))
            end
        end
    end
    q
end

@generated function reverse_cholesky_grad!(
    S::AbstractMutableLowerTriangularMatrix{P,T},
    A::AbstractMutableLowerTriangularMatrix{P,T}, # original chol
    B::AbstractMutableLowerTriangularMatrix{P,T} # Adjoint
) where {P,T}
# ) where {T,P}
    q = reverse_chol_grad_expr(P, T, halve_diagonal = false)
    quote
        @fastmath @inbounds begin
            $q
        end
        S
    end
end
function reverse_cholesky_grad(
    A::AbstractMutableLowerTriangularMatrix{P,T}, # original chol
    B::AbstractMutableLowerTriangularMatrix{P,T} # Adjoint
) where {P,T}
    S = MutableLowerTriangularMatrix{P,T}(undef)
    reverse_cholesky_grad!(S, A, B)
end
@generated function reverse_cholesky_grad(
    sptr::StackPointer,
    A::AbstractMutableLowerTriangularMatrix{P,T}, # original chol
    B::AbstractMutableLowerTriangularMatrix{P,T} # Adjoint
) where {P,T}
    Wm1 = VectorizationBase.pick_vector_width(T) - 1
    N = PaddedMatrices.calc_padding(binomial2(P+1))
    quote
        S = PtrLowerTriangularMatrix{$P,$T,$N}(pointer(sptr,$T))
        sptr + $(VectorizationBase.align(sizeof(T) * N)), spreverse_cholesky_grad!(S, A, B)
    end
end

struct RankUpdateAdjoint{P,T,L,A<:AbstractMutableLowerTriangularMatrix{P,T,L}}
    data::A
end

function ∂rank_update_quote(P, T; track_L::Bool, track_x::Bool, xscalar::Bool = false, store_S::Bool = false)
    (track_L || track_x) || return quote end
    q = reverse_chol_grad_expr(P, T, store_S = store_S)
    if track_L
        ### for L
        # if s isa Number
        #     sv = fill!(FixedSizeVector{8,Float64}(undef), s);
        # else
        #     sv = s
        # end
        # Lu = StructuredMatrices.rank_update(L1, sv)
        # StructuredMatrices.reverse_cholesky_grad!(S, Lu, Ladjoint)
        # LowerTriangular(Symmetric(Array(S),:L)*L1)
        for k ∈ 1:P
            for c ∈ 1:k
                assignment = k == c ? :(=) : :(+=)
                L_k_c = Symbol(:L_,k,:_,c)
                if k == c
                    push!(q.args, Expr(:(=), L_k_c, :(Linput[$k])))
                else
                    push!(q.args, Expr(:(=), L_k_c, :(Linput[$(lt_sub2ind_fast(P,k,c))])))
                end
                for r ∈ c:P
                    pL_r_c = Symbol(:∂L_,r,:_,c)
                    cs, rs = minmax(r, k)
                    S_r_k = Symbol(:S_,rs,:_,cs)
                    push!(q.args, Expr(assignment, pL_r_c, :($S_r_k * $L_k_c)))
                end
            end
        end
        inc = P
        for c ∈ 1:P
            push!(q.args, :(∂L[$c] = $(Symbol(:∂L_,c,:_,c))))
            for r ∈ c+1:P
                pL_r_c = Symbol(:∂L_,r,:_,c)
                inc += 1
                push!(q.args, :(∂L[$inc] = $pL_r_c))
            end
        end
    end
    if track_x
        ### for xscalar
        # x2p3 = fill!(FixedSizeVector{8,Float64}(undef), s);
        # Lu = StructuredMatrices.rank_update(L1, x2p3)
        # StructuredMatrices.reverse_cholesky_grad!(S, Lu, Ladjoint)
        # sumSsub = zero(eltype(S))
        # sumSdiag = zero(eltype(S))
        # N = size(S,2)
        # for n ∈ 1:N
        #     sumSdiag += S[n,n]
        # end
        # for c ∈ 1:N
        #     for r ∈ c+1:N
        #         sumSsub += S[r,c]
        #     end
        # end
        # s*(sumSdiag + 2sumSsub)
        W, Wshift = VectorizationBase.pick_vector_width_shift(P,T)
        Wm1 = W - 1
        indsumsub = 0
        if xscalar
            for c ∈ 1:P
                sumSdiag = Symbol(:sumdiag_, (c-1) & Wm1)
                assigment = ( ( (c-1) >>> Wshift ) == 0 ) ? :(=) : :(+=)
                push!(q.args, Expr(assigment, sumSdiag, Symbol(:S_,c,:_,c)))
                for r ∈ c+1:P
                    sumSsub = Symbol(:sumsub_, indsumsub & Wm1)
                    assigment = ( ( indsumsub >>> Wshift ) == 0 ) ? :(=) : :(+=)
                    push!(q.args, Expr(assigment, sumSsub, Symbol(:S_,r,:_,c)))
#                    push!(q.args, :(@show $(string(assigment)), $sumSsub, $(Symbol(:S_,r,:_,c))))
                    indsumsub += 1
                end
            end
            WLd = min(Wm1, P-1)
            WLs = min(Wm1, indsumsub-1)
            push!(q.args, Expr(:(=), :sumSdiag, Expr(:call, :(+), [Symbol(:sumdiag_,i) for i ∈ 0:WLd]...)))
            push!(q.args, Expr(:(=), :sumSsub, Expr(:call, :(+), [Symbol(:sumsub_,i) for i ∈ 0:WLs]...)))
#            push!(q.args, :(@show sumSdiag, sumSsub))
            push!(q.args, Expr(:(=), :∂x, :(xinput * (sumSdiag + $(T(2))*sumSsub))))
        ### for xvector
        # Lu = StructuredMatrices.rank_update(L1, s)
        # StructuredMatrices.reverse_cholesky_grad!(S, Lu, Ladjoint)
        # grad = similar(s)
        # for i ∈ eachindex(grad)
        #     g = zero(eltype(grad))
        #     for j ∈ 1:i-1
        #         g += S[i,j] * s[j]
        #     end
        #     g += S[i,i] * s[i]
        #     for j ∈ i+1:length(grad)
        #         g += S[j,i] * s[j]
        #     end
        #     grad[i] = g
        # end
        # grad
        else
            for k in 1:P
                assignment = k == 1 ? :(=) : :(+=)
                push!(q.args, :(x_k = xinput[$k]))
                for i in 1:P
                    cs, rs = minmax(i,k)
                    S_i_k = Symbol(:S_,rs,:_,cs)
                    push!(q.args, Expr(assignment, Symbol(:x_,i), :($S_i_k * x_k)))
                end
            end
            for i in 1:P
                push!(q.args, :(∂x[$i] = $(Symbol(:x_,i))))
            end
        end
    end
    q
end

@generated function ∂rank_update(
    sptr::StackPointer,
    A::AbstractMutableLowerTriangularMatrix{P,T}, # original output argument
    B::AbstractMutableLowerTriangularMatrix{P,T}, # adjoint
    Linput::AbstractMutableLowerTriangularMatrix{P,T}, # original input argument
    xinput::T # original input argument
#) where {P,T}
) where {T,P}
    q = ∂rank_update_quote(P, T, track_L = true, track_x = true, xscalar = true)
    PL = PaddedMatrices.calc_padding(binomial2(P+1), T)
    sptroffset = VectorizationBase.align(sizeof(T) * PL)
    quote
        ∂L = PtrLowerTriangularMatrix{$P,$T,$PL}(pointer(sptr,$T))
        #∂x = FixedSizeVector{$P,$T}(undef)
        @inbounds @fastmath begin
            $q
        end
        sptr + $sptroffset, (∂L, ∂x)
    end
end
@generated function ∂rank_update(
    A::AbstractMutableLowerTriangularMatrix{P,T}, # original output argument
    B::AbstractMutableLowerTriangularMatrix{P,T}, # adjoint
    Linput::AbstractMutableLowerTriangularMatrix{P,T}, # original input argument
    xinput::T # original input argument
#) where {P,T}
) where {T,P}
    q = ∂rank_update_quote(P, T, track_L = true, track_x = true, xscalar = true)
    quote
        ∂L = MutableLowerTriangularMatrix{$P,$T}(undef)
        #∂x = FixedSizeVector{$P,$T}(undef)
        @inbounds @fastmath begin
            $q
        end
        ∂L, ∂x
    end
end

@generated function ∂rank_update(
    sptr::StackPointer,
    A::AbstractMutableLowerTriangularMatrix{P,T}, # original output argument
    B::AbstractMutableLowerTriangularMatrix{P,T}, # adjoint
    Linput::AbstractMutableLowerTriangularMatrix{P,T}, # original input argument
    xinput::PaddedMatrices.AbstractMutableFixedSizeVector{P,T} # original input argument
) where {P,T}
    q = ∂rank_update_quote(P, T, track_L = true, track_x = true, xscalar = false)
    PL = PaddedMatrices.calc_padding(binomial2(P+1), T)
    sptroffset1 = VectorizationBase.align(sizeof(T) * PL)
    sptroffset2 = VectorizationBase.align(sptroffset1 + P*sizeof(T))
    quote
        _sptr = pointer(sptr,$T)
        ∂L = PtrLowerTriangularMatrix{$P,$T,$PL}(_sptr)
        ∂x = PtrVector{$P,$T}(_sptr + $sptroffset1)
        @inbounds @fastmath begin
            $q
        end
        sptr + $sptroffset2, (∂L, ∂x)
    end
end
@generated function ∂rank_update(
    A::AbstractMutableLowerTriangularMatrix{P,T}, # original output argument
    B::AbstractMutableLowerTriangularMatrix{P,T}, # adjoint
    Linput::AbstractMutableLowerTriangularMatrix{P,T}, # original input argument
    xinput::PaddedMatrices.AbstractMutableFixedSizeVector{P,T} # original input argument
) where {P,T}
# ) where {T,P}
    q = ∂rank_update_quote(P, T, track_L = true, track_x = true, xscalar = false)#, store_S = true)
    quote
        ∂L = MutableLowerTriangularMatrix{$P,$T}(undef)
        ∂x = FixedSizeVector{$P,$T}(undef)
        # S = MutableLowerTriangularMatrix{$P,$T}(undef)
        @inbounds @fastmath begin
            $q
        end
        # S, ∂L, ∂x
        ∂L, ∂x
    end
end

@generated function ∂rank_update(
    sptr::StackPointer,
    A::AbstractMutableLowerTriangularMatrix{P,T}, # original output argument
    B::AbstractMutableLowerTriangularMatrix{P,T}, # adjoint
    xinput::T  # original input argument
) where {P,T}
    q = ∂rank_update_quote(P, T, track_L = false, track_x = true, xscalar = true)
    quote
        @inbounds @fastmath begin
            $q
        end
        sptr, ∂x
    end
end
@generated function ∂rank_update(
    A::AbstractMutableLowerTriangularMatrix{P,T}, # original output argument
    B::AbstractMutableLowerTriangularMatrix{P,T}, # adjoint
    xinput::T # original input argument
) where {P,T}
    q = ∂rank_update_quote(P, T, track_L = false, track_x = true, xscalar = true)
    quote
        @inbounds @fastmath begin
            $q
        end
        ∂x
    end
end

@generated function ∂rank_update(
    sptr::StackPointer,
    A::AbstractMutableLowerTriangularMatrix{P,T}, # original output argument
    B::AbstractMutableLowerTriangularMatrix{P,T}, # adjoint
    xinput::PaddedMatrices.AbstractMutableFixedSizeVector{P,T} # original input argument
) where {P,T}
    q = ∂rank_update_quote(P, T, track_L = false, track_x = true, xscalar = false)
    sptroffset1 = VectorizationBase.align(P*sizeof(T))
    quote
        ∂x = PtrVector{$P,$T}(pointer(sptr,$T))
        @inbounds @fastmath begin
            $q
        end
        sptr + $sptroffset1, ∂x
    end
end
@generated function ∂rank_update(
    A::AbstractMutableLowerTriangularMatrix{P,T}, # original output argument
    B::AbstractMutableLowerTriangularMatrix{P,T}, # adjoint
    xinput::PaddedMatrices.AbstractMutableFixedSizeVector{P,T} # original input argument
) where {P,T}
    q = ∂rank_update_quote(P, T, track_L = false, track_x = true, xscalar = false)
    quote
        ∂x = FixedSizeVector{$P,$T}(undef)
        @inbounds @fastmath begin
            $q
        end
        ∂x
    end
end
@generated function ∂rank_update(
    sptr::StackPointer,
    A::AbstractMutableLowerTriangularMatrix{P,T}, # original output argument
    B::AbstractMutableLowerTriangularMatrix{P,T}, # adjoint
    Linput::AbstractMutableLowerTriangularMatrix{P,T} # original input argument
) where {P,T}
    q = ∂rank_update_quote(P, T, track_L = true, track_x = false)
    PL = PaddedMatrices.calc_padding(binomial2(P+1), T)
    quote
        ∂L = PtrLowerTriangularMatrix{$P,$T,$PL}(pointer(sptr,$T))
        @inbounds @fastmath begin
            $q
        end
        sptr + $(VectorizationBase.align(PL*sizeof(T))), ∂L
    end
end
@generated function ∂rank_update(
    A::AbstractMutableLowerTriangularMatrix{P,T}, # original output argument
    B::AbstractMutableLowerTriangularMatrix{P,T}, # adjoint
    Linput::AbstractMutableLowerTriangularMatrix{P,T} # original input argument
) where {P,T}
    q = ∂rank_update_quote(P, T, track_L = true, track_x = false)
    quote
        ∂L = MutableLowerTriangularMatrix{$P,$T}(undef)
        @inbounds @fastmath begin
            $q
        end
        ∂L
    end
end

# PaddedMatrices.@support_stack_pointer ∂rank_update!

#@generated function rank_update()
#
#end


text:
# This file was generated, do not modify it. # hide
function MLJ.fit(model::KNNRidgeBlend, verbosity::Int, X, y)
    Xs = source(X)
    ys = source(y)
    hot = machine(OneHotEncoder(), Xs)
    W = transform(hot, Xs)
    z = log(ys)
    ridge_model = model.ridge_model
    knn_model = model.knn_model
    ridge = machine(ridge_model, W, z)
    knn = machine(knn_model, W, z)
    # and finally
    ẑ = model.knn_weight * predict(knn, W) + (1.0 - model.knn_weight) * predict(ridge, W)
    ŷ = exp(ẑ)

    mach = machine(Deterministic(), Xs, ys; predict=ŷ)
    fit!(mach, verbosity=verbosity - 1)
    return mach()
end

ext: jl


text:
struct OrderedStronglyConvexError <: Exception
end

function Base.showerror(io::IO, e::OrderedStronglyConvexError)
    print(io, "Points are not ordered or do not represent a strongly convex set")
end


text:
include("./utilities.jl")
include("../../src/CGPO.jl")

# ==============================================================================
# Floats on the exact matrix.

function FloatsOnExactMatrix()
    N = 64
    function PerformForRho(rho)
        Eigens = EigDistribution(N, rho=rho, eigen_min=1e-5)
        A, _ = TinyIntervalTestMatrices(
            Eigens,
            scale=N*1e-5(maximum(Eigens)/minimum(Eigens))
        )
        b = ones(size(A, 1))
        RelErrors = PerformCGFor(A, b, exact=false, epsilon=1e-10)
    return RelErrors end

    RelErrors = PerformForRho(0.8)
    fig = plot(
        RelErrors,
        yaxis=:log10,
        xlabel="iterations: k",
        ylabel="\$\\frac{\\Vert e_k \\Vert_A}{\\Vert e_0\\Vert_A}\$",
        label="\$\\rho=0.8\$",
        left_margin = 10Plots.mm,
        size=(750, 500), dpi=300,
    )
    RelErrors = PerformForRho(0.9)
    plot!(
        fig,
        RelErrors,
        label="\$\\rho=0.9\$",
        linestyle=:dot
    )
    RelErrors = PerformForRho(0.7)
    plot!(
        fig,
        RelErrors,
        label="\$\\rho=0.7\$",
        linestyle=:dashdot
    )
    RelErrors = PerformForRho(1)
    plot!(
        fig,
        RelErrors,
        label="\$\\rho=1\$",
        linestyle=:dash
    )
    title!(fig, "Float CG on Exact Matrix")
    fig |> display
    savefig(fig, "$(@__DIR__)/plots/floats_cg_on_exact_matrix")
return nothing end

FloatsOnExactMatrix()

# ==============================================================================
# Exact on the Smeared out matrix.

function ExactOnSmearedMatrix()
    N = 64
    function PerformForRho(rho)
        Eigens = EigDistribution(N, rho=rho, eigen_min=1e-5)
        ConditionNumb = maximum(Eigens)/minimum(Eigens)
        _, A = TinyIntervalTestMatrices(
            Eigens,
            scale=2e-5ConditionNumb,
            m=100
        )
        b = ones(size(A, 1))
        RelErrors = PerformCGFor(A, b, exact=true, epsilon=1e-10)
    return RelErrors end

    RelErrors = PerformForRho(0.8)
    fig = plot(
        RelErrors,
        yaxis=:log10,
        xlabel="iterations: k",
        ylabel="\$\\frac{\\Vert e_k \\Vert_A}{\\Vert e_0\\Vert_A}\$",
        label="\$\\rho=0.8\$",
        left_margin = 10Plots.mm,
        size=(750, 500), dpi=300,
    )
    RelErrors = PerformForRho(0.9)
    plot!(
        fig,
        RelErrors,
        label="\$\\rho=0.9\$",
        linestyle=:dot
    )
    RelErrors = PerformForRho(0.7)
    plot!(
        fig,
        RelErrors,
        label="\$\\rho=0.7\$",
        linestyle=:dashdot
    )
    RelErrors = PerformForRho(1)
    plot!(
        fig,
        RelErrors,
        label="\$\\rho=1\$",
        linestyle=:dash
    )
    title!(fig, "Exact CG on Smeared Matrix")
    savefig(fig, "$(@__DIR__)/plots/exact_cg_on_smeared_matrix")
    fig |> display
return nothing end

ExactOnSmearedMatrix()


text:
include(joinpath(dirname(@__FILE__), "..", "src", "ImmutableNestedSetModel.jl"))
using ImmutableNestedSetModel
#using Base.Test

# write your own tests here
#@test 1 == 2

#=
table = ATree()

function pl()
  for node = table
    println(repeat("-", parentscount(node)) * repr(node))
  end
  println("")
end

pl()

append!(table, Node("--"))
append!(table, Node("-- 2"))
pl()

=#
#=
r = [
  Node("ELECTRONICS"),
  Node("TELEVISIONS"),
  Node("TUBE"),
  Node("LCD"),
  Node("PLASMA"),
  Node("PORTABLE ELECTRONICS"),
  Node("MP3 PLAYERS"),
  Node("FLASH"),
  Node("CD PLAYERS"),
  Node("2 WAY RADIOS"),
]

push!(r[1], r[2])
push!(r[2], r[3])
push!(r[2], r[4])
push!(r[2], r[5])
push!(r[1], r[6])
push!(r[6], r[7])
push!(r[7], r[8])
push!(r[6], r[9])
push!(r[6], r[10])

push!(tree.root, r[1])
prepare(tree)
=#
t = Tree()


text:
include("Thermo/Thermo.jl")
using .Thermo.UniversalConstants ,
      .Thermo.Compounds

testSize = 50
testResults=Array{Bool,1}(undef, testSize)
testErrors=Array{String,1}(undef, testSize)
testCounter = 0

println("")
println("-----------------------------------------------------------")
println("Iniciating tests:")
println("-----------------------------------------------------------")
println("")
testName= "UniversalConstants.getUniversalConstant"
testCounter+=1
try
    R = UniversalConstants.getUniversalConstant("molarGasConstant")
    testResults[testCounter]=true
    println(string("Test"," N°",testCounter,", ",testName,": Passed"))
catch error
    testResults[testCounter]=false
    println("******************************")
    println(string("Test"," N°",testCounter,", ",testName,": Failed"))
    testErrors[testCounter] = string(error)
    println(testErrors[testCounter])
    println("******************************")
end

testName= "UniversalConstants.getUniversalConstantList"
testCounter+=1
try
    R = UniversalConstants.getUniversalConstantList()
    testResults[testCounter]=true
    println(string("Test"," N°",testCounter,", ",testName,": Passed"))
catch error;
    testResults[testCounter]=false
    println("******************************")
    println(string("Test"," N°",testCounter,", ",testName,": Failed"))
    testErrors[testCounter] = string(error)
    println(testErrors[testCounter])
    println("******************************")
end

testName= "Compounds.getConstPropList"
testCounter+=1
try
    R = Compounds.getConstPropList()
    testResults[testCounter]=true
    println(string("Test"," N°",testCounter,", ",testName,": Passed"))
catch error;
    testResults[testCounter]=false
    println("******************************")
    println(string("Test"," N°",testCounter,", ",testName,": Failed"))
    testErrors[testCounter] = string(error)
    println(testErrors[testCounter])
    println("******************************")
end

testName= "Compounds.newCompound"
testCounter+=1
try
    x1 = Dict{String,String}(
    "ID"=>"compound1",
    "iupacName"=>"sdasd",
    "casRegistryNumber"=>"sdasd112312")

    x2 = Dict{String,Float64}(
        "molecularWeight"=>18.0,
        "criticalTemperature"=>647,
        "criticalPressure"=>220000)

    compound1 = Compound(x1)
    compound2 = Compound(x1,x2,ID="compound2")
    compound3 = Compound(x1,ID="compound3")
    #println(compound1.ID)
    #println(compound2.ID)
    #compoundlist1 = Compounds.newCompoundList([compound1,compound2],["11","22"]) #not valid anymore
    Compoundsz = [compound1,compound2,compound3]
    IDs = map(x->x.StringProperties["ID"],Compoundsz)
    println(IDs)
    compoundlist2 = Compounds.CompoundList(Compoundsz)
    #println(Compounds.GetCompoundList(compoundlist1))
   # println(Compounds.GetNumCompounds(compoundlist2))
    #println(compound1.StringProperties["ID"])
    #println(comp
    println(compoundlist2["molecularWeight"])
    #c2 = Compounds.getCompound(c1,1)

    testResults[testCounter]=true
    println(string("Test"," N°",testCounter,", ",testName,": Passed"))
catch error;
    testResults[testCounter]=false
  println("******************************")
    println(string("Test"," N°",testCounter,", ",testName,": Failed"))
    testErrors[testCounter] = string(error)
    println(testErrors[testCounter])
    println("******************************")
end

mutable struct po
    x::Real
    y::Real
end

p1 = po(1,1)
p1.x = 40

passedResults = length(filter(x ->  x== true, testResults)) #usage of filter function
coverage = 100*passedResults/testCounter

println("")
println("-----------------------------------------------------------")
println(string("Tests results:"))
println(string("Total executed tests = ",testCounter))
println(string("Coverage = ",round(coverage,digits=2),"%"))
println("-----------------------------------------------------------")
println("")

#x = newCompound()


text:
"""
	rounded_star(r, d, num, N)

Return a `ShapeParams` object containing the shape parametrized by
\$(x(θ),y(θ)) = (r + d*cos(θ*num))*(cos(θ),sin(θ))\$ with 2`N` nodes.
"""
function rounded_star(r, d, num, N)
    t = Float64[pi*j/N for j=0:(2*N-1)]
    Rt = r .+ d*cos.(num*t)
    dRt = (-d*num)*sin.(num*t)
    ft = [Rt.*cos.(t) Rt.*sin.(t)]
    dft = [dRt.*cos.(t) - Rt.*sin.(t)     dRt.*sin.(t) + Rt.*cos.(t)]
    return ShapeParams(t, ft, dft)
end

"""
	squircle(r, N)

Return a `ShapeParams` object containing the shape parametrized by
\$x(θ)^4 + y(θ)^4 = r^4\$ with 2`N` nodes.
"""
function squircle(r, N)
    #substituting x = rho*cost, y = rho*sint in x^4+y^4=R^4
    t = Float64[pi*j/N for j=0:(2*N-1)]
    cost = cos.(t); cos4t = cos.(4*t);
    sint = sin.(t); sin4t = sin.(4*t);
    rho = (sqrt(2)*r)./((3 .+ cos4t).^0.25)
    ft = [rho.*cost rho.*sint]
    drho = rho.*sin4t./(3 .+ cos4t)
    dx = drho.*cost - rho.*sint
    dy = drho.*sint + rho.*cost
    dft = [dx dy]
    return ShapeParams(t, ft, dft)
end

"""
	ellipse(r1, r2, N)

Return a `ShapeParams` object containing the shape parametrized by
`(x/r1)^2 + (y/r2)^2 = 1` with 2`N` nodes.
"""
function ellipse(r1, r2, N)
	t = Float64[pi*j/N for j=0:(2*N-1)]
    ft = [r1*cos.(t) r2*sin.(t)]
    dft = [(-r1)*sin.(t) r2*cos.(t)]
    return ShapeParams(t, ft, dft)
end

"""
	square_grid(a::Integer, d)

Return `centers`, an `(a^2,2)` array containing the points on an `a` by `a`
grid of points distanced `d`.
"""
function square_grid(a::Integer, d)
    offsetx = -0.5*(a-1)
    offsety = -0.5*(a-1)
    centers = d*[mod.(0:a^2-1, a) .+ offsetx   div.(0:a^2-1, a) .+ offsety]
end

"""
	rect_grid(a::Integer, b::Integer, dx, dy)

Return `centers`, an `(a*b,2)` array containing the points spanned by `a` points
distanced `dx` and `b` points distanced `dy`, in the x and y directions,
respectively.
"""
function rect_grid(a::Integer, b::Integer, dx, dy)
    offsetx = -0.5*(a-1)
    offsety = -0.5*(b-1)
	xpoints = dx*((0:a-1) .+ offsetx)
	ypoints = dy*((0:b-1) .+ offsety)

    centers = [repeat(xpoints, outer=[b]) 	repeat(ypoints, inner=[a])]
end

"""
	hex_grid(a::Integer, rows::Integer, d; minus1 = false)

Return `centers`, an `(M,2)` array containing points on a hexagonal lattice
with horizontal rows, with `a` points distanced `d` in each row and `rows` rows.
If `minus1` is true, the last point in every odd row is omitted.
"""
function hex_grid(a::Integer, rows::Integer, d; minus1 = false)
	h = d*sqrt(0.75) #row height
	M = minus1 ? a*rows - div(rows,2) : a*rows
	centers = Array{Float64}(undef, M, 2)
	ind = 1
	for r = 0:rows-1
		if mod(r,2) == 0
			for i = 1:a
				centers[ind,1] = i*d
				centers[ind,2] = (r - 1)*h
				ind += 1
			end
		else
			for i = 1:a-1
				centers[ind,1] = (i + 0.5)*d
				centers[ind,2] = (r - 1)*h
				ind += 1
			end
			if !minus1
				centers[ind,1] = (a + 0.5)*d
				centers[ind,2] = (r - 1)*h
				ind += 1
			end
		end
	end
	offset = mean(centers, dims=1)
    centers .-= offset
end

"""
	randpoints(M, dmin, width, height; failures = 100)

Return `centers`, an `(M,2)` array containing `M` points distanced at least
`dmin` in a `width` by `height` box. Fails `failures` times successively before
giving up.
"""
function randpoints(M, dmin, width, height; failures = 100)
    x_res = [rand(Float64)*width]
    y_res = [rand(Float64)*height]
    fail = 0
    dmin2 = dmin^2

    while fail < failures && length(x_res) < M
        accepted = true
        x_try = rand(Float64)*width
        y_try = rand(Float64)*height

        for i = 1:length(x_res)
            dist2 = (x_res[i] - x_try)^2 + (y_res[i] - y_try)^2
            if dist2 <= dmin2
                accepted = false
                break
            end
        end
        if accepted
            push!(x_res,x_try)
            push!(y_res,y_try)
            fail = 0
        else
            fail = fail + 1
        end
    end
    if length(x_res) < M
        error("""randpoints: Could not generate random points.
                Try increasing domain, decreasing minimum distance,
                or increasing subsequent failures allowed""")
    end
    centers = [x_res y_res]
end

"""
	randpoints(M, dmin, width, height, points; failures = 100)

Same as randpoints(M, dmin, width, height; failures = 100) but also
requires centers to be distanced at least `dmin` from `points`.
"""
function randpoints(M, dmin, width, height, points; failures = 100)
    dmin2 = dmin^2
    for i = 2:size(points,1), j = 1:i-1
        dist2 = sum(x -> x^2, points[i,:] - points[j,:])
        dist2 <= dmin2 && error("randpoints: given points have distance <= dmin")
    end

    x_res = Array{Float64}(undef, 0)
    y_res = Array{Float64}(undef, 0)
    fail = 0
    while fail < failures && length(x_res) < M
        accepted = true
        x_try = rand(Float64)*width
        y_try = rand(Float64)*height

        for i = 1:size(points,1)
            dist2 = (points[i,1] - x_try)^2 + (points[i,2] - y_try)^2
            if dist2 <= dmin2
                accepted = false
                break
            end
        end
        accepted && for i = 1:length(x_res)
            dist2 = (x_res[i] - x_try)^2 + (y_res[i] - y_try)^2
            if dist2 <= dmin2
                accepted = false
                break
            end
        end
        if accepted
            push!(x_res,x_try)
            push!(y_res,y_try)
            fail = 0
        else
            fail = fail + 1
        end
    end
    if length(x_res) < M
        error("""randpoints: Could not generate random points.
                Try increasing domain, decreasing minimum distance,
                or increasing subsequent failures allowed""")
    end
    centers = [x_res y_res]
end

"""
    verify_min_distance(shapes, centers::Array{Float64,2}, ids)
    verify_min_distance(sp::ScatteringProblem)
Returns `true` if the shapes placed at `centers` are properly distanced
(non-intersecting scattering disks).
"""
function verify_min_distance(shapes, centers::Array{Float64,2}, ids)
    Ns = size(centers,1)
    for ic = 1:Ns, ic2 = ic+1:Ns
        d = sqrt((centers[ic,1] - centers[ic2,1])^2 +
                (centers[ic,2] - centers[ic2,2])^2)
        d > shapes[ids[ic]].R + shapes[ids[ic2]].R || (return false)
    end
    true
end

"""
    verify_min_distance(shapes, centers::Array{Float64,2}, ids, points::Array{Float64,2})
    verify_min_distance(sp::ScatteringProblem, points)
Returns `true` if the shapes placed at `centers` are properly distanced
(non-intersecting scattering disks), and all `points` are outside the scattering
disks.
"""
function verify_min_distance(shapes, centers::Array{Float64,2}, ids, points::Array{Float64,2})
    Ns = size(centers,1)
    for ic = 1:Ns
        for ip = 1:size(points,1)
            d = sqrt((points[ip,1] - centers[ic,1])^2 +
                    (points[ip,2] - centers[ic,2])^2)
            d > shapes[ids[ic]].R || (return false)
        end
        for ic2 = ic+1:Ns
            d = sqrt((centers[ic,1] - centers[ic2,1])^2 +
                    (centers[ic,2] - centers[ic2,2])^2)
            d > (shapes[ids[ic]].R + shapes[ids[ic2]].R) || (return false)
        end
    end
    true
end

verify_min_distance(sp::ScatteringProblem, points) =
	verify_min_distance(sp.shapes, sp.centers, sp.ids, points)
verify_min_distance(sp::ScatteringProblem) =
	verify_min_distance(sp.shapes, sp.centers, sp.ids)

"""
    luneburg_grid(R_lens, N_cells, er; levels = 0, TM = true) -> centers, ids, rs

Returns the coordinates and radii of the circular inclusions in a Luneburg lens
device of radius `R_lens` with `N_cells` unit cells across its diameter. Radii
are determined by averaging over cell permittivity, assuming air outside and
relative permittivity `er` in the rods, and depends on incident field polarization
(TM/TE with respect to z-axis).
If `levels` == 0, groups identical radii together, such that rs[ids[n]]
is the radius of the rod centered at `(center[n,1],center[n,2])`. Otherwise
quantizes the radii to uniformly spaced levels.
"""
function luneburg_grid(R_lens, N_cells, er; levels = 0, TM = true)
    r_cell = Array{Float64}(undef, N_cells^2)
    centers = Array{Float64}(undef, N_cells^2, 2)
    flag_outside = Array{Bool}(undef, N_cells^2)
    a = 2*R_lens/N_cells #cell dimension
    for ix = 1:N_cells, iy = 1:N_cells
        ind = (ix-1)*N_cells + iy
        centers[ind,1] = a*(ix - N_cells/2.0 - 0.5)
        centers[ind,2] = a*(iy - N_cells/2.0 - 0.5)
        rho2 = sum(abs2,centers[ind,:])
        if rho2 > R_lens^2 #for now
            flag_outside[ind] = false
        else
            flag_outside[ind] = true
            r_cell[ind] = getLuneburgRadius(rho2, R_lens, a, er, TM)
        end
    end
    r_cell = r_cell[flag_outside]
    centers = centers[flag_outside,:]

    if levels > 0
        #uniform quantization
        rs = collect(range(minimum(r_cell), stop=maximum(r_cell), length=levels))
        ddd = (r_cell - rs[1])/(rs[2] - rs[1])
        ids = convert(Array{Int,1},round(ddd)) + 1
    else
        #only store unique vals
        ids,rs = uniqueind(r_cell)
    end
    return centers, ids, rs
end

function getLuneburgRadius(rho2, R_lens, a, er, TM)
    er_host = 1
    n_0 = 1
    n_eff = n_0*sqrt(2-rho2/R_lens^2)
    xx = a^2*(er_host-n_eff^2)/(pi*(er_host-er))
    if TM
        r = sqrt(xx)
    else #TE
        r = sqrt(xx*(er_host + er)/(er_host + n_eff^2))
    end
    return r
end


text:
using ACE, StaticArrays, BenchmarkTools

using Zygote: gradient

Nprop = 2
wL = 1.5
deg = 10
ord = 3

Bsel = SparseBasis(; maxorder = ord, p = 1, default_maxdeg = deg,
                        weight = Dict(:n => 1.0, :l => wL))
B1p = ACE.Utils.RnYlm_1pbasis(maxdeg = deg,
                                 maxL = ceil(Int, deg / wL),
                                 Bsel = Bsel)

φ = ACE.Invariant()
bsis = ACE.SymmetricBasis(φ, B1p, ACE.O3(), Bsel)

#create a multiple property model
W = rand(Nprop, length(bsis))
c = [SVector{size(W)[1]}(W[:,i]) for i in 1:size(W)[2]]
LM = ACE.LinearACEModel(bsis, c, evaluator = :standard)

# generate a configuration
nX = 100
Xs = rand(PositionState{Float64}, B1p.bases[1], nX)
cfg = ACEConfig(Xs)

@info("evaluation")
@btime ACE.evaluate(LM, cfg);

@info("der parameters")
@btime ACE.grad_params(LM, cfg);
@info("calculate gradient")
@btime ACE.grad_config(LM, cfg);
@info("Adjoint gradient")
@btime ACE._rrule_evaluate(ones(Nprop), LM, cfg);
@info("der grad according to params")
@btime ACE.grad_params_config(LM, cfg)
@info("Adjoint der grad according to params")
adj = [ACE.DState(rr=rand(SVector{3, Float64})) for _ = 1:length(cfg)]
@btime ACE.adjoint_EVAL_D(LM, LM.evaluator, cfg, adj)

@info("Zygote calls")

site_energy(m,x) = sum(ACE.val.(ACE.evaluate(m, x)))
@info("sum over properties evaluate")
@btime site_energy(LM, cfg)
@info("der parameters")
@btime gradient(x->site_energy(x,cfg), LM)
@info("calculate gradient")
@btime gradient(x->site_energy(LM,x), cfg)
@info("der params of gradient, sum for adjoint")
@btime gradient(m->sum(sum(gradient(x->site_energy(m,x), cfg)[1]).rr), LM)

@info("define a loss with gradients")
FS = props -> sum( (1 .+ ACE.val.(props).^2).^0.5 )
sqr(x) = x.rr .^ 2

loss(m, x) = (FS(ACE.evaluate(m,x)))^2 + sum(sum(sqr.(gradient(tx->FS(ACE.evaluate(m,tx)), x)[1])))
lossE(m, x) = (FS(ACE.evaluate(m,x)))^2

@info("evaluate loss full loss")
@btime loss(LM, cfg)

@info("only energies")
@btime lossE(LM, cfg)

@info("der full loss")
@btime gradient(m->loss(m,cfg), LM)

@info("der only energies")
@btime gradient(m->lossE(m,cfg), LM);


text:
#
println("2) Test example for One_Particle_Scaler operator.")

configs = [Configuration("1s^2 2s"), Configuration("1s^2 2p")]

# Generate a list of relativistic configurations and determine an ordered list of subshells for these configurations
relconfList = ConfigurationR[]
for  conf in configs    wa = Basics.generate("configuration list: relativistic", conf);     append!( relconfList, wa)   end
for  i = 1:length(relconfList)    println(">> include ", relconfList[i])    end
subshellList = Basics.generate("subshells: ordered list for relativistic configurations", relconfList)
Defaults.setDefaults("relativistic subshell list", subshellList; printout=false)

# Generate the relativistic CSF's for the given subshell list
csfList = CsfR[]
for  relconf in relconfList
    newCsfs = Basics.generate("CSF list: from single ConfigurationR", relconf, subshellList);   append!( csfList, newCsfs)
end

# Calculate angular coefficients for a scalar one- or two-particle operator
op = SpinAngular.OneParticleOperator(0, plus, true)
for  leftCsf in csfList
    for rightCsf in csfList
        coeffs = SpinAngular.computeCoefficients(op, leftCsf, rightCsf, subshellList)
        println(">> Angular coeffients for <$leftCsf || O^($(op.rank))  || $rightCsf> = \n $coeffs ")
        println("..................................................................................")
        println("")
     end
end


text:
using Test
using LinearAlgebra
using BQCESubroutine
using BQCESubroutine.Utilities
using YaoLocations: YaoLocations, plain, Locations, CtrlLocations, merge_locations
using BQCESubroutine: ctrl_offset
using YaoArrayRegister
using BenchmarkTools
using LoopVectorization
using BQCESubroutine: perm_kernel, BitRoutine, multi_broutine2x2!, subspace_mul_kernel, BitRoutine
using MatchCore

@smatch :(A = B) begin
    :($D = $D) => D
    :($B = $C) => (B, C)
    _ => nothing
end


text:
## test/hbdict.jl --- check if bond dict functions properly/as expected

@testset "HBondDict (2 element) (get/set/iterate)" begin
    ## declare a throwaway bonddict
    dict = HBondDict(2)             # check if the constructor works
    @test_throws ErrorException HBondDict(0)

    ## type and iteration
    @test dict isa HBondDict{4}
    @test length(dict) == 4
    @test size(dict) == (4,)
    keys_ = keys(dict)
    values_ = values(dict)

    # the values inside
    @test all(∈(1:4), keys_)
    @test all(==((0, 0)), values_)

    # getindex
    @test dict[1] == dict[2] == dict[3] == dict[4] == (0, 0)
    @test_throws BoundsError dict[5]

    # setindex
    @test_throws BoundsError dict[5] = (3, 4)
    dict[1] = (3, 4)
    @test (dict.interacting_residues[1] == 3 && dict.energies[1] == 4)
    dict[3] = (3, 4)
    @test (dict.interacting_residues[3] == 3 && dict.energies[3] == 4)

    ## donor/acceptor
    @test 3 in acceptorindices(dict)
    @test !(4 in acceptorindices(dict))
    @test 3 in donorindices(dict)
    @test !(4 in acceptorindices(dict))

    @testset "hbonded function" begin
        ## test vector containing 2 dicts
        vec = [HBondDict(2), HBondDict(2)]
        ## initialize the testable parts of dictionaries
        ## format is (index, energy)
        vec[1][1] = (2, -0.9)
        vec[1][3] = (3, -1.5)
        vec[2][3] = (1, -0.9)
        vec[2][2] = (3, -6.9)

        @test hbonded(vec, 1, 2)
        @test_throws BoundsError hbonded(vec, 3, 1)
        @test_throws BoundsError hbonded(vec, 2, 3)
    end
end

@testset "HBondDict (1 element) (get/set/iterate)" begin
    ## declare a throwaway bonddict
    dict = HBondDict(1)             # check if the constructor works

    ## type and iteration
    @test dict isa HBondDict{2}
    @test length(dict) == 2
    @test size(dict) == (2,)
    keys_ = keys(dict)
    values_ = values(dict)

    # the values inside
    @test all(∈(1:2), keys_)
    @test all(==((0, 0)), values_)

    # getindex
    @test dict[1] == dict[2] == (0, 0)
    @test_throws BoundsError dict[3]

    # setindex
    @test_throws BoundsError dict[3] = (3, 4)
    dict[1] = (3, 4)
    @test (dict.interacting_residues[1] == 3 && dict.energies[1] == 4)
    dict[2] = (3, 4)
    @test (dict.interacting_residues[2] == 3 && dict.energies[2] == 4)

    ## donor/acceptor
    @test 3 in acceptorindices(dict)
    @test !(4 in acceptorindices(dict))
    @test 3 in donorindices(dict)
    @test !(4 in acceptorindices(dict))

    @testset "hbonded function" begin
        vec = [HBondDict(1), HBondDict(1)]
        vec[2][2] = (1, -0.9)
        vec[1][1] = (2, -0.9)

        @test hbonded(vec, 1, 2)
        @test_throws BoundsError hbonded(vec, 3, 2)
        @test_throws BoundsError hbonded(vec, 1, 3)
    end
end

@testset "Record functions" begin
    ## FIRST AND FOREMOST: check if they are defined
    @test isdefined(ProtStructs, :recorddonor!)
    @test isdefined(ProtStructs, :recordacceptor!)

    # check for the M ∈ [1, 2, 3] cases
    @testset "Donor recording" begin
        for M in 1:5
            dict = HBondDict(M)
            # fill in with random energies and indices
            for k in 1:10
                # no chance that this is 0
                recorddonor!(dict, k, -(rand() + 0.1) * 10)
            end

            # all energies are negative
            @test let dnrs = donors(dict)
                all(dnrs) do (_, energy)
                    energy < 0
                end
            end

            # energies should be in increasing order (decreasing bond strength)
            @test let dnrs = donors(dict)
                Q = true; prev = -20.0
                for (_, energy) in dnrs
                    Q &= energy > prev
                    Q || break
                    prev = energy
                end
                Q
            end
        end
    end

    @testset "Acceptor recording" begin
        for M in 1:5
            dict = HBondDict(M)

            for k in 1:10
                # no chance this is 0
                recordacceptor!(dict, k, -(rand() + 0.1) * 10)
            end

            # all energies are negative
            @test let accs = acceptors(dict)
                all(accs) do (_, energy)
                    energy < 0
                end
            end

            # energies should be in increasing order (decreasing bond strength)
            @test let accs = acceptors(dict)
                Q = true; prev = -20.0
                for (_, energy) in accs
                    Q &= energy > prev
                    Q || break
                    prev = energy
                end
                Q
            end
        end
    end
end


text:
function DataFrames.rename!(df, ::DataFile{Val{:CV}})
    namemap = Dict("WE(1).Current (A)"=>"Current (mA)",
                   "WE(1).Potential (V)"=>"Potential (V)")
    rename!(df, namemap)
end

function process_data(::Val{:CV}, data; select)
    col, op, val = select
    for f in data["CV"]
        df = read_file(f)
        c = getproperty(df, col)
        fd = getindex(df, op.(c, val), :)
        push!(fd, fd[1, :])

        write_file(f, fd, ';')
    end
end