text:
abstract type AbstractBandMatrix{T,B} <: AbstractMatrix{T} end

struct BandMatrix{T,B} <: AbstractBandMatrix{T,B}
    ncol::Base.RefValue{Int}
    data::Vector{T}
end
struct DualBandMatrix{T,B,P} <: AbstractBandMatrix{ForwardDiff.Dual{Nothing,T,P},B}
    bm::BandMatrix{T,B}
    partials::Vector{ForwardDiff.ForwardDiff.Partials{P,T}}
end

function BandMatrix{T,B}(::UndefInitializer, ncols) where {T,B}
    ncol = Ref{Int}(ncols)
    data = Vector{T}(undef, ncols*(B+1))
    # @show ncol, data
    BandMatrix{T,B}(ncol, data)
end
function DualBandMatrix{T,B,P}(::UndefInitializer, ncols) where {T,B,P}
    DualBandMatrix{T,B,P}(
        BandMatrix{T,B}(undef, ncols),
        Vector{ForwardDiff.ForwardDiff.Partials{P,T}}(undef, ncols*(B+1))
    )
end

Random.rand!(bm::BandMatrix) = (rand!(bm.data); bm)
Random.randn!(bm::BandMatrix) = (randn!(bm.data); bm)
Random.randexp!(bm::BandMatrix) = (randexp!(bm.data); bm)

ncol(A::BandMatrix) = A.ncol[]
ncol(A::DualBandMatrix) = A.bm.ncol[]
data(A::BandMatrix) = A.data
data(A::DualBandMatrix) = A.bm.data

Base.size(A::AbstractBandMatrix) = (nc = ncol(A); (nc, nc))
Base.length(A::AbstractBandMatrix) = (ncol(A)^2)

function Base.Array(A::AbstractBandMatrix{T,B}) where {T,B}
    N = ncol(A)
    out = zeros(T, N, N)
    for i ∈ 1:N-B-1
        for b ∈ i:i+B
            out[b,i] = A[b,i]
        end
    end
    for i ∈ N-B:N
        for j ∈ i:N
            out[j,i] = A[j,i]
        end
    end
    out
end

function Base.resize!(bm::BandMatrix{T,B}, N) where {T,B}
    L = (B+1)*N
    bm.ncol[] = N
    resize!(bm.data, L)
    bm
end
function Base.resize!(bm::DualBandMatrix{T,B}, N) where {T,B}
    L = (B+1)*N
    bm.bm.ncol[] = N
    resize!(bm.bm.data, L)
    resize!(bm.partials, L)
    bm
end

function band_sub2ind(::Val{B}, i, j) where {B}
    diff = i - j
    Bp1 = B + 1
    Bp1 - diff + Bp1 * (j - 1 + diff)
    # 1 + diff + Bp1 * (j - 1)
end
function Base.getindex(A::BandMatrix{T,B}, i, j) where {T,B}
    diff = i - j
    if diff < 0
        return zero(T)
    elseif diff > B
        return zero(T)
    end
    A.data[band_sub2ind(Val{B}(), i, j)]
end
function Base.getindex(A::DualBandMatrix{T,B,P}, i, j) where {T,B,P}
    diff = i - j
    if diff < 0
        return zero(T)
    elseif diff > B
        return zero(T)
    end
    ind = band_sub2ind(Val{B}(), i, j)
    ForwardDiff.Dual{Nothing,T,P}(A.bm.data[ind], A.partials[ind])
end
@noinline throwboundserror(a, i) = throw(BoundsError(a, i))
@inline function Base.getindex(A::BandMatrix{T,B}, ind) where {T,B}
    @boundscheck ind > length(A.data) && throwboundserror(A.data, ind)
    @inbounds A.data[ind]
end
@inline function Base.getindex(A::DualBandMatrix{T,B,P}, ind) where {T,B,P}
    @boundscheck ind > length(A.partials) && throwboundserror(A.bm.data, ind)
    @inbounds ForwardDiff.Dual{Nothing,T,P}(A.bm.data[ind], A.partials[ind])
end
@inline function Base.setindex!(A::BandMatrix{T,B}, v, ind) where {T,B}
    @boundscheck ind > length(A.data) && throwboundserror(A.data, ind)
    A.data[ind] = v
end
@inline function Base.setindex!(A::DualBandMatrix{T,B,P}, v::ForwardDiff.Dual{Nothing,T,P}, ind) where {T,B,P}
    @boundscheck ind > length(A.partials) && throwboundserror(A.bm.data, ind)
    @inbounds A.bm.data[ind] = v.value
    @inbounds A.partials[ind] = v.partials
    v
end
@inline function Base.setindex!(A::DualBandMatrix{T,B,P}, v::T, ind) where {T,B,P}
    @boundscheck ind > length(A.partials) && throwboundserror(A.bm.data, ind)
    @inbounds A.bm.data[ind] = v
    # @inbounds A.partials[ind] = v.partials
    v
end
@inline function Base.setindex!(A::DualBandMatrix{T,B,P}, v::ForwardDiff.Partials{P,T}, ind) where {T,B,P}
    @boundscheck ind > length(A.partials) && throwboundserror(A.bm.data, ind)
    @inbounds A.partials[ind] = v.partials
    v
end

@inline function Base.setindex!(A::BandMatrix{T,B}, v, i, j) where {T,B}
    ind = band_sub2ind(Val{B}(), i, j)
    @boundscheck ind > length(A.data) && throwboundserror(A.data, ind)
    A.data[ind] = v
end
@inline function Base.setindex!(A::DualBandMatrix{T,B,P}, v::ForwardDiff.Dual{Nothing,T,P}, i, j) where {T,B,P}
    ind = band_sub2ind(Val{B}(), i, j)
    @boundscheck ind > length(A.partials) && throwboundserror(A.bm.data, ind)
    @inbounds A.bm.data[ind] = v.value
    @inbounds A.partials[ind] = v.partials
    v
end
@inline function Base.setindex!(A::DualBandMatrix{T,B,P}, v::T, i, j) where {T,B,P}
    ind = band_sub2ind(Val{B}(), i, j)
    @boundscheck ind > length(A.partials) && throwboundserror(A.bm.data, ind)
    @inbounds A.bm.data[ind] = v
    # @inbounds A.partials[ind] = v.partials
    v
end
@inline function Base.setindex!(A::DualBandMatrix{T,B,P}, v::ForwardDiff.Partials{P,T}, i, j) where {T,B,P}
    ind = band_sub2ind(Val{B}(), i, j)
    @boundscheck ind > length(A.partials) && throwboundserror(A.bm.data, ind)
    @inbounds A.partials[ind] = v.partials
    v
end

struct DualVector{T,P} <: AbstractVector{ForwardDiff.Dual{Nothing,T,P}}
    data::Vector{T}
    partials::Vector{ForwardDiff.ForwardDiff.Partials{P,T}}
end
function DualVector{T,P}(::UndefInitializer, L) where {T,P}
    DualVector{T,P}(
        Vector{T}(undef, L),
        Vector{ForwardDiff.Partials{P,T}}(undef, L)
    )
end
@inline Base.size(v::DualVector) = size(v.data)
@inline Base.length(v::DualVector) = length(v.data)
function Base.resize!(v::DualVector, N)
    resize!(v.data, N)
    resize!(v.partials, N)
    v
end

@inline function Base.getindex(A::DualVector{T,P}, ind) where {T,P}
    @boundscheck ind > length(A.partials) && throwboundserror(A.bm.data, ind)
    @inbounds ForwardDiff.Dual{Nothing,T,P}(A.data[ind], A.partials[ind])
end
@inline function Base.setindex!(A::DualVector{T,P}, v::ForwardDiff.Dual{Nothing,T,P}, ind) where {T,P}
    @boundscheck ind > length(A.partials) && throwboundserror(A.data, ind)
    @inbounds A.data[ind] = v.value
    @inbounds A.partials[ind] = v.partials
    v
end
@inline function Base.setindex!(A::DualVector{T,P}, v::T, ind) where {T,P}
    @boundscheck ind > length(A.partials) && throwboundserror(A.data, ind)
    @inbounds A.data[ind] = v
    # @inbounds A.partials[ind] = v.partials
    v
end
@inline function Base.setindex!(A::DualVector{T,P}, v::ForwardDiff.Partials{P,T}, ind) where {T,P}
    @boundscheck ind > length(A.partials) && throwboundserror(A.data, ind)
    @inbounds A.partials[ind] = v.partials
    v
end

Base.IndexStyle(::Type{<:DualVector}) = IndexLinear()

function self_dot(v::AbstractVector{T}) where {T}
    out = zero(eltype(v))
    @inbounds @simd for i ∈ eachindex(v)
        out += v[i] * v[i]
    end
    out
end
@generated function self_dot(v::DualVector{T,P}) where {T,P}
    quote
        real = self_dot(v.data)
        Base.Cartesian.@nexprs $P p -> d_p = zero(T)
        # L = length(v.partials)
        # @inbounds @fastmath for i ∈ 0:4:L-4
        #     Base.Cartesian.@nexprs 4 j -> begin
        #         scale_j = 2v.data[i+j]
        #         Base.Cartesian.@nexprs $P p -> d_p += v.partials[i+j][p] * scale_j
        #     end
        # end
        # @inbounds @fastmath for i ∈ L+1-(L&3):L
        #     scale = 2v.data[i]
        #     Base.Cartesian.@nexprs $P p -> d_p += v.partials[i][p] * scale
        # end
        @inbounds @fastmath for i ∈ eachindex(v.partials)
            scale = 2v.data[i]
            Base.Cartesian.@nexprs $P p -> d_p += v.partials[i][p] * scale
        end
        ForwardDiff.Dual(real, ForwardDiff.Partials((Base.Cartesian.@ntuple $P p -> d_p)))
    end
end

Rsym(i, j) = Symbol(:R_, i, :_, j)
Lsym(i, j) = Symbol(:L_, i, :_, j)
Ssym(i, j) = Symbol(:S_, i, :_, j)
sym(S, i, j) = Symbol(S, :_, i, :_, j)

function precision_band_iteration(B, Ldata = :L, ind_expr = :($(B+1)*n), preloop = false, T = Float64)
    Bp1 = B + 1
    q = quote
        $(sym(:R,Bp1,Bp1)) = one($T)
    end
    # Load last row
    for b ∈ 1:B
        push!(q.args, :( $(sym(:L,Bp1, b)) = $Ldata[$ind_expr + $b] ))
    end
    # Calculate that last row from the inverse using preceding rows of inverse
    push!(q.args, :( $(sym(:R, Bp1, B)) = - $(sym(:L,Bp1,B)) * $(sym(:R,B,B)) ))
    for b ∈ 2:B
        push!(q.args, :( $(sym(:R, Bp1, Bp1-b)) = - (
            $(Expr(:call, :+, [:( $(sym(:L,Bp1,Bp1-bᵢ))*$(sym(:R,Bp1-bᵢ,Bp1-b)) ) for bᵢ ∈ 1:b]...))
        ) ))
    end
    # first colum of S equals first column of R
    for b ∈ 2:Bp1
        push!(q.args, :( $(sym(:S,b,1)) = $(sym(:R,b,1)) ))
    end
    for bcol ∈ 2:B # col
        for brow ∈ bcol+1:B+1 # row
            push!(q.args, :( $(sym(:S, brow, bcol)) =
                $(Expr(:call, :+, [:($(sym(:R,brow,bᵢ))*$(sym(:R,bcol,bᵢ))) for bᵢ ∈ 1:bcol]...))
            ))
        end
    end
    # push!(q.args, :( @show $(Expr(:tuple, [ :(abs2($(sym(:R,Bp1,b)))) for b ∈ 1:Bp1 ]...)) ))
    push!(q.args, :( $(sym(:S,Bp1,Bp1)) = $(Expr(:call,:+,1,[ :(abs2($(sym(:R,Bp1,b)))) for b ∈ 1:B ]...)) ) )
    push!(q.args, :( $(sym(:rS,Bp1,Bp1)) = inv(sqrt($(sym(:S,Bp1,Bp1))))))
    for b ∈ 1:B
        push!(q.args, :( $(sym(:S,Bp1,b)) *= $(sym(:rS,Bp1,Bp1))))
    end
    if preloop
        B += 1
        adjust = -1
    else
        adjust = 0
        Rtemp = :_R
        Rtempinv = :_iR
        for b ∈ 2:Bp1
            push!(q.args, :( $(sym(Rtemp,b,1)) = $(sym(:S,b,1)) ))
        end
        push!(q.args, :($(sym(Rtemp,2,2)) = sqrt( one($T) - abs2($(sym(Rtemp,2,1)) ))))
        push!(q.args, :($(sym(Rtempinv,2,2)) = inv($(sym(Rtemp,2,2)))))
        for b ∈ 3:Bp1
            push!(q.args, :($(sym(Rtemp,b,2)) = ( $(sym(:S,b,2)) - $(sym(Rtemp,b,1))*$(sym(Rtemp,2,1)) ) * $(sym(Rtempinv,2,2)) ))
        end
        for b ∈ 3:B
            # push!(q.args, :(@show $(Expr(:tuple, [:($(sym(Rtemp,b,bᵢ))*$(sym(Rtemp,b,bᵢ))) for bᵢ ∈ 1:b-1]...))))
            push!(q.args, :(
                $(sym(Rtemp,b,b)) = sqrt( one($T) - $(Expr(:call, :+, [:($(sym(Rtemp,b,bᵢ))*$(sym(Rtemp,b,bᵢ))) for bᵢ ∈ 1:b-1]...)))
            ))
            push!(q.args, :(
                $(sym(Rtempinv,b,b)) = inv($(sym(Rtemp,b,b)))
            ))
            for brow ∈ b+1:Bp1
                push!(q.args, :($(sym(Rtemp,brow,b)) = ( $(sym(:S,brow,b)) -
                    $(Expr(:call, :+, [:($(sym(Rtemp,brow,bᵢ))*$(sym(Rtemp,b,bᵢ))) for bᵢ ∈ 1:b-1]...)))
                    * $(sym(Rtempinv,b,b))
                ) )
            end
        end
        # for bᵢ ∈ 1:B
        #     push!(q.args, :(@show $(sym(Rtemp,Bp1,bᵢ))))
        # end
        push!(q.args, :(
            $(sym(Rtemp,Bp1,Bp1)) = sqrt( one($T) - $(Expr(:call, :+, [:(abs2($(sym(Rtemp,Bp1,bᵢ)))) for bᵢ ∈ 1:B]...)))
        ))
        push!(q.args, :(
            $(sym(:L,Bp1,Bp1)) = $Ldata[$ind_expr + $Bp1] = last_diag = $(sym(Rtempinv,Bp1,Bp1)) = inv($(sym(Rtemp,Bp1,Bp1)))
        ))
    end

    for b ∈ 2:B
        push!(q.args, :( $(sym(:R,b,1)) = $(sym(:S,b+1+adjust,2+adjust)) ))
    end

    # push!(q.args, :(@show $(sym(:R,2,1))))

    push!(q.args, :( $(sym(:R,2,2)) = sqrt( one($T) - abs2($(sym(:R,2,1))) ) ))
    push!(q.args, :( $(sym(:iR,2,2)) = inv($(sym(:R,2,2))) ))
    for b ∈ 3:B
        push!(q.args, :( $(sym(:R,b,2)) = ($(sym(:S,b+1+adjust,3+adjust)) - $(sym(:R,2,1))*$(sym(:R,b,1))) * $(sym(:iR,2,2)) ))
    end
    for bcol ∈ 3:B
        push!(q.args, :( $(sym(:R,bcol,bcol)) = sqrt(one($T) - $(Expr(:call,:+,[:(abs2($(sym(:R,bcol,bᵢ)))) for bᵢ ∈ 1:bcol-1]...))) ) )
        if bcol < B
            push!(q.args, :( $(sym(:iR,bcol,bcol)) = inv($(sym(:R,bcol,bcol))) ))
        end
        for brow ∈ bcol+1:B
            push!(q.args, :(
                $(sym(:R,brow,bcol)) = ( $(sym(:S,brow+1+adjust,bcol+1+adjust)) - $(Expr(:call,:+,[:($(sym(:R,bcol,bᵢ))*$(sym(:R,brow,bᵢ))) for bᵢ ∈ 1:bcol-1]...)) ) * $(sym(:iR,bcol,bcol))
            ) )
        end
    end

    if preloop
        push!(q.args, :(
            $(sym(:L,Bp1,Bp1)) = $Ldata[$ind_expr + $B] = last_diag = inv($(sym(:R,B,B)))
        ))
    end
    q
end

# Assumes R is
# function fill_precision_factor!(L::BandMatrix{T,B}) where {T,B}
@generated function fill_precision_factor!(L::AbstractBandMatrix{T,B}) where {T,B}
    Bp1 = B + 1
    q = quote
        L[$Bp1] = one($T)
        R_1_1 = 1
        L_2_1 = L[$(Bp1 + B)]
        iR_2_2 = sqrt(1 + L_2_1^2)
        L[$(2Bp1)] = iR_2_2
        R_2_2 = 1 / iR_2_2
        R_2_1 = - L_2_1 * R_2_2
        # @show R_2_1
        # @show L_2_1
        # println("hi")
    end
    for b ∈ 2:B-1
        push!(q.args, precision_band_iteration(b, :L, Bp1*b + Bp1 - b - 1, true, T))
        # push!(q.args,
        #     quote
        #         Smat = SMatrix{3,3}(1.0, S_2_1, S_3_1,
        #                             S_2_1, (R_2_1^2 + R_2_2^2), S_3_2,
        #                             S_3_1, S_3_2, 1.0)
        #         SmatChol = cholesky(Smat).L
        #         Rmat = SMatrix{3,3}(1.0, 0.0, 0.0,
        #                             R_2_1, R_2_2, 0.0,
        #                             R_3_1, R_3_2, R_3_3)'
        #         println("Smat: ")
        #         display(Smat)
        #         println("Smatchol: ")
        #         display(SmatChol)
        #         println("Rmat: ")
        #         display(Rmat)
        #     end)
    end
    push!(q.args,
        quote
            for n ∈ $(B):L.ncol[]-1
                $(precision_band_iteration(B, :L, :($Bp1*n), false, T))
                # Smat = SMatrix{4,4}(1.0, S_2_1, S_3_1, S_4_1,
                #                     S_2_1, (R_2_1^2 + R_2_2^2), S_3_2, S_4_2,
                #                     S_3_1, S_3_2, (R_3_1^2 + R_3_2^2 + R_3_3^2), S_4_3,
                #                     S_4_1, S_4_2, S_4_3, 1.0)
                # SmatChol = cholesky(Smat).L
                # Rmat = SMatrix{4,4}(1.0, 0.0, 0.0, 0.0,
                #                     _R_2_1, _R_2_2, 0.0, 0.0,
                #                     _R_3_1, _R_3_2, _R_3_3, 0.0,
                #                     _R_4_1, _R_4_2, _R_4_3, _R_4_4)'
                # println("Smat: ")
                # display(Smat)
                # println("Smatchol: ")
                # display(SmatChol)
                # println("Rmat: ")
                # display(Rmat)
            end
        end
    )

    quote
        @fastmath @inbounds begin
        # begin
            $q
        end
    end
end

# Assumes R is
# function fill_precision_factor!(L::BandMatrix{T,B}) where {T,B}
@generated function fill_precision_factor_det!(L::AbstractBandMatrix{T,B}) where {T,B}
    Bp1 = B + 1
    q = quote
        L[$Bp1] = one($T)
        R_1_1 = 1
        L_2_1 = L[$(Bp1 + B)]
        iR_2_2 = sqrt(1 + L_2_1^2)
        L[$(2Bp1)] = iR_2_2
        R_2_2 = 1 / iR_2_2
        R_2_1 = - L_2_1 * R_2_2
        log_det = log(iR_2_2)
    end
    for b ∈ 2:B-1
        push!(q.args, precision_band_iteration(b, :L, Bp1*b + Bp1 - b - 1, true, T))
        push!(q.args, :(log_det += log(last_diag)))
    end
    push!(q.args,
        quote
            for n ∈ $(B):L.ncol[]-1
                $(precision_band_iteration(B, :L, :($Bp1*n), false, T))
                log_det += log(last_diag)
            end
        end
    )

    quote
        @fastmath @inbounds begin
        # begin
            $q
        end
        log_det
    end
end

# Assumes R is
# function fill_precision_factor!(L::BandMatrix{T,B}) where {T,B}
@generated function fill_precision_factor!(y, L::AbstractBandMatrix{T,B}, x) where {T,B}
    Bp1 = B + 1
    q = quote
        L[$Bp1] = one($T)
        R_1_1 = 1
        L_2_1 = L[$(Bp1 + B)]
        iR_2_2 = sqrt(1 + L_2_1^2)
        L[$(2Bp1)] = iR_2_2
        R_2_2 = 1 / iR_2_2
        R_2_1 = - L_2_1 * R_2_2
        # det = iR_2_2
        y[1] = (eltype(y))(x[1])
        y[2] = x[2] * iR_2_2 + x[1] * L_2_1
    end
    for b ∈ 2:B-1
        push!(q.args, precision_band_iteration(b, :L, Bp1*b + Bp1 - b - 1, true, T))
        # push!(q.args, :(det += last_diag))
        push!(q.args, :(y[$(b+1)] = $(Expr(:call,:+,[:($(sym(:L,b+1,bᵢ))*x[$(bᵢ)]) for bᵢ ∈ 1:b+1]...))))
    end
    push!(q.args,
        quote
            for n ∈ $(B):L.ncol[]-1
                $(precision_band_iteration(B, :L, :($Bp1*n), false, T))
                # det += last_diag
                push!(q.args, :(y[$(b+1)] = $(Expr(:call,:+,[:($(sym(:L,Bp1,bᵢ))*x[n+$(bᵢ - B)]) for bᵢ ∈ 1:Bp1]...))))
            end
        end
    )

    quote
        @fastmath @inbounds begin
        # begin
            $q
        end
    end
end

# Assumes R is
# function fill_precision_factor!(L::BandMatrix{T,B}) where {T,B}
@generated function fill_precision_factor_det!(y, L::AbstractBandMatrix{T,B}, x) where {B,T}
    Bp1 = B + 1
    q = quote
        L[$Bp1] = one($T)
        R_1_1 = 1
        L_2_1 = L[$(Bp1 + B)]
        iR_2_2 = sqrt(1 + L_2_1^2)
        L[$(2Bp1)] = iR_2_2
        R_2_2 = 1 / iR_2_2
        R_2_1 = - L_2_1 * R_2_2
        log_det = log(iR_2_2)
        y[1] = (eltype(y))(x[1])
        y[2] = x[2] * iR_2_2 + x[1] * L_2_1
    end
    for b ∈ 2:B-1
        push!(q.args, precision_band_iteration(b, :L, Bp1*b + Bp1 - b - 1, true, T))
        push!(q.args, :(log_det += log(last_diag)))
        push!(q.args, :(y[$(b+1)] = $(Expr(:call,:+,[:($(sym(:L,b+1,bᵢ))*x[$(bᵢ)]) for bᵢ ∈ 1:b+1]...))))
    end
    push!(q.args,
        quote
            for n ∈ $(B):ncol(L)-1
                $(precision_band_iteration(B, :L, :($Bp1*n), false, T))
                log_det += log(last_diag)
                # if any(isnan, log_det.partials.values)
                #     @show last_diag
                #     @show log_det
                # end
                y[n+1] = $(Expr(:call,:+,[:($(sym(:L,Bp1,bᵢ))*x[n+$(bᵢ - B)]) for bᵢ ∈ 1:Bp1]...))
            end
        end
    )

    quote
        @fastmath @inbounds begin
        # begin
            $q
        end
        log_det
    end
end

mutable struct InvCholCovar{T,B,P,R,L,O} <: DifferentiableObjects.AbstractDifferentiableObject{P,T}
    bandmat::DualBandMatrix{T,B,P}
    y::DualVector{T,P}
    x::Vector{T}
    δt::Vector{T}
    state::DifferentiableObjects.BFGSState{P,T,R,L}
    initial_x::PaddedMatrices.MutableFixedSizeVector{P,T,R,R}
    ls::DifferentiableObjects.BackTracking2{O,T,Int}
    ∇::PaddedMatrices.MutableFixedSizeVector{P,T,R,R}
end

@generated function InvCholCovar{T,B}(N::Integer) where {T,B}
    P = 2B
    quote
        InvCholCovar(
            DualBandMatrix{$T,$B,$P}(undef, N),
            DualVector{$T,$P}(undef, N),
            Vector{$T}(undef, $B*N),
            Vector{$T}(undef, $B*N),
            DifferentiableObjects.BFGSState(Val($P), $T),
            fill!(PaddedMatrices.MutableFixedSizeVector{$P,$T}(undef), zero($T)),
            DifferentiableObjects.BackTracking2{$T}(Val(2)),
            PaddedMatrices.MutableFixedSizeVector{$P,$T}(undef)
        )
    end
end
@generated function InvCholCovar{T,B}(x) where {T,B}
    P = 2B
    quote
        N = length(x)
        InvCholCovar(
            DualBandMatrix{$T,$B,$P}(undef, N),
            DualVector{$T,$P}(undef, N),
            x,
            Vector{$T}(undef, $B*N),
            DifferentiableObjects.BFGSState(Val($P), $T),
            fill!(PaddedMatrices.MutableFixedSizeVector{P,T,R,R}{$P,$T}(undef), zero($T)),
            DifferentiableObjects.BackTracking2{$T}(Val(2)),
            PaddedMatrices.MutableFixedSizeVector{P,T,R,R}{$P,$T}(undef)
        )
    end
end
function InvCholCovar{T,B}(x, t) where {T,B}
    icc = InvCholCovar{T,B}(x)
    fill_δt!(icc.δt, t, Val(B))
    icc
end

function Base.resize!(icc::InvCholCovar{T,B,P}, N) where {T,B,P}
    resize!(icc.bandmat, N)
    resize!(icc.y, N)
    resize!(icc.x, N)
    resize!(icc.δt, N*B)
    icc
end

function fit!(icc::InvCholCovar{T}) where {T}
    DifferentiableObjects.optimize_scale!(icc.state, icc, icc.initial_x, icc.ls, T(10), T(1e-3))
end

@generated function fill_δt!(δt, t, ::Val{B}) where B
    Bp1 = B + 1
    quote
        Base.Cartesian.@nexprs $B i -> begin
            tᵢ = t[i+1]
            Base.Cartesian.@nexprs i b -> (δt[i*$B + $B - i + b] = t[b] - tᵢ)
        end
        for i ∈ 1:length(t) - $Bp1
            tᵢ = t[i+$Bp1]
            Base.Cartesian.@nexprs $B b -> (δt[i*$B + $(B*B) + b] = t[i+b] - tᵢ)
        end
    end
end

function decorrelate_data!(icc::InvCholCovar{T,B,P}, X::AbstractMatrix, t) where {T,B,P}
    N, M = size(X)
    resize!(icc, N)
    fill_δt!(icc.δt, t, Val(B))
    for m ∈ 1:M
        allfinite = false
        for i ∈ 1:100
        # icc.x .= @view X[:,m]
            stdx = std(@view(X[:,m]))
            invstdx = 1 / stdx
            @inbounds @simd ivdep for i in 1:size(X,1)
                icc.x[i] = X[i,m] * invstdx
            end
            DifferentiableObjects.optimize_scale!(icc.state, icc, icc.initial_x, icc.ls, T(10), 1e-3)
            allfinite = all(isfinite, icc.y.data)
            if !allfinite
                Random.randn!(icc.initial_x)
                continue
            end
            @inbounds @simd ivdep for i in 1:size(X,1)
                X[i,m] = icc.y.data[i] * stdx
            end
            icc.initial_x .= zero(T)
            break
        end
        allfinite || throw("Not all finite after 100 attempts.")
    end
end
function decorrelate_data!(icc::InvCholCovar{T,B,P}, X::AbstractVector, t) where {T,B,P}
    N = length(X)
    resize!(icc, N)
    fill_δt!(icc.δt, t, Val(B))
    stdx = std(X)
    invstdx = 1 / stdx
    if isa(X, Vector{T})
        icc.x = X
        icc.x .*= invstdx
    else
        icc.x .= X .* invstdx
    end
    DifferentiableObjects.optimize_scale!(icc.state, icc, icc.initial_x, icc.ls, T(10), 1e-3)
    X .= icc.y.data .* stdx
end

@generated zero_tuple(::Val{N}, ::Type{T}) where {N,T} = Expr(:tuple, [zero(T) for n ∈ 1:N]...)
@inline function inv_sigmoid(x)
    expm1x = expm1(x)
    expm1x/(2+expm1x)
end
@inline function ∂inv_sigmoid(x)
    expx = exp(x)
    isx = (expx-1)/(expx+1)
    ∂isx = 2expx/(expx+1)^2
    isx, ∂isx
end

@generated function DifferentiableObjects.fdf(icc::InvCholCovar{T,B,P}, θ::AbstractVector{T}) where {T,B,P}
    Bp1 = B+1
    quote
        # Σθ² = zero(T)
        Base.Cartesian.@nexprs $B b -> begin
            @inbounds θ_b = θ[b]
            @inbounds θ_{$B+b} = exp(θ[$B+b])
            # signθ_b = sign(θ_b)
            # absθ_b, ∂absθ_b = ∂inv_sigmoid(abs(θ_b))
            # Σθ² = fma(θ_b, θ_b, Σθ²)
            # icc.∇[b] = θ_b + θ_b
        end
        δta = icc.δt
        @inbounds  begin
            Base.Cartesian.@nexprs $B i -> begin
                # tᵢ = t[i]
                Base.Cartesian.@nexprs i b -> begin
                    δt = δta[i*$B + $B - i + b]
                    # v1 = absθ_{$B - i + b}^(δt-1)
                    # v = copysign(absθ_{$B - i + b} * v1, signθ_{$B - i + b})
                    # p = δt * v1 * ∂absθ_{$B - i + b}
                    # icc.bandmat[i*$Bp1 + $B - i + b] = ForwardDiff.Dual(v,ForwardDiff.Partials(Base.setindex(zero_tuple(Val(P), T), p, $B - i + b)))

                    ∂ϕ_b = exp(δt * θ_{$P - i + b})
                    ϕ_b = θ_{$B - i + b} * ∂ϕ_b
                    ∂ϕ_{$B+b} = ϕ_b * δt * θ_{$P - i + b}

                    ∂ϕ_btup = Base.setindex(Base.setindex(zero_tuple(Val(P), T), ∂ϕ_b, $B - i + b), ∂ϕ_{$B+b}, $P - i + b)
                    icc.bandmat[i*$Bp1 + $B - i + b] = ForwardDiff.Dual(ϕ_b, ForwardDiff.Partials(∂ϕ_btup))
                end
            end
        end
        @inbounds for i ∈ 1:length(icc.x) - $Bp1
            # tᵢ = t[i+$Bp1]
            Base.Cartesian.@nexprs $B b -> begin
                δt = δta[i*$B + $(B*B) + b]
                # v1 = absθ_b^(δt-1)
                # v = copysign(absθ_b * v1, signθ_b)
                # p = δt * v1 * ∂absθ_b
                # icc.bandmat[i*$Bp1 + $(Bp1*B) + b] = ForwardDiff.Dual(v,ForwardDiff.Partials(Base.setindex(zero_tuple(Val(P), T), p, b)))
                # ϕ_b, ∂ϕ_b = ∂inv_sigmoid( θ_b*δt )
                # icc.bandmat[i*$Bp1 + $(Bp1*B) + b] = ForwardDiff.Dual(ϕ_b,ForwardDiff.Partials(Base.setindex(zero_tuple(Val(P), T), ∂ϕ_b * δt, b)))
                ∂ϕ_b = exp(δt * θ_{$B + b})
                ϕ_b = θ_b * ∂ϕ_b
                ∂ϕ_{$B+b} = ϕ_b * δt * θ_{$B + b}
                ∂ϕ_btup = Base.setindex(Base.setindex(zero_tuple(Val(P), T), ∂ϕ_b, b), ∂ϕ_{$B+b}, $B + b)
                icc.bandmat[i*$Bp1 + $(Bp1*B) + b] = ForwardDiff.Dual(ϕ_b,ForwardDiff.Partials(∂ϕ_btup))
            end
        end
        logdet = fill_precision_factor_det!(icc.y, icc.bandmat, icc.x)
        out = self_dot(icc.y) - logdet
        Base.Cartesian.@nexprs $P p -> @inbounds icc.∇[p] = out.partials[p]
        out.value #+ Σθ²
    end
end
@generated function DifferentiableObjects.f(icc::InvCholCovar{T,B,P}, θ::AbstractVector{T}) where {T,B,P}
    Bp1 = B+1
    quote
        # Σθ² = zero(T)
        Base.Cartesian.@nexprs $B b -> begin
            @inbounds θ_b = θ[b]
            @inbounds θ_{$B+b} = exp(θ[$B+b])
            # signθ_b = sign(θ_b)
            # absθ_b = abs(inv_sigmoid(θ_b))
            # Σθ² = fma(θ_b, θ_b, Σθ²)
        end
        # @show length(icc.bandmat.bm.data)
        # t = icc.t
        δta = icc.δt
        @inbounds begin
            Base.Cartesian.@nexprs $B i -> begin
                # tᵢ = t[i]
                Base.Cartesian.@nexprs i b -> begin
                    δt = δta[i*$B + $B - i + b]
                    # v = copysign(absθ_{$B - i + b}^δt, signθ_{$B - i + b})
                    # icc.bandmat.bm[i*$Bp1 + $B - i + b] = v
                    v = θ_{$B - i + b} * exp(δt * θ_{$P - i + b})
                    icc.bandmat.bm[i*$Bp1 + $B - i + b] = v# inv_sigmoid(θ_{$B - i + b} * δt)
                end
            end
        end
        @inbounds for i ∈ 1:length(icc.x) - $Bp1
            # tᵢ = t[i+$Bp1]
            Base.Cartesian.@nexprs $B b -> begin
                δt = δta[i*$B + $(B*B) + b]
                # v = copysign(absθ_b^δt, signθ_b)
                # icc.bandmat.bm[i*$Bp1 + $(Bp1*B) + b] = v
                v = θ_b * exp(δt * θ_{$B + b})
                icc.bandmat.bm[i*$Bp1 + $(Bp1*B) + b] = v #inv_sigmoid(θ_b * δt)
            end
        end
        logdet = fill_precision_factor_det!(icc.y.data, icc.bandmat.bm, icc.x)
        self_dot(icc.y.data) - logdet
    end
end

function DifferentiableObjects.scale_fdf(icc::InvCholCovar{T}, θ::AbstractVector{T}, scale_target) where {T}
    fval = DifferentiableObjects.fdf(icc, θ)
    scale = min(one(T), scale_target / norm(icc.∇))
    @inbounds @simd for i ∈ eachindex(icc.∇)
        icc.∇[i] *= scale
    end
    fval * scale, scale
end
function DifferentiableObjects.scaled_fdf(icc::InvCholCovar{T}, θ::AbstractVector{T}, scale) where {T}
    fval = DifferentiableObjects.fdf(icc, θ)
    @inbounds @simd for i ∈ eachindex(icc.∇)
        icc.∇[i] *= scale
    end
    fval * scale
end
DifferentiableObjects.gradient(icc::InvCholCovar) = icc.∇

@generated function logdettest(icc::InvCholCovar{T,B,P}, θ::AbstractVector{T}) where {T,B,P}
    Bp1 = B+1
    quote
        Base.Cartesian.@nexprs $B b -> begin
            @inbounds θ_b = θ[b]
            signθ_b = sign(θ_b)
            absθ_b, ∂absθ_b = ∂inv_sigmoid(abs(θ_b))
        end
        δta = icc.δt
        @inbounds  begin
            Base.Cartesian.@nexprs $B i -> begin
                # tᵢ = t[i]
                Base.Cartesian.@nexprs i b -> begin
                    # δt = tᵢ - t[b]
                    δt = δta[i*$B + $B - i + b]
                    v1 = absθ_{$B - i + b}^(δt-1)
                    v = copysign(absθ_{$B - i + b} * v1, signθ_{$B - i + b})
                    p = δt * v1 * ∂absθ_{$B - i + b}
                    # v = - exp(θ_b * δt)
                    # p = v * δt
                    # dual =
                    # @show dual
                    # ForwardDiff.Dual(v,ForwardDiff.Partials(Base.setindex(zero_tuple(Val(P), T), p, $B - i + b)))
                    v =
                    ∂ =
                    dualv = ForwardDiff.Dual(v,ForwardDiff.Partials(Base.setindex(zero_tuple(Val(P), T), p, $B - i + b)))
                    icc.bandmat[i*$Bp1 + $B - i + b] = dualv
                end
            end
        end
        @inbounds for i ∈ 1:length(icc.x) - $Bp1
            # tᵢ = t[i+$Bp1]
            Base.Cartesian.@nexprs $B b -> begin
                # δt = tᵢ - t[i+b]
                δt = δta[i*$B + $(B*B) + b]
                v1 = absθ_b^(δt-1)
                v = copysign(absθ_b * v1, signθ_b)
                p = δt * v1 * ∂absθ_b
                # v = - exp(θ_b * δt)
                # p = v * δt
                icc.bandmat[i*$Bp1 + $(Bp1*B) + b] = ForwardDiff.Dual(v,ForwardDiff.Partials(Base.setindex(zero_tuple(Val(P), T), p, b)))
            end
        end
        out = - fill_precision_factor_det!(icc.y, icc.bandmat, icc.x)
        Base.Cartesian.@nexprs $P p -> @inbounds icc.∇[p] = out.partials[p]
        out.value
    end
end

@generated function selfdottest(icc::InvCholCovar{T,B,P}, θ::AbstractVector{T}) where {T,B,P}
    Bp1 = B+1
    quote
        Base.Cartesian.@nexprs $B b -> begin
            @inbounds θ_b = θ[b]
            signθ_b = sign(θ_b)
            absθ_b, ∂absθ_b = ∂inv_sigmoid(abs(θ_b))
        end
        δta = icc.δt
        @inbounds  begin
            Base.Cartesian.@nexprs $B i -> begin
                # tᵢ = t[i]
                Base.Cartesian.@nexprs i b -> begin
                    # δt = tᵢ - t[b]
                    δt = δta[i*$B + $B - i + b]
                    v1 = absθ_{$B - i + b}^(δt-1)
                    v = copysign(absθ_{$B - i + b} * v1, signθ_{$B - i + b})
                    p = δt * v1 * ∂absθ_{$B - i + b}
                    # v = - exp(θ_b * δt)
                    # p = v * δt
                    # dual =
                    # @show dual
                    icc.bandmat[i*$Bp1 + $B - i + b] = ForwardDiff.Dual(v,ForwardDiff.Partials(Base.setindex(zero_tuple(Val(P), T), p, $B - i + b)))
                end
            end
        end
        @inbounds for i ∈ 1:length(icc.x) - $Bp1
            # tᵢ = t[i+$Bp1]
            Base.Cartesian.@nexprs $B b -> begin
                # δt = tᵢ - t[i+b]
                δt = δta[i*$B + $(B*B) + b]
                v1 = absθ_b^(δt-1)
                v = copysign(absθ_b * v1, signθ_b)
                p = δt * v1 * ∂absθ_b
                # v = - exp(θ_b * δt)
                # p = v * δt
                icc.bandmat[i*$Bp1 + $(Bp1*B) + b] = ForwardDiff.Dual(v,ForwardDiff.Partials(Base.setindex(zero_tuple(Val(P), T), p, b)))
            end
        end
        logdet = fill_precision_factor_det!(icc.y, icc.bandmat, icc.x)
        out = self_dot(icc.y)
        Base.Cartesian.@nexprs $P p -> @inbounds icc.∇[p] = out.partials[p]
        out.value
    end
end


text:
using Printf
try
  GaussQuadrature.neither
catch
  include("quadrature.jl")
end

# From the VTK source:
# https://gitlab.kitware.com/nick.laurenson/vtk/blob/d8aa5b89c622cf04b3322112fda420afc9d2a16d/Common/DataModel/vtkLagrangeHexahedron.cxx#L729-795
function PointIndexFromIJK(i, j, k, N)
  ibdy = (i == 0 || i == N)
  jbdy = (j == 0 || j == N)
  kbdy = (k == 0 || k == N)
  # How many boundaries do we lie on at once?
  nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0) + (kbdy ? 1 : 0)

  if (nbdy == 3) # Vertex DOF
    return (i != 0 ? (j != 0 ? 2 : 1) : (j != 0 ? 3 : 0)) + (k != 0 ? 4 : 0)
  end

  offset = 8
  if (nbdy == 2) # Edge DOF
    if (!ibdy)
      # On i axis
      return (i - 1) +
      (j != 0 ? N - 1 + N - 1 : 0) +
      (k != 0 ? 2 * (N - 1 + N - 1) : 0) +
      offset
    end
    if (!jbdy)
      # On j axis
      return (j - 1) +
      (i != 0 ? N - 1 : 2 * (N - 1) + N - 1) +
      (k != 0 ? 2 * (N - 1 + N - 1) : 0) +
      offset
    end
    # !kbdy, On k axis
    offset += 4 * (N - 1) + 4 * (N - 1)
    return (k - 1) + (N - 1) * (i != 0 ? (j != 0 ? 3 : 1) : (j != 0 ? 2 : 0)) + offset
  end

  offset += 4 * (N - 1 + N - 1 + N - 1)
  if (nbdy == 1) # Face DOF
    if (ibdy) # On i-normal face
      return (j - 1) + ((N - 1) * (k - 1)) + (i != 0 ? (N - 1) * (N - 1) : 0) + offset
    end
    offset += 2 * (N - 1) * (N - 1)
    if (jbdy) # On j-normal face
      return (i - 1) + ((N - 1) * (k - 1)) + (j != 0 ? (N - 1) * (N - 1) : 0) + offset
    end
    offset += 2 * (N - 1) * (N - 1)
    # kbdy, On k-normal face
    return (i - 1) + ((N - 1) * (j - 1)) + (k != 0 ? (N - 1) * (N - 1) : 0) + offset
  end

  # nbdy == 0: Body DOF
  offset += 2 * ( (N - 1) * (N - 1) +
                 (N - 1) * (N - 1) +
                 (N - 1) * (N - 1))
  return offset + (i - 1) + (N - 1) * ((j - 1) + (N - 1) * ((k - 1)))
end

let
  N = 5
  nelem = 4
  T = Float64
  dim = 3

  io = open("data$(dim)d.vtu", "w")
  Nq = N+1
  (r, _) = GaussQuadrature.legendre(T, Nq, GaussQuadrature.both)
  s = r
  t = dim == 2 ? [0] : r
  Nqk = length(t)
  Np = Nq * Nq * Nqk

  x = Array{T, 4}(undef, Nq, Nq, Nqk, nelem)
  y = Array{T, 4}(undef, Nq, Nq, Nqk, nelem)
  z = Array{T, 4}(undef, Nq, Nq, Nqk, nelem)

  for e = 1:nelem, k = 1:Nqk, j = 1:Nq, i = 1:Nq
    xoffset = 2e - 1 - nelem
    x[i, j, k, e], y[i, j, k, e], z[i, j, k, e] = r[i]-xoffset, s[j], t[k]
  end
  if dim == 2
    x, y = x + sin.(π * y) / 5, y + exp.(-x.^2)
  else
    x, y, z = x + sin.(π * y) / 5, y + exp.(-hypot.(x, z).^2), z + sin.(π * x) / 5
  end
  d = exp.(sin.(hypot.(x, y, x)))

  write(io,
        """
        <?xml version="1.0"?>
        <VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
          <UnstructuredGrid>
            <Piece NumberOfPoints="$(Np * nelem)" NumberOfCells="$nelem">
              <PointData>
                <DataArray type="Float32" Name="element_number" format="ascii">
        """
       )
  for n = 1:length(d)
    @printf(io, "           %16e\n", div(n-1, Np))
  end
  write(io,
        """
                </DataArray>
                <DataArray type="Float32" Name="data" format="ascii">
        """
       )
  for n = 1:length(d)
    @printf(io, "           %16e\n", d[n])
  end
  write(io,
        """
                </DataArray>
              </PointData>
              <CellData>
              </CellData>
              <Points>
                <DataArray type="Float32" Name="Points" NumberOfComponents="3" format="ascii">
        """
       )
  # Fill coordinates
  for n = 1:length(x)
    @printf(io, "           %16e %16e %16e\n", x[n], y[n], z[n])
  end

  write(io,
        """
                </DataArray>
              </Points>
              <Cells>
                <DataArray type="Int64" Name="connectivity" format="ascii">
        """
       )
  # Fill coordinates
  L = LinearIndices((1:Nq, 1:Nq, 1:Nqk))
  for e = 1:nelem
    o = (e-1) * Np - 1
    if dim == 2
      # corners
      @printf(io, "           %d %d %d %d\n",
              o + L[1, 1, 1],
              o + L[end, 1, 1],
              o + L[end, end, 1],
              o + L[1, end, 1])
      # edges
      @printf(io, "          ")
      for n = 2:Nq-1
        @printf(io, " %d", o + L[n, 1, 1])
      end
      @printf(io, "\n")
      @printf(io, "          ")
      for n = 2:Nq-1
        @printf(io, " %d", o + L[end, n, 1])
      end
      @printf(io, "\n")
      @printf(io, "          ")
      for n = 2:Nq-1
        @printf(io, " %d", o + L[n, end, 1])
      end
      @printf(io, "\n")
      @printf(io, "          ")
      for n = 2:Nq-1
        @printf(io, " %d", o + L[1, n, 1])
      end
      @printf(io, "\n")
      # interior
      for j = 2:Nq-1
        @printf(io, "          ")
        for i = 2:Nq-1
          @printf(io, " %d", o + L[i, j, 1])
        end
        j != Nq-1 && @printf(io, "\n")
      end
      @printf(io, "\n")
    elseif dim == 3
      # corners
      @printf(io, "           %d %d %d %d %d %d %d %d\n",
              o + L[1, 1, 1],
              o + L[end, 1, 1],
              o + L[end, end, 1],
              o + L[1, end, 1],
              o + L[1, 1, end],
              o + L[end, 1, end],
              o + L[end, end, end],
              o + L[1, end, end])
      tmp = 0
      @assert tmp == PointIndexFromIJK(0, 0, 0, N); tmp += 1
      @assert tmp == PointIndexFromIJK(N, 0, 0, N); tmp += 1
      @assert tmp == PointIndexFromIJK(N, N, 0, N); tmp += 1
      @assert tmp == PointIndexFromIJK(0, N, 0, N); tmp += 1
      @assert tmp == PointIndexFromIJK(0, 0, N, N); tmp += 1
      @assert tmp == PointIndexFromIJK(N, 0, N, N); tmp += 1
      @assert tmp == PointIndexFromIJK(N, N, N, N); tmp += 1
      @assert tmp == PointIndexFromIJK(0, N, N, N); tmp += 1
      # edges
      @printf(io, "          ")
      j = 1
      k = 1
      for i = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      @printf(io, "          ")
      i = Nq
      k = 1
      for j = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      @printf(io, "          ")
      j = Nq
      k = 1
      for i = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      @printf(io, "          ")
      i = 1
      k = 1
      for j = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")

      @printf(io, "          ")
      j = 1
      k = Nq
      for i = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      @printf(io, "          ")
      i = Nq
      k = Nq
      for j = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      @printf(io, "          ")
      j = Nq
      k = Nq
      for i = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      @printf(io, "          ")
      i = 1
      k = Nq
      for j = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")

      @printf(io, "          ")
      i = 1
      j = 1
      for k = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      @printf(io, "          ")
      i = Nq
      j = 1
      for k = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      @printf(io, "          ")
      i = 1
      j = Nq
      for k = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      @printf(io, "          ")
      i = Nq
      j = Nq
      for k = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      # faces
      @printf(io, "          ")
      i = 1
      for k = 2:Nq-1, j = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      @printf(io, "          ")
      i = Nq
      for k = 2:Nq-1, j = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      @printf(io, "          ")
      j = 1
      for k = 2:Nq-1, i = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      @printf(io, "          ")
      j = Nq
      for k = 2:Nq-1, i = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      @printf(io, "          ")
      k = 1
      for j = 2:Nq-1, i = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      @printf(io, "          ")
      k = Nq
      for j = 2:Nq-1, i = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
      # interior
      @printf(io, "          ")
      for k = 2:Nq-1, j = 2:Nq-1, i = 2:Nq-1
        @printf(io, " %d", o + L[i, j, k])
        @assert tmp == PointIndexFromIJK(i-1, j-1, k-1, N); tmp += 1
      end
      @printf(io, "\n")
    else
      error("invalid dim")
    end
  end
  write(io,
        """
                </DataArray>
                <DataArray type="Int64" Name="offsets" format="ascii">
        """
       )
  for e = 1:nelem
    @printf(io, "           %d\n", Np * e)
  end
  write(io,
        """
                </DataArray>
                <DataArray type="UInt8" Name="types" format="ascii">
        """
       )
  if dim == 2
    for e = 1:nelem
      @printf(io, "           70\n")
    end
  elseif dim == 3
    for e = 1:nelem
      @printf(io, "           72\n")
    end
  else
    error("invalid dim")
  end
  write(io,
        """
                </DataArray>
              </Cells>
            </Piece>
          </UnstructuredGrid>
        </VTKFile>
        """
       )

  close(io)
end

ext: jl


text:
"""
    findall_nan_chunks(x)

Finds the `(start, stop)` indices of each `NaN` chunk in `x` and returns a vector
of index tuples for the `NaN` ranges.

See also: [`findall_nan_chunks!`](@ref)

## Examples

```julia
x = [NaN, NaN, 2.3, NaN, 5.6, NaN, NaN, NaN]
findall_nan_chunks(x)
```
"""
function findall_nan_chunks(x)
    N = length(x)
    x_isnans = isnan.(x)

    # (start_idx, stop_idx) for each chunk of nans
    nan_chunks = Vector{Tuple{Int, Int}}(undef, 0)

    i = 1
    @inbounds while i <= N
        if x_isnans[i] == true
            # Now we know that the Nan range starts at position i
            j = i

            # Update j until there are no more nans
            while j < N && x_isnans[j + 1]
                j += 1
            end

            # So now we know the end of the range is at position j
            push!(nan_chunks, (i, j))

            # Update current position based on how many NaNs we found
            i = j + 1
        end

        i += 1
    end

    return nan_chunks
end

"""
    findall_nan_chunks!(v, x)

Finds the (start,a stop) indices of each `NaN` chunk in `x` and returns a vector
of those index tuples, using a preallocated boolean vector `v`, where
`length(x) == length(v)`, to keep track of `NaN` positions.

See also: [`findall_nan_chunks`](@ref)

## Example

```julia
x = [NaN, NaN, 2.3, NaN, 5.6, NaN, NaN, NaN]
v = zeros(Bool, length(x))
findall_nan_chunks!(v, x)
```
"""
function findall_nan_chunks!(v::AbstractVector{Bool}, x)
    v .= isnan.(x)

    # (start_idx, stop_idx) for each chunk of nans
    nan_chunks = Vector{Tuple{Int, Int}}(undef, 0)

    i = 1
    N = length(x)

    @inbounds while i <= N
        if v[i] == true
            # Now we know that the Nan range starts at position i
            j = i

            # Update j until there are no more nans
            while j < N && v[j + 1]
                j += 1
            end

            # So now we know the end of the range is at position j
            push!(nan_chunks, (i, j))

            # Update current position based on how many NaNs we found
            i = j + 1
        end

        i += 1
    end

    return nan_chunks
end

export findall_nan_chunks, findall_nan_chunks!


text:
function SetGen(generators, nodes, bus_gep)
   setN = 1:length(nodes)
   setI = 1:length(generators)

   wind_LBset = Int64[]
   wind_OSset = Int64[]
   pvset = Int64[]
   pv_BTMset = Int64[]
   hydroset = Int64[]
   nucset = Int64[]
   coalset = Int64[]
   oilset = Int64[]
   gasset = Int64[]

   for g in setI
      if generators[g].gtype == "Wind_LB"
         push!(wind_LBset, generators[g].gindex)
      elseif generators[g].gtype == "Wind_OS"
         push!(wind_OSset, generators[g].gindex)
      elseif generators[g].gtype == "PV"
         push!(pvset, generators[g].gindex)
      elseif generators[g].gtype == "PV_BTM"
         push!(pv_BTMset, generators[g].gindex)
      elseif generators[g].gtype == "Hydro"
         push!(hydroset, generators[g].gindex)
      elseif generators[g].gtype == "Nuclear"
         push!(nucset, generators[g].gindex)
      elseif generators[g].gtype == "Coal"
         push!(coalset, generators[g].gindex)
      elseif generators[g].gtype == "Oil"
         push!(oilset, generators[g].gindex)
      elseif generators[g].gtype == "Gas"
         push!(gasset, generators[g].gindex)
      end
   end

   # Avoid MIP (Not implementing UC problem)
   for g in setI
      generators[g].Pmin = 0
   end

   setI_hat = setI[end]+1 : setI[end]+5
   push!(wind_LBset, setI[end]+1)
   push!(wind_OSset, setI[end]+2)
   push!(pvset, setI[end]+3)
   push!(pv_BTMset, setI[end]+4)
   push!(gasset, setI[end]+5)

   B_g = []
   for g in setI
      push!(B_g, generators[g].location)
   end

   for jj in 1:5
      push!(B_g, bus_gep)
   end

   B_gn = Array{Array{Int64},1}(undef,length(nodes))
   for ii in 1:length(nodes)
      B_gn[ii] = Int64[]
   end
   for g in union(setI,setI_hat)
      push!(B_gn[B_g[g]], g)
   end

   setI = union(setI, setI_hat)

   I_R = intersect(setI, union(wind_LBset,wind_OSset,pvset,pv_BTMset,hydroset))
   I_R_hat = intersect(setI_hat, union(wind_LBset,wind_OSset,pvset,pv_BTMset,hydroset))
   I_C = setdiff(setI, union(wind_LBset,wind_OSset,pvset,pv_BTMset,hydroset))
   I_C_hat = setdiff(setI_hat,union(wind_LBset,wind_OSset,pvset,pv_BTMset,hydroset))

   return wind_LBset, wind_OSset, pvset, pv_BTMset, hydroset, nucset, coalset, oilset, gasset, B_g, B_gn, setI, setI_hat, I_R, I_R_hat, I_C, I_C_hat
end

max_stars_repo_name: JipKim/IA

max_issues_repo_name: JipKim/IA

max_forks_repo_name: JipKim/IA


text:
using Distances

using mjWrap
using MuJoCo

################################################## Experiment Centric Functions
## accepts model and trajectory data
## returns sum of trajectory rewards for each traj; saves per timestep reward
function rewardfunction(x::mjWrap.mjSet, tid::Int,
                        s0::AbstractVector{Float64}, # pass in states
                        s1::AbstractVector{Float64},
                        o0::AbstractVector{Float64}, # observations
                        o1::AbstractVector{Float64},
                        ctrl::AbstractVector{Float64})::Float64

   d = x.datas[tid]
   #hand = d.site_xpos[:,1]
   #goal = d.site_xpos[:,2]
   #ball = d.site_xpos[:,3]

   d_hand = o1[27] #norm(hand-ball)^2
   d_goal = o1[28] #norm(goal-ball)^2

   #reward = 1.0 - 3.0*d_goal - d_hand - 1e-1*norm(ctrl)^2
   # kendall
   #if d_hand >= 0.02
   #   reward = 1.0 - 1.0*d_hand
   #else
   #   reward = 2.0 - 1.0*d_goal
   #   if d_goal < 0.02
   #      reward += 3.0
   #   end
   #end
   ##reward -= 1e-1*norm(ctrl)^2
   #reward -= 1e-2*norm(d.qvel)^2

   #reward = 5 - 1*d_hand - 10*d_goal
   reward = -1*d_hand - 1*d_goal
   #if d_hand < 0.005
   #   reward += 5
   #end
   #if d_goal < 0.005
   #   reward += 10
   #end
   #reward -= 1e-2*norm(d.qvel)^2
   reward -= 1e-2*norm(ctrl)^2

   #mark = mj.MARKSTACK(x.datas[tid])
   #My = @view d.stack[mark+1:(mark+x.nv)]
   #mj.solveM(x.m.m, x.datas[tid].d, My, x.datas[tid].qfrc_actuator, x.nv)
   #reward -= 1e-4*(x.datas[tid].qfrc_actuator'*My)
   #mj.FREESTACK(x.datas[tid], Int32(mark))

   return reward #max(reward, -3.0)
end

function obsfunction(x::mjWrap.mjSet, tid::Int,
                     s::AbstractVector{Float64},
                     o::AbstractVector{Float64})
   hand = x.datas[tid].site_xpos[:,1]
   goal = x.datas[tid].site_xpos[:,2]
   ball = x.datas[tid].site_xpos[:,3]

   d_hand = norm(hand-ball)^2
   d_goal = norm(goal-ball)^2

   # magnet for ball to hand
   o[31:33] .= 0.0 #cosine_dist(hand, goal)
   if d_hand < 0.005
      x.datas[tid].qpos[11:13] .= hand # set ball qpos to be hand's
      s[11:13] .= hand # ALSO MODIFY THE STATE
      o[31:33] .= 1.0
   end

   o[1:13] = x.datas[tid].qpos
   o[14:26] .= x.datas[tid].qvel./10.0
   o[27] = d_hand
   o[28] = d_goal
   o[29] = 0.0 #cosine_dist(hand, ball)
   o[30] = 0.0 #cosine_dist(hand, goal)
end

## accepts model,
## default state,
## preallocated matrix of states (nqnv x numT)
function initfunction!(x::mjWrap.mjSet,
                       s::mjWrap.TrajSamples)
   s0 = s.s0
   s0[1:x.nq] = [0.369814, 0.40224, -1.18184, -1.50206, 1.46349, -0.215438, -0.000427812, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
   mag = 0.2
   z = 0.6
   rot = pi #2pi #2*pi
   init_state = view( s.state, :, 1, : )
   init_state .= 0.0
   for t=1:s.numT
      init_state[1:x.nq,t] = s0[1:x.nq] + rand(x.nq) * mag - mag/2.0

      # goal
      init_state[ 8:10,t] = z*rand(3) - z/2
      init_state[   10,t] += 0.3

      # ball
      init_state[11:13,t] = z*rand(3) - z/2
      init_state[   13,t] += 0.3
   end
end

function ctrlfunction!(x::mjWrap.mjSet, ctrl::Array{Float64})
   randn!(ctrl)
end

# modelfunc should takes in time index (iteration index) to know change curric
function modelfunction!(x::mjWrap.mjSet, iter::Int=0)
end

## evaluation function: reward function shapes exploration, evaluation
## function defines success or failure
function evalfunction(x::mjWrap.mjSet, d::mjWrap.TrajSamples)
   return mean(d.obs[28,end,:]) # x position at time T, for all numT
end

using ExpFunctions
myfuncs = FunctionSet(modelfunction!,
                      initfunction!,
                      ctrlfunction!,
                      obsfunction,
                      rewardfunction,
                      evalfunction)


text:
(function (main)
  (for ((line (eachline "data/words.txt")))
    (if (> (length line) 20)
        (println line))))

(main)


text:
@testset "maximization wrapper" begin
    @testset "univariate" begin
        resmax = maximize(x->x^3, -1, 9)
        resmin = optimize(x->-x^3, -1, 9)
        @test Optim.maximum(resmax) == -Optim.minimum(resmin)
        @test resmax.res.minimum == resmin.minimum
        for meth in (Brent(), GoldenSection())
            resmax = maximize(x->x^3, -1, 9, meth)
            resmin = optimize(x->-x^3, -1, 9, meth)
            @test Optim.maximum(resmax) == -Optim.minimum(resmin)
            @test resmax.res.minimum == resmin.minimum
        end
    end
    @testset "multivariate" begin
        resmax = maximize(x->x[1]^3+x[2]^2, [3.0, 0.0])
        resmin = optimize(x->-x[1]^3-x[2]^2, [3.0, 0.0])
        @test Optim.maximum(resmax) == -Optim.minimum(resmin)
        @test resmax.res.minimum == resmin.minimum
        for meth in (NelderMead(), BFGS(), LBFGS(), GradientDescent(), Newton(), NewtonTrustRegion(), SimulatedAnnealing())
            resmax = maximize(x->x[1]^3+x[2]^2, [3.0, 0.0])
            resmin = optimize(x->-x[1]^3-x[2]^2, [3.0, 0.0])
            @test Optim.maximum(resmax) == -Optim.minimum(resmin)
            @test resmax.res.minimum == resmin.minimum
        end
    end
end


text:
# ---
# title: 1392. Longest Happy Prefix
# id: problem1392
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String
# link: <https://leetcode.com/problems/longest-happy-prefix/description/>
# hidden: true
# ---
#
# A string is called a  _happy prefix_  if is a **non-empty** prefix which is
# also a suffix (excluding itself).
#
# Given a string `s`. Return the **longest happy prefix**  of `s` .
#
# Return an empty string if no such prefix exists.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "level"
#     Output: "l"
#     Explanation: s contains 4 prefix excluding itself ("l", "le", "lev", "leve"), and suffix ("l", "el", "vel", "evel"). The largest prefix which is also suffix is given by "l".
#
#
# **Example 2:**
#
#
#
#     Input: s = "ababab"
#     Output: "abab"
#     Explanation: "abab" is the largest prefix which is also suffix. They can overlap in the original string.
#
#
# **Example 3:**
#
#
#
#     Input: s = "leetcodeleet"
#     Output: "leet"
#
#
# **Example 4:**
#
#
#
#     Input: s = "a"
#     Output: ""
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 10^5`
#   * `s` contains only lowercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@discord_object struct MessageActivity
    type::MessageActivityType.MessageActivityTypeEnum
    party_id::String
end


text:
using BinaryProvider
using Compat
using Compat.Test

function pull_latest(url, dir)
    # Get the repo that we've already cloned, or clone a new one
    repo = if !isdir(dir)
        Compat.@info("Cloning $(basename(url))")
        LibGit2.clone(url, dir)
    else
        Compat.@info("Updating $(basename(url))")
        LibGit2.GitRepo(dir)
    end

    # Fetch latest changes
    LibGit2.fetch(repo)

    # Reset ourselves onto the latest `master`
    master_oid = LibGit2.GitCommit(repo, "remotes/origin/master")
    LibGit2.reset!(repo, master_oid, LibGit2.Consts.RESET_HARD)
    return nothing
end

# live_laugh_love(;irony=false)
function clone_build_test(builder_url, package_url, package_deps)
    # First, check out the builder and start building it
    builder_dir = abspath(joinpath("builder_repos", basename(builder_url)))
    pull_latest(builder_url, builder_dir)

    # Build for the current platform
    name, version, product_hashes = try
        cd(builder_dir) do
            Compat.@info("Building $(basename(builder_url))")
            m = Module(:__anon__)
            eval(m, quote
                ARGS = [$(triplet(platform_key()))]
                product_hashes = include(joinpath($(builder_dir), "build_tarballs.jl"))

                # Write out a build.jl file that points to this tarball
                bin_path = joinpath($(builder_dir), "products")
                BinaryBuilder.print_buildjl($(builder_dir), name, version, products(Prefix(bin_path)), product_hashes, bin_path)

                # Return back these three pieces of information
                return name, version, product_hashes
            end)
        end
    catch e
        display(e)
        warn("Building $(basename(builder_url)) failed.")
        return false
    end

    # Next, check out the package and its deps
    pkg_src_dirs = String[]
    for dep_url in package_deps
        dep_dir = abspath(joinpath("package_repos", basename(dep_url)))
        pull_latest(dep_url, dep_dir)
        push!(pkg_src_dirs, joinpath(dep_dir, "src"))
    end
    package_dir = abspath(joinpath("package_repos", basename(package_url)))
    pull_latest(package_url, package_dir)
    push!(pkg_src_dirs, joinpath(package_dir, "src"))
    pkg_env = merge(ENV, Dict("JULIA_LOAD_PATH" => join(pkg_src_dirs, ':')))

    # Copy over the new build.jl file and build it
    Compat.cp(joinpath(builder_dir, "products", "build_$(name).v$(version).jl"), joinpath(package_dir, "deps", "build.jl"); force=true)

    try
        Compat.@info("Building $(basename(package_url))")
        run(setenv(`$(Base.julia_cmd()) $(joinpath(package_dir, "deps", "build.jl"))`, pkg_env))
    catch e
       display(e)
        warn("Building $(basename(package_url)) failed.")
        return false
    end

    # Finally, test that package!
    try
        Compat.@info("Testing $(basename(package_url))")
        cd(joinpath(package_dir, "test")) do
            run(setenv(`$(Base.julia_cmd()) runtests.jl`, pkg_env))
        end
    catch e
        display(e)
        warn("Testing $(basename(package_url)) failed.")
        return false
    end
    return true
end

# Some test cases that ensure we can build and pass tests on a few different packages
test_cases = [
    ("https://github.com/staticfloat/OggBuilder", "https://github.com/staticfloat/Ogg.jl", ["https://github.com/JuliaIO/FileIO.jl"]),
    ("https://github.com/staticfloat/NettleBuilder", "https://github.com/staticfloat/Nettle.jl", []),
    # These disabled until they get updated with versioned outputs
    #("https://github.com/bicycle1885/ZlibBuilder", "https://github.com/bicycle1885/CodecZlib.jl", ["https://github.com/bicycle1885/TranscodingStreams.jl"]),
    #("https://github.com/JuliaWeb/MbedTLSBuilder", "https://github.com/JuliaWeb/MbedTLS.jl", []),
    # This segfaults for some reason.
    #("https://github.com/staticfloat/IpoptBuilder", "https://github.com/JuliaOpt/Ipopt.jl", ["https://github.com/JuliaOpt/MathProgBase.jl", "https://github.com/JuliaOpt/MathOptInterface.jl"]),
    # This awaiting the merging of https://github.com/davidanthoff/SnappyBuilder/pull/1
    #("https://github.com/davidanthoff/SnappyBuilder", "https://github.com/bicycle1885/Snappy.jl", []),
    # This awaiting https://github.com/dancasimiro/WAV.jl/pull/59
    #("https://github.com/staticfloat/FLACBuilder", "https://github.com/dmbates/FLAC.jl", ["https://github.com/staticfloat/Ogg.jl", "https://github.com/JuliaIO/FileIO.jl", "https://github.com/dancasimiro/WAV.jl"]),
]

mkpath("builder_repos")
mkpath("package_repos")

# Check everything out
@testset "Ecosystem tests" begin
    for (builder_url, package_url, package_deps) in test_cases
        @testset "$(basename(package_url))" begin
            @test clone_build_test(builder_url, package_url, package_deps)
        end
    end
end


text:
mutable struct Bits{N}
  offset::Int
  state::State
end

const Bit = Bits{1}

const states = WeakKeyDict()

function Bits{N}(vs::Integer...) where N
  b = Bits{N}(0, State(N, vs))
  states[b.state] = WeakKeyDict(b=>nothing)
  return b
end

Bits{N}() where N = Bits{N}(0)

bits(b::Bits) = b

nbits(::Bits{N}) where N = N

function Base.show(io::IO, b::Bits)
  println(io, "$(nbits(b)) bits with offset $(b.offset)")
  show(io, b.state)
end

function combine(a::State, b::State)
  a === b && return a
  of = nbits(a)
  a.ψ = a.ψ ⊗ b.ψ
  for buf in keys(states[b])
    buf.state = a
    buf.offset += of
    states[a][buf] = nothing
  end
  return a
end

combine(b::Bits) = b.state
combine(bs::Bits...) = reduce(combine, map(x -> x.state, bs))

bitindices(b::Bits{N}) where N = ntuple(i -> b.offset+i, Val(N))

tcat(x::Tuple) = x
tcat(x::Tuple, y::Tuple, ys::Tuple...) = tcat((x..., y...), ys...)

function apply!(U, bs::Bits...)
  s = combine(bs...)
  is = tcat(bitindices.(bs)...)
  apply!(U, s, is...)
  return
end

apply!(U, bs...) = apply!(U, bits.(bs)...)

measure(b::Bits) = measure!(b.state, bitindices(b)...)

measure(x) = repr(typeof(x),measure(bits(x)))


text:
"""
    weights = projection_weights(simplex)

This function implements Johnson's distance subalgorithm, as described in
E. G. Gilbert, D. W. Johnson, and S. S. Keerthi, “A fast procedure for
computing the distance between complex objects in three-dimensional space,”
1988. Given a simplex (a length N+1 vector of points of dimension N), it
returns weights such that dot(weights, simplex) yields the point in the convex
hull of the simplex which is closest to the origin.

This is the critical loop of the GJK algorithm, so it has been heavily optimized
to precompute, inline, and unroll as much of the algorithm as possible. For a
more readable (and much slower) implementation, see
projection_weights_reference()
"""
@generated function projection_weights(simplex::SVector{M, SVector{N, T}}) where {M,N,T}
    projection_weights_impl(simplex)
end

num_johnson_subsets(simplex_length::Integer) = 2^simplex_length - 1

"""
Compute all subsets of the points in the simplex in a reliable order. The
order is arbitrary, but was chosen to match the implemenation in
S. Cameron, “Enhancing GJK: computing minimum and penetration distances between
convex polyhedra,”. Specifically, subset i contains point j iff the binary
representation of i has a one at the jth bit.
"""
function johnson_subsets(simplex_length::Integer)
    num_subsets = num_johnson_subsets(simplex_length)
    subsets = falses(simplex_length, num_subsets)
    for i in 1:num_subsets
        for j in 1:simplex_length
            subsets[j, i] = (i >> (j - 1)) & 1
        end
    end
    subsets
end

function initial_deltas_impl(simplex_length, T)
    num_subsets = num_johnson_subsets(simplex_length)

    expr = quote
        deltas = zero(SVector{$num_subsets, SVector{$simplex_length, $T}})
    end

    # Set the weight of every singleton subset to 1.
    for i in 1:simplex_length
        elements = [:(zero($T)) for j in 1:simplex_length]
        elements[i] = :(one($T))
        arg = Expr(:call, :(SVector{$simplex_length, $T}), elements...)

        push!(expr.args, quote
            deltas = setindex(deltas, $(arg), $(2^(i - 1)))
        end)
    end
    expr
end

"""
This is the function which actually computes the weight values.
It returns an *expression* to compute that result, with all loops
unrolled and all indices pre-computed based on the size and dimension
of the simplex. For a more readable version of the same function, check
out projection_weights_reference().
"""
function projection_weights_impl(::Type{SVector{M, SVector{N, T}}}) where {M,N,T}
    num_subsets = num_johnson_subsets(M)
    subsets = johnson_subsets(M)
    complements = .!subsets

    expr = quote
        deltas = $(initial_deltas_impl(M, T))
    end

    for s in 1:(num_subsets - 1)
        k = first((1:M)[subsets[:,s]])
        push!(expr.args, quote
            viable = true
        end)
        for j in (1:M)[complements[:,s]]
            s2 = s + (1 << (j - 1))
            push!(expr.args, quote
                d = deltas[$s2][$j]
            end)
            for i in (1:M)[subsets[:,s]]
                push!(expr.args, quote
                    d += deltas[$s][$i] *
                    (dot(simplex[$i], simplex[$k]) - dot(simplex[$i], simplex[$j]))
                end)
            end
            push!(expr.args, quote
                deltas = setindex(deltas, setindex(deltas[$s2], d, $j), $s2)
            end)
            push!(expr.args, quote
                if d > 0
                    viable = false
                end
            end)
        end
        push!(expr.args, quote
            viable = viable && all($(Expr(:tuple,
            [:(deltas[$s][$i] >= 0) for i in (1:M)[subsets[:,s]]]...)))
            if viable
                return deltas[$s] ./ sum(deltas[$s])
            end
        end)
    end
    push!(expr.args, quote
        return deltas[$num_subsets] ./ sum(deltas[$num_subsets])
    end)
    return expr
end

function projection_weights_reference(simplex::SVector{M, SVector{N, T}}) where {M,N,T}
    num_subsets = 2^M - 1
    complements = falses(M, num_subsets)
    subsets = falses(M, num_subsets)
    for i in 1:num_subsets
        for j in 1:M
            subsets[j, i] = (i >> (j - 1)) & 1
            complements[j, i] = !subsets[j, i]
        end
    end

    deltas = zeros(MMatrix{M, num_subsets, T})
    viable = ones(MVector{num_subsets, Bool})

    for i in 1:M
        deltas[i, 2^(i - 1)] = 1
    end

    for s in 1:(num_subsets - 1)
        k = first((1:M)[subsets[:,s]])
        for j in (1:M)[complements[:,s]]
            s2 = s + (1 << (j - 1))
            delta_j_s2 = 0
            for i in (1:M)[subsets[:,s]]
                delta_j_s2 += deltas[i, s] * (dot(simplex[i], simplex[k]) - dot(simplex[i], simplex[j]))
            end
            if delta_j_s2 > 0
                viable[s] = false
            elseif delta_j_s2 < 0
                viable[s2] = false
            end
            deltas[j, s2] = delta_j_s2
        end
        if viable[s]
            return deltas[:,s] ./ sum(deltas[:,s])
        end
    end
    return deltas[:,end] ./ sum(deltas[:,end])
end


text:
# Copied from
# https://github.com/JuliaFolds/FLoops.jl/blob/master/test/test_executors.jl
# (for now).  TODO: Create FLoopsTesting.jl so that writing something
# like this easier?
module TestExecutors

using FLoops
using FoldsDagger
using Test

function f_sum(executor)
    @floop executor for x in 1:10
        @reduce(s += x)
    end
    return s
end

function f_filter_sum(executor)
    @floop executor for x in 1:10
        if isodd(x)
            @reduce(s += x)
        end
    end
    return s
end

function f_sum_nested_loop(executor)
    @floop executor for x in 1:10
        for y in 1:x
            @reduce(s += y)
        end
    end
    return s
end

function f_sum_update(executor)
    @floop executor for x in 1:10
        if isodd(x)
            @reduce(s += 2x)
        end
    end
    return s
end

function f_sum_op_init(executor)
    @floop executor for x in 1:10
        if isodd(x)
            @reduce(s = 0 + 2x)
        end
    end
    return s
end

function f_count_update(executor)
    @floop executor for x in 1:10
        if isodd(x)
            @reduce(s += 1)
        end
    end
    return s
end

function f_count_op_init(executor)
    @floop executor for x in 1:10
        if isodd(x)
            @reduce(s = 0 + 1)
        end
    end
    return s
end

function f_sum_continue(executor)
    @floop executor for x in 1:10
        x > 4 && continue
        @reduce(s += x)
    end
    return s
end

function f_sum_break(executor)
    @floop executor for x in 1:10
        @reduce(s += x)
        x == 3 && break
    end
    return s
end

function f_find_return(executor)
    @floop executor for x in 1:10
        @reduce(s += x)
        x == 3 && return (:found, x)
    end
    return s
end

function f_find_goto(executor)
    @floop executor for x in 1:10
        @reduce() do (s; x)
            s = x
        end
        x == 3 && @goto FOUND
    end
    return s
    @label FOUND
    return (:found, s)
end

function test_executors()
    @testset "$f" for (f, desired) in [
        (f_sum, 55),
        (f_filter_sum, 25),
        (f_sum_nested_loop, 220),
        (f_sum_update, 50),
        (f_sum_op_init, 50),
        (f_count_update, 5),
        (f_count_op_init, 5),
        (f_sum_continue, 10),
        (f_sum_break, 6),
        (f_find_return, (:found, 3)),
        (f_find_goto, (:found, 3)),
    ]
        @test f(SequentialEx()) === desired
        @test f(DaggerEx()) === desired
        @testset for basesize in 2:10
            @test f(DaggerEx(basesize = 2)) === desired
        end
    end
end

end  # module


text:
using Distributions
import Random: rand!, rand
import Distributions: logpdf
import Base: eltype, length
using GaussianDistributions

#===============================================================================
                        One dimensional random walk
===============================================================================#
abstract type RandomWalk end

mutable struct UniformRandomWalk{T} <: RandomWalk
    ϵ::T
    pos::Bool

    UniformRandomWalk(ϵ::T, pos=false) where T <: Number = new{T}(ϵ, pos)
end

eltype(::UniformRandomWalk{T}) where T = T
length(::UniformRandomWalk) = 1

function _rand(rw::UniformRandomWalk, θ, coord_idx)
    @assert truelength(coord_idx) == 1
    ϑ = thetainc(coord_idx, θ)[1]
    δ = rand(Uniform(-rw.ϵ, rw.ϵ))
    ϑ += rw.pos ? 0.0 : δ
    ϑ *= rw.pos ? exp(δ) : 1.0
    ϑ
end

function rand(rw::UniformRandomWalk, θ, coord_idx)
    ϑ = _rand(rw, θ, coord_idx)
    move_to_proper_place(ϑ, θ, coord_idx)
end

function rand!(rw::UniformRandomWalk, θ, coord_idx)
    ϑ = _rand(rw, θ, coord_idx)
    move_to_proper_place!(ϑ, θ, coord_idx)
end

function logpdf(rw::UniformRandomWalk, coord_idx, θ, θᵒ)
    @assert truelength(coord_idx) == 1
    !rw.pos && return 0.0
    ϑᵒ = thetainc(coord_idx, θᵒ)[1]
    -log(2.0*rw.ϵ)-log(ϑᵒ)
end

function readjust!(rw::UniformRandomWalk, accpt_track, param, ::Any, mcmc_iter, ::Any)
    δ = compute_δ(param, mcmc_iter)
    a_r = acceptance_rate(accpt_track)
    ϵ = compute_ϵ(rw.ϵ, param, a_r, δ)
    print("Updating random walker...\n")
    print("acceptance rate: ", round(a_r, digits=2),
          ", previous ϵ: ", round(rw.ϵ, digits=3),
          ", new ϵ: ", round(ϵ, digits=3), "\n")
    rw.ϵ = ϵ
    ϵ
end

get_dispersion(rw::UniformRandomWalk) = copy(rw.ϵ)
#===============================================================================
                    Multidimensional Gaussian random walk
===============================================================================#
mutable struct GaussianRandomWalk{T} <: RandomWalk
    Σ::Array{T,2}
    pos::Vector{Bool}

    function GaussianRandomWalk(Σ::Array{T,2}, pos=nothing) where T
        @assert size(Σ)[1] == size(Σ)[2]
        pos = (pos === nothing) ? fill(false, size(Σ)[1]) : pos
        new{T}(Σ, pos)
    end
end

eltype(::GaussianRandomWalk{T}) where T = T
length(rw::GaussianRandomWalk) = length(rw.pos)
remove_constraints!(rw::GaussianRandomWalk, ϑ) = (ϑ[rw.pos] = log.(ϑ[rw.pos]))
reimpose_constraints!(rw::GaussianRandomWalk, ϑ) = (ϑ[rw.pos] = exp.(ϑ[rw.pos]))

function _rand(rw::GaussianRandomWalk, θ, coord_idx)
    ϑ = [thetainc(coord_idx, θ)...]
    remove_constraints!(rw, ϑ)
    ϑᵒ = rand(Gaussian(ϑ, rw.Σ))
    reimpose_constraints!(rw, ϑᵒ)
    ϑᵒ
end

function rand(rw::GaussianRandomWalk, θ, coord_idx)
    ϑ = _rand(rw, θ, coord_idx)
    θᵒ = copy(θ)
    move_to_proper_place!(ϑ, θᵒ, coord_idx)
end

function rand!(rw::GaussianRandomWalk, θ, coord_idx)
    ϑ = _rand(rw, θ, coord_idx)
    move_to_proper_place!(ϑ, θ, coord_idx)
end

_logjacobian(rw::GaussianRandomWalk, ϑ, coord_idx) = -sum(log.(ϑ[rw.pos]))

function logpdf(rw::GaussianRandomWalk, coord_idx, θ, θᵒ)
    ϑ = [thetainc(coord_idx, θ)...]
    ϑᵒ = [thetainc(coord_idx, θᵒ)...]
    logJ = _logjacobian(rw, ϑᵒ, coord_idx)
    remove_constraints!(rw, ϑ)
    remove_constraints!(rw, ϑᵒ)
    logpdf(Gaussian(ϑ, rw.Σ), ϑᵒ) + logJ
end

function readjust!(rw::GaussianRandomWalk, ::Any, ::Any, cov_mat, ::Any, coord_idx)
    ρ = matinc(cov_mat, coord_idx)
    Σ = 2.38^2/length(rw)*ρ
    print("Updating multivariate random walker...\n")
    print("covariance: ", round.(ρ, digits=2),
          ", previous Σ: ", round.(rw.Σ, digits=3),
          ", new ϵ: ", round.(Σ, digits=3), "\n")
    rw.Σ = Σ
end

#===============================================================================
            Multidimensional mixture of Gaussian random walks
===============================================================================#

struct GaussianRandomWalkMix{T} <: RandomWalk
    gsn_A::GaussianRandomWalk{T}
    gsn_B::GaussianRandomWalk{T}
    pos::Vector{Bool}
    λ::Float64

    function GaussianRandomWalkMix(Σ_A::Array{T,2}, Σ_B::Array{T,2}, λ=0.5,
                                   pos=nothing) where T
        @assert 0.0 <= λ <= 1.0
        gsn_A = GaussianRandomWalk(Σ_A, pos)
        gsn_B = GaussianRandomWalk(Σ_B, pos)
        new{T}(gsn_A, gsn_B, pos, λ)
    end
end

eltype(::GaussianRandomWalkMix{T}) where T = T
length(rw::GaussianRandomWalkMix) = length(rw.pos)

function rand(rw::GaussianRandomWalkMix, θ, coord_idx)
    rw_i = pick_kernel(rw)
    rand(rw_i, θ, coord_idx)
end

function rand!(rw::GaussianRandomWalkMix, θ, coord_idx)
    rw_i = pick_kernel(rw)
    rand!(rw_i, θ, coord_idx)
end

function pick_kernel(rw::GaussianRandomWalkMix)
    rand(Bernoulli(rw.λ)) ? rw.gsn_B : rw.gsn_A
end

function logpdf(rw::GaussianRandomWalkMix, coord_idx, θ, θᵒ)
    log( (1-rw.λ)*exp(logpdf(rw.gsn_A, coord_idx, θ, θᵒ))
          + rw.λ *exp(logpdf(rw.gsn_B, coord_idx, θ, θᵒ)) )
end

function readjust!(rw::GaussianRandomWalkMix, ::Any, ::Any, cov_mat, ::Any, coord_idx)
    readjust!(rw.gsn_A, nothing, nothing, cov_mat, nothing, coord_idx)
    0.0 # this is just a temporary solution to prevent some unnecessary errors
end


text:
### A Pluto.jl notebook ###
# v0.11.13

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 7978f646-f25e-11ea-052d-2b6a98c1e3c9
begin
    using Markdown
    hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]));
    tip(text="md(これはTipです)") = Markdown.MD(Markdown.Admonition("correct", "Tip", [text]));
    md"下記のヒントを読んでね"
end

# ╔═╡ eaba0bfe-f253-11ea-0f40-b51436fe83ce
begin
    using Plots
    using PlutoUI
    using LaTeXStrings
end

# ╔═╡ 1a82d8e4-f255-11ea-013b-a12ca030a545
using Images, TestImages

# ╔═╡ f111e400-f25c-11ea-39c9-2ff328de00b0
md"# Pluto.jl 環境上での対話的ノートブック"

# ╔═╡ 90b3cc2e-f25f-11ea-2991-818b8dddbe36
md"Pluto.jl を用いたサンプルを紹介する．Jupyter notebook でも Interact.jl で実装はできるが， Pluto.jl によるノートブックの方が滑らかにスライダーの変化を楽しむことができる．"

# ╔═╡ 504a76a4-f25f-11ea-2d2b-7f7e2be7fbd0
tip(md"こんな感じのTipを載せることもできる. 次のHintの上にマウスをあててみよう. ぼやけた文字をくっきり写すことができるはずだ")

# ╔═╡ a7d73f04-f25e-11ea-0e3d-3b352eba207f
hint(md"Jupyter notebook のようにコードをまとめたブロックをセル (cell) と呼ぶが, セルに複数行のコードを記述するときは `begin ... end` ブロックで囲む必要がある. またこのHintはマークダウンで記述できる．これ自体の記述をどうしているかを見たい場合はこのセルの左に出てくる目マークを表しているアイコンをクリックすると表示できる．")

# ╔═╡ 7e366740-f260-11ea-131a-a7084a072758
tip(md"セルの右側にある再生ボタンを順次クリックするまたは Shift+Enter で進む. おそらくそんなことをしなくてもPlutoが全てのセルを実行してくれるはずではあるが")

# ╔═╡ 2fd22542-f25d-11ea-18bb-1713026f5026
md"## Plots"

# ╔═╡ e9f43d5e-f25d-11ea-0581-9d370db65eb7
md"""
### Slider を使ってグリグリする

```julia
@bind <変数名> <PlutoUIのオブジェクト>
```

のようにして UI を構築できる UI のコントロール（例えばスライダー）を変化させるとそれに対応する値が<変数名>に反映される.
"""

# ╔═╡ 7e8b6224-f254-11ea-090e-e1c6bd90bb33
begin
    t_slider = @bind t Slider(-π:0.01:π)
    md"t: $(t_slider) <- touch me"
end

# ╔═╡ a5246d68-f254-11ea-018e-79939c9689f5
begin
    x = range(-π,π,step=0.01)
    y = @. sin(x + 2t)
    plot(x, y, title= "sin(x+t)" * " where t=$t")
    xticks!(
        [-π,-π/2,0,π/2,π],
        [L"-\pi", L"-\pi/2", L"0", L"\pi/2", L"\pi"]
    )
end

# ╔═╡ 627a1d1a-f25d-11ea-38ea-19f2a7ee5bde
md"## アニメーションを gif で表示"

# ╔═╡ d4996e9c-f25c-11ea-3ad1-6b184dcdb79a
begin
    anim = @animate for i in 1:10 scatter(rand(10), rand(10), title="i=$i") end
    gif(anim, fps=2)
end

# ╔═╡ 00b739ec-f261-11ea-251e-e929f665952c
md"# 画像アプリケーションの例"

# ╔═╡ 416d498e-f259-11ea-3930-79118483e2c1
begin
    degree_slider = @bind degree Slider(0:360)
    md"degree: $(degree_slider)"
end

# ╔═╡ 0e957840-f25b-11ea-3cb1-334f405c5803
@bind name Select(TestImages.remotefiles)

# ╔═╡ 0dcce14c-f255-11ea-0348-ffddd69c792b
begin
    img = testimage(name)
    img1 = imrotate(img, deg2rad(degree)) # without cropping
    img2 = imrotate(img, deg2rad(degree), axes(img)) # with cropping
    img1 = imresize(img1, img2|>size)
    hcat(img1, img2)
end

# ╔═╡ Cell order:
# ╟─f111e400-f25c-11ea-39c9-2ff328de00b0
# ╟─90b3cc2e-f25f-11ea-2991-818b8dddbe36
# ╟─7978f646-f25e-11ea-052d-2b6a98c1e3c9
# ╟─504a76a4-f25f-11ea-2d2b-7f7e2be7fbd0
# ╟─a7d73f04-f25e-11ea-0e3d-3b352eba207f
# ╟─7e366740-f260-11ea-131a-a7084a072758
# ╟─2fd22542-f25d-11ea-18bb-1713026f5026
# ╠═eaba0bfe-f253-11ea-0f40-b51436fe83ce
# ╟─e9f43d5e-f25d-11ea-0581-9d370db65eb7
# ╠═7e8b6224-f254-11ea-090e-e1c6bd90bb33
# ╠═a5246d68-f254-11ea-018e-79939c9689f5
# ╟─627a1d1a-f25d-11ea-38ea-19f2a7ee5bde
# ╠═d4996e9c-f25c-11ea-3ad1-6b184dcdb79a
# ╟─00b739ec-f261-11ea-251e-e929f665952c
# ╠═1a82d8e4-f255-11ea-013b-a12ca030a545
# ╠═416d498e-f259-11ea-3930-79118483e2c1
# ╠═0e957840-f25b-11ea-3cb1-334f405c5803
# ╠═0dcce14c-f255-11ea-0348-ffddd69c792b


text:
srand(843)
using IterativeSolvers
include("mybicgstab.jl")

n=1000

#generate eigenvectors:
R=rand(Complex128,n,n)
(X1,R)=qr(R)

#generate eigenvalues
L1=zeros(n,n)
for i in 1:n
    L1[i,i]=i
end

#generate eigenvalue problem matrices
A=X1*L1*X1'

z=2.5+1.0*im
M=(z*eye(n)-A)

x=rand(n,1)
b=M*x

x0=zeros(n,1)
maxit=1000
eps=1e-5

xest=zbicgstabBlock(M,b,x0,maxit,eps)
#xest=bicgstabl(M,b,2,max_mv_products=500,tol=1e-5,initial_zero=true) #zbicgstabBlock(A,b,x0,maxit,eps)

res=norm(b-M*xest)/norm(b)
println("res=$res")

max_stars_repo_name: rleegates/feastjl

max_issues_repo_name: rleegates/feastjl

max_forks_repo_name: rleegates/feastjl


text:
export FELT, newFELT, NLDID, newNLDID

# Each application is characterized by:
# 1. (y,X) (well, formula, data, :i, : for the user)
# 2. Which cut points to use
# 3. How to interpolate the estimated function in between the cut points.
# If the outcome variable is discrete, then generally 2. will correspond to all points in the support of the dependent variable, and 3. will be a step function.
# If the outcome variable is continuous, then generally 2. will consist of some cleverly chosen points, and 3. will be a linear interpolation, or a B-spline, or something like that. I like linear interpolation because it is consistent with censoring.
mutable struct FELT
    formula
    data::DataFrame
    isymbol::Symbol
    tsymbol::Symbol
    y::Array{Float64,2}
    X::Array{Float64,2}
    y1s::Vector{Float64}
    y2s::Vector{Float64}
    b_hat::Vector{Float64}
    gamma_1_hat::Vector{Float64}
    gamma_2_hat::Vector{Float64}
    discrete::Bool
end

# Building on FELT, a nonlinear difference-in-differences structure
#   consists of:
#     1. a FELT object, containing control outcomes and estimation results
#     2. treatment indicator
#     3. treatment data
#     4. counterfactual outcomes.
#     5. Linear DiD estimate

mutable struct NLDID
    treatsymbol::Symbol
    felt::FELT
    t::Vector{Float64}
    y_treat::Vector{Float64}
    X_treat::Array{Float64,2}
    linearATT
    nonlinearATT
end

function newFELT(formula,data,isymbol,tsymbol; ys = false, discrete = true, levels = 10)

    # This only works for discrete!

    y, X = ConvertPanelToDiffCS_2(formula,data,isymbol,tsymbol)
    n,K = size(X)

    if ys !== false

        cuts1 = ys
        cuts2 = ys

    else

        if discrete
            cuts1 = sort(unique(y[:,1]))[2:end] #2:end: skip the first one
            cuts2 = sort(unique(y[:,2]))[2:end]

        else  #well then it must be continuous.

            # Compute `levels` equi-spaced quantiles.
            q_points = collect(0:(1/(levels+1)):1)[2:(end-1)]

            y1low = findmin(y[:,1])[1]
            y2low = findmin(y[:,2])[1]

            cuts1 = setdiff(quantile(y[:,1], q_points)[:,1],y1low)
            cuts2 = setdiff(quantile(y[:,2], q_points)[:,1],y2low)

        end

    end
    # Need these for extracting b, gamma1, gamma2
    n_1 = length(cuts1) - 1
    n_2 = length(cuts2)

    theta_hat = FixedEffectsOrderedLogit_2(y,X,cuts1,cuts2)
    b_hat = theta_hat[1:K]
    gamma_1_hat = [0;theta_hat[(K+1):(K+n_1)]]
    gamma_2_hat = theta_hat[(K+n_1+1):(K+n_1+n_2)]

    return FELT(formula,data,isymbol,tsymbol,
    y,X,cuts1,cuts2,
    b_hat,gamma_1_hat,gamma_2_hat,
    discrete)

end

function newNLDID(formula,data,isymbol,tsymbol,treatsymbol; ys = false, discrete = true, levels = 10)

    # First split
    data_0 = data[data[treatsymbol].==0,:]
    data_1 = data[data[treatsymbol].==1,:]

    # Create the FELT instance
    f = newFELT(formula,data_0,isymbol,tsymbol; ys = ys, discrete = discrete, levels = levels)

    # Extract the treat data.
    t = Vector(data_1[tsymbol])
    MF = ModelFrame(formula,data_1)
    X_treat = ModelMatrix(MF).m[:,2:end]
    y_treat = model_response(MF)

    y0 = y_treat[t .== 0]
    y1 = y_treat[t .== 1]
    X0 = X_treat[t .== 0,:]
    X1 = X_treat[t .== 1,:]

    if discrete
        yLO = minimum(y1)
        yHI = maximum(y1)
        y1_cf = counterfactual_discrete(y0,yLO,yHI,
                                        f.gamma_1_hat,f.y1s,
                                        f.gamma_2_hat,f.y2s)
        nlATT = [mean(y1)-y1_cf[2], mean(y1)-y1_cf[1]]
    else

        y1_cf = counterfactual_continuous(y0,X0,X1,f.b_hat,f.y1s,f.gamma_1_hat,f.y2s,f.gamma_2_hat,minimum(y_treat),maximum(y_treat))
        nlATT = mean(y1) - mean(y1_cf)
    end

    return NLDID(treatsymbol,f,t,y_treat,X_treat,"TBD",nlATT)
end


text:
# updateGammaz does not recognize this as bad diamond II
# Claudia April 2015

if !isdefined(:individualtest) individualtest = false; end

if(individualtest)
    include("../src/types.jl")
    include("../src/functions.jl")
end

if isdefined(:PhyloNetworks)
    PhyloNetworks.CHECKNET || error("need CHECKNET==true in PhyloNetworks to test snaq in test_correctLik.jl")
else
    CHECKNET || error("need CHECKNET==true to test snaq in test_correctLik.jl")
end

@testset "test: bad diamond, max pseudo lik" begin

tree = "(6,(5,#H7:0.0):9.970714072991349,(3,(((2,1):0.2950382234364404,4):0.036924483697671304)#H7:0.00926495670648208):1.1071489442240392);"
net = readTopologyLevel1(tree);
checkNet(net)
#printNodes(net)
#printEdges(net)
@test net.node[10].number == 3 # or: wrong hybrid
@test net.node[10].hybrid # or: does not know it is hybrid
@test net.node[10].isBadDiamondII # or: does not know it is bad diamond II
##plot(net,showEdgeNumber=true)
@test [e.inCycle for e in net.edge] == [-1,-1,3,3,-1,-1,-1,-1,-1,-1,3,3] # or: error in incycle
@test [e.containRoot for e in net.edge] == [true,true,false,true,true,false,false,false,false,false,false,true] # or: error in contain root
@test [e.istIdentifiable for e in net.edge] == [false,false,true,true,false,false,false,true,false,false,true,true] # or: istIdentifiable not correct
@test (net.edge[3].hybrid && net.edge[11].hybrid) # or: hybrid edges wrong")

df=DataFrame(t1=["1","1","2","2","1","2","2","2","2","2","1","2","2","3","2"],
             t2=["3","3","3","3","3","1","5","1","1","1","5","1","1","5","3"],
             t3=["5","5","5","5","6","5","6","3","6","5","6","3","3","6","6"],
             t4=["6","4","4","6","4","6","4","5","4","4","4","6","4","4","4"],
             CF1234=[0.565,0.0005,0.0005,0.565,0.00003,0.99986,0.0410167,1,0.99987,1,0.040167,0.998667,1,0.073167,0.00003],
             CF1324=[0.0903,0.8599,0.8599,0.0903,0.8885,0.00006,0.263,0,0.00006,0,0.2630,0.00006,0,0.0424667,0.8885])
df[:CF1423] = 1-df[:CF1234]-df[:CF1324]
d = readTableCF(df)

net2 = topologyMaxQPseudolik!(net,d, ftolRel=1e-5,ftolAbs=1e-6,xtolRel=1e-3,xtolAbs=1e-4)
@test 340 < net2.loglik < 340.5 # or: wrong loglik
@test net2.edge[3].istIdentifiable # or: wrong hybrid is t identifiable
@test 9.95 < net2.edge[3].length < 9.99 # or: wrong bl estimated
@test net2.edge[11].istIdentifiable # or: wrong hybrid is t identifiable
@test net2.edge[11].length < 0.01 # or: wrong bl estimated
@test net2.edge[10].length == 0 # or: tree edge in bad diamond II not 0
#printEdges(net2)

@test_nowarn show(DevNull, net2)
@test_nowarn [show(DevNull, net2.node[i]) for i in [1,3,10]];
@test_nowarn [show(DevNull, net2.edge[i]) for i in [1,3,11]];
@test_nowarn show(DevNull, d)
@test_nowarn show(DevNull, d.quartet[1])
@test_nowarn show(DevNull, d.quartet[1].qnet)
@test tipLabels(d.quartet) == ["1","2","3","4","5","6"]
a = (@test_nowarn fittedQuartetCF(d, :wide));
@test size(a) == (15,10)
a = (@test_nowarn fittedQuartetCF(d, :long));
@test size(a) == (45,7)

end

## testing readTopology----------------------------------------------------------------
## will remove this from the test because readTopology should not have to worry about istIdentifiable
## we move into updateGammaz
## tree = "(6,(5,#H7:0.0):9.970714072991349,(3,(((2,1):0.2950382234364404,4):0.036924483697671304)#H7:0.00926495670648208):1.1071489442240392);"
## net2 = readTopology(tree)
## printEdges(net2)

## for e in net2.edge
##     if(e.node[1].leaf || e.node[2].leaf)
##         !e.istIdentifiable || error("ext edge should not identifiable")
##     else
##         e.istIdentifiable || error("int edge should not identifiable")
##     end
## end

## ## this case works fine
## tree = "((((8,10))#H1,7),6,(4,#H1));" # Case I Bad diamond I
## net = readTopologyLevel1(tree)
## checkNet(net)
## net2 = readTopology(tree)
## printEdges(net2)

## for e in net2.edge
##     if(e.node[1].leaf || e.node[2].leaf)
##         !e.istIdentifiable || error("ext edge should not identifiable")
##     else
##         e.istIdentifiable || error("int edge should not identifiable")
##     end
## end


text:
module TestCategoricalAlgebra

using Test

@testset "ShapeDiagrams" begin
  include("ShapeDiagrams.jl")
end

@testset "FinSets" begin
  include("FinSets.jl")
end

@testset "Permutations" begin
  include("Permutations.jl")
end

end