text:
using IterTools, StatsBase
export imcmc_multi_split_involution!, imcmc_multi_merge_involution!
export imcmc_multi_split_prop_sample!, imcmc_multi_merge_prop_sample!
export imcmc_special_multi_split_prop_sample!, imcmc_special_multi_merge_prop_sample!
export log_binomial_ratio_split_merge

function imcmc_multi_split_involution!(
    S_curr::InteractionSequence{T},
    S_prop::InteractionSequence{T},
    ind_split::Vector{Int},
    ind_loc::Vector{Int},
    k_vec::Vector{Int}, # Split locations (in each path),
    curr_pointers::InteractionSequence{T},
    prop_pointers::InteractionSequence{T}
    ) where {T}

    store = Vector{T}[]
    log_ratio_inc = 0.0 # Increment to log_ratio
    for (i,k) in Iterators.reverse(zip(ind_split, k_vec))
        tmp1 = popat!(S_prop, i)
        n = length(tmp1)
        @assert n>1 "Cannot split path of length one."
        log_ratio_inc += log(n-1)
        tmp2 = popfirst!(prop_pointers)
        copy!(tmp2, tmp1[(k+1):end])
        resize!(tmp1, k)
        pushfirst!(store, tmp2)
        pushfirst!(store, tmp1)
    end

    for i in ind_loc
        tmp = popfirst!(store)
        insert!(S_prop, i, tmp)
    end

    return log_ratio_inc

end

function imcmc_multi_split_prop_sample!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    mcmc::T
    ) where {T<:SimPosteriorSampler}

    ν = mcmc.ν_td
    N = length(S_curr)
    ε = rand(1:min(N, ν))
    ind_split = view(mcmc.ind_td_split, 1:ε)
    ind_loc = view(mcmc.ind_td_merge, 1:2ε)
    prop_pointers = mcmc.curr_pointers

    store = mcmc.store
    empty!(store)

    log_ratio = 0.0

    StatsBase.seqsample_a!(1:N, ind_split)
    StatsBase.seqsample_a!(1:(N+ε), ind_loc)

    for i in Iterators.reverse(ind_split)
        @assert length(S_prop[i])>1 "Cannot split path of length one."
        tmp1 = popat!(S_prop, i)
        n = length(tmp1)
        k = rand(1:(n-1)) # Split location
        log_ratio += log(n-1)
        tmp2 = popfirst!(prop_pointers)
        copy!(tmp2, tmp1[(k+1):end])
        resize!(tmp1, k)
        pushfirst!(store, tmp2)
        pushfirst!(store, tmp1)
    end

    for i in ind_loc
        tmp = popfirst!(store)
        insert!(S_prop, i, tmp)
    end

    log_ratio += log(min(N, ν)) - log(min(floor(length(S_prop/2)), ν))

    return log_ratio

end

function imcmc_multi_merge_involution!(
    S_curr::InteractionSequence{T},
    S_prop::InteractionSequence{T},
    ind_merge::Vector{Int},
    ind_loc::Vector{Int},
    curr_pointers::InteractionSequence{T},
    prop_pointers::InteractionSequence{T}
    ) where {T}

    store = Vector{T}[]
    log_ratio_inc = 0.0
    i = length(ind_merge)
    while i > 0
        # Merge i-1 and i
        tmp2 = popat!(S_prop, ind_merge[i])
        tmp1 = popat!(S_prop, ind_merge[i-1])
        append!(tmp1, tmp2)
        pushfirst!(prop_pointers, tmp2)
        pushfirst!(store, tmp1)
        log_ratio_inc -= log(length(tmp1)-1)
        i -= 2
    end

    for i in ind_loc
        tmp = popfirst!(store)
        insert!(S_prop, i, tmp)
    end

    return log_ratio_inc
end

function imcmc_multi_merge_prop_sample!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    mcmc::T
    ) where {T<:SimPosteriorSampler}

    ν = mcmc.ν_td
    N = length(S_curr)
    ε = rand(1:min(floor(Int, N/2), ν))

    ind_merge = view(mcmc.ind_td_merge, 1:2ε)
    ind_loc = view(mcmc.ind_td_split, 1:ε)

    store = mcmc.store
    empty!(store)
    prop_pointers = mcmc.prop_pointers

    log_ratio = 0.0

    StatsBase.seqsample_a!(1:N, ind_merge)
    StatsBase.seqsample_a!(1:(N-ε), ind_loc)

    i = length(ind_merge)
    while i > 0
        # Merge i-1 and i
        tmp2 = popat!(S_prop, ind_merge[i])
        tmp1 = popat!(S_prop, ind_merge[i-1])
        append!(tmp1, tmp2)
        pushfirst!(prop_pointers, tmp2)
        pushfirst!(store, tmp1)
        log_ratio -= log(length(tmp1)-1)
        i -= 2
    end

    for i in ind_loc
        tmp = popfirst!(store)
        insert!(S_prop, i, tmp)
    end

    log_ratio += log(min(floor(Int, N/2), ν)) - log(min(length(S_prop), ν))

    return log_ratio
end

function log_binomial_ratio_split_merge(N, τ, ε)
    # \binom{N-τ}{ε}/\binom{N}{ε}
    out = 0.0
    if τ < ε
        for j in 0:(τ-1)
            out += log(N-ε-j) - log(N-j)
        end
    else
        for j in 0:(ε-1)
            out += log(N-τ-j) - log(N-j)
        end
    end
    return out
end

# Special case (when we have paths of length one present)
function imcmc_special_multi_split_prop_sample!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    mcmc::T
    ) where {T<:SimPosteriorSampler}

    ν = mcmc.ν_td
    N = length(S_curr)
    ind_ok_len = findall(x->length(x)>1, S_curr) # Good paths (length > 1)
    τ = N-length(ind_ok_len)                     # Number of bad paths
    ε = rand(1:min(N-τ, ν))
    ind_split = view(mcmc.ind_td_split, 1:ε)
    ind_loc = view(mcmc.ind_td_merge, 1:2ε)
    prop_pointers = mcmc.curr_pointers

    store = mcmc.store
    empty!(store)

    log_ratio = 0.0

    StatsBase.seqsample_a!(ind_ok_len, ind_split)  # Sample subsequence of good paths
    StatsBase.seqsample_a!(1:(N+ε), ind_loc)

    for i in Iterators.reverse(ind_split)
        @assert length(S_prop[i])>1 "Cannot split path of length one."
        tmp1 = popat!(S_prop, i)
        n = length(tmp1)
        k = rand(1:(n-1)) # Split location
        log_ratio += log(n-1)
        tmp2 = popfirst!(prop_pointers)
        copy!(tmp2, tmp1[(k+1):end])
        resize!(tmp1, k)
        pushfirst!(store, tmp2)
        pushfirst!(store, tmp1)
    end

    for i in ind_loc
        tmp = popfirst!(store)
        insert!(S_prop, i, tmp)
    end

    log_ratio += log(min(N-τ, ν)) - log(min(floor(length(S_prop/2)), ν))

    log_ratio += log_binomial_ratio_split_merge(N, τ, ε)

    return log_ratio

end

function imcmc_special_multi_merge_prop_sample!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    mcmc::T
    ) where {T<:SimPosteriorSampler}

    ν = mcmc.ν_td
    N = length(S_curr)
    ε = rand(1:min(floor(Int, N/2), ν))

    ind_merge = view(mcmc.ind_td_merge, 1:2ε)
    ind_loc = view(mcmc.ind_td_split, 1:ε)

    store = mcmc.store
    empty!(store)
    prop_pointers = mcmc.prop_pointers

    log_ratio = 0.0

    StatsBase.seqsample_a!(1:N, ind_merge)
    StatsBase.seqsample_a!(1:(N-ε), ind_loc)

    i = length(ind_merge)
    while i > 0
        # Merge i-1 and i
        tmp2 = popat!(S_prop, ind_merge[i])
        tmp1 = popat!(S_prop, ind_merge[i-1])
        append!(tmp1, tmp2)
        pushfirst!(prop_pointers, tmp2)
        pushfirst!(store, tmp1)
        log_ratio -= log(length(tmp1)-1)
        i -= 2
    end

    for i in ind_loc
        tmp = popfirst!(store)
        insert!(S_prop, i, tmp)
    end

    M = length(S_prop)
    τ = sum(length(x)==1 for x in S_prop) # Find number of length one paths
    log_ratio -= log_binomial_ratio_split_merge(M, τ, ε)

    log_ratio += log(min(floor(Int, N/2), ν)) - log(min(M-τ, ν))

    return log_ratio
end

# Adding necessary functions for mcmc samplers to work
# We first need the accept reject function for this move
function double_iex_trans_dim_accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SimPosterior,
    γ_curr::Float64,
    mcmc::SimIexSplitMerge,
    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_lb, K_in_ub = (K_inner.l, K_inner.u)
    data = posterior.data
    dist = posterior.dist
    V = posterior.V
    γ_prior = posterior.S_prior.γ
    mode_prior = posterior.S_prior.mode

    curr_pointers = mcmc.curr_pointers
    prop_pointers = mcmc.prop_pointers
    aux_mcmc = mcmc.aux_mcmc

    log_ratio = 0.0

    do_split = false

    τ = sum(length(x)==1 for x in S_curr)
    N = length(S_curr)
    # If current state has one path of length one, nothing can be done
    if (N==τ==1)
        return 0, suff_stat_curr
    # If one path of length > 1 we split with probability one
    elseif N==1
        do_split = true
        log_ratio += imcmc_multi_split_prop_sample!(S_curr, S_prop, mcmc)
        log_ratio += all(length(x)==1 for x in S_prop) ? 0.0 : log(0.5)
    # If no length one paths, use std proposal fn
    elseif τ==0
        do_split = rand(Bool)
        if do_split
            log_ratio += imcmc_multi_split_prop_sample!(S_curr, S_prop, mcmc)
        else
            log_ratio += imcmc_multi_merge_prop_sample!(S_curr, S_prop, mcmc)
        end
    # If all paths are length one, do merge with probability one
    elseif τ==N
        do_split = false
        log_ratio += imcmc_special_multi_merge_prop_sample!(S_curr, S_prop, mcmc)
        log_ratio += log(0.5) # Accounts for uneven choice between this and split move
    # Else to use special proposal function
    else
        do_split = rand(Bool)
        if do_split
            log_ratio += imcmc_special_multi_split_prop_sample!(S_curr, S_prop, mcmc)
        else
            log_ratio += imcmc_special_multi_merge_prop_sample!(S_curr, S_prop, mcmc)
        end
    end

    # Account for dimension bounds (make acc_prob 0 if outside bounds)
    if !(K_out_lb ≤ length(S_prop) ≤ K_out_ub)
        log_ratio += -Inf
    end
    if any(!(K_in_lb ≤ length(x) ≤ K_in_ub) for x in S_prop)
        log_ratio += -Inf
    end

    # Now do accept-reject step (**THIS IS WHERE WE DIFFER FROM MODEL SAMPLER***)
    aux_model = SIM(
        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_multinom_term = log_multinomial_ratio(S_curr, S_prop)

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

    ε = abs(length(S_curr)-length(S_prop))
    if log(rand()) < log_α
        # We accept!
        if do_split
            # We've increased dimension
            # Bring in pointers from curr_pointers
            for i in 1:ε
                tmp = pop!(curr_pointers)
                push!(S_curr,tmp)
            end
        else
            # We've decreased dimension
            # Send pointers back to curr_pointers
            for i in 1:ε
                tmp = pop!(S_curr)
                push!(curr_pointers, tmp)
            end
        end
        # Copy across from S_prop
        for (x,y) in zip(S_curr,S_prop)
            copy!(x, y)
        end
        return 1, suff_stat_prop
    else
        # We reject!
        if do_split
            # We've increased dimension
            # Send back pointers to prop_pointers
            for i in 1:ε
                tmp = pop!(S_prop)
                push!(prop_pointers, tmp)
            end
        else
            # We've decreased dimension
            # Bring in pointers from prop_pointers
            for i in 1:ε
                tmp = pop!(prop_pointers)
                push!(S_prop, tmp)
            end
        end
        # Copy across from S_curr
        for (x,y) in zip(S_prop,S_curr)
            copy!(x, y)
        end
        return 0, suff_stat_curr
    end

end

function draw_sample_mode!(
    sample_out::Union{InteractionSequenceSample{Int}, SubArray},
    mcmc::SimIexSplitMerge,
    posterior::SimPosterior,
    γ_fixed::Float64;
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::InteractionSequence{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) * lag + burn_in, # 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 = SIM(
        S_curr, γ_curr,
        posterior.dist,
        posterior.V,
        posterior.K_inner,
        posterior.K_outer)
    draw_sample!(aux_data, aux_mcmc, aux_model, burn_in=10000)

    # 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.α)
    while sample_count ≤ length(sample_out)
        i += 1
        # Store value
        if (i > burn_in) & (((i - burn_in - 1) % lag)==0)
            @inbounds 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_accept_reject!(
                S_curr, S_prop,
                posterior, γ_curr,
                mcmc,
                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

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

# Mode accept-reject

function accept_reject_mode!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SimPosterior,
    γ_curr::Float64,
    mcmc::SimIexSplitMerge,
    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=false
    )

    β = 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_accept_reject!(
            S_curr, S_prop,
            posterior, γ_curr,
            mcmc,
            aux_data,
            suff_stat_curr,
            aux_init_at_prev
        )
        acc_count[3] += was_accepted
        count[3] += 1
    end
    return suff_stat_curr
end


text:
# script file to reproduce benchmarks in the analysis at
# TODO http ref
using SearchBenchmarks
using Test

path = "" # working directory including path delimiter at end
source = path*"7z.dll"

benchmark(200, 1000000, source, [naivesearch,bloom_v0], path*"(1) BIN julia impl much better than naive search",false)
benchmark(200, 1000000, source, [bloom_v0,bloom_v1,bloom_v2,bloom_v3], path*"(2) BIN julia bloom impl optimizations",false)
benchmark(200, 1100000, source, [bloom_v0,bloom_v3,bloom_v4], path*"(3) BIN do bloom test first",false)
benchmark(2000, 1000000, source, [naivesearch,bloom_v0,bloom_v3,bloom_v4], path*"(4) BIN degeneration on large patterns",false)
benchmark(2000, 1000000, source, [naivesearch,bloom_v0,bloom_v3,bloom_best], path*"(5) BIN limited bits in bloom filter",false)
benchmark(200, 1000000, source, [bloom_v0,bloom_v3,bloom_best,bmsearch], path*"(6) BIN boyer moore alternative",false)
benchmark(2000, 1000000, source, [naivesearch,bloom_v0,bloom_v3,bloom_best,bmsearch], path*"(7) BIN boyer moore on large patterns",false)

ext: jl


text:
#============================================

Profiling for Skiplist generation.

============================================#

pushfirst!(LOAD_PATH, joinpath(@__DIR__, ".."))

using Logging
using Profile
using Random
using SkipLists

Profile.init(; delay=1e-4)

const OUTPUT_PATH = joinpath("results", "profile_generation.txt")
const N = 1_000_000

# Pre-compile code
list = SkipList{Int64}()
for ii = 1:100
    insert!(list, 0)
end

Profile.clear_malloc_data()

# Perform profiling
for ii = 1:5
    @info "Iteration $ii"
    list = SkipList{Int64}()
    for jj = shuffle(1:N)
        @profile insert!(list, jj)
    end
end

open(OUTPUT_PATH, "w") do f
    for outfile in (f, stdout)
        Profile.print(IOContext(outfile, :displaysize => (24,500)); maxdepth=40)
    end
end

println("Output written to $OUTPUT_PATH")


text:
module TestLocalization

using GeoIP
using Sockets: IPv4, @ip_str
using Test

TEST_DATADIR = joinpath(dirname(@__FILE__), "data")

ENV["GEOIP_DATADIR"] = TEST_DATADIR
ENV["GEOIP_ZIPFILE"] = "GeoLite2-City-CSV.zip"

@testset "Locales setup" begin
    db = load(locales = [:en, :ru])
    geoip1 = db[ip"1.0.8.1"]
    @test geoip1["locale_code"] == "en"
    @test geoip1["continent_name"] == "Asia"

    db2 = setlocale(db, :ru)
    geoip1 = db2[ip"1.0.8.1"]
    @test geoip1["locale_code"] == "ru"
    @test geoip1["continent_name"] == "Азия"
end

@testset "Missed locales" begin
    db = load(locales = [:en, :ru], deflocale = :ru)
    geoip1 = db[ip"1.0.1.1"]
    @test !haskey(geoip1, "locale_code")

    db2 = setlocale(db, :en)
    geoip1 = db2[ip"1.0.1.1"]
    @test geoip1["locale_code"] == "en"
end

end # module


text:
using ImageFiltering: centered, KernelFactors.gaussian, kernelfactors, imfilter
using ImageDistances
using TimerOutputs

"""
    move_landmarks(locations, flow)

Move each `location::Tuple(T,T) where T <: Real` from the `locations` vector using the `flow`.
"""
function move_landmarks(locations, flow)
    itp = Interpolations.interpolate(flow, BSpline(Linear()))
    shifted_locations = Array{Float64}(undef, size(locations))
    for k in 1:size(locations, 1)
        shifted_locations[k,:] .= locations[k,:] .+ [imag(itp(locations[k,:]...)), real(itp(locations[k,:]...))]
    end
    return shifted_locations
end

"""
    pad_images(image1::Image, image2::Image)

Adds zeros to the right and bottom of `image1` and `image2` to make them the same size.
"""
function pad_images(image1::Image, image2::Image)

    (a, b) = size(image1)
    (c, d) = size(image2)

    # if there is nothing to do, end
    if (a, b) == (c, d)
        return image1, image2
    end

    if (a < c)
        image1 = [image1; zeros(eltype(image1), c - a, b)];
        a = c
    elseif (a > c)
        image2 = [image2; zeros(eltype(image2), a - c, d)];
        c = a
    end

    if (b < d)
        image1 = [image1 zeros(eltype(image1), a, d - b)];
    elseif (b > d)
        image2 = [image2 zeros(eltype(image2), c, b - d)];
    end
    return image1, image2
end

"""
    pad_images(image1::Image, image2::Image)

Adds zeros to the right and bottom of `image1` and `image2` to make them the same size. Also create a mask,
which has zeros at the edges of added zeros.
"""
function pad_images_mask(target::Image, source::Image)

    (a, b) = size(target)
    (c, d) = size(source)

    mask = trues(a > c ? a : c, b > d ? b : d)

    # if there is nothing to do, end
    if (a, b) == (c, d)
        return target, source, mask
    end

    if (a < c)
        target = [target; zeros(eltype(target), c - a, b)];
        # add horizontal mask
        mask[a, :] .= false
        mask[a+1, :] .= false
        a = c
    elseif (a > c)
        source = [source; zeros(eltype(source), a - c, d)];
        c = a
    end

    if (b < d)
        target = [target zeros(eltype(target), a, d - b)];
        # add vertical mask
        mask[:, b] .= false
        mask[:, b+1] .= false
    elseif (b > d)
        source = [source zeros(eltype(source), c, b - d)];
    end
    return target, source, mask
end

"""
    normalize_to_zero_one(arr)

Normalize `arr` so that it has a maximum of 1 and minimum of 0.
"""
function normalize_to_zero_one(arr)
    normalized = arr .- minimum(arr)
    normalized = normalized ./ maximum(normalized)
    return normalized
end

"""
    rescale!(image1::Image, image2::Image)

Rescale `image1` and `image2` intensities to span the whole `[0, 1]`.
"""
function rescale!(image1::Image, image2::Image)

    max_intensity = maximum([image1 image2])
    min_intensity = minimum([image1 image2])

    image1 = (image1 .- min_intensity)./(max_intensity - min_intensity);
    image2 = (image2 .- min_intensity)./(max_intensity - min_intensity);

    return image1, image2
end

"""
    smooth_with_gaussian!(A::Matrix{<:Number}, window_half_size_one_dim::Integer)

Clean the Matrix `A` by smoothing using a square 2D Gaussian filter of size `2 * window_half_size_one_dim + 1` in each dimension.
"""
@inline function smooth_with_gaussian!(A::Matrix{<:Number}, window_half_size_one_dim::Integer)# where T<:Integer
    window_half_size = [window_half_size_one_dim, window_half_size_one_dim]
    return smooth_with_gaussian!(A, window_half_size)
end

"""
    smooth_with_gaussian!(A::Matrix{<:Number}, window_half_size)

Smooth the Matrix `A` by smoothing using a 2D Gaussian filter of size `2 * window_half_size + 1`.
"""
function smooth_with_gaussian!(A::Matrix{<:Number}, window_half_size; timer=TimerOutput("smoothing"))

    σ = 2 * window_half_size[1]

    @timeit_debug timer "create 1d gauss" begin
        gaus = gaussian(σ, 2 * σ + 1)
    end
    @timeit_debug timer "imfilter" begin
        imfilter!(A, A, kernelfactors((gaus, gaus)), "symmetric")
    end
    return A
end

function create_sparse_flow_from_sparse(flow_vectors, inds, flow_size)
    flow = zeros(flow_size) .+ zeros(flow_size) .* im
    for (vec, ind) in zip(flow_vectors, inds)
        flow[ind] .= vec
    end
    return flow
end

function create_sparse_flow_from_full(full_flow, inds)
    sparse_flow = zeros(size(full_flow)) .+ zeros(size(full_flow)) .* im

    for ind in inds
        sparse_flow[ind] .= full_flow[ind]
    end
    return sparse_flow
end

"""
    function angle_rmse(x, y)

Calculate the root mean square error in angle between `x` and `y`. Output in degrees.
"""
function angle_rmse(x, y)
    # @assert eltype(x) == eltype(y)
    @assert eltype(x) <: Complex

    return sqrt(mse(rad2deg.(angle.(x)), rad2deg.(angle.(y))))
end

"""
    function angle_mae(x, y)

Calculate the mean absolute error in angle between `x` and `y`. Output in degrees.
"""
function angle_mae(x, y)
    return mean(abs.(rad2deg.(angle.(x)) - rad2deg.(angle.(y))))
end

"""
    function vec_len(x)

Calculate the lenght of vector `x`. `x` is a complex number.
"""
function vec_len(x)
    return sqrt(real(x)^2 + imag(x)^2)
end

"""
    function mean(x)

Calculate the mean of `x`.
"""
function mean(x)
    sum(x)/length(x)
end

# TODO: sppeeudp
"""
    inds_to_points(inds::Array{CartesianIndex{2}, 1})

Transform an array of `CartersianIndexes` to an array of where each column is a vector of indices of the input array.
"""
function inds_to_points(inds::Array{CartesianIndex{2}, 1})
    pos_x = [ind[1] for ind in inds]
    pos_y = [ind[2] for ind in inds]
    return collect(transpose([pos_x pos_y]))
end

"""
    points_to_inds(points::Array{Int64,2})

Reverse of the `inds_to_points` function.
"""
function points_to_inds(points)
    return [CartesianIndex(points[1, k], points[2, k]) for k in 1:size(points,2)]
end

"""
    max_displacement(flow)

Find the maximum displacement of `flow`, ignoring `NaN`s.
"""
function max_displacement(flow)
    magnitudes = map(x -> vec_len(x), flow)
    max_mag = maximum(filter(!isnan, magnitudes))
    return max_mag
end

"""
    lap(args; kwargs)

Return a 2D `Flow` matrix of displacements that transforms `source` closer to `target` and `source` image transformed with this displacement field, `source_reg`.

# Details
Perform the classic Local All-Pass algorithm with post-proccessing (inpainting and smoothing), then warps the `source` image with the resulting estimate flow.

# Arguments:
- `target::Image`: target/fixed grayscale image.
- `source::Image`: source/moving grayscale image.
- `fhs::Integer`: the half size of the base of the gaussian filters used.
- `window_size=(2fhs+1, 2fhs+1)`: the size of the local window (tuple of 2 ints), usually same as filter_size.

# Keyword Arguments:
- `timer::TimerOutput=TimerOutput("lap")`: provide a timer which times certain blocks in the function.
- `display::Bool=false`: verbose and debug prints.

# Outputs:
- `flow_estim`: full estimated flow.
- `source_reg`: `source` image warped by `flow_estim`.

See also: [`single_lap`](@ref), [`inpaint_nans!`](@ref), [`smooth_with_gaussian!`](@ref), [`warp_img`](@ref).
"""
function lap(target::Image,
             source::Image,
             fhs,
             window_size=(2fhs+1, 2fhs+1);
             timer::TimerOutput=TimerOutput("lap"),
             display::Bool=false)

    @timeit_debug timer "lap" begin
        flow_estim = single_lap(target, source, fhs, window_size, timer=timer, display=display) # , new_feature=new_feature
    end
    @timeit_debug timer "inpainting" begin
        @timeit_debug timer "replicating borders" begin
            window_half_size = Int64.((window_size.-(1,1))./2)
            middle_vals = flow_estim[window_half_size[1]+1:end-window_half_size[1],
                                     window_half_size[2]+1:end-window_half_size[2]]
            flow_estim = parent(padarray(middle_vals, Pad(:replicate, window_half_size...)))
        end # "replication borders"
        @timeit_debug timer "inside" begin
            inpaint_nans!(flow_estim)
        end # "inside"
    end
    @timeit_debug timer "smoothing" begin
        smooth_with_gaussian!(flow_estim, window_half_size, timer=timer)
    end
    @timeit_debug timer "generate source_reg" begin
        source_reg = warp_img(source, -real(flow_estim), -imag(flow_estim))
        # source_reg = warp_img(source, -real(flow_estim), -imag(flow_estim), target)
    end
    # if display
    #     print_timer(timer)
    # end
    return flow_estim, source_reg
end

"""
    sparse_lap(args; kwargs)

Return a 2D `Flow` matrix of displacements that transforms `source` closer to `target`.

# Details
Uses the Local All-Pass algorithm at indices, `inds`, of the target image, with a high gradient magnitude, to estimate displacement vectors at these locations.
These estimated vectors are then fit into a parametric model (`flow_interpolation_method`) and a full displacemnt field/flow, or `full_flow_estim`, is made.
This flow is then used to warp the `source` image, creating the `source_reg` image.

# Arguments:
- `target::Image`: target/fixed grayscale image.
- `source::Image`: source/moving grayscale image.
- `fhs::Integer`: the half size of the base of the gaussian filters used.
- `window_size=(2fhs+1, 2fhs+1)`: the size of the local window (tuple of 2 ints), usually same as filter_size.

# Keyword Arguments:
- `spacing::Integer=25`: the smallest number of pixels that can separate two `inds`. See also: [`find_edge_points`](@ref)
- `point_count::Integer=100`: the number of `inds` we are looking for.
(The method will use the maximum amount up to this number if the this amount of `inds` isn't found.) See also: [`find_edge_points`](@ref)
- `timer::TimerOutput=TimerOutput("sparse_lap")`: provide a timer which times certain blocks in the function.
- `display::Bool=false`: verbose and debug prints.
- `flow_interpolation_method::Symbol=:quad`: choose which strategy to use for sparse flow interpolation.
(Choices: `:quad` -> fits to a global quadratic model, `:rbf` -> fits to a local rbf model.) See also: [`interpolate_flow`](@ref)
- `base_method_kwargs=Dict(:timer => timer, :display => display))`: keyword arguments passed to the base method.
(In this case: [`single_lap_at_points`](@ref)).

# Outputs:
- `full_flow_estim`: full estimated flow.
- `source_reg`: `source` image warped by `full_flow_estim`. See also: [`warp_img`](@ref).
- `flow_estim_at_inds`: estimated displacement vectors, at the indices `inds`. These are fitted into the `flow_interpolation_method` to
make the full estimated flow `full_flow_estim`.
- `inds`: indices chosen for the high gradient at which the `flow_estim_at_inds` is estimated.

See also: [`showflow`](@ref), [`imgshow`](@ref), [`find_edge_points`](@ref), [`interpolate_flow`](@ref)
"""
function sparse_lap(target,
                    source,
                    fhs,
                    window_size=(2fhs+1, 2fhs+1);
                    spacing::Integer=15,
                    point_count::Integer=200,
                    timer::TimerOutput=TimerOutput("sparse_lap"),
                    display::Bool=false,
                    flow_interpolation_method::Symbol=:quad,
                    base_method_kwargs=Dict(:timer => timer, :display => display))

    mask = parent(padarray(trues(size(target).-(2*fhs, 2*fhs)), Fill(false, (fhs, fhs), (fhs, fhs))))
    @timeit_debug timer "find edge points" begin
        inds = find_edge_points(target, spacing=spacing, point_count=point_count, mask=mask)
    end
    @timeit_debug timer "sparse lap" begin
        flow_estim_at_inds, inds = single_lap_at_points(target, source, fhs, window_size, inds; base_method_kwargs...)
    end
    if all(isnan, flow_estim_at_inds)
        @timeit_debug timer "interpolate flow" begin
            full_flow_estim = zeros(size(target)) .* im .+ zeros(size(target))
        end
    else
        @timeit_debug timer "interpolate flow" begin
            full_flow_estim = interpolate_flow(flow_estim_at_inds, inds, size(target), method=flow_interpolation_method)
        end
    end
    @timeit_debug timer "generate source_reg" begin
        # source_reg = warp_img(source, -real(full_flow_estim), -imag(full_flow_estim))
        source_reg = warp_img(source, -real(full_flow_estim), -imag(full_flow_estim), target)
    end
    return full_flow_estim, source_reg, flow_estim_at_inds, inds
end


text:
type RegularizedLTIModel <: DischargeModel
    M::Int
    λ::Float64
    β::Vector{Float64}
    H::Matrix{Float64}
    Q::Vector{Float64}
end

RegularizedLTIModel(M,N,λ=0.0) = RegularizedLTIModel(M,λ,zeros(M),zeros(M,N),zeros(N))

function StatsBase.fit(M::RegularizedLTIModel,H::Matrix{Float64},Q::Vector{Float64})
    U,s,V = svd(H')
    D = diagm(s./(s.^2+M.λ^2))
    RegularizedLTIModel(M.M,M.λ,V*D*U'Q,H,Q)
end

function StatsBase.fit!(M::RegularizedLTIModel,H::Matrix{Float64},Q::Vector{Float64})
    U,s,V = svd(H')
    D = diagm(s./(s.^2+M.λ^2))
    M.β = V*D*U'Q
    M.H = H
    M.Q = Q
    M
end

StatsBase.coef(M::RegularizedLTIModel) = M.β
StatsBase.nobs(M::RegularizedLTIModel) = length(M.Q)
StatsBase.model_response(M::RegularizedLTIModel) = M.Q
StatsBase.predict(M::RegularizedLTIModel) = M.H'M.β
StatsBase.predict(M::RegularizedLTIModel,H) = H'M.β
StatsBase.residuals(M::RegularizedLTIModel) = M.Q - predict(M)
StatsBase.residuals(M::RegularizedLTIModel,H,Q) = Q-predict(M,H)


text:
using Test
using LinearAlgebra
using GenericLinearAlgebra
using GenericLinearAlgebra: qrBlocked!

@testset "The QR decomposition" begin
    @testset "Problem dimension ($m,$n) with block size $bz" for
        (m, n) in (( 10,  5), ( 10,  10), ( 5,  10),
                   (100, 50), (100, 100), (50, 100)),
            bz in (1, 2, 3, 4, 7, 8, 9, 15, 16, 17, 31, 32, 33)

        A = randn(m, n)
        Aqr = qrBlocked!(copy(A), bz)
        AqrQ = Aqr[Tuple{:QBlocked}]
        if m >= n
            @test (AqrQ'A)[1:min(m,n),:] ≈ Aqr[Tuple{:R}]
        else # For type stbility getindex(,Tuple{:R}) throw when the output is trapezoid
            @test (AqrQ'A) ≈ triu(Aqr.factors)
        end
        @test AqrQ'*(AqrQ*A) ≈ A
    end

    @testset "Error paths" begin
        @test_throws DimensionMismatch LinearAlgebra.reflectorApply!(zeros(5, 5), zeros(4), 1.0)
        @test_throws ArgumentError     qrBlocked!(randn(5, 10))[Tuple{:R}]
    end
end


text:
__precompile__()

module LSODA

using Compat, DiffEqBase
import DiffEqBase: solve
using LinearAlgebra, Printf

const warnkeywords =
    (:save_idxs, :d_discontinuities, :isoutofdomain, :unstable_check,
     :calck, :progress, :timeseries_steps, :dense,
     :dtmin, :dtmax,
     :internalnorm, :gamma, :beta1, :beta2, :qmax, :qmin, :qoldinit)

function __init__()
    global warnlist = Set(warnkeywords)
end

abstract type LSODAAlgorithm <: DiffEqBase.AbstractODEAlgorithm end
abstract type LSODAAlgorithm <: DiffEqBase.AbstractODEAlgorithm end
struct lsoda <: LSODAAlgorithm end

const depsfile = joinpath(dirname(dirname(@__FILE__)),"deps","deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("LSODA is not properly installed. Please run Pkg.build(\"LSODA\")")
end

export lsoda, lsoda_0, lsoda_opt_t, lsoda_context_t, lsoda_prepare, lsoda_reset, lsoda_opt_t, lsoda_free, lsoda_evolve!, UserFunctionAndData

export LSODAAlgorithm, solve

include("types_and_consts.jl")
include("handle.jl")
include("solver.jl")
include("common.jl")

end # module


text:
using BijectiveHilbert
using Documenter

makedocs(;
    modules=[BijectiveHilbert],
    authors="Andrew Dolgert <adolgert@uw.edu>",
    repo="https://github.com/adolgert/BijectiveHilbert.jl/blob/{commit}{path}#L{line}",
    sitename="BijectiveHilbert.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://adolgert.github.io/BijectiveHilbert.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Usage" => "usage.md",
        "Algorithms" => [
            "hilbert.md",
            "simple2d.md",
            "globalgray.md",
            "compact.md",
            "facecontinuous.md"
        ]
    ],
)

deploydocs(;
    devbranch = "main",
    repo="github.com/adolgert/BijectiveHilbert.jl",
    deploy_config=Documenter.GitHubActions()
)


text:
using CPIMC
using CPIMC.PlaneWaves

import CPIMC: energy
import CPIMC.UniformElectronGas: λ, β

import LinearAlgebra: dot

struct FreeElectronGas <: Model end

energy(::FreeElectronGas, o::PlaneWave) = dot(o.vec, o.vec)

function move_particle(m::Model, e::Ensemble, c::Configuration)
    x = rand(c.occupations)
    oe = setdiff!(sphere_with_same_spin(x), c.occupations)

    if isempty(oe)
        return 1.0, Step()
    end
    y = rand(oe)

    # weight change
    dw = exp(-e.β*(energy(m,y)-energy(m,x)))

    # get orbitals for reverse update
    oe2 = setdiff!(sphere_with_same_spin(y), c.occupations)

    # quotient of proposal probabilities
    δv = length(oe)/length(oe2)

    δv * dw, Step((x,),(y,))
end

function Ekin(m::Model, e::Ensemble, c::Configuration) :: UInt
    sum(energy(m, n) for n in c.occupations)
end

function occupations(m::Model, e::Ensemble, c::Configuration, emax::Int=100) :: Array{UInt,1}
    nk = zeros(UInt, emax)

    ens = [ energy(m, n) for n in c.occupations ]

    for ε in ens[ens .< emax]
        nk[ε+1] = nk[ε+1] + 1
    end
    nk
end

function main()
    # MC options
    NMC = 10^6
    cyc = 10
    NEquil = 10^4
    # system parameters
    θ = 1.0
    rs = 0.5

    S = sphere_with_same_spin(PlaneWave((0,0,0),Up),dk=2)### use 33 particles
    N = length(S)
    ξ = fractional_spin_polarization(S)
    c = Configuration(S)
    d = dimension(c.occupations)

    println("θ: ", θ)
    println("rs: ", rs)
    println("N: ", N)
    println("d: ", d)

    e = CEnsemble(λ(N,rs,d), β(θ,N,ξ,d), N) # β only works for 3D

    updates = [ move_particle ]

    measurements = Dict(
      :Ekin => (Variance(UInt), Ekin)
    , :occs => (Group([Variance(UInt) for i in 1:100]), occupations))

    println("Start MC process ... ")
    d = sweep!(FreeElectronGas(), e, c, updates, measurements, NMC, cyc, NEquil)
    println(" finished.")

    print(d)

    print_results(measurements, e)

end

main()


text:
using Flux, HDF5, ProgressMeter
using MLDataUtils, BSON
using Statistics, Random
using Base.Iterators: partition
#regression neural net
function regressionNN(dimIn, dimOut, neuronas)
    Chain(
        Dense(dimIn, neuronas, Flux.relu),
        Dense(neuronas, neuronas, Flux.relu),
        Dense(neuronas, neuronas, Flux.relu),
        Dense(neuronas, neuronas, Flux.relu),
        Dense(neuronas, neuronas, Flux.relu),
        #Dense(neuronas, neuronas, Flux.relu),
        Dense(neuronas, neuronas, Flux.relu),
        Dense(neuronas, dimOut)) # |> gpu
end
function machine(xdata, ydata, xval, yval, seed;
        neurons=40, nepochs=100, bs=512, best_val = 5.0,
        last_improvement = 0,patience=5, gpu = false)
    #saving best model
    dict_models = Dict()
    dict_models["model_bm$(seed)"] = NaN
    #optimizer and loss function
    Random.seed!(seed)
    model = regressionNN(size(xdata)[1], size(ydata)[1], neurons)
    # if training on gpu's is an option
    model = gpu == true ? model |> gpu : model
    xdata = gpu == true ? xdata |> gpu : xdata
    ydata = gpu == true ? ydata |> gpu : ydata

    opt = ADAM()
    loss(x,y) = Flux.mse(model(x), y)
    #nepochs = 100
    lossVal = zeros(nepochs)
    #mini-batch
    #batch_size 512 is ok too, similar performance at test accuracy
    mb_idxs = partition(1:size(xdata)[2], bs)
    #xs, ys = shuffleobs((Float32.(xdatacoef), ycdata))
    p = Progress(nepochs, dt=1, barglyphs=BarGlyphs("[=> ]"), color = :yellow)
    for epoch_indx in 1:nepochs
        xs, ys = shuffleobs((Float32.(xdata), Float32.(ydata)))
        train_set = [(xs[:,p], ys[:,p]) for p in mb_idxs]
        Flux.train!(loss, params(model), train_set, opt)

        validation_loss = loss(xval, yval)
        lossVal[epoch_indx] = validation_loss
        # If this is the best val_loss we've seen so far, save the model out
        if validation_loss <= best_val
            #@info(" -> New best val_loss! Saving model out to model_bm$(seed).bson")
            last_improvement = epoch_indx
            dict_models["model_bm$(seed)"] = model # model_bestmodel
            best_val = validation_loss
        end
        if epoch_indx - last_improvement >= patience
            #@info(" -> Early-exiting iteration $(r) and epoch $(epoch_indx): no more patience")
            break
        end
        next!(p; showvalues = [(:epoch,epoch_indx), (:Model_seed, seed)])
        #ProgressMeter.ijulia_behavior(:append)
        ProgressMeter.ijulia_behavior(:clear)
    end
    model = dict_models["model_bm$(seed)"]
    #BSON.@save "models_bm$(seed).bson" dict_models lossVal
    model, lossVal
end

function ensembleMachines(xtrain, ytrain, xval, yval;
        nmodels = 10, nepoch = 100, neurons=40, bs = 512, best_val = 5.0,
        last_improvement = 0, patience=10, gpu = false)

    losses = zeros(nepoch, nmodels)
    machines = []
    dict_models = Dict()
    #pro = Progress(nmodels, dt=0.5, barglyphs=BarGlyphs("[=> ]"), color = :black)
    for seed in 1:nmodels
        modelo1, lossVal = machine(xtrain, ytrain, xval, yval,seed;
            neurons= neurons, nepochs= nepoch, bs=bs, best_val = best_val,
            last_improvement = last_improvement, patience=patience, gpu = gpu)

        losses[:,seed] = lossVal
        dict_models["model_bm$(seed)"] = modelo1
        #next!(pro; showvalues = [(:Model, seed)])
        #ProgressMeter.ijulia_behavior(:append)
        #ProgressMeter.ijulia_behavior(:clear)
    end
    #BSON.@save "models_bm$(nmodels).bson" dict_models losses
    dict_models, losses
end

function splitdata(xdata, ydata)
    #μnoi = mean(xdata)
    #stdnoi = std(xdata)
    #μref = mean(ydata_1)
    #stdref = std(ydata_1)
    #xdata = (xdata .- μnoi)./stdnoi
    #ydata = (ydata .- μref)./stdref

    xtrain = xdata[:,1:32000]
    ytrain = ydata[:,1:32000]
    xval = xdata[:,32001:36000]
    yval = ydata[:,32001:36000]
    xtest = xdata[:,36001:end]
    ytest = ydata[:,36001:end]
    xtrain, ytrain, xval, yval, xtest, ytest #, [μnoi, stdnoi, μref, stdref]
end

max_stars_repo_name: lazarusA/LightMatterInteraction

max_issues_repo_name: lazarusA/LightMatterInteraction

max_forks_repo_name: lazarusA/LightMatterInteraction


text:
@testset "Affine Parametric Array construction" begin
    @linvars x y z

    vars = [x, y, z]

    A = AffineParametricArray([x+1 y+2;x+y+z+1 x-z], vars)
    A1 = AffineParametricArray([1 2;3 4], vars)
    B = AffineParametricArray([x+1 y+2;x+y+z+1 x-z+1], vars)

    @test A.vars == vars
    @test A.coeffs == [[1 0;1 1], [0 1;1 0], [0 0;1 -1], [1 2;1 0]]
    @test A1.coeffs == [[0 0;0 0], [0 0;0 0], [0 0;0 0], [1 2;3 4]]

    @test A[1, 2] == y + 2
    @test A[:, 1] == [x+1, x+y+z+1]

    A1[1, 2] = x
    @test A1 == AffineParametricArray([1 x;3 4], vars)
    A1[:, 1] = [x+1, z-1]
    @test A1 == AffineParametricArray([x+1 x;z-1 4], vars)

    @test A() == [x+1 y+2;x+y+z+1 x-z]
    @test A([1..2, 2..3, 3..4]) == [2..3 4..5; 7..10 -3.. -1]
end

@testset "Affine parametric array operations" begin
    @linvars x y z
    vars = [x, y, z]

    A = AffineParametricArray([x+1 y+2;x+y+z+1 x-z], vars)
    B = AffineParametricArray([x+1 y+2;x+y+z+1 x-z+1], vars)

    @test A == A
    @test A != B

    @test +A == A
    @test -A == AffineParametricArray([-x-1 -y-2; -x-y-z-1 -x+z], vars)
    @test A + B == AffineParametricArray([2x+2 2y+4; 2x+2y+2z+2 2x-2z+1], vars)
    @test A - B == AffineParametricArray([0 0;0 -1], vars)

    C = [2 1;1 1]
    @test A + C == AffineParametricArray([x+3 y+3;x+y+z+2 x-z+1], vars)
    @test C + A == AffineParametricArray([x+3 y+3;x+y+z+2 x-z+1], vars)

    @test A * C == AffineParametricArray([2x+y+4 x+y+3;3x+2y+z+2 2x+y+1], vars)
    @test C * A == AffineParametricArray([3x+y+z+3 x+2y-z+4;2x+y+z+2 x+y-z+2], vars)
    @test C \ A == inv(C) * A

    @test A * 2 == AffineParametricArray([2x+2 2y+4;2x+2y+2z+2 2x-2z], vars)
    @test 2 * A == AffineParametricArray([2x+2 2y+4;2x+2y+2z+2 2x-2z], vars)
    @test A / 2 == AffineParametricArray([0.5x+0.5 0.5y+1;0.5x+0.5y+0.5z+0.5 0.5x-0.5z], vars)
end

max_stars_repo_name: lucaferranti/IntervalLinearAlgebra

max_issues_repo_name: lucaferranti/IntervalLinearAlgebra

max_forks_repo_name: lucaferranti/IntervalLinearAlgebra


text:
module BetweenFlags

include("types.jl")
include("tokenize.jl")

end


text:
"""
    $(SIGNATURES)
Create a new CloudGraphs-based DFG factor graph using a Neo4j.Connection.
"""
function CloudGraphsDFG{T}(neo4jConnection::Neo4j.Connection, userId::String, robotId::String, sessionId::String, encodePackedTypeFunc, getPackedTypeFunc, decodePackedTypeFunc, rebuildFactorMetadata!; description::String="CloudGraphs DFG", solverParams::T=NoSolverParams(), useCache::Bool=false) where T <: AbstractParams
    graph = Neo4j.getgraph(neo4jConnection)
    neo4jInstance = Neo4jInstance(neo4jConnection, graph)
    return CloudGraphsDFG{T}(neo4jInstance, description, userId, robotId, sessionId, encodePackedTypeFunc, getPackedTypeFunc, decodePackedTypeFunc, rebuildFactorMetadata!, Dict{Symbol, Int64}(), Dict{Symbol, DFGVariable}(), Dict{Symbol, DFGFactor}(), Symbol[], solverParams, useCache)
end
"""
    $(SIGNATURES)
Create a new CloudGraphs-based DFG factor graph by specifying the Neo4j connection information.
"""
function CloudGraphsDFG{T}(host::String, port::Int, dbUser::String, dbPassword::String, userId::String, robotId::String, sessionId::String, encodePackedTypeFunc, getPackedTypeFunc, decodePackedTypeFunc, rebuildFactorMetadata!; description::String="CloudGraphs DFG", solverParams::T=NoSolverParams(), useCache::Bool=false) where T <: AbstractParams
    neo4jConnection = Neo4j.Connection(host, port=port, user=dbUser, password=dbPassword);
    return CloudGraphsDFG{T}(neo4jConnection, userId, robotId, sessionId, encodePackedTypeFunc, getPackedTypeFunc, decodePackedTypeFunc, rebuildFactorMetadata!, description=description, solverParams=solverParams, useCache=useCache)
end

"""
    $(SIGNATURES)
Gets an empty and unique CloudGraphsDFG derived from an existing DFG.
"""
function _getDuplicatedEmptyDFG(dfg::CloudGraphsDFG)::CloudGraphsDFG
    count = 0
    sessionId = dfg.sessionId*"_$count"
    while true #do..while loop
        count += 1
        sessionId = dfg.sessionId*"_$count"
        length(_getLabelsFromCyphonQuery(dfg.neo4jInstance, "(node:$(dfg.userId):$(dfg.robotId):$(sessionId))")) == 0 && break
    end
    @debug "Unique+empty copy session name: $sessionId"
    return CloudGraphsDFG{typeof(dfg.solverParams)}(dfg.neo4jInstance.connection, dfg.userId, dfg.robotId, sessionId, dfg.encodePackedTypeFunc, dfg.getPackedTypeFunc, dfg.decodePackedTypeFunc, dfg.rebuildFactorMetadata!, solverParams=deepcopy(dfg.solverParams), description="(Copy of) $(dfg.description)", useCache=dfg.useCache)
end

# Accessors
getLabelDict(dfg::CloudGraphsDFG) = dfg.labelDict
getDescription(dfg::CloudGraphsDFG) = dfg.description
setDescription(dfg::CloudGraphsDFG, description::String) = dfg.description = description
getAddHistory(dfg::CloudGraphsDFG) = dfg.addHistory
getSolverParams(dfg::CloudGraphsDFG) = dfg.solverParams
function setSolverParams(dfg::CloudGraphsDFG, solverParams::T)::T where T <: AbstractParams
    return dfg.solverParams = solverParams
end

function getSerializationModule(dfg::CloudGraphsDFG)::Module where G <: AbstractDFG
    # TODO: If we need to specialize this for RoME etc, here is where we'll change it.
    return Main
end

"""
    $(SIGNATURES)
True if the variable or factor exists in the graph.
"""
function exists(dfg::CloudGraphsDFG, nId::Symbol)
    # If in the dictionary, then shortcut return true
    dfg.useCache && haskey(dfg.labelDict, nId) && return true
    # Otherwise try get it
    nodeId = _tryGetNeoNodeIdFromNodeLabel(dfg.neo4jInstance, dfg.userId, dfg.robotId, dfg.sessionId, nId)
    if nodeId != nothing
        push!(dfg.labelDict, nId=>nodeId)
        return true
    end
    return false
end
function exists(dfg::CloudGraphsDFG, node::N) where N <: DFGNode
    return exists(dfg, node.label)
end

"""
    $SIGNATURES

Return whether `sym::Symbol` represents a variable vertex in the graph.
"""
isVariable(dfg::CloudGraphsDFG, sym::Symbol)::Bool =
    "VARIABLE" in _getNodeTags(dfg.neo4jInstance, [dfg.userId, dfg.robotId, dfg.sessionId, String(sym)])

"""
    $SIGNATURES

Return whether `sym::Symbol` represents a factor vertex in the graph.
"""
isFactor(dfg::CloudGraphsDFG, sym::Symbol)::Bool =
    "FACTOR" in _getNodeTags(dfg.neo4jInstance, [dfg.userId, dfg.robotId, dfg.sessionId, String(sym)])

"""
    $(SIGNATURES)
Add a DFGVariable to a DFG.
"""
function addVariable!(dfg::CloudGraphsDFG, variable::DFGVariable)::Bool
    if exists(dfg, variable)
        error("Variable '$(variable.label)' already exists in the factor graph")
    end
    props = Dict{String, Any}()
    props["label"] = string(variable.label)
    props["timestamp"] = string(variable.timestamp)
    props["tags"] = JSON2.write(variable.tags)
    props["estimateDict"] = JSON2.write(variable.estimateDict)
    props["solverDataDict"] = JSON2.write(Dict(keys(variable.solverDataDict) .=> map(vnd -> pack(dfg, vnd), values(variable.solverDataDict))))
    props["smallData"] = JSON2.write(variable.smallData)
    props["ready"] = variable.ready
    props["backendset"] = variable.backendset
    # Don't handle big data at the moment.

    neo4jNode = Neo4j.createnode(dfg.neo4jInstance.graph, props);
    variable._internalId = neo4jNode.id
    Neo4j.updatenodelabels(neo4jNode, union([string(variable.label), "VARIABLE", dfg.userId, dfg.robotId, dfg.sessionId], variable.tags))

    # Make sure that if there exists a SESSION sentinel that it is attached.
    # TODO: Optimize this.
    _bindSessionNodeToInitialVariable(dfg.neo4jInstance, dfg.userId, dfg.robotId, dfg.sessionId, string(variable.label))

    # Update our internal dictionaries.
    push!(dfg.labelDict, variable.label=>variable._internalId)
    push!(dfg.variableCache, variable.label=>variable)
    # Track insertion
    push!(dfg.addHistory, variable.label)

    return true
end

"""
    $(SIGNATURES)
Add a DFGFactor to a DFG.
"""
function addFactor!(dfg::CloudGraphsDFG, variables::Vector{DFGVariable}, factor::DFGFactor)::Bool
    if exists(dfg, factor)
        error("Factor '$(factor.label)' already exists in the factor graph")
    end

    # Update the variable ordering
    factor._variableOrderSymbols = map(v->v.label, variables)

    # Construct the properties to save
    props = Dict{String, Any}()
    props["label"] = string(factor.label)
    props["tags"] = JSON2.write(factor.tags)
    # Pack the node data
    fnctype = factor.data.fnc.usrfnc!
    packtype = getfield(_getmodule(fnctype), Symbol("Packed$(_getname(fnctype))"))
    packed = convert(PackedFunctionNodeData{packtype}, factor.data)
    props["data"] = JSON2.write(packed)
    # Include the type
    props["fnctype"] = String(_getname(fnctype))
    props["_variableOrderSymbols"] = JSON2.write(factor._variableOrderSymbols)
    props["backendset"] = factor.backendset
    props["ready"] = factor.ready
    # Don't handle big data at the moment.

    neo4jNode = Neo4j.createnode(dfg.neo4jInstance.graph, props);
    factor._internalId = neo4jNode.id
    Neo4j.updatenodelabels(neo4jNode, union([string(factor.label), "FACTOR", dfg.userId, dfg.robotId, dfg.sessionId], factor.tags))

    # Add all the relationships - get them to cache them + make sure the links are correct
    for variable in variables
        v = getVariable(dfg, variable.label)
        vNode = Neo4j.getnode(dfg.neo4jInstance.graph, v._internalId)
        Neo4j.createrel(neo4jNode, vNode, "FACTORGRAPH")
    end

    # Graphs.add_vertex!(dfg.g, v)
    push!(dfg.labelDict, factor.label=>factor._internalId)
    push!(dfg.factorCache, factor.label=>factor)
    # Track insertion only for variables
    # push!(dfg.addHistory, factor.label

    return true
end

"""
    $(SIGNATURES)
Add a DFGFactor to a DFG.
"""
function addFactor!(dfg::CloudGraphsDFG, variableIds::Vector{Symbol}, factor::DFGFactor)::Bool
    variables = map(vId -> getVariable(dfg, vId), variableIds)
    return addFactor!(dfg, variables, factor)
end

"""
    $(SIGNATURES)
Get a DFGVariable from a DFG using its underlying integer ID.
"""
function getVariable(dfg::CloudGraphsDFG, variableId::Int64)::DFGVariable
    props = getnodeproperties(dfg.neo4jInstance.graph, variableId)
    # Time to do deserialization
    # props["label"] = Symbol(variable.label)
    timestamp = DateTime(props["timestamp"])
    tags =  JSON2.read(props["tags"], Vector{Symbol})
    #TODO this will work for some time, but unpacking in an <: AbstractPointParametricEst would be lekker.
    estimateDict = JSON2.read(props["estimateDict"], Dict{Symbol, MeanMaxPPE})
    smallData = nothing
    smallData = JSON2.read(props["smallData"], Dict{String, String})

    packed = JSON2.read(props["solverDataDict"], Dict{String, PackedVariableNodeData})
    solverData = Dict(Symbol.(keys(packed)) .=> map(p -> unpack(dfg, p), values(packed)))

    # Rebuild DFGVariable
    variable = DFGVariable(Symbol(props["label"]), variableId)
    variable.timestamp = timestamp
    variable.tags = tags
    variable.estimateDict = estimateDict
    variable.solverDataDict = solverData
    variable.smallData = smallData
    variable.ready = props["ready"]
    variable.backendset = props["backendset"]

    # Add to cache
    push!(dfg.variableCache, variable.label=>variable)

    return variable
end

"""
    $(SIGNATURES)
Get a DFGVariable from a DFG using its label.
"""
function getVariable(dfg::CloudGraphsDFG, label::Union{Symbol, String})::DFGVariable
    if typeof(label) == String
        label = Symbol(label)
    end
    dfg.useCache && haskey(dfg.variableCache, label) && return dfg.variableCache[label]
    # Else try get it
    nodeId = _tryGetNeoNodeIdFromNodeLabel(dfg.neo4jInstance, dfg.userId, dfg.robotId, dfg.sessionId, label)
    if nodeId == nothing
        error("Unable to retrieve the ID for variable '$label'. Please check your connection to the database and that the variable exists.")
    end

    return getVariable(dfg, nodeId)
end

"""
    $(SIGNATURES)
Get a DFGFactor from a DFG using its underlying integer ID.
"""
function getFactor(dfg::CloudGraphsDFG, factorId::Int64)::DFGFactor
    props = getnodeproperties(dfg.neo4jInstance.graph, factorId)

    label = props["label"]
    tags = JSON2.read(props["tags"], Vector{Symbol})

    data = props["data"]
    datatype = props["fnctype"]
    # fulltype = getfield(Main, Symbol(datatype))
    packtype = getfield(Main, Symbol("Packed"*datatype))
    packed = JSON2.read(data, GenericFunctionNodeData{packtype,String})
    fullFactor = dfg.decodePackedTypeFunc(dfg, packed)

    # Include the type
    _variableOrderSymbols = JSON2.read(props["_variableOrderSymbols"], Vector{Symbol})
    backendset = props["backendset"]
    ready = props["ready"]

    # Rebuild DFGVariable
    factor = DFGFactor{typeof(fullFactor.fnc), Symbol}(Symbol(label), factorId)
    factor.tags = tags
    factor.data = fullFactor
    factor._variableOrderSymbols = _variableOrderSymbols
    factor.ready = ready
    factor.backendset = backendset

    # Lastly, rebuild the metadata
    factor = dfg.rebuildFactorMetadata!(dfg, factor)
    # GUARANTEED never to bite us in the future...
    # ... TODO: refactor if changed: https://github.com/JuliaRobotics/IncrementalInference.jl/issues/350
    solverData(factor).fncargvID = _variableOrderSymbols

    # Add to cache
    push!(dfg.factorCache, factor.label=>factor)

    return factor
end

"""
    $(SIGNATURES)
Get a DFGFactor from a DFG using its label.
"""
function getFactor(dfg::CloudGraphsDFG, label::Union{Symbol, String})::DFGFactor
    if typeof(label) == String
        label = Symbol(label)
    end
    dfg.useCache && haskey(dfg.factorCache, label) && return dfg.factorCache[label]
    # Else try get it
    nodeId = _tryGetNeoNodeIdFromNodeLabel(dfg.neo4jInstance, dfg.userId, dfg.robotId, dfg.sessionId, label)
    if nodeId == nothing
        error("Unable to retrieve the ID for factor '$label'. Please check your connection to the database and that the factor exists.")
    end

    return getFactor(dfg, nodeId)
end

"""
    $(SIGNATURES)
Update a complete DFGVariable in the DFG.
"""
function updateVariable!(dfg::CloudGraphsDFG, variable::DFGVariable)::DFGVariable
    if !exists(dfg, variable)
        error("Variable label '$(variable.label)' does not exist in the factor graph")
    end
    nodeId = _tryGetNeoNodeIdFromNodeLabel(dfg.neo4jInstance, dfg.userId, dfg.robotId, dfg.sessionId, variable.label)
    # Update the node ID
    variable._internalId = nodeId

    neo4jNode = Neo4j.getnode(dfg.neo4jInstance.graph, nodeId)
    props = getnodeproperties(dfg.neo4jInstance.graph, nodeId)

    props["label"] = string(variable.label)
    props["timestamp"] = string(variable.timestamp)
    props["tags"] = JSON2.write(variable.tags)
    props["estimateDict"] = JSON2.write(variable.estimateDict)
    props["solverDataDict"] = JSON2.write(Dict(keys(variable.solverDataDict) .=> map(vnd -> pack(dfg, vnd), values(variable.solverDataDict))))
    props["smallData"] = JSON2.write(variable.smallData)
    props["ready"] = variable.ready
    props["backendset"] = variable.backendset
    # Don't handle big data at the moment.

    Neo4j.updatenodeproperties(neo4jNode, props)
    Neo4j.updatenodelabels(neo4jNode, union([string(variable.label), "VARIABLE", dfg.userId, dfg.robotId, dfg.sessionId], variable.tags))
    return variable
end

"""
    $(SIGNATURES)
Update solver and estimate data for a variable (variable can be from another graph).
"""
function mergeUpdateVariableSolverData!(dfg::CloudGraphsDFG, sourceVariable::DFGVariable)::DFGVariable
    if !exists(dfg, sourceVariable)
        error("Source variable '$(sourceVariable.label)' doesn't exist in the graph.")
    end
    var = getVariable(dfg, sourceVariable.label)
    newEsts = merge!(var.estimateDict, deepcopy(sourceVariable.estimateDict))
    newSolveData = merge!(var.solverDataDict, sourceVariable.solverDataDict)
    Neo4j.setnodeproperty(dfg.neo4jInstance.graph, var._internalId, "estimateDict",
        JSON2.write(newEsts))
    Neo4j.setnodeproperty(dfg.neo4jInstance.graph, var._internalId, "solverDataDict",
        JSON2.write(Dict(keys(newSolveData) .=> map(vnd -> pack(dfg, vnd), values(newSolveData)))))
    return sourceVariable
end

"""
    $(SIGNATURES)
Update a complete DFGFactor in the DFG.
"""
function updateFactor!(dfg::CloudGraphsDFG, factor::DFGFactor)::DFGFactor
    if !exists(dfg, factor)
        error("Factor label '$(factor.label)' does not exist in the factor graph")
    end
    nodeId = _tryGetNeoNodeIdFromNodeLabel(dfg.neo4jInstance, dfg.userId, dfg.robotId, dfg.sessionId, factor.label)
    # Update the _internalId
    factor._internalId = nodeId
    neo4jNode = Neo4j.getnode(dfg.neo4jInstance.graph, nodeId)
    props = getnodeproperties(dfg.neo4jInstance.graph, nodeId)

    props["label"] = string(factor.label)
    props["tags"] = JSON2.write(factor.tags)
    # Pack the node data
    fnctype = factor.data.fnc.usrfnc!
    packtype = getfield(_getmodule(fnctype), Symbol("Packed$(_getname(fnctype))"))
    packed = convert(PackedFunctionNodeData{packtype}, factor.data)
    props["data"] = JSON2.write(packed)
    # Include the type
    props["fnctype"] = String(_getname(fnctype))
    props["_variableOrderSymbols"] = JSON2.write(factor._variableOrderSymbols)
    props["backendset"] = factor.backendset
    props["ready"] = factor.ready
    # Don't handle big data at the moment.

    Neo4j.updatenodeproperties(neo4jNode, props)
    Neo4j.updatenodelabels(neo4jNode, union([string(factor.label), "FACTOR", dfg.userId, dfg.robotId, dfg.sessionId], factor.tags))

    return factor
end

"""
    $(SIGNATURES)
Update a complete DFGFactor in the DFG and update it's relationships.
"""
function updateFactor!(dfg::CloudGraphsDFG, variables::Vector{DFGVariable}, factor::DFGFactor)::DFGFactor
    # Update the body
    factor = updateFactor!(dfg, factor)

    # Now update the relationships
    existingNeighbors = getNeighbors(dfg, factor)
    if symdiff(existingNeighbors, map(v->v.label, variables)) == []
        # Done, otherwise we need to remake the edges.
        return factor
    end
    # Delete existing relationships
    fNode = Neo4j.getnode(dfg.neo4jInstance.graph, factor._internalId)
    for relationship in Neo4j.getrels(fNode)
        if relationship.reltype == "FACTORGRAPH"
            Neo4j.deleterel(relationship)
        end
    end
    # Make new ones
    for variable in variables
        v = getVariable(dfg, variable.label)
        vNode = Neo4j.getnode(dfg.neo4jInstance.graph, v._internalId)
        @info "Creating REL between factor $(fNode) and variable $(vNode)"
        Neo4j.createrel(fNode, vNode, "FACTORGRAPH")
    end

    return factor
end

"""
    $(SIGNATURES)
Update a complete DFGFactor in the DFG and update it's relationships.
"""
function updateFactor!(dfg::CloudGraphsDFG, variableIds::Vector{Symbol}, factor::DFGFactor)::DFGFactor
    variables = map(vId -> getVariable(dfg, vId), variableIds)
    return updateFactor!(dfg, variables, factor)
end

"""
    $(SIGNATURES)
Delete a DFGVariable from the DFG using its label.
"""
function deleteVariable!(dfg::CloudGraphsDFG, label::Symbol)::DFGVariable
    variable = nothing
    if dfg.useCache && haskey(dfg.variableCache, label)
        variable = dfg.variableCache[label]
    else
        # Else try get it
        variable = getVariable(dfg, label)
    end
    if variable == nothing
        error("Unable to retrieve the ID for variable '$label'. Please check your connection to the database and that the variable exists.")
    end

    # Perform detach+deletion
    _getNeoNodesFromCyphonQuery(dfg.neo4jInstance, "(node:VARIABLE) where id(node)=$(variable._internalId) detach delete node ")

    # Clearing history
    dfg.addHistory = symdiff(dfg.addHistory, [label])
    haskey(dfg.variableCache, label) && delete!(dfg.variableCache, label)
    haskey(dfg.labelDict, label) && delete!(dfg.labelDict, label)
    return variable
end

#Alias
"""
    $(SIGNATURES)
Delete a referenced DFGVariable from the DFG.
"""
deleteVariable!(dfg::CloudGraphsDFG, variable::DFGVariable)::DFGVariable = deleteVariable!(dfg, variable.label)

"""
    $(SIGNATURES)
Delete a DFGFactor from the DFG using its label.
"""
function deleteFactor!(dfg::CloudGraphsDFG, label::Symbol)::DFGFactor
    factor = nothing
    if dfg.useCache && haskey(dfg.factoreCache, label)
        factor = dfg.factorCache[label]
    else
        # Else try get it
        factor = getFactor(dfg, label)
    end
    if factor == nothing
        error("Unable to retrieve the ID for factor '$label'. Please check your connection to the database and that the factor exists.")
    end

    # Perform detach+deletion
    _getNeoNodesFromCyphonQuery(dfg.neo4jInstance, "(node:FACTOR) where id(node)=$(factor._internalId) detach delete node ")

    # Clearing history
    dfg.addHistory = symdiff(dfg.addHistory, [label])
    haskey(dfg.factorCache, label) && delete!(dfg.factorCache, label)
    haskey(dfg.labelDict, label) && delete!(dfg.labelDict, label)
    return factor
end

# Alias
"""
    $(SIGNATURES)
Delete the referened DFGFactor from the DFG.
"""
deleteFactor!(dfg::CloudGraphsDFG, factor::DFGFactor)::DFGFactor = deleteFactor!(dfg, factor.label)

"""
    $(SIGNATURES)
List the DFGVariables in the DFG.
Optionally specify a label regular expression to retrieves a subset of the variables.
"""
function getVariables(dfg::CloudGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[])::Vector{DFGVariable}
    variableIds = getVariableIds(dfg, regexFilter)
    # TODO: Optimize to use tags in query here!
    variables = map(vId->getVariable(dfg, vId), variableIds)
    if length(tags) > 0
        mask = map(v -> length(intersect(v.tags, tags)) > 0, variables )
        return variables[mask]
    end
    return variables
end

"""
    $(SIGNATURES)
Get a list of IDs of the DFGVariables in the DFG.
Optionally specify a label regular expression to retrieves a subset of the variables.
"""
function getVariableIds(dfg::CloudGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing)::Vector{Symbol}
    # Optimized for DB call
    if regexFilter == nothing
        return _getLabelsFromCyphonQuery(dfg.neo4jInstance, "(node:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId):VARIABLE)")
    else
        return _getLabelsFromCyphonQuery(dfg.neo4jInstance, "(node:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId):VARIABLE) where node.label =~ '$(regexFilter.pattern)'")
    end
end

# Alias
"""
    $(SIGNATURES)
List the DFGVariables in the DFG.
Optionally specify a label regular expression to retrieves a subset of the variables.
"""
ls(dfg::CloudGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing)::Vector{Symbol} = getVariableIds(dfg, regexFilter)

"""
    $(SIGNATURES)
List the DFGFactors in the DFG.
Optionally specify a label regular expression to retrieves a subset of the factors.
"""
function getFactors(dfg::CloudGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing)::Vector{DFGFactor}
    factorIds = getFactorIds(dfg, regexFilter)
    return map(vId->getFactor(dfg, vId), factorIds)
end

"""
    $(SIGNATURES)
Get a list of the IDs of the DFGFactors in the DFG.
Optionally specify a label regular expression to retrieves a subset of the factors.
"""
function getFactorIds(dfg::CloudGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing)::Vector{Symbol}
    # Optimized for DB call
    if regexFilter == nothing
        return _getLabelsFromCyphonQuery(dfg.neo4jInstance, "(node:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId):FACTOR)")
    else
        return _getLabelsFromCyphonQuery(dfg.neo4jInstance, "(node:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId):FACTOR) where node.label =~ '$(regexFilter.pattern)'")
    end
end

# Alias
"""
    $(SIGNATURES)
List the DFGFactors in the DFG.
Optionally specify a label regular expression to retrieves a subset of the factors.
"""
lsf(dfg::CloudGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing)::Vector{Symbol} = getFactorIds(dfg, regexFilter)

# Alias - getNeighbors
"""
    $(SIGNATURES)
Get neighbors around a given node. TODO: Refactor this
"""
function lsf(dfg::CloudGraphsDFG, label::Symbol)::Vector{Symbol}
  return getNeighbors(dfg, label)
end

"""
    $(SIGNATURES)
Checks if the graph is fully connected, returns true if so.
"""
function isFullyConnected(dfg::CloudGraphsDFG)::Bool
    # If the total number of nodes == total number of distinct connected nodes, then it is fully connected
    # Total nodes
    varIds = getVariableIds(dfg)
    factIds = getFactorIds(dfg)
    length(varIds) + length(factIds) == 0 && return false

    # Total connected nodes - thank you Neo4j for 0..* awesomeness!!
    query = """
        match (n:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId):$(varIds[1]))-[FACTORGRAPH*]-(node:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId))
        WHERE (n:VARIABLE OR n:FACTOR OR node:VARIABLE OR node:FACTOR) and not (node:SESSION)
        WITH collect(n)+collect(node) as nodelist
        unwind nodelist as nodes
        return count(distinct nodes)"""
    @debug "[Querying] $query"
    result = _queryNeo4j(dfg.neo4jInstance, query)
    # Neo4j.jl data structure sometimes feels brittle... like below
    return result.results[1]["data"][1]["row"][1] == length(varIds) + length(factIds)
end

#Alias
"""
    $(SIGNATURES)
Checks if the graph is not fully connected, returns true if it is not contiguous.
"""
hasOrphans(dfg::CloudGraphsDFG)::Bool = !isFullyConnected(dfg)

"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor.
"""
function getNeighbors(dfg::CloudGraphsDFG, node::T; ready::Union{Nothing, Int}=nothing, backendset::Union{Nothing, Int}=nothing)::Vector{Symbol}  where T <: DFGNode
    query = "(n:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId):$(node.label))--(node) where (node:VARIABLE or node:FACTOR) "
    if ready != nothing || backendset != nothing
        if ready != nothing
            query = query * " and node.ready = $(ready)"
        end
        if backendset != nothing
            query = query * " and node.backendset = $(backendset)"
        end
    end
    @debug "[Query] $query"
    neighbors = _getLabelsFromCyphonQuery(dfg.neo4jInstance, query)
    # If factor, need to do variable ordering
    if T <: DFGFactor
        neighbors = intersect(node._variableOrderSymbols, neighbors)
    end
    return neighbors
end

"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor specified by its label.
"""
function getNeighbors(dfg::CloudGraphsDFG, label::Symbol; ready::Union{Nothing, Int}=nothing, backendset::Union{Nothing, Int}=nothing)::Vector{Symbol}
    query = "(n:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId):$(label))--(node) where (node:VARIABLE or node:FACTOR) "
    if ready != nothing || backendset != nothing
        if ready != nothing
            query = query * " and node.ready = $(ready)"
        end
        if backendset != nothing
            query = query * " and node.backendset = $(backendset)"
        end
    end
    @debug "[Query] $query"
    neighbors = _getLabelsFromCyphonQuery(dfg.neo4jInstance, query)
    # If factor, need to do variable ordering
    if isFactor(dfg, label)
        # Server is authority
        serverOrder = Symbol.(JSON2.read(_getNodeProperty(dfg.neo4jInstance, [dfg.userId, dfg.robotId, dfg.sessionId, String(label)], "_variableOrderSymbols")))
        neighbors = intersect(serverOrder, neighbors)
    end
    return neighbors
end

# Aliases
"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor.
"""
function ls(dfg::CloudGraphsDFG, node::T)::Vector{Symbol} where T <: DFGNode
    return getNeighbors(dfg, node)
end
"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor specified by its label.
"""
function ls(dfg::CloudGraphsDFG, label::Symbol)::Vector{Symbol} where T <: DFGNode
    return getNeighbors(dfg, label)
end

## This is moved to services/AbstractDFG.jl
# function _copyIntoGraph!(sourceDFG::CloudGraphsDFG, destDFG::CloudGraphsDFG, variableFactorLabels::Vector{Symbol}, includeOrphanFactors::Bool=false)::Nothing

"""
    $(SIGNATURES)
Retrieve a deep subgraph copy around a given variable or factor.
Optionally provide a distance to specify the number of edges should be followed.
Optionally provide an existing subgraph addToDFG, the extracted nodes will be copied into this graph. By default a new subgraph will be created.
Note: By default orphaned factors (where the subgraph does not contain all the related variables) are not returned. Set includeOrphanFactors to return the orphans irrespective of whether the subgraph contains all the variables.
"""
function getSubgraphAroundNode(dfg::CloudGraphsDFG, node::DFGNode, distance::Int64=1, includeOrphanFactors::Bool=false, addToDFG::AbstractDFG=_getDuplicatedEmptyDFG(dfg))::AbstractDFG
    distance < 1 && error("getSubgraphAroundNode() only works for distance > 0")

    # Thank you Neo4j for 0..* awesomeness!!
    neighborList = _getLabelsFromCyphonQuery(dfg.neo4jInstance,
        "(n:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId):$(node.label))-[FACTORGRAPH*0..$distance]-(node:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId)) WHERE (n:VARIABLE OR n:FACTOR OR node:VARIABLE OR node:FACTOR) and not (node:SESSION)")

    # Copy the section of graph we want
    _copyIntoGraph!(dfg, addToDFG, neighborList, includeOrphanFactors)
    return addToDFG
end

"""
    $(SIGNATURES)
Get a deep subgraph copy from the DFG given a list of variables and factors.
Optionally provide an existing subgraph addToDFG, the extracted nodes will be copied into this graph. By default a new subgraph will be created.
Note: By default orphaned factors (where the subgraph does not contain all the related variables) are not returned. Set includeOrphanFactors to return the orphans irrespective of whether the subgraph contains all the variables.
"""
function getSubgraph(dfg::CloudGraphsDFG,
                     variableFactorLabels::Vector{Symbol},
                     includeOrphanFactors::Bool=false,
                     addToDFG::G=_getDuplicatedEmptyDFG(dfg) )::G where {G <: AbstractDFG}
    # Making a copy session if not specified

    _copyIntoGraph!(dfg, addToDFG, variableFactorLabels, includeOrphanFactors)

    return addToDFG
end

"""
    $(SIGNATURES)
Get an adjacency matrix for the DFG, returned as a Matrix{Union{Nothing, Symbol}}.
Rows are all factors, columns are all variables, and each cell contains either nothing or the symbol of the relating factor.
The first row and first column are factor and variable headings respectively.
This is optimized for database usage.
"""
function getAdjacencyMatrix(dfg::CloudGraphsDFG)::Matrix{Union{Nothing, Symbol}}
    varLabels = sort(getVariableIds(dfg))
    factLabels = sort(getFactorIds(dfg))
    vDict = Dict(varLabels .=> [1:length(varLabels)...].+1)
    fDict = Dict(factLabels .=> [1:length(factLabels)...].+1)

    adjMat = Matrix{Union{Nothing, Symbol}}(nothing, length(factLabels)+1, length(varLabels)+1)
    # Set row/col headings
    adjMat[2:end, 1] = factLabels
    adjMat[1, 2:end] = varLabels

    # Now ask for all relationships for this session graph
    loadtx = transaction(dfg.neo4jInstance.connection)
    query = "START n=node(*) MATCH (n:VARIABLE:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId))-[r:FACTORGRAPH]-(m:FACTOR:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId)) RETURN n.label as variable, m.label as factor;"
    nodes = loadtx(query; submit=true)
    # Have to finish the transaction
    commit(loadtx)
    if length(nodes.errors) > 0
        error(string(nodes.errors))
    end
    # Add in the relationships
    varRel = Symbol.(map(node -> node["row"][1], nodes.results[1]["data"]))
    factRel = Symbol.(map(node -> node["row"][2], nodes.results[1]["data"]))
    for i = 1:length(varRel)
        adjMat[fDict[factRel[i]], vDict[varRel[i]]] = factRel[i]
    end

    return adjMat
end

function getAdjacencyMatrixSparse(dfg::CloudGraphsDFG)::Tuple{SparseMatrixCSC, Vector{Symbol}, Vector{Symbol}}
    varLabels = getVariableIds(dfg)
    factLabels = getFactorIds(dfg)
    vDict = Dict(varLabels .=> [1:length(varLabels)...])
    fDict = Dict(factLabels .=> [1:length(factLabels)...])

    adjMat = spzeros(Int, length(factLabels), length(varLabels))

    # Now ask for all relationships for this session graph
    loadtx = transaction(dfg.neo4jInstance.connection)
    query = "START n=node(*) MATCH (n:VARIABLE:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId))-[r:FACTORGRAPH]-(m:FACTOR:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId)) RETURN n.label as variable, m.label as factor;"
    nodes = loadtx(query; submit=true)
    # Have to finish the transaction
    commit(loadtx)
    if length(nodes.errors) > 0
        error(string(nodes.errors))
    end
    # Add in the relationships
    varRel = Symbol.(map(node -> node["row"][1], nodes.results[1]["data"]))
    factRel = Symbol.(map(node -> node["row"][2], nodes.results[1]["data"]))
    for i = 1:length(varRel)
        adjMat[fDict[factRel[i]], vDict[varRel[i]]] = 1
    end

    return adjMat, varLabels, factLabels
end

# """
#     $(SIGNATURES)
# Produces a dot-format of the graph for visualization.
# """
# function toDot(dfg::CloudGraphsDFG)::String
#     m = PipeBuffer()
#     write(m,Graphs.to_dot(dfg.g))
#     data = take!(m)
#     close(m)
#     return String(data)
# end
#
# """
#     $(SIGNATURES)
# Produces a dot file of the graph for visualization.
# Download XDot to see the data
#
# Note
# - Default location "/tmp/dfg.dot" -- MIGHT BE REMOVED
# - Can be viewed with the `xdot` system application.
# - Based on graphviz.org
# """
# function toDotFile(dfg::CloudGraphsDFG, fileName::String="/tmp/dfg.dot")::Nothing
#     open(fileName, "w") do fid
#         write(fid,Graphs.to_dot(dfg.g))
#     end
#     return nothing
# end


text:
@doc """
 COST 231- Cost-Hata Extension Model
 freq: signal frequency
 txH: base station's height
 rxH: mobile's height
 cityKind : cyte type(1-small, 2-medium e 3-large).
"""
type Cost231HataModel<:Model
  freq::Float64
  txH::Float64
  rxH::Float64
  areaKind::Int
  checkFreqRange::Bool
  Cost231HataModel() = new(1650,50.0,1.5,CityKind.Medium,true)
end

"Calculate Path loss for distance dist (km) with COST 231 Hata Extension Model"
function pathloss_op(model::Cost231HataModel, dist::AbstractFloat)
  #150 a 2000 MHz
  hm::Float64 = model.rxH
  hb::Float64 = model.txH
  d::Float64 = dist
  f::Float64 = model.freq
  C ::Float64 = model.areaKind==AreaKind.Urban ? 3 : 0
  ar::Float64 =(1.1*log10(f)-0.7)*hm-(1.56*log(f)-0.8)

  # Checking the frequency range
  if (model.checkFreqRange)
    if (f < 1500 || f > 2000)
     error("The frequency range for COST213 Hata Extension Model is 1500MHz-2000Mhz")
    end
  end
  #Pathloss
  L::Float64 = 46.3 +33.9*log10(f)-13.82*log10(hb)-ar+(44.9-6.55*log(hb))*log(d)+C
  return L

end


text:
"""
    EncodingDirection

Possible values: `ipos`, `jpos`, `kpos`, `ineg, `jneg`, `kneg` (the axis of the NIfTI data along which
slices were acquired, and the direction in which SliceTiming is defined with
respect to). `ipos`, `jpos`, `kpos` identifiers correspond to the first, second and third axis
of the data in the NIfTI file. `*neg` indicates that the contents of
SliceTiming are defined in reverse order - that is, the first entry corresponds
to the slice with the largest index, and the final entry corresponds to slice
index zero.
"""
@enum EncodingDirection begin
    ipos = 1
    ineg = -1
    jpos = 2
    jneg = -2
    kpos = 3
    kneg = -3
end

EncodingDirection(e::AbstractString) = EncodingDirection(Symbol(e))

function EncodingDirection(e::Symbol)
    if  (e === Symbol("i")) | (e === Symbol("ipos"))
        return ipos
    elseif (e === Symbol("i-")) | (e === Symbol("ineg"))
        return ineg
    elseif  (e === Symbol("j")) | (e === Symbol("jpos"))
        return jpos
    elseif (e === Symbol("j-")) | (e === Symbol("jneg"))
        return jneg
    elseif  (e === Symbol("k")) | (e === Symbol("kpos"))
        return kpos
    elseif (e === Symbol("k-")) | (e === Symbol("kneg"))
        return kneg
    else
        error("$e is not a supported encoding direction.")
    end
end
Base.String(e::EncodingDirection) = String(Symbol(e))

@inline function encoding_name(x::Int)
    if x === 1
        return :left
    elseif x === -1
        return :right
    elseif x === 2
        return :posterior
    elseif x === -2
        return :anterior
    elseif x === 3
        return :inferior
    elseif x === -3
        return :superior
    else
        error("$x does not map to a dimension name.")
    end
end

encoding_name(x) = encoding_name(Int(x))


text:
# push!(LOAD_PATH,"../src/")
using Documenter
include("../src/simulate.jl")

makedocs(sitename="ISRL-MBD Documentation",
         doctest = false,
         format = :html,
         pages = ["Introduction" => "index.md",
                  "User Documentation" => [
                  "Simulate" => "simulate.md"
                  "Rigid Body" => "rigidbody.md"
                  "Joint" => "joint.md"
                  "Constraint Forces" => "constraintUK.md"
                  ]])

# using Documenter, ForwardDiff
#
# makedocs(modules=[ForwardDiff],
#          doctest = false,
#          format = :html,
#          sitename = "ForwardDiff",
#          pages = ["Introduction" => "index.md",
#                   "User Documentation" => [
#                     "Limitations of ForwardDiff" => "user/limitations.md",
#                     "Differentiation API" => "user/api.md",
#                     "Advanced Usage Guide" => "user/advanced.md",
#                     "Upgrading from Older Versions" => "user/upgrade.md"],
#                   "Developer Documentation" => [
#                     "How ForwardDiff Works" => "dev/how_it_works.md",
#                     "How to Contribute" => "dev/contributing.md"]])
#
# deploydocs(repo = "github.com/JuliaDiff/ForwardDiff.jl.git",
#            osname = "linux",
#            julia = "0.7",
#            target = "build",
#            deps = nothing,
#            make = nothing)

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
#=
NOTES:
- Stepsize will be updated in each iteration in phases Warmup, Adaptionˢˡᵒʷ, Adaptionᶠᵃˢᵗ.
- Within phases updates are the un-smoothed stepsizes.
- At the end of each iteration, the smoothed stepsize is returned. Hence, in Exploration
    the smoothed stepsize from phase Adaptionᶠᵃˢᵗ is used and kept constant.
- After the first proposal step in new iteration in phases Warmup, Adaptionˢˡᵒʷ, Adaptionᶠᵃˢᵗ,
    DualAverageParameter is updated and stepsize multiplier will only come into effect at
    NEXT iteration with un-smoothed epsilon, in case adaption is still performed.
=#

############################################################################################
"""
$(TYPEDEF)

Default configuration for stepsize adaption.

# Fields
$(TYPEDFIELDS)
"""
struct ConfigStepsize{A<:BaytesCore.UpdateBool,B<:BaytesCore.UpdateBool}
    "Initial Discretization size"
    ϵ::Float64
    "Step size adaption"
    stepsizeadaption::A
    "Boolean if initial stepsize should be estimated"
    initialstepsize::B
    function ConfigStepsize(;
        ϵ = 0.1,
        stepsizeadaption = BaytesCore.UpdateTrue(),
        initialstepsize = BaytesCore.UpdateFalse()
        )
        @argcheck ϵ > 0.0 "Discretization size has to be positive"
        return new{typeof(stepsizeadaption), typeof(initialstepsize)}(
            ϵ, stepsizeadaption, initialstepsize
            )
    end
end

############################################################################################
"""
$(TYPEDEF)

Contains information for default Dual Averaging algorithm.

# Fields
$(TYPEDFIELDS)
"""
struct DualAverageParameter{T<:AbstractFloat}
    "Target acceptance rate"
    δ::T
    "Regularization scale"
    γ::T
    "Relaxation exponent - for Average log step size"
    κ::T
    "Offset"
    t₀::Int64
    function DualAverageParameter(δ::T, γ::T, κ::T, t₀::Int64) where {T<:AbstractFloat}
        @argcheck 0 < δ < 1
        @argcheck γ > 0
        @argcheck 0.5 < κ ≤ 1
        @argcheck t₀ ≥ 0
        return new{T}(δ, γ, κ, t₀)
    end
end
function DualAverageParameter(
    ::Type{T}; δ=0.80, γ=0.05, κ=0.75, t₀=10
) where {T<:AbstractFloat}
    return DualAverageParameter(T(δ), T(γ), T(κ), t₀)
end

############################################################################################
"""
$(TYPEDEF)

Contains DualAverage tuning information and runtime parameter.

# Fields
$(TYPEDFIELDS)
"""
struct DualAverage{T<:AbstractFloat}
    #!NOTE: Updated each iteration
    adaption::DualAverageParameter{T}
    "Upwards bias for target acceptance rate - proposals are biased upwards to stay away from 0."
    μ::T
    "Time update, starts with 0"
    t::Int64
    "Average part of first equation in Hoffman(2014), p 1607, (6)"
    H̄::T
    "Log step"
    logϵ::T
    "AVERAGED log step"
    logϵ̄::T
    function DualAverage(
        adaption::DualAverageParameter{T}, μ::A, t::I, H̄::B, logϵ::C, logϵ̄::D
    ) where {T<:Real,I<:Integer, A<:Real, B<:Real, C<:Real, D<:Real}
        return new{T}(adaption, μ, t, H̄, logϵ, logϵ̄)
    end
end

############################################################################################
"""
$(SIGNATURES)
Initialize new DualAverage struct.

# Examples
```julia
```

"""
function init(::Type{DualAverage}, δ, ϵ::T; multiplier=T(10.0)) where {T<:AbstractFloat}
    @argcheck ϵ > 0
    logϵ = log(ϵ)
    return DualAverage(
        DualAverageParameter(T; δ=δ),
        log(multiplier) + logϵ,
        0,
        zero(logϵ),
        logϵ,
        zero(logϵ),
    )
end

############################################################################################
"""
$(SIGNATURES)
Update DualAverage struct.

# Examples
```julia
```

"""
function update(dualaverage::DualAverage, acceptratio::T) where {T<:AbstractFloat}
    @argcheck 0 ≤ acceptratio ≤ 1
    @unpack adaption, μ, t, H̄, logϵ, logϵ̄ = dualaverage
    @unpack δ, γ, κ, t₀ = adaption
    ## Update t and H
    t += 1
    H̄ += (δ - acceptratio - H̄) / (t + t₀)
    ## Calculate logϵ and MA
    logϵ = μ - T(√t) / γ * H̄
    logϵ̄ += t^(-κ) * (logϵ - logϵ̄)
    ## Pack container
    return DualAverage(adaption, μ, t, H̄, logϵ, logϵ̄) #T(logϵ), T(logϵ̄ ))
end

############################################################################################
"""
$(TYPEDEF)

Initial Step Size Search factors.

# Fields
$(TYPEDFIELDS)
"""
struct InitialStepsizeSearch{T<:AbstractFloat}
    "Lowest local acceptance rate."
    a_min::T
    "Highest local acceptance rate."
    a_max::T
    "Initial stepsize."
    ϵ₀::T
    "Scale factor for initial bracketing, > 1."
    C::T
    "Maximum number of iterations for initial bracketing."
    maxiter_crossing::Int64
    "Maximum number of iterations for bisection."
    maxiter_bisect::Int64
    function InitialStepsizeSearch(
        ::Type{T};
        a_min=0.25,
        a_max=0.75,
        ϵ₀=0.5,
        C=2.0,
        maxiter_crossing=1000,
        maxiter_bisect=1000,
    ) where {T}
        @argcheck 0 < a_min < a_max < 1
        @argcheck 0 < ϵ₀
        @argcheck 1 < C
        @argcheck maxiter_crossing ≥ 50
        @argcheck maxiter_bisect ≥ 50
        return new{T}(a_min, a_max, ϵ₀, C, maxiter_crossing, maxiter_bisect)
    end
end

############################################################################################
"""
$(SIGNATURES)
Find the stepsize for which the local acceptance rate `A(ϵ)` crosses `a`.

# Examples
```julia
```

"""
function find_crossing_stepsize(parameters, A, ϵ₀, Aϵ₀=A(ϵ₀))
    @unpack a_min, a_max, C, maxiter_crossing = parameters
    s, a = Aϵ₀ > a_max ? (1.0, a_max) : (-1.0, a_min)
    if s < 0                    # when A(ϵ) < a,
        C = 1 / C                 # decrease ϵ
    end
    for _ in 1:maxiter_crossing
        ϵ = ϵ₀ * C
        Aϵ = A(ϵ)
        if s * (Aϵ - a) ≤ 0
            return ϵ₀, Aϵ₀, ϵ, Aϵ
        else
            ϵ₀ = ϵ
            Aϵ₀ = Aϵ
        end
    end
    dir = s > 0 ? "below" : "above"
    return error("Reached maximum number of iterations searching for ϵ from $(dir).")
end

############################################################################################
"""
$(SIGNATURES)
Return the desired stepsize `ϵ` by bisection.

# Examples
```julia
```

"""
function bisect_stepsize(parameters, A, ϵ₀, ϵ₁, Aϵ₀=A(ϵ₀), Aϵ₁=A(ϵ₁))
    @unpack a_min, a_max, maxiter_bisect = parameters
    @argcheck ϵ₀ < ϵ₁
    @argcheck Aϵ₀ > a_max && Aϵ₁ < a_min
    for _ in 1:maxiter_bisect
        ϵₘ = Statistics.middle(ϵ₀, ϵ₁) #0.5*(ϵ₀ + ϵ₁)
        Aϵₘ = A(ϵₘ)
        if a_min ≤ Aϵₘ ≤ a_max  # in
            return ϵₘ
        elseif Aϵₘ < a_min      # above
            ϵ₁ = ϵₘ
            Aϵ₁ = Aϵₘ
        else                    # below
            ϵ₀ = ϵₘ
            Aϵ₀ = Aϵₘ
        end
    end
    return error("Reached maximum number of iterations while bisecting interval for ϵ.")
end

############################################################################################
"""
$(SIGNATURES)
Find an initial stepsize that matches the conditions of `parameters` (see
[`InitialStepsizeSearch`](@ref)). `A` is the local acceptance ratio (unbounded).

# Examples
```julia
```

"""
function find_initial_stepsize(parameters::InitialStepsizeSearch, A)
    @unpack a_min, a_max, ϵ₀ = parameters
    #!NOTE: A is a closure for the MCMC step to calculate acceptance rate for 1 MCMC step
    Aϵ₀ = A(ϵ₀)
    if a_min ≤ Aϵ₀ ≤ a_max
        ϵ₀
    else
        ϵ₀, Aϵ₀, ϵ₁, Aϵ₁ = find_crossing_stepsize(parameters, A, ϵ₀, Aϵ₀)
        #!NOTE: in interval
        if a_min ≤ Aϵ₁ ≤ a_max
            ϵ₁
            #!NOTE: order as necessary
        elseif ϵ₀ < ϵ₁
            bisect_stepsize(parameters, A, ϵ₀, ϵ₁, Aϵ₀, Aϵ₁)
        else
            bisect_stepsize(parameters, A, ϵ₁, ϵ₀, Aϵ₁, Aϵ₀)
        end
    end
end

############################################################################################
"""
$(TYPEDEF)

Tuning stepsize parameter for MCMC algorithm.

# Fields
$(TYPEDFIELDS)
"""
mutable struct StepSizeTune{A<:BaytesCore.UpdateBool,F<:AbstractFloat}
    "If true, stepsize will be adapted."
    adaption::A
    #!NOTE: isbits(dualaverage) == true
    "Dualaverage struct"
    dualaverage::DualAverage{F}
    "Current stepsize"
    ϵ::F
    function StepSizeTune(
        adaption::A, dualaverage::DualAverage{F}, ϵ::F
    ) where {A<:BaytesCore.UpdateBool,F<:AbstractFloat}
        @argcheck 0.0 < ϵ
        return new{A,F}(adaption, dualaverage, ϵ)
    end
end

"""
$(SIGNATURES)
Update stepsize based on acceptance ratio α.

# Examples
```julia
```

"""
function update!(
    stepsize::StepSizeTune,
    stepsizeadaption::BaytesCore.UpdateTrue,
    α::T,
    samplingphase::N,
    iterationupdate::Val{false},
) where {N<:Union{Warmup,Adaptionˢˡᵒʷ,Adaptionᶠᵃˢᵗ},T<:AbstractFloat}
    @unpack dualaverage = stepsize
    ## Update stepsize
    dualaverage = update(dualaverage, α)
    ϵ = exp(dualaverage.logϵ)
    @pack! stepsize = dualaverage, ϵ
    return nothing
end
function update!(
    stepsize::StepSizeTune,
    stepsizeadaption::BaytesCore.UpdateTrue,
    α::T,
    samplingphase::Exploration,
    iterationupdate::Val{false},
) where {T<:AbstractFloat}
    return nothing
end
function update!(
    stepsize::StepSizeTune,
    stepsizeadaption::BaytesCore.UpdateTrue,
    α::T,
    samplingphase::S,
    iterationupdate::Val{true},
) where {S<:SamplingPhase,T<:AbstractFloat}
    @unpack dualaverage = stepsize
    ## Update stepsize
    dualaverage = update(dualaverage, α)
    ϵ = exp(dualaverage.logϵ̄) #!NOTE: SMOOTHED VERSION
    ## Assign new dualaverage parameter with smoothed log epsilon
    #= This will be Dualaveraging Parameter for first update! in new samplingphase, hence in final sampling phase only the smoothed ϵ is returned,
    and dualaverage has no impact cause it wont compute anymore when Exploration+update=false. Hence, the stepsize multiplier will only come
    into effect at NEXT iteration with un-smoothed epsilon, in case adaption is still performed. =#
    dualaverage = init(DualAverage, dualaverage.adaption.δ, ϵ)
    @pack! stepsize = dualaverage, ϵ
    return nothing
end
function update!(
    stepsize::StepSizeTune, stepsizeadaption::BaytesCore.UpdateFalse, α, samplingphase, iterationupdate
)
    return nothing
end
function update!(stepsize::StepSizeTune, α::T, phase::PhaseTune) where {T<:AbstractFloat}
    return update!(
        stepsize,
        stepsize.adaption,
        α,
        phase.name[phase.counter.current],
        Val(phase.update.current),
    )
end

############################################################################################
# Export
export StepSizeTune, ConfigStepsize, DualAverage, DualAverageParameter


text:
### log pdf of a multivariate normal
function logpfun_mv(x, μ, σ, w)
    logw = log.(w)
    K = length(logw)
    N = size(x, 2)
    sum(logsumexp(logw[k] + logpdf(MvNormal(μ[k], σ), x[:,n]) for k in 1:K) for n in 1:N)
end

@model function (x, nc, nsites, σ,  p,  wc=ones(nc)/nc, wp=ones(nsites-1)/(nsites-1), w=wc.*wp', q = 2)

    n = size(x, 2)
    startpos ~ filldist(Uniform(-1500, 1500), nc)
    delta    ~ filldist(Gamma(q*p, 1/q), nc, nsites)
    configs = startpos .+ cumsum(delta, dims=2)

    μ = [[configs[i, j], configs[i, j+1]] for i = 1:nc, j=1:(nsites-1)]

    #x ~ filldist(MixtureModel(Normal.(vec(μ), w)), n)
    Turing.@addlogprob! logpfun_mv(x, vec(μ), σ, w)
end


text:
# # Online Gaussian Process
#
# ### Loading necessary packages
# ## Preliminary steps
# ### Load the necessary packages
using Plots
using AugmentedGaussianProcesses
using MLDataUtils, Distributions

# ### We create a toy dataset with a noisy sinus
N = 2000
σ = 0.1
X, y = noisy_sin(N, 0, 20; noise=σ)
X_train = X[1:2:end];
y_train = y[1:2:end]; # We split the data equally
X_test = X[2:2:end];
y_test = y[2:2:end];
scatter(X_train, y_train)

# ### Plot model at each step
function plot_model(model, X, X_test, X_train, y_train)
    y_pred, sig_pred = proba_y(model, X_test)
    plot(X, sin; lab="f", color=:black, lw=3.0, ylims=(-2, 2))
    plot!(X_test, y_pred; ribbon=sqrt.(sig_pred), lab="Prediction", lw=3.0)
    scatter!(X_train, y_train; msw=0.0, alpha=0.5, lab="Data")
    return scatter!(first.(model[1].Z), mean(model[1]); lab="IP")
end

# ## Model training
# ### Create a kernel
k = SqExponentialKernel();
# ### Create an inducing point selection method
IP_alg = OIPS(0.8);
# ### Create the model and stream the data
model = OnlineSVGP(k, GaussianLikelihood(σ), AnalyticVI(), IP_alg; optimiser=false)
anim = Animation()
size_batch = 100
for (i, (X_batch, y_batch)) in enumerate(eachbatch((X_train, y_train); obsdim=1, size=size_batch))
    train!(model, X_batch, y_batch; iterations=3)
    plot_model(model, X, X_test, X_train[1:(i * size_batch)], y_train[1:(i * size_batch)])
    frame(anim)
end
gif(anim; fps=4)
# This works just as well with any likelihood! Just try it out!