text:
export get_F_all

function get_F_all(x::Array{Float64,2})
    dhp = zeros(Float64, 3, 3*(N+2))
    #Units are STILL not clear
    F = zeros(Float64, N, 3)
    get_f = get_F_au_au(x, F)
    F_au_au = [zeros(Float64, 2, 3); get_f]

    F_au_n_n = get_F_au_n_neutral(x)
    F_au_o_n = get_F_au_o_neutral(x)
    F_n_o_n = [get_F_n_o_neutral(x[1:2,:]); zeros(Float64, N, 3)]
    F_neutral = @. F_au_n_n + F_au_o_n + F_n_o_n + F_au_au #add half force of gold-gold?

    F_au_n_i = get_F_au_n_ion(x)
    F_au_o_i= get_F_au_o_ion(x)
    F_n_o_i = [get_F_n_o_ion(x[1:2,:]); zeros(Float64, N, 3)]
    zcom = zcom_NO(x[1:2, :])
    F_image_i = [get_F_image_ion(zcom); zeros(Float64, N, 3)]
    F_ion = @. F_au_n_i + F_au_o_i + F_n_o_i + F_image_i + F_au_au  #add half force of gold-gold?

    F_au_n_c = get_F_au_n_coup(x)
    F_au_o_c = get_F_au_o_coup(x)
    F_coupling = @. F_au_n_c + F_au_o_c

    #diabatic matrix elements of NO-Au Hamiltonian
    dhp[1, :] = vec(F_neutral')
    dhp[2, :] = vec(F_ion')
    dhp[3, :] = vec(F_coupling')
    return -dhp
end
#---------------------------Au-O Coupling----------------------------------
function get_F_au_n_coup_loop!(x::Array{Float64,2}, F::Array{Float64,2}, i::Int)
    delta_xn = @. x[i+2, :] - x[1, :]
    delta_xn[1] = @.delta_xn[1] - cell[1]*round(delta_xn[1]/cell[1])
    delta_xn[2] = @.delta_xn[2] - cell[2]*round(delta_xn[2]/cell[2])

    cutoff_test = @. abs(delta_xn) - cutoff
    bool_cutoff_test = @. cutoff_test <= 0.0
    if sum(bool_cutoff_test) != 0.0
        norm_delta_xn = sqrt(delta_xn' * delta_xn)
        if norm_delta_xn <= cutoff
            term1 = exp(coup_β_N*norm_delta_xn)
            term2 = (1.0 + coup_b_N*term1)^2 * norm_delta_xn
            term3 = coup_a_N * coup_b_N * coup_β_N/term2

            dV_dri = @. -term1*delta_xn*term3
            dV_drn = -dV_dri
            F[1, :] = F[1, :] - dV_drn
            F[i+2, :] = F[i+2, :] - dV_dri
        end
    end
    return F
end

function get_F_au_n_coup(x::Array{Float64,2})
    F = zeros(Float64, N+2, 3)
    for i in 1:N
        F = get_F_au_n_coup_loop!(x, F, i)
    end
    return F
end
#---------------------------Au-O Coupling----------------------------------
function get_F_au_o_coup_loop!(x::Array{Float64,2}, F::Array{Float64,2}, i::Int)
    delta_xo =  x[i+2, :] .- x[2, :]
    delta_xo[1] = @.delta_xo[1] - cell[1]*round(delta_xo[1]/cell[1])
    delta_xo[2] = @.delta_xo[2] - cell[2]*round(delta_xo[2]/cell[2])

    cutoff_test = @. abs(delta_xo) - cutoff
    bool_cutoff_test =  cutoff_test .<= 0.0
    if sum(bool_cutoff_test) != 0.0
        norm_delta_xo = sqrt(delta_xo' * delta_xo)
        if norm_delta_xo <= cutoff
            term1 = exp(coup_β_O*norm_delta_xo)
            term2 = (1.0 + coup_b_O*term1)^2 * norm_delta_xo
            term3 = coup_a_O * coup_b_O * coup_β_O/term2

            dV_dri = @. -term1*delta_xo*term3
            dV_dro = -dV_dri
            F[2, :] = F[2, :] - dV_dro
            F[i+2, :] = F[i+2, :] - dV_dri
        end
    end
    return F
end

function get_F_au_o_coup(x::Array{Float64,2})
    F = zeros(Float64, N+2, 3)
    for i in 1:N
        F = get_F_au_o_coup_loop!(x, F, i)
    end
    return F
end

#---------------------------Ionic N O ----------------------------------

function get_F_n_o_ion(x_no::Array{Float64,2})
    F = zeros(Float64, 2, 3)
    r_no = x_no[1, :] - x_no[2,:]
    norm_r_no = sqrt(r_no' * r_no)
    term = 2*F_i *δ_i * (exp(-2.0 * δ_i * (norm_r_no - rNO_e_i)) - exp(-δ_i*(norm_r_no - rNO_e_i)))
    F[1, :] = @. term * r_no/norm_r_no
    F[2, :] = - F[1, :]
    return F
end

#---------------------------Ionic Au O ----------------------------------
function get_F_au_o_ion_loop!(x::Array{Float64,2}, F::Array{Float64,2}, i::Int)
    delta_xo =  x[i+2, :] .- x[2, :]
    delta_xo[1] = @.delta_xo[1] - cell[1]*round(delta_xo[1]/cell[1])
    delta_xo[2] = @.delta_xo[2] - cell[2]*round(delta_xo[2]/cell[2])

    cutoff_test = @. abs(delta_xo) - cutoff
    bool_cutoff_test = @. cutoff_test <= 0.0
    if sum(bool_cutoff_test) != 0.0
        norm_delta_xo = sqrt(delta_xo' * delta_xo)
        if norm_delta_xo <= cutoff
            term = -A_n * exp(-α_n*norm_delta_xo) * α_n/norm_delta_xo
            dV_dri = term * delta_xo
            dV_dro = -dV_dri
            F[2, :] = @. F[2, :] - dV_dro
            F[i+2, :] = @. F[i+2, :] - dV_dri
        end
    end
    return F
end

function get_F_au_o_ion(x::Array{Float64,2})
    F = zeros(Float64, N+2, 3)
    for i in 1:N
        F = get_F_au_o_ion_loop!(x, F, i)
    end
    return F
end

#---------------------------Ionic Au N ----------------------------------

function get_F_au_n_ion_loop!(x::Array{Float64,2}, F::Array{Float64,2}, i::Int)
    delta_xn =  x[i+2, :] .- x[1, :]
    delta_xn[1] = @.delta_xn[1] - cell[1]*round(delta_xn[1]/cell[1])
    delta_xn[2] = @.delta_xn[2] - cell[2]*round(delta_xn[2]/cell[2])

    cutoff_test = @. abs(delta_xn) - cutoff
    bool_cutoff_test =  cutoff_test .<= 0.0
    if sum(bool_cutoff_test) != 0.0
        norm_delta_xn = sqrt(delta_xn' * delta_xn)
        if norm_delta_xn <= cutoff
            r_no = x[1, :] - x[2,:]
            norm_r_no = sqrt(r_no' * r_no)
            r_diff = norm_delta_xn - rN_sur_e_i
            term1 = r_no[3]
            term2 = exp(-β_i*r_diff)
            term3 = exp(-2.0 * β_i * r_diff)
            term4 = norm_r_no^4

            dV_xN_1 = 4.0*term2*(x[1, 1]-x[2, 1])*term1*term1
            dV_xN_1 = dV_xN_1/term4
            dV_xN_4 = 2.0*term3*delta_xn[1]*β_i
            dV_xN_4 = dV_xN_4/norm_delta_xn
            dV_xN_5 = 2.0*term2*delta_xn[1]*term1*term1*β_i
            dV_xN_5 = dV_xN_5/(norm_delta_xn*norm_r_no*norm_r_no)
            dV_xN = B_i*(dV_xN_1+dV_xN_4-dV_xN_5)

            dV_yN_1 = 4.0*term2*(x[1, 2]-x[2, 2])*term1*term1
            dV_yN_1 = dV_yN_1/term4
            dV_yN_4 = 2.0*term3*delta_xn[2]*β_i
            dV_yN_4 = dV_yN_4/norm_delta_xn
            dV_yN_5 = 2.0*term2*delta_xn[2]*term1*term1*β_i
            dV_yN_5 = dV_yN_5/(norm_delta_xn*norm_r_no*norm_r_no)
            dV_yN = B_i*(dV_yN_1+dV_yN_4-dV_yN_5)

            dV_zN_1 = 4.0*term2*(x[1, 3]-x[2, 3])*term1*term1
            dV_zN_1 = dV_zN_1/term4
            dV_zN_2 = 4.0*term2*term1
            dV_zN_2 = dV_zN_2/(norm_r_no*norm_r_no)
            dV_zN_4 = 2.0*term3*delta_xn[3]*β_i
            dV_zN_4 = dV_zN_4/norm_delta_xn
            dV_zN_5 = 2.0*term2*delta_xn[3]*term1*term1*β_i
            dV_zN_5 = dV_zN_5/(norm_delta_xn*norm_r_no*norm_r_no)
            dV_zN = B_i*(dV_zN_1+dV_zN_2+dV_zN_4-dV_zN_5)

            dV_xO_1 = 4.0*term2*(x[1, 1]-x[1, 2])*term1*term1
            dV_xO_1 = dV_xO_1/term4
            dV_xO = B_i*(-dV_xO_1)

            dV_yO_1 = 4.0*term2*(x[1, 1]-x[2, 2])*term1*term1
            dV_yO_1 = dV_yO_1/term4
            dV_yO = B_i*(-dV_yO_1)

            dV_zO_1 = 4.0*term2*(x[1, 3]-x[2, 3])*term1*term1
            dV_zO_1 = dV_zO_1/term4
            dV_zO_2 = 4.0*term2*term1
            dV_zO_2 = dV_zO_2/(norm_r_no*norm_r_no)
            dV_zO = B_i*(-dV_zO_1-dV_zO_2)

            dV_xi_3 = 2.0*term3*delta_xn[1]*β_i
            dV_xi_3 = dV_xi_3/norm_delta_xn
            dV_xi_4 = 2.0*term2*delta_xn[1]*term1*term1*β_i
            dV_xi_4 = dV_xi_4/(norm_delta_xn*norm_r_no*norm_r_no)
            dV_xi = B_i*(-dV_xi_3+dV_xi_4)

            dV_yi_3 = 2.0*term3*delta_xn[2]*β_i
            dV_yi_3 = dV_yi_3/norm_delta_xn
            dV_yi_4 = 2.0*term2*delta_xn[2]*term1*term1*β_i
            dV_yi_4 = dV_yi_4/(norm_delta_xn*norm_r_no*norm_r_no)
            dV_yi = B_i*(-dV_yi_3+dV_yi_4)

            dV_zi_3 = 2.0*term3*delta_xn[3]*β_i
            dV_zi_3 = dV_zi_3/norm_delta_xn
            dV_zi_4 = 2.0*term2*delta_xn[3]*term1*term1*β_i
            dV_zi_4 = dV_zi_4/(norm_delta_xn*norm_r_no*norm_r_no)
            dV_zi = B_i*(-dV_zi_3+dV_zi_4)

            #------   Calculating derivatives of cutoff function

            r_diff = cutoff - rN_sur_e_i

            term2 = exp(-β_i*r_diff)
            term3 = exp(-2.0*β_i*r_diff)

            dVc_xN_1 = 4.0*term2*(x[1, 1]-x[2, 1])*term1*term1
            dVc_xN_1 = dVc_xN_1/term4
            dVc_xN = B_i*(dVc_xN_1)

            dVc_yN_1 = 4.0*term2*(x[1, 2]-x[2, 2])*term1*term1
            dVc_yN_1 = dVc_yN_1/term4
            dVc_yN = B_i*(dVc_yN_1)

            dVc_zN_1 = 4.0*term2*(x[1, 3]-x[2, 3])*term1*term1
            dVc_zN_1 = dVc_zN_1/term4
            dVc_zN_2 = 4.0*term2*term1
            dVc_zN_2 = dVc_zN_2/(norm_r_no*norm_r_no)
            dVc_zN = B_i*(dVc_zN_1+dVc_zN_2)

            dVc_xO_1 = 4.0*term2*(x[1, 1]-x[2, 1])*term1*term1
            dVc_xO_1 = dVc_xO_1/term4
            dVc_xO = B_i*(-dVc_xO_1)

            dVc_yO_1 = 4.0*term2*(x[1, 2]-x[2, 2])*term1*term1
            dVc_yO_1 = dVc_yO_1/term4
            dVc_yO = B_i*(-dVc_yO_1)

            dVc_zO_1 = 4.0*term2*(x[1, 3]-x[2, 3])*term1*term1
            dVc_zO_1 = dVc_zO_1/term4
            dVc_zO_2 = 4.0*term2*term1
            dVc_zO_2 = dVc_zO_2/(norm_r_no*norm_r_no)
            dVc_zO = B_i*(-dVc_zO_1-dVc_zO_2)

            dVc_xi = 0.0
            dVc_yi = 0.0
            dVc_zi = 0.0

            F[1, 1] = F[1, 1] - dV_xN + dVc_xN
            F[1, 2] = F[1, 2] - dV_yN + dVc_yN
            F[1, 3] = F[1, 3] - dV_zN + dVc_zN

            F[2, 1] = F[2, 1] - dV_xO + dVc_xO
            F[2, 2] = F[2, 2] - dV_yO + dVc_yO
            F[2, 3] = F[2, 3] - dV_zO + dVc_zO

            F[i+2, 1] = F[i+2, 1] - dV_xi + dVc_xi
            F[i+2, 2] = F[i+2, 2] - dV_yi + dVc_yi
            F[i+2, 3] = F[i+2, 3] - dV_zi + dVc_zi
        end
    end
    return F
end

function get_F_au_n_ion(x::Array{Float64,2})
    F = zeros(Float64, N+2, 3)
    for i in 1:N
        F = get_F_au_n_ion_loop!(x, F, i)
    end
    return F
end

#---------------------------Ionic Image Potential----------------------------------
function get_F_image_ion(zcom::Float64)
    F = zeros(Float64, 2, 3)
    term = (D_i*(zcom - z0_i))/(C_i^2 + (zcom-z0_i)^2)^(1.5)
    F[1, 3] = term * mass_arr[1]/(mass_arr[1] + mass_arr[2])
    F[2, 3] = term * mass_arr[2]/(mass_arr[1] + mass_arr[2])
    return -F
end
#---------------------------Ionic N-O----------------------------------
function get_F_n_o_ion(x_no::Array{Float64,2})
    F = zeros(Float64, 2, 3)
    r_no = x_no[1, :] - x_no[2,:]
    norm_r_no = sqrt(r_no' * r_no)
    term = 2*F_i*δ_i*exp(-δ_i*(norm_r_no - rNO_e_i))*(1-exp(-δ_i*(norm_r_no - rNO_e_i)))
    F[1, :] = @. term * r_no/norm_r_no
    F[2, :] = -F[1, :]
    return -F
end

#---------------------------Neutral N-O----------------------------------

function get_F_n_o_neutral(x_no::Array{Float64,2})
    F = zeros(Float64, 2, 3)
    r_no = x_no[1, :] - x_no[2,:]
    norm_r_no = sqrt(r_no' * r_no)
    term = 2*F_n*δ_n*exp(-δ_n*(norm_r_no - r_0_NO))*(1-exp(-δ_n*(norm_r_no - r_0_NO)))
    F[1, :] = @. term * r_no/norm_r_no
    F[2, :] = -F[1, :]
    return -F
end

#---------------------------Neutral AU_O----------------------------------
function get_F_au_o_neutral_loop!(x_all::Array{Float64,2}, i::Int, F::Array{Float64,2})
    delta_xo =  x_all[i+2, :] .- x_all[2, :]
    delta_xo[1] = @. delta_xo[1] - cell[1]*round(delta_xo[1]/cell[1])
    delta_xo[2] = @. delta_xo[2] - cell[2]*round(delta_xo[2]/cell[2])

    cutoff_test = @. abs(delta_xo) - cutoff
    bool_cutoff_test = @. cutoff_test <= 0.0
    if sum(bool_cutoff_test) != 0.0
        norm_delta_xo = sqrt(delta_xo' * delta_xo)
        if norm_delta_xo <= cutoff
            term = -A_n * exp(-α_n * norm_delta_xo) * α_n/norm_delta_xo
            dV_dxi = @. term * delta_xo
            dV_dxO = -dV_dxi

            F[2,:] = F[2,:] - dV_dxO
            F[i+2,:] = F[i+2,:] - dV_dxi
        end
    end
    return F
end

function get_F_au_o_neutral(x_all::Array{Float64,2})
    F::Array{Float64,2} = zeros(Float64, N+2, 3) #Force acting on N and all Au atoms
    for i in 1:N
        F = get_F_au_o_neutral_loop!(x_all, i, F)
    end
    return F
end

#---------------------------Neutral AU_N----------------------------------
function get_F_au_n_neutral_loop!(x_all::Array{Float64,2}, i::Int, F::Array{Float64,2})
    delta_xn =  x_all[i+2, :] .- x_all[1, :]
    delta_xn[1] = @. delta_xn[1] - cell[1]*round(delta_xn[1]/cell[1])
    delta_xn[2] = @. delta_xn[2] - cell[2]*round(delta_xn[2]/cell[2])

    cutoff_test = @. abs(delta_xn) - cutoff
    bool_cutoff_test = @. cutoff_test <= 0.0
    if sum(bool_cutoff_test) != 0.0
        norm_delta_xn = sqrt(delta_xn' * delta_xn)
        if norm_delta_xn <= cutoff
            term = -B_n * exp(-β_n * norm_delta_xn) * β_n/norm_delta_xn
            dV_dxi = @. term * delta_xn
            dV_dxN = -dV_dxi

            F[1,:] = F[1,:] + dV_dxN
            F[i+2,:] = F[i+2,:] + dV_dxi
        end
    end
    return F
end

function get_F_au_n_neutral(x_all::Array{Float64,2})
    F::Array{Float64,2} = zeros(Float64, N+2, 3) #Force acting on N and all Au atoms
    for i in 1:N
        F = get_F_au_n_neutral_loop!(x_all, i, F)
    end
    return -F
end

#---------------------------GOLD AU AU LATTICE----------------------------------
function F_au_au_loop_if!(x::Array{Float64,2}, F::Array{Float64,2}, i::Int, j::Int)
    if nn[i, j] != 0
        m::Int = nn[i, j]+2
        xm = x[m, :]
        xi = x[i+2, :]
        temp = xm .- xi
        temp .= temp./cell
        temp_floor = floor.(temp .+ 0.5)
        temp .=  temp .- temp_floor
        r = temp .* cell .- r0[:, j]

        if j == 1 || j== 7
            F[i,:] = F[i,:]  .+ (d1_new + d1_new') * r
        elseif j == 2 || j== 8
            F[i,:] = F[i,:]  .+  (d2_new + d2_new') * r
        elseif j == 3 || j== 9
            F[i,:] = F[i,:]  .+  (d3_new + d3_new') * r
        elseif j == 4 || j== 10
            F[i,:] = F[i,:]  .+  (d4_new + d4_new') * r
        elseif j == 5|| j== 11
            F[i,:] = F[i,:]  .+  (d5_new + d5_new') * r
        elseif j == 6 || j== 12
            F[i,:] = F[i,:]  .+  (d6_new + d6_new') * r
        end
    end
    return F
end

function F_au_au_loop!(x_au::Array{Float64,2}, F::Array{Float64,2})
    for i in 1:N
        for j in 1:12
            F = F_au_au_loop_if!(x_au, F, i, j)
        end
    end
    return F
end

function get_F_au_au(x_au::Array{Float64,2}, F::Array{Float64,2})
    #Force acting on all N AU atoms
    F = F_au_au_loop!(x_au, F)
    return -F
end


text:
struct ReplacingTraces <: AbstractTraceUpdate end
struct AccumulatingTraces <: AbstractTraceUpdate end

_accumulate_trace(::AccumulatingTraces, e::AbstractMatrix, s::AbstractVector, pred_inds::Nothing) =
    e .+= s'
_accumulate_trace(::AccumulatingTraces, e::AbstractMatrix, s::AbstractVector, pred_inds) =
    e[pred_inds,:] .+= s'

_accumulate_trace(::AccumulatingTraces, e::AbstractMatrix, s::AbstractMatrix, pred_inds::Nothing) =
    e .+= s
_accumulate_trace(::AccumulatingTraces, e::AbstractMatrix, s::AbstractMatrix, pred_inds) = begin
    e[pred_inds,:] .+= s
end

_accumulate_trace(::ReplacingTraces, e::AbstractMatrix, s::AbstractVector, pred_inds::Nothing) =
    e .= s'
_accumulate_trace(::ReplacingTraces, e::AbstractMatrix, s::AbstractVector, pred_inds) =
    e[pred_inds,:] .= s'

_accumulate_trace(::ReplacingTraces, e::AbstractMatrix, s::AbstractMatrix, pred_inds::Nothing) =
    e .= s
_accumulate_trace(::ReplacingTraces, e::AbstractMatrix, s::AbstractMatrix, pred_inds) =
    e[pred_inds,:] .= s

function update_trace!(t::AbstractTraceUpdate, e::AbstractMatrix, s, λ, discounts, ρ, pred_inds=nothing; emphasis=nothing)
    e .*= λ * discounts .* ρ

    if !isnothing(emphasis)
        emphasized_s = emphasis * s'
        _accumulate_trace(t, e, emphasized_s, pred_inds)
    else
        _accumulate_trace(t, e, s, pred_inds)
    end
end

function get_demon_pis(horde::GVFSRHordes.SRHorde, num_actions, state, obs)
    target_pis = zeros(length(horde), num_actions)
    for i in 1:length(horde.PredHorde)
        for a in 1:num_actions
            target_pis[i, a] = get(GVFHordes.policy(horde.PredHorde.gvfs[i]); state_t = obs, action_t = a)
        end
    end

    state_action_feature_length = Int( length(horde.SFHorde) / horde.num_SFs)
    inds_per_task = collect(1:state_action_feature_length:length(horde.SFHorde))

    for a in 1:num_actions
        unique_pi_SF = map(gvf -> get(GVFHordes.policy(gvf); state_t=obs, action_t=a), horde.SFHorde.gvfs[inds_per_task])
        upsf = repeat(unique_pi_SF, inner = state_action_feature_length)
        target_pis[length(horde.PredHorde)+1:end, a] .=  upsf
    end

    return target_pis
end

function get_demon_pis(horde, num_actions, state, obs)
    target_pis = zeros(length(horde), num_actions)
    for i in 1:length(horde)
        for a in 1:num_actions
            target_pis[i, a] = get(GVFHordes.policy(horde.gvfs[i]); state_t = obs, action_t = a)
        end
    end
    return target_pis
end

function get_demon_parameters(lu::LearningUpdate, learner, demons, obs, state, action, next_obs, next_state, next_action, env_reward)
    C, next_discounts, _ = get(demons; state_t = obs, action_t = action, state_tp1 = next_obs, action_tp1 = next_action, reward = env_reward)
    target_pis = get_demon_pis(demons, learner.num_actions, state, obs)
    next_target_pis = get_demon_pis(demons, learner.num_actions, next_state, next_obs)
    C, next_discounts, target_pis, next_target_pis
end

function get_interest(learner, lu, obs, action)
    # An easy way to swap around the interest function easily.
    # This might want to be included into the agent config at some point,
    # but currently we can just leave it as is.

    # 1 for all states
    return ones(learner.num_demons)

    # if (obs[1] == 14.0)
    #     return ones(learner.num_demons)
    # else
    #     return zeros(learner.num_demons)
    # end
end

ext: jl

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
# This requires a little more understanding of acting along certain axes of an
# array.  The only useful references I can find for that are
# https://julialang.org/blog/2016/02/iteration/#filtering_along_a_specified_dimension_exploiting_multiple_indexes
# and
# https://discourse.julialang.org/t/crazy-allocations-using-cartesianindices/42262/2

const Rotator = Union{Quaternion, Rotor}

function _unflip!(q, Rpre, R, Rpost)
    @inbounds for Ipost in Rpost
        for i in R
            for Ipre in Rpre
                if q[Ipre, i-1, Ipost] ⋅ q[Ipre, i, Ipost] < 0
                    q[Ipre, i, Ipost] *= -1
                end
            end
        end
    end
end

function _unflip!(q, p, Rpre, R, Rpost)
    @inbounds for Ipost in Rpost
        let i = 1
            for Ipre in Rpre
                p[Ipre, i, Ipost] = q[Ipre, i, Ipost]
            end
        end
        for i in R
            for Ipre in Rpre
                if p[Ipre, i-1, Ipost] ⋅ q[Ipre, i, Ipost] < 0
                    p[Ipre, i, Ipost] = -q[Ipre, i, Ipost]
                else
                    p[Ipre, i, Ipost] = q[Ipre, i, Ipost]
                end
            end
        end
    end
end

function unflip!(q::AbstractArray{<:AbstractQuaternion}; dim::Integer=1)
    Rpre = CartesianIndices(size(q)[1:dim-1])
    R = 2:size(q, dim)
    Rpost = CartesianIndices(size(q)[dim+1:end])
    _unflip!(q, Rpre, R, Rpost)
end

"""
    unflip(q, [dim=1])
    unflip!(q, [dim=1])

Flip the signs of successive quaternions along dimension `dim` so that they are
as continuous as possible.

If `q` represents a series of rotations, the sign of each element is arbitrary.
However, for certain purposes — such as interpolation and differentiation — the
continuity of the quaternions matters, and so we want the *quaternions* to be
as continuous as possible without changing the *rotations* that they represent.

# Examples
```jldoctest
julia> q = [imx, -imx, imx, -imx];

julia> unflip(q)
4-element Vector{QuatVec{Int64}}:
 0 + 1𝐢 + 0𝐣 + 0𝐤
 0 + 1𝐢 + 0𝐣 + 0𝐤
 0 + 1𝐢 + 0𝐣 + 0𝐤
 0 + 1𝐢 + 0𝐣 + 0𝐤
```
"""
function unflip(q::AbstractArray{<:AbstractQuaternion}; dim::Integer=1)
    p = similar(q)
    Rpre = CartesianIndices(size(q)[1:dim-1])
    R = 2:size(q, dim)
    Rpost = CartesianIndices(size(q)[dim+1:end])
    _unflip!(q, p, Rpre, R, Rpost)
    p
end

"""
    slerp(q₁, q₂, τ; [unflip=false])

"Spherical Linear intERPolation" of a pair of quaternions.

The result of a "slerp" is given by

        (q₂ / q₁)^τ * q₁

When `τ` is 0, this evaluates to `q₁`; when `τ` is 1, this evaluates to `q₂`;
for any other values the result varies between the two.

Note that applying this to successive pairs of quaternions as in `slerp(q₁, q₂,
τₐ)` and `slerp(q₂, q₃, τᵦ)` will be continuous, but the derivative will be
discontinuous when moving from the first pair to the second.  See
[`squad`](@ref) for a more continuous curve.

If `unflip=true` is passed as a keyword, and the input quaternions are more
anti-parallel than parallel, the sign of `q₂` will be flipped before the result
is computed.

See also [`slerp∂slerp∂τ`](@ref), to simultaneously evaluate this function and
its derivative with respect to `τ`, or [`slerp∂slerp`](@ref) to evaluate this
function and its derivative with respect to each parameter of the input.

"""
function slerp(q₁::R1, q₂::R2, τ::Real; unflip::Bool=false) where {R1<:Rotator, R2<:Rotator}
    if unflip && q₁⋅q₂ < 0
        return (-q₂ / q₁)^τ * q₁
    end
    (q₂ / q₁)^τ * q₁
end

@doc raw"""
    squad_control_points(R::AbstractVector{Rotor}, t::AbstractVector{<:AbstractFloat}, i::Int)

This is a helper function for the `squad` routines, returning the control
points between one pair of input rotors.

The expressions for ``A`` and ``B`` (assuming all indices are valid) are
```math
\begin{aligned}
A_{i} &= R_{i}\, \exp\left\{
  \frac{1}{4}
  \left[
    \log\left(\bar{R}_{i-1}\, R_i\right) \frac{t_{i+1} - t_{i}} {t_{i} - t_{i-1}}
    - \log\left(\bar{R}_{i}\, R_{i+1}\right)
  \right]
\right\},
\\
B_{i} &= R_{i+1}\, \exp\left\{
  -\frac{1}{4}
  \left[
    \log\left(\bar{R}_{i+1}\, R_{i+2}\right) \frac{t_{i+1} - t_{i}} {t_{i+2} - t_{i+1}}
    - \log\left(\bar{R}_{i}\, R_{i+1}\right)
  \right]
\right\}.
\end{aligned}
```
These expressions will be invalid for ``A_{\mathrm{begin}}``, ``A_{\mathrm{end}}``,
``B_{\mathrm{end-1}}``, and ``B_{\mathrm{end}}``, because they all involve out-of-bounds indices of
``R_i``.  We can simply extend the input `R` values by linearly extrapolating, which results in the
following simplified results:
```math
\begin{aligned}
A_{\mathrm{begin}} &= R_{\mathrm{begin}} \\
A_{\mathrm{end}} &= R_{\mathrm{end}} \\
B_{\mathrm{end-1}} &= R_{\mathrm{end}} \\
B_{\mathrm{end}} &= R_{\mathrm{end}}\, \bar{R}_{\mathrm{end-1}}\, R_{\mathrm{end}} \\
                 &= 2\left(R_{\mathrm{end}}\cdot R_{\mathrm{end-1}}\right)\, R_{\mathrm{end}} - R_{\mathrm{end-1}}.
\end{aligned}
```
"""
function squad_control_points(R::AbstractVector{<:Rotor}, t::AbstractVector{<:AbstractFloat}, i::Int)
    if i==1
        A = R[1]
    elseif i==length(R)
        A = R[end]  # COV_EXCL_LINE
    else
        A = R[i] * exp(
            (
                log(conj(R[i-1]) * R[i]) * ((t[i+1] - t[i]) / (t[i] - t[i-1]))
                - log(conj(R[i]) * R[i+1])
            ) / 4
        )
    end
    if i<length(R)-1
        B = R[i+1] * exp(
            (
                log(conj(R[i+1]) * R[i+2]) * ((t[i+1] - t[i]) / (t[i+2] - t[i+1]))
                - log(conj(R[i]) * R[i+1])
            ) / -4
        )
    elseif i==length(R)-1
        B = R[i+1]
    else # i==length(R)
        B = Rotor{eltype(R)}((2*(R[i] ⋅ R[i-1]) * R[i] - R[i-1]).components...)  # COV_EXCL_LINE
    end
    A, B
end

const unflip_func = unflip  # `unflip` will be used as a local variable in
                            # `squad`, but the function will also be needed

"""
    squad!(Rout, Ω⃗out, Ṙout, Rin, tin, tout; [unflip=false], [validate=false])
    squad!(Rout, Rin, tin, tout; [unflip=false], [validate=false])

In-place evaluation of "Spherical QUADrangular interpolation".  Note that this
is intended mostly as a utility function; the [`squad`](@ref) is more
user-friendly.  However, for efficiency, this function may be preferable.

The first three arrays will be modified in place, and must have the same length
as `tout`.  Their elements must be `Rotor`, `QuatVec`, and `Quaternion`,
respectively.  Optionally, either or both of `Ω⃗out` and `Ṙout` maybe `nothing`,
in which case they will not be computed.

See also [`squad`](@ref).

"""
function squad!(
    Rout::AbstractVector{<:Rotor}, Ω⃗out::Union{Nothing, AbstractVector{<:QuatVec}}, Ṙout::Union{Nothing, AbstractVector{<:Quaternion}},
    Rin::AbstractVector{<:Rotor}, tin::AbstractVector{<:Real}, tout::AbstractVector{<:Real}; unflip=false, validate=false
)
    if length(tout) == 0
        return
    end
    t_begin, t_end = extrema(tin)
    @assert t_begin < t_end  # Proves that there are at least 2 tin
    @assert length(Rin) == length(tin)  # Proves that there are at least 2 Rin
    @assert length(Rout) == length(tout)
    evaluate_Ω⃗ = (Ω⃗out !== nothing)
    evaluate_Ṙ = (Ṙout !== nothing)
    if evaluate_Ω⃗
        @assert length(Ω⃗out) == length(tout)
    end
    if evaluate_Ṙ
        @assert length(Ṙout) == length(tout)
    end
    if validate
        @assert minimum(diff(tin)) > 0
        if length(tout) > 1
            @assert minimum(diff(tout)) > 0
            tout_begin, tout_end = extrema(tout)
            @assert t_begin ≤ tout_begin
            @assert tout_end ≤ t_end
        elseif length(tout) == 1
            @assert t_begin ≤ tout[1] ≤ t_end
        end
    end
    if unflip
        Rin = unflip_func(Rin)
    end
    j = 1
    toutj = tout[j]
    if toutj == tin[1]
        if evaluate_Ω⃗ || evaluate_Ṙ
            i = 1
            A, B = squad_control_points(Rin, tin, i)
            ta, tb = tin[i], tin[i+1]
            qᵢ, qᵢ₊₁ = Rin[i], Rin[i+1]
            s, ∂s∂t = squad∂squad∂t(qᵢ, A, B, qᵢ₊₁, ta, tb, tout[j])
            if evaluate_Ω⃗
                Ω⃗out[j] = 2 * eltype(Ω⃗out)(∂s∂t / s)
            end
            if evaluate_Ṙ
                Ṙout[j] = ∂s∂t
            end
        end
        Rout[1] = Rin[1]
        j += 1
    end
    while j ≤ length(Rout)
        toutj = tout[j]
        i = searchsortedfirst(tin, toutj) - 1
        if i == length(tin)
            error("Searching for $toutj went out of range [$t_begin, $t_end]")  # COV_EXCL_LINE
        end
        A, B = squad_control_points(Rin, tin, i)
        ta, tb = tin[i], tin[i+1]
        qᵢ, qᵢ₊₁ = Rin[i], Rin[i+1]
        while j ≤ length(Rout) && tb ≥ tout[j]
            if evaluate_Ω⃗ || evaluate_Ṙ
                s, ∂s∂t = squad∂squad∂t(qᵢ, A, B, qᵢ₊₁, ta, tb, tout[j])
                Rout[j] = s
                if evaluate_Ω⃗
                    Ω⃗out[j] = 2 * eltype(Ω⃗out)(∂s∂t / s)
                end
                if evaluate_Ṙ
                    Ṙout[j] = ∂s∂t
                end
            else
                τ = (tout[j] - ta) / (tb - ta)
                Rout[j] = slerp(
                    slerp(qᵢ, qᵢ₊₁, τ),
                    slerp(A, B, τ),
                    2τ*(1-τ)
                )
            end
            j += 1
        end
    end
end

"""
    squad(Rin, tin, tout; [kwargs...])

"Spherical QUADrangle interpolation" of the input `Rotor`s `Rin` with
corresponding times `tin`, to the output times `tout`.

This is a slightly generalized version of [Shoemake's "spherical Bézier
curves"](https://doi.org/10.1145/325165.325242), to allow for time steps of
varying sizes.

The input `Rin` and `tin` must be vectors of the same length.  The output
`tout` may be either a single real number or a vector of real numbers.  Both
`tin` and `tout` are assumed to be sorted, and `tout` is assumed to be
contained entirely within `tin`; no extrapolation will be done.

See also [`squad!`](@ref) for in-place versions of this function.

# Keyword arguments

If `unflip=true` is passed as a keyword, the [`unflip`](@ref) function will be
applied to `Rin`.

If `validate=true` is passed as a keyword, the time ordering of the input `tin`
and `tout` will be tested to ensure that no extrapolation will be done.

If `compute_angular_velocity=true` is passed as a keyword, the return value
will be a tuple.  The first element of the tuple will be a vector of `Rotor`s
as before, but the second element will be a vector of `QuatVec`s representing
the angular velocity.

If `compute_derivative=true` is passed as a keyword, the return value will be
a tuple.  The first element of the tuple will be a vector of `Rotor`s as
before, but the last element will be a vector of `Quaternion`s representing the
time-derivative of the rotors.  Note that if `compute_angular_velocity=true`,
this tuple will have three elements.

"""
function squad(
        Rin::AbstractVector{Rotor{T}}, tin::AbstractVector{<:Real}, tout::AbstractVector{<:Real};
        unflip=false, validate=false, compute_angular_velocity=false, compute_derivative=false
) where {T}
    Rout_eltype = promote_type(eltype(tin), eltype(tout), T)
    Rout = similar(Rin, Rotor{Rout_eltype}, length(tout))
    Ω⃗out = compute_angular_velocity ? similar(Rin, QuatVec{Rout_eltype}, length(tout)) : nothing
    Ṙout = compute_derivative ? similar(Rin, Quaternion{Rout_eltype}, length(tout)) : nothing
    squad!(Rout, Ω⃗out, Ṙout, Rin, tin, tout; unflip=unflip, validate=validate)
    if compute_angular_velocity && compute_derivative
        return (Rout, Ω⃗out, Ṙout)
    elseif compute_angular_velocity
        return (Rout, Ω⃗out)
    elseif compute_derivative
        return (Rout, Ṙout)
    end
    return Rout
end

function squad(
        Rin::AbstractVector{Rotor{T}}, tin::AbstractVector{<:Real}, tout::Real;
        unflip=false, validate=false, compute_angular_velocity=false, compute_derivative=false
) where {T}
    result = squad(Rin, tin, [tout]; unflip, validate, compute_angular_velocity, compute_derivative)
    if compute_angular_velocity && compute_derivative
        return (result[1][1], result[2][1], result[3][1])
    elseif compute_angular_velocity
        return (result[1][1], result[2][1])
    elseif compute_derivative
        return (result[1][1], result[2][1])
    end
    return result[1]
end


text:
# This experiment will compare the behavior of pvalue, effect size, and
# product over a wide range of screen and library designs

function compare_methods(filepath; debug=false, quiet=false)

    runs = []

    if !debug
        num_runs = 25
        levels = (10, 100, 1000)
    else
        num_runs = 1
        levels = (10)
    end

    screen_combos = (((FacsScreen, NaN), (GrowthScreen, 10), (GrowthScreen, 20)),
                    levels,
                    (CRISPRi(), CRISPRn()),
                    (Linear(), Sigmoidal()))

    for (screentype, representation, crisprtype, genetype) in IterTools.product(screen_combos...)
        screen = screentype[1]()
        screen.representation = representation
        screen.bottleneck_representation = representation
        screen.seq_depth = representation
        if screentype[1] == FacsScreen
            screen.bin_info = OrderedDict{Symbol, Tuple{Float64, Float64}}(:bin1 => (0.0, 0.25), :bin2 => (0.75, 1.0))
            max_phenotype_dists = Dict{Symbol, Tuple{Float64, Sampleable}}(
                :inactive => (0.75, Delta(0.0)),
                :negcontrol => (0.05, Delta(0.0)),
                :increasing => (0.1, TruncatedNormal(0.55, 0.2, 0.1, 1)),
                :decreasing => (0.1, TruncatedNormal(-0.55, 0.2, -1, -0.1))
            )
        else
            screen.num_bottlenecks = screentype[2]
            max_phenotype_dists = Dict{Symbol, Tuple{Float64, Sampleable}}(
                :inactive => (0.85, Delta(0.0)),
                :negcontrol => (0.05, Delta(0.0)),
                :decreasing => (0.1, TruncatedNormal(-0.55, 0.2, -1, -0.1))
            )
        end

        kd_phenotype_relationships = Dict{Symbol, Tuple{Float64, KDPhenotypeRelationship}}(
            Symbol(typeof(genetype)) => (1.0, genetype)
        )

        if typeof(crisprtype) == CRISPRn
            knockdown_dist = Dict{Symbol, Tuple{Float64, Sampleable}}(
                    :high => (0.9, Delta(1.0)),
                    :low => (0.1, TruncatedNormal(0.05, 0.07, 0, 1))
                )
        else
            knockdown_dist = Dict{Symbol, Tuple{Float64, Sampleable}}(
                :high => (0.9, TruncatedNormal(0.90, 0.1, 0, 1)),
                :low => (0.1, TruncatedNormal(0.05, 0.07, 0, 1))
            )
        end

        lib = Library(knockdown_dist, max_phenotype_dists, kd_phenotype_relationships, crisprtype)

        for i in 1:num_runs
            push!(runs, (screen, lib, i))
        end
    end

    const overlap = intersect(fieldnames(FacsScreen), fieldnames(GrowthScreen))
    # custom function for handling both growth and a facs screen in the
    # same relational datastructure
    flatten_overlap = (setup, lib) -> begin
        local results = Any[typeof(setup)]
        for name in overlap
            push!(results, getfield(setup, name))
        end
        push!(results, typeof(setup) == GrowthScreen ? getfield(setup, :num_bottlenecks) : NaN)
        push!(results, as_array(lib)...)
        results
    end

    compare_reduction_methods = (bc_counts, genes) -> begin
        sig, noi = signal(bc_counts), noise(bc_counts)

        pvalmeanprod = auprc(abs.(genes[:pvalmeanprod_bin2_div_bin1]), genes[:class], Set([:increasing, :decreasing]))[1]
        pvalue = auprc(genes[:pvalue_bin2_div_bin1], genes[:class], Set([:increasing, :decreasing]))[1]
        effectsize = auprc(genes[:absmean_bin2_div_bin1], genes[:class], Set([:increasing, :decreasing]))[1]

        (sig, noi, pvalmeanprod, pvalue, effectsize)
    end

    before = time()
    results = pmap(args -> run_exp(args[1],
                                   args[2],
                                   compare_reduction_methods;
                                   flatten_func=flatten_overlap,
                                   run_idx=args[3]), runs)
    (!quiet) && println("$(time() - before) seconds")
    results = DataFrame(permutedims(hcat(results...), [2,1]))
    new_names = [:method; :score; :screentype; overlap...; :num_bottlenecks; array_names(Library)...; :run_idx]
    hierarchy = reshape([:signal, :noise, :product, :pvalue, :effectsize], 5, 1)
    results = construct_hierarchical_label(hierarchy, results, new_names)
    CSV.write(filepath, results)
end


text:
##############################
#
# obj-Files
#
##############################

function load(io::Stream{format"OBJ"}; facetype=GLTriangleFace,
              pointtype=Point3f0, normaltype=Vec3f0, uvtype=Vec2f0)

    points, v_normals, uv, faces = pointtype[], normaltype[], uvtype[], facetype[]
    f_uv_n_faces = (faces, GLTriangleFace[], GLTriangleFace[])
    last_command = ""
    attrib_type  = nothing
    for full_line in eachline(stream(io))
        # read a line, remove newline and leading/trailing whitespaces
        line = strip(chomp(full_line))
        !isascii(line) && error("non valid ascii in obj")

        if !startswith(line, "#") && !isempty(line) && !all(iscntrl, line) #ignore comments
            lines = split(line)
            command = popfirst!(lines) #first is the command, rest the data

            if "v" == command # mesh always has vertices
                push!(points, Point{3, Float32}(parse.(Float32, lines))) # should automatically convert to the right type in vertices(mesh)
            elseif "vn" == command
                push!(v_normals, Vec3f0(parse.(Float32, lines)))
            elseif "vt" == command
                if length(lines) == 2
                    push!(uv, Vec2f0(parse.(Float32, lines)))
                elseif length(lines) == 3
                    push!(uv, Vec3f0(parse.(Float32, lines)))
                else
                    error("Unknown UVW coordinate: $lines")
                end
            elseif "f" == command # mesh always has faces
                if any(x-> occursin("//", x), lines)
                    fs = process_face_normal(lines)
                elseif any(x-> occursin("/", x), lines)
                    fs = process_face_uv_or_normal(lines)
                else
                    append!(faces, triangulated_faces(facetype, lines))
                    continue
                end
                for i = 1:length(first(fs))
                    append!(f_uv_n_faces[i], triangulated_faces(facetype, getindex.(fs, i)))
                end
            else
                #TODO
            end
        end
    end
    point_attributes = Dict{Symbol, Any}()
    uv_faces = f_uv_n_faces[2]
    normal_faces = f_uv_n_faces[3]
    if !isempty(v_normals)
        if !isempty(normal_faces)
            normals_remapped = similar(points, eltype(v_normals))
            for (vf, nf) in zip(faces, normal_faces)
                for (vidx, nidx) in zip(vf, nf)
                    normals_remapped[vidx] = v_normals[nidx]
                end
            end
            v_normals = normals_remapped
        else
            # these are not per vertex normals, which we
            # can't deal with at the moment, so we just generate our own!
            if length(points) != length(v_normals)
                v_normals = normals(points, faces, normaltype=normaltype)
            end
        end
        point_attributes[:normals] = v_normals
    end
    if !isempty(uv)
        if !isempty(uv_faces)
            uv_remapped = similar(points, eltype(uv))
            for (vf, uvf) in zip(faces, uv_faces)
                for (vidx, uvidx) in zip(vf, uvf)
                    uv_remapped[vidx] = uv[uvidx]
                end
            end
            uv = uv_remapped
        end
        point_attributes[:uv] = uv
    end

    return Mesh(meta(points; point_attributes...), faces)
end

# of form "faces v1 v2 v3 ....""
process_face(lines::Vector) = (lines,) # just put it in the same format as the others
# of form "faces v1//vn1 v2//vn2 v3//vn3 ..."
process_face_normal(lines::Vector) = split.(lines, "//")
# of form "faces v1/vt1 v2/vt2 v3/vt3 ..." or of form "faces v1/vt1/vn1 v2/vt2/vn2 v3/vt3/vn3 ...."
process_face_uv_or_normal(lines::Vector) = split.(lines, ('/',))

function triangulated_faces(::Type{Tf}, vertex_indices::Vector) where {Tf}
    poly_face = NgonFace{length(vertex_indices), UInt32}(parse.(UInt32, vertex_indices))
    return convert_simplex(Tf, poly_face)
end


text:
# using colorview function to convert an array to image
using Images, ImageView
random_img_array = rand(3, 8, 8); # channel, height, width
img = colorview(RGB, random_img_array);
imshow(img)

# changing order of dimensions prior to conversion
using Images, ImageView
random_img_array = rand(40, 100, 3); # height, width, channel
img_perm = permuteddimsview(random_img_array, (3, 1, 2))
img = colorview(RGB, Float16.(img_perm))
imshow(img)

# the alternative way of writing the same code is as follows:
using Images, ImageView
random_img_array = rand(40, 100, 3); # height, width, channel
img_perm = permutedims(random_img_array, [3, 1, 2])
img = colorview(RGB, img_perm)
imshow(img)


text:
struct UnboundedFluidQueue <: Model
    T::Array{Float64,2}
    S::PhaseSet
    function UnboundedFluidQueue(T::Array{Float64,2}, S::PhaseSet)
        DiscretisedFluidQueues._fluid_queue_checks(T,S)
        !(all(sum(T,dims=2).≈0.0))&&@warn "row sums of T should be 0"
        return new(T,S)
    end
end

"""
    UnboundedFluidQueue(T::Array{Float64,2},c::Array{Float64,1},b::Float64)

Alias to `UnboundedFluidQueue(T,PhaseSet(c))`.
"""
UnboundedFluidQueue(T::Array{Float64,2},c::Array{Float64,1}) = UnboundedFluidQueue(T,PhaseSet(c))

"""
    _model_dicts(model::Model)

    input: a Model object
outputs:
     - SDict: a dictionary with keys `"+","-","0","bullet"`
    and corresponding values `findall(model.C .> 0)`, `findall(model.C .< 0)`,
    `findall(model.C .== 0)`, `findall(model.C .!= 0)`, respectively.
     - TDict: a dictionary of submatrices of `T` with keys
    `"ℓm"` with ``ℓ,m∈{+,-,0,bullet}`` and corresponding values
    `model.T[S[ℓ],S[m]]`.
"""
function _model_dicts(model::Model)
    nPhases = n_phases(model)
    SDict = Dict{String,Array}("S" => 1:nPhases)
    SDict["+"] = findall(rates(model) .> 0)
    SDict["-"] = findall(rates(model) .< 0)
    SDict["0"] = findall(rates(model) .== 0)
    SDict["bullet"] = findall(rates(model) .!= 0)

    TDict = Dict{String,Array}("T" => model.T)
    for ℓ in ["+" "-" "0" "bullet"], m in ["+" "-" "0" "bullet"]
        TDict[ℓ*m] = model.T[SDict[ℓ], SDict[m]]
    end

    return SDict, TDict
end

"""
Construct and evaluate ``Ψ(s)`` for a triditional SFM.

Uses newtons method to solve the Ricatti equation
``D⁺⁻(s) + Ψ(s)D⁻⁺(s)Ψ(s) + Ψ(s)D⁻⁻(s) + D⁺⁺(s)Ψ(s) = 0.``

    psi_fun_x( model::Model; s = 0, MaxIters = 1000, err = 1e-8)

# Arguments
- `model`: a Model object
- `s::Real`: a value to evaluate the LST at
- `MaxIters::Int`: the maximum number of iterations of newtons method
- `err::Float64`: an error tolerance for terminating newtons method. Terminates
    when `max(Ψ_{n} - Ψ{n-1}) .< eps`.

# Output
- `Ψ(s)::Array{Float64,2}` the matrix ``Ψ``
"""
function psi_fun_x( model::Model; s = 0.0, MaxIters = 1000, err = 1e-8)
    SDict, TDict = _model_dicts(model)

    T00inv = inv(TDict["00"] - s * LinearAlgebra.I)
    # construct the generator Q(s)
    Q =
        (1 ./ abs.(rates(model)[SDict["bullet"]])) .* (
            TDict["bulletbullet"] - s * LinearAlgebra.I -
            TDict["bullet0"] * T00inv * TDict["0bullet"]
        )

    # model_without_zero_phases = @suppress_err BoundedFluidQueue(Q,model.S[SDict["bullet"]],model.b)

     ~, QDict = make_QDict(Q,model.S[SDict["bullet"]]) #_model_dicts(model_without_zero_phases)

    Ψ = zeros(Float64, length(SDict["+"]), length(SDict["-"]))
    A = QDict["++"]
    B = QDict["--"]
    D = QDict["+-"]
    # use netwons method to solve the Ricatti equation
    for n in 1:MaxIters
        Ψ = LinearAlgebra.sylvester(A,B,D)
        if maximum(abs.(sum(Ψ,dims=2).-1)) < err
            break
        end
        A = QDict["++"] + Ψ * QDict["-+"]
        B = QDict["--"] + QDict["-+"] * Ψ
        D = QDict["+-"] - Ψ * QDict["-+"] * Ψ
    end

    return Ψ
end

function make_QDict(Q,S_without_zero)
    SDict = Dict{String,Array}("S" => 1:length(S_without_zero))
    SDict["+"] = findall(rates(S_without_zero) .> 0)
    SDict["-"] = findall(rates(S_without_zero) .< 0)
    SDict["0"] = findall(rates(S_without_zero) .== 0)
    SDict["bullet"] = findall(rates(S_without_zero) .!= 0)

    TDict = Dict{String,Array}("T" => Q)
    for ℓ in ["+" "-" "0" "bullet"], m in ["+" "-" "0" "bullet"]
        TDict[ℓ*m] = Q[SDict[ℓ], SDict[m]]
    end

    return SDict, TDict
end

"""
Construct the vector ``ξ`` containing the distribution of the phase at the time
when ``X(t)`` first hits `0`.

    xi_x( model::Model, Ψ::Array)

# Arguments
- `model`: a Model object
- `Ψ`: an array as output from `psi_fun_x`

# Output
- the vector `ξ`
"""
function xi_x( model::Model, Ψ::Array)
    # the system to solve is [ξ 0](-[B₋₋ B₋₀; B₀₋ B₀₀])⁻¹[B₋₊; B₀₊]Ψ = ξ
    # writing this out and using block inversion (as described on wikipedia)
    # we can solve this in the following way
    SDict, TDict = _model_dicts(model)

    T00inv = inv(TDict["00"])
    invT₋₋ =
        inv(TDict["--"] - TDict["-0"] * T00inv * TDict["0-"])
    invT₋₀ = -invT₋₋ * TDict["-0"] * T00inv

    A =
        -(
            invT₋₋ * TDict["-+"] * Ψ + invT₋₀ * TDict["0+"] * Ψ +
            LinearAlgebra.I
        )
    b = zeros(1, size(TDict["--"], 1))
    A[:, 1] .= 1.0 # normalisation conditions
    b[1] = 1.0 # normalisation conditions

    ξ = b / A

    return ξ
end

"""
Construct the stationary distribution of the SFM

    stationary_distribution_x( model::Model, Ψ::Array, ξ::Array)

# Arguments
- `model`: a Model object
- `Ψ`: an array as output from `psi_fun_x`
- `ξ`: an array as returned from `xi_x`

# Output
- `pₓ::Array{Float64,2}`: the point masses of the SFM
- `πₓ(x)` a function with two methods
    - `πₓ(x::Real)`: for scalar inputs, returns the stationary density evaluated
        at `x` in all phases.
    - `πₓ(x::Array)`: for array inputs, returns an array with the same shape
        as is output by Coeff2Dist.
- `K::Array{Float64,2}`: the matrix in the exponential of the density.
"""
function stationary_distribution_x( model::Model, Ψ::Array, ξ::Array)
    # using the same block inversion trick as in xi_x
    SDict, TDict = _model_dicts(model)

    T00inv = inv(TDict["00"])
    invT₋₋ =
        inv(TDict["--"] - TDict["-0"] * T00inv * TDict["0-"])
    invT₋₀ = -invT₋₋ * TDict["-0"] * T00inv

    Q =
        (1 ./ abs.(rates(model)[SDict["bullet"]])) .* (
            TDict["bulletbullet"] -
            TDict["bullet0"] * T00inv * TDict["0bullet"]
        )

    # model_without_zero_phases = @suppress_err BoundedFluidQueue(Q,model.S[SDict["bullet"]],model.b)

    # ~, QDict = _model_dicts(model_without_zero_phases)
     ~, QDict = make_QDict(Q,model.S[SDict["bullet"]])

    K = QDict["++"] + Ψ * QDict["-+"]

    A = -[invT₋₋ invT₋₀]

    # unnormalised values
    αpₓ = ξ * A

    απₓ = αpₓ *
        [TDict["-+"]; TDict["0+"]] *
        -inv(K) *
        [LinearAlgebra.I(length(SDict["+"])) Ψ] *
        LinearAlgebra.diagm(1 ./ abs.(rates(model)[SDict["bullet"]]))

    απₓ0 = -απₓ * [TDict["+0"];TDict["-0"]] * T00inv

    # normalising constant
    α = sum(αpₓ) + sum(απₓ) + sum(απₓ0)

    # normalised values
    # point masses
    pₓ = αpₓ/α
    # density method for scalar x-values
    idx = [findall(rates(model).>0);
        findall(rates(model).<0);
        findall(rates(model).==0)]
    function πₓ(x::Real)
        out = zeros(n_phases(model))
        out[idx] = (pₓ *
        [TDict["-+"]; TDict["0+"]] *
        exp(K*x) *
        [LinearAlgebra.I(length(SDict["+"])) Ψ] *
        LinearAlgebra.diagm(1 ./ abs.(rates(model)[SDict["bullet"]])) *
        [LinearAlgebra.I(sum(rates(model) .!= 0)) [TDict["+0"];TDict["-0"]] * -T00inv])
        return out
    end
    # density method for arrays so that πₓ returns an array with the same shape
    # as is output by Coeff2Dist
    function πₓ(x::Array)
        temp = πₓ.(x)
        Evalπₓ = zeros(Float64, size(x,1), size(x,2), n_phases(model))
        for cell in 1:size(x,2)
            for basis in 1:size(x,1)
                Evalπₓ[basis,cell,:] = temp[basis,cell]
            end
        end
        return Evalπₓ
    end

    # CDF method for scalar x-values
    function Πₓ(x::Real)
        out = zeros(n_phases(model))
        out[idx] = [zeros(1,sum(rates(model).>0)) pₓ] .+
        pₓ *
        [TDict["-+"]; TDict["0+"]] *
        (exp(K*x) - LinearAlgebra.I) / K *
        [LinearAlgebra.I(length(SDict["+"])) Ψ] *
        LinearAlgebra.diagm(1 ./ abs.(rates(model)[SDict["bullet"]])) *
        [LinearAlgebra.I(sum(rates(model) .!= 0)) [TDict["+0"];TDict["-0"]] * -T00inv]
        return out
    end
    # CDF method for arrays so that Πₓ returns an array with the same shape
    # as is output by Coeff2Dist
    function Πₓ(x::Array)
        temp = Πₓ.(x)
        Evalπₓ = zeros(Float64, size(x,1), size(x,2), n_phases(model))
        for cell in 1:size(x,2)
            for basis in 1:size(x,1)
                Evalπₓ[basis,cell,:] = temp[basis,cell]
            end
        end
        return Evalπₓ
    end

    return pₓ, πₓ, Πₓ, K
end

function stationary_distribution_x( model::Model)
    Ψ = psi_fun_x( model)
    ξ = xi_x( model, Ψ)
    pₓ, πₓ, Πₓ, K = stationary_distribution_x(model,Ψ,ξ)
    return pₓ, πₓ, Πₓ, K
end
rate_reverse(model::BoundedFluidQueue) = BoundedFluidQueue(model.T,-rates(model),model.P_lwr,model.P_upr,model.b)

function expected_crossings(model::BoundedFluidQueue)
    rate_reverse_model = rate_reverse(model)
    ~, TDict = _model_dicts(model)
    ~, T̃Dict = _model_dicts(rate_reverse_model)
    Ψ = psi_fun_x(model)
    Ψ̃ = psi_fun_x(rate_reverse_model)
    K = TDict["++"] + Ψ*TDict["-+"]
    K̃ = T̃Dict["++"] + Ψ̃*T̃Dict["-+"]
    n₊ = size(TDict["++"],1)
    n₋ = size(TDict["--"],1)
    I₊ = LinearAlgebra.I(n₊)
    I₋ = LinearAlgebra.I(n₋)
    N(x) = [I₊ exp(K*model.b)*Ψ; exp(K̃*model.b)*Ψ̃ I₋]\[exp(K*x) zeros(n₊,n₋); zeros(n₋,n₊) exp(K̃*(model.b-x))]*
        [I₊ Ψ; Ψ̃ I₋]
    return N
end

function int_expected_crossings(model::BoundedFluidQueue)
    rate_reverse_model = rate_reverse(model)
    ~, TDict = _model_dicts(model)
    ~, T̃Dict = _model_dicts(rate_reverse_model)
    Ψ = psi_fun_x(model)
    Ψ̃ = psi_fun_x(rate_reverse_model)
    K = TDict["++"] + Ψ*TDict["-+"]
    K̃ = T̃Dict["++"] + Ψ̃*T̃Dict["-+"]
    n₊ = size(TDict["++"],1)
    n₋ = size(TDict["--"],1)
    I₊ = LinearAlgebra.I(n₊)
    I₋ = LinearAlgebra.I(n₋)
    N(x) = [I₊ exp(K*model.b)*Ψ; exp(K̃*model.b)*Ψ̃ I₋]\[K\(exp(K*x)-I₊) zeros(n₊,n₋); zeros(n₋,n₊) K̃\(exp(K̃*model.b)-exp(K̃*(model.b-x)))]*
        [I₊ Ψ; Ψ̃ I₋]
    return N
end

function hat_matrices(model::BoundedFluidQueue)
    rate_reverse_model = rate_reverse(model)
    ~, TDict = _model_dicts(model)
    # ~, T̃Dict = _model_dicts(rate_reverse_model)
    Ψ = psi_fun_x(model)
    Ψ̃ = psi_fun_x(rate_reverse_model)
    U = TDict["--"] + TDict["-+"]*Ψ
    Ũ = TDict["++"] + TDict["+-"]*Ψ̃
    n₊ = size(TDict["++"],1)
    n₋ = size(TDict["--"],1)
    I₊ = LinearAlgebra.I(n₊)
    I₋ = LinearAlgebra.I(n₋)
    b = model.b
    Λ̂₁₁ = (I₊ - Ψ*Ψ̃)*exp(Ũ*b)/(I₊ - Ψ*exp(U*b)*Ψ̃*exp(Ũ*b))
    Ψ̂₁₂ = (Ψ - exp(Ũ*b)*Ψ*exp(U*b))/(I₋ - Ψ̃*exp(Ũ*b)*Ψ*exp(U*b))
    Ψ̂₂₁ = (Ψ̃ - exp(U*b)*Ψ̃*exp(Ũ*b))/(I₊ - Ψ*exp(U*b)*Ψ̃*exp(Ũ*b))
    Λ̂₂₂ = (I₋ - Ψ̃*Ψ)*exp(U*b)/(I₋ - Ψ̃*exp(Ũ*b)*Ψ*exp(U*b))
    return Λ̂₁₁, Ψ̂₁₂, Ψ̂₂₁, Λ̂₂₂
end

function H_matrix(model)
    Λ̂₁₁, Ψ̂₁₂, Ψ̂₂₁, Λ̂₂₂ = hat_matrices(model)
    hat_matrix = [Λ̂₁₁ Ψ̂₁₂; Ψ̂₂₁ Λ̂₂₂]
    S, TDict = _model_dicts(model)
    n₊ = size(TDict["++"],1)
    n₋ = size(TDict["--"],1)
    n = length(model.S)
    Ĥ = hat_matrix*[zeros(n₊,n₊) model.P_upr[:,S["-"]]; model.P_lwr[:,S["+"]] zeros(n₋,n₋)] +
        hat_matrix*[model.P_upr[:,[i∉S["-"] for i in 1:n]] zeros(n₊,n-n₊); zeros(n₋,n-n₋) model.P_lwr[:,[i∉S["+"] for i in 1:n]]]*
            [-inv(model.T[[i∉S["-"] for i in 1:n],[i∉S["-"] for i in 1:n]]) zeros(n-n₋,n-n₊);
             zeros(n-n₊,n-n₋) -inv(model.T[[i∉S["+"] for i in 1:n],[i∉S["+"] for i in 1:n]])]*
                [zeros(n-n₋,n₊) model.T[[i∉S["-"] for i in 1:n],S["-"]];
                model.T[[i∉S["+"] for i in 1:n],S["+"]] zeros(n-n₊,n₋)]
    return Ĥ
end

function nu(model)
    A = H_matrix(model) - LinearAlgebra.I
    A[:,1] .= 1.0
    b = zeros(1,size(A,1))
    b[1] = 1.0
    return b/A
end

function normalising_constant_matrix(model)
    rate_reverse_model = rate_reverse(model)
    ~, TDict = _model_dicts(model)
    ~, T̃Dict = _model_dicts(rate_reverse_model)
    Ψ = psi_fun_x(model)
    Ψ̃ = psi_fun_x(rate_reverse_model)
    K = TDict["++"] + Ψ*TDict["-+"]
    K̃ = T̃Dict["++"] + Ψ̃*T̃Dict["-+"]
    n₊ = size(TDict["++"],1)
    n₋ = size(TDict["--"],1)
    I₊ = LinearAlgebra.I(n₊)
    I₋ = LinearAlgebra.I(n₋)
    if abs(LinearAlgebra.det(K))>sqrt(eps()) && abs(LinearAlgebra.det(K̃))>sqrt(eps())
        N = [I₊ exp(K*model.b)*Ψ; exp(K̃*model.b)*Ψ̃ I₋]\[Matrix(K)\(exp(K*model.b)-I₊) zeros(n₊,n₋); zeros(n₋,n₊) Matrix(K̃)\(exp(K̃*model.b)-I₋)]*
            [I₊ Ψ; Ψ̃ I₋]
    elseif abs(LinearAlgebra.det(K))>sqrt(eps())
        N = Matrix([I₊ exp(K*model.b)*Ψ; exp(K̃*model.b)*Ψ̃ I₋])\Matrix([Matrix(K)\(exp(K*model.b)-I₊) zeros(n₊,n₋); zeros(n₋,n₊) model.b*I₋])*
            [I₊ Ψ; Ψ̃ I₋]
    elseif abs(LinearAlgebra.det(K̃))>sqrt(eps())
        N = [I₊ exp(K*model.b)*Ψ; exp(K̃*model.b)*Ψ̃ I₋]\[model.b*I₊ zeros(n₊,n₋); zeros(n₋,n₊) Matrix(K̃)\(exp(K̃*model.b)-I₋)]*
            [I₊ Ψ; Ψ̃ I₋]
    else
        N = [I₊ exp(K*model.b)*Ψ; exp(K̃*model.b)*Ψ̃ I₋]\[model.b*I₊ zeros(n₊,n₋); zeros(n₋,n₊) model.b*I₋]*
            [I₊ Ψ; Ψ̃ I₋]
    end
    return N
end

function B_matrices(model::BoundedFluidQueue)
    Λ̂₁₁, Ψ̂₁₂, Ψ̂₂₁, Λ̂₂₂ = hat_matrices(model)
    S, ~ = _model_dicts(model)
    n = n_phases(model)
    B₁ₒ = (LinearAlgebra.I - Ψ̂₁₂*model.P_lwr[:,S["+"]])\Ψ̂₁₂*model.P_lwr[:,[i∉S["+"] for i in 1:n]]
    B₂ₒ = (LinearAlgebra.I - Ψ̂₂₁*model.P_upr[:,S["-"]])\Ψ̂₂₁*model.P_upr[:,[i∉S["-"] for i in 1:n]]
    B₁₁ = (LinearAlgebra.I - Ψ̂₁₂*model.P_lwr[:,S["+"]])\Λ̂₁₁
    A₂₂ = (LinearAlgebra.I - Ψ̂₂₁*model.P_upr[:,S["-"]])\Λ̂₂₂
    return B₁ₒ, B₂ₒ, B₁₁, A₂₂
end

function omega_matrices(model::BoundedFluidQueue)
    B₁ₒ, B₂ₒ, B₁₁, A₂₂ = B_matrices(model)
    n = n_phases(model)
    S, ~ = _model_dicts(model)
    Ω̂₁ₒ = (LinearAlgebra.I - B₁₁*model.P_upr[:,S["-"]]*A₂₂*model.P_lwr[:,S["+"]])\
        (B₁₁*model.P_upr[:,[i∉S["-"] for i in 1:n]] + B₁₁*model.P_upr[:,S["-"]]*B₂ₒ)
    Ω₂ₒ = (LinearAlgebra.I - A₂₂*model.P_lwr[:,S["+"]]*B₁₁*model.P_upr[:,S["-"]])\
        (A₂₂*model.P_lwr[:,[i∉S["+"] for i in 1:n]] + A₂₂*model.P_lwr[:,S["+"]]*B₁ₒ)
    Ω₁ₒ = (LinearAlgebra.I - B₁₁*model.P_upr[:,S["-"]]*A₂₂*model.P_lwr[:,S["+"]])\
        (B₁₁*model.P_upr[:,S["-"]]*A₂₂*model.P_lwr[:,S["+"]]*B₁ₒ +B₁₁*model.P_upr[:,S["-"]]*A₂₂*model.P_lwr[:,[i∉S["+"] for i in 1:n]])
    Ω̂₂ₒ = (LinearAlgebra.I - A₂₂*model.P_lwr[:,S["+"]]*B₁₁*model.P_upr[:,S["-"]])\
        (A₂₂*model.P_lwr[:,S["+"]]*B₁₁*model.P_upr[:,S["-"]]*B₂ₒ + A₂₂*model.P_lwr[:,S["+"]]*B₁₁*model.P_upr[:,[i∉S["-"] for i in 1:n]])
    return Ω̂₁ₒ, Ω₂ₒ, Ω₁ₒ, Ω̂₂ₒ
end

function point_mass(model)
    B₁ₒ, B₂ₒ, B₁₁, A₂₂ = B_matrices(model)
    Ω̂₁ₒ, Ω₂ₒ, Ω₁ₒ, Ω̂₂ₒ = omega_matrices(model)
    n = n_phases(model)
    S, ~ = _model_dicts(model)
    Z₁₁ = model.T[[i∉S["-"] for i in 1:n],[i∉S["-"] for i in 1:n]] + model.T[[i∉S["-"] for i in 1:n],S["-"]]*(B₂ₒ+Ω̂₂ₒ)
    Z₁₂ = model.T[[i∉S["-"] for i in 1:n],S["-"]]*Ω₂ₒ
    Z₂₁ = model.T[[i∉S["+"] for i in 1:n],S["+"]]*Ω̂₁ₒ
    Z₂₂ = model.T[[i∉S["+"] for i in 1:n],[i∉S["+"] for i in 1:n]] + model.T[[i∉S["+"] for i in 1:n],S["+"]]*(B₁ₒ+Ω₁ₒ)
    Z = [Z₁₁ Z₁₂; Z₂₁ Z₂₂]
    Z[:,1] .= 1.0
    b = zeros(1,size(Z,1))
    b[1] = 1.0
    return b/Z
end

function preprocess(model::Model)
    S,T = _model_dicts(model)
    not_zero = [i∉S["0"] for i in 1:n_phases(model)]
    C = rates(model)[not_zero]
    Cinv = 1.0./abs.(C)
    Q = Cinv.*model.T[not_zero,not_zero] +
        Cinv.*model.T[not_zero,S["0"]]*(-model.T[S["0"],S["0"]])^-1*model.T[S["0"],not_zero]
    P_lwr = model.P_lwr[:,not_zero]
    P_upr = model.P_lwr[:,not_zero]
    return BoundedFluidQueue(Q,C,P_lwr,P_upr,model.b)
end

"""
This might be the least robust piece of code ever written. Use with caution. Answers you get might
just be straight up wrong without warning. Be careful.
"""
function stationary_distribution_x(model::BoundedFluidQueue)
    ν = nu(model)
    S, T = _model_dicts(model)
    Cinv = LinearAlgebra.diagm(1.0./[rates(model)[S["+"]];abs.(rates(model)[S["-"]])])
    n = n_phases(model)

    # density matrix
    N = expected_crossings(model)

    # cdf matrix
    int_N = int_expected_crossings(model)

    # boundary probabilities
    p = point_mass(model)
    n_upr = n-length(S["-"])
    p_upr = p[1:n_upr]'
    p_lwr = p[n_upr+1:end]'
    # p_to_density = [p_lwr*model.T[[i∉S["+"] for i in 1:n],S["+"]] p_upr*model.T[[i∉S["-"] for i in 1:n],S["-"]]]
    k = sum(p)

    # normalising constant
    # integral over denisty matrix
    N_const = int_N(model.b)# DiscretisedFluidQueues.normalising_constant_matrix(model)
    # integrals over densities for + and -
    plus_minus_m_norm_vec = ν*N_const*Cinv
    # intetegral over zero phases
    zero_norm_vec = plus_minus_m_norm_vec*[T["+0"];T["-0"]]*(-inv(T["00"]))
    # the constant
    z = sum(plus_minus_m_norm_vec) + sum(zero_norm_vec) #+ sum(p)

    if (sum(p_lwr)==0.0)&&(sum(p_upr)==0.0)
        cz = z
        F0 = zeros(n)
        Fb = zeros(n)
    elseif (sum(p_upr)==0.0)
        vec0₋ = (ν*N(0)*Cinv)[length(S["+"])+1:end]
        π0₋ = vec0₋.*model.P_lwr[:,[i∉S["+"] for i in 1:n]]
        z0 = sum(π0₋)

        k0 = sum(p_lwr*model.T[[i∉S["+"] for i in 1:n],S["+"]])

        r0 = z0/k0

        cz = z + k*r0
        cp = cz/r0

        p /= cp
        p_upr = p[1:n_upr]'
        p_lwr = p[n_upr+1:end]'

        F0 = zeros(n)
        Fb = zeros(n)
        counter0 = 0
        counterb = 0
        for i in 1:n
            if rates(model)[i]<=0
                counter0 += 1
                F0[i] = p_lwr[counter0]
            end
            if rates(model)[i]>=0
                counterb += 1
                Fb[i] = p_upr[counterb]
            end
        end
    else
        vec0₊ = (ν*N(model.b)*Cinv)[1:length(S["+"])]
        π0₊ = vec0₊.*model.P_upr[:,[i∉S["-"] for i in 1:n]]
        z0 = sum(π0₊)

        k0 = sum(p_upr*model.T[[i∉S["-"] for i in 1:n],S["-"]])

        r0 = z0/k0

        cz = z + k*r0
        cp = cz/r0

        p /= cp
        p_upr = p[1:n_upr]'
        p_lwr = p[n_upr+1:end]'

        F0 = zeros(n)
        Fb = zeros(n)
        counter0 = 0
        counterb = 0
        for i in 1:n
            if rates(model)[i]<=0
                counter0 += 1
                F0[i] = p_lwr[counter0]
            end
            if rates(model)[i]>=0
                counterb += 1
                Fb[i] = p_upr[counterb]
            end
        end
    end
    # p /= z
    # p_upr = p[1:n_upr]'
    # p_lwr = p[n_upr+1:end]'
    # p_to_density = [p_lwr*model.T[[i∉S["+"] for i in 1:n],S["+"]] p_upr*model.T[[i∉S["-"] for i in 1:n],S["-"]]]

    idx = [findall(rates(model).>0);
        findall(rates(model).<0);
        findall(rates(model).==0)]
    function πₓ(x)
        vec = ν*N(x)*Cinv
        ϕ̂ = vec*[T["+0"];T["-0"]]*(-inv(T["00"]))
        out = [vec ϕ̂]./cz
        return out[idx]
    end
    function Πₓ(x)
        vec = ν*int_N(x)*Cinv
        ϕ̂ = vec*[T["+0"];T["-0"]]*(-inv(T["00"]))
        out = [vec ϕ̂]./cz
        return (F0 + out[idx])*(zero(x)<=x) + Fb*(model.b<=x)
    end
    return πₓ, Πₓ
end

## HITTING TIMES ##
_error_on_neg(x) = (x>=0.0) ? true : throw(DomainError("x must be positive"))
"""
Returns G₋₋(s,x), G₊₋(s,x)
"""
function draining_times_lst(model::Model)
    Ψ(s::ComplexF64) = psi_fun_x(model,s=s)
    S,T = _model_dicts(model)
    B(s::ComplexF64) = (-1.0./rates(model)[S["-"]])*(T["--"]-s*LinearAlgebra.I + T["-+"]*Ψ(s) - T["-0"]*inv(T["00"]-s*LinearAlgebra.I)*(T["0-"] + T["0+"]*Ψ(s)))
    G₋₋(s::ComplexF64,x::Float64) = _error_on_neg(x)&&exp(B(s)*x)
    G₊₋(s::ComplexF64,x::Float64) = Ψ(s)*G₋₋(s,x)
    G₋₋(s::Real,x::Real) = G₋₋(Complex{Float64}(s),Float64(x))
    G₊₋(s::Real,x::Real) = G₊₋(Complex{Float64}(s),Float64(x))
    return G₋₋, G₊₋
end
function draining_times_lst(model::Model,s,x)
    _error_on_neg(x)
    Ψ = psi_fun_x(model,s=s)
    S,T = _model_dicts(model)
    B = (-1.0./rates(model)[S["-"]])*(T["--"]-s*LinearAlgebra.I + T["-+"]*Ψ - T["-0"]*inv(T["00"]-s*LinearAlgebra.I)*(T["0-"] + T["0+"]*Ψ))
    G₋₋ = exp(B*x)
    G₊₋ = Ψ*G₋₋
    return G₋₋, G₊₋
end
"""
Returns H₊₊(s,x), H₋₊(s,x)
"""
filling_times_lst(model::Model) = draining_times_lst(rate_reverse(model))
filling_times_lst(model::Model,s,x) = draining_times_lst(rate_reverse(model),s,x)

function draining_filling_times_lst(model::Model,s,x,y)
    _error_on_neg(x)
    _error_on_neg(y)
    _error_on_neg(y-x)
    Ψ = psi_fun_x(model,s=s)
    S,T = _model_dicts(model)
    the_inv = inv(T["00"]-s*LinearAlgebra.I)
    B = (-1.0./rates(model)[S["-"]])*(T["--"]-s*LinearAlgebra.I + T["-+"]*Ψ - T["-0"]*the_inv*(T["0-"] + T["0+"]*Ψ))
    G₋₋y = exp(B*y)
    G₊₋y = Ψ*G₋₋y
    G₋₋x = exp(B*x)
    G₊₋x = Ψ*G₋₋x

    Ξ = psi_fun_x(rate_reverse(model),s=s)
    # S,T = _model_dicts(model)
    A = (1.0./rates(model)[S["+"]])*(T["++"]-s*LinearAlgebra.I + T["+-"]*Ξ - T["+0"]*the_inv*(T["0+"] + T["0-"]*Ξ))
    H₊₊y = exp(A*y)
    H₋₊y = Ξ*H₊₊y
    H₊₊z = exp(A*(y-x))
    H₋₊z = Ξ*H₊₊z

    n₊, n₋ = size(G₊₋y)

    Gy = [zeros(n₊,n₊) G₊₋y; zeros(n₋,n₊) G₋₋y]
    Gx = [zeros(n₊,n₊) G₊₋x; zeros(n₋,n₊) G₋₋x]
    Hy= [H₊₊y zeros(n₊,n₋); H₋₊y zeros(n₋,n₋)]
    Hz= [H₊₊z zeros(n₊,n₋); H₋₊z zeros(n₋,n₋)]
    lst = [Gx Hz]/[LinearAlgebra.I Hy; Gy LinearAlgebra.I]
    return lst
end

function build_H_lst(model::Model)
    H₊₊, H₋₊ = filling_times_lst(model)
    S, ~ = _model_dicts(model)
    n₋, n₊ = length(S["-"]), length(S["+"])
    H(s::ComplexF64,x::Float64) = _error_on_neg(x)&&[H₊₊(s,x) zeros(n₊,n₋); H₋₊(s,x) zeros(n₋,n₋)]
    H(s::Real,x::Real) = H(Complex{Float64}(s),Float64(x))
    return H
end
function build_H_lst(model::Model,s,x)
    H₊₊, H₋₊ = filling_times_lst(model,s,x)
    S, ~ = _model_dicts(model)
    n₋, n₊ = length(S["-"]), length(S["+"])
    H = [H₊₊ zeros(n₊,n₋); H₋₊ zeros(n₋,n₋)]
    return H
end

function build_G_lst(model::Model)
    G₋₋, G₊₋ = draining_times_lst(model)
    S, ~ = _model_dicts(model)
    n₋, n₊ = length(S["-"]), length(S["+"])
    G(s::ComplexF64,x::Float64) = [zeros(n₊,n₊) G₊₋(s,x); zeros(n₋,n₊) G₋₋(s,x)]
    G(s::Real,x::Real) = G(Complex{Float64}(s),Float64(x))
    return G
end
function build_G_lst(model::Model,s,x)
    G₋₋, G₊₋ = draining_times_lst(model,s,x)
    S, ~ = _model_dicts(model)
    n₋, n₊ = length(S["-"]), length(S["+"])
    G = [zeros(n₊,n₊) G₊₋; zeros(n₋,n₊) G₋₋]
    return G
end

function hitting_times_lst(model::Model)
    G = build_G_lst(model)
    H = build_H_lst(model)
    lst(s::ComplexF64,x::Float64,y::Float64) = [G(s,x) H(s,y-x)]/[LinearAlgebra.I H(s,y); G(s,y) LinearAlgebra.I]
    lst(s::Real,x::Real,y::Real) = lst(Complex{Float64}(s),Float64(x),Float64(y))
    return lst
end
function hitting_times_lst(model::Model,s,x,y)
    Gx = build_G_lst(model,s,x)
    Gy = build_G_lst(model,s,y)
    Hz = build_H_lst(model,s,y-x)
    Hy = build_H_lst(model,s,y)
    lst = [Gx Hz]/[LinearAlgebra.I Hy; Gy LinearAlgebra.I]
    return lst
end

for fun_name in (:filling, :draining)
    @eval function $(Symbol(fun_name,"_times_fun"))(model::Model,t,x)
        fun = ilt(cme,s->$(Symbol(fun_name,"_times_lst"))(model,s,x),t,n)
        return fun
    end
end

function hitting_times_cdf(model,t,x,y,n=21)
    return ilt(cme,s->draining_filling_times_lst(model,s,x,y)./s,t,n)
end


text:
@testset "Access and Iterations" begin
    dna_kmer = mer"ACTG"dna
    rna_kmer = mer"ACUG"rna
    aa_kmer  = mer"MVXN"aa

    @testset "Access DNA Kmer" begin
        @test dna_kmer[1] == DNA_A
        @test dna_kmer[2] == DNA_C
        @test dna_kmer[3] == DNA_T
        @test dna_kmer[4] == DNA_G

        # Access indexes out of bounds
        @test_throws BoundsError dna_kmer[-1]
        @test_throws BoundsError dna_kmer[0]
        @test_throws BoundsError dna_kmer[5]
        @test_throws BoundsError getindex(dna_kmer,-1)
        @test_throws BoundsError getindex(dna_kmer, 0)
        @test_throws BoundsError getindex(dna_kmer, 5)
    end

    @testset "Iteration through DNA Kmer" begin
        @test iterate(DNAKmer("ACTG")) == (DNA_A, 2)

        @test iterate(DNAKmer("ACTG"), 1) == (DNA_A, 2)
        @test iterate(DNAKmer("ACTG"), 4) == (DNA_G, 5)

        @test iterate(DNAKmer("ACTG"), 1)  !== nothing
        @test iterate(DNAKmer("ACTG"), 4)  !== nothing
        @test iterate(DNAKmer("ACTG"), 5)  === nothing
        @test_throws BoundsError iterate(DNAKmer("ACTG"), -1)

        dna_vec = [DNA_A, DNA_C, DNA_T, DNA_G]
        @test all([nt === dna_vec[i] for (i, nt) in enumerate(dna_kmer)])
    end

    @testset "Access RNA Kmer" begin
        @test rna_kmer[1] == RNA_A
        @test rna_kmer[2] == RNA_C
        @test rna_kmer[3] == RNA_U
        @test rna_kmer[4] == RNA_G

        # Access indexes out of bounds
        @test_throws BoundsError rna_kmer[-1]
        @test_throws BoundsError rna_kmer[0]
        @test_throws BoundsError rna_kmer[5]
        @test_throws BoundsError getindex(rna_kmer, -1)
        @test_throws BoundsError getindex(rna_kmer, 0)
        @test_throws BoundsError getindex(rna_kmer, 5)
    end

    @testset "Iteration through RNA Kmer" begin
        @test iterate(RNAKmer("ACUG")) == (RNA_A, 2)

        @test iterate(RNAKmer("ACUG"), 1) == (RNA_A, 2)
        @test iterate(RNAKmer("ACUG"), 4) == (RNA_G, 5)

        @test iterate(RNAKmer("ACUG"), 1)  !== nothing
        @test iterate(RNAKmer("ACUG"), 4)  !== nothing
        @test iterate(RNAKmer("ACUG"), 5)  === nothing
        @test_throws BoundsError iterate(RNAKmer("ACUG"), -1)

        rna_vec = [RNA_A, RNA_C, RNA_U, RNA_G]
        @test all([nt === rna_vec[i] for (i, nt) in enumerate(rna_kmer)])
    end

    @testset "Access AA Kmer" begin
        @test aa_kmer[1] == AA_M
        @test aa_kmer[2] == AA_V
        @test aa_kmer[3] == AA_X
        @test aa_kmer[4] == AA_N

        # Access indexes out of bounds
        @test_throws BoundsError aa_kmer[-1]
        @test_throws BoundsError aa_kmer[0]
        @test_throws BoundsError aa_kmer[5]
        @test_throws BoundsError getindex(aa_kmer,-1)
        @test_throws BoundsError getindex(aa_kmer, 0)
        @test_throws BoundsError getindex(aa_kmer, 5)
    end
end


text:
module RegistryIndexer
using Base: UUID, SHA1
using TOML
using Tar
using LazilyInitializedFields

export registry_info

include("Versions.jl")
# The content of a registry is assumed to be constant during the
# lifetime of a `Registry`. Create a new `Registry` if you want to have
# a new view on the current registry.

function to_tar_path_format(file::AbstractString)
    @static if Sys.iswindows()
        file = replace(file, "\\" => "/")
    end
    return file
end

# See loading.jl
const TOML_CACHE = Base.TOMLCache(TOML.Parser(), Dict{String,Dict{String,Any}}())
const TOML_LOCK = ReentrantLock()
_parsefile(toml_file::AbstractString) = Base.parsed_toml(toml_file, TOML_CACHE, TOML_LOCK)
function parsefile(in_memory_registry::Union{Dict,Nothing}, folder::AbstractString, file::AbstractString)
    if in_memory_registry === nothing
        return _parsefile(joinpath(folder, file))
    else
        content = in_memory_registry[to_tar_path_format(file)]
        return TOML.Internals.parse(TOML.Parser(content; filepath=file))
    end
end

custom_isfile(in_memory_registry::Union{Dict,Nothing}, folder::AbstractString, file::AbstractString) =
    in_memory_registry === nothing ? isfile(joinpath(folder, file)) : haskey(in_memory_registry, to_tar_path_format(file))

# Info about each version of a package
@lazy mutable struct VersionInfo
    git_tree_sha1::Base.SHA1
    yanked::Bool
    @lazy uncompressed_compat::Union{Dict{UUID,Versions.VersionSpec}}
end
VersionInfo(git_tree_sha1::Base.SHA1, yanked::Bool) = VersionInfo(git_tree_sha1, yanked, uninit)

# This is the information that exists in e.g. General/A/ACME
struct PkgInfo
    # Package.toml:
    repo::Union{String,Nothing}
    subdir::Union{String,Nothing}

    # Versions.toml:
    version_info::Dict{VersionNumber,VersionInfo}

    # Compat.toml
    compat::Dict{Versions.VersionRange,Dict{String,Versions.VersionSpec}}

    # Deps.toml
    deps::Dict{Versions.VersionRange,Dict{String,UUID}}
end

isyanked(pkg::PkgInfo, v::VersionNumber) = pkg.version_info[v].yanked
treehash(pkg::PkgInfo, v::VersionNumber) = pkg.version_info[v].git_tree_sha1

function uncompress(compressed::Dict{Versions.VersionRange,Dict{String,T}}, vsorted::Vector{VersionNumber}) where {T}
    @assert issorted(vsorted)
    uncompressed = Dict{VersionNumber,Dict{String,T}}()
    for v in vsorted
        uncompressed[v] = Dict{String,T}()
    end
    for (vs, data) in compressed
        first = length(vsorted) + 1
        # We find the first and last version that are in the range
        # and since the versions are sorted, all versions in between are sorted
        for i in eachindex(vsorted)
            v = vsorted[i]
            v in vs && (first = i; break)
        end
        last = 0
        for i in reverse(eachindex(vsorted))
            v = vsorted[i]
            v in vs && (last = i; break)
        end
        for i in first:last
            v = vsorted[i]
            uv = uncompressed[v]
            for (key, value) in data
                if haskey(uv, key)
                    # Change to an error?
                    error("Overlapping ranges for $(key) in $(repr(path)) for version $v.")
                else
                    uv[key] = value
                end
            end
        end
    end
    return uncompressed
end

const JULIA_UUID = UUID("1222c4b2-2114-5bfd-aeef-88e4692bbb3e")
function initialize_uncompressed!(pkg::PkgInfo, versions=keys(pkg.version_info))
    # Only valid to call this with existing versions of the package
    # Remove all versions we have already uncompressed
    versions = filter!(v -> !isinit(pkg.version_info[v], :uncompressed_compat), collect(versions))

    sort!(versions)

    uncompressed_compat = uncompress(pkg.compat, versions)
    uncompressed_deps   = uncompress(pkg.deps,   versions)

    for v in versions
        vinfo = pkg.version_info[v]
        compat = Dict{UUID,Versions.VersionSpec}()
        uncompressed_deps_v = uncompressed_deps[v]
        # Everything depends on Julia
        uncompressed_deps_v["julia"] = JULIA_UUID
        uncompressed_compat_v = uncompressed_compat[v]
        for (pkg, uuid) in uncompressed_deps_v
            vspec = get(uncompressed_compat_v, pkg, nothing)
            compat[uuid] = vspec === nothing ? Versions.VersionSpec() : vspec
        end
        @init! vinfo.uncompressed_compat = compat
    end
    return pkg
end

function compat_info(pkg::PkgInfo)
    initialize_uncompressed!(pkg)
    return Dict(v => info.uncompressed_compat for (v, info) in pkg.version_info)
end

@lazy struct PkgEntry
    # Registry.toml:
    path::String
    registry_path::String
    name::String
    uuid::UUID

    in_memory_registry::Union{Dict{String,String},Nothing}
    # Version.toml / (Compat.toml / Deps.toml):
    @lazy info::PkgInfo
end

registry_info(pkg::PkgEntry) = init_package_info!(pkg)

function init_package_info!(pkg::PkgEntry)
    # Already uncompressed the info for this package, return early
    @isinit(pkg.info) && return pkg.info
    path = pkg.registry_path

    d_p = parsefile(pkg.in_memory_registry, pkg.registry_path, joinpath(pkg.path, "Package.toml"))
    name = d_p["name"]::String
    name != pkg.name && error("inconsistend name in Registry.toml and Package.toml for pkg at $(path)")
    repo = get(d_p, "repo", nothing)::Union{Nothing,String}
    subdir = get(d_p, "subdir", nothing)::Union{Nothing,String}

    # Versions.toml
    d_v = custom_isfile(pkg.in_memory_registry, pkg.registry_path, joinpath(pkg.path, "Versions.toml")) ?
        parsefile(pkg.in_memory_registry, pkg.registry_path, joinpath(pkg.path, "Versions.toml")) : Dict{String,Any}()
    version_info = Dict{VersionNumber,VersionInfo}(VersionNumber(k) =>
        VersionInfo(SHA1(v["git-tree-sha1"]::String), get(v, "yanked", false)::Bool) for (k, v) in d_v)

    # Compat.toml
    compat_data_toml = custom_isfile(pkg.in_memory_registry, pkg.registry_path, joinpath(pkg.path, "Compat.toml")) ?
        parsefile(pkg.in_memory_registry, pkg.registry_path, joinpath(pkg.path, "Compat.toml")) : Dict{String,Any}()
    # The Compat.toml file might have string or vector values
    compat_data_toml = convert(Dict{String,Dict{String,Union{String,Vector{String}}}}, compat_data_toml)
    compat = Dict{Versions.VersionRange,Dict{String,Versions.VersionSpec}}()
    for (v, data) in compat_data_toml
        vr = Versions.VersionRange(v)
        d = Dict{String,Versions.VersionSpec}(dep => Versions.VersionSpec(vr_dep) for (dep, vr_dep) in data)
        compat[vr] = d
    end

    # Deps.toml
    deps_data_toml = custom_isfile(pkg.in_memory_registry, pkg.registry_path, joinpath(pkg.path, "Deps.toml")) ?
        parsefile(pkg.in_memory_registry, pkg.registry_path, joinpath(pkg.path, "Deps.toml")) : Dict{String,Any}()
    # But the Deps.toml only have strings as values
    deps_data_toml = convert(Dict{String,Dict{String,String}}, deps_data_toml)
    deps = Dict{Versions.VersionRange,Dict{String,UUID}}()
    for (v, data) in deps_data_toml
        vr = Versions.VersionRange(v)
        d = Dict{String,UUID}(dep => UUID(uuid) for (dep, uuid) in data)
        deps[vr] = d
    end
    # All packages depend on julia
    deps[Versions.VersionRange()] = Dict("julia" => JULIA_UUID)

    @init! pkg.info = PkgInfo(repo, subdir, version_info, compat, deps)

    return pkg.info
end

function uncompress_registry(tar_gz::AbstractString)
                if !isfile(tar_gz)
        error("$(repr(tar_gz)): No such file")
    end
    data = Dict{String,String}()
    buf = Vector{UInt8}(undef, Tar.DEFAULT_BUFFER_SIZE)
    io = IOBuffer()
    open(`$(exe7z()) x $tar_gz -so`) do tar
        Tar.read_tarball(x -> true, tar; buf=buf) do hdr, _
            if hdr.type == :file
                Tar.read_data(tar, io; size=hdr.size, buf=buf)
                data[hdr.path] = String(take!(io))
            end
        end
    end
    return data
end

struct RegistryInstance
    path::String
    name::String
    uuid::UUID
    repo::Union{String,Nothing}
    description::Union{String,Nothing}
    pkgs::Dict{UUID,PkgEntry}
    tree_info::Union{Base.SHA1,Nothing}
    in_memory_registry::Union{Nothing,Dict{String,String}}
    # various caches
    name_to_uuids::Dict{String,Vector{UUID}}
end

const REGISTRY_CACHE = Dict{String,Tuple{Base.SHA1,Bool,RegistryInstance}}()

function get_cached_registry(path, tree_info::Base.SHA1, compressed::Bool)
    if !ispath(path)
        delete!(REGISTRY_CACHE, path)
        return nothing
    end
    v = get(REGISTRY_CACHE, path, nothing)
    if v !== nothing
        cached_tree_info, cached_compressed, reg = v
        if cached_tree_info == tree_info && cached_compressed == compressed
            return reg
        end
    end
    # Prevent hogging up memory indefinitely
    length(REGISTRY_CACHE) > 20 && empty!(REGISTRY_CACHE)
    return nothing
end

function RegistryInstance(path::AbstractString)
    compressed_file = nothing
    if isfile(path)
        @assert splitext(path)[2] == ".toml"
        d_reg_info = parsefile(nothing, dirname(path), basename(path))
        compressed_file = d_reg_info["path"]::String
        tree_info = Base.SHA1(d_reg_info["git-tree-sha1"]::String)
    else
        tree_info_file = joinpath(path, ".tree_info.toml")
        tree_info = if isfile(tree_info_file)
            Base.SHA1(parsefile(nothing, path, ".tree_info.toml")["git-tree-sha1"]::String)
        else
            nothing
        end
    end
    # Reuse an existing cached registry if it exists for this content
    if tree_info !== nothing
        reg = get_cached_registry(path, tree_info, compressed_file !== nothing)
        if reg isa RegistryInstance
            return reg
        end
    end

    in_memory_registry = if compressed_file !== nothing
        uncompress_registry(joinpath(dirname(path), compressed_file))
    else
        nothing
    end

    d = parsefile(in_memory_registry, path, "Registry.toml")
    pkgs = Dict{UUID,PkgEntry}()
    for (uuid, info) in d["packages"]::Dict{String,Any}
        uuid = UUID(uuid::String)
        info::Dict{String,Any}
        name = info["name"]::String
        pkgpath = info["path"]::String
        pkg = PkgEntry(pkgpath, path, name, uuid, in_memory_registry, uninit)
        pkgs[uuid] = pkg
    end
    reg = RegistryInstance(
        path,
        d["name"]::String,
        UUID(d["uuid"]::String),
get(d, "repo", nothing)::Union{String,Nothing},
        get(d, "description", nothing)::Union{String,Nothing},
        pkgs,
        tree_info,
        in_memory_registry,
        Dict{String,UUID}(),
    )
    if tree_info !== nothing
        REGISTRY_CACHE[path] = (tree_info, compressed_file !== nothing, reg)
    end
    return reg
end

function Base.show(io::IO, ::MIME"text/plain", r::RegistryInstance)
    println(io, "Registry: $(repr(r.name)) at $(repr(r.path)):")
    println(io, "  uuid: ", r.uuid)
    println(io, "  repo: ", r.repo)
    if r.tree_info !== nothing
        println(io, "  git-tree-sha1: ", r.tree_info)
    end
    println(io, "  packages: ", length(r.pkgs))
end

function uuids_from_name(r::RegistryInstance, name::String)
    create_name_uuid_mapping!(r)
    return get(Vector{UUID}, r.name_to_uuids, name)
end

function create_name_uuid_mapping!(r::RegistryInstance)
    isempty(r.name_to_uuids) || return
    for (uuid, pkg) in r.pkgs
        uuids = get!(Vector{UUID}, r.name_to_uuids, pkg.name)
        push!(uuids, pkg.uuid)
            end
    return
end

function registry_use_pkg_server()
    get(ENV, "JULIA_PKG_SERVER", nothing) !== ""
end

registry_read_from_tarball() =
    registry_use_pkg_server() && !(get(ENV, "JULIA_PKG_UNPACK_REGISTRY", "") == "true")

function verify_compressed_registry_toml(path::String)
    d = TOML.tryparsefile(path)
    if d isa TOML.ParserError
        @warn "Failed to parse registry TOML file at $(repr(path))" exception = d
        return false
    end
    for key in ("git-tree-sha1", "uuid", "path")
        if !haskey(d, key)
            @warn "Expected key $(repr(key)) to exist in registry TOML file at $(repr(path))"
            return false
        end
    end
    compressed_file = joinpath(dirname(path), d["path"])
    if !isfile(compressed_file)
        @warn "Expected the compressed registry for $(repr(path)) to exist at $(repr(compressed_file))"
        return false
    end
    return true
end

function reachable_registries(; depots::Union{String,Vector{String}}=Base.DEPOT_PATH)
    # collect registries
    if depots isa String
        depots = [depots]
    end
    registries = RegistryInstance[]
    for d in depots
        isdir(d) || continue
        reg_dir = joinpath(d, "registries")
        isdir(reg_dir) || continue
        reg_paths = readdir(reg_dir; join=true)
        candidate_registries = String[]
        # All folders could be registries
        append!(candidate_registries, filter(isdir, reg_paths))
        if registry_read_from_tarball()
            compressed_registries = filter(endswith(".toml"), reg_paths)
            # if we are reading compressed registries, ignore compressed registries
            # with the same name
            compressed_registry_names = Set([splitext(basename(file))[1] for file in compressed_registries])
            filter!(x -> !(basename(x) in compressed_registry_names), candidate_registries)
            for compressed_registry in compressed_registries
                if verify_compressed_registry_toml(compressed_registry)
                    push!(candidate_registries, compressed_registry)
                end
            end
        end

        for candidate in candidate_registries
            # candidate can be either a folder or a TOML file
            if isfile(joinpath(candidate, "Registry.toml")) || isfile(candidate)
                push!(registries, RegistryInstance(candidate))
            end
        end
    end
    return registries
end

# Dict interface

function Base.haskey(r::RegistryInstance, uuid::UUID)
    return haskey(r.pkgs, uuid)
end

function Base.keys(r::RegistryInstance)
    return keys(r.pkgs)
end

Base.getindex(r::RegistryInstance, uuid::UUID) = r.pkgs[uuid]

Base.get(r::RegistryInstance, uuid::UUID, default) = get(r.pkgs, uuid, default)

Base.iterate(r::RegistryInstance) = iterate(r.pkgs)
Base.iterate(r::RegistryInstance, state) = iterate(r.pkgs, state)

end


text:
@static if VERSION < v"0.7-"
    srand(12345)
else
    Random.seed!(12345)
end

@testset "gemv!" begin
    for elty in eltypes
        A = rand(elty, m_L2, n_L2)
        A_cl = cl.CLArray(queue, A)
        x = rand(elty, n_L2)
        x_cl = cl.CLArray(queue, x)
        y = rand(elty, m_L2)
        y_cl = cl.CLArray(queue, y)
        α = rand(elty)
        β = rand(elty)

        @test_throws DimensionMismatch CLBlast.gemv!('T', α, A_cl, x_cl, β, y_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.gemv!('C', α, A_cl, x_cl, β, y_cl, queue=queue)
        CLBlast.gemv!('N', α, A_cl, x_cl, β, y_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.gemv!('N', α, A, x, β, y)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y

        @test_throws DimensionMismatch CLBlast.gemv!('N', α, A_cl, y_cl, β, x_cl, queue=queue)
        CLBlast.gemv!('T', α, A_cl, y_cl, β, x_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.gemv!('T', α, A, y, β, x)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y

        @test_throws DimensionMismatch CLBlast.gemv!('N', α, A_cl, y_cl, β, x_cl, queue=queue)
        CLBlast.gemv!('C', α, A_cl, y_cl, β, x_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.gemv!('C', α, A, y, β, x)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y

        @test_throws ArgumentError CLBlast.gemv!('A', α, A_cl, y_cl, β, x_cl, queue=queue)
    end
end

@testset "gbmv!" begin
    for elty in eltypes
        A = rand(elty, kl+ku+1, n_L2)
        A_cl = cl.CLArray(queue, A)
        x = rand(elty, n_L2)
        x_cl = cl.CLArray(queue, x)
        y = rand(elty, m_L2)
        y_cl = cl.CLArray(queue, y)
        α = rand(elty)
        β = rand(elty)

        CLBlast.gbmv!('N', m_L2, kl, ku, α, A_cl, x_cl, β, y_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.gbmv!('N', m_L2, kl, ku, α, A, x, β, y)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y

        CLBlast.gbmv!('T', m_L2, kl, ku, α, A_cl, y_cl, β, x_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.gbmv!('T', m_L2, kl, ku, α, A, y, β, x)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y

        CLBlast.gbmv!('C', m_L2, kl, ku, α, A_cl, y_cl, β, x_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.gbmv!('C', m_L2, kl, ku, α, A, y, β, x)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y

        @test_throws ArgumentError CLBlast.gbmv!('A', m_L2, kl, ku, α, A_cl, y_cl, β, x_cl, queue=queue)
    end
end

@testset "hemv!" begin
    for elty in eltypes
        elty <: Complex || continue

        A = rand(elty, n_L2, n_L2)
        A_cl = cl.CLArray(queue, A)
        x = rand(elty, n_L2)
        x_cl = cl.CLArray(queue, x)
        y = rand(elty, n_L2)
        y_cl = cl.CLArray(queue, y)
        α = rand(elty)
        β = rand(elty)

        CLBlast.hemv!('U', α, A_cl, x_cl, β, y_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.hemv!('U', α, A, x, β, y)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y

        CLBlast.hemv!('L', α, A_cl, x_cl, β, y_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.hemv!('L', α, A, x, β, y)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y

        @test_throws ArgumentError CLBlast.hemv!('A', α, A_cl, y_cl, β, x_cl, queue=queue)

        y = rand(elty, m_L2)
        y_cl = cl.CLArray(queue, y)
        @test_throws DimensionMismatch CLBlast.hemv!('U', α, A_cl, x_cl, β, y_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.hemv!('U', α, A_cl, y_cl, β, x_cl, queue=queue)
    end
end

@testset "hbmv!" begin
    for elty in eltypes
        elty <: Complex || continue

        A = rand(elty, ku+1, n_L2)
        A_cl = cl.CLArray(queue, A)
        x = rand(elty, n_L2)
        x_cl = cl.CLArray(queue, x)
        y = rand(elty, n_L2)
        y_cl = cl.CLArray(queue, y)
        α = rand(elty)
        β = rand(elty)

        CLBlast.hbmv!('U', ku, α, A_cl, x_cl, β, y_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.hbmv!('U', ku, α, A, x, β, y)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y

        CLBlast.hbmv!('L', ku, α, A_cl, x_cl, β, y_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.hbmv!('L', ku, α, A, x, β, y)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y

        @test_throws ArgumentError CLBlast.hbmv!('A', ku, α, A_cl, x_cl, β, y_cl, queue=queue)

        y = rand(elty, m_L2)
        y_cl = cl.CLArray(queue, y)
        @test_throws DimensionMismatch CLBlast.hbmv!('U', ku, α, A_cl, x_cl, β, y_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.hbmv!('U', ku, α, A_cl, y_cl, β, x_cl, queue=queue)
    end
end

@testset "symv!" begin
    for elty in eltypes
        elty <: Real || continue

        A = rand(elty, n_L2, n_L2)
        A_cl = cl.CLArray(queue, A)
        x = rand(elty, n_L2)
        x_cl = cl.CLArray(queue, x)
        y = rand(elty, n_L2)
        y_cl = cl.CLArray(queue, y)
        α = rand(elty)
        β = rand(elty)

        CLBlast.symv!('U', α, A_cl, x_cl, β, y_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.symv!('U', α, A, x, β, y)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y

        CLBlast.symv!('L', α, A_cl, x_cl, β, y_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.symv!('L', α, A, x, β, y)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y

        @test_throws ArgumentError CLBlast.symv!('A', α, A_cl, y_cl, β, x_cl, queue=queue)

        y = rand(elty, m_L2)
        y_cl = cl.CLArray(queue, y)
        @test_throws DimensionMismatch CLBlast.symv!('U', α, A_cl, x_cl, β, y_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.symv!('U', α, A_cl, y_cl, β, x_cl, queue=queue)
    end
end

@testset "sbmv!" begin
    for elty in eltypes
        elty <: Real || continue

        A = rand(elty, ku+1, n_L2)
        A_cl = cl.CLArray(queue, A)
        x = rand(elty, n_L2)
        x_cl = cl.CLArray(queue, x)
        y = rand(elty, n_L2)
        y_cl = cl.CLArray(queue, y)
        α = rand(elty)
        β = rand(elty)

        CLBlast.sbmv!('U', ku, α, A_cl, x_cl, β, y_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.sbmv!('U', ku, α, A, x, β, y)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y

        CLBlast.sbmv!('L', ku, α, A_cl, x_cl, β, y_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.sbmv!('L', ku, α, A, x, β, y)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y

        @test_throws ArgumentError CLBlast.sbmv!('A', ku, α, A_cl, x_cl, β, y_cl, queue=queue)

        y = rand(elty, m_L2)
        y_cl = cl.CLArray(queue, y)
        @test_throws DimensionMismatch CLBlast.sbmv!('U', ku, α, A_cl, x_cl, β, y_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.sbmv!('U', ku, α, A_cl, y_cl, β, x_cl, queue=queue)
    end
end

@testset "trmv!" begin
    for elty in eltypes
        A = rand(elty, n_L2, n_L2)
        A_cl = cl.CLArray(queue, A)
        x = rand(elty, n_L2)
        x_cl = cl.CLArray(queue, x)

        for uplo in ['U','L'], trans in ['N','T','C'], diag in ['N','U']
            CLBlast.trmv!(uplo, trans, diag, A_cl, x_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.trmv!(uplo, trans, diag, A, x)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(x_cl, queue=queue) ≈ x
        end

        @test_throws ArgumentError CLBlast.trmv!('A', 'N', 'N', A_cl, x_cl, queue=queue)
        @test_throws ArgumentError CLBlast.trmv!('U', 'A', 'N', A_cl, x_cl, queue=queue)
        @test_throws ArgumentError CLBlast.trmv!('U', 'N', 'A', A_cl, x_cl, queue=queue)

        y = rand(elty, m_L2)
        y_cl = cl.CLArray(queue, y)
        @test_throws DimensionMismatch CLBlast.trmv!('U', 'N', 'N', A_cl, y_cl, queue=queue)

        A = rand(elty, m_L2, n_L2)
        A_cl = cl.CLArray(queue, A)
        @test_throws DimensionMismatch CLBlast.trmv!('U', 'N', 'N', A_cl, x_cl, queue=queue)
    end
end

@testset "trsv!" begin
    for elty in eltypes
        A = rand(elty, n_L2, n_L2)
        for i in 1:n_L2
            A[i,i] = i
        end
        A_cl = cl.CLArray(queue, A)
        x = rand(elty, n_L2)
        x_cl = cl.CLArray(queue, x)

        for uplo in ['U','L'], trans in ['N','T','C'], diag in ['N','U']
            CLBlast.trsv!(uplo, trans, diag, A_cl, x_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.trsv!(uplo, trans, diag, A, x)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(x_cl, queue=queue) ≈ x
        end

        @test_throws ArgumentError CLBlast.trsv!('A', 'N', 'N', A_cl, x_cl, queue=queue)
        @test_throws ArgumentError CLBlast.trsv!('U', 'A', 'N', A_cl, x_cl, queue=queue)
        @test_throws ArgumentError CLBlast.trsv!('U', 'N', 'A', A_cl, x_cl, queue=queue)

        y = rand(elty, m_L2)
        y_cl = cl.CLArray(queue, y)
        @test_throws DimensionMismatch CLBlast.trsv!('U', 'N', 'N', A_cl, y_cl, queue=queue)

        A = rand(elty, m_L2, n_L2)
        A_cl = cl.CLArray(queue, A)
        @test_throws DimensionMismatch CLBlast.trsv!('U', 'N', 'N', A_cl, x_cl, queue=queue)
    end
end

@testset "ger!" begin
    for elty in eltypes
        A = rand(elty, m_L2, n_L2)
        A_cl = cl.CLArray(queue, A)
        x = rand(elty, m_L2)
        x_cl = cl.CLArray(queue, x)
        y = rand(elty, n_L2)
        y_cl = cl.CLArray(queue, y)
        α = rand(elty)

        CLBlast.ger!(α, x_cl, y_cl, A_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.ger!(α, x, y, A)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(x_cl, queue=queue) ≈ x
        @test cl.to_host(y_cl, queue=queue) ≈ y

        @test_throws DimensionMismatch CLBlast.ger!(α, x_cl, x_cl, A_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.ger!(α, y_cl, y_cl, A_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.ger!(α, y_cl, x_cl, A_cl, queue=queue)
    end
end

@testset "her!" begin
    for elty in eltypes
        elty <: Complex || continue

        A = rand(elty, n_L2, n_L2)
        A_cl = cl.CLArray(queue, A)
        x = rand(elty, n_L2)
        x_cl = cl.CLArray(queue, x)
        α = real(rand(elty))

        for uplo in ['U', 'L']
            CLBlast.her!(uplo, α, x_cl, A_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.her!(uplo, α, x, A)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(x_cl, queue=queue) ≈ x
        end

        y = rand(elty, m_L2)
        y_cl = cl.CLArray(queue, y)
        @test_throws DimensionMismatch CLBlast.her!('U', α, y_cl, A_cl, queue=queue)
    end
end

@testset "syr!" begin
    for elty in eltypes
        elty <: Real || continue

        A = rand(elty, n_L2, n_L2)
        A_cl = cl.CLArray(queue, A)
        x = rand(elty, n_L2)
        x_cl = cl.CLArray(queue, x)
        α = real(rand(elty))

        for uplo in ['U', 'L']
            CLBlast.syr!(uplo, α, x_cl, A_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.syr!(uplo, α, x, A)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(x_cl, queue=queue) ≈ x
        end

        y = rand(elty, m_L2)
        y_cl = cl.CLArray(queue, y)
        @test_throws DimensionMismatch CLBlast.syr!('U', α, y_cl, A_cl, queue=queue)
    end
end


text:
module StringBenchmarks

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))

using .RandUtils
using BenchmarkTools

const SUITE = BenchmarkGroup()

####################
# replace (#12224) #
####################

str = join(samerand('a':'d', 10^4))

SUITE["replace"] = @benchmarkable replace($str, "a" => "b")
SUITE["join"] = @benchmarkable join($str, $str) time_tolerance=0.40

str = "Gf6FPPWevqer3di13haDSzrRrSiThqmV3k02dALLu7OHdYRR0dfrKf4iCMcDvgZBawx"

# The searchindex deprecation target makes the updated searchindex group
# redundant with the updated search group, so they've been combined here
g = addgroup!(SUITE, "findfirst")
g["Char"] = @benchmarkable findfirst(isequal($('x')), $str)
g["String"] = @benchmarkable findfirst($("x"), $str)

######################
# readuntil (#20621) #
######################

buffer = IOBuffer("A" ^ 50000)

g = addgroup!(SUITE, "readuntil")
for len in (1, 2, 1000, 50000)
    g["target length $len"] = @benchmarkable readuntil(seekstart($buffer), $("A" ^ len))
end

buffer = IOBuffer(("A" ^ 50000) * "B")
target = ("A" ^ 5000) * "Z"
g["backtracking"] = @benchmarkable readuntil(seekstart($buffer), $target)

buffer = IOBuffer(String(rand(RandUtils.SEED, 'A':'X', 40000)) * target)
target = "Y" * ("Z" ^ 999)
g["no backtracking"] = @benchmarkable readuntil(seekstart($buffer), $target)

buffer = IOBuffer(("bar" ^ 20000) * "ians")
target = ("bar" ^ 300) * "ian"
g["barbarian backtrack"] = @benchmarkable readuntil(seekstart($buffer), $target)

#################
# repeat #22462 #
#################

g = addgroup!(SUITE, "repeat")
g["repeat str len 1"] = @benchmarkable repeat(" ", 500)
g["repeat str len 16"] = @benchmarkable repeat("repeatmerepeatme", 500)
g["repeat char 1"] = @benchmarkable repeat(' ', 500)
g["repeat char 2"] = @benchmarkable repeat('α', 500)

end # module


text:
# this script generates the SnoopCompile precompile list

using SnoopCompile

### Log the compiles
# This only needs to be run once (to generate "/tmp/approxfun_compiles.csv")

SnoopCompile.@snoop "/tmp/approxfun_compiles.csv" begin
    using ApproxFun
    a=Fun()
    b=Fun(exp)
    a+b
    a*b
    a/b
    a-b

    # following causes SEGFAULT
    # S = Chebyshev()
    # D = Derivative(S)
    # B = dirichlet(S)
    # A = [B;D^2+I]
    # u = A\[1.0,0.0,b]
end

### Parse the compiles and generate precompilation scripts
# This can be run repeatedly to tweak the scripts

# IMPORTANT: we must have the module(s) defined for the parcelation
# step, otherwise we will get no precompiles for the ApproxFun module
using ApproxFun

data = SnoopCompile.read("/tmp/approxfun_compiles.csv")

# Blacklist helps fix problems:
# - MIME uses type-parameters with symbols like :image/png, which is
#   not parseable
blacklist = ["MIME"]

# Use these two lines if you want to create precompile functions for
# individual packages
pc, discards = SnoopCompile.parcel(data[end:-1:1,2], blacklist=blacklist)
SnoopCompile.write("/tmp/precompile", pc)


text:
"""
    Wing(left :: HalfWing, right :: HalfWing)

A composite type consisting of two `HalfWing`s with fields `left` and `right` for constructing a wing.
"""
struct Wing{T <: Real} <: Aircraft
    left :: HalfWing{T}
    right :: HalfWing{T}
end

# function Wing(left :: HalfWing{T}, right :: HalfWing{T}) where T <: Real
#     @assert first(left.chords) == first(right.chords) "Properties of first section must match."
#     # Add more conditions...
#     Wing{T}(left, right)
# end

left(wing :: Wing)      = wing.left
right(wing :: Wing)     = wing.right
foils(wing :: Wing)     = @views [ (foils ∘ left)(wing)[end:-1:2]	  ; (foils ∘ right)(wing)     ]
chords(wing :: Wing)    = @views [ (chords ∘ left)(wing)[end:-1:2]	  ; (chords ∘ right)(wing)    ]
twists(wing :: Wing)    = @views [ (twists ∘ left)(wing)[end:-1:2]	  ; (twists ∘ right)(wing)   ]
spans(wing :: Wing)     = @views [ (reverse ∘ spans ∘ left)(wing)	  ; (spans ∘ right)(wing)     ]
dihedrals(wing :: Wing) = @views [ (reverse ∘ dihedrals ∘ left)(wing) ; (dihedrals ∘ right)(wing) ]
sweeps(wing :: Wing)    = @views [ (reverse ∘ sweeps ∘ left)(wing)    ; (sweeps ∘ right)(wing)   ]

# Symmetric wing
Wing(; chords :: Vector{T}, twists :: Vector{T}, spans :: Vector{T}, dihedrals :: Vector{T}, sweep_LEs :: Vector{T}, foils :: Vector{Foil{T}} = fill(Foil(naca4((0,0,1,2))), length(chords))) where T <: Real = let w = HalfWing(foils = foils, chords = chords, twists = twists, spans = spans, dihedrals = dihedrals, sweep_LEs = sweep_LEs); Wing(w, w) end

# Single section for convenience
WingSection(; span = 1., dihedral = 0., sweep_LE = 0., taper = 1., root_chord = 1., root_twist = 0., tip_twist = 0., root_foil = naca4((0,0,1,2)), tip_foil = naca4((0,0,1,2))) = let w = HalfWingSection(span = span, dihedral = dihedral, sweep_LE = sweep_LE, taper = taper, root_chord = root_chord, root_twist = root_twist, tip_twist = tip_twist, root_foil = root_foil, tip_foil = tip_foil); Wing(w, w) end

"""
    span(wing :: Wing)

Compute the span of a `Wing`.
"""
span(wing :: Wing) = (span ∘ left)(wing) + (span ∘ right)(wing)

taper_ratio(wing :: Wing) = (taper_ratio ∘ left)(wing), (taper_ratio ∘ right)(wing)

"""
    projected_area(wing :: Wing)

Compute the projected_area of a `Wing`.
"""
projected_area(wing :: Wing) = (projected_area ∘ left)(wing) + (projected_area ∘ right)(wing)

"""
    mean_aerodynamic_chord(wing :: Wing)

Compute the mean aerodynamic chord of a `Wing`.
"""
mean_aerodynamic_chord(wing :: Wing) = ((mean_aerodynamic_chord ∘ left)(wing) + (mean_aerodynamic_chord ∘ right)(wing)) / 2

"""
    wing_bounds(wing :: Wing)

Return the leading and trailing edge coordinates of a `Wing`.
"""
function wing_bounds(wing :: Wing)
    left_lead, left_trail 	= wing_bounds(left(wing), true)
    right_lead, right_trail = wing_bounds(right(wing))

    leading  = @views [ left_lead[1:end-1] ; right_lead  ]
    trailing = @views [ left_trail[1:end-1]; right_trail ]

    leading, trailing
end

leading_edge(wing :: Wing)  = @views [ leading_edge(left(wing), true)[1:end-1] ; leading_edge(right(wing))  ]
trailing_edge(wing :: Wing) = @views [ trailing_edge(left(wing), true)[1:end-1]; trailing_edge(right(wing)) ]

mean_aerodynamic_center(wing :: Wing, factor = 0.25) = (mean_aerodynamic_center(right(wing), factor) .+ mean_aerodynamic_center(left(wing), factor)) ./ 2

max_forks_repo_name: GodotMisogi/AeroMDAO