text:
abstract type AbstractAutoregressiveMatrix{T,V,S} <: AbstractMatrix{T} end
abstract type AbstractAutoregressiveMatrixAdjoint{T,V,S} <: AbstractMatrix{T} end

abstract type AbstractInvervalSpacing end
abstract type AbstractEvenSpacing <: AbstractInvervalSpacing end
abstract type AbstractUnevenSpacing <: AbstractInvervalSpacing end

struct EvenSpacing{T1,T2} <: AbstractEvenSpacing
    ρᵗ::T1
    invOmρ²ᵗ::T2
    rinvOmρ²ᵗ::T2
    # nρᵗrinvOmρ²ᵗ::T2
end
# We don't precompute, because storage cost is high..
struct UnevenSpacing <: AbstractUnevenSpacing end
struct CachedUnevenSpacing{T} <: AbstractUnevenSpacing
    ρᵗ::T
    invOmρ²ᵗ::T
    rinvOmρ²ᵗ::T
    # nρᵗrinvOmρ²ᵗ::T
end

"""
ρ is the autoregressive parameter.
τ is either:
1. A range of times, eg 1:10, -5:16, or 0.12:0.23:15. Use this for regularly spaced time intervals.
2. A vector of time deltas, eg [0.2, 0.15, 0.35, 0.3] would correspond to times [0.0, 0.2, 0.35, 0.7, 1.0].
    Use this for irregularly spaced time intervals. Note that the size of the matrix is 1 greater than the
    length of these deltas (eg, length(time deltas) == 4 in this example, corresponding to 5 times, or a 5x5 matrix).

"""
struct AutoregressiveMatrixLowerCholeskyInverse{T,V <: AbstractVector, S <: AbstractInvervalSpacing} <: AbstractAutoregressiveMatrix{T,V,S}
    ρ::T
    τ::V
    spacing::S
end
# struct AutoregressiveMatrixLowerCholeskyInverseAdjoint{T,V <: AbstractVector} <: AbstractAutoregressiveMatrixAdjoint{T,V}
#     ρ::T
#     ρᵗ::TV
#     invOmρ²ᵗ::TV
#     rinvOmρ²ᵗ::TV
#     τ::V
# end
struct AutoregressiveMatrix{T,V <: AbstractVector, S <: AbstractInvervalSpacing} <: AbstractAutoregressiveMatrix{T,V,S}
    ρ::T
    τ::V
    spacing::S
end
# struct AutoregressiveMatrixAdjoint{T,V <: AbstractVector} <: AbstractAutoregressiveMatrixAdjoint{T,V}
#     ρ::T
#     ρᵗ::TV
#     invOmρ²ᵗ::TV
#     rinvOmρ²ᵗ::TV
#     τ::V
# end

function AutoregressiveMatrixLowerCholeskyInverse(ρ::T, τ::AbstractUnitRange) where {T}
    invOmρ²ᵗ = 1 / ( 1 - ρ*ρ )
    rinvOmρ²ᵗ = sqrt(invOmρ²ᵗ)
    AutoregressiveMatrixLowerCholeskyInverse(
        ρ, τ, EvenSpacing(nothing, invOmρ²ᵗ, rinvOmρ²ᵗ)#, - ρ * rinvOmρ²ᵗ)
    )
end
function AutoregressiveMatrixLowerCholeskyInverse(ρ::T, τ::AbstractRange) where {T}
    ρᵗ = copysign(abs(ρ)^(step(τ)), ρ)
    invOmρ²ᵗ = 1 / ( 1 - ρᵗ*ρᵗ )
    rinvOmρ²ᵗ = sqrt(invOmρ²ᵗ)
    AutoregressiveMatrixLowerCholeskyInverse(
        ρ, τ, EvenSpacing(ρᵗ, invOmρ²ᵗ, rinvOmρ²ᵗ)#, - ρᵗ * rinvOmρ²ᵗ)
    )
end
# @generated function AutoregressiveMatrixLowerCholeskyInverse(ρ::T, τ::AbstractFixedSizeVector{M,T,L}) where {M,L,T}
#     quote
#         ρᵗ = FixedSizeVector{L,T}(undef)
#         invOmρ²ᵗ = = FixedSizeVector{L,T}(undef)
#         rinvOmρ²ᵗ = = FixedSizeVector{L,T}(undef)
#         @vectorize $T for i ∈ 1:$L
#             ρᵗ[i] = SIMDPirates.vcopysign(SIMDPirates.vpow(SIMDPirates.vabs(ρ),  τ[i] ), ρ[i])
#             vinvOmρ²ᵗ = 1 / (1 - ρᵗ[i]*ρᵗ[i])
#             invOmρ²ᵗ[i] = vinvOmρ²ᵗ
#             rinvOmρ²ᵗ[i] = sqrt(vinvOmρ²ᵗ)
#         end
#         AutoregressiveMatrixLowerCholeskyInverse(
#             ρ, ConstantFixedSizeVector(ρᵗ), ConstantFixedSizeVector(invOmρ²ᵗ), ConstantFixedSizeVector(rinvOmρ²ᵗ), τ
#         )
#     end
# end
"""
τ must be a vector of differences.
"""
function AutoregressiveMatrixLowerCholeskyInverse(ρ::T, τ::AbstractVector) where {T}
    cache(AutoregressiveMatrixLowerCholeskyInverse(
        ρ, τ, UnevenSpacing()
    ))
end
@inline function AutoregressiveMatrix(ρ::T, τ::AbstractFixedSizeVector{nT,T}) where {nT,T}
    AutoregressiveMatrix(ρ, τ, UnevenSpacing())
end
function AutoregressiveMatrix(ρ::T, τ::AbstractVector) where {T}
    ar = cache(AutoregressiveMatrix(
        ρ, τ, UnevenSpacing()
    ))
    AutoregressiveMatrix(ar.ρ, ar.τ, ar.spacing)
end
cache(A::AbstractAutoregressiveMatrix) = A
@generated function cache(A::AbstractAutoregressiveMatrix{T,V,UnevenSpacing}) where {M,T,L,V <: AbstractFixedSizeVector{M,T,L}}
    quote
        ρ = A.ρ
        τ = A.τ
        ρᵗ = FixedSizeVector{$M,$T}(undef)
        invOmρ²ᵗ = FixedSizeVector{$M,$T}(undef)
        rinvOmρ²ᵗ = FixedSizeVector{$M,$T}(undef)
        # nρᵗrinvOmρ²ᵗ = FixedSizeVector{$M,$T}(undef)
        if ρ != 0
            absρ = abs(ρ)
            @vectorize $T for i ∈ 1:$L
                ρᵗ[i] = SIMDPirates.vcopysign(SIMDPirates.vpow(absρ,  τ[i]), ρ)
                vinvOmρ²ᵗ = 1 / (1 - ρᵗ[i]*ρᵗ[i])
                invOmρ²ᵗ[i] = vinvOmρ²ᵗ
                rinvOmρ²ᵗ[i] = sqrt(vinvOmρ²ᵗ)
                # nρᵗrinvOmρ²ᵗ[i] = -ρᵗ[i] * rinvOmρ²ᵗ[i]
            end
        # elseif ρ < 0
        #     @vectorize $T for i ∈ 1:$M
        #         ρᵗ[i] = SIMDPirates.vcopysign(SIMDPirates.vpow(SIMDPirates.vabs(ρ),  τ[i] ), ρ)
        #         vinvOmρ²ᵗ = 1 / (1 - ρᵗ[i]*ρᵗ[i])
        #         invOmρ²ᵗ[i] = vinvOmρ²ᵗ
        #         rinvOmρ²ᵗ[i] = sqrt(vinvOmρ²ᵗ)
        #         # nρᵗrinvOmρ²ᵗ[i] = -ρᵗ[i] * rinvOmρ²ᵗ[i]
        #     end
        else # ρ == 0
            fill!(ρᵗ, zero(T))
            fill!(invOmρ²ᵗ, one(T))
            fill!(rinvOmρ²ᵗ, one(T))
        end
        AutoregressiveMatrixLowerCholeskyInverse(
            ρ, τ, CachedUnevenSpacing(
                ConstantFixedSizeVector(ρᵗ),
                ConstantFixedSizeVector(invOmρ²ᵗ),
                ConstantFixedSizeVector(rinvOmρ²ᵗ)#,
                # ConstantFixedSizeVector(nρᵗrinvOmρ²ᵗ)
            )
        )
    end
end

# function AutoregressiveMatrix(ρ::AutoregressiveParameter{T}, τ::AbstractRange) where {T}
#     absρ = abs(ρ)

# end
# function AutoregressiveMatrix(ρ::AutoregressiveParameter{T}, τ::AbstractUnitRange) where {T}
#     absρ = abs(ρ)
#     ptm1 = copysign(absρ, ρ)

# end
# function AutoregressiveMatrix(ρ::AutoregressiveParameter{T}, τ::AbstractVector) where {T}
#     AutoregressiveMatrix(ρ, τ)
# end

@inline Base.size(AR::AbstractAutoregressiveMatrix{T,V}) where {T,V<:AbstractRange} = (t = length(AR.τ); (t,t))
@inline Base.size(AR::AbstractAutoregressiveMatrix) = (t = length(AR.τ); (t,t))
# @inline Base.size(AR::AbstractAutoregressiveMatrixAdjoint{T,V}) where {T,V<:AbstractRange} = (t = length(AR.τ); (t,t))
# @inline Base.size(AR::AbstractAutoregressiveMatrixAdjoint) = (t = length(AR.τ)+1; (t,t))

function Base.getindex(AR::AutoregressiveMatrixLowerCholeskyInverse{T,V,S}, i, j) where {T,V<:AbstractUnitRange,S <: AbstractEvenSpacing}
    @boundscheck max(i,j) > length(AR.τ) || PaddedMatrices.ThrowBoundsError("max(i,j) = $(max(i,j)) > length(AR.τ) = $(length(AR.τ)).")
    if i == j
        if i == 1
            return one(T)
        else
            return AR.spacing.rinvOmρ²ᵗ
        end
    end
    # i - j  == 1 ? AR.spacing.nρᵗrinvOmρ²ᵗ : zero(T)
    i - j  == 1 ? - AR.ρ * AR.spacing.rinvOmρ²ᵗ : zero(T)
end

function Base.getindex(AR::AutoregressiveMatrixLowerCholeskyInverse{T,V,S}, i, j) where {T,V<:AbstractRange,S <: AbstractEvenSpacing}
    @boundscheck max(i,j) > length(AR.τ) || PaddedMatrices.ThrowBoundsError("max(i,j) = $(max(i,j)) > length(AR.τ) = $(length(AR.τ)).")
    if i == j
        if i == 1
            return one(T)
        else
            return AR.spacing.rinvOmρ²ᵗ
        end
    end
    # i - j  == 1 ? AR.spacing.nρᵗrinvOmρ²ᵗ : zero(T)
    i - j  == 1 ? - AR.spacing.ρᵗ * AR.spacing.rinvOmρ²ᵗ : zero(T)
end

function Base.getindex(AR::AutoregressiveMatrixLowerCholeskyInverse{T,V,S}, i, j) where {T,V,S <: UnevenSpacing}
    @boundscheck max(i,j)-1 > length(AR.τ) && PaddedMatrices.ThrowBoundsError("max(i,j) = $(max(i,j)) > length(AR.τ) = $(length(AR.τ)).")
    if i == j
        if i == 1
            return one(T)
        else
            rhot = AR.ρ ^ (AR.τ[i-1])
            return 1 / sqrt(1 - rhot * rhot)
        end
    end
    δij = i - j
    δij == 1 || return zero(T)
    rhot = AR.ρ ^ (AR.τ[j])
    - rhot / sqrt(1 - rhot*rhot)
end
function Base.getindex(AR::AutoregressiveMatrixLowerCholeskyInverse{T,V,S}, i, j) where {T,V,S <: CachedUnevenSpacing}
    @boundscheck max(i,j)-1 > length(AR.τ) && PaddedMatrices.ThrowBoundsError("max(i,j) = $(max(i,j)) > length(AR.τ) = $(length(AR.τ)).")
    if i == j
        if i == 1
            return one(T)
        else
            return AR.spacing.rinvOmρ²ᵗ[i-1]
        end
    end
    δij = i - j
    δij == 1 ? -AR.spacing.ρᵗ[j] * AR.spacing.rinvOmρ²ᵗ[j] : zero(T)
end

@inline function Base.getindex(AR::AutoregressiveMatrix{T,V,S}, i, j) where {T,V<:AbstractUnitRange,S <: AbstractEvenSpacing}
    @boundscheck max(i,j)-1 > length(AR.τ) && PaddedMatrices.ThrowBoundsError("max(i,j) = $(max(i,j)) > length(AR.τ) = $(length(AR.τ)).")
    SIMDPirates.vcopysign(SLEEFPirates.power(AR.ρ, SIMDPirates.abs(i - j), AR.ρ))
end

@inline function Base.getindex(AR::AutoregressiveMatrix{T,V,S}, i, j) where {T,V<:AbstractRange,S <: AbstractEvenSpacing}
    @boundscheck max(i,j)-1 > length(AR.τ) && PaddedMatrices.ThrowBoundsError("max(i,j) = $(max(i,j)) > length(AR.τ) = $(length(AR.τ)).")
    SLEEFPirates.power(AR.ρ, step(AR.τ) * SIMDPirates.abs(i - j))
end

function Base.getindex(AR::AutoregressiveMatrix{T,V,S}, i, j) where {T,V <: AbstractVector,S <: AbstractUnevenSpacing}
    @boundscheck max(i,j)-1 > length(AR.τ) && PaddedMatrices.ThrowBoundsError("max(i,j) = $(max(i,j)) > length(AR.τ) = $(length(AR.τ)).")
    i == j && return one(T)
    δt = zero(T)
    i, j = minmax(i, j)
    for k ∈ i:j-1
        δt += AR.τ[k]
    end
    ρᵗ = abs(AR.ρ)^δt
    iseven(j - i) ? ρᵗ : copysign(ρᵗ, AR.ρ)
end

@inline function PaddedMatrices.FixedSizeMatrix(
        AR::Union{AutoregressiveMatrix{T,<:AbstractFixedSizeVector{M}},AutoregressiveMatrixLowerCholeskyInverse{T,<:AbstractFixedSizeVector{M}}}) where {T,M}
    pAR = FixedSizeMatrix{M+1,M+1,T}(undef)
    @inbounds for mc ∈ 1:M+1
        for mr ∈ 1:M+1
            pAR[mr,mc] = AR[mr,mc]
        end
    end
    pAR
end
function PaddedMatrices.ConstantFixedSizeMatrix(
        AR::Union{AutoregressiveMatrix{T,<:AbstractFixedSizeVector{M}},AutoregressiveMatrixLowerCholeskyInverse{T,<:AbstractFixedSizeVector{M}}}) where {T,M}
    ConstantFixedSizeMatrix(FixedSizeMatrix(AR))
end

@generated function Base.:*(
            A::AbstractAutoregressiveMatrix{T,R,S},
            B::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}}
        ) where {M,N,W,T,R,S}

    # register_count = VectorizationBase.REGISTER_COUNT

    V = NTuple{W,Core.VecElement{T}}
    # if R <: AbstractUnitRange
    #     loop_body = :( $vinvOmρ² * ( B[ 1+m + M*n ] - A.ρ * B[ m + M*n ] ) )
    if S <: AbstractEvenSpacing
        loop_body = :( AB[ 1+m + $M*n ] = SIMDPirates.vmul( vinvOmρ², SIMDPirates.vfnmadd( vρ, B[ m + $M*n ], B[ 1+m + $M*n ] ) ) )
    elseif S == UnevenSpacing
        loop_body = quote
            ρᵗ = A.ρ ^ (A.τ[m])
            rinvOmρ²ᵗ = 1 / sqrt(1 - ρᵗ*ρᵗ)
            AB[ 1+m + $M*n ] = SIMDPirates.vmul( SIMDPirates.vbroadcast($V, rinvOmρ²ᵗ), SIMDPirates.vfnmadd( SIMDPirates.vbroadcast($V, ρᵗ), B[ m + M*n ], B[ 1+m + M*n ] ) )
        end
    elseif S == CachedUnevenSpacing
        loop_body = quote
            vinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ[m])
            vρᵗ = SIMDPirates.vbroadcast($V, A.ρᵗ[m])
            AB[ 1+m + $M*n ] = SIMDPirates.vmul( vinvOmρ², SIMDPirates.vfnmadd( vρᵗ, B[ m + M*n ], B[ 1+m + M*n ] ) )
        end
    end

    q = quote
        AB = FixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}}(undef)
        triind = N
        @inbounds for n ∈ 0:N-1
            AB[ 1 + M*n ] = B[ 1 + M*n ]
            for m ∈ 1:M-1
                $loop_body
            end
        end
        ConstantFixedSizeMatrix(AB)
    end
    if R <: AbstractUnitRange
        pushfirst!(q.args, quote
            vinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.ρ)
        end)
    elseif R <: AbstractRange
        pushfirst!(q.args, quote
            vinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ)
        end)
    end
    q
end

function LinearAlgebra.det(A::AbstractAutoregressiveMatrix{T,R,S}) where {T,R,S<:EvenSpacing}
    A.spacing.rinvOmρ²ᵗ ^ (length(A.τ) - 1)
end
function LinearAlgebra.logdet(A::AbstractAutoregressiveMatrix{T,R,S}) where {T,R,S<:EvenSpacing}
    (length(A.τ) - 1) * log(A.spacing.rinvOmρ²ᵗ)
end

function ∂logdet(A::AbstractAutoregressiveMatrix{T,R,S}) where {T,R <: AbstractUnitRange,S <: EvenSpacing}
    N = (length(A.τ) - 1)
    N * log(A.spacing.rinvOmρ²ᵗ), N * A.τ * A.ρ * A.spacing.invOmρ²ᵗ[i]
end
function ∂logdet(A::AbstractAutoregressiveMatrix{T,R,S}) where {T,R <: AbstractRange,S <: EvenSpacing}
    N = (length(A.τ) - 1)
    N * log(A.spacing.rinvOmρ²ᵗ), N * A.τ * A.spacing.ρᵗ * A.spacing.ρᵗ * A.spacing.invOmρ²ᵗ / A.ρ
end

@generated function LinearAlgebra.logdet(A::AbstractAutoregressiveMatrix{T,R,S}) where {T,R<:AbstractVector,S <: CachedUnevenSpacing}
    quote
        out = zero(T)
        rinvOmρ²ᵗ = A.spacing.rinvOmρ²ᵗ
        @vectorize $T for i ∈ eachindex(A.τ)
            out += log(rinvOmρ²ᵗ[i])
        end
        out
    end
end
@generated function ∂logdet(A::AbstractAutoregressiveMatrix{T,R,S}) where {T,R<:AbstractVector,S <: CachedUnevenSpacing}
    quote
        out = zero(T)
        ∂out = zero(T)
        rinvOmρ²ᵗ = A.spacing.rinvOmρ²ᵗ
        invOmρ²ᵗ = A.spacing.invOmρ²ᵗ
        ρᵗ = A.spacing.ρᵗ
        ρ = A.ρ
        τ = A.τ
        # if ρ != 0
        @vectorize $T for i ∈ eachindex(τ)
            out += log(rinvOmρ²ᵗ[i])
            ∂out += τ[i] * ρᵗ[i] * ρᵗ[i] * invOmρ²ᵗ[i] / ρ
        end
        return out, ∂out
        # else
        #     @vectorize $T for i ∈ eachindex(τ)
        #         ∂out += τ[i] * ρᵗ[i] * SIMDPirates.vabs(ρ[i-1] * invOmρ²ᵗ[i] / ρ
        #     end
        #     return out, ∂out
        # end
        return out, ∂out
    end
end
function ∂det(A::AbstractAutoregressiveMatrix)
    logdetA, ∂logdetA = ∂logdet(A)
    detA = exp(logdetA)
    ∂detA = ∂logdetA * detA
    detA, ∂detA
end

@generated function quadform(
            A::AbstractAutoregressiveMatrix{T,R,S},
            B::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}}
        # ) where {T,R,M,N,W,S}
        ) where {T,R,M,N,W,S}
    V = NTuple{W,Core.VecElement{T}}
    # broadcastB = VoT == T
    q = quote
        $(Expr(:meta, :inline))
        qf = SIMDPirates.vbroadcast($V, zero($T))
    end
    if R <: AbstractUnitRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.ρ)
            ρ⁻¹2 = SIMDPirates.vbroadcast($V, one($T) / A.ρ)
            ρ²ᵗ = SIMDPirates.vbroadcast($V, A.ρ*A.ρ)
        end)
    elseif R <: AbstractRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ)
            ρ⁻¹2 = SIMDPirates.vbroadcast($V, $T(step(A.τ)) / A.ρ)
            ρ²ᵗ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ*A.spacing.ρᵗ)
        end)
    # else
    end
    if S <: AbstractEvenSpacing
        push!(q.args, quote
            @inbounds for n ∈ 0:$(N-1)
                δBC_last = B[1 + $M*n]
                A_1 = δBC_last
                # product[1 + $M*n] = A_m
                δ_last = A_1

                # δ_previous = one($T)
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                for m ∈ 1:$(M-1)
                    δBC_next = B[m+1 + $M*n]
                    ρᵗB = SIMDPirates.vmul( vρ, δBC_last)
                    A_m = SIMDPirates.vmul( vrinvOmρ²ᵗ, SIMDPirates.vsub(δBC_next, ρᵗB ) )

                    δ_next = SIMDPirates.vmul( A_m, vrinvOmρ²ᵗ )
                    qf = SIMDPirates.vmuladd(A_m, A_m, qf)
                    # product[m + $M*n] = SIMDPirates.vfnmadd(vρ, δ_next, δ_last)

                    δBC_last = δBC_next
                    # δ_last = δ_next
                end
                # product[$M + $M*n] = δ_last
            end
        end)
    elseif S <: CachedUnevenSpacing
        push!(q.args, quote
            # $([ :($(Symbol(:riOmρ²_,m)) = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ[$m])) for m ∈ 2:M]...)
            # $([ :($(Symbol(:ρᵗ_,    m)) = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ[$m]))        for m ∈ 2:M]...)
            # $([ :($(Symbol(:ρ²ᵗ_,   m)) = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ[$m]*A.spacing.ρᵗ[$m])) for m ∈ 2:M]...)
            # $([ :($(Symbol(:ρ⁻¹2_, m)) = SIMDPirates.vbroadcast($V, AR.τ[$m] / A.ρ)) for m ∈ 2:M]...)
            @inbounds ρ²ᵗtup = $(Expr(:tuple, [ :(A.spacing.ρᵗ[$m]*A.spacing.ρᵗ[$m]) for m ∈ 1:M-1]...))
            @inbounds ρ⁻¹2tup = $(Expr(:tuple, [ :(A.τ[$m] / A.ρ) for m ∈ 1:M-1]...))

            @inbounds for n ∈ 0:$(N-1)
                δBC_last = B[1 + $M*n]
                A_1 = δBC_last
                # product[1 + $M*n] = A_m
                δ_last = A_1

                # δ_previous = one($T)
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                for m ∈ 1:$(M-1)
                    δBC_next = B[m+1 + $M*n]
                    ρᵗB = SIMDPirates.vmul( A.spacing.ρᵗ[m], δBC_last)
                    A_m = SIMDPirates.vmul( A.spacing.rinvOmρ²ᵗ[m], SIMDPirates.vsub(δBC_next, ρᵗB ) )

                    δ_next = SIMDPirates.vmul( A_m, A.spacing.rinvOmρ²ᵗ[m] )
                    qf = SIMDPirates.vmuladd(A_m, A_m, qf)
                    # product[m + $M*n] = SIMDPirates.vfnmadd(A.spacing.ρᵗ[m], δ_next, δ_last)

                    δBC_last = δBC_next
                    # δ_last = δ_next
                end
                # product[$M + $M*n] = δ_last
            end
        end)
    else
        throw("Spacing type $S not yet supported. Perhaps try caching intermediate results: `cache(AR_matrix)`")
    end
    push!(q.args, :qf)

    q

end

@generated function quadformdiff(
            A::AbstractAutoregressiveMatrix{T,R},
            B::AbstractFixedSizeMatrix{M,N,VoT},
            C::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}}
        ) where {T,R,M,N,W,VoT}
    V = NTuple{W,Core.VecElement{T}}
    q = quote
        $(Expr(:meta, :inline))
        qf = SIMDPirates.vbroadcast($V, zero($T))
    end
    if R <: AbstractUnitRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.ρ)
        end)
    elseif R <: AbstractRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ)
        end)
    end
    if R <: AbstractRange
        push!(q.args, quote
            @inbounds for n ∈ 0:$(N-1)
                δBC_last = SIMDPirates.vsub(B[1 + $M*n], C[1 + $M*n])
                A_1 = δBC_last
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                for m ∈ 2:$M
                    δBC_next = SIMDPirates.vsub(B[m + $M*n], C[m + $M*n])
                    A_m = SIMDPirates.vmul( vrinvOmρ², SIMDPirates.vfnmadd( vρ, δBC_last, δBC_next ) )
                    qf = SIMDPirates.vmuladd(A_m, A_m, qf)
                    δBC_last = δBC_next
                end
            end
        end)
    else
    # loop_body = :( SIMDPirates.vmul( , SIMDPirates.vfnmadd( , B[ m + M*n ], B[ 1+m + M*n ] ) ) )
        push!(q.args, quote
            # @inbounds ρᵗtup = $(Expr(:tuple, [:(A.spacing.ρᵗ[$m]) for m ∈ 1:M-1]...))
            # @inbounds riOmρ²tup = $(Expr(:tuple, [:(A.spacing.rinvOmρ²ᵗ[$m]) for m ∈ 1:M-1]...))
            # $([ :($(Symbol(:ρᵗ_,    m)) = SIMDPirates.vbroadcast($V, A.ρᵗ[$m]))        for m ∈ 2:M]...)
            # $([ :($(Symbol(:riOmρ²_,m)) = SIMDPirates.vbroadcast($V, A.rinvOmρ²ᵗ[$m])) for m ∈ 2:M]...)
            @inbounds for n ∈ 0:$(N-1)
                δBC_last = SIMDPirates.vsub(B[1 + $M*n], C[1 + $M*n])
                A_1 = δBC_last
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                for m ∈ 1:$(M-1)
                    δBC_next = SIMDPirates.vsub(B[m+1 + $M*n], C[m+1 + $M*n])
                    A_m = SIMDPirates.vmul( A.spacing.rinvOmρ²ᵗ[m], SIMDPirates.vfnmadd( A.spacing.ρᵗ[m], δBC_last, δBC_next ) )
                    qf = SIMDPirates.vmuladd(A_m, A_m, qf)
                    δBC_last = δBC_next
                end
            end
        end)
    end

    push!(q.args, :qf)
    q

end

@generated function mul_and_quadform(
            A::AbstractAutoregressiveMatrix{T,R},
            B::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}}
        ) where {T,R,M,N,W}
    V = NTuple{W,Core.VecElement{T}}
    q = quote
        $(Expr(:meta, :inline))
        qf = SIMDPirates.vbroadcast($V, zero($T))
        product = FixedSizeMatrix{$M,$N,$V}(undef)
    end
    if R <: AbstractUnitRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.ρ)
        end)
    elseif R <: AbstractRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ)
        end)
    # else
    end
    if R <: AbstractRange
        push!(q.args, quote
            for n ∈ 0:$(N-1)
                A_1 = B[1 + $M*n]
                product[1 + $M*n] = A_1
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                $([quote
                        $(Symbol(:A_,m)) = SIMDPirates.vmul( vrinvOmρ², SIMDPirates.vfnmadd( vρ, B[ $(m-1) + $M*n ], B[ $m + $M*n ] ) )
                        product[$m + $M*n] = $(Symbol(:A_,m))
                        qf = SIMDPirates.vmuladd($(Symbol(:A_,m)), $(Symbol(:A_,m)), qf)
                    end for m ∈ 2:M]...)
            end
        end)
    else
        # loop_body = :( SIMDPirates.vmul( , SIMDPirates.vfnmadd( , B[ m + M*n ], B[ 1+m + M*n ] ) ) )
        push!(q.args, quote
            $([ :($(Symbol(:ρᵗ_,    m)) = SIMDPirates.vbroadcast($V, A.ρᵗ[$m]))        for m ∈ 2:M]...)
            $([ :($(Symbol(:riOmρ²_,m)) = SIMDPirates.vbroadcast($V, A.rinvOmρ²ᵗ[$m])) for m ∈ 2:M]...)
            for n ∈ 0:$(N-1)
                A_1 = B[1 + $M*n]
                product[1 + $M*n] = A_1
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                $([quote
                        $(Symbol(:A_,m)) = SIMDPirates.vmul( $(Symbol(:riOmρ²_,m)), SIMDPirates.vfnmadd( $(Symbol(:ρᵗ,m)), B[ $(m-1) + $M*n ], B[ $m + $M*n ] ) )
                        product[$m + $M*n] = $(Symbol(:A_,m))
                        qf = SIMDPirates.vmuladd($(Symbol(:A_,m)), $(Symbol(:A_,m)), qf)
                    end for m ∈ 2:M]...)
            end
        end)
    end

    push!(q.args, :(qf, ConstantFixedSizeMatrix(product)))

    q

end
function ∂mul_and_quadform(
            A::AbstractAutoregressiveMatrix{T,R,S},
            B::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}}
        ) where {T,R,M,N,W,S}
    V = NTuple{W,Core.VecElement{T}}
    q = quote
        $(Expr(:meta, :inline))
        qf = SIMDPirates.vbroadcast($V, zero($T))
        product = FixedSizeMatrix{$M,$N,$V}(undef)
    end
    if R <: AbstractUnitRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.ρ)
            ρ⁻¹2 = SIMDPirates.vbroadcast($V, 2 / A.ρ)
        end)
    elseif R <: AbstractRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ)
            ρ⁻¹2 = SIMDPirates.vbroadcast($V, T(2*step(AR.τ)) / A.ρ)
        end)
    # else
    end
    if S <: AbstractEvenSpacing
        push!(q.args, quote
            ρ²ᵗ = SIMDPirates.vbroadcast($V, $( R <: AbstractUnitRange ? :(A.ρ*A.ρ) : :(A.spacing.ρᵗ*A.spacing.ρᵗ)  ) )
            for n ∈ 0:$(N-1)
                A_1 = B[1 + $M*n]
                product[1 + $M*n] = A_1
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                $([quote
                        ρᵗB = SIMDPirates.vmul(vρ, B[$(m-1)+$M*n ])
                        $(Symbol(:A_,m)) = SIMDPirates.vmul( vrinvOmρ², SIMDPirates.vsub(B[$m+$M*n], ρᵗB ) )
                        product[$m + $M*n] = $(Symbol(:A_,m))

                        qf = SIMDPirates.vmuladd($(Symbol(:A_,m)), $(Symbol(:A_,m)), qf)

                        δ = SIMDPirates.vmul( $(Symbol(:A_,m)), vrinvOmρ² )
                        ∂out = SIMDPirates.vmuladd(SIMDPirates.vmul(δ, ρ⁻¹2), SIMDPirates.vfmsub(δ, ρ²ᵗ, ρᵗB), ∂out)

                    end for m ∈ 2:M]...)
            end
        end)
    elseif S == CachedUnevenSpacing
        push!(q.args, quote
            $([ :($(Symbol(:riOmρ²_,m)) = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ[$m])) for m ∈ 2:M]...)
            $([ :($(Symbol(:ρᵗ_,    m)) = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ[$m]))        for m ∈ 2:M]...)
            $([ :($(Symbol(:ρ²ᵗ_,   m)) = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ[$m]*A.spacing.ρᵗ[$m])) for m ∈ 2:M]...)
            $([ :($(Symbol(:ρ⁻¹2_, m)) = SIMDPirates.vbroadcast($V, T(2*AR.τ[$m]) / A.ρ)) for m ∈ 2:M]...)
            for n ∈ 0:$(N-1)
                A_1 = B[1 + $M*n]
                product[1 + $M*n] = A_1
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                $([quote
                        ρᵗB = SIMDPirates.vmul($(Symbol(:ρᵗ,     m)), B[$(m-1)+$M*n ])
                        $(Symbol(:A_,m)) = SIMDPirates.vmul( $(Symbol(:riOmρ²_,m)), SIMDPirates.vsub(B[$m+$M*n], ρᵗB ) )
                        product[$m + $M*n] = $(Symbol(:A_,m))

                        qf = SIMDPirates.vmuladd($(Symbol(:A_,m)), $(Symbol(:A_,m)), qf)

                        δ = SIMDPirates.vmul( $(Symbol(:A_,m)), $(Symbol(:riOmρ²_,m)) )
                        ∂out = SIMDPirates.vmuladd(SIMDPirates.vmul(δ, $(Symbol(:ρ⁻¹2_, m))), SIMDPirates.vfmsub(δ, $(Symbol(:ρ²ᵗ_,   m)), ρᵗB), ∂out)

                    end for m ∈ 2:M]...)
            end
        end)
    else
        throw("Spacing type $S not yet supported. Perhaps try caching intermediate results: `cache(AR_matrix)`")
    end

    # push!(q.args, :(qf, ConstantFixedSizeMatrix(product)))
    # relying on inline to avoid allocations
    push!(q.args, :(qf, product))

    q

end

function ∂quadform_quote(M,N,W,T,R,S)
    V = NTuple{W,Core.VecElement{T}}
    # broadcastB = VoT == T
    q = quote
        # $(Expr(:meta, :inline))
        qf = SIMDPirates.vbroadcast($V, zero($T))
        # product = FixedSizeMatrix{$M,$N,$V}(undef)
        ∂out = SIMDPirates.vbroadcast($V, zero($T))
    end
    if R <: AbstractUnitRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.ρ)
            ρ⁻¹2 = SIMDPirates.vbroadcast($V, one($T) / A.ρ)
            ρ²ᵗ = SIMDPirates.vbroadcast($V, A.ρ*A.ρ)
        end)
    elseif R <: AbstractRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ)
            ρ⁻¹2 = SIMDPirates.vbroadcast($V, $T(step(A.τ)) / A.ρ)
            ρ²ᵗ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ*A.spacing.ρᵗ)
        end)
    # else
    end
    if S <: AbstractEvenSpacing
        push!(q.args, quote
            @inbounds for n ∈ 0:$(N-1)
                δBC_last = B[1 + $M*n]
                A_1 = δBC_last
                # product[1 + $M*n] = A_m
                δ_last = A_1

                # δ_previous = one($T)
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                for m ∈ 1:$(M-1)
                    δBC_next = B[m+1 + $M*n]
                    ρᵗB = SIMDPirates.vmul( vρ, δBC_last)
                    A_m = SIMDPirates.vmul( vrinvOmρ²ᵗ, SIMDPirates.vsub(δBC_next, ρᵗB ) )

                    δ_next = SIMDPirates.vmul( A_m, vrinvOmρ²ᵗ )
                    qf = SIMDPirates.vmuladd(A_m, A_m, qf)
                    # product[m + $M*n] = SIMDPirates.vfnmadd(vρ, δ_next, δ_last)

                    ∂out = SIMDPirates.vmuladd(SIMDPirates.vmul(δ_next, ρ⁻¹2), SIMDPirates.vfnmadd(δ_next, ρ²ᵗ, ρᵗB), ∂out)
                    δBC_last = δBC_next
                    # δ_last = δ_next
                end
                # product[$M + $M*n] = δ_last
            end
        end)
    elseif S <: CachedUnevenSpacing
        push!(q.args, quote
            # $([ :($(Symbol(:riOmρ²_,m)) = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ[$m])) for m ∈ 2:M]...)
            # $([ :($(Symbol(:ρᵗ_,    m)) = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ[$m]))        for m ∈ 2:M]...)
            # $([ :($(Symbol(:ρ²ᵗ_,   m)) = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ[$m]*A.spacing.ρᵗ[$m])) for m ∈ 2:M]...)
            # $([ :($(Symbol(:ρ⁻¹2_, m)) = SIMDPirates.vbroadcast($V, AR.τ[$m] / A.ρ)) for m ∈ 2:M]...)
            @inbounds ρ²ᵗtup = $(Expr(:tuple, [ :(A.spacing.ρᵗ[$m]*A.spacing.ρᵗ[$m]) for m ∈ 1:M-1]...))
            @inbounds ρ⁻¹2tup = $(Expr(:tuple, [ :(A.τ[$m] / A.ρ) for m ∈ 1:M-1]...))

            @inbounds for n ∈ 0:$(N-1)
                δBC_last = B[1 + $M*n]
                A_1 = δBC_last
                # product[1 + $M*n] = A_m
                δ_last = A_1

                # δ_previous = one($T)
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                for m ∈ 1:$(M-1)
                    δBC_next = B[m+1 + $M*n]
                    ρᵗB = SIMDPirates.vmul( A.spacing.ρᵗ[m], δBC_last)
                    A_m = SIMDPirates.vmul( A.spacing.rinvOmρ²ᵗ[m], SIMDPirates.vsub(δBC_next, ρᵗB ) )

                    δ_next = SIMDPirates.vmul( A_m, A.spacing.rinvOmρ²ᵗ[m] )
                    qf = SIMDPirates.vmuladd(A_m, A_m, qf)
                    # product[m + $M*n] = SIMDPirates.vfnmadd(A.spacing.ρᵗ[m], δ_next, δ_last)

                    ∂out = SIMDPirates.vmuladd(SIMDPirates.vmul(δ_next, ρ⁻¹2tup[m]), SIMDPirates.vfnmadd(δ_next, ρ²ᵗtup[m], ρᵗB), ∂out)
                    δBC_last = δBC_next
                    # δ_last = δ_next
                end
                # product[$M + $M*n] = δ_last
            end
        end)
    else
        throw("Spacing type $S not yet supported. Perhaps try caching intermediate results: `cache(AR_matrix)`")
    end

    # push!(q.args, :(qf, -0.5*(SIMDPirates.vsum(∂out)), ConstantFixedSizeMatrix(product)))
    # relying on inlining to avoid allocations.
    # push!(q.args, :(qf, ∂out, product))

    q
end

@generated function ∂quadform(
            A::AbstractAutoregressiveMatrix{T,R,S},
            B::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}}
        # ) where {T,R,M,N,W,S}
        ) where {T,R,M,N,W,S}

    quote
        $(Expr(:meta, :inline))
        $(∂quadform_quote(M,N,W,T,R,S))
        qf, ∂out#, product
    end
end

function selfcrossmul_and_quadform_quote(M,N,W,T,R,S)
    V = NTuple{W,Core.VecElement{T}}
    # broadcastB = VoT == T
    q = quote
        # $(Expr(:meta, :inline))
        qf = SIMDPirates.vbroadcast($V, zero($T))
        # product = FixedSizeMatrix{$M,$N,$V}(undef)
        # ∂out = SIMDPirates.vbroadcast($V, zero($T))
    end
    if R <: AbstractUnitRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.ρ)
            ρ⁻¹2 = SIMDPirates.vbroadcast($V, one($T) / A.ρ)
            ρ²ᵗ = SIMDPirates.vbroadcast($V, A.ρ*A.ρ)
        end)
    elseif R <: AbstractRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ)
            ρ⁻¹2 = SIMDPirates.vbroadcast($V, $T(step(A.τ)) / A.ρ)
            ρ²ᵗ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ*A.spacing.ρᵗ)
        end)
    # else
    end
    if S <: AbstractEvenSpacing
        push!(q.args, quote
            @inbounds for n ∈ 0:$(N-1)
                δBC_last = B[1 + $M*n]
                A_1 = δBC_last
                # product[1 + $M*n] = A_m
                δ_last = A_1

                # δ_previous = one($T)
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                for m ∈ 1:$(M-1)
                    δBC_next = B[m+1 + $M*n]
                    ρᵗB = SIMDPirates.vmul( vρ, δBC_last)
                    A_m = SIMDPirates.vmul( vrinvOmρ²ᵗ, SIMDPirates.vsub(δBC_next, ρᵗB ) )

                    δ_next = SIMDPirates.vmul( A_m, vrinvOmρ²ᵗ )
                    qf = SIMDPirates.vmuladd(A_m, A_m, qf)
                    product[m + $M*n] = SIMDPirates.vfnmadd(vρ, δ_next, δ_last)

                    # ∂out = SIMDPirates.vmuladd(SIMDPirates.vmul(δ_next, ρ⁻¹2), SIMDPirates.vfnmadd(δ_next, ρ²ᵗ, ρᵗB), ∂out)
                    δBC_last = δBC_next
                    δ_last = δ_next
                end
                product[$M + $M*n] = δ_last
            end
        end)
    elseif S <: CachedUnevenSpacing
        push!(q.args, quote
            # $([ :($(Symbol(:riOmρ²_,m)) = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ[$m])) for m ∈ 2:M]...)
            # $([ :($(Symbol(:ρᵗ_,    m)) = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ[$m]))        for m ∈ 2:M]...)
            # $([ :($(Symbol(:ρ²ᵗ_,   m)) = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ[$m]*A.spacing.ρᵗ[$m])) for m ∈ 2:M]...)
            # $([ :($(Symbol(:ρ⁻¹2_, m)) = SIMDPirates.vbroadcast($V, AR.τ[$m] / A.ρ)) for m ∈ 2:M]...)
            @inbounds ρ²ᵗtup = $(Expr(:tuple, [ :(A.spacing.ρᵗ[$m]*A.spacing.ρᵗ[$m]) for m ∈ 1:M-1]...))
            @inbounds ρ⁻¹2tup = $(Expr(:tuple, [ :(A.τ[$m] / A.ρ) for m ∈ 1:M-1]...))

            @inbounds for n ∈ 0:$(N-1)
                δBC_last = B[1 + $M*n]
                A_1 = δBC_last
                # product[1 + $M*n] = A_m
                δ_last = A_1

                # δ_previous = one($T)
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                for m ∈ 1:$(M-1)
                    δBC_next = B[m+1 + $M*n]
                    ρᵗB = SIMDPirates.vmul( A.spacing.ρᵗ[m], δBC_last)
                    A_m = SIMDPirates.vmul( A.spacing.rinvOmρ²ᵗ[m], SIMDPirates.vsub(δBC_next, ρᵗB ) )

                    δ_next = SIMDPirates.vmul( A_m, A.spacing.rinvOmρ²ᵗ[m] )
                    qf = SIMDPirates.vmuladd(A_m, A_m, qf)
                    product[m + $M*n] = SIMDPirates.vfnmadd(A.spacing.ρᵗ[m], δ_next, δ_last)

                    # ∂out = SIMDPirates.vmuladd(SIMDPirates.vmul(δ_next, ρ⁻¹2tup[m]), SIMDPirates.vfnmadd(δ_next, ρ²ᵗtup[m], ρᵗB), ∂out)
                    δBC_last = δBC_next
                    δ_last = δ_next
                end
                product[$M + $M*n] = δ_last
            end
        end)
    else
        throw("Spacing type $S not yet supported. Perhaps try caching intermediate results: `cache(AR_matrix)`")
    end

    # push!(q.args, :(qf, -0.5*(SIMDPirates.vsum(∂out)), ConstantFixedSizeMatrix(product)))
    # relying on inlining to avoid allocations.
    # push!(q.args, :(qf, ∂out, product))

    q
end

"""
D = A(ρ) * (B - C)
qf = self_dot(D)

returns: qf, -0.5* ∂qf/∂ρ, A(ρ)' * D
"""
@generated function selfcrossmul_and_quadform(
            A::AbstractAutoregressiveMatrix{T,R,S},
            B::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}}
        # ) where {T,R,M,N,W,S}
        ) where {T,R,M,N,W,S}

    quote
        $(Expr(:meta, :inline))
        product = FixedSizeMatrix{$M,$N,NTuple{$W,Core.VecElement{$T}}}(undef)
        $(selfcrossmul_and_quadform_quote(M,N,W,T,R,S))
        qf, product
    end
end
@generated function selfcrossmul_and_quadform!(
            product::PaddedMatrices.AbstractMutableFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            A::AbstractAutoregressiveMatrix{T,R,S},
            B::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}}
        # ) where {T,R,M,N,W,S}
        ) where {T,R,M,N,W,S}

    quote
        $(Expr(:meta, :inline))
        $(selfcrossmul_and_quadform_quote(M,N,W,T,R,S))
        qf#, product
    end
end

function ∂selfcrossmul_and_quadform_quote(M,N,W,T,R,S)
    V = NTuple{W,Core.VecElement{T}}
    # broadcastB = VoT == T
    q = quote
        # $(Expr(:meta, :inline))
        qf = SIMDPirates.vbroadcast($V, zero($T))
        # product = FixedSizeMatrix{$M,$N,$V}(undef)
        ∂out = SIMDPirates.vbroadcast($V, zero($T))
    end
    if R <: AbstractUnitRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.ρ)
            ρ⁻¹2 = SIMDPirates.vbroadcast($V, one($T) / A.ρ)
            ρ²ᵗ = SIMDPirates.vbroadcast($V, A.ρ*A.ρ)
        end)
    elseif R <: AbstractRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ)
            ρ⁻¹2 = SIMDPirates.vbroadcast($V, $T(step(A.τ)) / A.ρ)
            ρ²ᵗ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ*A.spacing.ρᵗ)
        end)
    # else
    end
    if S <: AbstractEvenSpacing
        push!(q.args, quote
            @inbounds for n ∈ 0:$(N-1)
                δBC_last = B[1 + $M*n]
                A_1 = δBC_last
                # product[1 + $M*n] = A_m
                δ_last = A_1

                # δ_previous = one($T)
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                for m ∈ 1:$(M-1)
                    δBC_next = B[m+1 + $M*n]
                    ρᵗB = SIMDPirates.vmul( vρ, δBC_last)
                    A_m = SIMDPirates.vmul( vrinvOmρ²ᵗ, SIMDPirates.vsub(δBC_next, ρᵗB ) )

                    δ_next = SIMDPirates.vmul( A_m, vrinvOmρ²ᵗ )
                    qf = SIMDPirates.vmuladd(A_m, A_m, qf)
                    product[m + $M*n] = SIMDPirates.vfnmadd(vρ, δ_next, δ_last)

                    ∂out = SIMDPirates.vmuladd(SIMDPirates.vmul(δ_next, ρ⁻¹2), SIMDPirates.vfnmadd(δ_next, ρ²ᵗ, ρᵗB), ∂out)
                    δBC_last = δBC_next
                    δ_last = δ_next
                end
                product[$M + $M*n] = δ_last
            end
        end)
    elseif S <: CachedUnevenSpacing
        push!(q.args, quote
            # $([ :($(Symbol(:riOmρ²_,m)) = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ[$m])) for m ∈ 2:M]...)
            # $([ :($(Symbol(:ρᵗ_,    m)) = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ[$m]))        for m ∈ 2:M]...)
            # $([ :($(Symbol(:ρ²ᵗ_,   m)) = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ[$m]*A.spacing.ρᵗ[$m])) for m ∈ 2:M]...)
            # $([ :($(Symbol(:ρ⁻¹2_, m)) = SIMDPirates.vbroadcast($V, AR.τ[$m] / A.ρ)) for m ∈ 2:M]...)
            @inbounds ρ²ᵗtup = $(Expr(:tuple, [ :(A.spacing.ρᵗ[$m]*A.spacing.ρᵗ[$m]) for m ∈ 1:M-1]...))
            @inbounds ρ⁻¹2tup = $(Expr(:tuple, [ :(A.τ[$m] / A.ρ) for m ∈ 1:M-1]...))

            @inbounds for n ∈ 0:$(N-1)
                δBC_last = B[1 + $M*n]
                A_1 = δBC_last
                # product[1 + $M*n] = A_m
                δ_last = A_1

                # δ_previous = one($T)
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                for m ∈ 1:$(M-1)
                    δBC_next = B[m+1 + $M*n]
                    ρᵗB = SIMDPirates.vmul( A.spacing.ρᵗ[m], δBC_last)
                    A_m = SIMDPirates.vmul( A.spacing.rinvOmρ²ᵗ[m], SIMDPirates.vsub(δBC_next, ρᵗB ) )

                    δ_next = SIMDPirates.vmul( A_m, A.spacing.rinvOmρ²ᵗ[m] )
                    qf = SIMDPirates.vmuladd(A_m, A_m, qf)
                    product[m + $M*n] = SIMDPirates.vfnmadd(A.spacing.ρᵗ[m], δ_next, δ_last)

                    ∂out = SIMDPirates.vmuladd(SIMDPirates.vmul(δ_next, ρ⁻¹2tup[m]), SIMDPirates.vfnmadd(δ_next, ρ²ᵗtup[m], ρᵗB), ∂out)
                    δBC_last = δBC_next
                    δ_last = δ_next
                end
                product[$M + $M*n] = δ_last
            end
        end)
    else
        throw("Spacing type $S not yet supported. Perhaps try caching intermediate results: `cache(AR_matrix)`")
    end

    # push!(q.args, :(qf, -0.5*(SIMDPirates.vsum(∂out)), ConstantFixedSizeMatrix(product)))
    # relying on inlining to avoid allocations.
    # push!(q.args, :(qf, ∂out, product))

    q
end

"""
D = A(ρ) * (B - C)
qf = self_dot(D)

returns: qf, -0.5* ∂qf/∂ρ, A(ρ)' * D
"""
@generated function ∂selfcrossmul_and_quadform(
            A::AbstractAutoregressiveMatrix{T,R,S},
            B::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}}
        # ) where {T,R,M,N,W,S}
        ) where {T,R,M,N,W,S}

    quote
        $(Expr(:meta, :inline))
        product = FixedSizeMatrix{$M,$N,NTuple{$W,Core.VecElement{$T}}}(undef)
        $(∂selfcrossmul_and_quadform_quote(M,N,W,T,R,S))
        qf, ∂out, product
    end
end
@generated function ∂selfcrossmul_and_quadform!(
            product::PaddedMatrices.AbstractMutableFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            A::AbstractAutoregressiveMatrix{T,R,S},
            B::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}}
        # ) where {T,R,M,N,W,S}
        ) where {T,R,M,N,W,S}

    quote
        $(Expr(:meta, :inline))
        $(∂selfcrossmul_and_quadform_quote(M,N,W,T,R,S))
        qf, ∂out#, product
    end
end

function ∂selfcrossmuldiff_and_quadform_quote(M,N,W,T,R,S)
    V = NTuple{W,Core.VecElement{T}}
    q = quote
        # $(Expr(:meta, :inline))
        qf = SIMDPirates.vbroadcast($V, zero($T))
        # product = FixedSizeMatrix{$M,$N,$V}(undef)
        ∂out = SIMDPirates.vbroadcast($V, zero($T))
    end
    if R <: AbstractUnitRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.ρ)
            ρ⁻¹2 = SIMDPirates.vbroadcast($V, one($T) / A.ρ)
            ρ²ᵗ = SIMDPirates.vbroadcast($V, A.ρ*A.ρ)
        end)
    elseif R <: AbstractRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ)
            ρ⁻¹2 = SIMDPirates.vbroadcast($V, $T(step(A.τ)) / A.ρ)
            ρ²ᵗ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ*A.spacing.ρᵗ)
        end)
    # else
    end
    if S <: AbstractEvenSpacing
        push!(q.args, quote
            @inbounds for n ∈ 0:$(N-1)
                δBC_last = SIMDPirates.vsub(B[1 + $M*n], C[1 + $M*n])
                A_1 = δBC_last
                # product[1 + $M*n] = A_m
                δ_last = A_1

                # δ_previous = one($T)
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                for m ∈ 1:$(M-1)
                    δBC_next = SIMDPirates.vsub(B[m+1 + $M*n], C[m+1 + $M*n])
                    ρᵗB = SIMDPirates.vmul( vρ, δBC_last)
                    A_m = SIMDPirates.vmul( vrinvOmρ²ᵗ, SIMDPirates.vsub(δBC_next, ρᵗB ) )

                    δ_next = SIMDPirates.vmul( A_m, vrinvOmρ²ᵗ )
                    qf = SIMDPirates.vmuladd(A_m, A_m, qf)
                    product[m + $M*n] = SIMDPirates.vfnmadd(vρ, δ_next, δ_last)

                    ∂out = SIMDPirates.vmuladd(SIMDPirates.vmul(δ_next, ρ⁻¹2), SIMDPirates.vfnmadd(δ_next, ρ²ᵗ, ρᵗB), ∂out)
                    δBC_last = δBC_next
                    δ_last = δ_next
                end
                product[$M + $M*n] = δ_last
            end
        end)
    elseif S <: CachedUnevenSpacing
        push!(q.args, quote
            # $([ :($(Symbol(:riOmρ²_,m)) = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ[$m])) for m ∈ 2:M]...)
            # $([ :($(Symbol(:ρᵗ_,    m)) = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ[$m]))        for m ∈ 2:M]...)
            # $([ :($(Symbol(:ρ²ᵗ_,   m)) = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ[$m]*A.spacing.ρᵗ[$m])) for m ∈ 2:M]...)
            # $([ :($(Symbol(:ρ⁻¹2_, m)) = SIMDPirates.vbroadcast($V, AR.τ[$m] / A.ρ)) for m ∈ 2:M]...)
            @inbounds ρ²ᵗtup = $(Expr(:tuple, [ :(A.spacing.ρᵗ[$m]*A.spacing.ρᵗ[$m]) for m ∈ 1:M-1]...))
            @inbounds ρ⁻¹2tup = $(Expr(:tuple, [ :(A.τ[$m] / A.ρ) for m ∈ 1:M-1]...))

            @inbounds for n ∈ 0:$(N-1)
                δBC_last = SIMDPirates.vsub(B[1 + $M*n], C[1 + $M*n])
                A_1 = δBC_last
                # product[1 + $M*n] = A_m
                δ_last = A_1

                # δ_previous = one($T)
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                for m ∈ 1:$(M-1)
                    δBC_next = SIMDPirates.vsub(B[m+1 + $M*n], C[m+1 + $M*n])
                    ρᵗB = SIMDPirates.vmul( A.spacing.ρᵗ[m], δBC_last)
                    A_m = SIMDPirates.vmul( A.spacing.rinvOmρ²ᵗ[m], SIMDPirates.vsub(δBC_next, ρᵗB ) )

                    δ_next = SIMDPirates.vmul( A_m, A.spacing.rinvOmρ²ᵗ[m] )
                    qf = SIMDPirates.vmuladd(A_m, A_m, qf)
                    product[m + $M*n] = SIMDPirates.vfnmadd(A.spacing.ρᵗ[m], δ_next, δ_last)

                    ∂out = SIMDPirates.vmuladd(SIMDPirates.vmul(δ_next, ρ⁻¹2tup[m]), SIMDPirates.vfnmadd(δ_next, ρ²ᵗtup[m], ρᵗB), ∂out)
                    δBC_last = δBC_next
                    δ_last = δ_next
                end
                product[$M + $M*n] = δ_last
            end
        end)
    else
        throw("Spacing type $S not yet supported. Perhaps try caching intermediate results: `cache(AR_matrix)`")
    end

    # push!(q.args, :(qf, -0.5*(SIMDPirates.vsum(∂out)), ConstantFixedSizeMatrix(product)))
    # relying on inlining to avoid allocations.
    # push!(q.args, :(qf, ∂out, product))

    q
end

"""
D = A(ρ) * (B - C)
qf = self_dot(D)

returns: qf, -0.5* ∂qf/∂ρ, A(ρ)' * D
"""
@generated function ∂selfcrossmuldiff_and_quadform(
            A::AbstractAutoregressiveMatrix{T,R,S},
            B::AbstractFixedSizeMatrix{M,N,VoT},
            C::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}}
        # ) where {T,R,M,N,W,VoT,S}
        ) where {T,R,M,N,W,S,VoT}

    quote
        $(Expr(:meta, :inline))
        product = FixedSizeMatrix{$M,$N,NTuple{$W,Core.VecElement{$T}}}(undef)
        $(∂selfcrossmuldiff_and_quadform_quote(M,N,W,T,R,S))
        qf, ∂out, product
    end
end
@generated function ∂selfcrossmuldiff_and_quadform!(
            product::PaddedMatrices.AbstractMutableFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            A::AbstractAutoregressiveMatrix{T,R,S},
            B::AbstractFixedSizeMatrix{M,N,VoT},
            C::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}}
        # ) where {T,R,M,N,W,S,VoT}
        ) where {T,R,M,N,W,VoT,S}

    quote
        $(Expr(:meta, :inline))
        $(∂selfcrossmuldiff_and_quadform_quote(M,N,W,T,R,S))
        qf, ∂out#, product
    end
end

@generated function mul_and_quadform(
            A::AbstractAutoregressiveMatrix{T,R,S},
            B::AbstractFixedSizeMatrix{M,N,T}
        ) where {T <: Number,R,M,N,S}

    W, Wshift = VectorizationBase.pick_vector_width_shift(M-1, T)
    V = NTuple{W,Core.VecElement{T}}
    q = quote
        $(Expr(:meta, :inline))
        Base.Cartesian.@nexprs $N n -> qf_n = SIMDPirates.vbroadcast($V, zero($T))
        product = FixedSizeMatrix{$M,$N,$T}(undef)
    end
    if R <: AbstractUnitRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.ρ)
        end)
    elseif R <: AbstractRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ)
        end)
    # else
    end
    if S <: AbstractEvenSpacing
        push!(q.args, quote
            A_1 = B[1]
            product[1] = A_1
            qf = A_1 * A_1
            Base.Cartesian.@nexprs $(N-1) n -> begin
                A_1 = B[1 + $M*n]
                product[1 + $M*n] = A_1
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
            end
            # Base.Cartesian.@nexprs $(N-1) n -> qf_n = SIMDPirates.vmul(A_1, A_1)
            @vectorize $T for m ∈ 1:$(M-1)
                $([quote
                        A_m = SIMDPirates.vmul( vrinvOmρ², SIMDPirates.vfnmadd( vρ, B[ m + $(M*n-1) ], B[ m + $(M*n) ] ) )
                        product[m + $(M*n)] = A_m
                        $(Symbol(:qf_,n)) = SIMDPirates.vmuladd(A_m, A_m, $(Symbol(:qf_,n)))
                    end for n ∈ 1:N]...)
            end
            # for n ∈ 0:$(N-1)
            #     A_1 = B[1 + $M*n]
            #     product[1 + $M*n] = A_1
            #     qf = SIMDPirates.vmuladd(A_1, A_1, qf)
            #     @vectorize $T for m ∈ 1:$(M-1)
            #         A_m = SIMDPirates.vmul( vrinvOmρ², SIMDPirates.vfnmadd( vρ, B[ m - 1 + $M*n ], B[ m + $M*n ] ) )
            #         product[m + $M*n] = A_m
            #         qf = SIMDPirates.vmuladd(A_m, A_m, qf)
            #     end
            # end
        end)
    elseif S == CachedUnevenSpacing
        push!(q.args, quote
            A_1 = B[1]
            product[1] = A_1
            qf = A_1 * A_1
            Base.Cartesian.@nexprs $(N-1) n -> begin
                A_1 = B[1 + $M*n]
                product[1 + $M*n] = A_1
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
            end
            rinvOmρ²ᵗ = A.spacing.rinvOmρ²ᵗ
            ρᵗ = A.spacing.ρᵗ
            @vectorize $T for m ∈ 1:$(M-1)
                vρ = ρᵗ[m]
                vrinvOmρ² = rinvOmρ²ᵗ[m]
                $([quote
                    A_m = SIMDPirates.vmul( vrinvOmρ², SIMDPirates.vfnmadd( vρ, B[ m + $(M*n-1) ], B[ m + $(M*n) ] ) )
                    product[m + $(M*n)] = A_m
                    $(Symbol(:qf_,n)) = SIMDPirates.vmuladd(A_m, A_m, $(Symbol(:qf_,n)))
                end for n ∈ 1:N]...)
            end
        end)
    end

    push!(q.args, :(qf, ConstantFixedSizeMatrix(product)))

    q

end
function ∂mul_and_quadform(
            A::AbstractAutoregressiveMatrix{T,R,S},
            B::AbstractFixedSizeMatrix{M,N,T}
        ) where {T <: Number,R,M,N,W,S}
    V = NTuple{W,Core.VecElement{T}}
    q = quote
        $(Expr(:meta, :inline))
        qf = SIMDPirates.vbroadcast($V, zero($T))
        product = FixedSizeMatrix{$M,$N,$V}(undef)
    end
    if R <: AbstractUnitRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.ρ)
            ρ⁻¹2 = SIMDPirates.vbroadcast($V, 2 / A.ρ)
        end)
    elseif R <: AbstractRange
        push!(q.args, quote
            vrinvOmρ² = SIMDPirates.vbroadcast($V, A.spacing.rinvOmρ²ᵗ)
            vρ = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ)
            ρ⁻¹2 = SIMDPirates.vbroadcast($V, T(2*step(AR.τ)) / A.ρ)
        end)
    # else
    end
    if S <: AbstractEvenSpacing
        push!(q.args, quote
            ρ²ᵗ = SIMDPirates.vbroadcast($V, $( R <: AbstractUnitRange ? :(A.ρ*A.ρ) : :(A.spacing.ρᵗ*A.spacing.ρᵗ)  ) )
            for n ∈ 0:$(N-1)
                A_1 = B[1 + $M*n]
                product[1 + $M*n] = A_1
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                $([quote
                        ρᵗB = SIMDPirates.vmul(vρ, B[$(m-1)+$M*n ])
                        $(Symbol(:A_,m)) = SIMDPirates.vmul( vrinvOmρ², SIMDPirates.vsub(B[$m+$M*n], ρᵗB ) )
                        product[$m + $M*n] = $(Symbol(:A_,m))

                        qf = SIMDPirates.vmuladd($(Symbol(:A_,m)), $(Symbol(:A_,m)), qf)

                        δ = SIMDPirates.vmul( $(Symbol(:A_,m)), vrinvOmρ² )
                        ∂out = SIMDPirates.vmuladd(SIMDPirates.vmul(δ, ρ⁻¹2), SIMDPirates.vfmsub(δ, ρ²ᵗ, ρᵗB), ∂out)

                    end for m ∈ 2:M]...)
            end
        end)
    elseif S == CachedUnevenSpacing
        push!(q.args, quote
            $([ :($(Symbol(:riOmρ²_,m)) = SIMDPirates.vbroadcast($V, A.rinvOmρ²ᵗ[$m])) for m ∈ 2:M]...)
            $([ :($(Symbol(:ρᵗ_,    m)) = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ[$m]))        for m ∈ 2:M]...)
            $([ :($(Symbol(:ρ²ᵗ_,   m)) = SIMDPirates.vbroadcast($V, A.spacing.ρᵗ[$m]*A.spacing.ρᵗ[$m])) for m ∈ 2:M]...)
            $([ :($(Symbol(:ρ⁻¹2_, m)) = SIMDPirates.vbroadcast($V, T(2*AR.τ[$m]) / A.ρ)) for m ∈ 2:M]...)
            for n ∈ 0:$(N-1)
                A_1 = B[1 + $M*n]
                product[1 + $M*n] = A_1
                qf = SIMDPirates.vmuladd(A_1, A_1, qf)
                $([quote
                        ρᵗB = SIMDPirates.vmul($(Symbol(:ρᵗ, m)), B[$(m-1)+$M*n ])
                        $(Symbol(:A_,m)) = SIMDPirates.vmul( $(Symbol(:riOmρ²_,m)), SIMDPirates.vsub(B[$m+$M*n], ρᵗB ) )
                        product[$m + $M*n] = $(Symbol(:A_,m))

                        qf = SIMDPirates.vmuladd($(Symbol(:A_,m)), $(Symbol(:A_,m)), qf)

                        δ = SIMDPirates.vmul( $(Symbol(:A_,m)), $(Symbol(:riOmρ²_,m)) )
                        ∂out = SIMDPirates.vmuladd(SIMDPirates.vmul(δ, $(Symbol(:ρ⁻¹2_, m))), SIMDPirates.vfmsub(δ, $(Symbol(:ρ²ᵗ_,   m)), ρᵗB), ∂out)

                    end for m ∈ 2:M]...)
            end
        end)
    else
        throw("Spacing type $S not yet supported. Perhaps cache(AR_matrix) will work.")
    end

    push!(q.args, :(qf, ConstantFixedSizeMatrix(product)))

    q

end

@generated function Base.copyto!(
    A::PaddedMatrices.AbstractMutableFixedSizeMatrix{nT,nT,T,R1},
    AR::AutoregressiveMatrix{T,VT}
) where {T,nT,R1,R2,VT <: PaddedMatrices.AbstractFixedSizeVector{nT,T,R2}}
    # We will assume rho > 0
    R = min(R1,R2)
    W, Wshift = VectorizationBase.pick_vector_width_shift(nT, T)
    Wm1 = W - 1
    V = Vec{W,T}
    nTl = (nT + Wm1) & ~Wm1
    # We want it to be a multiple of 8, and we don't want nTl bleeding over the edge
    T_size = sizeof(T)
    q = quote
        times = AR.τ
        ptr_A = pointer(A)
        ρ = AR.ρ
        # We want to use Base.log, and not fastmath log.
        logρ = $(T == Float64 ? :(ccall(:log,Float64,(Float64,),ρ))  :  :(Base.log(ρ)))
#        logρ = log2(ρ)
#        vρ = SIMDPirates.vbroadcast($V, ρ)
    end
    reps = (nTl >> Wshift) - 1
#    if false#reps > 10
    if true#reps > 10
        rem = nTl & Wm1
        inner_loop = quote
            for i ∈ tcouter:$((nTl>>Wshift)-1)
                vtimes_tr = SIMDPirates.vload($V, ptr_time + i * $(W))
                vδt = SIMDPirates.vabs(SIMDPirates.vsub(vtimes_tr, vtime_tc))
                AR_offset = (tc-1)*$(T_size*R) + i*$(W*T_size)
                SIMDPirates.vstore!(ptr_A + AR_offset, SLEEFPirates.exp(SIMDPirates.vmul(vδt,vlogρ)))
            end
        end
        if rem > 0
            mask = VectorizationBase.mask(Val(W), rem)
            push!(inner_loop.args, quote
                  vtimes_tr = SIMDPirates.vload($V, ptr_time + i * $(W), $mask)
                  vδt = SIMDPirates.vabs(SIMDPirates.vsub(vtimes_tr, vtime_tc))
                  AR_offset = (tc-1)*$(T_size*R) + i*$(W*T_size)
                  SIMDPirates.vstore!(ptr_A + AR_offset, SLEEFPirates.exp(SIMDPirates.vmul(vδt,vlogρ)), $mask)
                  end)
        end
        loop_q = quote
            ptr_time = VectorizationBase.vpointer(times)
            vlogρ = SIMDPirates.vbroadcast($V, logρ)
            @inbounds for tcouter ∈ 0:$((nTl>>Wshift)-1)
                for tcinner ∈ 1:$W
                    Wtcouter = $W*tcouter
                    tc = tcinner + Wtcouter
                    tc > $nT && break
                    for tr ∈ 1:Wtcouter
                        A[tr,tc] = A[tc,tr]
                    end
                    vtime_tc = SIMDPirates.vbroadcast($V, times[tc])
                    $inner_loop
                end
            end
        end
    else
        loop_q = quote
            @inbounds for c ∈ 0:$(nT-1)
                time_c = times[c+1]
                LoopVectorization.@vvectorize $T for r ∈ 1:$R
                    times_r = times[r]
                    δt = abs(times_r - time_c)
                    A[r + $R1*c] = exp(δt * logρ)
                end
            end
        end
    end
    push!(q.args, loop_q)
    q
end


text:
#using Knet,ArgParse,Main

for p in ("Knet","ArgParse")
    Pkg.installed(p) == nothing && Pkg.add(p)
end

using Knet,ArgParse,Main
#some global functions
z(P,i,x) = P[1][i]*x + P[2][i]
sig_der(x) = sigm(x)*(1-sigm(x))

#Analytical Solution of the Eq
diff_sol(x) = exp(-(x^2)/2)/(1+x+(x)^3) + x^2

#The phi trial solution
phi(P,x) = 1 + x*(predict(P,x))

#The Differential Equation function for Example 1
f(P,x) = x^3 + 2*x + (x^2)*((1+3*(x^2))/(1+x+(x^3))) - (phi(P,x))*(x + ((1+3*(x^2))/(1+x+x^3)))

 x_tr = linspace(0.0,1,10)
 x_tst = linspace(0.1,1,10)

function train(P, prms, data; epochs =100, iters=10000)
    #print(size(P),size(g),size(prms))
      for epoch=1:epochs
      	      #i = 0
        for (x,y) in data

          #i = i+1
	        g = lossgradient(P, x,y)
          #print(size(P),size(g),size(prms))
          update!(P, g, prms)
          #println(P[1][1],P[2][1],P[3][1],P[4][1])
          if (iters -= 1) <= 0
            return P
          end
        end
      end

    return P
end

function predict(P,x)
  w1, b1, w2, b2 = P
  h = sigm(w1 * x .+ b1)
  return w2 * h .+ b2
end

# function loss(P,x)
#
#     dNx = sum([P[3][i]*P[1][i]*sig_der(P[1][i]*x + P[2][i]) for i = 1:10])
#     ((predict(P,x) + (x*dNx) - f(P,x))[1])^2
# end

function loss(w,x,ygold)
    sumabs2(ygold - predict(w,x)) / size(ygold,2)
end

lossgradient = grad(loss)

function init_params(;ftype=Float32,atype=KnetArray)
    #P = Vector{Vector{Float32}}(4)
    P = Array{Any}(4)
    P[1] = randn(Float32,10,1)
    P[2] = zeros(Float32,10,1)
    P[3] = randn(Float32,1,10)
    P[4] = zeros(Float32,1,1)
    P = map(x -> convert(atype, x), P)

    return P
end

function accuracy(w, dtst, pred=predict)
  ninstance = nerror = 0
  for (x, ygold) in dtst
      ypred = pred(w, x)[1]
      #println(" X: ",x," Prediction: ",ypred[1]," True: ",ygold)
      nerror += abs(ygold - ypred[1])
      ninstance += 1
  end
  return (nerror/ninstance)
end

function params(ws, o)
	prms = Any[]

	for i=1:length(ws)
		w = ws[i]
		if o[:optim] == "Sgd"
			prm = Sgd(;lr=o[:lr])
		elseif o[:optim] == "Momentum"
			prm = Momentum(lr=o[:lr], gamma=o[:gamma])
		elseif o[:optim] == "Adam"
			prm = Adam(lr=o[:lr], beta1=o[:beta1], beta2=o[:beta2], eps=o[:eps])
		elseif o[:optim] == "Adagrad"
			prm = Adagrad(lr=o[:lr], eps=o[:eps])
		elseif o[:optim] == "Adadelta"
			prm = Adadelta(lr=o[:lr], rho=o[:rho], eps=o[:eps])
		elseif o[:optim] == "Rmsprop"
			prm = Rmsprop(lr=o[:lr], rho=o[:rho], eps=o[:eps])
		else
			error("Unknown optimization method!")
		end
		push!(prms, prm)
	end

	return prms
end

function batch(x, y, batchsize; atype=Array{Float32})
    data = Any[]
    x_ = convert(atype,x)
    y_ = convert(atype,y)
    for i=1:10
        push!(data, (x_[i],y_[i]))
    end
    return data
end

function main(args=ARGS)
    s = ArgParseSettings()
    s.description="optimizers.jl (c) Ozan Arkan Can and Deniz Yuret, 2016. Demonstration of different sgd based optimization methods using LeNet."
    s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        ("--seed"; arg_type=Int; default=-1; help="random number seed: use a nonnegative int for repeatable results")
        ("--batchsize"; arg_type=Int; default=1; help="minibatch size")
        ("--lr"; arg_type=Float64; default=0.1; help="learning rate")
	("--eps"; arg_type=Float64; default=1e-6; help="epsilon parameter used in adam, adagrad, adadelta")
	("--gamma"; arg_type=Float64; default=0.95; help="gamma parameter used in momentum")
	("--rho"; arg_type=Float64; default=0.9; help="rho parameter used in adadelta and rmsprop")
	("--beta1"; arg_type=Float64; default=0.9; help="beta1 parameter used in adam")
	("--beta2"; arg_type=Float64; default=0.95; help="beta2 parameter used in adam")
        ("--epochs"; arg_type=Int; default=10; help="number of epochs for training")
        ("--iters"; arg_type=Int; default=6000; help="number of updates for training")
	("--optim"; default="Sgd"; help="optimization method (Sgd, Momentum, Adam, Adagrad, Adadelta, Rmsprop)")
    end
    println(s.description)
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    println("opts=",[(k,v) for (k,v) in o]...)
    o[:seed] > 0 && srand(o[:seed])
    gpu() >= 0 || error("LeNet only works on GPU machines.")

    #isdefined(MNIST,:xtrn) || MNIST.loaddata()
    global dtrn = batch(x_tr, map(diff_sol,x_tr), o[:batchsize])
    global dtst = batch(x_tst, map(diff_sol,x_tst), o[:batchsize])
    # dtrn = Any[]
    # dtst = Any[]
    # for i=1:10
    #     push!(dtrn, (x_tr[i],map(diff_sol,x_tr)[i]))
    #     push!(dtst, (x_tst[i],map(diff_sol,x_tst)[i]))
    # end
    # #print("params initiated nad trn, tst pushed")
    w = init_params()
    prms = params(w, o)

    # log = Any[]
    # report(epoch)=push!(log, (:epoch,epoch,:trn,accuracy(w,dtrn,predict),:tst,accuracy(w,dtst,predict)))
    report(epoch)=println((:epoch,epoch,:trn,accuracy(w,dtrn,predict),:tst,accuracy(w,dtst,predict)))
    report(0)
    iters = o[:iters]
    #print(size(w),size(g),size(prms))

    @time for epoch=1:o[:epochs]

	    train(w, prms, dtrn; epochs=100, iters=1000)
	    report(epoch)
      # println(w[1][1])
      # println(w[2][1])
      # println(w[3][1])
      # println(w[4][1])
        (iters -= length(dtrn)) <= 0 && break
    end

    # for t in log; println(t); end
    return w
end

main()

ext: jl


text:
using SupervisedLearning
using UnicodePlots

import EmpiricalRisks
import Regression

# using RDatasets

# dataFrame = dataset("datasets", "mtcars")
# formula = AM ~ DRat + WT

# data = dataSource(formula, dataFrame, SignedClassEncoding)
# model = Classifier.LogisticRegression(l2_coef = .001)
# solver = Solver.GradientDescent()

# train!(model, data, solver, max_iter = 1000, break_every = 10) do
#   println("Iterations: $(iterations(model))")
#   println("Cost: $(cost(model))")
# end

# print(lineplot(trainingCurve(model)...))

# d = 100
# n = 10000
# w = randn(d+1)
# X = randn(d, n)
# t = sign(X'w[1:d] + w[d+1] + 0.01 * randn(n))

# f1 = Regression.RegRiskFun(EmpiricalRisks.riskmodel(EmpiricalRisks.AffinePred(d, 1.), EmpiricalRisks.LogisticLoss()),
#                           Regression.SqrL2Reg(1.0e-4),
#                           X, t)

# buf = zeros(size(t))
# f(θ) = Regression.value!(buf, f1, θ)
# function g!{T<:FloatingPoint}(θ::StridedArray{T}, storage::StridedArray{T})
#   v_risk, _ = Regression.value_and_addgrad!(buf, f1.rmodel, zero(T), storage, one(T), θ, f1.X, f1.Y)
#   v_regr, _ = Regression.value_and_addgrad!(f1.reg, one(T), storage, one(T), θ)
#   storage
# end
# function f_and_g!{T<:FloatingPoint}(θ::StridedArray{T}, storage::StridedArray{T})
#   v_risk, _ = Regression.value_and_addgrad!(buf, f1.rmodel, zero(T), storage, one(T), θ, f1.X, f1.Y)
#   v_regr, _ = Regression.value_and_addgrad!(f1.reg, one(T), storage, one(T), θ)
#   v_risk + v_regr
# end

# using Optim

# dall = DifferentiableFunction(f, g!, f_and_g!)

# println("Optim")
# theta1 = zeros(d+1)
# @time theta1 = optimize(dall, theta1, method = :gradient_descent, iterations = 100, ftol = 1.0e-6, xtol = 1.0e-8, grtol = 1.0e-8)
# theta1 = zeros(d+1)
# @time theta1 = optimize(dall, theta1, method = :gradient_descent, iterations = 10000, ftol = 1.0e-6, xtol = 1.0e-8, grtol = 1.0e-8)

# println("Regression")
# theta2 = zeros(d+1)
# @time sol = Regression.solve!(Regression.GD(), f1, theta2,
#                     Regression.Options(maxiter = 100, ftol = 1.0e-6, xtol = 1.0e-8, grtol = 1.0e-8),
#                     Regression.no_op)
# theta2 = zeros(d+1)
# @time sol = Regression.solve!(Regression.GD(), f1, theta2,
#                     Regression.Options(maxiter = 10000, ftol = 1.0e-6, xtol = 1.0e-8, grtol = 1.0e-8),
#                     Regression.no_op)

# println()
# println(theta1)
# println()
# println(theta2)
# println(sol)
# data = dataSource(X, t, SignedClassEncoding(["no","yes"]))
# model = Classifier.LogisticRegression(l2_coef = .01)
# solver = Solver.GradientDescent()

# println("State: $(state(model))")
# println("Iterations: $(iterations(model))")
# println("Cost: $(cost(model))")

# @time train!(model, data, solver, max_iter = 100, break_every = 500)
# @time train!(model, data, solver, max_iter = 10000, break_every = 500)
# train!(model, data, solver, max_iter = 10000, break_every = 500) do
#   println("Iterations: $(iterations(model))")
#   println("Cost: $(cost(model))")
# end

# train!(model, data, solver, max_iter = 1000, break_every = 100) do
#   println("Iterations: $(iterations(model))")
# end
# #lineplot(trainingCurve(model)...)
# println()
# print(lineplot(trainingCurve(model)...))
# println("State: $(state(model))")
# println("Iterations: $(iterations(model))")
# println("Fitness: $(fitness(model))")
# println("Cost: $(cost(model))")
#@test_approx_eq w model._risk_learner._coefficients

# using EmpiricalRisks
# using ArrayViews

# pm = AffinePred(d, 1.)
# @time s = predict(pm, w, X)
# @time s = predict(pm, w, X)

# r = zeros(n)
# @time predict!(r, pm, w, X)
# @time predict!(r, pm, w, X)

# println("Loop with predict")
# @time for i = 1:1000
#   s = predict(pm, w, X)
# end

# println("Loop with predict!")
# @time for i = 1:1000
#   s = predict!(r, pm, w, X)
# end

# s = predict(pm, w, X)
# @test_approx_eq s r

# println("Risk models")
# rmodel = riskmodel(pm, LogisticLoss())
# value(rmodel, w, X, t)
# @time value(rmodel, w, X, t)

# rmodel = riskmodel(pm, LogisticLoss())
# value!(r, rmodel, w, X, t)
# @time value!(r, rmodel, w, X, t)

# println("Loop with value")
# @time for i = 1:100
#   s = value(rmodel, w, X, t)
# end

# println("Loop with value!")
# @time for i = 1:100
#   s = value!(r, rmodel, w, X, t)
# end

# data = dataSource(X, t, SignedClassEncoding(["no","yes"]))
# model = Classifier.LogisticRegression(l2_coef = .001)
# solver = Solver.GradientDescent()

# println("State: $(state(model))")
# println("Iterations: $(iterations(model))")
# println("Cost: $(cost(model))")

# train!(model, data, solver, max_iter = 1000, break_every = 100)
# @time train!(model, data, solver, max_iter = 1000, break_every = 100)

# # train!(model, data, solver, max_iter = 1000, break_every = 100) do
# #   println("Iterations: $(iterations(model))")
# # end
# #lineplot(trainingCurve(model)...)
# #@time print(lineplot(trainingCurve(model)...))
# println("State: $(state(model))")
# println("Iterations: $(iterations(model))")
# println("Fitness: $(fitness(model))")
# println("Cost: $(cost(model))")


text:
# Create example data
angle = LinRange(0, 2π, 40)
radius = LinRange(0, 10, 20)
z = sin.(angle') .* cos.(radius)
polarheatmap(z)


text:
__precompile__()
module BoundaryIntegralMethod
using Distributed

export Boundary, SimpleBoundary, Circle, DiscretizedBoundary, Point2D, nExt,
      nInt, isInside, nIdxRefr, boundaryPoint, boundaryPointDeriv, normVector, boundaryIdx,
      scaledParameter, absDiff, boundaryPointVector, A, A′, residualInverseIteration, NoeckelQuadrupoleBoundary,
      exportTotalField, ψField, G, ∂νG, arcLen, curvature, unitNormVector

include("boundaries.jl")
include("bimCalculations.jl")
include("modeSearch.jl")
include("totalField.jl")
include("DirectSweep.jl")
include("Follow.jl")
include("Husimi.jl")
include("Utilities.jl")

end # module
# @everywhere using BoundaryIntegralMethod
# @everywhere using DistributedArrays


text:
if Pkg.installed("JuMP") == nothing
    Pkg.add("JuMP")
end

if Pkg.installed("NLPModels") == nothing
    Pkg.add("NLPModels")
end

if Pkg.installed("Ipopt") == nothing
    Pkg.add("Ipopt") #nécessite un compilateur fortran (gfortran)
end

#if Pkg.installed("PyPlot") == nothing
#    Pkg.add("PyPlot") #nécessite Python d'installer
#end

if Pkg.installed("BenchmarkProfiles") == nothing
    Pkg.add("BenchmarkProfiles")
end

if Pkg.installed("CUTEst") == nothing
    Pkg.add("CUTEst") #nécessite gfortran, wget, gsl
end

if Pkg.installed("OptimizationProblems") == nothing
 #using OptimizationProblems
 Pkg.clone("https://github.com/JuliaSmoothOptimizers/OptimizationProblems.jl.git")
end

if Pkg.installed("AmplNLReader") == nothing
 Pkg.clone("https://github.com/JuliaSmoothOptimizers/AmplNLReader.jl.git")
 Pkg.build("AmplNLReader") #nécessite cmake
end


text:
#############################################################################
# dotsort.jl
# dotsort(a,b) computes dot(sort(a), sort(b))
# All expressions and atoms are subtpyes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################

# This atom computes dot(sort(x), sort(w)), where w is constant
# for example, if w = [1 1 1 0 0 0 ... 0], it computes the sum of the 3 largest elements of x
struct DotSortAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int,Int}
    w::Value

    function DotSortAtom(x::AbstractExpr, w::Value)
        if sign(x) == ComplexSign()
            error("Argument should be real instead it is $(sign(x))")
        else
            if !(length(w) == length(x))
                error("x and w must be the same size")
            end
            children = (x,)
            vecw = reshape(w, length(x))
            return new(:dotsort, hash((children, vecw)), children, (1, 1), vecw)
        end
    end
end

function sign(x::DotSortAtom)
    if all(x.w .>= 0)
        return sign(x.children[1])
    elseif all(x.w .<= 0)
        return sign(x.children[1])
    else
        return NoSign()
    end
end

function monotonicity(x::DotSortAtom)
    if all(x.w .>= 0)
        return (Nondecreasing(),)
    else
        return (NoMonotonicity(),)
    end
end

function curvature(x::DotSortAtom)
    return ConvexVexity()
end

function evaluate(x::DotSortAtom)
    return sum(
        sort(vec(evaluate(x.children[1])), rev = true) .*
        sort(vec(x.w), rev = true),
    )
end

function conic_form!(x::DotSortAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, x)
        y = x.children[1]
        w = x.w
        sy = size(y)
        if sy[1] > 1 && sy[2] > 1
            y = vec(y)
        end
        mu = Variable(size(y))
        nu = Variable(size(y))
        onesvec = ones(size(y))
        # given by the solution to
        # minimize sum(mu) + sum(nu)
        # subject to y*w' <= onesvec*nu' + mu*onesvec'
        objective = conic_form!(sum(mu) + sum(nu), unique_conic_forms)
        conic_form!(y * w' <= onesvec * nu' + mu * onesvec', unique_conic_forms)
        cache_conic_form!(unique_conic_forms, x, objective)
    end
    return get_conic_form(unique_conic_forms, x)
end

dotsort(a::AbstractExpr, b::Value) = DotSortAtom(a, b)
dotsort(b::Value, a::AbstractExpr) = DotSortAtom(a, b)


text:
# Base.summary(io::IO, x::CachedNode) = (print(io, "(Tracked) "); summary(io, x.output))
# Base.show(io::IO, m::MIME"text/plain", x::CachedNode{<:Node, <:AbstractArray}) = (print(io, "(tracked) "); show(io, m, x.output))
# Base.show(io::IO, x::CachedNode) = show(io, x.output)

# printing
"""
    PrintTrait(node) -> Trait
"""
function PrintTrait end

struct ScalarVariable{T}
    value::T
end

# a trait type to forward variable to array's type tree
struct ArrayVariable{T, N, AT <: AbstractArray{T, N}} <: AbstractArray{T, N}
    value::AT
end

struct Other{NT}
    node::NT
end

# forward to default
PrintTrait(var::Number) = ScalarVariable(var)
PrintTrait(var::Value{<:Number}) = ScalarVariable(value(var))
PrintTrait(var::Value{AT}) where {T, N, AT <: AbstractArray{T, N}} = ArrayVariable{T, N, AT}(value(var))
PrintTrait(var::Value{T}) where {T} = Other(var)

Base.size(x::ArrayVariable) = size(x.value)
Base.getindex(x::ArrayVariable, i...) = getindex(x.value, i...)

Base.show(io::IO, x::Value) = show(io, "text/plain", PrintTrait(x))
Base.show(io::IO, m::MIME"text/plain", x::Value) = show(io, m, PrintTrait(x))
Base.summary(io::IO, x::Union{Variable, CachedNode, Node}) = summary(io, PrintTrait(x))
Base.show(io::IO, x::ScalarVariable) = print(io, "(Tracked) ",x.value)
Base.summary(io::IO, x::ArrayVariable) = (print(io, "(Tracked) "); summary(io, x.value))

function Base.show(io::IO, ::MIME"text/plain", x::Other{<:CachedNode})
    println(io, "CachedNode:")
    println(io, "       f: ", operator(x.node))
    println(io, "    args: ", args(x.node))
    println(io, "  kwargs: ", kwargs(x.node))
      print(io, "  output: ", x.node.output)
end

function Base.show(io::IO, ::MIME"text/plain", x::Other{<:Variable})
    println(io, "Variable:")
    println(io, "  value: ", x.node.value)

    if isdefined(x.node, :grad)
        println(io, "   grad: ", x.node.grad)
    end
end


text:
export SCS_init, SCS_solve, SCS_finish, SCS_version

function SCS_version()
    return unsafe_string(ccall((:scs_version, SCS.scs), Cstring, ()))
end

function SCS_init(data::SCSData, cone::SCSCone)
    # Initialize the info struct
    info = SCSInfo(0, convert(Int128, 0), convert(Int128, 0), 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

    p_work = ccall((:scs_init, SCS.scs), Ptr{Void},
        (Ptr{SCSData}, Ptr{SCSCone}, Ptr{SCSInfo}),
        &data, &cone, &info)

    return p_work, info
end

function SCS_solve(p_work::Ptr{Void}, data::SCSData, cone::SCSCone, info::SCSInfo, solution::SCSSolution)
    solution_ptr = pointer([solution])

    info_ptr = pointer([info])

    status = ccall((:scs_solve, SCS.scs), Int,
        (Ptr{Void}, Ptr{SCSData}, Ptr{SCSCone}, Ptr{SCSSolution}, Ptr{SCSInfo}),
        p_work, &data, &cone, solution_ptr, info_ptr)

    solution = unsafe_load(solution_ptr)
    info = unsafe_load(info_ptr)

    return status, solution, info, p_work
end

function SCS_solve(p_work::Ptr{Void}, data::SCSData, cone::SCSCone, info::SCSInfo)
    solution = SCSSolution(pointer(zeros(data.n)), pointer(zeros(data.m)), pointer(zeros(data.m)))
    return SCS_solve(p_work, data, cone, info, solution)
end

function SCS_solve(data::SCSData, cone::SCSCone, solution::SCSSolution)
    p_work, info = SCS_init(data, cone)
    return SCS_solve(p_work, data, cone, info, solution)
end

function SCS_finish(p_work::Ptr{Void})
    ccall((:scs_finish, SCS.scs), Void,
        (Ptr{Void}, ),
        p_work)
end


text:
using Base: RefValue
using LoopVectorization: parents, name, parent

struct Func
    instr::Instruction
    output::RefValue{Int} # Should support unpacking of returned objecs, whether homogonous arrays or heterogenous tuples.
    vparents::Vector{Int} # should these be Vector{vparent}, or Vector{Int} of ids ?
    loopsetid::Int#index into Model's ::Vector{LoopSet}; 0 indicates no ls
    probdistapi::Bool
    lowered::RefValue{Bool}
end
function Func(instr::Instruction, probdistapi::Bool = false, loopsetid::Int = 0)
    Func(instr, Ref{Int}(), Int[], loopsetid, probdistapi, Ref(false))
end

function Base.hash(f::Func, u::UInt)
    @unpack instr, vparents = f
    u = hash(instr, u)
    for p ∈ vparents
        u = hash(p, u)
    end
    u
end
Base.isequal(f1::Func, f2::Func) = f1.instr == f2.instr && f1.vparents == f2.vparents

function uses!(f::Func, v::Variable)
    push!(parents(f), v.varid)
    v
end
function returns!(f::Func, v::Variable)
    @assert name(v) != Symbol("##TARGET####BAR##")
    f.output[] = v.varid
    v
end
islowered(f::Func) = f.lowered[]
LoopVectorization.parents(f::Func) = f.vparents
LoopVectorization.instruction(f::Func) = f.instr
ReverseDiffExpressionsBase.InstructionArgs(f::Func) = InstructionArgs(instruction(f), length(parents(f)))

stackpointercall_expr(mod) = Expr(:(.), Expr(:(.), mod, QuoteNode(:ReverseDiffExpressions)), QuoteNote(:stack_pointer_call))

function isindexfunc(func::Func)
    func.instr.instr === :view || func.instr.instr === :getindex || func.instr.instr ∈ (:first, :second, :third, :fourth, :fifth, :sixth, :seventh, :eigth, :ninth, :last)
end

function isconstrainttransform(func::Func)
    (func.instr.instr === :constrain) || (func.instr.instr === :constrain_pullback!)
end
function notupdating(func::Func)
    func ∉ keys(ReverseDiffExpressionsBase.MAKEUPDATING)
end

# function getconstindex!(m::Model, index::Integer, vin::Variable)
# end


text:
# Diag can have either Vector storage, in which case
# it is a general Diag tensor, or scalar storage,
# in which case the diagonal has a uniform value
struct Diag{ElT,VecT} <: TensorStorage{ElT}
  data::VecT
  function Diag{ElT,VecT}(data) where {ElT,VecT<:AbstractVector{ElT}}
    return new{ElT,VecT}(data)
  end
  function Diag{ElT,ElT}(data) where {ElT}
    return new{ElT,ElT}(data)
  end
end

Diag(data::VecT) where {VecT<:AbstractVector{ElT}} where {ElT} = Diag{ElT,VecT}(data)

Diag(data::ElT) where {ElT<:Number} = Diag{ElT,ElT}(data)

function Diag{ElR}(data::AbstractVector{ElT}) where {ElR<:Number,ElT<:Number}
  return ElT == ElR ? Diag(data) : Diag(ElR.(data))
end

Diag(::Type{ElT}, n::Integer) where {ElT<:Number} = Diag(zeros(ElT, n))

Diag(x::ElT, n::Integer) where {ElT<:Number} = Diag(fill(x, n))

copy(D::Diag) = Diag(copy(data(D)))

const NonuniformDiag{ElT,VecT} = Diag{ElT,VecT} where {VecT<:AbstractVector}

const UniformDiag{ElT,VecT} = Diag{ElT,VecT} where {VecT<:Number}

getindex(D::UniformDiag, i::Int) = data(D)

function setindex!(D::UniformDiag, val, i::Int)
  return error("Cannot set elements of a uniform Diag storage")
end

Base.real(::Type{Diag{ElT,Vector{ElT}}}) where {ElT} = Diag{real(ElT),Vector{real(ElT)}}
Base.real(::Type{Diag{ElT,ElT}}) where {ElT} = Diag{real(ElT),real(ElT)}

complex(::Type{Diag{ElT,Vector{ElT}}}) where {ElT} = Diag{complex(ElT),Vector{complex(ElT)}}
complex(::Type{Diag{ElT,ElT}}) where {ElT} = Diag{complex(ElT),complex(ElT)}

# Deal with uniform Diag conversion
convert(::Type{<:Diag{ElT,VecT}}, D::Diag) where {ElT,VecT} = Diag(convert(VecT, data(D)))

# TODO: write in terms of ::Int, not inds
similar(D::NonuniformDiag) = Diag(similar(data(D)))
#similar(D::NonuniformDiag,inds) = Diag(similar(data(D),minimum(dims(inds))))
#function similar(D::Type{<:NonuniformDiag{ElT,VecT}},inds) where {ElT,VecT}
#  return Diag(similar(VecT,diaglength(inds)))
#end

similar(D::UniformDiag{ElT}) where {ElT} = Diag(zero(ElT))
similar(D::UniformDiag, inds) = similar(D)
similar(::Type{<:UniformDiag{ElT}}, inds) where {ElT} = Diag(zero(ElT))

similar(D::Diag, n::Int) = Diag(similar(data(D), n))

similar(D::Diag, ::Type{ElR}, n::Int) where {ElR} = Diag(similar(data(D), ElR, n))

# TODO: make this work for other storage besides Vector
zeros(::Type{<:NonuniformDiag{ElT}}, dim::Int64) where {ElT} = Diag(zeros(ElT, dim))
zeros(::Type{<:UniformDiag{ElT}}, dim::Int64) where {ElT} = Diag(zero(ElT))

setdata(D::Diag, ndata) = Diag(ndata)

#
# Type promotions involving Diag
# Useful for knowing how conversions should work when adding and contracting
#

function promote_rule(
  ::Type{<:UniformDiag{ElT1}}, ::Type{<:UniformDiag{ElT2}}
) where {ElT1,ElT2}
  ElR = promote_type(ElT1, ElT2)
  return Diag{ElR,ElR}
end

function promote_rule(
  ::Type{<:NonuniformDiag{ElT1,VecT1}}, ::Type{<:NonuniformDiag{ElT2,VecT2}}
) where {ElT1,VecT1<:AbstractVector,ElT2,VecT2<:AbstractVector}
  ElR = promote_type(ElT1, ElT2)
  VecR = promote_type(VecT1, VecT2)
  return Diag{ElR,VecR}
end

# This is an internal definition, is there a more general way?
#promote_type(::Type{Vector{ElT1}},
#                  ::Type{ElT2}) where {ElT1<:Number,
#                                       ElT2<:Number} = Vector{promote_type(ElT1,ElT2)}
#
#promote_type(::Type{ElT1},
#                  ::Type{Vector{ElT2}}) where {ElT1<:Number,
#                                               ElT2<:Number} = promote_type(Vector{ElT2},ElT1)

# TODO: how do we make this work more generally for T2<:AbstractVector{S2}?
# Make a similartype(AbstractVector{S2},T1) -> AbstractVector{T1} function?
function promote_rule(
  ::Type{<:UniformDiag{ElT1,VecT1}}, ::Type{<:NonuniformDiag{ElT2,Vector{ElT2}}}
) where {ElT1,VecT1<:Number,ElT2}
  ElR = promote_type(ElT1, ElT2)
  VecR = Vector{ElR}
  return Diag{ElR,VecR}
end

function promote_rule(
  ::Type{DenseT1}, ::Type{<:NonuniformDiag{ElT2,VecT2}}
) where {DenseT1<:Dense,ElT2,VecT2<:AbstractVector}
  return promote_type(DenseT1, Dense{ElT2,VecT2})
end

function promote_rule(
  ::Type{DenseT1}, ::Type{<:UniformDiag{ElT2,VecT2}}
) where {DenseT1<:Dense,ElT2,VecT2<:Number}
  return promote_type(DenseT1, ElT2)
end

# Convert a Diag storage type to the closest Dense storage type
dense(::Type{<:NonuniformDiag{ElT,VecT}}) where {ElT,VecT} = Dense{ElT,VecT}
dense(::Type{<:UniformDiag{ElT}}) where {ElT} = Dense{ElT,Vector{ElT}}

const DiagTensor{ElT,N,StoreT,IndsT} = Tensor{ElT,N,StoreT,IndsT} where {StoreT<:Diag}
const NonuniformDiagTensor{ElT,N,StoreT,IndsT} =
  Tensor{ElT,N,StoreT,IndsT} where {StoreT<:NonuniformDiag}
const UniformDiagTensor{ElT,N,StoreT,IndsT} =
  Tensor{ElT,N,StoreT,IndsT} where {StoreT<:UniformDiag}

IndexStyle(::Type{<:DiagTensor}) = IndexCartesian()

# TODO: this needs to be better (promote element type, check order compatibility,
# etc.
function convert(::Type{<:DenseTensor{ElT,N}}, T::DiagTensor{ElT,N}) where {ElT<:Number,N}
  return dense(T)
end

convert(::Type{Diagonal}, D::DiagTensor{<:Number,2}) = Diagonal(data(D))

# These are rules for determining the output of a pairwise contraction of NDTensors
# (given the indices of the output tensors)
function contraction_output_type(
  TensorT1::Type{<:DiagTensor}, TensorT2::Type{<:DenseTensor}, IndsR::Type
)
  return similartype(promote_type(TensorT1, TensorT2), IndsR)
end
function contraction_output_type(
  TensorT1::Type{<:DenseTensor}, TensorT2::Type{<:DiagTensor}, IndsR::Type
)
  return contraction_output_type(TensorT2, TensorT1, IndsR)
end

# This performs the logic that DiagTensor*DiagTensor -> DiagTensor if it is not an outer
# product but -> DenseTensor if it is
# TODO: if the tensors are both order 2 (or less), or if there is an Index replacement,
# then they remain diagonal. Should we limit DiagTensor*DiagTensor to cases that
# result in a DiagTensor, for efficiency and type stability? What about a general
# SparseTensor result?
function contraction_output_type(
  TensorT1::Type{<:DiagTensor{<:Number,N1}},
  TensorT2::Type{<:DiagTensor{<:Number,N2}},
  IndsR::Type,
) where {N1,N2}
  if ValLength(IndsR) === Val{N1 + N2}
    # Turn into is_outer(inds1,inds2,indsR) function?
    # How does type inference work with arithmatic of compile time values?
    return similartype(dense(promote_type(TensorT1, TensorT2)), IndsR)
  end
  return similartype(promote_type(TensorT1, TensorT2), IndsR)
end

# The output must be initialized as zero since it is sparse, cannot be undefined
function contraction_output(T1::DiagTensor, T2::Tensor, indsR)
  return zero_contraction_output(T1, T2, indsR)
end
contraction_output(T1::Tensor, T2::DiagTensor, indsR) = contraction_output(T2, T1, indsR)

function contraction_output(T1::DiagTensor, T2::DiagTensor, indsR)
  return zero_contraction_output(T1, T2, indsR)
end

function array(T::DiagTensor{ElT,N}) where {ElT,N}
  return array(dense(T))
end
matrix(T::DiagTensor{<:Number,2}) = array(T)
vector(T::DiagTensor{<:Number,1}) = array(T)

function Array{ElT,N}(T::DiagTensor{ElT,N}) where {ElT,N}
  return array(T)
end

function Array(T::DiagTensor{ElT,N}) where {ElT,N}
  return Array{ElT,N}(T)
end

function zeros(TensorT::Type{<:DiagTensor}, inds)
  return tensor(zeros(storagetype(TensorT), mindim(inds)), inds)
end

function zeros(TensorT::Type{<:DiagTensor}, inds::Tuple{})
  return tensor(zeros(storagetype(TensorT), mindim(inds)), inds)
end

# Needed to get slice of DiagTensor like T[1:3,1:3]
function similar(
  T::DiagTensor{<:Number,N}, ::Type{ElR}, inds::Dims{N}
) where {ElR<:Number,N}
  return tensor(similar(storage(T), ElR, minimum(inds)), inds)
end

"""
getdiagindex(T::DiagTensor,i::Int)

Get the ith value along the diagonal of the tensor.
"""
getdiagindex(T::DiagTensor{<:Number}, ind::Int) = storage(T)[ind]

"""
setdiagindex!(T::DiagTensor,i::Int)

Set the ith value along the diagonal of the tensor.
"""
setdiagindex!(T::DiagTensor{<:Number}, val, ind::Int) = (storage(T)[ind] = val)

"""
setdiag(T::UniformDiagTensor,val)

Set the entire diagonal of a uniform DiagTensor.
"""
setdiag(T::UniformDiagTensor, val) = tensor(Diag(val), inds(T))

@propagate_inbounds function getindex(
  T::DiagTensor{ElT,N}, inds::Vararg{Int,N}
) where {ElT,N}
  if all(==(inds[1]), inds)
    return getdiagindex(T, inds[1])
  else
    return zero(eltype(ElT))
  end
end
@propagate_inbounds getindex(T::DiagTensor{<:Number,1}, ind::Int) = storage(T)[ind]
@propagate_inbounds getindex(T::DiagTensor{<:Number,0}) = storage(T)[1]

# Set diagonal elements
# Throw error for off-diagonal
@propagate_inbounds function setindex!(
  T::DiagTensor{<:Number,N}, val, inds::Vararg{Int,N}
) where {N}
  all(==(inds[1]), inds) || error("Cannot set off-diagonal element of Diag storage")
  setdiagindex!(T, val, inds[1])
  return T
end
@propagate_inbounds function setindex!(T::DiagTensor{<:Number,1}, val, ind::Int)
  return (storage(T)[ind] = val)
end
@propagate_inbounds setindex!(T::DiagTensor{<:Number,0}, val) = (storage(T)[1] = val)

function setindex!(T::UniformDiagTensor{<:Number,N}, val, inds::Vararg{Int,N}) where {N}
  return error("Cannot set elements of a uniform Diag storage")
end

# TODO: make a fill!! that works for uniform and non-uniform
#fill!(T::DiagTensor,v) = fill!(storage(T),v)

function dense(::Type{<:Tensor{ElT,N,StoreT,IndsT}}) where {ElT,N,StoreT<:Diag,IndsT}
  return Tensor{ElT,N,dense(StoreT),IndsT}
end

# convert to Dense
function dense(T::TensorT) where {TensorT<:DiagTensor}
  R = zeros(dense(TensorT), inds(T))
  for i in 1:diaglength(T)
    setdiagindex!(R, getdiagindex(T, i), i)
  end
  return R
end

denseblocks(T::DiagTensor) = dense(T)

function outer!(
  R::DenseTensor{<:Number,NR}, T1::DiagTensor{<:Number,N1}, T2::DiagTensor{<:Number,N2}
) where {NR,N1,N2}
  for i1 in 1:diaglength(T1), i2 in 1:diaglength(T2)
    indsR = CartesianIndex{NR}(ntuple(r -> r ≤ N1 ? i1 : i2, Val(NR)))
    R[indsR] = getdiagindex(T1, i1) * getdiagindex(T2, i2)
  end
  return R
end

# TODO: write an optimized version of this?
function outer!(R::DenseTensor{ElR}, T1::DenseTensor, T2::DiagTensor) where {ElR}
  R .= zero(ElR)
  outer!(R, T1, dense(T2))
  return R
end

function outer!(R::DenseTensor{ElR}, T1::DiagTensor, T2::DenseTensor) where {ElR}
  R .= zero(ElR)
  outer!(R, dense(T1), T2)
  return R
end

# Right an in-place version
function outer(T1::DiagTensor{ElT1,N1}, T2::DiagTensor{ElT2,N2}) where {ElT1,ElT2,N1,N2}
  indsR = unioninds(inds(T1), inds(T2))
  R = tensor(Dense(zeros(promote_type(ElT1, ElT2), dim(indsR))), indsR)
  outer!(R, T1, T2)
  return R
end

function permutedims!(
  R::DiagTensor{<:Number,N},
  T::DiagTensor{<:Number,N},
  perm::NTuple{N,Int},
  f::Function=(r, t) -> t,
) where {N}
  # TODO: check that inds(R)==permute(inds(T),perm)?
  for i in 1:diaglength(R)
    @inbounds setdiagindex!(R, f(getdiagindex(R, i), getdiagindex(T, i)), i)
  end
  return R
end

function permutedims(
  T::DiagTensor{<:Number,N}, perm::NTuple{N,Int}, f::Function=identity
) where {N}
  R = similar(T, permute(inds(T), perm))
  g(r, t) = f(t)
  permutedims!(R, T, perm, g)
  return R
end

function permutedims(
  T::UniformDiagTensor{<:Number,N}, perm::NTuple{N,Int}, f::Function=identity
) where {N}
  R = tensor(Diag(f(getdiagindex(T, 1))), permute(inds(T), perm))
  return R
end

# Version that may overwrite in-place or may return the result
function permutedims!!(
  R::NonuniformDiagTensor{<:Number,N},
  T::NonuniformDiagTensor{<:Number,N},
  perm::NTuple{N,Int},
  f::Function=(r, t) -> t,
) where {N}
  R = convert(promote_type(typeof(R), typeof(T)), R)
  permutedims!(R, T, perm, f)
  return R
end

function permutedims!!(
  R::UniformDiagTensor{ElR,N},
  T::UniformDiagTensor{ElT,N},
  perm::NTuple{N,Int},
  f::Function=(r, t) -> t,
) where {ElR,ElT,N}
  R = convert(promote_type(typeof(R), typeof(T)), R)
  R = tensor(Diag(f(getdiagindex(R, 1), getdiagindex(T, 1))), inds(R))
  return R
end

function permutedims!(
  R::DenseTensor{ElR,N}, T::DiagTensor{ElT,N}, perm::NTuple{N,Int}, f::Function=(r, t) -> t
) where {ElR,ElT,N}
  for i in 1:diaglength(T)
    @inbounds setdiagindex!(R, f(getdiagindex(R, i), getdiagindex(T, i)), i)
  end
  return R
end

function permutedims!!(
  R::DenseTensor{ElR,N}, T::DiagTensor{ElT,N}, perm::NTuple{N,Int}, f::Function=(r, t) -> t
) where {ElR,ElT,N}
  RR = convert(promote_type(typeof(R), typeof(T)), R)
  permutedims!(RR, T, perm, f)
  return RR
end

# TODO: make a single implementation since this is
# the same as the version with the input types
# swapped.
function permutedims!!(
  R::DiagTensor{ElR,N}, T::DenseTensor{ElT,N}, perm::NTuple{N,Int}, f::Function=(r, t) -> t
) where {ElR,ElT,N}
  RR = convert(promote_type(typeof(R), typeof(T)), R)
  permutedims!(RR, T, perm, f)
  return RR
end

function _contract!!(
  R::UniformDiagTensor{ElR,NR},
  labelsR,
  T1::UniformDiagTensor{<:Number,N1},
  labelsT1,
  T2::UniformDiagTensor{<:Number,N2},
  labelsT2,
) where {ElR,NR,N1,N2}
  if NR == 0  # If all indices of A and B are contracted
    # all indices are summed over, just add the product of the diagonal
    # elements of A and B
    R = setdiag(R, diaglength(T1) * getdiagindex(T1, 1) * getdiagindex(T2, 1))
  else
    # not all indices are summed over, set the diagonals of the result
    # to the product of the diagonals of A and B
    R = setdiag(R, getdiagindex(T1, 1) * getdiagindex(T2, 1))
  end
  return R
end

function contract!(
  R::DiagTensor{ElR,NR},
  labelsR,
  T1::DiagTensor{<:Number,N1},
  labelsT1,
  T2::DiagTensor{<:Number,N2},
  labelsT2,
) where {ElR,NR,N1,N2}
  if NR == 0  # If all indices of A and B are contracted
    # all indices are summed over, just add the product of the diagonal
    # elements of A and B
    Rdiag = zero(ElR)
    for i in 1:diaglength(T1)
      Rdiag += getdiagindex(T1, i) * getdiagindex(T2, i)
    end
    setdiagindex!(R, Rdiag, 1)
  else
    min_dim = min(diaglength(T1), diaglength(T2))
    # not all indices are summed over, set the diagonals of the result
    # to the product of the diagonals of A and B
    for i in 1:min_dim
      setdiagindex!(R, getdiagindex(T1, i) * getdiagindex(T2, i), i)
    end
  end
  return R
end

function contract!(
  C::DenseTensor{ElC,NC},
  Clabels,
  A::DiagTensor{ElA,NA},
  Alabels,
  B::DenseTensor{ElB,NB},
  Blabels,
  α::Number=one(ElC),
  β::Number=zero(ElC);
  convert_to_dense::Bool=true,
) where {ElA,NA,ElB,NB,ElC,NC}
  #@timeit_debug timer "diag-dense contract!" begin
  if all(i -> i < 0, Blabels)
    # If all of B is contracted
    # TODO: can also check NC+NB==NA
    min_dim = minimum(dims(B))
    if length(Clabels) == 0
      # all indices are summed over, just add the product of the diagonal
      # elements of A and B
      # Assumes C starts set to 0
      c₁ = zero(ElC)
      for i in 1:min_dim
        c₁ += getdiagindex(A, i) * getdiagindex(B, i)
      end
      setdiagindex!(C, α * c₁ + β * getdiagindex(C, 1), 1)
    else
      # not all indices are summed over, set the diagonals of the result
      # to the product of the diagonals of A and B
      # TODO: should we make this return a Diag storage?
      for i in 1:min_dim
        setdiagindex!(
          C, α * getdiagindex(A, i) * getdiagindex(B, i) + β * getdiagindex(C, i), i
        )
      end
    end
  else
    # Most general contraction
    if convert_to_dense
      contract!(C, Clabels, dense(A), Alabels, B, Blabels, α, β)
    else
      if !isone(α) || !iszero(β)
        error(
          "`contract!(::DenseTensor, ::DiagTensor, ::DenseTensor, α, β; convert_to_dense = false)` with `α ≠ 1` or `β ≠ 0` is not currently supported. You can call it with `convert_to_dense = true` instead.",
        )
      end
      astarts = zeros(Int, length(Alabels))
      bstart = 0
      cstart = 0
      b_cstride = 0
      nbu = 0
      for ib in 1:length(Blabels)
        ia = findfirst(==(Blabels[ib]), Alabels)
        if !isnothing(ia)
          b_cstride += stride(B, ib)
          bstart += astarts[ia] * stride(B, ib)
        else
          nbu += 1
        end
      end

      c_cstride = 0
      for ic in 1:length(Clabels)
        ia = findfirst(==(Clabels[ic]), Alabels)
        if !isnothing(ia)
          c_cstride += stride(C, ic)
          cstart += astarts[ia] * stride(C, ic)
        end
      end

      # strides of the uncontracted dimensions of
      # B
      bustride = zeros(Int, nbu)
      custride = zeros(Int, nbu)
      # size of the uncontracted dimensions of
      # B, to be used in CartesianIndices
      busize = zeros(Int, nbu)
      n = 1
      for ib in 1:length(Blabels)
        if Blabels[ib] > 0
          bustride[n] = stride(B, ib)
          busize[n] = size(B, ib)
          ic = findfirst(==(Blabels[ib]), Clabels)
          custride[n] = stride(C, ic)
          n += 1
        end
      end

      boffset_orig = 1 - sum(strides(B))
      coffset_orig = 1 - sum(strides(C))
      cartesian_inds = CartesianIndices(Tuple(busize))
      for inds in cartesian_inds
        boffset = boffset_orig
        coffset = coffset_orig
        for i in 1:nbu
          ii = inds[i]
          boffset += ii * bustride[i]
          coffset += ii * custride[i]
        end
        c = zero(ElC)
        for j in 1:diaglength(A)
          # With α == 0 && β == 1
          C[cstart + j * c_cstride + coffset] +=
            getdiagindex(A, j) * B[bstart + j * b_cstride + boffset]
          # XXX: not sure if this is correct
          #C[cstart+j*c_cstride+coffset] += α * getdiagindex(A, j)* B[bstart+j*b_cstride+boffset] + β * C[cstart+j*c_cstride+coffset]
        end
      end
    end
  end
  #end # @timeit
end

function contract!(
  C::DenseTensor,
  Clabels,
  A::DenseTensor,
  Alabels,
  B::DiagTensor,
  Blabels,
  α::Number=one(eltype(C)),
  β::Number=zero(eltype(C)),
)
  return contract!(C, Clabels, B, Blabels, A, Alabels, α, β)
end

function show(io::IO, mime::MIME"text/plain", T::DiagTensor)
  summary(io, T)
  return print_tensor(io, T)
end


text:
using BenchmarkTools
using MCLightTransport
using Plots

s = Ray()
# s = PointSource()

μ_a = 1.6  # cm^-1
μ_s = 6.9  # cm^-1
g   = 0.5
nt  = 1.33

# corner1 = [-2.0, -2.0, -2.0]
# corner2 = [2.0, 2.0, 2.0]
# d_xyz = [0.01, 0.01, 0.01]
# cvol = CatesianVolume(corner1, corner2, d_xyz, [μ_a, μ_s, g, nt]);
cvol = CatesianVolume(-2, 2, 0.05, [μ_a, μ_s, g, nt]);

## Running the Simulation
N = 100_000;
@btime simulate!(cvol, s, N);

## Plotting Max Projection
xx, yy, zz = [l:d:h for (l,d,h) in zip(cvol.corner1, cvol.d_xyz, cvol.corner2)]

vol_log = log10.(collect(cvol.data));
maxproj(v, d=1) = dropdims(maximum(v, dims=d), dims=d)

p1 = heatmap(xx, yy, maxproj(vol_log, 3), aspect_ratio=:equal, legend=:none)
p2 = heatmap(xx, zz, maxproj(vol_log, 2), aspect_ratio=:equal, legend=:none)
p3 = heatmap(yy, zz, maxproj(vol_log, 1), aspect_ratio=:equal, legend=:none)
plot(p1,p2,p3, layout=grid(3,1), size=(300,950))