text:
default_diff_fn(c1::CT1,c2::CT2) where {CT1<:Union{Colorant,Real}, CT2<:Union{Colorant,Real}} = sqrt(sum(abs2,(c1)-accum_type(CT2)(c2)))

"""
    seg_img = seeded_region_growing(img, seeds, [kernel_dim], [diff_fn])
    seg_img = seeded_region_growing(img, seeds, [neighbourhood], [diff_fn])

Segments the N-D image `img` using the seeded region growing algorithm
and returns a [`SegmentedImage`](@ref) containing information about the segments.

# Arguments:
* `img`             :  N-D image to be segmented (arbitrary axes are allowed)
* `seeds`           :  `Vector` containing seeds. Each seed is a Tuple of a
                       CartesianIndex{N} and a label. See below note for more
                       information on labels.
* `kernel_dim`      :  (Optional) `Vector{Int}` having length N or a `NTuple{N,Int}`
                       whose ith element is an odd positive integer representing
                       the length of the ith edge of the N-orthotopic neighbourhood
* `neighbourhood`   :  (Optional) Function taking CartesianIndex{N} as input and
                       returning the neighbourhood of that point.
* `diff_fn`         :  (Optional) Function that returns a difference measure(δ)
                       between the mean color of a region and color of a point

!!! note
    The labels attached to points must be positive integers, although multiple
    points can be assigned the same label. The output includes a labelled array
    that has same indexing as that of input image. Every index is assigned to
    either one of labels or a special label '0' indicating that the algorithm
    was unable to assign that index to a unique label.

# Examples

```jldoctest; setup = :(using ImageCore, ImageMorphology,ImageSegmentation)
julia> img = zeros(Gray{N0f8},4,4);

julia> img[2:4,2:4] .= 1;

julia> seeds = [(CartesianIndex(3,1),1),(CartesianIndex(2,2),2)];

julia> seg = seeded_region_growing(img, seeds);

julia> labels_map(seg)
4×4 $(Matrix{Int}):
 1  1  1  1
 1  2  2  2
 1  2  2  2
 1  2  2  2
```

# Citation:

Albert Mehnert, Paul Jackaway (1997), "An improved seeded region growing algorithm",
Pattern Recognition Letters 18 (1997), 1065-1071
"""
function seeded_region_growing(img::AbstractArray{CT,N}, seeds::AbstractVector{Tuple{CartesianIndex{N},Int}},
    kernel_dim::Union{Vector{Int}, NTuple{N, Int}} = ntuple(i->3,N), diff_fn::Function = default_diff_fn) where {CT<:Union{Colorant,Real}, N}
    length(kernel_dim) == N || error("Dimension count of image and kernel_dim do not match")
    for dim in kernel_dim
        dim > 0 || error("Dimensions of the kernel must be positive")
        isodd(dim) || error("Dimensions of the kernel must be odd")
    end
    pt = CartesianIndex(ntuple(i->kernel_dim[i]÷2, N))
    neighbourhood_gen(t) = c->CartesianIndices(_colon(c-t,c+t))
    seeded_region_growing(img, seeds, neighbourhood_gen(pt), diff_fn)
end

function seeded_region_growing(img::AbstractArray{CT,N}, seeds::AbstractVector{Tuple{CartesianIndex{N},Int}},
    neighbourhood::Function, diff_fn::Function = default_diff_fn) where {CT<:Union{Colorant,Real}, N}

    _QUEUE_SZ = 10

    # Check if labels are positive integers
    for seed in seeds
        seed[2] > 0 || error("Seed labels need to be positive integers!")
    end

    TM = meantype(CT)

    # Required data structures
    result              =   fill(-1, axes(img))                                     # Array to store labels
    nhq                 =   Queue{CartesianIndex{N}}(_QUEUE_SZ)                     # Neighbours holding queue
    pq                  =   PriorityQueue{CartesianIndex{N}, Float64}()             # Priority Queue to hold the pixels of same δ value
    labelsq             =   Queue{Int}(_QUEUE_SZ)                                   # Queue to hold labels
    holdingq            =   Queue{CartesianIndex{N}}(_QUEUE_SZ)                     # Queue to hold points corresponding to the labels in `labelsq`
    region_means        =   Dict{Int, TM}()                                         # A map containing (label, mean) pairs
    region_pix_count    =   Dict{Int, Int}()                                        # A map containing (label, pixel_count) pairs
    labels              =   Vector{Int}()                                           # A vector containing list of labels

    # Labelling initial seeds and initialising `region_means` and `region_pix_count`
    for seed in seeds
        result[seed[1]] = seed[2]
        region_pix_count[seed[2]] = get(region_pix_count, seed[2], 0) + 1
        region_means[seed[2]] = get(region_means, seed[2], zero(TM)) + (img[seed[1]] - get(region_means, seed[2], zero(TM)))/(region_pix_count[seed[2]])
        if ! (seed[2] in labels)
            push!(labels, seed[2])
        end
    end

    #=  Labeling scheme for the Array "result"-
            Unlabelled => -1
            In nhq => -2
            In pq => -3
            Tied => 0
            Labelled => Seed value
    =#

    # Enqueue all the neighbours of seeds into `nhq`
    for seed in seeds
        for point in neighbourhood(seed[1])
            if point != seed[1] && checkbounds(Bool, img, point) && result[point] == -1
                enqueue!(nhq, point)
                @inbounds result[point] = -2
            end
        end
    end

    while !isempty(pq) || !isempty(nhq)

        while !isempty(nhq)
            # For every neighbouring point, get the minimum δ
            p = dequeue!(nhq)
            δ = Inf
            for point in neighbourhood(p)
                if point != p && checkbounds(Bool, img, point) && result[point] > 0
                    curr_diff = diff_fn(region_means[result[point]], img[p])
                    if δ > curr_diff
                        δ = curr_diff
                    end
                end
            end
            pq[p] = δ
            @inbounds result[p] = -3
        end

        # Get the pixels with minimum δ from `pq` and add them to `holdingq` and their labels to `labelsq`
        if !isempty(pq)
            δ_min = peek(pq)[2]
        end
        while (!isempty(pq) && isapprox(peek(pq)[2], δ_min)) #, atol=1e-8))
            p = dequeue!(pq)
            mindifflabel = -1
            mindiff = Inf
            for point in neighbourhood(p)
                if point!=p && checkbounds(Bool, img, point) && result[point] > 0
                    if mindifflabel < 0
                        mindifflabel = result[point]
                        mindiff = diff_fn(region_means[result[point]], img[p])
                    elseif mindifflabel != result[point]
                        curr_diff = diff_fn(region_means[result[point]], img[p])
                        if curr_diff < mindiff
                            mindiff = curr_diff
                            mindifflabel = result[point]
                        elseif isapprox(curr_diff, mindiff) #, atol=1e-8)
                            mindifflabel = 0
                        end
                    end
                end
            end
            if mindifflabel != 0                # new labels, including resolved ties in last step
                enqueue!(labelsq, mindifflabel)
                enqueue!(holdingq, p)
                @inbounds result[p] = -2
            elseif result[p] == -3              # new ties are enqueued in pq
                pq[p] = Inf
                result[p] = 0
            end                                 # old ties are not enqueued again
        end

        # Add label to each point in `holdingq` and add their neighbours to `nhq`
        while !isempty(holdingq)
            label = dequeue!(labelsq)
            p = dequeue!(holdingq)
            result[p] = label
            region_pix_count[label] += 1
            region_means[label] += (img[p] - region_means[label])/(region_pix_count[label])
            for point in neighbourhood(p)
                if point!=p && checkbounds(Bool, img, point) && (result[point] == -1 || result[point] == -3)
                    enqueue!(nhq, point)
                    @inbounds result[point] = -2
                end
            end
        end

    end

    c0 = count(i->(i==0),result)
    if c0 != 0
        push!(labels, 0)
        region_pix_count[0] = c0
    end
    SegmentedImage(result, labels, region_means, region_pix_count)
end

"""
    seg_img = unseeded_region_growing(img, threshold, [kernel_dim], [diff_fn])
    seg_img = unseeded_region_growing(img, threshold, [neighbourhood], [diff_fn])

Segments the N-D image using automatic (unseeded) region growing algorithm
and returns a [`SegmentedImage`](@ref) containing information about the segments.

# Arguments:
* `img`             :  N-D image to be segmented (arbitrary axes are allowed)
* `threshold`       :  Upper bound of the difference measure (δ) for considering
                       pixel into same segment
* `kernel_dim`      :  (Optional) `Vector{Int}` having length N or a `NTuple{N,Int}`
                       whose ith element is an odd positive integer representing
                       the length of the ith edge of the N-orthotopic neighbourhood
* `neighbourhood`   :  (Optional) Function taking CartesianIndex{N} as input and
                       returning the neighbourhood of that point.
* `diff_fn`         :  (Optional) Function that returns a difference measure (δ)
                       between the mean color of a region and color of a point

# Examples

```jldoctest; setup = :(using ImageCore, ImageMorphology, ImageSegmentation)
julia> img = zeros(Gray{N0f8},4,4);

julia> img[2:4,2:4] .= 1;

julia> seg = unseeded_region_growing(img, 0.2);

julia> labels_map(seg)
4×4 $(Matrix{Int}):
 1  1  1  1
 1  2  2  2
 1  2  2  2
 1  2  2  2
```

"""
function unseeded_region_growing(img::AbstractArray{CT,N}, threshold::Real,
    kernel_dim::Union{Vector{Int}, NTuple{N, Int}} = ntuple(i->3,N), diff_fn::Function = default_diff_fn) where {CT<:Colorant, N}
    length(kernel_dim) == N || error("Dimension count of image and kernel_dim do not match")
    for dim in kernel_dim
        dim > 0 || error("Dimensions of the kernel must be positive")
        isodd(dim) || error("Dimensions of the kernel must be odd")
    end
    pt = CartesianIndex(ntuple(i->kernel_dim[i]÷2, N))
    neighbourhood_gen(t) = c->CartesianIndices(_colon(c-t,c+t))
    unseeded_region_growing(img, threshold, neighbourhood_gen(pt), diff_fn)
end

function unseeded_region_growing(img::AbstractArray{CT,N}, threshold::Real, neighbourhood::Function, diff_fn = default_diff_fn) where {CT<:Colorant,N}
    TM = meantype(CT)

    # Required data structures
    result                  =   fill(-1, axes(img))                             # Array to store labels
    neighbours              =   PriorityQueue{CartesianIndex{N},Float64}()      # Priority Queue containing boundary pixels with δ as the priority
    region_means            =   Dict{Int, TM}()                                 # A map containing (label, mean) pairs
    region_pix_count        =   Dict{Int, Int}()                                # A map containing (label, pixel_count) pairs
    labels                  =   Vector{Int}()                                   # Vector containing assigned labels

    # Initialize data structures
    start_point = first(CartesianIndices(axes(img)))
    result[start_point] = 1
    push!(labels, 1)
    region_means[1] = img[start_point]
    region_pix_count[1] = 1

    # Enqueue neighouring points of `start_point`
    for p in neighbourhood(start_point)
        if p != start_point && checkbounds(Bool, img, p) && result[p] == -1
            enqueue!(neighbours, p, diff_fn(region_means[result[start_point]], img[p]))
        end
    end

    while !isempty(neighbours)
        point = dequeue!(neighbours)
        δ = Inf
        minlabel = -1
        pixelval = img[point]
        for p in neighbourhood(point)
            if p != point && checkbounds(Bool, img, p) && result[p] > 0
                curδ = diff_fn(region_means[result[p]], pixelval)
                if curδ < δ
                    δ = curδ
                    minlabel = result[p]
                end
            end
        end

        if δ < threshold
            # Assign point to minlabel
            result[point] = minlabel
        else
            # Select most appropriate label
            δ = Inf
            minlabel = -1
            for label in labels
                curδ = diff_fn(region_means[label], pixelval)
                if curδ < δ
                    δ = curδ
                    minlabel = label
                end
            end

            if δ < threshold
                result[point] = minlabel
            else
                # Assign point to a new label
                minlabel = length(labels) + 1
                push!(labels, minlabel)
                result[point] = minlabel
            end
        end

        #Update region_means
        region_pix_count[minlabel] = get(region_pix_count, minlabel, 0) + 1
        region_means[minlabel] = get(region_means, minlabel, zero(TM)) + (pixelval - get(region_means, minlabel, zero(TM)))/(region_pix_count[minlabel])

        # Enqueue neighbours of `point`
        for p in neighbourhood(point)
            if checkbounds(Bool, img, p) && result[p] == -1
                if haskey(neighbours, p)
                    continue
                end
                δ = Inf
                for tp in neighbourhood(p)
                    if checkbounds(Bool, img, tp) && tp != p && result[tp] > 0
                        δ = min(δ, diff_fn(region_means[result[tp]], img[p]))
                    end
                end
                enqueue!(neighbours, p, δ)
            end
        end

    end
    SegmentedImage(result, labels, region_means, region_pix_count)
end


text:
function maxent_max_pol(method, model_key)

    ## -------------------------------------------------------------------
    # Monitor
    mon = SimT.OnDiskMonitor(cachedir(iJR, :ME_MONITOR))
    SimT.reset!(mon)

    # Feed jobs
    Ch = Channel(nthreads()) do ch
        cGLCs = Nd.val("cGLC")
        for (exp, cGLC)  in enumerate(cGLCs)
            put!(ch, (exp, cGLC))
        end
    end

    @threads for _ in 1:nthreads()
        thid = threadid()
        for (exp, cGLC) in Ch

            ## -------------------------------------------------------------------
            # handle cache
            is_cached(;method, exp) && continue

            ## -------------------------------------------------------------------
            # Excluded because of unfeasibility or convergence problems
            exp in IGNORE_EXPS && continue

            ## -------------------------------------------------------------------
            # SetUp
            model = iJR.load_model(model_key)
            M, N = size(model)
            biomidx = ChU.rxnindex(model, iJR.BIOMASS_IDER)
            glcidx = ChU.rxnindex(model, iJR.EX_GLC_IDER)
            exp_growth = Nd.val("D", exp) # experimental biom growth rate (equals D)
            # glc per biomass unit supply
            cgD_X = -Nd.cval(:GLC, exp) * Nd.val(:D, exp) / Nd.val(:X, exp)
            biom_beta = 0.0 # current biomass beta
            biom_betas = [biom_beta] # biomass beta round history
            vg_beta = 0.0 # current vg beta
            vg_betas = [vg_beta] # vg beta round history
            vg_avPME = 0.0 # current vg exchange average
            vg_avPME_vgb0 = 0.0 # vg exchange average at vg_beta = 0
            biom_avPME = 0.0 # current biom exchange average
            biom_avPME_vgb0 = 0.0 # biom exchange average at vg_beta = 0
            biom_diff = 0.0 # distance between exp_growth and biom_avPME
            vg_diff = 0.0 # distance between cgD_X and vg_diff
            beta_vec = zeros(N) # ep beta vector
            epouts = Dict() # epout pool for each round (will contain the solution)
            epout = nothing # current epout (used as seed)
            hasvalid_moments = false # a flag that indicate is the momentous are valid
            roundconv = false # global (round) converge flag

            ## ----------------------------------------------------
            # ep params
            me_params = lglob(iJR, :maxent, :params)
            @extract me_params: alpha epsconv epmaxiter=maxiter
            @extract me_params: damp maxvar minvar

            ## ----------------------------------------------------
            # other params
            gdmaxiter = 3000 # maxiter for each gradient descent
            gdth = 0.01  # th of each gradient descend
            roundth = 0.01 # th of the whole simulation
            smooth = 0.1 # gd smooth th

            # After a while without converge, accelerate
            turbo_iter0 = 100 # iter for initing turbo
            turbo_frec = 10 # iter frec for apply turbo
            turbo_factor = 2.0 # turbo strength

            # damp will be reduced when detected
            damp_factor = 0.5 # gd damp penalty factor
            biom_gddamp = 1.0 # biom gd current damp
            vg_gddamp = 1.0 # vg gd current damp

            gdit = -1 # current gd iter
            gderr = -1.0 # current gd error
            last_infotime = -1.0 # time to check if gd needs to update
            upfrec_time = 15 # update info frequency

            # the whole simulation converge as ~log,
            # so I force betas increment by stepping
            beta_scale_rounditer0 = 3 # starting round for beta stepping
            beta_scale_factor = 0.5 # stepping scaling factor

            rounditer = 1 # current round iter
            maxrounds = 50 # max no of rounds

            # monitor
            SimT.record!(mon) do dat
                tdat = get!(dat, exp, Dict())
                tdat[:cgD_X] = cgD_X
                tdat[:method] = method
                tdat[:exp_growth] = exp_growth
            end

            ## -------------------------------------------------------------------
            # closure functions
            function check_roundconv()
                hasvalid_biom_moment = abs(biom_avPME - exp_growth)/abs(exp_growth) <= roundth
                hasvalid_vg_moment = abs(vg_avPME) <= abs(cgD_X) ||
                    abs(vg_avPME - cgD_X)/abs(cgD_X) <= roundth
                hasvalid_moments = hasvalid_biom_moment && hasvalid_vg_moment
                return hasvalid_moments
            end

            ## -------------------------------------------------------------------
            function print_info(msg; varargs...)

                isinfotime = gdit == 1 || abs(last_infotime - time()) > upfrec_time ||
                    epout.status != ChEP.CONVERGED_STATUS

                !isinfotime && return

                lock(WLOCK) do
                    @info(msg,
                        varargs...,
                        epout.status, epout.iter,
                        (biom_avPME_vgb0, biom_avPME, exp_growth), biom_diff,
                        (vg_avPME_vgb0, vg_avPME, cgD_X), vg_diff,
                        (biom_beta, vg_beta),
                        thid
                    ); println()
                    last_infotime = time()
                end
            end

            function check_nan()
                params = [biom_beta, vg_beta, biom_avPME, vg_avPME]
                any(isnan.(params))
            end

            ## -------------------------------------------------------------------
            function gd_core_fun(gdmodel; msg)

                beta_vec[biomidx] = biom_beta
                beta_vec[glcidx] = vg_beta

                # MAXENT
                epout = ChEP.maxent_ep(model;
                    beta_vec,
                    alpha, epsconv,
                    maxvar, minvar, damp,
                    maxiter = epmaxiter,
                    verbose = false,
                    solution = epout
                )

                # update
                biom_avPME = ChU.av(model, epout, iJR.BIOMASS_IDER)
                vg_avPME = ChU.av(model, epout, iJR.EX_GLC_IDER)
                biom_diff = abs(biom_avPME - exp_growth)
                vg_diff = abs(vg_avPME - cgD_X)

                gderr = gdmodel.ϵi
                print_info(msg;
                    exp, rounditer, gdit, gderr,
                    vg_gddamp, biom_gddamp,
                )

                # MONITOR
                SimT.record!(mon) do dat
                    tdat = get!(dat, exp, Dict())
                    tdat[:live_prove] = rand()
                    gddat = get!(tdat, :gd, Dict())
                    SimT.get!push!(gddat;
                        vg_beta, biom_beta,
                        biom_avPME, vg_avPME
                    )
                end

                # RUN OUT OF PATIENT
                (gdit >= turbo_iter0 && rem(gdit, turbo_frec) == 0) &&
                    (gdmodel.maxΔx *= turbo_factor)

                gdit += 1
            end

            while true

                ## -------------------------------------------------------------------
                # BETA SCALING
                scalebeta = rounditer >= beta_scale_rounditer0
                scalebeta && let
                    biom_beta_step = biom_betas[end] - biom_betas[end - 1]
                    biom_beta += biom_beta_step * beta_scale_factor

                    vg_beta_step = vg_betas[end] - vg_betas[end - 1]
                    vg_beta += vg_beta_step * beta_scale_factor
                end

                ## -------------------------------------------------------------------
                # Z GRAD DESCEND: Match biomass momentums
                let
                    target = exp_growth
                    x0 = biom_beta
                    maxΔx = max(abs(biom_beta) * 0.05, 5e2)
                    minΔx = maxΔx * 0.001
                    x1 = x0 + maxΔx * 0.01

                    gdit = 1

                    ## -------------------------------------------------------------------
                    function z_fun(gdmodel)
                        biom_beta = SimT.gd_value(gdmodel)
                        gd_core_fun(gdmodel; msg = "z grad descent... ")
                        return biom_avPME
                    end

                    ## -------------------------------------------------------------------
                    function z_break_cond(gdmodel)
                        roundconv = check_roundconv()
                        zconv = abs(biom_avPME - exp_growth)/abs(exp_growth) <= gdth
                        roundconv || zconv
                    end

                    ## -------------------------------------------------------------------
                    gdmodel = SimT.grad_desc(z_fun;
                        x0, x1, gdth, minΔx, maxΔx, smooth,
                        target, maxiter = gdmaxiter,
                        damp_factor, damp = biom_gddamp,
                        break_cond = z_break_cond,
                        verbose = false
                    )
                    biom_beta = SimT.gd_value(gdmodel)
                    biom_gddamp = gdmodel.damp
                end

                ## -------------------------------------------------------------------
                # AT VG BETA 0 MOMENTS
                isfirstround = (rounditer == 1)
                isfirstround && let
                    biom_avPME_vgb0 = biom_avPME
                    vg_avPME_vgb0 = vg_avPME
                end

                ## -------------------------------------------------------------------
                # FORCE VG BOUNDARY
                let
                    ## -------------------------------------------------------------------
                    # VG GRAD DESCEND: Match biomass momentums
                    target = cgD_X * 0.99 # force to be inside
                    x0 = vg_beta
                    maxΔx = max(abs(vg_beta) * 0.05, 1e2)
                    minΔx = maxΔx * 0.001
                    x1 = x0 + maxΔx * 0.01

                    gdit = 1

                    ## -------------------------------------------------------------------
                    function vg_fun(gdmodel)
                        vg_beta = SimT.gd_value(gdmodel)
                        gd_core_fun(gdmodel; msg = "vg grad descent... ")
                        return vg_avPME
                    end

                    ## -------------------------------------------------------------------
                    function vg_break_cond(epmodel)
                        vgconv = abs(vg_avPME) <= abs(cgD_X)
                        roundconv = check_roundconv()
                        roundconv || vgconv
                    end

                    ## -------------------------------------------------------------------
                    gdmodel = SimT.grad_desc(vg_fun;
                        x0, x1, gdth, minΔx, maxΔx,
                        break_cond = vg_break_cond,
                        damp_factor, damp = vg_gddamp,
                        target, maxiter = gdmaxiter,
                        verbose = false
                    )

                    vg_beta = SimT.gd_value(gdmodel)
                    vg_gddamp = gdmodel.damp
                end

                ## -------------------------------------------------------------------
                # COLLECTING
                push!(biom_betas, biom_beta)
                push!(vg_betas, vg_beta)
                epouts[(biom_beta, vg_beta)] = epout

                ## -------------------------------------------------------------------
                # CHECK NAN
                if check_nan()
                    print_info("Nan detected (BIG PROBLEMS HERE)";
                        exp, rounditer
                    ); break
                end

                ## -------------------------------------------------------------------
                # MONITOR
                SimT.record!(mon) do dat
                    tdat = get!(dat, exp, Dict())
                    tdat[:live_prove] = rand()
                    rdat = get!(tdat, :round, Dict())
                    SimT.get!push!(rdat;
                        vg_beta, biom_beta,
                        biom_avPME, vg_avPME
                    )
                end

                ## -------------------------------------------------------------------
                # PRINT INFO
                print_info("Round Done"; exp, rounditer,
                    hasvalid_moments, roundconv
                )

                ## -------------------------------------------------------------------
                # BREAK
                roundconv && break
                rounditer += 1
                rounditer > maxrounds && break

            end # round while

            ## -------------------------------------------------------------------
            # Storing
            lock(WLOCK) do

                dat = Dict()
                dat[:exp_beta] = (biom_beta, vg_beta)
                dat[:epouts] = epouts
                dat[:model] = model |> ChU.compressed_model

                datfile = dat_file(;method, exp)
                sdat(iJR, dat, datfile)

                print_info("Finished "; exp, rounditer)
            end

        end # for exp, cGLC
    end # for thid

    SimT.reset!(mon)

end

ext: jl


text:
function solve(prob::NonlinearProblem{<:Number}, ::NewtonRaphson, args...; xatol = nothing, xrtol = nothing, maxiters = 1000, kwargs...)
  f = Base.Fix2(prob.f, prob.p)
  x = float(prob.u0)
  T = typeof(x)
  atol = xatol !== nothing ? xatol : oneunit(T) * (eps(one(T)))^(4//5)
  rtol = xrtol !== nothing ? xrtol : eps(one(T))^(4//5)

  xo = oftype(x, Inf)
  for i in 1:maxiters
    fx, dfx = value_derivative(f, x)
    iszero(fx) && return NewtonSolution(x, DEFAULT)
    Δx = dfx \ fx
    x -= Δx
    if isapprox(x, xo, atol=atol, rtol=rtol)
        return NewtonSolution(x, DEFAULT)
    end
    xo = x
  end
  return NewtonSolution(x, MAXITERS_EXCEED)
end

function scalar_nlsolve_ad(prob, alg, args...; kwargs...)
  f = prob.f
  p = value(prob.p)
  u0 = value(prob.u0)

  newprob = NonlinearProblem(f, u0, p; prob.kwargs...)
  sol = solve(newprob, alg, args...; kwargs...)

  uu = getsolution(sol)
  if p isa Number
    f_p = ForwardDiff.derivative(Base.Fix1(f, uu), p)
  else
    f_p = ForwardDiff.gradient(Base.Fix1(f, uu), p)
  end

  f_x = ForwardDiff.derivative(Base.Fix2(f, p), uu)
  pp = prob.p
  sumfun = let f_x′ = -f_x
    ((fp, p),) -> (fp / f_x′) * ForwardDiff.partials(p)
  end
  partials = sum(sumfun, zip(f_p, pp))
  return sol, partials
end

function solve(prob::NonlinearProblem{<:Number, iip, <:Dual{T,V,P}}, alg::NewtonRaphson, args...; kwargs...) where {iip, T, V, P}
  sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
  return NewtonSolution(Dual{T,V,P}(sol.u, partials), sol.retcode)
end
function solve(prob::NonlinearProblem{<:Number, iip, <:AbstractArray{<:Dual{T,V,P}}}, alg::NewtonRaphson, args...; kwargs...) where {iip, T, V, P}
  sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
  return NewtonSolution(Dual{T,V,P}(sol.u, partials), sol.retcode)
end

# avoid ambiguities
for Alg in [Bisection, Falsi]
  @eval function solve(prob::NonlinearProblem{uType, iip, <:Dual{T,V,P}}, alg::$Alg, args...; kwargs...) where {uType, iip, T, V, P}
    sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
    return BracketingSolution(Dual{T,V,P}(sol.left, partials), Dual{T,V,P}(sol.right, partials), sol.retcode)
  end
  @eval function solve(prob::NonlinearProblem{uType, iip, <:AbstractArray{<:Dual{T,V,P}}}, alg::$Alg, args...; kwargs...) where {uType, iip, T, V, P}
    sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
    return BracketingSolution(Dual{T,V,P}(sol.left, partials), Dual{T,V,P}(sol.right, partials), sol.retcode)
  end
end

function solve(prob::NonlinearProblem, ::Bisection, args...; maxiters = 1000, kwargs...)
  f = Base.Fix2(prob.f, prob.p)
  left, right = prob.u0
  fl, fr = f(left), f(right)

  if iszero(fl)
    return BracketingSolution(left, right, EXACT_SOLUTION_LEFT)
  end

  i = 1
  if !iszero(fr)
    while i < maxiters
      mid = (left + right) / 2
      (mid == left || mid == right) && return BracketingSolution(left, right, FLOATING_POINT_LIMIT)
      fm = f(mid)
      if iszero(fm)
        right = mid
        break
      end
      if sign(fl) == sign(fm)
        fl = fm
        left = mid
      else
        fr = fm
        right = mid
      end
      i += 1
    end
  end

  while i < maxiters
    mid = (left + right) / 2
    (mid == left || mid == right) && return BracketingSolution(left, right, FLOATING_POINT_LIMIT)
    fm = f(mid)
    if iszero(fm)
      right = mid
      fr = fm
    else
      left = mid
      fl = fm
    end
    i += 1
  end

  return BracketingSolution(left, right, MAXITERS_EXCEED)
end


text:
using Interpolations

export plot_h, plot_dist, plot_dist_control, plot_dist_treat, counterfactual_discrete, counterfactual_continuous

function plot_h(f::FELT)

    x1 = f.gamma_1_hat
    y1 = f.y1s
    x2 = f.gamma_2_hat
    y2 = f.y2s

    if f.discrete
        plot(x1,y1-0.1, linetype = :steppost, linewidth = 5)
        plot!(x2,y2+0.1, linetype = :steppost, linewidth = 5)
    else
        plot(x1,y1, linewidth = 2)
        plot!(x2,y2, linewidth = 2)
    end

    xmin = min(findmin(x1)[1],findmin(x2)[1])
    xmax = max(findmax(x1)[1],findmax(x2)[1])
    extend = 0.1*(xmax-xmin)
    plot!(xlims=(xmin-extend,xmax+extend))

    yvals = sort(union(y1,y2))
    ymin = findmin(yvals)[1]
    ymax = findmax(yvals)[1]
    #yticks!(yvals)
    ylims!(ymin-1,ymax+1)
    p = plot!(legend=:topleft)

    return(p)

end

function plot_dist(f::FELT)
    return histogram2d(f.y[:,1],f.y[:,2], bins = f.y1s)
end

function plot_h(f::NLDID)
    plot_h(f.felt)
    # Just redirect to the plot function for a FELT object.
end

function plot_dist_control(f::NLDID)
    # Delegate to the function for FELT.
    plot_dist(f.felt)
end

function plot_dist_treat(f::NLDID)
    y1 = f.y_treat[f.t.==0]
    y2 = f.y_treat[f.t.==1]
    histogram2d(y1,y2, bins = f.felt.y1s)
end

#############
# Get the treatment effect.
#############

function counterfactual_discrete(Y0,YLOW,YHIGH,gamma_1,y1,gamma_2,y2)

    # Y0 is a vector of observed values in period 1 for which we want to compute
    #       counterfactuals in time 2.
    #
    # Should be passed in as
    #     Y0 = f.y_treat[f.t .== 0]
    # where f is a NLDID object.

    # YLOW is the lowest value the variable can take
    # YHIGH is the highest value the variable can take

    # y1 is the set of cut points used in period 1. These will exclude YLOW.
    # gamma_1 is the estimated thresholds on the latent variable, e.g. y*>=gamma_1[2]
    #    means y>=y1[2].
    #
    # same for y2, gamma_2.
    #
    #
    # Start by computing bounds on the latent variable.

    # lb, ub represent the lower and upper bound on the latent variable
    #   for period 1. They are vectors of the same size as Y0.
    lb = similar(Y0)*0.
    ub = similar(Y0)*0.
    n_treat = length(Y0)

    for i in 1:n_treat

        index = sum((Y0[i] .>= y1))

        if index == 0
            lb[i] = -Inf
        else
            lb[i] = gamma_1[index]
        end

        if index == length(y1)
            ub[i] = Inf
        else
            ub[i] = gamma_1[index+1]
        end
    end

    # Now convert them into bounds on the outcome variable.

    # First the lower bounds.
    lb_counterfactual = similar(Y0)*0.
    ub_counterfactual = similar(Y0)*0.

    for i in 1:n_treat

        index = sum(lb[i] .>= gamma_2)
        if index == 0
            lb_counterfactual[i] = YLOW
        else
            lb_counterfactual[i] = y2[index]
        end

        index = sum(ub[i] .>= gamma_2)
        if index == 0
            ub_counterfactual[i] = YLOW
        else
            ub_counterfactual[i] = y2[index]
        end
    end

    return mean(ub_counterfactual), mean(lb_counterfactual)

end

# Write it up as a function
function counterfactual_continuous(Yobs,Xobs,xto,beta,y_from,c_from,y_to,c_to,YLO,YHI)

    n,K = size(Xobs)

    if K>0 #deal with situations where no X is present.
        dXb = (xto.-Xobs)*beta
    else
        dXb = 0
    end

    # The following interpolation objects will serve as functions.
    # h1inv is an interpolation object that is the inverse of
    #   the transformation function in the time period from which
    #   Yobs originates.
    h1inv = interpolate((y_from,),c_from, Gridded(Linear()))
    # h2 is an interpolation object that allows you to go
    #   from latent variable to outcome in period "to", for which we're
    #   constructing a counterfactual.
    h2 = interpolate((c_to,),y_to, Gridded(Linear()))

    # Map it!
    Y_cf = h2[h1inv[Yobs] + dXb]

    # Stay in bounds, for example while censoring.
    Y_cf[Y_cf .< YLO] = YLO
    Y_cf[Y_cf .> YHI] = YHI
    return Y_cf
end


text:
using Test

using LinkageDisequilibrium

samples = Bool[
    0 1 1 1 0 1 0 1
    0 1 0 1 0 1 0 1
    1 1 0 1 0 1 0 0
    1 1 0 1 1 1 0 1
    1 1 1 1 1 1 1 1
    0 0 1 1 1 1 0 0
    0 0 0 0 0 0 0 0
]

expected = [2, 5, 5, 6, 7, 7, 7]
actual = [
    length(ld_select(samples; r²_threshold = τ))
    for τ in 0.2:0.1:0.8
]

@test expected == actual


text:
using BioBackgroundModels,BioSequences,DataFrames,Distributed,Distributions,FASTX,ProgressMeter,Random,StatsFuns,Test
import BioBackgroundModels: autotransition_init, load_balancer, CompoundAlphabet, get_order_n_seqs, code_seqs, code_job_obs, make_padded_df, add_partition_masks!, partition_genome_coordinates, divide_partitions_by_scaffold, mask_sequence_by_partition, find_position_partition, setup_sample_jobs, rectify_identifier, add_metacoordinates!, meta_to_feature_coord, get_feature_row_index, get_strand_dict, build_scaffold_seq_dict, get_feature_params_from_metacoord, determine_sample_window, fetch_sequence, get_sample_set, get_sample, assert_hmm, linear_step, get_BGHMM_symbol_lh,make_padded_df,add_partition_masks!,BGHMM_likelihood_calc,linear_likelihood, t_add_categorical_counts!
import Serialization: deserialize
include("synthetic_sequence_gen.jl")
include("ref_fns.jl")

#JOB FILEPATHS
#GFF3 feature database, FASTA genome and index paths
Sys.islinux() ? genome =  (@__DIR__) * "/synthetic.fna" : genome = (@__DIR__) * "\\synthetic.fna"
Sys.islinux() ? index =  (@__DIR__) * "/synthetic.fna.fai" : index = (@__DIR__) * "\\synthetic.fna.fai"
Sys.islinux() ? gff =  (@__DIR__) * "/synthetic.gff3" : gff = (@__DIR__) * "\\synthetic.gff3"
Sys.islinux() ? posfasta =  (@__DIR__) * "/syntheticpos.fa" : posfasta = (@__DIR__) * "\\syntheticpos.fa"

!isfile(genome) && print_synthetic_fasta(genome)
!isfile(index) && print_synthetic_index(index)
!isfile(gff) && print_synthetic_gff(gff)
!isfile(posfasta) && print_synthetic_position(posfasta)

Random.seed!(1)

@testset "BHMM/BHMM.jl constructors" begin
    good_a = [.20,.30,.20,.30]
    bad_probvec_a = [.25,.25,.25,.27]
    good_A = fill(.25,4,4)
    bad_probvec_A = deepcopy(good_A)
    bad_probvec_A[:,3].=.27
    bad_size_A = fill(.25,5,4)
    good_D = Categorical([.25,.25,.25,.25])
    bad_size_D = Categorical([.2,.2,.2,.2,.2])
    good_Dvec = [good_D for i in 1:4]
    bad_size_Dvec =[good_D, good_D, good_D, bad_size_D]
    bad_length_Dvec=[good_D for i in 1:5]

    @test_throws ArgumentError BHMM(bad_probvec_a, good_A, good_Dvec)
    @test_throws ArgumentError BHMM(good_a, bad_probvec_A, good_Dvec)
    @test_throws ArgumentError BHMM(good_a, bad_size_A, good_Dvec)
    @test_throws ArgumentError BHMM(good_a, good_A, bad_size_Dvec)
    @test_throws ArgumentError BHMM(good_a, good_A, bad_length_Dvec)

    hmm=BHMM(good_a, good_A, good_Dvec)
    @test typeof(hmm)==BHMM{Float64}
    @test hmm.a==good_a
    @test hmm.A==good_A
    @test hmm.B==good_Dvec
    @test size(hmm) == (4,4)

    hmm_copy=copy(hmm)
    @test hmm.a==hmm_copy.a
    @test hmm.A==hmm_copy.A
    @test hmm.B==hmm_copy.B

    hmm2=BHMM(good_A, good_Dvec)
    @test hmm2.a==[.25,.25,.25,.25]
    @test hmm2.A==good_A
    @test hmm2.B==good_Dvec
end

@testset "BHMM/chain.jl Chain_ID and EM_step constructors" begin
    obs_id="test"
    K,zero_k,neg_K=4,0,-1
    order,badorder=0,-1
    replicate,zero_rep,neg_rep=1,0,-1

    @test_throws ArgumentError Chain_ID(obs_id, zero_k, order, replicate)
    @test_throws ArgumentError Chain_ID(obs_id, neg_K, order, replicate)
    @test_throws ArgumentError Chain_ID(obs_id, K, badorder, replicate)
    @test_throws ArgumentError Chain_ID(obs_id, K, order, zero_rep)
    @test_throws ArgumentError Chain_ID(obs_id, K, order, neg_rep)

    id=Chain_ID(obs_id, K, order, replicate)
    @test typeof(id)==Chain_ID
    @test id.obs_id==obs_id
    @test id.K==K
    @test id.order==order
    @test id.replicate==replicate

    good_A = fill(.25,4,4)
    good_D = Categorical([.25,.25,.25,.25])
    good_Dvec = [good_D for i in 1:4]
    hmm=BHMM(good_A,good_Dvec)
    log_p=-.001

    @test_throws ArgumentError EM_step(0,hmm,0.0,0.0,true)
    @test_throws ArgumentError EM_step(-1,hmm,0.0,0.0,true)
    @test_throws ArgumentError EM_step(1,hmm,1.0,0.0,true)

    steppe=EM_step(1,hmm,log_p,0.0,false)
    @test typeof(steppe)==EM_step
    @test steppe.iterate==1
    @test steppe.hmm==hmm
    @test steppe.log_p==log_p
    @test steppe.delta==0.0
    @test steppe.converged==false
end

@testset "utilities/log_prob_sum.jl Log probability summation functions" begin
    @test isapprox(lps([.1, .1, .1]), .3)
    @test lps([-Inf, .1, .1]) == -Inf
    @test isapprox(lps(.1, .1, .1), .3)
    @test lps(-Inf, .1, .1) == -Inf
    @test lps(.1, -Inf, .1) == -Inf
end

@testset "utilities/HMM_init.jl initialization functions" begin
    K=4
    order=2
    hmm=autotransition_init(K,order)
    @test typeof(hmm)==BHMM{Float64}
    @test size(hmm)==(4,64)
end

@testset "utilities/load_balancer.jl EM_converge load balancing structs and functions" begin
    @test_throws ArgumentError LoadConfig(0:5,1:5)
    @test_throws ArgumentError LoadConfig(1:5,-1:5)

    job_ids=[Chain_ID("test",6,2,1),Chain_ID("test",4,2,1),Chain_ID("test",2,0,1)]
    obs_sets=Dict("test"=>[LongSequence{DNAAlphabet{2}}("AAAAAAAAAAAAAA")])

    K_exclusive_config=LoadConfig(1:1,0:0)
    O_exclusive_config=LoadConfig(1:6,3:4)
    blacklist_config=LoadConfig(1:6,0:2,blacklist=job_ids)
    whitelist_config=LoadConfig(1:6,0:2,whitelist=[Chain_ID("test",6,2,2)])
    high_config=LoadConfig(4:6,0:2,whitelist=job_ids)

    no_input_hmms, chains, input_channel, output_channel = setup_EM_jobs!(job_ids,obs_sets)
    @test load_balancer(no_input_hmms, input_channel, K_exclusive_config, 3) == (0,0,0,0,0)

    no_input_hmms, chains, input_channel, output_channel = setup_EM_jobs!(job_ids,obs_sets)
    @test load_balancer(no_input_hmms, input_channel, O_exclusive_config, 3) == (0,0,0,0,0)

    no_input_hmms, chains, input_channel, output_channel = setup_EM_jobs!(job_ids,obs_sets)
    @test load_balancer(no_input_hmms, input_channel, blacklist_config, 3) == (0,0,0,0,0)

    no_input_hmms, chains, input_channel, output_channel = setup_EM_jobs!(job_ids,obs_sets)
    @test load_balancer(no_input_hmms, input_channel, whitelist_config, 3) == (0,0,0,0,0)

    no_input_hmms, chains, input_channel, output_channel = setup_EM_jobs!(job_ids,obs_sets)
    high_set=load_balancer(no_input_hmms, input_channel, high_config)

    @test high_set[1]==Chain_ID("test",6,2,1)
    @test high_set[2]==1
    @test typeof(high_set[3])==BHMM{Float64}
    @test high_set[4]==obs_lh_given_hmm(high_set[5],high_set[3],linear=false)
    @test typeof(high_set[5])==Matrix{UInt8}
end

@testset "utilities/observation_coding.jl Order coding structs and functions" begin
    compound_alphabet=CompoundAlphabet(ACGT, 2)
    @test typeof(compound_alphabet.symbols)==Dict{Mer{DNAAlphabet{2}},Integer}
    @test length(compound_alphabet.symbols)==64
    @test compound_alphabet.symbols[Mer{DNAAlphabet{2}}("AAA")]==1
    @test compound_alphabet.symbols[Mer{DNAAlphabet{2}}("TTT")]==64

    test_seqs = [BioSequences.LongSequence{DNAAlphabet{2}}("ACGTACGTACGTACGT"),BioSequences.LongSequence{DNAAlphabet{2}}("TTTTTTT")]
    target0= [1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 0; 4 4 4 4 4 4 4 0 0 0 0 0 0 0 0 0 0]
    target2= [37 58 15 20 37 58 15 20 37 58 15 20 37 58 0; 64 64 64 64 64 0 0 0 0 0 0 0 0 0 0]

    order0_seqs = get_order_n_seqs(test_seqs,0)
    @test typeof(order0_seqs)==BioBackgroundModels.N_Order_ntSequence
    for (mer, code) in CompoundAlphabet(ACGT,0).symbols
        @test order0_seqs.alphabet.symbols[mer]==code
    end
    @test order0_seqs.seq_lengths==[16,7]
    @test order0_seqs.order_kmers[1][1]==Mer{DNAAlphabet{2}}("A")
    code0_seqs = code_seqs(order0_seqs)
    @test target0 == code0_seqs

    order2_seqs = get_order_n_seqs(test_seqs,2)
    @test typeof(order2_seqs)==BioBackgroundModels.N_Order_ntSequence
    for (mer, code) in CompoundAlphabet(ACGT,2).symbols
        @test order2_seqs.alphabet.symbols[mer]==code
    end
    @test order2_seqs.seq_lengths==[16,7]
    @test order2_seqs.order_kmers[1][1]==Mer{DNAAlphabet{2}}("ACG")
    code2_seqs = code_seqs(order2_seqs)
    @test target2 == code2_seqs

    bw_o0_seqs = get_order_n_seqs(test_seqs[end:-1:1], 0)
    sorted_code_seqs = code_seqs(bw_o0_seqs, sorted=true)
    @test target0 == sorted_code_seqs

    obs_sets=Dict("test1"=>test_seqs,
                    "test2"=>test_seqs)
    job_ids=[Chain_ID("test1",4,0,1),Chain_ID("test2",4,0,1),Chain_ID("test1",2,2,1)]

    code_dict=code_job_obs(job_ids,obs_sets)
    @test ("test1",0) in keys(code_dict)
    @test ("test1",2) in keys(code_dict)
    @test ("test2",0) in keys(code_dict)
    @test !(("test2",2) in keys(code_dict))
    @test code_dict[("test1",0)]==target0
    @test code_dict[("test1",2)]==target2
    @test code_dict[("test2",0)]==target0

    #test type selection for high order numbers
    o5_nos=get_order_n_seqs(test_seqs,5)
    o5_codes=code_seqs(o5_nos)
    @test typeof(o5_codes)==Matrix{UInt16}

    o7_nos=get_order_n_seqs(test_seqs,7)
    o7_codes=code_seqs(o7_nos)
    @test typeof(o7_codes)==Matrix{UInt32}

end

@testset "genome_sampling/partition_masker.jl functions" begin
    synthetic_seq = BioSequences.LongSequence{DNAAlphabet{2}}(generate_synthetic_seq())
    position_start = 501
    position_length=141
    position_pad=350
    perigenic_pad = 250

    position_df = make_padded_df(posfasta, gff, genome, index, position_pad)
    add_partition_masks!(position_df, gff, perigenic_pad)

    @test position_df.SeqID[1]=="1"
    @test length(position_df.PadSeq[1]) == position_pad+position_length == length(position_df.PadStart[1]:position_df.PadStart[1]+position_length+position_pad-1) == position_df.End[1]-position_df.Start[1]+1+position_pad == size(position_df.MaskMatrix[1])[1]
    @test synthetic_seq[position_df.PadStart[1]:position_df.End[1]]==position_df.PadSeq[1]

    mm = position_df.MaskMatrix[1]
    idx=1#pos 151
    @test sum(mm[idx:idx+99,1] .== 1)==length(mm[idx:idx+99,1])
    idx+=100#pos 251
    @test sum(mm[idx:idx+259,1] .== 2)==length(mm[idx:idx+259,1])
    idx+=260 #pos 511
    @test sum(mm[idx:idx+59,1] .== 3)==length(mm[idx:idx+59,1])
    idx+=60#pos 571
    @test sum(mm[idx:idx+29,1] .== 2)==length(mm[idx:idx+29,1])
    idx+=30#pos601
    @test sum(mm[idx:idx+40,1] .== 3)==length(mm[idx:idx+40,1])

    #test exception
    position_df = make_padded_df(posfasta, gff, genome, index, position_pad)

    partition_coords_dict=partition_genome_coordinates(gff, perigenic_pad)
    partitioned_scaffolds=divide_partitions_by_scaffold(partition_coords_dict)
    scaffold_coords_dict = Dict{String,DataFrame}()
    scaffold="1"
    for ((partition, part_scaffold), df) in partitioned_scaffolds
        if scaffold == part_scaffold
            scaffold_coords_dict[partition] = df
        end
    end

    @test_throws DomainError mask_sequence_by_partition(1001,151,scaffold_coords_dict)
end

@testset "genome_sampling/sequence_sampler.jl functions & API/genome_sampling.jl interface" begin
    #rectifier tests
    @test rectify_identifier("1")=="1"

    partitions = 3 #exonic, periexonic, intragenic
    sample_set_length=100
    min_sample_window=5
    max_sample_window=25
    perigenic_pad=250
    syn_intergenic_starts = [1]
    syn_intergenic_ends = [250]
    syn_periexonic_starts = [251,571,661,761,861,961]
    syn_periexonic_ends = [510,600,700,800,900,1000]
    syn_periexonic_strands = ['-','-','-','-','-','-']
    syn_exonic_starts = [511,601,701,801,901]
    syn_exonic_ends = [570,660,760,860,960]
    syn_exonic_strands = ['+','-','-','-','-']
    partition_lengths = Dict("exon"=>5*60,"intergenic"=>250,"periexonic"=>450)

    # "Testing sequence sampler fns..."
    synthetic_seq = BioSequences.LongSequence{DNAAlphabet{2}}(generate_synthetic_seq())

    # "Partitioning synthetic genome coordinates..."
    coordinate_partitions = partition_genome_coordinates(gff, perigenic_pad)

    @test coordinate_partitions["intergenic"].Start == syn_intergenic_starts
    @test coordinate_partitions["intergenic"].End == syn_intergenic_ends

    @test coordinate_partitions["periexonic"].Start == syn_periexonic_starts
    @test coordinate_partitions["periexonic"].End == syn_periexonic_ends
    @test coordinate_partitions["periexonic"].Strand == syn_periexonic_strands

    @test coordinate_partitions["exon"].Start == syn_exonic_starts
    @test coordinate_partitions["exon"].End == syn_exonic_ends
    @test coordinate_partitions["exon"].Strand == syn_exonic_strands

    # "Checking sampling functions at all synthetic indices..."

    input_test_channel, completed_test_channel, progress_channel, setlength = setup_sample_jobs(genome, index, gff, sample_set_length, min_sample_window, max_sample_window, perigenic_pad; deterministic=true)
    while isready(input_test_channel)
        genome_path, genome_index_path, partition_df, partitionid, sample_set_length, sample_window_min, sample_window_max, deterministic = take!(input_test_channel)

        for feature in eachrow(partition_df)
            seqid, scaffold_start, scaffold_end = meta_to_feature_coord(feature.MetaStart, feature.MetaEnd, partition_df)
            @test scaffold_start == feature.Start && scaffold_end == feature.End
        end

        stranded = get_strand_dict()[partitionid]
        @test typeof(stranded) == Bool

        scaffold_sequence_record_dict = build_scaffold_seq_dict(genome, index)
        @test scaffold_sequence_record_dict["1"] == synthetic_seq

        partition_length = partition_df.MetaEnd[end]
        @test partition_length == partition_lengths[partitionid]

        metacoordinate_bitarray = trues(partition_df.MetaEnd[end])

        for bitindex in findall(metacoordinate_bitarray)
            feature_metaStart, feature_metaEnd, strand = get_feature_params_from_metacoord(bitindex, partition_df, stranded)
            @test 1 <= feature_metaStart < feature_metaEnd <= length(metacoordinate_bitarray)
            @test feature_metaStart in partition_df.MetaStart
            @test feature_metaEnd in partition_df.MetaEnd
            if partitionid == "periexonic"
                @test strand == '-'
            elseif partitionid == "exon"
                @test strand in ['-','+']
            end
            feature_length = length(feature_metaStart:feature_metaEnd)
            window = determine_sample_window(feature_metaStart, feature_metaEnd, bitindex, metacoordinate_bitarray, sample_window_min, sample_window_max) #get an appropriate sampling window around the selected index, given the feature boundaries and params
            @test 1 <= feature_metaStart <= window[1] < window[1]+sample_window_min-1 < window[1]+sample_window_max-1<= window[2] <= feature_metaEnd <= length(metacoordinate_bitarray)
            sample_scaffoldid, sample_scaffold_start, sample_scaffold_end = meta_to_feature_coord(window[1],window[2],partition_df)
            @test sample_scaffoldid == "1"
            @test 1 <= sample_scaffold_start <= sample_scaffold_start+sample_window_min <= sample_scaffold_end <= min(sample_scaffold_start+sample_window_max,1000) <= 1000

            strand == '-' ? target_seq=reverse_complement(synthetic_seq[sample_scaffold_start:sample_scaffold_end]) :
                target_seq = synthetic_seq[sample_scaffold_start:sample_scaffold_end]

            proposal_sequence = fetch_sequence(sample_scaffoldid, scaffold_sequence_record_dict, sample_scaffold_start, sample_scaffold_end, strand; deterministic=deterministic) #get the sequence associated with the sample window

            @test proposal_sequence == target_seq
        end
    end

    # "Verifying sampling channels..."

    input_sample_channel, completed_sample_channel, progress_channel, setlength = setup_sample_jobs(genome, index, gff, sample_set_length, min_sample_window, max_sample_window, perigenic_pad; deterministic=true)
    get_sample_set(input_sample_channel, completed_sample_channel, progress_channel)

    #collect sample dfs by partition id when ready
    collected_counter = 0
    sample_record_dfs = Dict{String,DataFrame}()
    while collected_counter < partitions
        wait(completed_sample_channel)
        partition_id, sample_df = take!(completed_sample_channel)
        sample_record_dfs[partition_id] = sample_df
        collected_counter += 1
    end

    #get_sample entire feature selection
    coordinate_partitions = partition_genome_coordinates(gff, perigenic_pad)
    partition_id="exon"
    partition=coordinate_partitions[partition_id]
    add_metacoordinates!(partition)
    metacoordinate_bitarray=trues(partition.MetaEnd[end])
    scaffold_sequence_record_dict = build_scaffold_seq_dict(genome, index)

    #test fetch_sequence strand char ArgumentError
    @test_throws ArgumentError fetch_sequence("1", Dict{String,LongSequence}("1"=>LongSequence{DNAAlphabet{2}}("AAAAAAAAAAAAA")), 1, 5, 'L')

    #test get_feature_row_index ArgumentError
    @test_throws DomainError get_feature_row_index(partition,0)

    @test length(get_sample(metacoordinate_bitarray, 1, 250, partition, scaffold_sequence_record_dict)[6])==60

    synthetic_seq = BioSequences.LongSequence{DNAAlphabet{2}}(generate_synthetic_seq())
    for (partid, df) in sample_record_dfs
        for sample in eachrow(df)
            target_seq = synthetic_seq[sample.SampleStart:sample.SampleEnd]
            strand = sample.Strand
            if sample.Strand == '-'
                target_seq = reverse_complement(target_seq)
            end
            @test sample.SampleSequence == target_seq
        end
    end

    #execute_sample_jobs API
    wkpool=addprocs(1)
    @everywhere using BioBackgroundModels

    channels = setup_sample_jobs(genome, index, gff, sample_set_length, min_sample_window, max_sample_window, perigenic_pad; deterministic=true)
    records=execute_sample_jobs(channels, wkpool)

    rmprocs(wkpool)
end

@testset "likelihood_funcs/hmm.jl & bg_lh_matrix.jl functions" begin
    pvec = [.4,.3,.2,.1]
    trans = ones(1,1)
    B = [Categorical(pvec)]
    hmm = BHMM(trans, B)

    testseq=zeros(Int64,1,5)
    testseq[1:4] = [1,2,3,4]
    @test isapprox(get_BGHMM_symbol_lh(testseq, hmm)[1], log.(pvec[1]))
    @test isapprox(obs_lh_given_hmm(testseq,hmm),obs_lh_given_hmm(testseq,hmm,linear=false))
    @test isapprox(sum(get_BGHMM_symbol_lh(testseq,hmm)), obs_lh_given_hmm(testseq,hmm))
    @test isapprox(sum(get_BGHMM_symbol_lh(testseq,hmm)), obs_lh_given_hmm(testseq,hmm,linear=false))

    pvec=[.25,.25,.25,.25]
    trans = [.9 .1
         .1 .9]
    B = [Categorical(pvec), Categorical(pvec)]
    hmm = BHMM(trans, B)

    @test isapprox(obs_lh_given_hmm(testseq,hmm),obs_lh_given_hmm(testseq,hmm,linear=false))
    @test isapprox(sum(get_BGHMM_symbol_lh(testseq,hmm)), obs_lh_given_hmm(testseq,hmm))
    @test isapprox(sum(get_BGHMM_symbol_lh(testseq,hmm)), obs_lh_given_hmm(testseq,hmm,linear=false))

    testseq=zeros(Int64,1,1001)
    testseq[1,1:1000]=rand(1:4,1000)

    @test isapprox(sum(get_BGHMM_symbol_lh(testseq, hmm)),obs_lh_given_hmm(testseq,hmm))

    Dex = [Categorical([.3, .1, .3, .3]),Categorical([.15, .35, .35, .15])]
    Dper = [Categorical([.15, .35, .35, .15]),Categorical([.4, .1, .1, .4])]
    Dint = [Categorical([.4, .1, .1, .4]),Categorical([.45, .05, .05, .45])]
    BGHMM_dict = Dict{String,BHMM}()
    BGHMM_dict["exon"] = BHMM(trans, Dex)
    BGHMM_dict["periexonic"] = BHMM(trans, Dper)
    BGHMM_dict["intergenic"] = BHMM(trans, Dint)

    position_length=141;perigenic_pad=250;
    position_df = make_padded_df(posfasta, gff, genome, index, position_length)

    reader=open(FASTA.Reader, genome, index=index)
    seq=FASTA.sequence(reader["CM002885.2.1"])
    @test position_df.Start[1]==501
    @test position_df.End[1]==641
    @test position_df.PadSeq[1]==seq[360:641]
    @test position_df.PadStart[1]==360
    @test position_df.RelStart[1]==141
    @test position_df.SeqOffset[1]==0

    offset_position_df = make_padded_df(posfasta, gff, genome, index,600)
    @test offset_position_df.Start[1]==501
    @test offset_position_df.End[1]==641
    @test offset_position_df.PadSeq[1]==seq[1:641]
    @test offset_position_df.PadStart[1]==1
    @test offset_position_df.RelStart[1]==500
    @test offset_position_df.SeqOffset[1]==100

    add_partition_masks!(position_df, gff, perigenic_pad)
    @test all(position_df.MaskMatrix[1][1:151,1].==2)
    @test all(position_df.MaskMatrix[1][152:211,1].==3)
    @test all(position_df.MaskMatrix[1][212:241,1].==2)
    @test all(position_df.MaskMatrix[1][242:end,1].==3)
    @test all(position_df.MaskMatrix[1][1:151,2].==-1)
    @test all(position_df.MaskMatrix[1][152:211,2].==1)
    @test all(position_df.MaskMatrix[1][212:end,2].==-1)

    lh_matrix=BGHMM_likelihood_calc(position_df,BGHMM_dict)
    @test size(lh_matrix)==(position_length*2,1)

    periexonic_frag=reverse_complement!(LongSequence{DNAAlphabet{2}}(seq[360:510]))
    pno=get_order_n_seqs([periexonic_frag],0)
    pcode=code_seqs(pno)
    plh=get_BGHMM_symbol_lh(pcode, BGHMM_dict["periexonic"])
    @test isapprox(lh_matrix[1:151,1], reverse(plh))

    exonic_frag=LongSequence{DNAAlphabet{2}}(seq[511:570])
    eno=get_order_n_seqs([exonic_frag],0)
    ecode=code_seqs(eno)
    elh=get_BGHMM_symbol_lh(ecode, BGHMM_dict["exon"])
    @test isapprox(lh_matrix[152:211,1],elh)
end

@testset "EM/baum-welch.jl MS_HMMBase equivalency and functions" begin
    A = [.5 .5
         .5 .5]
    B = [Categorical(ones(4)/4), Categorical([.7,.1,.1,.1])]
    hmm = BHMM(A, B)
    log_A = log.(hmm.A)

    obs = zeros(UInt8, 22,1)
    obs[1:21] = [4,3,3,2,3,2,1,1,2,3,3,3,4,4,2,3,2,3,4,3,2]
    obs_lengths=[21]

    lls = BioBackgroundModels.bw_llhs(hmm,obs)
    m_lls = mouchet_log_likelihoods(hmm,obs[1:21])

    K,Tmaxplus1,O = size(lls)
    T=Tmaxplus1-1

    #TEST LOG LIKELIHOOD FN
    @test transpose(lls[:,1:end-1,1]) == m_lls

    log_α = BioBackgroundModels.messages_forwards_log(hmm.a, hmm.A, lls,  obs_lengths)
    m_log_α = mouchet_messages_forwards_log(hmm.a, hmm.A, m_lls)

    #TEST FORWARD MESSAGES FN
    @test transpose(log_α[:,1:end-1,1]) == m_log_α

    log_β = BioBackgroundModels.messages_backwards_log(hmm.A, lls,  obs_lengths)
    m_log_β = mouchet_messages_backwards_log(hmm.A, m_lls)

    #TEST BACKWARDS  MESSAGES FN
    @test isapprox(transpose(log_β[:,1:end-1,1]),  m_log_β)

    lls = permutedims(lls, [2,3,1]) # from (K,T,O) to (T,O,K)
    log_α = permutedims(log_α, [2,3,1])
    log_β = permutedims(log_β, [2,3,1])

    # "Testing obs probability calcs..."
    #TEST OBSERVATION PROBABILITY CALC
    normalizer = logsumexp(m_log_α[1,:] + m_log_β[1,:])
    log_pobs = logsumexp(lps.(log_α[1,1,:], log_β[1,1,:]))

    @test isapprox(normalizer, log_pobs)

    # "Testing ξ, γ calcs..."

    log_ξ = fill(-Inf, Tmaxplus1,O,K,K)
    log_γ = fill(-Inf, Tmaxplus1,O,K)

    m_log_ξ = zeros(T, K, K)

    for t = 1:T-1, i = 1:K, j = 1:K
        m_log_ξ[t,i,j] = m_log_α[t,i] + log_A[i,j] + m_log_β[t+1,j] + m_lls[t+1,j] - normalizer
    end

    for i = 1:K, j = 1:K, o = 1:O
        obsl = obs_lengths[o]
        for t = 1:obsl-1 #log_ξ & log_γ calculated to T-1 for each o
           log_ξ[t,o,i,j] = lps(log_α[t,o,i], log_A[i,j], log_β[t+1,o,j], lls[t+1,o,j], -log_pobs[o])
           log_γ[t,o,i] = lps(log_α[t,o,i], log_β[t,o,i], -log_pobs[o])
        end
        t=obsl #log_ξ @ T = 0
        log_ξ[t,o,i,j] = 0
        log_γ[t,o,i] = lps(log_α[t,o,i], log_β[t,o,i], -log_pobs)
    end

    ξ = exp.(log_ξ)
    k_ξ = sum(ξ, dims=[3,4])
    nan_mask = k_ξ .== 0; k_ξ[nan_mask] .= Inf #prevent NaNs in dummy renorm arising from multiple sequences indexing
    ξ  ./=  k_ξ #dummy renorm across K to keep numerical creep from causing isprobvec to fail on new new_A during hmm creation

    m_ξ = exp.(m_log_ξ)
    m_ξ ./= sum(m_ξ, dims=[2,3])

    #check equivalency of xi calc methods
    @test isapprox(reshape(ξ,Tmaxplus1,K,K)[1:21,:,:],m_ξ)

    γ = exp.(log_γ)
    k_γ = sum(γ, dims=3); k_γ[nan_mask[:,:,:]] .= Inf #prevent NaNs in dummy renorm
    γ ./= k_γ

    m_γ = exp.((m_log_α .+ m_log_β) .- normalizer)

    #check equivalency of gamma calculations
    @test isapprox(reshape(γ,Tmaxplus1,K)[1:21,:,:],m_γ)

    # "Testing initial and transition matrix calcs..."

    m_new_A = sum(m_ξ[1:end-1,:,:], dims=1)[1,:,:]
    m_new_A ./= sum(m_new_A, dims=2)

    new_A = zeros(K,K)
    for i=1:K, j=1:K
        ∑otξ_vec = zeros(O)
        ∑otγ_vec = zeros(O)
        for o in 1:O
            ∑otξ_vec[o] = sum(ξ[1:obs_lengths[o]-1,o,i,j])
            ∑otγ_vec[o] = sum(γ[1:obs_lengths[o]-1,o,i])
        end
        new_A[i,j] = sum(∑otξ_vec) / sum(∑otγ_vec)
    end
    new_A ./= sum(new_A, dims=[2]) #dummy renorm

    #check equivalency of transition matrix calculations
    @test isapprox(m_new_A,new_A)

    m_new_a = exp.((m_log_α[1,:] + m_log_β[1,:]) .- normalizer)
    m_new_a ./= sum(m_new_a)

    new_a = (sum(γ[1,:,:], dims=1)./sum(sum(γ[1,:,:], dims=1)))[1,:]
    new_a ./= sum(new_a) #dummy renorm

    @test isapprox(m_new_a,new_a)

    # "Testing distribution calcs..."

    F=Univariate
    m_D = Distribution{F}[]
    for (i, d) in enumerate(hmm.B)
        # Super hacky...
        # https://github.com/JuliaStats/Distributions.jl/issues/809
        push!(m_D, fit_mle(eval(typeof(d).name.name), permutedims(obs[1:21]), m_γ[:,i]))
    end

    obs_mask = .!nan_mask
    obs_collection = obs[obs_mask[:,:]]

    B = Distribution{F}[]
    @inbounds for (i, d) in enumerate(hmm.B)
        # Super hacky...
        # https://github.com/JuliaStats/Distributions.jl/issues/809
        γ_d = γ[:,:,i]
        push!(B, fit_mle(eval(typeof(d).name.name), obs_collection, γ_d[obs_mask[:,:]]))
        #slowest call by orders of magnitude
    end

    #test maximization equivalency
    for (d, dist) in enumerate(B)
        @test isapprox(dist.support, m_D[d].support)
        @test isapprox(dist.p, m_D[d].p)
    end

    # "Testing mle_step..."

    #verify that above methods independently produce equivalent output, and that this is true of multiple identical obs, but not true of different obs sets
    mouchet_hmm = mouchet_mle_step(hmm, obs[1:21])

    new_hmm = BioBackgroundModels.bw_step(hmm, obs, obs_lengths)

    dblobs = zeros(UInt8, 22,2)
    dblobs[1:21,1] = [4,3,3,2,3,2,1,1,2,3,3,3,4,4,2,3,2,3,4,3,2]
    dblobs[1:21,2] = [4,3,3,2,3,2,1,1,2,3,3,3,4,4,2,3,2,3,4,3,2]
    dblobs_lengths=[21,21]
    dbl_hmm = BioBackgroundModels.bw_step(hmm, dblobs, dblobs_lengths)

    otherobs = dblobs
    otherobs[1:21,2] = [1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1]
    other_hmm = BioBackgroundModels.bw_step(hmm, otherobs, dblobs_lengths)

    for n in fieldnames(typeof(new_hmm[1]))
        if n == :B
            for (d, dist) in enumerate(new_hmm[1].B)
                @test dist.support == mouchet_hmm[1].B[d].support
                @test isapprox(dist.p,mouchet_hmm[1].B[d].p)
                @test dist.support == dbl_hmm[1].B[d].support
                @test isapprox(dist.p,dbl_hmm[1].B[d].p)
                @test dist.support == other_hmm[1].B[d].support
                @test !isapprox(dist.p, other_hmm[1].B[d].p)
            end
        elseif n == :partition
            @test getfield(new_hmm[1],n) == getfield(mouchet_hmm[1],n)
            @test getfield(new_hmm[1],n) == getfield(dbl_hmm[1],n)
            @test getfield(new_hmm[1],n) == getfield(other_hmm[1],n)
        else
            @test isapprox(getfield(new_hmm[1],n), getfield(mouchet_hmm[1],n))
            @test isapprox(getfield(new_hmm[1],n), getfield(dbl_hmm[1],n))
            @test !isapprox(getfield(new_hmm[1],n), getfield(other_hmm[1],n))
        end
    end

    @test new_hmm[2] == mouchet_hmm[2] != dbl_hmm[2] != other_hmm[2]

    # "Testing fit_mle!..."

    # "test fit_mle! function"
    input_hmms= RemoteChannel(()->Channel{Tuple}(1))
    output_hmms = RemoteChannel(()->Channel{Tuple}(30))
    chainid=Chain_ID("Test",2,0,1)
    put!(input_hmms, (chainid, 2, hmm, 0.0, transpose(obs)))
    BioBackgroundModels.EM_converge!(input_hmms, output_hmms, 1; EM_func=BioBackgroundModels.bw_step, max_iterates=4, verbose=true)
    wait(output_hmms)
    workerid, jobid, iterate, hmm3, log_p, epsilon, converged = take!(output_hmms)
    @test jobid == chainid
    @test iterate == 3
    @test assert_hmm(hmm3.a, hmm3.A, hmm3.B)
    @test size(hmm3) == size(hmm) == (2,4)
    @test log_p < 1
    @test log_p == obs_lh_given_hmm(transpose(obs),hmm, linear=false)
    @test converged == false
    wait(output_hmms)
    workerid, jobid, iterate, hmm4, log_p, epsilon, converged = take!(output_hmms)
    @test jobid == chainid
    @test iterate == 4
    @test assert_hmm(hmm4.a, hmm4.A, hmm4.B)
    @test size(hmm4) == size(hmm) == (2,4)
    @test log_p < 1
    @test log_p == obs_lh_given_hmm(transpose(obs),hmm3,linear=false)
    @test converged == false

    # "Test convergence.."
    obs=zeros(UInt8, 101, 2)
    for i in 1:size(obs)[2]
        obs[1:100,i]=rand(1:4,100)
    end
    input_hmms= RemoteChannel(()->Channel{Tuple}(1))
    output_hmms = RemoteChannel(()->Channel{Tuple}(30))
    put!(input_hmms, (chainid, 2, hmm, 0.0, transpose(obs)))
    BioBackgroundModels.EM_converge!(input_hmms, output_hmms, 1; EM_func=BioBackgroundModels.bw_step, delta_thresh=.05, max_iterates=100, verbose=true)
    wait(output_hmms)
    workerid, jobid, iterate, hmm4, log_p, epsilon, converged = take!(output_hmms)
    while isready(output_hmms)
        workerid, jobid, iterate, hmm4, log_p, epsilon, converged = take!(output_hmms)
    end
    @test converged==1

    #test t_add_categorical_counts DimensionMismatch
    @test_throws DimensionMismatch t_add_categorical_counts!(zeros(4), zeros(UInt8,2,3), zeros(2,2))

end

@testset "EM/churbanov.jl Baum-Welch equivalency and functions" begin
    # "Setting up for MLE function tests.."
    A = fill((1/6),6,6)
    B = [Categorical(ones(4)/4), Categorical([.7,.05,.15,.1]),Categorical([.15,.35,.4,.1]), Categorical([.6,.15,.15,.1]),Categorical([.1,.4,.4,.1]), Categorical([.2,.2,.3,.3])]
    hmm = BHMM(A, B)
    log_A = log.(hmm.A)

    obs = zeros(Int64,1,250)
    obs[1:249] = rand(1:4,249)
    obs_lengths=[249]
    # "Testing mle_step..."

    #verify that above methods independently produce equivalent output, and that this is true of multiple identical obs, but not true of different obs sets
    mouchet_hmm = mouchet_mle_step(hmm, obs[1:249])

    new_hmm = linear_step(hmm, obs, obs_lengths)

    ms_sng = BioBackgroundModels.bw_step(hmm, Array(transpose(obs)), obs_lengths)

    dblobs = zeros(Int64, 2,250)
    dblobs[1,1:249] = obs[1:249]
    dblobs[2,1:249] = obs[1:249]
    dblobs_lengths=[249,249]
    dbl_hmm = linear_step(hmm, dblobs, dblobs_lengths)

    ms_dbl =BioBackgroundModels.bw_step(hmm, Array(transpose(dblobs)),dblobs_lengths)

    otherobs = deepcopy(dblobs)
    otherobs[2,1:249] = rand(1:4,249)
    if otherobs[1,1]==otherobs[2,1]
        otherobs[1,1]==1&&(otherobs[2,1]=2)
        otherobs[1,1]==2&&(otherobs[2,1]=3)
        otherobs[1,1]==3&&(otherobs[2,1]=4)
        otherobs[1,1]==4&&(otherobs[2,1]=1)
    end

    other_hmm = linear_step(hmm, otherobs, dblobs_lengths)

    ms_hmm = BioBackgroundModels.bw_step(hmm, Array(transpose(otherobs)),dblobs_lengths)

    for n in fieldnames(typeof(new_hmm[1]))
        if n == :B
            for (d, dist) in enumerate(new_hmm[1].B)
            @test dist.support==mouchet_hmm[1].B[d].support
            @test isapprox(dist.p,mouchet_hmm[1].B[d].p)
            @test dist.support==ms_sng[1].B[d].support
            @test isapprox(dist.p,ms_sng[1].B[d].p)
            @test dist.support==dbl_hmm[1].B[d].support
            @test isapprox(dist.p,dbl_hmm[1].B[d].p)
            @test dist.support==ms_dbl[1].B[d].support
            @test isapprox(dist.p,ms_dbl[1].B[d].p)
            @test dist.support==other_hmm[1].B[d].support
            @test ms_hmm[1].B[d].support==other_hmm[1].B[d].support
            @test !isapprox(dist.p, other_hmm[1].B[d].p)
            @test isapprox(ms_hmm[1].B[d].p, other_hmm[1].B[d].p)

            end
        elseif n == :partition
            @test getfield(new_hmm[1],n) == getfield(mouchet_hmm[1],n)
            @test getfield(new_hmm[1],n) == getfield(dbl_hmm[1],n)
            @test getfield(new_hmm[1],n) == getfield(other_hmm[1],n)
        else
            @test isapprox(getfield(new_hmm[1],n), getfield(mouchet_hmm[1],n))
            @test isapprox(getfield(new_hmm[1],n), getfield(ms_sng[1],n))

            @test isapprox(getfield(new_hmm[1],n), getfield(dbl_hmm[1],n))
            @test isapprox(getfield(new_hmm[1],n), getfield(ms_dbl[1],n))

            @test !isapprox(getfield(new_hmm[1],n), getfield(other_hmm[1],n))
            @test isapprox(getfield(ms_hmm[1],n), getfield(other_hmm[1],n))

        end
    end

    @test ms_sng[2] == new_hmm[2] == mouchet_hmm[2] != dbl_hmm[2] == ms_dbl[2] != other_hmm[2] == ms_hmm[2]

    # "Testing fit_mle!..."

    #test fit_mle! function
    input_hmms= RemoteChannel(()->Channel{Tuple}(1))
    output_hmms = RemoteChannel(()->Channel{Tuple}(30))
    chainid=Chain_ID("Test",6,0,1)
    put!(input_hmms, (chainid, 2, hmm, 0.0, obs))
    BioBackgroundModels.EM_converge!(input_hmms, output_hmms, 1; max_iterates=4, verbose=true)
    wait(output_hmms)
    workerid, jobid, iterate, hmm3, log_p, delta, converged = take!(output_hmms)
    @test jobid == chainid
    @test iterate == 3
    @test BioBackgroundModels.assert_hmm(hmm3.a, hmm3.A, hmm3.B)
    @test size(hmm3) == size(hmm) == (6,4)
    @test log_p < 1
    @test log_p == linear_likelihood(obs, hmm)
    @test converged == false
    wait(output_hmms)
    workerid, jobid, iterate, hmm4, log_p, delta, converged = take!(output_hmms)
    @test jobid == chainid
    @test iterate == 4
    @test BioBackgroundModels.assert_hmm(hmm4.a, hmm4.A, hmm4.B)
    @test size(hmm4) == size(hmm) == (6,4)
    @test log_p < 1
    @test log_p == linear_likelihood(obs, hmm3)
    @test converged == false

    # "Test convergence.."
    obs=zeros(Int64, 4, 1001)
    for o in 1:size(obs)[1]
        obsl=rand(100:1000)
        obs[o,1:obsl]=rand(1:4,obsl)
    end
    input_hmms= RemoteChannel(()->Channel{Tuple}(1))
    output_hmms = RemoteChannel(()->Channel{Tuple}(Inf))
    put!(input_hmms, (chainid, 2, hmm, 0.0, obs))
    BioBackgroundModels.EM_converge!(input_hmms, output_hmms, 1; delta_thresh=.05, max_iterates=100, verbose=true)
    wait(output_hmms)
    workerid, jobid, iterate, hmm4, log_p, delta, converged = take!(output_hmms)
    while isready(output_hmms)
        workerid, jobid, iterate, hmm4, log_p, delta, converged = take!(output_hmms)
    end
    @test converged==1
end

@testset "API/EM_master.jl and reports.jl API tests" begin
    #CONSTANTS FOR BHMM LEARNING
    replicates = 4 #repeat optimisation from this many seperately initialised samples from the prior
    Ks = [1,2,4,6] #mosaic class #s to test
    order_nos = [0,1,2] #DNA kmer order #s to test
    sample_set_length=100
    min_sample_window=5
    max_sample_window=25
    perigenic_pad=250

    wkpool=addprocs(2)
    @everywhere using BioBackgroundModels

    channels = setup_sample_jobs(genome, index, gff, sample_set_length, min_sample_window, max_sample_window, perigenic_pad; deterministic=true)
    sample_record_dfs=execute_sample_jobs(channels, wkpool)
    training_sets, test_sets = split_obs_sets(sample_record_dfs)

    rmprocs(wkpool)

    job_ids=Vector{Chain_ID}()
    for (obs_id, obs) in training_sets, K in Ks, order in order_nos, rep in 1:replicates
        push!(job_ids, Chain_ID(obs_id, K, order, rep))
    end

    no_input_hmms, chains, input_hmms, output_hmms = setup_EM_jobs!(job_ids, training_sets)
    @test no_input_hmms == replicates*length(Ks)*length(order_nos)*length(training_sets)

    while isready(input_hmms)
        jobid, start_iterate, hmm, last_norm, observations = take!(input_hmms)
        @test last_norm == linear_likelihood(observations, hmm)
        obs_lengths = [findfirst(iszero,observations[o,:])-1 for o in 1:size(observations)[1]]
        #make sure input HMMs are valid and try to mle_step them and ensure their 1-step children are valid
        @test assert_hmm(hmm.a, hmm.A, hmm.B)
        new_hmm, prob = linear_step(hmm,observations,obs_lengths)
        @test assert_hmm(hmm.a, hmm.A, hmm.B)
        @test prob < 0
    end

    genome_reader = open(FASTA.Reader, genome, index=index)
    seq=LongSequence{DNAAlphabet{2}}(FASTA.sequence(genome_reader["CM002885.2.1"]))
    seqdict = Dict("test"=>[seq for i in 1:3])
    seqdict["blacklist"]=seqdict["test"]

    job_ids=Vector{Chain_ID}()
    push!(job_ids, Chain_ID("blacklist", 1, 0, 1))
    push!(job_ids, Chain_ID("blacklist", 1, 0, 2))
    Ks=[1,2]; order_nos=[0,1]
    for K in Ks, order in order_nos, replicate in 1:2
        push!(job_ids, Chain_ID("test", K, order, replicate))
    end

    wkpool=addprocs(2, topology=:master_worker)
    @everywhere using BioBackgroundModels

    load_dict=Dict{Int64,LoadConfig}()
    for (n,wk) in enumerate(wkpool)
        n==1 && (load_dict[wk]=LoadConfig(1:2,0:1,blacklist=[Chain_ID("blacklist",1,0,1)]))
        n==2 && (load_dict[wk]=LoadConfig(1:2,0:1))
    end
    #test work resumption, load configs

    em_jobset=setup_EM_jobs!(job_ids, seqdict, delta_thresh=10000.)
    execute_EM_jobs!(wkpool, em_jobset..., "testchains", delta_thresh=10000., verbose=true, load_dict=load_dict)

    for (chain_id,chain) in em_jobset[2]
        @test chain[end].converged == true
        @test chain[end].delta <= 10000
        @test typeof(chain[end].hmm)==BHMM{Float64}
        @test 1 <= chain[end].iterate <= 5000
    end

    wkpool=addprocs(2, topology=:master_worker)
    @everywhere using BioBackgroundModels

    em_jobset=setup_EM_jobs!(job_ids, seqdict, chains=deserialize("testchains"), delta_thresh=.1)
    execute_EM_jobs!(wkpool, em_jobset..., "testchains", delta_thresh=.1, verbose=true)

    for (chain_id,chain) in em_jobset[2]
        @test chain[end].converged == true
        @test chain[end].delta <= .1
        @test typeof(chain[end].hmm)==BHMM{Float64}
        @test 1 <= chain[end].iterate <= 5000
    end

    #test for already converged warning
    wkpool=addprocs(2, topology=:master_worker)
    @everywhere using BioBackgroundModels
    @test_throws ArgumentError execute_EM_jobs!(wkpool, em_jobset..., "testchains", delta_thresh=.01, verbose=true)

    rm("testchains")

    @test_throws ArgumentError generate_reports(Dict{Chain_ID,Vector{EM_step}}(),seqdict)
    @test_throws ArgumentError generate_reports(em_jobset[2],Dict{String,Vector{LongSequence{DNAAlphabet{2}}}}())

    @test_throws ArgumentError BioBackgroundModels.report_chains(em_jobset[2],Dict{String,Vector{LongSequence{DNAAlphabet{2}}}}())

    report_folders=generate_reports(em_jobset[2],seqdict)
    folder=report_folders["test"]
    @test "test"==folder.partition_id
    show(folder.partition_report)
    show(folder.replicate_report)
    for id in folder.partition_report.best_repset
        show(folder.chain_reports[id])
    end
end

rm(genome)
rm(index)
rm(gff)
rm(posfasta)


text:
function render_template( filename, args = dict{String, String}())
  f = open(filename)
# to be finished soon sorry :(
end


text:
###
### Initialization
###

# Initialize particle swarm state

function initialize(rng::AbstractRNG, r::Union{ParamRange, Tuple{ParamRange, Any}}, n::Int)
    return initialize(rng, [r], n)
end

function initialize(rng::AbstractRNG, rs::AbstractVector, n::Int)
    # `length` is 1 for `NumericRange` and the number of categories for `NominalRange`
    ranges, parameters, lengths, Xᵢ = zip(_initialize.(rng, rs, n)...)
    indices = _to_indices(lengths)
    X = hcat(Xᵢ...)
    V = zero(X)
    pbest_X = copy(X)
    gbest_X = copy(X)
    pbest = fill(eltype(X)(Inf), n)
    gbest = similar(pbest)
    return ParticleSwarmState(
        ranges, parameters, indices, X, V, pbest_X, gbest_X, pbest, gbest
    )
end

# Unpack tuple of range and distribution

function _initialize(rng, t::Tuple{ParamRange, Any}, n)
    return _initialize(rng, t[1], t[2], n)
end

# Initialize parameters with default distributions

function _initialize(rng, r::NominalRange{T, N}, n) where {T, N}
    d = Dirichlet(ones(N))
    return _initialize(rng, r, d, n)
end

function _initialize(rng, r::NumericRange, n)
    D = _initializer(MLJTuning.boundedness(r))
    return _initialize(rng, r, D, n)
end

_initializer(::Type{MLJBase.Bounded}) = Uniform

_initializer(::Type{MLJTuning.PositiveUnbounded}) = Gamma

_initializer(::Type{MLJTuning.Other}) = Normal

# Fit distributions and initialize parameters

function _initialize(rng, r::ParamRange, D::Type{<:Distribution}, n)
    throw(ArgumentError("$D distribution is unsupported for $(typeof(r))."))
end

function _initialize(rng, r::NumericRange, D::Type{<:UnivariateDistribution}, n)
    d = Distributions.fit(D, r)
    return _initialize(rng, r, d, n)
end

# Initialize parameters with fitted/provided distributions

function _initialize(rng, r::ParamRange, d::Distribution, n)
    throw(ArgumentError("$(typeof(d)) distribution is unsupported for $(typeof(r))."))
end

function _initialize(rng, r::NominalRange{T, N}, d::Dirichlet, n) where {T, N}
    N != d.alpha0 &&
        throw(ArgumentError("Provided distribution's number of categories don't match $r."))
    p = Vector{T}(undef, n)
    X = rand(rng, d, n)'
    return r, p, N, X
end

function _initialize(rng, r::NumericRange{T}, d::UnivariateDistribution, n) where {T}
    p = Vector{T}(undef, n)
    X = rand(rng, d, n)
    return r, p, 1, X
end

# Helper function to get ranges' corresponding indices
# E.g., `_to_indices((1, 2, 1, 3))` returns `(1, 2:3, 4, 5:7)`

function _to_indices(lengths)
    curr = 1
    return map(lengths) do length
        start = curr
        stop = start + length - 1
        curr = stop + 1
        start == stop ? stop : (start:stop)
    end
end

###
### Retrieval
###

# For updating `state.parameters`, the model hyperparameters to be returned, from their
# internal representation `state.X`

function retrieve!(rng::AbstractRNG, state::ParticleSwarmState)
    ranges, params, indices, X = state.ranges, state.parameters, state.indices, state.X
    for (r, p, i) in zip(ranges, params, indices)
        _retrieve!(rng, p, r, view(X, :, i))
    end
    return state
end

function _retrieve!(rng, p, r::NominalRange, X)
    return p .= getindex.(
        Ref(r.values),
        rand.(rng, Categorical.(X[i,:] for i in axes(X, 1)))
    )
end

function _retrieve!(rng, p, r::NumericRange{T}, X) where {T<:Integer}
    return @. p = round(T, _transform(r.scale, X))
end

function _retrieve!(rng, p, r::NumericRange{T}, X) where {T<:Real}
    return @. p = _transform(r.scale, X)
end

_transform(::Symbol, X) = X

_transform(scale, X) = scale(X)


text:
module Utils

using LinearAlgebra
export ref_to_pointer, indexOfField, memmove, sizeof_obj

include("MathUtils.jl")
include("StringUtils.jl")

memmove(dest, src, n) = ccall(:memmove, Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Int), dest, src, n)

function sizeof_obj(x)
    total = 0;
    fieldNames = fieldnames(typeof(x));
    if length(fieldNames) == 0
        return sizeof(x);
    else
        for fieldName in fieldNames
            total += sizeof_obj(getfield(x, fieldName));
        end
        return total;
    end
end

function indexOfField(type, field)
    names = fieldnames(type)
    count = length(names)
    for i = 1 : count
        if names[i] == field
            return i
        end
    end
    return 0
end

function ref_to_pointer(type::Type, ref)
    ref == nothing ? C_NULL : Base.unsafe_convert(Ptr{type}, ref)
end

end #module


text:
# # Pinning Julia Threads
#
# In the following, we demonstrate how to use LIKWID.jl to pin
# Julia threads to specific cores.

#md # !!! note
#md #     Note that Julia implements task-based multithreading where `N`
#md #     tasks get mapped onto `M` OS threads (M:N hybrid threading).
#md #     We will pin the Julia (p)threads and not the tasks.
#md #     Depending on how the latter are started/configured, tasks may
#md #     migrate between Julia threads!

# First, we load the packages and standard libraries that we'll use.
using LIKWID
using Base.Threads: nthreads, @threads
using Test
using Random

# Note that I have started julia with multiple threads (`julia -t N`), concretely
NT = nthreads()

# Let's find out on which cores are running our Julia threads.
# We will use [`LIKWID.get_processor_id`](@ref) in combination with
# `@threads :static for`, which guarantees that the associated tasks
# are themselves pinned to Julia threads.
coreids = zeros(Int, NT)
@threads :static for i in 1:NT
    coreids[i] = LIKWID.get_processor_id()
end
println("Cores: ", coreids)

# Since querying all core ids is a common operation, we provide [`LIKWID.get_processor_ids`](@ref) which returns all core ids right away.

# To pin a thread to a specific core, there is [`LIKWID.pinthread`](@ref). Using `@threads :static for` like above, we can, for example, pin the `NT` Julia threads to the first `NT` cores:
## core numbering starts with zero
cores_firstNT = 0:NT-1
@threads :static for i in 1:NT
    LIKWID.pinthread(cores_firstNT[i])
end
@test LIKWID.get_processor_ids() == cores_firstNT

# To avoid the explicit for-loop, we can directly use [`LIKWID.pinthreads`](@ref) to pin all Julia threads.

# Let's try a less trivial shuffled mapping
cores_firstNT_shuffeled = shuffle(cores_firstNT)
LIKWID.pinthreads(cores_firstNT_shuffeled)
@test LIKWID.get_processor_ids() == cores_firstNT_shuffeled

# Finally, let's choose the cores at random
topo = LIKWID.get_cpu_topology()
ncores = topo.numSockets * topo.numCoresPerSocket
cores_rand = shuffle(0:ncores-1)[1:NT]
LIKWID.pinthreads(cores_rand)
@test LIKWID.get_processor_ids() == cores_rand


text:
# bitmap

const RES_PATH = isassigned(ARGS, 1) ? ARGS[1] : "./res"

struct Bitmap
    width::UInt
    height::UInt
    num_channels::UInt8
    max_color::UInt16
    colors::Vector{UInt16}
end

function read_pbm_file(name::String, path::String, verbose::Bool)::Bitmap
    fullpath = abspath(joinpath(path, "$name.pbm"))
    if verbose
        println("path: $fullpath")
    end

    num_channels = max_color = width = height = colors = nothing
    open(fullpath, "r") do file_io
        # header
        magic_number = file_io |> readline |> strip
        @assert startswith(magic_number, "P") "first line must start with 'P'. - $magic_number"
        if magic_number == "P1"
            # 0-1 (white & black)
            num_channels = 1
            max_color = 1
        else
            error("not supported magic number. - P$magic_number")
        end

        # color stream
        width, height = parse.([UInt], split(readline(file_io), " "))
        colors = zeros(Int16, width * height * num_channels)
        color_index = 1
        for line in eachline(file_io)
            for color in parse.([UInt16], split(line, " "))
                colors[color_index] = color
                color_index += 1
            end
        end
    end

    bitmap = Bitmap(width, height, num_channels, max_color, colors)
    if verbose
        println("size: ($width, $height), num of channels: $num_channels, max color: $max_color")
        println("length: ", length(colors))
        println("content:")
        show(bitmap)
    end
    return bitmap
end

function float_colors(bitmap::Bitmap)::Vector{Float32}
    return bitmap.colors ./ bitmap.max_color
end

function pick_color(bitmap::Bitmap, x::UInt, y::UInt)::Vector{UInt16}
    offset = (y-1)*bitmap.width + x*bitmap.num_channels
    limit = offset + bitmap.num_channels - 1
    return bitmap.colors[offset:limit]
end

function pick_color_code(bitmap::Bitmap, x::UInt, y::UInt)::String
    pad = 1
    if bitmap.max_color >= 0xffffff
        pad = 4
    elseif bitmap.max_color >= 0xff
        pad = 2
    end
    return join([string(color, base=16, pad=pad) for color in pick_color(bitmap, x, y)], "")
end

function show(bitmap::Bitmap)
    for y = 1:bitmap.height
        line = join([pick_color_code(bitmap, x, y) for x = 1:bitmap.width], " ")
        println("$line")
    end
end

Bitmap(name::String; path::String=RES_PATH, verbose::Bool=false) = read_pbm_file(name, path, verbose)

# bias

const TRAINING_ALPHA = 0.9f0
const TRAINING_ETA = 0.1f0

mutable struct Bias
    value::Float32
    momentum::Float32

    Bias() = new(zero(Float32), zero(Float32))
end

function random_biases(num::Int)::Vector{Bias}
    biases = fill(Bias(), num)
    for i = 1:num
        biases[i].value = rand() - 0.5  # -0.5 ~ 0.5
        biases[i].momentum = zero(Float32)
    end
    return biases
end

function update!(bias::Bias, delta::Float32)
    delta_bias = TRAINING_ALPHA*bias.momentum + TRAINING_ETA*delta
    bias.value += delta_bias
    bias.momentum = delta_bias
end

# weight

struct Weight
    values::Vector{Float32}
    momentums::Vector{Float32}

    Weight(num::Int) = new(Vector{Float32}(undef, num), zeros(Float32, num))
end

function random_weight(num::Int)::Weight
    weight = Weight(num)
    for i = 1:num
        weight.values[i] = rand() - 0.5 # -0.5 ~ 0.5
    end
    return weight
end

function update!(weight::Weight, delta::Float32, prev_outputs::Vector{Float32})
    for i = 1:length(weight.values)
        delta_weight = TRAINING_ALPHA*weight.momentums[i] + TRAINING_ETA*prev_outputs[i]*delta
        weight.values[i] += delta_weight
        weight.momentums[i] = delta_weight
    end
end

# layer

abstract type Layer end

sigmoid(value::Float32)::Float32 = one(Float32) / (one(Float32) + exp(-value))

function input!(layer::Layer, prev_layer::Layer)
    num_prev_outputs = length(prev_layer.outputs)

    for i = 1:length(layer.outputs)
        if !isassigned(layer.weights, i)
            layer.weights[i] = random_weight(num_prev_outputs)
        end

        sum_nets = zero(Float32)
        for j = 1:num_prev_outputs
            sum_nets += layer.weights[i].values[j] * prev_layer.outputs[j]
        end
        layer.outputs[i] = sigmoid(sum_nets + layer.biases[i].value)
    end
end

struct InputLayer <: Layer
    outputs::Vector{Float32}

    InputLayer(num::Int) = new(zeros(Float32, num))
end

function input!(layer::InputLayer, bitmaps::Bitmap)
    copyto!(layer.outputs, float_colors(bitmaps))
end

struct HiddenLayer <: Layer
    biases::Vector{Bias}
    weights::Vector{Weight}
    deltas::Vector{Float32}
    outputs::Vector{Float32}

    HiddenLayer(num::Int) = new(
        random_biases(num), Vector{Weight}(undef, num), zeros(Float32, num), zeros(Float32, num))
end

function update!(layer::HiddenLayer, prev_layer::Layer, next_layer::Layer)
    copyto!(layer.deltas, next_layer.deltas)

    for i = 1:length(layer.outputs)
        layer.deltas[i] *= layer.outputs[i] * (one(Float32) - layer.outputs[i]) # delta=o*(1-o)*SUM(u_delta*Woh)
        delta = layer.deltas[i]

        update!(layer.biases[i], delta)
        update!(layer.weights[i], delta, prev_layer.outputs)
    end
end

mutable struct OutputLayer <: Layer
    biases::Vector{Bias}
    weights::Vector{Weight}
    deltas::Vector{Float32}
    outputs::Vector{Float32}

    OutputLayer(num::Int) = new(
        random_biases(num), Vector{Weight}(undef, num), zeros(Float32, num), zeros(Float32, num))
end

function diff_target(layer::OutputLayer, index::Int, is_correct::Bool)::Float32
    return (is_correct ? one(Float32) : zero(Float32)) - layer.outputs[index]
end

function update!(layer::OutputLayer, prev_layer::Layer, target_index::Int)::Float64
    num_prev_outputs = length(prev_layer.outputs)
    layer.deltas = zeros(Float32, num_prev_outputs)

    num_outputs = length(layer.outputs)
    sum_errors::Float64 = zero(Float64)
    for i = 1:num_outputs
        diff = diff_target(layer, i, i == target_index)
        sum_errors += 0.5 * diff^2.0                                        # E=1/2(d-o)^2
        delta = diff * layer.outputs[i] * (one(Float32) - layer.outputs[i]) # delta=(d-o)o(1-o)
        for j = 1:num_prev_outputs
            layer.deltas[j] += delta * layer.weights[i].values[j]
        end

        update!(layer.biases[i], delta)
        update!(layer.weights[i], delta, prev_layer.outputs)
    end
    return sum_errors / num_outputs
end

# recognizer

abstract type NeuralNetwork end

struct Recognizer <: NeuralNetwork
    input_layer::InputLayer
    hidden_layer::HiddenLayer
    output_layer::OutputLayer

    Recognizer(input_size::Int, output_size::Int) = new(
        InputLayer(input_size), HiddenLayer(output_size * 3), OutputLayer(output_size))
end

function input!(recognizer::Recognizer, bitmap::Bitmap)
    input!(recognizer.input_layer, bitmap)
    input!(recognizer.hidden_layer, recognizer.input_layer)
    input!(recognizer.output_layer, recognizer.hidden_layer)
end

function update!(recognizer::Recognizer, target_index::Int)::Float64    # reverse propagation
    err = update!(recognizer.output_layer, recognizer.hidden_layer, target_index)
    update!(recognizer.hidden_layer, recognizer.input_layer, recognizer.output_layer)
    return err
end

function train!(nn::NeuralNetwork, bitmaps::Vector{Bitmap})::Float64
    num_bitmaps = length(bitmaps)
    sum_errors = zero(Float64)
    for i = 1:num_bitmaps
        input!(nn, bitmaps[i])
        sum_errors += update!(nn, i)
    end
    return sum_errors / num_bitmaps
end

function classify(nn::NeuralNetwork, bitmap::Bitmap)::Int
    input!(nn, bitmap)
    return argmax(nn.output_layer.outputs)
end

# main

const BITMAP_SIZE = 100         # width * height
const NUM_CHARS = 'Z' - 'A' + 1 # number of alphabet
const TRAINING_ERROR_THRESHOLD = 0.0001
const TRAINING_COUNT_LIMIT = 10000

to_char(index::Int)::Char = index + 'A' - 1

println("[[ file loading ]]")
bitmaps = Vector{Bitmap}([ch |> string |> Bitmap for ch in "ABCDEFGHIJKLMNOPQRSTUVWXYZ"])
recognizer = Recognizer(BITMAP_SIZE, NUM_CHARS)
println("")

println("[[ pattern training ]]")
for cnt = 1:TRAINING_COUNT_LIMIT
    print("training pattern $cnt...")
    err = train!(recognizer, bitmaps)
    println("done. (error: $err)")
    if err < TRAINING_ERROR_THRESHOLD
        break
    end
end
println("")

println("[[ classification testing ]]")
for i = 1:NUM_CHARS
    ch = to_char(i)
    print("test character '$ch'...")
    classified = classify(recognizer, bitmaps[i]) |> to_char
    if classified == ch
        println("correct!")
    else
        println("incorrect! (classified $classified)")
    end
end
println("")


text:
using FractionalSystems, ControlSystems
using Test

@testset "Test Lyapunov exponents" begin
    function RF(du, u, t)
        du[1] = u[2]*(u[3]-1+u[1]*u[1])+0.1*u[1]
        du[2] = u[1]*(3*u[3]+1-u[1]*u[1])+0.1*u[2]
        du[3] = -2*u[3]*(0.98+u[1]*u[2])
    end

    (LE, tspan)=FOLyapunov(RF, 0.999, 0, 0.02, 300, [0.1; 0.1; 0.1], 0.005, 1000)

    # Test the latest computed(final) Lyapunov exponents
    @test isapprox(LE[end-2:end], [0.06111650166568285
    0.0038981396237095034
   -1.8324646820425692]; atol=1e-3)
end


text:
module AdaptiveSampling
#using Plots
include("anchor_point.jl")

end # module


text:
using NBInclude
nbinclude(joinpath(@__DIR__, "..", "demo.ipynb"))


text:
# Functions for exporting molecular dynamics simulation results to visualization tools

"""
    write_nbs_animation(result::NBSResult, filename::String)

Animate an `NBSResult` and store the result in a .gif file.
"""
function write_nbs_animation(result::NBSResult, filename::String)
    animate(result.result, filename)
end

"""
    write_ase_trajectory(result::MolecularDynamicsResult, element::DFTK.Element, lattice, filename::String)

Write the trajectory of a `MolecularDynamicsResult` to a .traj file.

The file can be visualized by running `ase gui <filename>` on the command line.
"""
# TODO: support multiple species using atom data
function write_ase_trajectory(result::MolecularDynamicsResult, element::DFTK.Element, lattice, filename::String)
    ase = pyimport_e("ase")
    if ispynull(ase)
        @warn "ASE is not installed, skipping write_trajectory"
    else
        filename = abspath(expanduser(filename))
        if isfile(filename)
            @info "overwriting trajectory at" filename
            rm(filename)
        else
            @info "writing trajectory to" filename
        end
        traj = pyimport("ase.io.trajectory").Trajectory(filename, "a")
        for t ∈ 1:length(get_time_range(result))
            traj.write(ase_atoms(lattice, dftk_atoms(get_system(result, t), element)))
        end
        traj.close()
    end
end


text:
# Since Kronecker does not implement `TensorCore.:⊗` but instead exports its own function
# `Kronecker.:⊗`, only the module is imported and Kronecker.:⊗ and Kronecker.kronecker are
# called explicitly.
using .Kronecker: Kronecker

export kernelkronmat
export kronecker_kernelmatrix

function kernelkronmat(κ::Kernel, X::AbstractVector, dims::Int)
    @assert iskroncompatible(κ) "The chosen kernel is not compatible for kroenecker matrices (see [`iskroncompatible`](@ref))"
    k = kernelmatrix(κ, X)
    return Kronecker.kronecker(k, dims)
end

function kernelkronmat(
    κ::Kernel, X::AbstractVector{<:AbstractVector}; obsdim::Int=defaultobs
)
    @assert iskroncompatible(κ) "The chosen kernel is not compatible for Kronecker matrices"
    Ks = kernelmatrix.(κ, X)
    return K = reduce(Kronecker.:⊗, Ks)
end

"""
    To be compatible with kroenecker constructions the kernel must satisfy
    the property : for x,x' ∈ ℜᴰ
    k(x,x') = ∏ᵢᴰ k(xᵢ,x'ᵢ)
"""
@inline iskroncompatible(κ::Kernel) = false # Default return for kernels

function _kernelmatrix_kroneckerjl_helper(
    ::Type{<:MOInputIsotopicByFeatures}, Kfeatures, Koutputs
)
    return Kronecker.kronecker(Kfeatures, Koutputs)
end

function _kernelmatrix_kroneckerjl_helper(
    ::Type{<:MOInputIsotopicByOutputs}, Kfeatures, Koutputs
)
    return Kronecker.kronecker(Koutputs, Kfeatures)
end

"""
    kronecker_kernelmatrix(
        k::Union{IndependentMOKernel,IntrinsicCoregionMOKernel}, x::MOI, y::MOI
    ) where {MOI<:IsotopicMOInputsUnion}

Requires Kronecker.jl: Computes the `kernelmatrix` for the `IndependentMOKernel` and the
`IntrinsicCoregionMOKernel`, but returns a lazy kronecker product. This object can be very
efficiently inverted or decomposed. See also [`kernelmatrix`](@ref).
"""
function kronecker_kernelmatrix(
    k::Union{IndependentMOKernel,IntrinsicCoregionMOKernel}, x::MOI, y::MOI
) where {MOI<:IsotopicMOInputsUnion}
    x.out_dim == y.out_dim ||
        throw(DimensionMismatch("`x` and `y` must have the same `out_dim`"))
    Kfeatures = kernelmatrix(k.kernel, x.x, y.x)
    Koutputs = _mo_output_covariance(k, x.out_dim)
    return _kernelmatrix_kroneckerjl_helper(MOI, Kfeatures, Koutputs)
end

function kronecker_kernelmatrix(
    k::Union{IndependentMOKernel,IntrinsicCoregionMOKernel}, x::MOI
) where {MOI<:IsotopicMOInputsUnion}
    Kfeatures = kernelmatrix(k.kernel, x.x)
    Koutputs = _mo_output_covariance(k, x.out_dim)
    return _kernelmatrix_kroneckerjl_helper(MOI, Kfeatures, Koutputs)
end

function kronecker_kernelmatrix(
    k::MOKernel, x::IsotopicMOInputsUnion, y::IsotopicMOInputsUnion
)
    return throw(
        ArgumentError("This kernel does not support a lazy kronecker kernelmatrix.")
    )
end

function kronecker_kernelmatrix(k::MOKernel, x::IsotopicMOInputsUnion)
    return kronecker_kernelmatrix(k, x, x)
end


text:
using MyPackage
using Test

@testset "my_func.jl" begin
    # Write your own tests here.
    @test my_func(1,2) == 4
    @test my_func(2,2) == 6
    @test my_func(2,3) == 7
end

@testset "derivative tests" begin
    @test derivative_of_my_func(2,1) == 2
end


text:
SRC_DIR = joinpath(Base.source_dir(), "..", "..", "src")
OUT_DIR = joinpath(Base.source_dir(), "..", "..", "output")
mkpath(OUT_DIR)

#include(joinpath(SRC_DIR, "BridgeSDEInference.jl"))
#using Main.BridgeSDEInference
include(joinpath(SRC_DIR, "BridgeSDEInference_for_tests.jl"))

using Distributions
using Random
using DataFrames
using CSV

DIR = "auxiliary"
include(joinpath(SRC_DIR, DIR, "read_and_write_data.jl"))
include(joinpath(SRC_DIR, DIR, "transforms.jl"))
include(joinpath(SRC_DIR, DIR, "utility_functions.jl"))

# decide if first passage time observations or partially observed diffusion
fptObsFlag = false

# pick dataset
filename = "path_part_obs_conj.csv"

# fetch the data
(df, x0, obs, obs_time, fpt,
      fptOrPartObs) = readData(Val(fptObsFlag), joinpath(OUT_DIR, filename))

param = :complexConjug
θ_init = [10.0, -8.0, 15.0, 0.0, 3.0]
P˟ = FitzhughDiffusion(param, θ_init...)
P̃ = std_aux_laws(FitzhughDiffusionAux, param, θ_init, obs, obs_time, 1)
L = @SMatrix [1. 0.]
Σ = @SMatrix [10^(-10)]

setup = MCMCSetup(P˟, P̃, fptOrPartObs)
set_observations!(setup, [L for _ in P̃], [Σ for _ in P̃], obs, obs_time, fpt)
set_imputation_grid!(setup, 1/1000)
set_transition_kernels!(setup,
                        [RandomWalk([],[]),
                         RandomWalk([0.0, 0.0, 0.0, 0.0, 0.5], 5)],
                        0.975, true, ((1,2,3),(5,)), # the first is memory param for pCN
                        (ConjugateUpdt(),
                         MetropolisHastingsUpdt(),
                         ))
set_priors!(setup,
            Priors((MvNormal([0.0,0.0,0.0], diagm(0=>[1000.0, 1000.0, 1000.0])),
                    ImproperPrior(),
                    )),
            GsnStartingPt(x0, @SMatrix [20. 0; 0 20.]),
            x0)
# 1*10^3, 3*10^2, 10^2, 10^0, 0
# num_mcmc_steps, save_iter, verb_iter, skip_for_save, warm_up
set_mcmc_params!(setup, 1*10^4, 3*10^2, 10^2, 10^0, 100) #1*10^4, 3*10^2, 10^2, 10^0, 100
set_blocking!(setup)
set_solver!(setup, Vern7(), NoChangePt())
initialise!(eltype(x0), setup)

Random.seed!(4)
out, elapsed = @timeit mcmc(setup)
display(out.accpt_tracker)

include(joinpath(SRC_DIR, DIR, "plotting_fns.jl"))
plot_chains(out; truth=[10.0, -8.0, 15.0, 0.0, 3.0])

plot_paths(out; transf=[(x,θ)->x, (x,θ)->conjugToRegular(x, θ[1], 0)],
           obs=(times=obs_time[2:end], vals=obs[2:end], indices=1))