text:
using Distributions, StatsBase, ProgressMeter, InvertedIndices, IterTools

export draw_sample_mode!, draw_sample_mode
export draw_sample_gamma!, draw_sample_gamma
export rand_multivariate_bernoulli, rand_multinomial_dict
export flip!, multinomial_flip!, flip_informed!
export rand_restr_bins, flip_informed_excl!

# NOTE μ_cusum must have 0.0 first entry, and for n bins μ_cusum must be of size n+1
function rand_multivariate_bernoulli(μ_cusum::Vector{Float64})
    @assert μ_cusum[1] ≈ 0 "First entry must be 0.0 (for differencing to find probabilities)."
    β = rand()
    for i in 1:length(μ_cusum)
        if β < μ_cusum[i]
            return i-1, μ_cusum[i]-μ_cusum[i-1]
        else
            continue
        end
    end
end

function rand_multinomial_dict(μ_cusum::Vector{Float64}, ntrials::Int)
    @assert μ_cusum[1] ≈ 0 "First entry must be 0.0 (for differencing to find probabilities)."
    out = Dict{Int,Int}()
    for i in 1:ntrials
        β = rand()
        j = findfirst(x->x>β, μ_cusum)
        out[j-1] = get(out, j-1, 0) + 1
    end
    return out
end

function rand_restr_bins(
    bins::Vector{Int},
    n::Int
    )
    cum_bins = cumsum(bins)
    @assert n <= cum_bins[end] "Cannot sample more than capacity of bins."
    ind_flat = zeros(Int, n)
    StatsBase.seqsample_a!(1:cum_bins[end], ind_flat)
    out = Dict{Int,Int}()
    for i in ind_flat
        j = findfirst(x->x>=i, cum_bins)
        out[j] = get(out,j,0) + 1
    end
    return out
end

function rand_restr_cum_bins(
    cum_bins::Vector{Int},
    n::Int
    )
    @assert n <= cum_bins[end] "Cannot sample more than capacity of bins."
    ind_flat = zeros(Int, n)
    StatsBase.seqsample_a!(1:cum_bins[end], ind_flat)
    out = Dict{Int,Int}()
    for i in ind_flat
        j = findfirst(x->x>=i, cum_bins)
        out[j] = get(out,j,0) + 1
    end
    return out
end

function delete_insert_informed!(
    x::Path,
    ind_del::AbstractArray{Int},
    ind_add::AbstractArray{Int},
    vals_del::AbstractArray{Int},
    P::CumCondProbMatrix
    )

    @views for (i, index) in enumerate(ind_del)
        tmp_ind = index - i + 1  # Because size is changing must adapt index
        vals_del[i] = x[tmp_ind]
        deleteat!(x, tmp_ind)
    end

    # Now find distrbution for insertions via co-occurence
    curr_vertices = unique(x)
    if length(curr_vertices) == 0
        V = size(P)[2]
        tmp = fill(1/V, V)
        pushfirst!(tmp, 0.0)
        μ_cusum = cumsum(tmp)
    else
        μ_cusum = sum(P[:,curr_vertices], dims=2)[:] ./ length(curr_vertices)
    end

    # @show μ_cusum

    log_ratio_entries = 0.0
    # Add probability of values deleted
    for v in vals_del
        log_ratio_entries += log(μ_cusum[v+1]-μ_cusum[v])
    end

    # Sample new entries and add probabilility
    @views for index in ind_add
        # @show i, index, val
        val, prob = rand_multivariate_bernoulli(μ_cusum)
        insert!(x, index, val)
        log_ratio_entries += -log(prob)
    end

    return log_ratio_entries
end

function flip!(
    x::Path,
    ind_flip::AbstractArray{Int},
    V::Vector{Int}
    )

    # Get unique vertices not being flipped
    for i in ind_flip
        x[i] = curr_val
        curr_vert_ind = findfirst(x->x==curr_val, V)
        tmp = rand(1:(length(V)-1))

        if tmp >= curr_vert_ind
            x[i] = V[tmp+1]
        else
            x[i] = V[tmp]
        end

    end

end

function flip!(
    x::Path,
    ind_flip::AbstractArray{Int},
    V::UnitRange{Int}
    )

    # Get unique vertices not being flipped
    for i in ind_flip
        x[i] = curr_val
        tmp = rand(1:(length(V)-1))
        if tmp >= curr_val
            x[i] = tmp+1
        else
            x[i] = tmp
        end
    end

end

function flip_informed!(
    x::Path,
    ind_flip::AbstractArray{Int},
    P::CumCondProbMatrix
    )

    curr_vertices = unique(view(x, Not(ind_flip)))
    if length(curr_vertices) == 0
        V = size(P)[2]
        tmp = fill(1/V, V)
        pushfirst!(tmp, 0.0)
        μ_cusum = cumsum(tmp)
    else
        μ_cusum = sum(P[:,curr_vertices], dims=2)[:] ./ length(curr_vertices)
    end
    log_ratio = 0.0
    for index in ind_flip
        curr_val = x[index]
        x[index], prob = rand_multivariate_bernoulli(μ_cusum)
        log_ratio += (
            log(μ_cusum[curr_val+1]-μ_cusum[curr_val]) - log(prob)
        )
    end
    return log_ratio
end

function flip_informed!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    mcmc::Union{SisIexInsertDelete,SimIexInsertDelete},
    P::CumCondProbMatrix
    )

    δ = rand(1:mcmc.ν_ed)
    alloc = rand_restr_bins(length.(S_curr), δ)
    ind = mcmc.ind_add
    log_ratio = 0.0

    for (key,val) in pairs(alloc)
        ind_flip = view(ind, 1:val)
        log_ratio += flip_informed!(
            S_prop[key],
            ind_flip,
            P
        )
    end
    return log_ratio

end

function flip_informed_excl!(
    x::Path,
    ind_flip::AbstractArray{Int},
    P::CumCondProbMatrix
    )

    curr_vertices = unique(view(x, Not(ind_flip)))
    if length(curr_vertices) == 0
        V = size(P)[2]
        tmp = fill(1/V, V)
        pushfirst!(tmp, 0.0)
        μ_cusum = cumsum(tmp)
    else
        μ_cusum = sum(P[:,curr_vertices], dims=2)[:] ./ length(curr_vertices)
    end
    log_ratio = 0.0
    for index in ind_flip
        curr_val = x[index]
        prop_val = curr_val
        prop_prob = 0.0
        while prop_val == curr_val
            prop_val, prop_prob = rand_multivariate_bernoulli(μ_cusum)
        end
        x[index] = prop_val
        curr_prob = μ_cusum[curr_val+1]-μ_cusum[curr_val]
        log_ratio += (
            log(curr_prob) + log(1-curr_prob)
            - log(prop_prob) - log(1-prop_prob)
        )
    end
    return log_ratio
end

function flip_informed_excl!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    mcmc::Union{SisIexInsertDelete,SimIexInsertDelete},
    P::CumCondProbMatrix
    )

    δ = rand(1:mcmc.ν_ed)
    alloc = rand_restr_bins(length.(S_curr), δ)
    ind = mcmc.ind_add
    log_ratio = 0.0

    for (key,val) in pairs(alloc)
        ind_flip = view(ind, 1:val)
        log_ratio += flip_informed_excl!(
            S_prop[key],
            ind_flip,
            P
        )
    end
    return log_ratio

end

function double_iex_multinomial_edit_accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SisPosterior,
    γ_curr::Float64,
    mcmc::SisIexInsertDelete,
    P::CumCondProbMatrix,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool
    )

    N = length(S_curr)
    dist = posterior.dist
    V = posterior.V
    K_inner, K_outer = (posterior.K_inner, posterior.K_outer)
    K_in_lb, K_in_ub = (K_inner.l, K_inner.u)
    data = posterior.data
    γ_prior = posterior.S_prior.γ
    mode_prior = posterior.S_prior.mode

    aux_mcmc = mcmc.aux_mcmc

    δ = rand(1:mcmc.ν_ed)  # Number of edits to enact
    rem_edits = δ # Remaining edits to allocate
    j = 0 # Keeps track how many interaction have been edited
    log_ratio = 0.0

    for i = 1:N

        # If at end we just assign all remaining edits to final interaction
        if i == N
            δ_tmp = rem_edits
        # Otherwise we sample the number of edits via rescaled Binomial
        else
            p = 1/(N-i+1)
            δ_tmp = rand(Binomial(rem_edits, p)) # Number of edits to perform on ith interaction
        end

        # println("   Index $i getting $δ_tmp edits...")
        # If we sampled zero edits we skip to next iteration
        if δ_tmp == 0
            continue
        else
            j += 1 # Increment j
            # Make edits ....
            @inbounds n = length(S_curr[i])
            # a,b = (lb(n, δ_tmp, model), ub(n, δ_tmp))
            d = rand(0:min(n-K_in_lb, δ_tmp))
            m = n + δ_tmp - 2*d

            # Catch invalid proposals
            if  (m > K_in_ub)
                # Here we just reject the proposal
                for i in 1:N
                    copy!(S_prop[i], S_curr[i])
                end
                return 0, suff_stat_curr
            end

            # tot_dels += d
            # println("       Deleting $d and adding $(δ_tmp-d)")
            ind_del = view(mcmc.ind_del, 1:d)
            ind_add = view(mcmc.ind_add, 1:(δ_tmp-d))
            vals_del = view(mcmc.vals, 1:d)

            # println("           ind_del: $ind_del ; ind_add: $ind_add")

            # Sample indexing info and new entries (all in-place)
            StatsBase.seqsample_a!(1:n, ind_del)
            StatsBase.seqsample_a!(1:m, ind_add)

            # *** HERE IS DIFFERENT FROM MODEL SAMPLER ***
            # The delete_insert_informed() function does the sampling + editing
            log_ratio += delete_insert_informed!(
                S_prop[i],
                ind_del, ind_add, vals_del,
                P)

            mcmc.ind_update[j] = i # Store which interaction was updated

            # Add to log_ratio
            # log_prod_term += log(b - a + 1) - log(ub(m, δ_tmp) - lb(m, δ_tmp, model) +1)
            log_ratio += log(min(n-K_in_lb, δ_tmp)+1) - log(min(m-K_in_lb, δ_tmp)+1)

        end

        # Update rem_edits
        rem_edits -= δ_tmp

        # If no more left terminate
        if rem_edits == 0
            break
        end

    end

    aux_model = SIS(
        S_prop, γ_curr,
        dist,
        V,
        K_inner,
        K_outer
        )

    if aux_init_at_prev
        @inbounds tmp = deepcopy(aux_data[end])
        draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux_data, aux_mcmc, aux_model)
    end

    aux_log_lik_ratio = -γ_curr * (
        mapreduce(x -> dist(x, S_curr), + , aux_data)
        - mapreduce(x -> dist(x, S_prop), +, aux_data)
    )

    suff_stat_prop = mapreduce(x -> dist(x, S_prop), + , data)
    log_lik_ratio = -γ_curr * (
        suff_stat_prop - suff_stat_curr
    )

    log_prior_ratio = -γ_prior * (
        dist(S_prop, mode_prior) - dist(S_curr, mode_prior)
    )

    # Log acceptance probability
    log_α = log_lik_ratio + log_prior_ratio + aux_log_lik_ratio + log_ratio

    # Accept-reject step. Use info in mcmc.ind_update to know which interaction are to be copied over
    if log(rand()) < log_α
        for i in view(mcmc.ind_update, 1:j)
            copy!(S_curr[i], S_prop[i])
        end
        return 1, suff_stat_prop
    else
        for i in view(mcmc.ind_update, 1:j)
            copy!(S_prop[i], S_curr[i])
        end
        return 0, suff_stat_curr
    end
end

function double_iex_flip_accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SisPosterior,
    γ_curr::Float64,
    mcmc::SisIexInsertDelete,
    P::CumCondProbMatrix,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool
    )

    dist = posterior.dist
    V = posterior.V
    data = posterior.data
    γ_prior = posterior.S_prior.γ
    mode_prior = posterior.S_prior.mode
    K_inner, K_outer = (posterior.K_inner, posterior.K_outer)

    aux_mcmc = mcmc.aux_mcmc
    lengths = length.(S_curr)

    δ = rand(1:min(mcmc.ν_ed,sum(lengths)))
    # @show δ, lengths
    alloc = rand_restr_bins(lengths, δ)
    ind = mcmc.ind_add
    log_ratio = 0.0

    for (key,val) in pairs(alloc)
        ind_flip = view(ind, 1:val)
        StatsBase.seqsample_a!(1:lengths[key], ind_flip)
        log_ratio += flip_informed_excl!(
            S_prop[key],
            ind_flip,
            P
        )
    end

    aux_model = SIS(
        S_prop, γ_curr,
        dist,
        V,
        K_inner,
        K_outer
        )

    if aux_init_at_prev
        tmp = deepcopy(aux_data[end])
        draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux_data, aux_mcmc, aux_model)
    end

    aux_log_lik_ratio = -γ_curr * (
        mapreduce(x -> dist(x, S_curr), + , aux_data)
        - mapreduce(x -> dist(x, S_prop), +, aux_data)
    )

    suff_stat_prop = mapreduce(x -> dist(x, S_prop), + , data)
    log_lik_ratio = -γ_curr * (
        suff_stat_prop - suff_stat_curr
    )

    log_prior_ratio = -γ_prior * (
        dist(S_prop, mode_prior) - dist(S_curr, mode_prior)
    )

    # Log acceptance probability
    log_α = log_lik_ratio + log_prior_ratio + aux_log_lik_ratio + log_ratio

    # Accept-reject step. Use info in mcmc.ind_update to know which interaction are to be copied over
    if log(rand()) < log_α
        for i in keys(alloc)
            copy!(S_curr[i], S_prop[i])
        end
        return 1, suff_stat_prop
    else
        for i in keys(alloc)
            copy!(S_prop[i], S_curr[i])
        end
        return 0, suff_stat_curr
    end

end

# Insertion/deletion proposal generation functions
# ================================================

# With vector of indices (uninformed)
# -----------------------------------
function imcmc_multi_insert_prop_sample!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    mcmc::T,
    ind::AbstractVector{Int},
    V::UnitRange,
    K_in_ub::Int
    ) where {T<:Union{SisPosteriorSampler,SimPosteriorSampler}}

    prop_pointers = mcmc.prop_pointers
    ν_td = mcmc.ν_td
    N = length(S_curr)
    len_dist = mcmc.len_dist

    log_ratio = 0.0
    for i in ind
        tmp = popfirst!(prop_pointers)
        m = rand(len_dist)
        resize!(tmp, m)
        sample!(V, tmp)
        insert!(S_prop, i, tmp)
        log_ratio += - logpdf(len_dist, m) + m*log(length(V)) - Inf * (m > K_in_ub)
    end
    return log_ratio

end

function imcmc_multi_delete_prop_sample!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    mcmc::T,
    ind::AbstractVector{Int},
    V::UnitRange
    ) where {T<:Union{SisPosteriorSampler,SimPosteriorSampler}}

    prop_pointers = mcmc.prop_pointers
    ν_td = mcmc.ν_td
    len_dist = mcmc.len_dist
    N = length(S_curr)

    log_ratio = 0.0

    for i in Iterators.reverse(ind)
        tmp = popat!(S_prop, i)
        pushfirst!(prop_pointers, tmp)
        m = length(tmp)
        log_ratio += logpdf(len_dist, m) - m * log(length(V))
    end
    return log_ratio

end

# With number of insertions/deletions
# -----------------------------------
function imcmc_multi_insert_prop_sample!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    mcmc::T,
    ε::Int,
    V::UnitRange,
    K_in_ub::Int
    ) where {T<:Union{SisPosteriorSampler,SimPosteriorSampler}}

    prop_pointers = mcmc.prop_pointers
    len_dist = mcmc.len_dist
    log_ratio = 0.0
    ind = mcmc.ind_td_add

    n = length(S_prop)+ε
    k = ε
    i = 0
    j = 0
    while k > 0
        u = rand()
        q = (n - k) / n
        while q > u  # skip
            i += 1
            n -= 1
            q *= (n - k) / n
        end
        i += 1
        j += 1
        # Insert path at index i
        # i is now index to insert
        ind[j] = i
        tmp = popfirst!(prop_pointers) # Get storage
        m = rand(len_dist)  # Sample length
        resize!(tmp, m) # Resize
        sample!(V, tmp) # Sample new entries uniformly
        @inbounds insert!(S_prop, i, tmp) # Insert path into S_prop
        log_ratio += - logpdf(len_dist, m) + m*log(length(V)) - Inf * (m > K_in_ub)  # Add to log_ratio term
        n -= 1
        k -= 1
    end
    return log_ratio

end

function imcmc_multi_delete_prop_sample!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    mcmc::T,
    ε::Int,
    V::UnitRange
    ) where {T<:Union{SisPosteriorSampler,SimPosteriorSampler}}

    prop_pointers = mcmc.prop_pointers
    len_dist = mcmc.len_dist
    log_ratio = 0.0
    ind = mcmc.ind_td_del # Store which entries were deleted

    n = length(S_curr)
    k = ε
    i = 0
    j = 0
    live_index = 0
    while k > 0
        u = rand()
        q = (n - k) / n
        while q > u  # skip
            i += 1
            n -= 1
            q *= (n - k) / n
        end
        i+=1
        j+=1
        # Delete path
        # i is now index to delete
        @inbounds ind[j] = i
        @inbounds tmp = popat!(S_prop, i - live_index)
        pushfirst!(prop_pointers, tmp)
        m = length(tmp)
        log_ratio += logpdf(len_dist, m) - m * log(length(V))
        live_index += 1
        n -= 1
        k -= 1
    end
    return log_ratio
end

# Informed insertions/deletions
# -----------------------------

function imcmc_multi_insert_prop_sample_informed!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    mcmc::T,
    ε::Int,
    p_ins::Categorical,
    K_in_ub::Int
    ) where {T<:Union{SisPosteriorSampler,SimPosteriorSampler}}

    prop_pointers = mcmc.prop_pointers
    len_dist = mcmc.len_dist
    log_ratio = 0.0
    ind = mcmc.ind_td_add

    n = length(S_curr)+ε
    k = ε
    i = 0
    j = 0
    while k > 0
        u = rand()
        q = (n - k) / n
        while q > u  # skip
            i += 1
            n -= 1
            q *= (n - k) / n
        end
        i += 1
        j += 1
        # Insert path at index i
        # i is now index to insert
        ind[j] = i
        tmp = popfirst!(prop_pointers) # Get storage
        m = rand(len_dist)  # Sample length
        resize!(tmp, m) # Resize
        for i in eachindex(tmp) # Sample new entries
            v = rand(p_ins)
            tmp[i] = v
            log_ratio += - logpdf(p_ins, v)
        end
        @inbounds insert!(S_prop, i, tmp) # Insert path into S_prop
        log_ratio += - logpdf(len_dist, m) - Inf * (m > K_in_ub)  # Add to log_ratio term
        n -= 1
        k -= 1
    end
    return log_ratio

end

function imcmc_multi_delete_prop_sample_informed!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    mcmc::T,
    ε::Int,
    p_ins::Categorical
    ) where {T<:Union{SisPosteriorSampler,SimPosteriorSampler}}

    prop_pointers = mcmc.prop_pointers
    len_dist = mcmc.len_dist
    log_ratio = 0.0
    ind = mcmc.ind_td_del # Store which entries were deleted

    n = length(S_curr)
    k = ε
    i = 0
    j = 0
    live_index = 0
    while k > 0
        u = rand()
        q = (n - k) / n
        while q > u  # skip
            i += 1
            n -= 1
            q *= (n - k) / n
        end
        i+=1
        j+=1
        # Delete path
        # i is now index to delete
        @inbounds ind[j] = i
        @inbounds tmp = popat!(S_prop, i - live_index)
        pushfirst!(prop_pointers, tmp)
        m = length(tmp)
        for v in tmp
            log_ratio += logpdf(p_ins,v)
        end
        log_ratio += logpdf(len_dist, m)
        # Update counters
        live_index += 1
        n -= 1
        k -= 1
    end
    return log_ratio
end

function double_iex_trans_dim_accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SisPosterior,
    γ_curr::Float64,
    mcmc::SisIexInsertDelete,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool
    )
    K_inner, K_outer = (posterior.K_inner, posterior.K_outer)
    K_out_lb, K_out_ub = (K_outer.l, K_outer.u)
    K_in_ub = K_inner.u
    data = posterior.data
    dist = posterior.dist
    V = posterior.V
    γ_prior = posterior.S_prior.γ
    mode_prior = posterior.S_prior.mode

    ν_td = mcmc.ν_td
    curr_pointers = mcmc.curr_pointers
    prop_pointers = mcmc.prop_pointers
    aux_mcmc = mcmc.aux_mcmc

    log_ratio = 0.0

    # Enact insertion / deletion
    N = length(S_curr)
    ε = rand(1:ν_td)
    d = rand(0:min(ν_td, N))
    a = ε - d
    # Catch invalid proposal (outside dimension bounds)
    M = N - d + a  # Number of paths in proposal
    if (M < K_out_lb) | (M > K_out_ub)
        return 0, suff_stat_curr
    end
    if d > 0
        log_ratio += imcmc_multi_delete_prop_sample!(
            S_curr, S_prop,
            mcmc,
            d,
            V
        )
    end
    if a > 0
        log_ratio += imcmc_multi_insert_prop_sample!(
            S_curr, S_prop,
            mcmc,
            a,
            V, K_in_ub
        )
    end

    log_ratio += log(min(ν_td, N) + 1) - log(min(ν_td, M) + 1)

    # Now do accept-reject step (**THIS IS WHERE WE DIFFER FROM MODEL SAMPLER***)
    aux_model = SIS(
        S_prop, γ_curr,
        dist,
        V,
        K_inner,
        K_outer
        )

    if aux_init_at_prev
        tmp = deepcopy(aux_data[end])
        draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux_data, aux_mcmc, aux_model)
    end

    aux_log_lik_ratio = -γ_curr * (
        mapreduce(x -> dist(x, S_curr), + , aux_data)
        - mapreduce(x -> dist(x, S_prop), +, aux_data)
    )
    suff_stat_prop = mapreduce(x -> dist(x, S_prop), + , data)
    log_lik_ratio = -γ_curr * (
        suff_stat_prop - suff_stat_curr
    )

    log_prior_ratio = -γ_prior * (
        dist(S_prop, mode_prior) - dist(S_curr, mode_prior)
    )

    # Log acceptance probability
    log_α = log_lik_ratio + log_prior_ratio + aux_log_lik_ratio + log_ratio

    # Note that we copy interactions between S_prop (resp. S_curr) and prop_pointers (resp .curr_pointers) by hand.
    if log(rand()) < log_α
        # Do deletions to S_curr
        ind = view(mcmc.ind_td_del, 1:d)
        for i in Iterators.reverse(ind) # If we didnt do reverse would have to update indices
            @inbounds tmp = popat!(S_curr, i)
            pushfirst!(curr_pointers, tmp)
        end
        # Do insertions
        ind = view(mcmc.ind_td_add, 1:a)
        for i in ind
            tmp = popfirst!(curr_pointers)
            copy!(tmp, S_prop[i])
            insert!(S_curr, i, tmp)
        end
        return 1, suff_stat_prop
    else
        # Remove insertions and return to prop_pointers
        ind = view(mcmc.ind_td_add, 1:a)
        for i in Iterators.reverse(ind)
            @inbounds tmp = popat!(S_prop, i)
            pushfirst!(prop_pointers, tmp)
        end
        # Re-insert deletions
        ind = view(mcmc.ind_td_del, 1:d)
        for i in ind
            @inbounds tmp = popfirst!(prop_pointers)
            copy!(tmp, S_curr[i])
            insert!(S_prop, i, tmp)
        end
        return 0, suff_stat_curr
    end
end

function double_iex_trans_dim_informed_accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SisPosterior,
    γ_curr::Float64,
    mcmc::SisIexInsertDelete,
    p_ins::Categorical,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool
    )
    K_inner, K_outer = (posterior.K_inner, posterior.K_outer)
    K_out_lb, K_out_ub = (K_outer.l, K_outer.u)
    K_in_ub = K_inner.u
    data = posterior.data
    dist = posterior.dist
    V = posterior.V
    γ_prior = posterior.S_prior.γ
    mode_prior = posterior.S_prior.mode

    ν_td = mcmc.ν_td
    curr_pointers = mcmc.curr_pointers
    prop_pointers = mcmc.prop_pointers
    aux_mcmc = mcmc.aux_mcmc

    log_ratio = 0.0

    # Enact insertion / deletion
    N = length(S_curr)
    ε = rand(1:ν_td)
    d = rand(0:min(ν_td, N))
    a = ε - d
    # Catch invalid proposal (outside dimension bounds)
    M = N - d + a  # Number of paths in proposal
    if (M < K_out_lb) | (M > K_out_ub)
        return 0, suff_stat_curr
    end
    if d > 0
        log_ratio += imcmc_multi_delete_prop_sample_informed!(
            S_curr, S_prop,
            mcmc,
            d,
            p_ins
        )
    end
    if a > 0
        log_ratio += imcmc_multi_insert_prop_sample_informed!(
            S_curr, S_prop,
            mcmc,
            a,
            p_ins,
            K_in_ub
        )
    end

    log_ratio += log(min(ν_td, N) + 1) - log(min(ν_td, M) + 1)

    # Now do accept-reject step (**THIS IS WHERE WE DIFFER FROM MODEL SAMPLER***)
    aux_model = SIS(
        S_prop, γ_curr,
        dist,
        V,
        K_inner,
        K_outer
        )

    if aux_init_at_prev
        tmp = deepcopy(aux_data[end])
        draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux_data, aux_mcmc, aux_model)
    end

    aux_log_lik_ratio = -γ_curr * (
        mapreduce(x -> dist(x, S_curr), + , aux_data)
        - mapreduce(x -> dist(x, S_prop), +, aux_data)
    )
    suff_stat_prop = mapreduce(x -> dist(x, S_prop), + , data)
    log_lik_ratio = -γ_curr * (
        suff_stat_prop - suff_stat_curr
    )

    log_prior_ratio = -γ_prior * (
        dist(S_prop, mode_prior) - dist(S_curr, mode_prior)
    )

    # Log acceptance probability
    log_α = log_lik_ratio + log_prior_ratio + aux_log_lik_ratio + log_ratio

    # @show is_insert, suff_stat_curr, suff_stat_prop, exp(log_α)
    # Note that we copy interactions between S_prop (resp. S_curr) and prop_pointers (resp .curr_pointers) by hand.
    # Note that we copy interactions between S_prop (resp. S_curr) and prop_pointers (resp .curr_pointers) by hand.
    if log(rand()) < log_α
        # Do deletions to S_curr
        ind = view(mcmc.ind_td_del, 1:d)
        for i in Iterators.reverse(ind) # If we didnt do reverse would have to update indices
            @inbounds tmp = popat!(S_curr, i)
            pushfirst!(curr_pointers, tmp)
        end
        # Do insertions
        ind = view(mcmc.ind_td_add, 1:a)
        for i in ind
            tmp = popfirst!(curr_pointers)
            copy!(tmp, S_prop[i])
            insert!(S_curr, i, tmp)
        end
        return 1, suff_stat_prop
    else
        # Remove insertions and return to prop_pointers
        ind = view(mcmc.ind_td_add, 1:a)
        for i in Iterators.reverse(ind)
            @inbounds tmp = popat!(S_prop, i)
            pushfirst!(prop_pointers, tmp)
        end
        # Re-insert deletions
        ind = view(mcmc.ind_td_del, 1:d)
        for i in ind
            @inbounds tmp = popfirst!(prop_pointers)
            copy!(tmp, S_curr[i])
            insert!(S_prop, i, tmp)
        end
        return 0, suff_stat_curr
    end
end

# =========================
# Samplers
# =========================

# Mode conditional
# ----------------

function draw_sample_mode!(
    sample_out::Union{InteractionSequenceSample{Int}, SubArray},
    mcmc::SisIexInsertDelete,
    posterior::SisPosterior,
    γ_fixed::Float64;
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::Vector{Path{Int}}=sample_frechet_mean(posterior.data, posterior.dist),
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    )

    if loading_bar
        iter = Progress(
            length(sample_out), # How many iters
            1,  # At which granularity to update loading bar
            "Chain for γ = $(γ_fixed) and n = $(posterior.sample_size) (mode conditional)....")  # Loading bar. Minimum update interval: 1 second
    end

    # Define aliases for pointers to the storage of current vals and proposals
    curr_pointers = mcmc.curr_pointers
    prop_pointers = mcmc.prop_pointers
    β = mcmc.β
    aux_mcmc = mcmc.aux_mcmc

    S_curr = InteractionSequence{Int}()
    S_prop = InteractionSequence{Int}()
    for i in 1:length(S_init)
        migrate!(S_curr, curr_pointers, i, 1)
        migrate!(S_prop, prop_pointers, i, 1)
        copy!(S_curr[i], S_init[i])
        copy!(S_prop[i], S_init[i])
    end

    γ_curr = γ_fixed

    sample_count = 1 # Keeps which sample to be stored we are working to get
    i = 0 # Keeps track all samples (included lags and burn_ins)

    tr_dim_count = 0
    tr_dim_acc_count = 0
    ed_count = 0
    ed_acc_count = 0
    flp_count = 0
    flp_acc_count = 0

    aux_data = [[Int[]] for i in 1:posterior.sample_size]
    # Initialise the aux_data
    aux_model = SIS(
        S_curr, γ_curr,
        posterior.dist,
        posterior.V,
        posterior.K_inner,
        posterior.K_outer)
    draw_sample!(aux_data, aux_mcmc, aux_model)

    # Evaluate sufficient statistic
    suff_stat_curr = mapreduce(
        x -> posterior.dist(S_curr, x),
        +,
        posterior.data
        )
    suff_stats = Float64[suff_stat_curr] # Storage for all sufficient stats (for diagnostics)

    P, vmap, vmap_inv = get_informed_proposal_matrix(posterior, mcmc.α)
    p_ins = get_informed_insertion_dist(posterior, mcmc.α)
    while sample_count ≤ length(sample_out)
        i += 1
        # Store value
        if (i > burn_in) & (((i-1) % lag)==0)
            sample_out[sample_count] = deepcopy(S_curr)
            push!(suff_stats, suff_stat_curr)
            sample_count += 1
        end
        # W.P. do update move (accept-reject done internally by function call)
        if rand() < β
            if rand() < 0.5
                was_acc, suff_stat_curr = double_iex_multinomial_edit_accept_reject!(
                    S_curr, S_prop,
                    posterior, γ_curr,
                    mcmc, P,
                    aux_data,
                    suff_stat_curr,
                    aux_init_at_prev
                )
                ed_acc_count += was_acc
                ed_count += 1
            else
                was_acc, suff_stat_curr = double_iex_flip_accept_reject!(
                    S_curr, S_prop,
                    posterior, γ_curr,
                    mcmc, P,
                    aux_data,
                    suff_stat_curr,
                    aux_init_at_prev
                )
                flp_acc_count += was_acc
                flp_count += 1
            end
        # Else do trans-dim move. We will do accept-reject move here
        else
            was_acc, suff_stat_curr = double_iex_trans_dim_informed_accept_reject!(
                S_curr, S_prop,
                posterior, γ_curr,
                mcmc,
                p_ins,
                aux_data,
                suff_stat_curr,
                aux_init_at_prev
            )
            tr_dim_acc_count += was_acc
            tr_dim_count += 1
        end
        if loading_bar
            next!(iter)
        end
    end
    for i in 1:length(S_curr)
        migrate!(curr_pointers, S_curr, 1, 1)
        migrate!(prop_pointers, S_prop, 1, 1)
    end
    return (
                ed_count, ed_acc_count,
                flp_count, flp_acc_count,
                tr_dim_count, tr_dim_acc_count,
                suff_stats
            )
end

function draw_sample_mode(
    mcmc::SisIexInsertDelete,
    posterior::SisPosterior,
    γ_fixed::Float64;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::Vector{Path{Int}}=sample_frechet_mean(posterior.data, posterior.dist),
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    )

    sample_out = Vector{Vector{Path{Int}}}(undef, desired_samples)
    draw_sample_mode!(
        sample_out,
        mcmc, posterior,
        γ_fixed,
        burn_in=burn_in, lag=lag, S_init=S_init,
        loading_bar=loading_bar,
        aux_init_at_prev=aux_init_at_prev
        )
    return sample_out

end

function (mcmc::SisIexInsertDelete)(
    posterior::SisPosterior,
    γ_fixed::Float64;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::Vector{Path{Int}}=sample_frechet_mean(posterior.data, posterior.dist),
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    )
    sample_out = Vector{Vector{Path{Int}}}(undef, desired_samples)

    (
        edit_count, edit_acc_count,
        flip_count, flip_acc_count,
        trans_dim_count, trans_dim_acc_count,
        suff_stats
        ) = draw_sample_mode!(
            sample_out,
            mcmc,
            posterior, γ_fixed,
            burn_in=burn_in,
            lag=lag,
            S_init=S_init,
            loading_bar=loading_bar,
            aux_init_at_prev=aux_init_at_prev
            )

    p_measures = Dict(
            "Proportion Update Moves" => (edit_count+flip_count)/(edit_count+flip_count+trans_dim_count),
            "Edit Alloc Move Acceptance Probability" => edit_acc_count / edit_count,
            "Flip Alloc Acceptance Probability" => flip_acc_count / flip_count,
            "Trans-Dimensional Move Acceptance Probability" => trans_dim_acc_count / trans_dim_count
        )
    output = SisPosteriorModeConditionalMcmcOutput(
            γ_fixed,
            sample_out,
            posterior,
            suff_stats,
            p_measures
            )

    return output

end

# Dispersion Conditional
# ----------------------

function draw_sample_gamma!(
    sample_out::Union{Vector{Float64}, SubArray},
    mcmc::SisIexInsertDelete,
    posterior::SisPosterior,
    S_fixed::InteractionSequence{Int};
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    γ_init::Float64=4.0,
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    )

    if loading_bar
        iter = Progress(
            length(sample_out), # How many iters
            1,  # At which granularity to update loading bar
            "Chain for n = $(posterior.sample_size) (dispersion conditional)....")  # Loading bar. Minimum update interval: 1 second
    end

    # Define aliases for pointers to the storage of current vals and proposals
    ε = mcmc.ε
    aux_mcmc = mcmc.aux_mcmc

    acc_count = 0
    i = 1 # Which iteration we are on
    sample_count = 1  # Which sample we are working to get

    S_curr = deepcopy(S_fixed)
    γ_curr = γ_init
    aux_data = [[Int[]] for i in 1:posterior.sample_size]

    # Evaluate sufficient statistic
    suff_stat = sum(
        x -> posterior.dist(S_curr, x),
        posterior.data
        )
    # Initialise the aux_data
    aux_model = SIS(
        S_curr, γ_curr,
        posterior.dist,
        posterior.V,
        posterior.K_inner,
        posterior.K_outer)
    draw_sample!(aux_data, aux_mcmc, aux_model)

    while sample_count ≤ length(sample_out)
        # Store value
        if (i > burn_in) & (((i-1) % lag)==0)
            sample_out[sample_count] = γ_curr
            sample_count += 1
        end

        γ_prop = rand_reflect(γ_curr, ε, 0.0, Inf)

        aux_model = SIS(
            S_curr, γ_prop,
            posterior.dist,
            posterior.V,
            posterior.K_inner, posterior.K_outer
            )
        if aux_init_at_prev
            tmp = deepcopy(aux_data[end])
            draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
        else
            draw_sample!(aux_data, aux_mcmc, aux_model)
        end

        # Accept reject

        log_lik_ratio = (γ_curr - γ_prop) * suff_stat

        dist = posterior.dist
        aux_log_lik_ratio = (γ_prop - γ_curr) * sum(x->dist(x,S_curr),aux_data)

        log_α = (
            logpdf(posterior.γ_prior, γ_prop)
            - logpdf(posterior.γ_prior, γ_curr)
            + log_lik_ratio + aux_log_lik_ratio
        )
        if log(rand()) < log_α
            γ_curr = γ_prop
            acc_count += 1
        end
        if loading_bar
            next!(iter)
        end
        i += 1

    end
    return acc_count
end

function draw_sample_gamma(
    mcmc::SisIexInsertDelete,
    posterior::SisPosterior,
    S_fixed::InteractionSequence{Int};
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    γ_init::Float64,
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    )

    sample_out = Vector{Float64}(undef, desired_samples)
    draw_sample_gamme!(
        sample_out,
        mcmc, posterior,
        S_fixed,
        burn_in=burn_in, lag=lag, γ_init=γ_init,
        loading_bar=loading_bar,
        aux_init_at_prev=aux_init_at_prev
        )
    return sample_out

end

function (mcmc::SisIexInsertDelete)(
    posterior::SisPosterior,
    S_fixed::InteractionSequence{Int};
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    γ_init::Float64=5.0,
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    )

    sample_out = Vector{Float64}(undef, desired_samples)

    acc_count = draw_sample_gamma!(
            sample_out,
            mcmc,
            posterior, S_fixed,
            burn_in=burn_in,
            lag=lag,
            γ_init=γ_init,
            loading_bar=loading_bar,
            aux_init_at_prev=aux_init_at_prev
            )

    p_measures = Dict(
            "Acceptance Probability" => acc_count/desired_samples
        )

    output = SisPosteriorDispersionConditionalMcmcOutput(
            S_fixed,
            sample_out,
            posterior,
            p_measures
            )

    return output

end

# Joint Distribution
# ------------------

# Mode accept-reject

function accept_reject_mode!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SisPosterior,
    γ_curr::Float64,
    mcmc::SisIexInsertDelete,
    P::CumCondProbMatrix,
    p_ins::Categorical,
    aux_data::InteractionSequenceSample{Int},
    acc_count::Vector{Int},
    count::Vector{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool
    )

    β = mcmc.β
    if rand() < β
        if rand() < 0.5
            was_accepted, suff_stat_curr = double_iex_multinomial_edit_accept_reject!(
                S_curr, S_prop,
                posterior, γ_curr,
                mcmc, P,
                aux_data,
                suff_stat_curr,
                aux_init_at_prev
            )
            acc_count[1] += was_accepted
            count[1] += 1
        else
            was_accepted, suff_stat_curr = double_iex_flip_accept_reject!(
                S_curr, S_prop,
                posterior, γ_curr,
                mcmc, P,
                aux_data,
                suff_stat_curr,
                aux_init_at_prev
            )
            acc_count[2] += was_accepted
            count[2] += 1
        end

    else
        was_accepted, suff_stat_curr = double_iex_trans_dim_informed_accept_reject!(
            S_curr, S_prop,
            posterior,
            γ_curr,
            mcmc,
            p_ins,
            aux_data,
            suff_stat_curr,
            aux_init_at_prev
        )

        acc_count[3] += was_accepted
        count[3] += 1
    end
    return suff_stat_curr
end

function accept_reject_gamma!(
    γ_curr::Float64,
    S_curr::InteractionSequence{Int},
    posterior::SisPosterior,
    mcmc::SisIexInsertDelete,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool
    )

    ε = mcmc.ε
    aux_mcmc = mcmc.aux_mcmc

    γ_prop = rand_reflect(γ_curr, ε, 0.0, Inf)

    aux_model = SIS(
        S_curr, γ_prop,
        posterior.dist,
        posterior.V,
        posterior.K_inner, posterior.K_outer
        )
    if aux_init_at_prev
        tmp = deepcopy(aux_data[end])
        draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux_data, aux_mcmc, aux_model)
    end

    # Accept reject

    log_lik_ratio = (γ_curr - γ_prop) * suff_stat_curr
    dist = posterior.dist
    aux_log_lik_ratio = (γ_prop - γ_curr) * sum(x->dist(x,S_curr),aux_data)

    log_α = (
        logpdf(posterior.γ_prior, γ_prop)
        - logpdf(posterior.γ_prior, γ_curr)
        + log_lik_ratio + aux_log_lik_ratio
    )
    if log(rand()) < log_α
        return γ_prop, 1
    else
        return γ_curr, 0
    end
end

function draw_sample!(
    sample_out_S::Union{InteractionSequenceSample{Int},SubArray},
    sample_out_gamma::Union{Vector{Float64},SubArray},
    mcmc::SisIexInsertDelete,
    posterior::SisPosterior;
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::Vector{Path{Int}}=sample_frechet_mean(posterior.data, posterior.dist),
    γ_init::Float64=5.0,
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    )

    if loading_bar
        iter = Progress(
            length(sample_out_S), # How many iters
            1,  # At which granularity to update loading bar
            "Chain for n = $(posterior.sample_size) (joint)....")  # Loading bar. Minimum update interval: 1 second
    end

    # Define aliases for pointers to the storage of current vals and proposals
    curr_pointers = mcmc.curr_pointers
    prop_pointers = mcmc.prop_pointers
    aux_mcmc = mcmc.aux_mcmc

    S_curr = InteractionSequence{Int}()
    S_prop = InteractionSequence{Int}()
    for i in 1:length(S_init)
        migrate!(S_curr, curr_pointers, i, 1)
        migrate!(S_prop, prop_pointers, i, 1)
        copy!(S_curr[i], S_init[i])
        copy!(S_prop[i], S_init[i])
    end
    γ_curr = copy(γ_init)

    sample_count = 1 # Keeps which sample to be stored we are working to get
    i = 1 # Keeps track all samples (included lags and burn_ins)

    acc_count = [0,0,0]
    count = [0,0,0]
    γ_acc_count = 0

    aux_data = [[Int[]] for i in 1:posterior.sample_size]
    # Initialise the aux_data
    aux_model = SIS(
        S_curr, γ_curr,
        posterior.dist,
        posterior.V,
        posterior.K_inner,
        posterior.K_outer)
    draw_sample!(aux_data, aux_mcmc, aux_model)
    # Initialise sufficient statistic
    suff_stat_curr = mapreduce(
        x -> posterior.dist(S_curr, x),
        +,
        posterior.data
    )
    suff_stats = Float64[suff_stat_curr]
    # Get informed proposal matrix
    P, vmap, vmap_inv = get_informed_proposal_matrix(posterior, mcmc.α)
    p_ins = get_informed_insertion_dist(posterior, mcmc.α)
    while sample_count ≤ length(sample_out_S)
        # Store values
        if (i > burn_in) & (((i-1) % lag)==0)
            sample_out_S[sample_count] = deepcopy(S_curr)
            sample_out_gamma[sample_count] = copy(γ_curr)
            push!(suff_stats, suff_stat_curr)
            sample_count += 1
        end

        # Update mode
        # -----------
        suff_stat_curr = accept_reject_mode!(
            S_curr, S_prop,
            posterior, γ_curr,
            mcmc, P, p_ins,
            aux_data,
            acc_count, count,
            suff_stat_curr,
            aux_init_at_prev
        )
        # Update gamma
        # ------------
        γ_curr, tmp =  accept_reject_gamma!(
            γ_curr,
            S_curr,
            posterior,
            mcmc,
            aux_data,
            suff_stat_curr,
            aux_init_at_prev
        )
        γ_acc_count += tmp
        if loading_bar
            next!(iter)
        end
        i += 1
    end

    for i in 1:length(S_curr)
        migrate!(curr_pointers, S_curr, 1, 1)
        migrate!(prop_pointers, S_prop, 1, 1)
    end

    ed_acc_prob = acc_count[1]/count[1]
    flip_acc_prob = acc_count[2]/count[2]
    td_acc_prob = acc_count[3]/count[3]
    γ_acc_prob = γ_acc_count / sum(count)
    return ed_acc_prob, flip_acc_prob, td_acc_prob, γ_acc_prob, suff_stats
end

function draw_sample(
    mcmc::SisIexInsertDelete,
    posterior::SisPosterior;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::InteractionSequence{Int}=sample_frechet_mean(posterior.data, posterior.dist),
    γ_init::Float64=5.0,
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    )

    sample_out_S = Vector{InteractionSequence}(undef, desired_samples)
    sample_out_gamma = Vector{Float64}(undef, desired_samples)

    draw_sample!(
        sample_out_S,
        sample_out_gamma,
        mcmc,
        posterior,
        burn_in=burn_in, lag=lag,
        S_init=S_init, γ_init=γ_init,
        loading_bar=loading_bar,
        aux_init_at_prev=aux_init_at_prev
    )

    return (S=sample_out_S, gamma=sample_out_gamma)
end

function (mcmc::SisIexInsertDelete)(
    posterior::SisPosterior;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::InteractionSequence{Int}=sample_frechet_mean(posterior.data, posterior.dist),
    γ_init::Float64=5.0,
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    )

    sample_out_S = Vector{InteractionSequence{Int}}(undef, desired_samples)
    sample_out_gamma = Vector{Float64}(undef, desired_samples)

    ed_acc_prob, flip_acc_prob, td_acc_prob, γ_acc_prob, suff_stats = draw_sample!(
        sample_out_S,
        sample_out_gamma,
        mcmc,
        posterior,
        burn_in=burn_in, lag=lag,
        S_init=S_init, γ_init=γ_init,
        loading_bar=loading_bar,
        aux_init_at_prev=aux_init_at_prev
    )

    p_measures = Dict(
            "Dipsersion acceptance probability" => γ_acc_prob,
            "Edit allocation acceptance probability" => ed_acc_prob,
            "Flip allocation acceptance probability" => flip_acc_prob,
            "Trans-dimensional acceptance probability" => td_acc_prob
        )

    return SisPosteriorMcmcOutput(
        sample_out_S,
        sample_out_gamma,
        posterior,
        suff_stats,
        p_measures
    )

end


text:
import Statistics

lines = readlines("input.txt")
matrix = parse.(Bool, hcat(collect.(lines)...))

oxygen = co2 = matrix'
for idx = axes(matrix, 1)
    if size(oxygen, 1) > 1
        mean = Statistics.mean(oxygen[:,idx]) .>= 0.5
        global oxygen = oxygen[findall(x -> x[idx] == mean, eachrow(oxygen)), :]
    end

    if size(co2, 1) > 1
        mean = Statistics.mean(co2[:,idx]) .>= 0.5
        global co2 = co2[findall(x -> x[idx] != mean, eachrow(co2)), :]
    end
end

function bitarr_to_int(arr)
    return sum(arr .* (2 .^ collect(length(arr)-1:-1:0)))
end

println(bitarr_to_int(oxygen') * bitarr_to_int(co2'))

ext: jl


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Types and structures.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    struct PrintInfo{Td,Th}

This structure stores the information required so that the backends can print
the tables.

"""
@with_kw struct PrintInfo{Td,Th}
    data::Td
    header::Th
    id_cols::Vector{Int}
    id_rows::Vector{Int}
    num_rows::Int
    num_cols::Int
    num_printed_cols::Int
    num_printed_rows::Int
    header_num_rows::Int
    header_num_cols::Int
    alignment::Vector{Symbol}
end


text:
################################################################################
### Type and general utility
################################################################################

"""
    BlockDiagonal(blocks...)
    BlockDiagonal{eltype, N_blocks, MatrixType}(init, n, m)

Creates a block diagonal matrix from a collection of blocks or an initializer,
i.e. `undef` or `I`. In the latter case `n, m` specifies the size of the full
block diagonal matrix.

Notes:
- `vmul!` and friends are implemented for real and complex `BlockDiagonal`
- some base methods are implemented - `Matrix`, `getidnex`, `setindex`, `copyto!`,
`size`, `copy`, `log`, `exp`, `det`, `*`
- BlockDiagonal of 2 equally sized square matrices is heavily tested and used
- BlockDiagonal of more equally sized square matrices should work
- irregularly sized and non-square matrices are not supported
- methods generally assume correct size and use @inbounds for performance.
"""
struct BlockDiagonal{T, N, AT <: AbstractMatrix{T}} <: AbstractMatrix{T}
    blocks::NTuple{N, AT}

    function BlockDiagonal(blocks::NTuple{N, AT}) where {T, N, AT <: AbstractMatrix{T}}
        n, m = size(first(blocks))
        n == m || throw(ErrorException("Each block must be square."))
        all(b -> size(b) == (n, n), blocks) || throw(ErrorException(
            "All Blocks must have the same size."
        ))
        # One matrix with zeros is faster than multiple 3x3 matrices on my
        # machine (avx2). This will most likely depend on AVX version...
        # if n < 4
        #     full = zeros(T, n*N, n*N)
        #     for i in 1:N
        #         @views copyto!(full[(i-1)*n+1 : i*n, (i-1)*n+1 : i*n], blocks[i])
        #     end
        #     # new{T, 1, AT}((full,))
        #     # Is this bad? - Yes, because the rest of the code doesn't know about Matrix returns
        #     return full
        # else
            return new{T, N, AT}(blocks)
        # end
    end
end

BlockDiagonal(blocks::AT...) where {T, AT <: AbstractMatrix{T}} = BlockDiagonal(blocks)
function BlockDiagonal{T, N, AT}(::UndefInitializer, n, m) where {T, N, AT}
    n == m || throw(ErrorException("Each block must be square."))
    (n, m) = divrem(n, N)
    m == 0 || throw(ErrorException(
        "The size of the Matrix is incompatible with the blocking ($n % $N != 0)"
    ))
    BlockDiagonal([AT(undef, n, n) for _ in 1:N]...)
end
function BlockDiagonal{T, N, AT}(::UniformScaling, n, m) where {T, N, AT}
    n == m || throw(ErrorException("Each block must be square."))
    (n, m) = divrem(n, N)
    m == 0 || throw(ErrorException(
        "The size of the Matrix is incompatible with the blocking ($n % $N != 0)"
    ))
    BlockDiagonal([AT(I, n, n) for _ in 1:N]...)
end

# Shorthand for typing functions
const BD{N} = BlockDiagonal{<: Number, N}

function Base.Matrix(B::BlockDiagonal{T, N}) where {T, N}
    n = size(B.blocks[1], 1)
    output = zeros(T, N * n, N * n)
    copyto!(output, B)
    output
end

function Base.checkbounds(B::BlockDiagonal, i::Integer, j::Integer)
    N = M = 0
    @inbounds for b in B.blocks
        N += size(b, 1)
        M += size(b, 2)
    end
    (1 ≤ i ≤ N) && (1 ≤ j ≤ M)
end

function Base.getindex(B::BlockDiagonal{T, N}, i, j) where {T, N}
    @boundscheck checkbounds(B, i, j)
    @inbounds n = size(B.blocks[1], 1)

    b, k = divrem(i-1, n)
    _b, l = divrem(j-1, n)

    if b != _b
        return zero(T)
    else
        @inbounds return B.blocks[b+1][k+1, l+1]
    end
end

function Base.setindex!(B::BlockDiagonal, val, i, j)
    i0 = j0 = 0
    for block in B.blocks
        i1 = size(block, 1) + i0
        j1 = size(block, 2) + j0
        if (i0 < i <= i1) && (j0 < j <= j1)
            @inbounds block[i - i0, j - j0] = val
            return B
        end
        i0 = i1; j0 = j1
    end
    throw(BoundsError(B, (i, j)))
end

function Base.copyto!(B1::BD{N}, B2::BD{N}) where {N}
    @inbounds for i in 1:N
        copyto!(B1.blocks[i], B2.blocks[i])
    end
    return B1
end
function Base.copyto!(B::BlockDiagonal, ::UniformScaling)
    @inbounds for block in B.blocks
        copyto!(block, I)
    end
    return B
end
function Base.copyto!(B::BlockDiagonal, D::Diagonal)
    n = 0
    @inbounds for block in B.blocks
        s = 1 + n
        n += size(block, 1)
        @views copyto!(block, Diagonal(D.diag[s:n]))
    end
    return B
end

function Base.copyto!(output::Matrix, B::BlockDiagonal{T, N}) where {T, N}
    n = size(B.blocks[1], 1)
    for i in 1:N
        @views copyto!(output[(i-1)*n+1 : i*n, (i-1)*n+1 : i*n], B.blocks[i])
    end
    output
end

# for Base.show
Base.size(B::BlockDiagonal) = mapreduce(b -> size(b), (a, b) -> a .+ b, B.blocks)

# Needed for tests
Base.copy(B::BlockDiagonal) = deepcopy(B)

# Needed for stack building (shouldn't be performance critical)
function Base.:(*)(x::Number, B::BlockDiagonal)
    BlockDiagonal(map(block -> x * block, B.blocks)...)
end
Base.:(*)(B::BlockDiagonal, x::Number) = x * B
function Base.:(*)(B1::BlockDiagonal{T1, N}, B2::BlockDiagonal{T2, N}) where {T1, T2, N}
    BlockDiagonal(map(*, B1.blocks, B2.blocks)...)
end

function Base.exp(B::BlockDiagonal)
    BlockDiagonal(map(block -> exp(block), B.blocks))
end

# This takes super long to compile...?
# Dot syntax is another 6x slower
# function Base.log(B::BlockDiagonal)
#     BlockDiagonal(map(block -> log(block), B.blocks))
# end

# function LinearAlgebra.det(B::BlockDiagonal{T}) where {T}
#     output = T(1)
#     for b in B.blocks
#         output *= det(b)
#     end
#     output
# end

# I thought this would be needed for greens(k, l), but it's not?
function LinearAlgebra.transpose!(A::BlockDiagonal{T, N}, B::BlockDiagonal{T, N}) where {T, N}
    @inbounds for i in 1:N
        transpose!(A.blocks[i], B.blocks[i])
    end
    A
end
# This however is used for greens(k, l)
function LinearAlgebra.rmul!(B::BlockDiagonal, f::Number)
    @inbounds for block in B.blocks
        rmul!(block, f)
    end
    B
end

################################################################################
### AVX-powered matrix multiplications
################################################################################

function vmul!(C::BD{N}, A::BD{N}, B::BD{N}) where {N}
    @inbounds for i in 1:N
        vmul!(C.blocks[i], A.blocks[i], B.blocks[i])
    end
    nothing
end
function vmul!(C::BD{N}, A::BD{N}, B::Diagonal) where {N}
    n = size(C.blocks[1], 1)
    @inbounds for i in 1:N
        vmul!(C.blocks[i], A.blocks[i], B, (i-1)*n+1 : i*n)
    end
    nothing
end
function vmul!(C::BD{N}, A::Diagonal, B::BD{N}) where {N}
    # Assuming correct size
    n = size(C.blocks[1], 1)
    @inbounds for i in 1:N
        vmul!(C.blocks[i], A, B.blocks[i], (i-1)*n+1 : i*n)
    end
    nothing
end
function vmul!(C::BD{N}, A::BD{N}, B::Adjoint) where {N}
    @inbounds for i in 1:N
        vmul!(C.blocks[i], A.blocks[i], Adjoint(B.parent.blocks[i]))
    end
end
function vmul!(C::BD{N}, A::Adjoint, B::BD{N}) where {N}
    @inbounds for i in 1:N
        vmul!(C.blocks[i], Adjoint(A.parent.blocks[i]), B.blocks[i])
    end
end
function vmul!(C::BD{N}, A::Adjoint, B::Adjoint) where {N}
    @inbounds for i in 1:N
        vmul!(C.blocks[i], Adjoint(A.parent.blocks[i]), Adjoint(B.parent.blocks[i]))
    end
end
function rvmul!(A::BD{N}, B::Diagonal) where {N}
    n = size(A.blocks[1], 1)
    @inbounds for i in 1:N
        @views rvmul!(A.blocks[i], Diagonal(B.diag[(i-1)*n+1 : i*n]))
    end
end
function lvmul!(A::Diagonal, B::BD{N}) where {N}
    n = size(B.blocks[1], 1)
    @inbounds for i in 1:N
        @views lvmul!(Diagonal(A.diag[(i-1)*n+1 : i*n]), B.blocks[i])
    end
end

# used in greens(k, l)
function rvadd!(A::BD{N}, B::BD{N}) where {N}
    @inbounds for i in 1:N
        rvadd!(A.blocks[i], B.blocks[i])
    end
end
# used in equal time greens
function rvadd!(B::BD{N}, D::Diagonal) where {N}
    # Assuming correct size
    n = size(B.blocks[1], 1)
    @inbounds for i in 1:N
        @views rvadd!(B.blocks[i], Diagonal(D.diag[(i-1)*n+1 : i*n]))
    end
end
# used in CombinedGreensIterator
function vsub!(O::BD{N}, A::BD{N}, ::UniformScaling) where {N}
    @inbounds for i in 1:N
        vsub!(O.blocks[i], A.blocks[i], I)
    end
end

function rdivp!(A::BD, T::BD, O::BD, pivot) where {ET<:Real, N, BD <: BlockDiagonal{ET, N}}
    # assume Diagonal is ±1!
    @inbounds begin
        n = size(A.blocks[1], 1)

        # Apply pivot
        for b in 1:N
            o = O.blocks[b]
            a = A.blocks[b]
            t = T.blocks[b]
            offset = (b-1)*n

            for j in 1:n
                p = pivot[j + offset]
                @turbo for i in 1:n
                    o[i, j] = a[i, p]
                end
            end

            # do the rdiv
            # @turbo will segfault on `k in 1:0`, so pull out first loop
            invt11 = 1.0 / t[1, 1]
            @turbo for i in 1:n
                a[i, 1] = o[i, 1] * invt11
            end
            for j in 2:n
                invtjj = 1.0 / t[j, j]
                @turbo for i in 1:n
                    x = o[i, j]
                    for k in 1:j-1
                        x -= a[i, k] * t[k, j]
                    end
                    a[i, j] = x * invtjj
                end
            end
        end
    end
    A
end

################################################################################
### Fallbacks
################################################################################

function rdivp!(A::BD, T::BD, O::BD, pivot) where {ET, N, BD <: BlockDiagonal{ET, N}}
    # Assuming correct size
    n = size(A.blocks[1], 1)
    @inbounds for i in 1:N
        @views rdivp!(A.blocks[i], T.blocks[i], O.blocks[i], pivot[(i-1)*n+1 : i*n])
    end
    A
end

################################################################################
### UDT
################################################################################

function udt_AVX_pivot!(
        U::BlockDiagonal{Float64, N, AT},
        D::AbstractArray{Float64, 1},
        input::BlockDiagonal{Float64, N, AT},
        pivot::AbstractArray{Int64, 1} = Vector(UnitRange(1:size(input, 1))),
        temp::AbstractArray{Float64, 1} = Vector{Float64}(undef, length(D)),
        apply_pivot::Val = Val(true)
    ) where {N, AT <: AbstractMatrix{Float64}}
    # Assuming correct size
    n = 0
    @inbounds for i in 1:N
        s = 1+n
        n += size(U.blocks[i], 1)
        @views udt_AVX_pivot!(
            U.blocks[i], D[s:n], input.blocks[i],
            pivot[s:n], temp[s:n], apply_pivot
        )
    end
end
function udt_AVX_pivot!(
        U::BlockDiagonal{ComplexF64, N, AT},
        D::AbstractArray{Float64, 1},
        input::BlockDiagonal{ComplexF64, N, AT},
        pivot::AbstractArray{Int64, 1} = Vector(UnitRange(1:size(input, 1))),
        temp::AbstractArray{ComplexF64, 1} = Vector{ComplexF64}(undef, length(D)),
        apply_pivot::Val = Val(true)
    ) where {N, AT <: AbstractMatrix{ComplexF64}}
    # Assuming correct size
    n = 0
    @inbounds for i in 1:N
        s = 1+n
        n += size(U.blocks[i], 1)
        @views udt_AVX_pivot!(
            U.blocks[i], D[s:n], input.blocks[i],
            pivot[s:n], temp[s:n], apply_pivot
        )
    end
end


text:
###
#   An implementation of Transitional Markov Chain Monte Carlo in Julia
#
#            Institute for Risk and Uncertainty, Uni of Liverpool
#
#                       Authors: Ander Gray, Adolphus Lye
#
#                       Email: Ander.Gray@liverpool.ac.uk,
#                              Adolphus.Lye@liverpool.ac.uk
#
#
#   This Transitional MCMC algorithm is inspirted by OpenCOSSAN:
#                   https://github.com/cossan-working-group/OpenCossan
#
#
#   Algorithm originally proposed by:
#           J. Ching, and Y. Chen (2007). Transitional Markov Chain Monte Carlo method
#           for Bayesian model updating, Model class selection, and Model averaging.
#           Journal of Engineering Mechanics, 133(7), 816-832.
#           doi:10.1061/(asce)0733-9399(2007)133:7(816)
#
###

function tmcmc(log_fD_T, fT, sample_fT, Nsamples, burnin=20, thin=3, beta2=0.01)

    j1 = 0;                     # Iteration number
    βj = 0;                     # Tempering parameter
    θ_j = sample_fT(Nsamples);  # Samples of prior
    Lp_j = zeros(Nsamples, 1);   # Log liklihood of first iteration

    Log_ev = 0                  # Log Evidence

    Ndims = size(θ_j, 2)         # Number of dimensions (input)

    while βj < 1

        j1 = j1 + 1
        println()
        println("Beginning iteration $j1")

        ###
        # Parallel evaluation of the likelihood
        ###
        print("Computing likelihood with $(nworkers()) workers....")
        # Lp_j = log_fD_T(θ_j')
        ins = [θ_j[i,:] for i in 1:size(θ_j, 1)]
        Lp_j = pmap(log_fD_T, ins)
        Lp_j = reduce(vcat, Lp_j)
        # Lp_j = pamp(log_fD_T,θ_j')
        println("Done!")

        ###
        #   Computing new βj
        #   Uses bisection method
        ###
        println("Computing Bj")

        low_β = βj; hi_β = 2; Lp_adjust = maximum(Lp_j);
        x1 = (hi_β + low_β) / 2;

        while (hi_β - low_β) / ((hi_β + low_β) / 2) > 1e-6
            x1 = (hi_β + low_β) / 2;
            wj_test = exp.((x1 .- βj ) .* (Lp_j .- Lp_adjust));
            cov_w   = std(wj_test) / mean(wj_test);
            if cov_w > 1; hi_β = x1; else; low_β = x1; end
        end

        βj1 = min(1, x1)
        println("B_$(j1) = $(βj1)")

        ###
        #   Computation of normalised weights
        ###
        println("Computing weights")

        w_j = exp.((βj1 - βj) .* (Lp_j .- Lp_adjust))       # Nominal weights from likilhood and βjs

        Log_ev = log(mean(w_j)) + (βj1 - βj) * Lp_adjust + Log_ev   # Log evidence in current iteration

        # Normalised weights
        wn_j = w_j ./ sum(w_j);

        Th_wm = θ_j .* wn_j                 # Weighted mean of samples

        ###
        #   Calculation of COV matrix of proposal
        ###
        SIGMA_j = zeros(Ndims, Ndims)

        for l = 1:Nsamples
            SIGMA_j = SIGMA_j + beta2 .* wn_j[l] .* (θ_j[l,:]' .- Th_wm)' * (θ_j[l,:]' .- Th_wm)
        end

        # Ensure that cov is symetric
        SIGMA_j = (SIGMA_j' + SIGMA_j) / 2

        prop = mu -> proprnd(mu, SIGMA_j, fT)           # Anonymous function for proposal

        target = x -> log_fD_T(x) .* βj1 .+ log.(fT(x)) # Anonymous function for transitional distribution

        # Weighted resampling of θj (indecies with replacement)
        randIndex = sample(1:Nsamples, Weights(wn_j), Nsamples, replace=true)

        ins = [θ_j[randIndex[i], :] for i = 1:Nsamples]

        print("Markov chains with $(nworkers()) workers...")
        # pmap is a parallel map
        θ_j1 = pmap(x -> runChains(target, prop, x, burnin, thin), ins)

        θ_j1 = reduce(vcat, θ_j1)

        println("Done!")

        # println("Mean α = $(meanα)")

        βj = βj1
        θ_j = θ_j1
    end
    return θ_j, Log_ev
end

function proprnd(mu, covMat, prior)

    samp = rand(MvNormal(mu, covMat), 1)
    while iszero(prior(samp))
        samp = rand(MvNormal(mu, covMat),1)
    end
    return samp[:]
end

function runChains(target, prop, θ_js, burnin, thin)
    samps, α  = MHsampleSimple(target, prop, θ_js, 1, burnin, thin)
    return samps
end

function runChains2(target, prop, θ_js, burnin, thin)
    Nsamples = size(θ_js, 1)
    Ndims    = size(θ_js, 2)
    θ_j1 = zeros(Nsamples, Ndims)
    α = zeros(Nsamples, 1)
    for i = 1:Nsamples
        samps, αs = MHsampleSimple(target, prop, θ_j[i], 1, burnin, thin)
        θ_j1[i,:] = samps
        α[i] = αs
    end

    return θ_j1, α

end


text:
close all
n=10;
Param.nPanel=n;
Param.Grav=9.81;
Param.Omega=2*pi/(24*3600);
Param.RadEarth=6.37122d+6;
Param.Coriolis=true;
Param.HyperDM=0.0; %1.e15*(n/30)^3.2;
Param.HyperD=1.e15*(n/30)^3.2;
Param.Flat=true;

fig=1;
level=1;
OrdPoly=4;
CG.OrdPoly=OrdPoly;
Param.Grid=CubedGrid(n,OrientFaceSphere,Param);
nz=2;
z=zeros(nz,1);
dz=500;
z(1)=dz/2;
for i=2:nz
  z(i)=z(i-1)+dz;
end
Param.Grid.nz=nz;
Param.Grid.z=z;
Param.Grid.dz=dz;
[CG.Faces,CG.NumG,CG.NumI,CG.Glob,CG.FaceGlob,CG.Stencil]...
  =NumberingFemCG(Param.Grid,OrdPoly);
Param.J=zeros(OrdPoly+1,OrdPoly+1,Param.Grid.NumFaces,nz);
Param.X=zeros(OrdPoly+1,OrdPoly+1,3,Param.Grid.NumFaces,nz);
Param.dXdx=zeros(OrdPoly+1,OrdPoly+1,Param.Grid.NumFaces,nz,3,3);
Param.lat=zeros(OrdPoly+1,OrdPoly+1,Param.Grid.NumFaces);
[CG.w,CG.xw]=GaussLobattoQuad(CG.OrdPoly);
[CG.wX,CG.xwX]=GaussLobattoQuad(CG.OrdPoly);
[CG.wY,CG.xwY]=GaussLobattoQuad(CG.OrdPoly);
[CG.DW,CG.DS]=DerivativeMatrixSingle(CG);

for iF=1:Param.Grid.NumFaces
  for iz=1:nz
  for j=1:OrdPoly+1
    for i=1:OrdPoly+1
      [Param.X(i,j,:,iF),Param.J(i,j,iF,iz),Param.dXdx(i,j,iF,iz,:,:),Param.lat(i,j,iF)]...
        =JacobiSphere3(CG.xw(i),CG.xw(j),Param.Grid.Faces(iF),z(iz),Param.Grid);
    end
  end
  end
end
[CG.M]=MassCG(CG,Param);
U=zeros(CG.NumG,nz,2);
W=zeros(CG.NumG,nz+1,1);
Profh='Galewsky';
ProfVel='Galewsky';
Param.ProfVel=ProfVel;
Param.Profh=Profh;
Rho=Project(@fh,CG,@JacobiSphere3,Param);
[U(:,:,2),U(:,:,3)]=ProjectVec(@fVel,CG,@JacobiSphere3,Param);
Param.hS=zeros(size(U(:,:,1)));
Th=Rho;
CFL=0.125;
dtau=100;
time=0;

IntMethod='RungeKutta';
nIter=200;
%v = VideoWriter ('Galewsky.avi');
%open (v);
switch IntMethod
  case 'RungeKutta'
    for i=1:nIter
      [FRho,FU,FW,FTh]=FcnNHCurlVec(Rho,U,W,Th,CG,Param);
      RhoNeu=Rho+1/3*dtau*FRho;
      UNeu=U+1/3*dtau*FU;
      WNeu=W+1/3*dtau*FW;
      ThNeu=Th+1/3*dtau*FTh;
      F=FcnShallowCurlVec(UNeu,CG,Param);
      UNeu=U+1/2*dtau*F;
      F=FcnShallowCurlVec(UNeu,CG,Param);
      U=U+dtau*F;
      time=time+dtau;
      if mod(i,200)==0
        Vort=FCurl2(U(:,:,2),U(:,:,3),CG,Param);
        PlotCG(Vort,CG,@JacobiSphere3,Param,fig,level);
%         frame=getframe(gcf);
%         writeVideo (v, frame);
      end
    end
end
%fig=PlotCG(U(:,1),CG,@JacobiSphere2,Param,fig);
%close(v)

for i=1:200
GradVec=FGrad2Vec(U(:,1),CG,Param);
Grad=FGrad2(U(:,1),CG,Param);
DivVec=FDiv2Vec(U(:,2),U(:,3),CG,Param);
Div=FDiv2(U(:,2),U(:,3),CG,Param);
FVec=FCurlNon2Vec(U(:,2),U(:,3),CG,Param);
F=FCurlNon2(U(:,2),U(:,3),CG,Param);
end
bb=3;


text:
module TestMultiobjective

using Gurobi
using Test

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

const MOI = Gurobi.MOI

const GRB_ENV = isdefined(Main, :GRB_ENV) ? Main.GRB_ENV : Gurobi.Env()

function test_multiobjective()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    MOI.Utilities.loadfromstring!(
        model,
        """
variables: x, y
minobjective: 2x + y
c1: x + y >= 1.0
c2: 0.5 * x + 1.0 * y >= 0.75
c3: x >= 0.0
c4: y >= 0.25
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    y = MOI.get(model, MOI.VariableIndex, "y")

    f = MOI.ScalarAffineFunction(
        [MOI.ScalarAffineTerm(1.0, x), MOI.ScalarAffineTerm(3.0, y)],
        0.0,
    )

    MOI.set(model, Gurobi.MultiObjectiveFunction(2), f)

    @test MOI.get(model, Gurobi.MultiObjectiveWeight(1)) == 1.0
    @test MOI.get(model, Gurobi.MultiObjectiveWeight(2)) == 1.0
    @test MOI.get(model, Gurobi.MultiObjectivePriority(1)) == 0
    @test MOI.get(model, Gurobi.MultiObjectivePriority(2)) == 0

    MOI.optimize!(model)
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1.5
    @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 0.5
    @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 0.5

    BFS = [
        (x = 1.0, y = 0.25, f1 = 2.25, f2 = 1.75),
        (x = 0.5, y = 0.5, f1 = 1.5, f2 = 2.0),
        (x = 0.0, y = 1.0, f1 = 1.0, f2 = 3.0),
    ]
    for (i, λ) in enumerate([0.2, 0.5, 0.8])
        MOI.set(model, Gurobi.MultiObjectiveWeight(1), λ)
        MOI.set(model, Gurobi.MultiObjectiveWeight(2), 1 - λ)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ BFS[i].x
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ BFS[i].y
        @test MOI.get(model, Gurobi.MultiObjectiveValue(1)) ≈ BFS[i].f1
        @test MOI.get(model, Gurobi.MultiObjectiveValue(2)) ≈ BFS[i].f2
    end

    MOI.set(model, Gurobi.MultiObjectiveWeight(1), 1.0)
    MOI.set(model, Gurobi.MultiObjectiveWeight(2), 1.0)
    MOI.set(model, Gurobi.MultiObjectivePriority(1), 1)
    MOI.set(model, Gurobi.MultiObjectivePriority(2), 2)
    @test MOI.get(model, Gurobi.MultiObjectivePriority(1)) == 1
    @test MOI.get(model, Gurobi.MultiObjectivePriority(2)) == 2

    MOI.optimize!(model)

    @test MOI.get(model, MOI.VariablePrimal(), x) ≈ BFS[1].x
    @test MOI.get(model, MOI.VariablePrimal(), y) ≈ BFS[1].y
    @test MOI.get(model, Gurobi.MultiObjectiveValue(1)) ≈ BFS[1].f1
    @test MOI.get(model, Gurobi.MultiObjectiveValue(2)) ≈ BFS[1].f2

    MOI.set(model, Gurobi.MultiObjectivePriority(1), 2)
    MOI.set(model, Gurobi.MultiObjectivePriority(2), 1)
    @test MOI.get(model, Gurobi.MultiObjectivePriority(1)) == 2
    @test MOI.get(model, Gurobi.MultiObjectivePriority(2)) == 1

    MOI.optimize!(model)

    @test MOI.get(model, MOI.VariablePrimal(), x) ≈ BFS[3].x
    @test MOI.get(model, MOI.VariablePrimal(), y) ≈ BFS[3].y
    @test MOI.get(model, Gurobi.MultiObjectiveValue(1)) ≈ BFS[3].f1
    @test MOI.get(model, Gurobi.MultiObjectiveValue(2)) ≈ BFS[3].f2
end

end

TestMultiobjective.runtests()


text:
using Glo
using Test

@testset "Glo.jl" begin
    @test true
    #TODO: how to do tests without sharing secrets?
end

##


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""
    BezierCurve{P,S,N,M} <: Spline{P,S,N,M}

`N` is the dimension of the curve, `M` is the curve order

```julia
BezierCurve{P,S,N,M}(controlpoints::AbstractArray{<:AbstractArray{S,1},1})
```
"""
struct BezierCurve{P,S,N,M} <: Spline{P,S,N,M} # this has a funny name because of name conflict with Plots.
    controlpolygon::Array{SVector{N,S},1}

    function BezierCurve{P,S,N,M}(controlpoints::AbstractArray{<:AbstractArray{S,1},1}) where {P<:CurveType,S<:Number,N,M}
        @assert length(controlpoints) == M + 1 && length(controlpoints[1]) == N
        return new{P,S,N,M}([SVector{N,S}(point) for point in controlpoints])
    end
end
export BezierCurve

# m is the order of the moving lines, n is the order of the curve, d is the spatial dimension of the curve
_arraysize(m, n, d) = (2m + n + 2, (m + 1) * (d + 1))
# computes the proper row index so control point Pj gets multiplied by correct non-redundant basis function.
# The m term is to offset all the Pj values below the first m+1 rows to allow for the correct summation of the final g(d+1,k) terms.
_rowindex(k, j, movinglineorder) = k + j + movinglineorder + 2

function movinglinearray(curve::BezierCurve{P,S,N,M}) where {P,S,N,M}
    # curveorder = M, spatialdimension = N
    movinglineorder = M - 1 # minimum order of moving line polynomial necessary to represent curve of order M.
    result = zeros(S, _arraysize(movinglineorder, M, N))
    ctrlpts = curve.controlpolygon

    # upper right corner will contain just Bc(k,m) coefficients that get multiplied by the moving line vector g[movinglineorder*(d+1) + k] term
    start = movinglineorder * (N + 1) + 1

    for index in start:(start + movinglineorder)
        k = index - start

        println("$k $movinglineorder $index")

        result[k + 1, index] = Bc(k, movinglineorder)
    end

    for spatialindex in 0:(N - 1)
        colblock = spatialindex * (movinglineorder + 1) + 1
        for j in 0:M
            for k in 0:movinglineorder
                row = _rowindex(k, j, movinglineorder)
                result[row, colblock + k] = Bc(j, M) * Bc(k, movinglineorder) * ctrlpts[j + 1][spatialindex + 1]
            end
        end
    end
    return result
end

###################################################################################################################################

"""
    BezierSurface{P,S,N,M} <: SplineSurface{P,S,N,M}

Bezier surface defined by grid of control points.

!!! danger
    This surface does not create a valid half-space, requires updates to function correctly.

```julia
BezierSurface{P,S,N,M}(controlpoints::AbstractArray{<:AbstractArray{S,1},2})
```
"""
struct BezierSurface{P,S,N,M} <: SplineSurface{P,S,N,M}
    controlpolygon::Array{SVector{N,S},2}

    function BezierSurface{P,S,N,M}(controlpoints::AbstractArray{<:AbstractArray{S,1},2}) where {P<:CurveType,S<:Real,N,M}
        @assert (M + 1, M + 1) == size(controlpoints) && length(controlpoints[1, 1]) == N
        return new([SVector{N,S}(point) for point in controlpoints])
    end

    # function BezierSurface{P,S,N,M}(controlpoints::Array{SVector{N,S},2}) where {P<:CurveType,S<:Number,N,M}
    #     @assert (M + 1, M + 1) == size(controlpoints)
    #     newpoints = [copy(point) for point in controlpoints]
    #     # temp1 = SVector{M + 1,SVector{N,S}}(undef, M + 1)
    #     # temp2 = SVector{M + 1,SVector{N,S}}(undef, M + 1)
    #     return new(newpoints)
    # end
end
export BezierSurface

# """this function transforms the control points of the Bezier patch before ray tracing. Probably a better idea to use transforms in CSG operators,
# although this could be a little less computation at run time."""
# function Base.:*(a::Transform{T}, surf::BezierSurface{P,T,N,M}) where {P,T,N,M}
#     result = deepcopy(surf)
#     for index in CartesianIndices(surf.controlpolygon)
#         result.controlpolygon[index] = a * result.controlpolygon[index]
#     end
#     return result
# end

deepcopy(a::BezierSurface{P,S,N,M}) where {P,S,N,M} = BezierSurface{P,S,N,M}(a.controlpolygon)
uvrange(::Type{BezierSurface{P,S}}) where {P,S<:Real} = ((zero(S), one(S)), (zero(S), one(S)))
uvrange(::BezierSurface{P,S}) where {P,S<:Real} = uvrange(BezierSurface{P,S})
onsurface(::BezierSurface{P,S,3,M}, ::SVector{3,S}) where {P,S<:Real,M} = false

# DECASTELJAU
@inline function decasteljau(controlpoints::AbstractArray{SVector{N,S},1}, ::Val{2}, u::T) where {T<:Real,S,N}
    w = 1 - u
    return (controlpoints[1] * w + controlpoints[2] * u) * w + (controlpoints[2] * w + controlpoints[3] * u) * u
end

@inline function decasteljau(controlpoints::AbstractArray{SVector{N,S},1}, ::Val{3}, u::T) where {T<:Real,S,N}
    w = 1 - u
    np1 = controlpoints[1] * w + controlpoints[2] * u
    np2 = controlpoints[2] * w + controlpoints[3] * u
    np3 = controlpoints[3] * w + controlpoints[4] * u
    np1 = np1 * w + np2 * u
    np2 = np2 * w + np3 * u
    np1 = np1 * w + np2 * u
    return np1
end

@inline function decasteljau(controlpoints::AbstractArray{SVector{N,S},1}, ::Val{4}, u::T) where {T<:Real,S,N}
    w = 1 - u
    np1 = controlpoints[1] * w + controlpoints[2] * u
    np2 = controlpoints[2] * w + controlpoints[3] * u
    np3 = controlpoints[3] * w + controlpoints[4] * u
    np4 = controlpoints[4] * w + controlpoints[5] * u
    np1 = np1 * w + np2 * u
    np2 = np2 * w + np3 * u
    np3 = np3 * w + np4 * u
    np1 = np1 * w + np2 * u
    np2 = np2 * w + np3 * u
    np1 = np1 * w + np2 * u
    return np1
end

@inline function decasteljau(controlpoints::AbstractArray{SVector{N,S},1}, ::Val{M}, u::T) where {T<:Real,S,N,M}
    # general solution for decasteljau for arbitrary curve order
    newpoints = MVector{M + 1,SVector{N,T}}(controlpoints)
    # do this to force newpoints to have the type of u rather than controlpoints. When u is a dual number this makes the array the correct type.
    @inbounds for i in 1:M
        @inbounds for k in 1:(M + 1 - i)
            newpoints[k] = newpoints[k] * (1 - u) + newpoints[k + 1] * u
            # this kills forwarddiff and zygote. ForwardDiff doesn't work because newpoints has the type of controlpoints, which is going to be Float64 most times,
            # when it needs to be of type dual number. Zygote doesn't support mutable arrays. Arrrrgggghhh! Useless!
        end
    end
    return newpoints[1]
end

# POINT
function point(surf::BezierSurface{P,S,N,3}, u::T, v::T) where {T<:Real,P,S,N}
    q1 = decasteljau(view(surf.controlpolygon, :, 1), Val(3), u)
    q2 = decasteljau(view(surf.controlpolygon, :, 2), Val(3), u)
    q3 = decasteljau(view(surf.controlpolygon, :, 3), Val(3), u)
    q4 = decasteljau(view(surf.controlpolygon, :, 4), Val(3), u)
    return decasteljau(SVector{4,SVector{N,S}}(q1, q2, q3, q4), Val(3), v)
end

function point(surf::BezierSurface{P,S,N,4}, u::T, v::T) where {T<:Real,P,S,N}
    q1 = decasteljau(view(surf.controlpolygon, :, 1), Val(4), u)
    q2 = decasteljau(view(surf.controlpolygon, :, 2), Val(4), u)
    q3 = decasteljau(view(surf.controlpolygon, :, 3), Val(4), u)
    q4 = decasteljau(view(surf.controlpolygon, :, 4), Val(4), u)
    q5 = decasteljau(view(surf.controlpolygon, :, 5), Val(4), u)
    return decasteljau(SVector{5,SVector{N,S}}(q1, q2, q3, q4, q5), Val(4), v)
end

function point(surf::BezierSurface{P,S,N,M}, u::T, v::T) where {T<:Real,P,S,N,M}
    qi = zeros(MVector{M + 1,SVector{N,S}})
    @inbounds for i in 1:(M + 1)
        qi[i] = decasteljau(view(surf.controlpolygon, :, i), Val(M), u)
    end
    return decasteljau(qi, Val(M), v)
end

# DERIVATIVE
@inline function derivative(controlpoints::SVector{4,SVector{N,S}}, u::T) where {T<:Real,N,S}
    d1 = controlpoints[2] .- controlpoints[1]
    d2 = controlpoints[3] .- controlpoints[2]
    d3 = controlpoints[4] .- controlpoints[3]
    return 3 * decasteljau(SVector{3,SVector{N,S}}(d1, d2, d3), Val(2), u)
end

@inline function derivative(controlpoints::SVector{5,SVector{N,S}}, u::T) where {T<:Real,N,S}
    d1 = controlpoints[2] .- controlpoints[1]
    d2 = controlpoints[3] .- controlpoints[2]
    d3 = controlpoints[4] .- controlpoints[3]
    d4 = controlpoints[5] .- controlpoints[4]
    return 4 * decasteljau(SVector{4,SVector{N,S}}(d1, d2, d3, d4), Val(3), u)
end

@inline function derivative(controlpoints::MVector{Q,SVector{N,S}}, u::T) where {T<:Real,N,S,Q}
    derivpoints = zeros(MVector{Q - 1,SVector{N,T}})
    @inbounds for i in 1:(Q - 1)
        derivpoints[i] = controlpoints[i + 1] .- controlpoints[i]
    end
    return (Q - 1) * decasteljau(derivpoints, Val(Q - 2), u)
end

# PARTIALS
function partials(surf::BezierSurface{Euclidean,S,N,3}, u::T, v::T) where {T<:Real,S,N}
    q1 = decasteljau(view(surf.controlpolygon, :, 1), Val(3), u)
    q2 = decasteljau(view(surf.controlpolygon, :, 2), Val(3), u)
    q3 = decasteljau(view(surf.controlpolygon, :, 3), Val(3), u)
    q4 = decasteljau(view(surf.controlpolygon, :, 4), Val(3), u)
    partialv = derivative(SVector{4,SVector{N,S}}(q1, q2, q3, q4), v)
    q1 = decasteljau(view(surf.controlpolygon, 1, :), Val(3), v)
    q2 = decasteljau(view(surf.controlpolygon, 2, :), Val(3), v)
    q3 = decasteljau(view(surf.controlpolygon, 3, :), Val(3), v)
    q4 = decasteljau(view(surf.controlpolygon, 4, :), Val(3), v)
    partialu = derivative(SVector{4,SVector{N,S}}(q1, q2, q3, q4), u)
    return (partialu, partialv)
end

function partials(surf::BezierSurface{Euclidean,S,N,4}, u::T, v::T) where {T<:Real,S,N}
    q1 = decasteljau(view(surf.controlpolygon, :, 1), Val(4), u)
    q2 = decasteljau(view(surf.controlpolygon, :, 2), Val(4), u)
    q3 = decasteljau(view(surf.controlpolygon, :, 3), Val(4), u)
    q4 = decasteljau(view(surf.controlpolygon, :, 4), Val(4), u)
    q5 = decasteljau(view(surf.controlpolygon, :, 5), Val(4), u)
    partialv = derivative(SVector{5,SVector{N,S}}(q1, q2, q3, q4, q5), v)
    q1 = decasteljau(view(surf.controlpolygon, 1, :), Val(4), v)
    q2 = decasteljau(view(surf.controlpolygon, 2, :), Val(4), v)
    q3 = decasteljau(view(surf.controlpolygon, 3, :), Val(4), v)
    q4 = decasteljau(view(surf.controlpolygon, 4, :), Val(4), v)
    q5 = decasteljau(view(surf.controlpolygon, 5, :), Val(4), v)
    partialu = derivative(SVector{5,SVector{N,S}}(q1, q2, q3, q4, q5), u)
    return (partialu, partialv)
end

function partials(surf::BezierSurface{Euclidean,S,N,M}, u::T, v::T) where {T<:Real,S,N,M}
    qi = zeros(MVector{M + 1,SVector{N,T}})
    # do this to force newpoints to have the type of u rather than controlpoints. When u is a dual number this makes the array the correct type.
    @inbounds for i in 1:(M + 1)
        qi[i] = decasteljau(view(surf.controlpolygon, :, i), Val(M), u)
    end
    partialv = derivative(qi, v)
    @inbounds for i in 1:(M + 1)
        qi[i] = decasteljau(view(surf.controlpolygon, i, :), Val(M), v)
    end
    partialu = derivative(qi, u)
    return (partialu, partialv)
end

# "computes partials and point simultaneously for Euclidean surface"
# function pointandpartials(surf::BezierSurface{Euclidean,S,N,M}, u::T, v::T) where {T<:Real,S,N,M}
#     qi = zeros(MVector{M + 1,SVector{N,T}})
#     # do this to force newpoints to have the type of u rather than controlpoints. When u is a dual number this makes the array the correct type.
#     @inbounds for i in 1:(M + 1)
#         qi[i] = decasteljau(view(surf.controlpolygon, :, i), M, u)
#     end
#     pt = decasteljau(qi, M, v)
#     partialv = derivative(qi, v)
#     @inbounds for i in 1:(M + 1)
#         qi[i] = decasteljau(view(surf.controlpolygon, i, :), M, v)
#     end
#     partialu = derivative(qi, u)
#     return (pt, partialu, partialv)
# end

function partials(surf::BezierSurface{Rational,S,N,M}, u::T, v::T) where {T<:Real,S,N,M}
    # computes Euclidean space partials for Rational surface. Not consistent with point, which has euclideanpoint function, but not sure why user would ever want Rational partials.
    pt = point(surf, u, v)
    du, dv = partials(surf, u, v)
    rdu = rationalcorrection(pt, du)
    rdv = rationalcorrection(pt, dv)
    return (rdu, rdv)
end

function normal(surf::BezierSurface, u::T, v::T) where {T<:Real}
    du, dv = partials(surf, u, v)
    return normalize(cross(du, dv))
end

function rationalcorrection(pt::SVector{N,S}, dpt::SVector{N,S}) where {N,S}
    # have form A/B where A = (x,y,z) and B = w component of vector. Want partial(A/B,u). Get (B*partial(A,u) - A*partial(B,u))/B*B
    A = pt[1:(N - 1)]
    B = pt[N]
    dB = dpt[N]
    dA = dpt[1:(N - 1)]
    return (B .* dA - A .* dB) / B^2
end

# """attempt to generate inline code so Zygote and Forwarddiff could differentatiate properly. Didn't work."""
# function inlinedecasteljau(curveorder, N)
#     a = "(u,$(reduce(*, ("newpoints$(k)_$j," for k in 1:curveorder,j in 1:N)))) -> let
#     "
#     for i in 1:curveorder
#         for k in 1:(curveorder + 1 - i)
#             for j in 1:N
#                 a *= "newpoints$(k)_$j = newpoints$(k)_$j * (1-u) + newpoints$(k + 1)_$j * u
#                  "
#             end
#         end
#     end
#     a *= "return ($(reduce(*, ("newpoints1_$j," for j in 1:N))))
#     end"

#     return eval(Meta.parse(a))
# end

# "Creates a differentiable function to evaluate a surface point. Zygote has trouble with mutable arrays so this creates a function
# which unrolls all loops and converts array indexing to variables. Stopped working on this because it became so diffult to make it work."
# function pointinline(surf::BezierSurface{P,S,N,M}, u::T, v::T) where {T<:Real,P,S,N,M}
#     controlpoints = surf.controlpolygon
#     qi = zeros(MVector{N * (M + 1),MVector{N,S},1})
#     pointinl = inlinedecasteljau(M, N)
#     @inbounds for i in 1:(M + 1)
#         start = (i - 1) * N + 1
#         @inbounds for j in 1:N
#             temp = pointinl(u, controlpoints[:, i]...)
#             qi[start + j - 1] = temp[j]
#         end
#     end
#     return pointinl(v, qi...)
# end

function point(curve::BezierCurve{P,S,N,M}, u::T) where {T<:Real,P,S,N,M}
    # Will return homogeneous curve point if curve is homogeneous and euclidean point if curve is Euclidean.
    return decasteljau(curve.controlpolygon, M, u)
end

euclideanpoint(curve::BezierCurve{Euclidean,S,N,M}, u::T) where {T<:Real,S,N,M} = point(curve, u)

function euclideanpoint(curve::BezierCurve{Rational,S,N,M}, u::T) where {T<:Real,S,N,M}
    temp = point(curve, u)
    return toeuclidean(temp)
end

#########################################################################################################################

function inside(surf::AcceleratedParametricSurface{S,3,BezierSurface{P,S,3,M}}, p::SVector{3,S}) where {P,S<:Real,M}
    # approximate the normal to the bezier surface as a whole by taking the average of the normals of the triangles
    # defined by the extreme points of the surface
    A = point(surf, 0.0, 0.0)
    B = point(surf, 0.0, 1.0)
    C = point(surf, 1.0, 0.0)
    D = point(surf, 1.0, 1.0)
    n1 = cross(C - A, B - A)
    n2 = cross(B - D, C - D)
    approx_normal = (n1 + n2) / (norm(n1) + norm(n2))
    r = Ray(p, approx_normal)
    @inbounds for i in 1:length(surf.triangles)
        intv = surfaceintersection(surf.triangles[i], r)
        if !(intv isa EmptyInterval)
            return true
        end
    end
    return false
end

function BoundingBox(surf::BezierSurface{P,S,N,M}) where {P,S,N,M}
    big = fill(typemin(S), MVector{N,S})
    small = fill(typemax(S), MVector{N,S})

    for pt in surf.controlpolygon
        small .= min.(small, pt)
        big .= max.(big, pt)
    end

    return BoundingBox(small[1], big[1], small[2], big[2], small[3], big[3])
end


text:
#
# SamplerChoiceModel
#

# define abstract sampler type, the type tree of which it is a root, and concrete implementations of these types
include(joinpath("samplers", "sampler.jl"))

# The sampler choice model is a set of samplers, one assigned to each choice point (indexed by choice point id)
type SamplerChoiceModel <: ChoiceModel
	samplers::Dict{UInt, Sampler}
	maxresamplings::Int
	function SamplerChoiceModel(cpinfo::Dict, choicepointmapping::Function=defaultchoicepointmapping, maxresamplings::Int=9)
		maxresamplings >= 0 || error("maxresamplings cannot be negative")
		samplers = Dict{UInt, Sampler}()
		for (cpid, info) in cpinfo # gets info from sub-generators also
			samplers[cpid] = choicepointmapping(info)
		end
		new(samplers, maxresamplings)
	end
end

function SamplerChoiceModel(g::Generator; choicepointmapping::Function=defaultchoicepointmapping, maxresamplings::Int=9)
	SamplerChoiceModel(choicepointinfo(g), choicepointmapping, maxresamplings)
end

function setsamplerchoicemodel!(g::Generator; choicepointmapping::Function=defaultchoicepointmapping, maxresamplings::Int=9)
	Base.depwarn("Use setchoicemodel!(g, SamplerChoiceModel(g, ...)) in place of setsamplerchoicemodel!(g, ...)", :setsamplerchoicemodel!)
	setchoicemodel!(g, SamplerChoiceModel(g, choicepointmapping=choicepointmapping, maxresamplings=maxresamplings))
end

function defaultchoicepointmapping(info::Dict)
	cptype = info[:type]
	if cptype == :rule
		sampler = CategoricalSampler(info[:max])
		# justification: expected be a small range, is guaranteed to be closed, and no meaning is attached to order
		# varying support: not possible
	elseif cptype == :sequence
		sampler = AlignMinimumSupportSampler(GeometricSampler())
		# justification: often builds recursive structures, so smaller number of repetitions should be more likely
		# varying support: lower handled by AlignMinimumSupportSampler; upper by rejection sampling
	elseif cptype == :value
		datatype = info[:datatype]
		if datatype <: Bool
			sampler = BernoulliSampler()
			# justification: two mutually exclusive outcomes, so natural choice
			# varying support: not possible
		elseif datatype <: Integer # order of if clauses matters here since Bool <: Integer
			sampler = AdjustParametersToSupportSampler(DiscreteUniformSampler())
			# justification: simplicity
			# varying support: range of distribution is dynamically modified using AdjustParametersToSupport
		else # floating point, but may also be a rational type
			sampler = AdjustParametersToSupportSampler(UniformSampler())
			# justification: simplicity
			# varying support: range of distribution is dynamically modified using AdjustParametersToSupport
		end
	else
		error("unrecognised choice point type when creating sampler choice model")
	end
	sampler
end

#
# DataGenerators interface
#

# retrieve the sampler for this choice point and then ask it to return a godel number
# resample if necessary - usually because upperbound cannot be enforced directly by sampler,
# but eventually fall back to uniform distribution
function godelnumber(cm::SamplerChoiceModel, cc::ChoiceContext)
	lowerbound = isfinite(cc.lowerbound) ? cc.lowerbound : sign(cc.lowerbound) * realmax(Float64)
	upperbound = isfinite(cc.upperbound) ? cc.upperbound : sign(cc.upperbound) * realmax(Float64)
	sampler = cm.samplers[cc.cpid]
	fallback = true
	x, trace = nothing, nothing
	samplecount = 0
	while (samplecount <= cm.maxresamplings) && (fallback)
		x, trace = sample(sampler, (lowerbound, upperbound), cc::ChoiceContext)
		if lowerbound <= x <= upperbound
			fallback = false
		end
		samplecount += 1
	end
	if fallback
		warn("falling back to a uniform distribution after too many resamplings in sampler choice model")
		if cc.datatype <: Integer # includes Bool
			fallbacksampler = DiscreteUniformSampler([float(lowerbound), float(upperbound)])
		elseif cc.datatype <: AbstractFloat
			fallbacksampler = UniformSampler([lowerbound, upperbound])
		else
			@assert false
		end
		x, fallbacktrace = sample(fallbacksampler, (lowerbound, upperbound), cc::ChoiceContext)
		# we use the last trace form the 'proper' sampler but place in it this fallback sampled value
		amendtrace(sampler, trace, x)
	end
	cptrace = Dict(:sam=>trace, :fbk=>fallback)
	# some overhead to recording both the following, so could control using model-level parameters
	cptrace[:acs] = getimmediateancestorcpseqnumber(cc.derivationstate)
	cptrace[:rcd] = getcurrentrecursiondepth(cc.derivationstate)
	return x, cptrace
end

# Note: paramranges, setparams and getparams the following assumption: to guarantee the ordering both within a session,
# and across all uses of the same generator, the samplers are accessed in sort order of the key.
# (Currently the key is a UInt that is assigned in a deterministic order in the generator macro.)

# get parameter ranges for all samplers
function paramranges(cm::SamplerChoiceModel)
	ranges = Tuple{Float64,Float64}[]
	for k in sort(collect(keys(cm.samplers))) # guarantees ordering
		sampler = cm.samplers[k]
		ranges = [ranges; paramranges(sampler)]
	end
	ranges
end

# set parameters of all samplers
function setparams!(cm::SamplerChoiceModel, params)
	params = convert(Vector{Float64},params)
	if length(params) != numparams(cm)
		error("expected $(numparams(cm)) model parameter(s), but got $(length(params))")
	end
	idx = 1
	for k in sort(collect(keys(cm.samplers))) # guarantees ordering
		sampler = cm.samplers[k]
		nparams = numparams(sampler)
		# TODO replace this with a check on getnummodelparams
		@assert (idx+nparams-1)<=length(params)
		setparams!(sampler, params[idx:(idx+nparams-1)])
		# note: it is important here that sampler works on a copy of the parameters (which is ensured by [range] syntax)
		# as the sampler may adjust its parameters
		idx += nparams
	end
end

# get parameters of all samplers
function getparams(cm::SamplerChoiceModel)
	params = (Float64)[]
	for k in sort(collect(keys(cm.samplers))) # guarantees ordering
		sampler = cm.samplers[k]
		params = [params; getparams(sampler)]
	end
	params
end

# extract dict of trace info indexed by cpid from a vector of cm traces
function extractsamplertracesbycp(cm::SamplerChoiceModel, cmtraces)
	samplertracesbycp = Dict{UInt,Vector{Dict}}()
	for cmtrace in cmtraces
		for (cpid, trace) in cmtrace
			if !haskey(samplertracesbycp, cpid)
				samplertracesbycp[cpid] = Dict[]
			end
			push!(samplertracesbycp[cpid], trace[:sam])
		end
	end
	samplertracesbycp
end

# estimate the parameters of each sampler based on vector of traces from generation states
function estimateparams!(cm::SamplerChoiceModel, cmtraces)
	samplertracesbycp = extractsamplertracesbycp(cm, cmtraces)
	for cpid in keys(samplertracesbycp)
		estimateparams!(cm.samplers[cpid], samplertracesbycp[cpid])
	end
end

# at each cp in every trace, record the most recent gn emitted by every choice point (i.e. a "history")
# store this by cpid
# for ease of use in re-estimating sampler params, we simultaneously extract sampler traces for each history
function extractgnhistoriesbycp(cm::SamplerChoiceModel, cmtraces, recencywindow::Int, distinguishparentsbyrecursiondepth::Bool, restricttoancestors::Bool)

	# current history values at recent choice point, distinguised by recency (up to recencywindow)
	# note: to avoid creating unused cplabels (cp instances distinguished by recency, and optionally recursion depth),
	# that the model estimation code must neverthless process, we build up the cplabels as we go, and
	# then later pad to histories to the full length once all traces have been processed
	# note also: the last entry in each history is the current godel number

	# tuples that label each index in the history vector: tuple is (cpid, recursiondepth, recency, restricttoancestors)
	cplabels = Vector{Tuple{UInt, Any, Int, Bool}}() # Any is type of recursion depth since this is set to nothing if it is not required
	# provides quick lookup of index given a cplabel
	cplabelidxlookup = Dict{Tuple{UInt, Any, Int, Bool}, Int}()

	# process each trace
	gnhistoriesbycp = Dict{UInt, Vector{Vector{Any}}}()
	samplertracesbycp = Dict{UInt, Vector{Dict}}()

	# ... we do this in different ways depending on whether we need to restrict the history to ancestors or not;
	# in the former case, we can't keep a "running" history, but instead must re-evaluate the history each time
	# in order to filter according to ancestry

	if restricttoancestors

		for cmtrace in cmtraces

			for i in 1:length(cmtrace) # consider each choice point in trace in turn

				(currentcpid, currenttrace) = cmtrace[i]
				history = convert(Vector{Any}, fill(nothing, length(cplabels)))

				# we can traverse back through ancestor choice points using the ancestor sequence numbers stored in the trace
				ancestorcpseqnumber = currenttrace[:acs]
				while ancestorcpseqnumber != nothing

					(earliercpid, earliertrace) = cmtrace[ancestorcpseqnumber]

					earlierrecursiondepth = distinguishparentsbyrecursiondepth ? earliertrace[:rcd] : nothing

					# now, look for recency entry that has not yet been filled
					recency, historyrecorded = 1, false
					while (recency <= recencywindow) && (!historyrecorded)
						hidx = get(cplabelidxlookup, (earliercpid, earlierrecursiondepth, recency, restricttoancestors), 0)
						if hidx > 0
							if history[hidx] == nothing
								history[hidx] = earliertrace[:gdl]
								historyrecorded = true
							end
						else
							push!(history, earliertrace[:gdl])
							historyrecorded = true
							push!(cplabels, (earliercpid, earlierrecursiondepth, recency, restricttoancestors))
							cplabelidxlookup[(earliercpid, earlierrecursiondepth, recency, restricttoancestors)] = length(cplabels)
						end
						recency += 1
					end

					ancestorcpseqnumber = earliertrace[:acs]

				end # ancestor choice points

				if !haskey(gnhistoriesbycp, currentcpid)
					gnhistoriesbycp[currentcpid] = Vector{Vector{Any}}()
				end
				# add current Gödel number as *last* element
				push!(gnhistoriesbycp[currentcpid], [history; currenttrace[:gdl]]) # this concatenation also causes implicit copy

				# also store trace for this sampler (so in the same order as the history)
				if !haskey(samplertracesbycp, currentcpid)
					samplertracesbycp[currentcpid] = Vector{Dict}()
				end
				push!(samplertracesbycp[currentcpid], currenttrace[:sam])

			end # each choice point

		end # each trace

	else

		for cmtrace in cmtraces

			# in this case (without filtering by ancestry) we can build up the history
			# as we process each choice point in chronoligical order
			runninghistory = convert(Vector{Any}, fill(nothing, length(cplabels)))

			for cptrace in cmtrace

				(cpid, trace) = cptrace

				# store (copy of) running history against this cpid
				if !haskey(gnhistoriesbycp, cpid)
					gnhistoriesbycp[cpid] = Vector{Vector{Any}}()
				end
				gn = trace[:gdl]
				# add current Gödel number as *last* element
				push!(gnhistoriesbycp[cpid], [runninghistory; gn]) # this concatenation also causes implicit copy

				recursiondepth = distinguishparentsbyrecursiondepth ? trace[:rcd] : nothing

				# now update history with this choice point's Gödel number
				# first shift down values at other recencies
				hidx = get(cplabelidxlookup, (cpid, recursiondepth, recencywindow, restricttoancestors), 0)
				for recency in (recencywindow-1):-1:1
					newerhidx = get(cplabelidxlookup, (cpid, recursiondepth, recency, restricttoancestors), 0)
					@assert (newerhidx != 0) || (hidx == 0) # if no label for current recency, then can't be one for an older recency
					if newerhidx > 0
						if hidx > 0
							runninghistory[hidx] = runninghistory[newerhidx]
						else
							push!(runninghistory, runninghistory[newerhidx])
							push!(cplabels, (cpid, recursiondepth, recency+1, restricttoancestors))
							cplabelidxlookup[(cpid, recursiondepth, recency+1, restricttoancestors)] = length(cplabels)
						end
					end
					hidx = newerhidx
				end
				# and then store current Gödel number at recency 1
				if hidx > 0
					runninghistory[hidx] = trace[:gdl]
				else
					push!(runninghistory, trace[:gdl])
					push!(cplabels, (cpid, recursiondepth, 1, restricttoancestors))
					cplabelidxlookup[(cpid, recursiondepth, 1, restricttoancestors)] = length(cplabels)
				end

				# also store trace for this sampler (so in the same order as the history)
				if !haskey(samplertracesbycp, cpid)
					samplertracesbycp[cpid] = Vector{Dict}()
				end
				push!(samplertracesbycp[cpid], trace[:sam])

			end # each choice point

		end # each trace

	end

	# now pad histories so they are all the same length (which should be cplabels + 1 since first entry is actually gn of cp at
	# time the history was taken
	targetlength = length(cplabels) + 1
	for gnhistories in values(gnhistoriesbycp)
		for gnhistory in gnhistories
			padlength = targetlength - length(gnhistory)
			@assert padlength >= 0
			if padlength > 0
				# current gn was recorded as the last entry at the time of location, so remove and then add again after padding
				currentgn = pop!(gnhistory)
				append!(gnhistory, fill(nothing, padlength))
				push!(gnhistory, currentgn)
			end
		end
	end

	cplabels, gnhistoriesbycp, samplertracesbycp

end

function estimateconditionalmodel!(cm::SamplerChoiceModel, cmtraces; recencywindow::Int = 1, distinguishparentsbyrecursiondepth::Bool = false, restricttoancestors::Bool=false)

	recencywindow >= 1 || error("recency window cannot be less than 1")

	# extract the history (the most recent value of all choice points) at each choice point invocation in all traces
	cplabels, gnhistoriesbycp, samplertracesbycp = extractgnhistoriesbycp(cm, cmtraces, recencywindow, distinguishparentsbyrecursiondepth, restricttoancestors)

	# process each cp
	# TODO PARALLELISE!!!
	for cpid in keys(samplertracesbycp)
		# note: if samplertracsbycp has an entry for cpid, then so will gnhistoriesbycp
		if supportsconditionalmodelestimation(cm.samplers[cpid])
			estimateconditionalmodel!(cm.samplers[cpid], cplabels, gnhistoriesbycp[cpid], samplertracesbycp[cpid])
		else
			estimateparams!(cm.samplers[cpid], samplertracesbycp[cpid])
		end
	end

end

# pretty print the model

function show(io::IO, cm::SamplerChoiceModel)

	print(io, "sampler choice model")

	# now changed this method to override show, can no longer pass generator (TODO: could include ref to generator in cm)
	# # first build slightly more informative names for choice points
	# cpnames = Dict{UInt, AbstractString}()
	# cpinfos = choicepointinfo(g)
	# orderedcpids = sort(collect(keys(cpinfos)))
	# for cpid in orderedcpids
	# 	cpinfo = cpinfos[cpid]
	# 	cpname = "($(cpid))"
	# 	cptype = cpinfo[:type]
	# 	if cptype == :rule
	# 		cpname = "Rule $(cpinfo[:rulename]) " * cpname
	# 	elseif cptype == :sequence
	# 		cpname = "Sequence " * cpname
	# 	elseif cptype == :value
	# 		cpname = "Value $(cpinfo[:datatype]) " * cpname
	# 	end
	# 	cpnames[cpid] = cpname
	# end

	# now print samplers in id order (which should (roughly) be lexical order of choice points?)
	# for cpid in sort(collect(keys(cm.samplers)))
	# 	print(io, "$(hex(cpid)): ")
	# 	show(io, cm.samplers[cpid])
	# end

end