text:
#import Random: rand!
import AQFED.Math: norminv

export AbstractSeq, ScrambledSobolSeq, next!, nextn!, scramble, skip, skipTo
export Owen, FaureTezuka, OwenFaureTezuka, NoScrambling

include("ssoboldata.jl") #loads `sobol_a` and `sobol_minit` from Joe & Kuo.
include("ssobolrng.jl") #original RNG and RNG adapter for scrambling

abstract type AbstractSeq{N} end

abstract type Scrambling end

struct NoScrambling <: Scrambling end

# Owen's recommendation is maxd >= L+10 where L is the log2 of maxNumSamples. maxd=31 scrambles all bits of an int.
struct Owen <: Scrambling
    maxd::Int
    rng::ScramblingRng
end

struct FaureTezuka <: Scrambling
    rng::ScramblingRng
end

struct OwenFaureTezuka <: Scrambling
    maxd::Int
    rng::ScramblingRng
end

unirand(s::Union{Owen,FaureTezuka,OwenFaureTezuka}) = unirand(s.rng)

# N iis the dimension of sequence being generated, S is scrambling type
mutable struct ScrambledSobolSeq{N,S} <: AbstractSeq{N}
    v::Array{UInt32,2} #array of size (sdim, log2(n))
    x::Array{UInt32,1}
    shift::Array{UInt32,1}
    l::UInt32
    counter::UInt32
    scrambling::S
end

ndims(s::AbstractSeq{N}) where {N} = N::Int

unirand(s::ScrambledSobolSeq) = unirand(s.scrambling)

@inline normalize(
    s::ScrambledSobolSeq{N,S},
    x::UInt32,
) where {N,S<:Union{NoScrambling,FaureTezuka}} = ldexp(Float64(x), -s.l % Int32)
@inline normalize(
    s::ScrambledSobolSeq{N,S},
    x::UInt32,
) where {N,S<:Union{Owen,OwenFaureTezuka}} = ldexp(Float64(x), -s.scrambling.maxd % Int32)

#fl1 find last 1 index in value
@inline function fl1(value::Int)
	return 32 - leading_zeros(value % UInt32)
end

@inline function fl1Slow(n::Int)
	l = 0
	i = n
	while (i != 0)
		i >>= 1
		l += 1
	end
	l = min(l, 31)
end

function ScrambledSobolSeq(dimension::Int, n::Int, scrambling::Scrambling)
    d = dimension
    (d < 0 || d > (length(sobol_a) + 1)) && error("invalid Sobol dimension")
    #special cases
    d == 0 && return (ScrambledSobolSeq{0})(v, UInt32[], UInt32[], zero(UInt32))
    l = 31
    if (n > 0)
    	l = fl1(n+1) #+1 as we skip the first point
    end
    x = Vector{UInt32}(undef, d)
    fill!(x, 0)
    counter = 0
    v = ones(UInt32, (d, l))
    for i = 1:l
        @inbounds v[1, i] = one(UInt32) << (l - i)
    end
    for j = 2:d
        a = sobol_a[j-1]
        s = floor(Int, log2(a)) #degree of poly
		a >>= 1
        #set initial values of m from table
        if l <= s
            for i = 1:l
                v[j, i] = sobol_minit[i, j-1] << (l - i)
            end
        else
            #fill in remaining values using recurrence
            for i = 1:s
                v[j, i] = sobol_minit[i, j-1] << (l - i)
            end
            for i = s+1:l
                @inbounds v[j, i] = v[j, i-s] ⊻ (v[j, i-s] >> s)
                for k = 1:s-1
                    @inbounds v[j, i] ⊻= (((a >> (s - k - 1)) & one(UInt32)) * v[j, i-k])
                end
            end
        end
    end
    vt = ones(UInt32, (l, d))
    shift = Vector{UInt32}(undef, d)
    ss = ScrambledSobolSeq{d,typeof(scrambling)}(v, x, shift, l, counter, scrambling)
    scramble(ss)
    ss
end

ScrambledSobolSeq(N::Integer) = ScrambledSobolSeq(Int(N), 0, NoScrambling())

function scramble(s::ScrambledSobolSeq{D,NoScrambling}) where {D}
    s.counter = 1 #skip first item for non-scrambled version
    fill!(s.shift, 0)
end

function scramble(s::ScrambledSobolSeq{D,Owen}) where {D}
    s.counter = 0
    fill!(s.shift, 0)
    scrambleOwen(s, s.scrambling.maxd)
end

function scramble(s::ScrambledSobolSeq{D,FaureTezuka}) where {D}
    s.counter = 0
    fill!(s.shift, 0)
    scrambleTezuka(s, s.l)
end

function scramble(s::ScrambledSobolSeq{D,OwenFaureTezuka}) where {D}
    s.counter = 0
    fill!(s.shift, 0)
    scrambleOwen(s, s.scrambling.maxd)
    scrambleTezuka(s, s.scrambling.maxd)
end

function scrambleOwen(s::ScrambledSobolSeq{D}, maxd::Integer) where {D}
    if maxd < s.l
        throw(DomainError(string("maxd must be >= l where l = log2(n) = ", s.l)))
    end
    scrambledV = zeros(UInt32, (D, s.l))
    #norm = 1/2^MAXD
    lsm = genscrml(s)
    @inbounds for i = 1:D
        @inbounds for j = 1:s.l
            l = 1
            temp2 = 0
            @inbounds for p = maxd:-1:1
                temp1 = 0
                @inbounds for k = 1:s.l
                    temp1 += lbitbits(lsm[i, p], k - 1, 1) * lbitbits(s.v[i, j], k - 1, 1)
                end
                temp1 %= 2
                temp2 += temp1 << (l - 1)
                l += 1
            end
            scrambledV[i, j] = temp2
        end
    end
    s.v .= scrambledV
    s.x .= s.shift
end

function scrambleTezuka(s::ScrambledSobolSeq{D}, maxx::Integer) where {D}
    ushift = Vector{UInt32}(undef, s.l)
    usm = genscrmu(s, ushift)
    tv = zeros(UInt32, (s.l, maxx, D))

    @inbounds for i = 1:D
        @inbounds for j = 1:s.l
            p = maxx
            @inbounds for k = 1:maxx
                tv[j, p, i] = lbitbits(s.v[i, j], k - 1, 1)
                p -= 1
            end
        end
        @inbounds for pp = 1:s.l
            temp2 = 0
            temp4 = 0
            l = 1
            @inbounds for j = maxx:-1:1
                temp1 = 0
                temp3 = 0
                @inbounds for p = 1:s.l
                    temp1 += tv[p, j, i] * usm[pp, p]
                    if (pp == 1)
                        temp3 += tv[p, j, i] * ushift[p]
                    end
                end
                temp1 %= 2
                temp2 += temp1 * l
                if (pp == 1)
                    temp3 %= 2
                    temp4 += temp3 * l
                end
                l <<= 1
            end
            s.v[i, pp] = temp2
            if (pp == 1)
                s.shift[i] ⊻= temp4
            end
        end
    end
    s.x .= s.shift
end

function genscrml(s::ScrambledSobolSeq{D,S}) where {D,S<:Union{Owen,OwenFaureTezuka}}
    maxd = s.scrambling.maxd
    lsm = zeros(UInt32, (D, maxd))
    @inbounds for p = 1:D
        s.shift[p] = 0
        l = 1
        @inbounds for i = maxd:-1:1
            lsm[p, i] = 0
            stemp = unirand(s)
            s.shift[p] += stemp << (l - 1)
            l += 1
            ll = 1
            @inbounds for j = s.l:-1:1
                local temp
                if (j == i)
                    temp = 1
                elseif (j < i)
                    temp = unirand(s)
                else
                    temp = 0
                end
                lsm[p, i] += temp << (ll - 1)
                ll += 1
            end
        end
    end
    return lsm
end

function genscrmu(
    s::ScrambledSobolSeq{D,S},
    ushift::Vector{UInt32},
) where {D,S<:Union{FaureTezuka,OwenFaureTezuka}}
    usm = zeros(UInt32, (s.l, s.l))
    @inbounds for i = 1:s.l
        stemp = unirand(s)
        ushift[i] = stemp
        @inbounds for j = 1:s.l
            local temp
            if (j == i)
                temp = 1
            elseif (j > i)
                temp = unirand(s)
            else
                temp = 0
            end
            usm[j, i] = temp
        end
    end
    return usm
end

@inline function lbitbits(a::UInt32, b::Int, len::Int)
    x = UInt64(a)
    y = UInt64(0xffffffffffffffff) #-1
    x >>= b
    y <<= len
    return (x & ~y) % UInt32
end

@inline function next(s::ScrambledSobolSeq, ::Type{UInt32})
    d = ndims(s)
    if s.counter == 0
        s.counter += one(s.counter)
        return s.x
    end
    c = ffz(s.counter)
    c > s.l &&
        throw(DomainError(string("counter larger than sequence length: ", s.counter)))
    s.counter += one(s.counter)
	#alternatively: @inbounds s.x .⊻= @view(s.v[:,c])
	sx = s.x
    sv = s.v
    @inbounds for j = 1:d
        sx[j] ⊻= sv[j, c]
    end
    return s.x
end

#next vector of uniform [0,1) numbers at counter containing all dimensions
@inline function next!(s::ScrambledSobolSeq, points::AbstractVector{<:AbstractFloat})
    next(s, UInt32)
    sx = s.x
    @inbounds for j = 1:ndims(s)
        # if sx[j] == 0 #may happen somewhere because of scrambling
        #     points[j] = normalize(s, one(UInt32)) / 2
        # else
            points[j] = normalize(s, sx[j])
        # end
    end
    return points
end

#next set of normal numbers
@inline function nextn!(s::ScrambledSobolSeq, points::AbstractVector{<:AbstractFloat})
	next!(s, points)
	@. points = norminv(points)
end

#next set of normal numbers
@inline function nextn!(s::ScrambledSobolSeq,  dim::Int, points::AbstractVector{<:AbstractFloat})
	next!(s, dim, points)
	@. points = norminv(points)
end

#next vector for a given dimension (vertical) from counter to counter + length(points)
@inline function next!(s::ScrambledSobolSeq, dim::Int, points::AbstractVector{<:AbstractFloat})
    j = dim
    sx = s.x
    sv = s.v
    @inbounds for i = 1:length(points)
        if s.counter != 0
            c = ffz(s.counter)
            sx[j] ⊻= sv[j, c]
        end
        points[i] =
        #    (sx[j] == 0) ? normalize(s, one(UInt32)) / 2 :
            normalize(s, sx[j])
        s.counter += one(s.counter)
    end
    return points
end

@inline function skip(s::ScrambledSobolSeq, n::Int)
    m = n + s.counter
    g = m ⊻ (m >> 1)
    d = ndims(s)
    sx = s.x
    sv = s.v
    @inbounds for j = 1:d
        @inbounds for index = 1:s.l
            if ((g >> index) & one(UInt32)) != 0
                sx[j] ⊻= sv[j, index]
            end
        end
    end
    s.counter += n
    s
end

@inline function skipTo(s::ScrambledSobolSeq, n::Int)
    # Convert to Gray code
    g = n ⊻ (n >> 1)
    d = ndims(s)
    sx = s.x
    sv = s.v
    @inbounds for j = 1:d
        s.x[j] = s.shift[j]
        @inbounds for index = 1:s.l
            if ((g >> index) & one(UInt32)) != 0
                sx[j] ⊻= sv[j, index]
            end
        end
    end
    s.counter = n
    s
end

#skip to position n for a given dimension dim
@inline function skipTo(s::ScrambledSobolSeq, dim::Int, n::Int)
    # Convert to Gray code
    g = n ⊻ (n >> 1)
    j = dim
    sx = s.x
    sv = s.v
    s.x[j] = s.shift[j]

    @inbounds for index = 1:s.l
        if ((g >> index) & one(UInt32)) != 0
            sx[j] ⊻= sv[j, index]
        end
    end
    s.counter = n
    s
end

@inline function ffz(x::Integer)
    # if x == 0 #not used as we check before if x is 0
    #     return 1
    # end
    return trailing_zeros(x) + 1
end


text:
@testset "SymmetricAllDifferent2AllDifferentInverse: $(fct_type), dimension $(dim), $(T)" for fct_type in ["vector of variables", "vector affine function"], dim in [2, 3], T in [Int, Float64]
    mock = MOIU.MockOptimizer(AllDifferentInverseModel{T}())
    model = COIB.SymmetricAllDifferent2AllDifferentInverse{T}(mock)

    if T == Int
        @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    end
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.AllDifferent,
    )
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Inverse,
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.SymmetricAllDifferent,
    )

    if T == Int
        x, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:dim])
    elseif T == Float64
        x = MOI.add_variables(model, dim)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables(x)
    elseif fct_type == "vector affine function"
        MOIU.vectorize(x)
    else
        @assert false
    end
    c = MOI.add_constraint(model, fct, CP.SymmetricAllDifferent(dim))

    for i in 1:dim
        @test MOI.is_valid(model, x[i])
    end
    @test MOI.is_valid(model, c)

    bridge = first(MOIBC.bridges(model))[2]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VariableIndex, CP.SymmetricAllDifferent) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == Tuple{Type}[]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.VectorAffineFunction{T}, CP.AllDifferent),
            (MOI.VectorAffineFunction{T}, CP.Inverse),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 0
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, CP.AllDifferent}()) == 1
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, CP.Inverse}()) == 1

        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, CP.AllDifferent}()) == [bridge.con_all_diff]
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, CP.Inverse}()) == [bridge.con_inverse]
    end

    @testset "All different" begin
        @test MOI.is_valid(model, bridge.con_all_diff)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_all_diff)
        @test length(f.terms) == dim
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_all_diff) == CP.AllDifferent(dim)

        for i in 1:dim
            t = f.terms[i]
            @test t.output_index == i
            @test t.scalar_term.coefficient === one(T)
            @test t.scalar_term.variable == x[i]
        end
    end

    @testset "Inverse" begin
        @test MOI.is_valid(model, bridge.con_inverse)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_inverse)
        @test length(f.terms) == 2 * dim
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_inverse) == CP.Inverse(dim)

        for i in 1:dim
            t = f.terms[i]
            @test t.output_index == i
            @test t.scalar_term.coefficient === one(T)
            @test t.scalar_term.variable == x[i]

            t = f.terms[dim + i]
            @test t.output_index == dim + i
            @test t.scalar_term.coefficient === one(T)
            @test t.scalar_term.variable == x[i]
        end
    end
end

ext: jl


text:
#!/usr/bin/julia
#
# Hamiltonian type definitions.
import LinearAlgebra: checksquare;
import Unitful: 𝐌, 𝐋;

export Hamiltonian, SeparableHamiltonian, set_min_zero;

"Abstract type of Hamiltonian matricies."
abstract type Hamiltonian{T<:Number,TM<:Number,M<:AbstractMatrix{TM}, S<:Space} <:
              AbstractMatrix{TM} end;
Base.size(H::Hamiltonian) = size(H.H);
Base.getindex(H::Hamiltonian, I::Int,J::Int) = H.H[I,J];
Base.setindex!(H::Hamiltonian, v, I::Int, J::Int) = (H.H[I,J] = v);

"Hamiltonian which can be separated into kinetic and potential parts."
struct SeparableHamiltonian{T,TM,M,S} <: Hamiltonian{T,TM,M,S}
  "Hamiltonian in the space `S`."
  H::M;

  "Diagonal components of potential part in configuration space."
  V::Vector{T};

  "Diagonal components on kinetic part in momentum space."
  T::Vector{T};

  """
      SeparableHamiltonian{T,TM,M,S}(H::M, V::AbstractVector,
                                    K::AbstractVector) where {T,TM,M,S}

  Re-implement default constructor with some error checking.
  """
  function SeparableHamiltonian{T,TM,M,S}(H::AbstractMatrix, V::AbstractVector,
                                          K::AbstractVector) where {T,TM,M,S}
    LinearAlgebra.checksquare(H);
    check_same_size(V, K);
    if size(H)[1] != length(V)
      throw(DimensionMismatch("passed matrix size does not match the vector length."));
    end
    new(H, V, K);
  end
end

"""
    SeparableHamiltonian{T}(g::Grid1D{S}) where {S,T}

Initialize blank Hamiltonian in the space `S` using given grid to determine
appropriate Hamiltonian size.
"""
function SeparableHamiltonian{T}(g::Grid1D{S}) where {T<:Number,S}
  len = length(g);
  SeparableHamiltonian{T,T,Matrix{T},S}(zeros(T, len, len),
                                      zeros(T, len), zeros(T, len));
end
SeparableHamiltonian(g::Grid1D) = SeparableHamiltonian{Complex{Float64}}(g);

"""
    SeparableHamiltonian(g::Grid1D{S}, H::AbstractMatrix{T}) where {S,T}

Initialise Hamiltonian in the space `S`, as the passed matrix.
"""
function SeparableHamiltonian(g::Grid1D{S}, H::AbstractMatrix{T}) where {S,T}
  LinearAlgebra.checksquare(H);
  check_same_size(g, H);
  len = length(g);
  SeparableHamiltonian{T,T,typeof(H),S}(H, zeros(T,len), zeros(T,len));
end

"""
    SeparableHamiltonian{T,,TM,M}(xgrid::Grid1D{ConfigurationSpace},
                        kgrid::Grid1D{MomentumSpace}) where M<:AbstractMatrix{T}

Construct Hamiltonian in configuration space with 'standard' KE component.

This assumes the mass of the particle and ħ are both 1.
"""
function SeparableHamiltonian{T,TM,M}(xgrid::Grid1D{C}, kgrid::Grid1D{S}) where
  {T,TM,M<:Matrix{TM},C<:ConfigurationSpace,S<:MomentumSpace}
  check_same_size(xgrid, kgrid);
  SeparableHamiltonian{T,TM,M,C}(-0.5sec_dv(xgrid), zeros(T,length(xgrid)), 0.5kgrid.^2);
end

"""
    function SeparableHamiltonian{T,TM,M}(xgrid::Grid1D{ConfigurationSpace,T1},
                                     kgrid::Grid1D{MomentumSpace,T2},
                                     mass::Number = 1) where
        {T,T1<:Unitful.Quantity,T2<:Unitful.Quantity}

Create Hamiltonian with energy units given by `T` in configuration space with
the 'standard' KE component.

This assumes ħ = 1.
"""
function SeparableHamiltonian{T,TM,M}(xgrid::Grid1D{ConfigurationSpace,T1},
                                 kgrid::Grid1D{MomentumSpace,T2}, m::Number=1) where
  {N<:Number,T<:Unitful.Quantity{N,Unitful.𝐌*Unitful.𝐋^2/Unitful.𝐓^2},TM<:Number,M<:Matrix{TM},
   T1<:Unitful.Quantity,T2<:Unitful.Quantity}
  check_same_size(xgrid, kgrid);

  SeparableHamiltonian{T,TM,M,ConfigurationSpace}(-0.5ustrip.(sec_dv(xgrid))/ustrip(m)*unit(T),
                                                zeros(T,length(xgrid)),
                                                0.5ustrip.(kgrid).^2/ustrip(m)*unit(T));
end

"Convenience constructor for setting default matrix types."
SeparableHamiltonian{T}(xgrid::Grid1D{ConfigurationSpace},
                        kgrid::Grid1D{MomentumSpace}) where T =
  SeparableHamiltonian{T,T,Matrix{T}}(xgrid, kgrid);
SeparableHamiltonian{T,TM}(xgrid::Grid1D{ConfigurationSpace},
                           kgrid::Grid1D{MomentumSpace}) where {T,TM} =
  SeparableHamiltonian{T,TM,Matrix{TM}}(xgrid, kgrid);
SeparableHamiltonian{T}(xgrid::Grid1D{ConfigurationSpace},
                        kgrid::Grid1D{MomentumSpace},m::Number) where T =
  SeparableHamiltonian{T,T,Matrix{T}}(xgrid, kgrid, m);
SeparableHamiltonian{T,TM}(xgrid::Grid1D{ConfigurationSpace},
                           kgrid::Grid1D{MomentumSpace},m::Number) where {T,TM} =
  SeparableHamiltonian{T,TM,Matrix{TM}}(xgrid, kgrid, m);
SeparableHamiltonian(xgrid::Grid1D{ConfigurationSpace},
                     kgrid::Grid1D{MomentumSpace}) =
  SeparableHamiltonian{Complex{Float64},Complex{Float64},Matrix{Complex{Float64}}}(xgrid, kgrid);

"""
    SeparableHamiltonian{T,TM,M}(xgrid::Grid1D{ConfigurationSpace},
                            kgrid::Grid1D{MomentumSpace},
                            V::AbstractVector)

Construct Hamiltonian in configuration space with 'standard' KE component and
a potential given by `V`.

This assumes the mass of the particle and ħ are 1.
"""
function SeparableHamiltonian{T,TM,M}(xgrid::Grid1D{C}, kgrid::Grid1D{S},
                                   V::AbstractVector) where
    {T,TM,M<:Matrix{TM},C<:ConfigurationSpace,S<:MomentumSpace}
  check_same_size(xgrid, kgrid);
  check_same_size(xgrid, V);

  SeparableHamiltonian{T,TM,M,C}(-0.5sec_dv(xgrid) + LinearAlgebra.Diagonal(V),
                              V, 0.5kgrid.^2);
end

SeparableHamiltonian{T}(xgrid::Grid1D{ConfigurationSpace},
                        kgrid::Grid1D{MomentumSpace}, V::AbstractVector) where T =
  SeparableHamiltonian{T,T,Matrix{T}}(xgrid, kgrid, V);
SeparableHamiltonian{T,TM}(xgrid::Grid1D{ConfigurationSpace},
                           kgrid::Grid1D{MomentumSpace}, V::AbstractVector) where {T,TM} =
  SeparableHamiltonian{T,TM,Matrix{TM}}(xgrid, kgrid, V);
SeparableHamiltonian(xgrid::Grid1D{ConfigurationSpace},
                     kgrid::Grid1D{MomentumSpace}, V::AbstractVector{T}) where T =
  SeparableHamiltonian{T,T,Matrix{T}}(xgrid, kgrid, V);

"""
    SeparableHamiltonian{T,TM,M}(g::Grid1D{ConfigurationSpace}, V::AbstractVector)

Construct Hamiltonian in configuration space with the potential passed by `V`.
"""
function SeparableHamiltonian{T,TM,M}(g::Grid1D{S}, V::AbstractVector) where
      {T,TM,M<:Matrix{TM},S<:ConfigurationSpace}
  check_same_size(g, V);
  SeparableHamiltonian{T,TM,M,S}(LinearAlgebra.Diagonal(V), V, zeros(T,length(g)));
end

"""
    SeparableHamiltonian{T,TM,M}(g::Grid1D{MomentumSpace}, K::AbstractVector)

Construct Hamiltonian in momentum space with the kinetic component passed by `K`.
"""
function SeparableHamiltonian{T,TM,M}(g::Grid1D{S}, K::AbstractVector) where
      {T,TM,M<:Matrix{TM},S<:MomentumSpace}
  check_same_size(g, K);
  SeparableHamiltonian{T,TM,M,S}(LinearAlgebra.Diagonal(K), zeros(T,length(grid)), K);
end

"""
    SeparableHamiltonian(g::Grid1D{S}, V::AbstractVector{T}) where S<:Space

Construct Hamiltonian, with elements of type `T`, in the space `S` with
the momentum/potential component passed by `V` (depending on the space).
"""
SeparableHamiltonian(g::Grid1D, V::AbstractVector{T}) where T =
  SeparableHamiltonian{T,T,Matrix{T}}(g, V);
SeparableHamiltonian{T}(g::Grid1D, V::AbstractVector) where T =
  SeparableHamiltonian{T,T,Matrix{T}}(g, V);
SeparableHamiltonian{T,TM}(g::Grid1D, V::AbstractVector) where {T,TM} =
  SeparableHamiltonian{T,TM,Matrix{TM}}(g, V);

"""
    Base.:+(H₁::SeparableHamiltonian, H₂::SeparableHamiltonian)

Adding two `SeparableHamiltonian` objects of the same type returns a new
`SeparableHamiltonian`.
"""
Base.:+(H₁::SeparableHamiltonian{S,T,TM,M},
        H₂::SeparableHamiltonian{S,T,TM,M}) where {S,T,TM,M} =
  SeparableHamiltonian{S,T,TM,M}(H₁.H + H₂.H, H₁.V + H₂.V, H₁.T + H₂.T);
"""
    Base.:-(H₁::SeparableHamiltonian, H₂::SeparableHamiltonian)

Subtracting two `SeparableHamiltonian` objects of the same type returns a new
`SeparableHamiltonian`.
"""
Base.:-(H₁::SeparableHamiltonian{S,T,TM,M},
        H₂::SeparableHamiltonian{S,T,TM,M}) where {S,T,TM,M} =
  SeparableHamiltonian{S,T,TM,M}(H₁.H - H₂.H, H₁.V - H₂.V, H₁.T - H₂.T);

"""
    Base.:*(H::SeparableHamiltonian, m::Real)

Multiplying a `SeparableHamiltonian` object by a real scalar returns a
`SeparableHamiltonian` object.
"""
Base.:*(H::SeparableHamiltonian{S,T,TM,M}, m::Real) where {S,T,TM,M} =
  SeparableHamiltonian{S,T,TM,M}(H.H * m, H.V * m, H.T * m);
Base.:*(m::Real, H::SeparableHamiltonian) = Base.:*(H, m);

"""
    Base.:/(H::SeparableHamiltonian, m::Real)

Dividing a `SeparableHamiltonian` object by a real scalar returns a
`SeparableHamiltonian` object.
"""
Base.:/(H::SeparableHamiltonian{S,T,TM,M}, m::Real) where {S,T,TM,M} =
  SeparableHamiltonian{S,T,TM,M}(H.H / m, H.V / m, H.T / m);
Base.:/(m::Real, H::SeparableHamiltonian) = Base.:/(H, m);

"""
Hamiltonian constructors alias to the `SeparableHamiltonian` constructors.

See also: [`SeparableHamiltonian`](@ref)
"""
Hamiltonian(g::Grid1D) = SeparableHamiltonian(g);
Hamiltonian(xgrid::Grid1D{ConfigurationSpace}, kgrid::Grid1D{MomentumSpace}) =
  SeparableHamiltonian(xgrid, kgrid);
Hamiltonian(xgrid::Grid1D{ConfigurationSpace}, kgrid::Grid1D{MomentumSpace},
            V::AbstractVector) = SeparableHamiltonian(xgrid, kgrid, V);
Hamiltonian{T}(xgrid::Grid1D{ConfigurationSpace}, kgrid::Grid1D{MomentumSpace},
               V::AbstractVector) where T =
  SeparableHamiltonian{T}(xgrid, kgrid, V);
Hamiltonian(g::Grid1D, V::AbstractVector) = SeparableHamiltonian(g, V);
Hamiltonian{T}(g::Grid1D, V::AbstractVector) where T = SeparableHamiltonian{T}(g, V);

"""
    check_same_size(g::Grid1D, H::AbstractMatrix)

Check the passed matrix has first dimension equal to grid size.
"""
function check_same_size(g::Grid1D, H::AbstractMatrix)
  if (length(H[1,:]) != length(g))
    throw(DimensionMismatch("passed matrix size does not match grid size."));
  end
end
function check_same_size(H::AbstractMatrix, g::Grid1D)
  check_same_size(g, H)
end
# Note: functions implementing check_same_size() for Grid and AbstractArray
# are implemented in waveobjects.jl

"""
    set_min_zero(H::SeparableHamiltonian)

Shift the potential part of the Hamiltonian such that the minimum potential
value is zero.
"""
function set_min_zero(H::SeparableHamiltonian{S,T,TM,M}) where {S,T,TM,M}
  min = minimum(real.(H.V));
  V = H.V .- min;
  return SeparableHamiltonian{S,T,TM,M}(H .- LinearAlgebra.Diagonal(V), V, H.T);
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export transform_neuron_fn
"""
Data used by `neuron_fn_transformer` during the ASTWalk process

    `args`  : Stores the neuron fields encountered in the body of the function
              which should be passed as arguments to the generated function
"""
type NeuronTransformerData
    ensemble     :: AbstractEnsemble
    args         :: Set
    NeuronTransformerData(ensemble::AbstractEnsemble) =
        new(ensemble, Set())
end

function is_param(name::Symbol, ensemble::AbstractEnsemble)
    for param in ensemble.params
        if param.name == name
            return true
        end
    end
    return false
end

@doc """
Process a neuron function replacing field references with array reference
expressions
"""
function transform_neuron_fn(fn, ensemble)
    function walker(node, cbdata, index, top_level, read)
        if !isa(node, Expr)
            return ASTWALK_RECURSE
        end
        if node.head == :(.) && node.args[1] == :neuron
            name = node.args[2].value
            N = ndims(cbdata.ensemble)
            if name == :index
                # N += 1  # batch_dim
                idx = :($(symbol(:_neuron_index_,N)) - 1)
                buffer = symbol(cbdata.ensemble.name, :value)
                for i in N-1:-1:1
                    size = :(size($buffer, $i))
                    idx = :($idx * $size + $(symbol(:_neuron_index_,i)) - 1)
                end
                return :($idx + 1)
            end
            if name in cbdata.ensemble.batch_fields
                N += 1
            end
            name = symbol(cbdata.ensemble.name,name)
            idx = Any[symbol(:_neuron_index_,i) for i in 1:N]
            str_name = string(name)
            if !LOSSY_GRADIENTS && contains(str_name, "∇")
                result = split(str_name, "∇")
                if result[2] != "" && !contains(result[2], "inputs")
                    push!(idx, :(_omp_get_thread_num() + 1))
                end
            end
            if !contains(str_name, "inputs")
                push!(cbdata.args, name)
            end
            return :($name[$(idx...)])
        elseif node.head == :ref
            for i in 2:length(node.args)
                node.args[i] = AstWalk(node.args[i], walker, cbdata)
            end
            result = AstWalk(node.args[1], walker, cbdata)
            str_target = string(result.args[1])
            if endswith(str_target, "inputs")
                @assert length(node.args[3:end]) == 0
                index = result.args[2:end]
                if isa(cbdata.ensemble, ActivationEnsemble)
                    if contains(str_target, "∇")
                        result.args[1] = symbol(cbdata.ensemble.name,:∇)
                    else
                        result.args[1] = symbol(cbdata.ensemble.name,:value)
                    end
                else
                    conn_index = node.args[2]
                    result.args[1] = symbol(result.args[1], conn_index)
                end
                push!(cbdata.args, result.args[1])
                node = Expr(:ref, result.args[1], index...)
            elseif contains(str_target, "inputs")
                conn_index = parse(Int, split(str_target, "inputs")[end])
                if cbdata.ensemble.connections[conn_index].is_one_to_one
                    # FIXME: Needs tiling enabled check
                    @assert length(node.args[2:end]) == 1 && node.args[2] == 1
                    if result.args[end] != :NOTILE
                        push!(result.args, :NOTILE)
                    end
                    return result
                else
                    node = Expr(:ref, result.args[1], node.args[2:end]..., result.args[2:end]...)
                end
            else
                node = Expr(:ref, result.args[1], node.args[2:end]..., result.args[2:end]...)
            end
            return node
        elseif node.head == :call && node.args[1] == :length
            node.args[2] = AstWalk(node.args[2], walker, cbdata)
            if isa(node.args[2], Expr) && node.args[2].head == :ref
                if contains(node.args[2].args[1], "inputs")
                    return Expr(:call, :size, node.args[2].args[1], 1)
                end
            end
        end
        ASTWALK_RECURSE
    end
    cbdata = NeuronTransformerData(ensemble)
    ast = AstWalk(fn, walker, cbdata)
    ast, cbdata.args
end


text:
include("potentially_intersecting_simplices.jl")
include("reshape_simplices.jl")
include("function_barriers.jl")

import ..AbstractTransferOperator
import StateSpaceReconstruction:
    DelaunayTriangulation,
    Simplex,
	point_contained
import StateSpaceReconstruction: subsample_coeffs
import StaticArrays: Size
import Statistics: det
import InplaceOps: @!

#####################################
# Exports (don't export this estimator, the other one is more efficient)
#####################################
#export transferoperator_triangulation_approx

#####################################
# Constructors
#####################################

function transferoperator_triangulation_approxSlow(inv_pts;
                            n_sample_pts::Int = 200,
                            sample_randomly::Bool = false)

    # Triangulate all points but the last (we need the last point
    # for the forward projection of the simplices)
    DT = DelaunayTriangulation(inv_pts[1:(end - 1)])

    # Some constants used throughout the function
    n_simplices = length(DT)
    dim = length(DT[1]) - 1

    #=
    # Prepare memory-efficient representations of the simplices
    =#
    simplices = reshape_simplices(inv_pts[1:(end- 1)], DT)
    image_simplices = reshape_simplices(inv_pts[2:end], DT)

    #=
    Compute the convex expansions coefficients needed to generate points
    within a simplex. Also, update number of points if a shape-preserving
    simplex subdivision algorithm was used (in that case, we get a
    *minimum* of `n_sample_pts` set of coefficients.
    =#
    convex_coeffs = subsample_coeffs(dim, n_sample_pts, sample_randomly)
    n_coeffs::Int = size(convex_coeffs, 2)
    convex_coeffs = [convex_coeffs[:, i] for i = 1:n_coeffs]

    # Pre-allocated arrays (SizedArrays, for efficiency)
    pt          = Size(1, dim)(zeros(Float64, 1, dim))
    s_arr       = Size((dim+1)^2)(zeros(Float64, (dim+1)^2))
    signs       = Size(dim + 1)(zeros(Float64, dim + 1))

    # Re-arrange simplices so that look-up is a bit more efficient
    simplex_arrs = Vector{Array{Float64, 2}}(undef, n_simplices)
    imsimplex_arrs = Vector{Array{Float64, 2}}(undef, n_simplices)

    for i in 1:n_simplices
        simplex_arrs[i] = hcat(inv_pts[DT[i]]...,)
        imsimplex_arrs[i] = hcat(inv_pts[DT[i] .+ 1]...,)
    end

    # The Markov matrix
    M = zeros(Float64, n_simplices, n_simplices)

    for i in 1:n_simplices
        # Find the indices of the simplices who potentially intersect
        # with the image of simplex #i
        idxs = idxs_potentially_intersecting_simplices(inv_pts, DT, i)
        Sj = Vector{AbstractArray}(undef, length(idxs))
        get_simplices_at_inds!(Sj, idxs, simplices)

        # Generate points within the image of simplex #i and check
        # what fraction of those points fall into the simplices.
        # We approximate the intersecting volume between the
        # image simplex and the individual simplices by those
        # fractions.
        @views is = imsimplex_arrs[i]
        for k in 1:n_coeffs
            @! pt = transpose(convex_coeffs[k]) * transpose(is)
            innerloop!(idxs, signs, s_arr, Sj, pt, dim, M, i)
        end
    end
    # Need to normalise, because all we have up until now is counts
    # of how many points inside the image simplex falls into
    # the simplices.
    return transferoperator_triangulation_approx(transpose(M) ./ n_coeffs)
end

"""
Approximate the transfer operator for a set of points by triangulating them,
then compute approximate simplex intersections.
"""

function transferoperator_triangulation_approxNaive(invariant_pts;
        n_sample_pts::Int = 100, sample_randomly::Bool = true)
    original_pts = invariant_pts[1:end-1]
    forward_pts =  invariant_pts[2:end]

    DT = DelaunayTriangulation(invariant_pts[1:end-1])
    n_simplices = length(DT)
    dim = length(DT[1]) - 1

    # Collect all simplices and image simplices before the loop.
    # For all but the smallest triangulations, this is way faster
    # than gathering the points and creating new simplices at every
    # step of the loop.
    simplices = [Simplex(original_pts[DT[i]]) for i = 1:n_simplices]
    image_simplices = [Simplex(forward_pts[DT[i]]) for i = 1:n_simplices]

    #=
    Convex coefficients for computing interior points within the simplices.
    =#
    convex_coeffs = subsample_coeffs(dim, n_sample_pts, sample_randomly)

    #=
    # update number of points in case the subsampled points were generated
    # using a shape-preserving splitting of the parent simplex
    # (if so, because the number of subsimplices generated by the
    # shape-preserving splitting depends only on the dimension of the space,
    # there will be more points than we asked for). (If points are sampled
    # randomly, we will have exactly as many points within the simplex as
    # we asked for, and the number of points doesn't need to be adjusted)
    =#
    n_sample_pts = size(convex_coeffs, 2)

    @show n_sample_pts

    TO = zeros(Float64, n_simplices, n_simplices)

    interior_pt = zeros(Float64, 3)
    image_simplex_arr = zeros(Float64, dim, dim + 1)
    temp_simplex = MutableSimplex(zeros(Float64, dim, dim + 1))

    for j in 1:n_simplices
        image_simplex_arr[:, :] = image_simplices[j][:, :]

        for i in 1:n_simplices
            simplex = simplices[i]

            # Generate points within the simplex and check how many of them
            # lie inside the current image simplex.
            for k = 1:n_sample_pts
                @! interior_pt = image_simplex_arr * convex_coeffs[:, k]

                if point_contained!(temp_simplex, simplex, interior_pt)
                    TO[j, i] += 1.0
                end
            end
        end
    end

    return TransferOperatorTriangulationApprox(TO ./ n_sample_pts)
end


text:
using DrWatson, Test, Parameters, Dates

# Define a type we use at experiments
abstract type Species end
struct Mouse <: Species end
struct Cat <: Species end

@with_kw struct Experiment{S<:Species}
    n::Int = 50
    c::Float64 = 10.0
    x::Float64 = 0.2
    date::Date = Date(1991, 04, 13)
    species::S = Mouse()
    scientist::String = "George"
end
# Make a couple of experiments
e1 = Experiment()
e2 = Experiment(species = Cat())

# Implement defaults
DrWatson.default_prefix(e::Experiment) = "Experiment_"*string(e.date)
DrWatson.allaccess(::Experiment) = (:n, :c, :x, :species)

# Extend `Base.string` for custom types (`savename` uses `string`)
Base.string(::Mouse) = "mouse"
Base.string(::Cat) = "cat"

@test savename(e1) == "Experiment_1991-04-13_c=10_n=50_x=0.2"
@test savename(e2) == savename(e1)

# Add the extra type:
DrWatson.default_allowed(::Experiment) = (Real, String, Species)

@test savename(e1) ≠ "Experiment_1991-04-13_c=10_n=50_x=0.2"
@test savename(e2) ≠ savename(e1)

@test savename(e1) == "Experiment_1991-04-13_c=10_n=50_species=mouse_x=0.2"
@test savename(e2) == "Experiment_1991-04-13_c=10_n=50_species=cat_x=0.2"


text:
"Parameters for right-hand side function of variational partitioned Runge-Kutta methods."
mutable struct ParametersVPRKpSymmetric{DT, TT, ET <: IODE{DT,TT}, D, S} <: AbstractParametersVPRK{DT,TT,ET,D,S}
    equ::ET
    tab::TableauVPRK{TT}
    Δt::TT

    R::Vector{TT}

    t̅::TT
    q̅::Vector{DT}
    p̅::Vector{DT}
end

function ParametersVPRKpSymmetric(equ::ET, tab::TableauVPRK{TT}, Δt::TT) where {DT, TT, ET <: IODE{DT,TT}}
    R = convert(Vector{TT}, [1, tab.R∞])
    q̅ = zeros(DT, ndims(equ))
    p̅ = zeros(DT, ndims(equ))

    ParametersVPRKpSymmetric{DT, TT, ET, ndims(equ), tab.s}(equ, tab, Δt, R, zero(TT), q̅, p̅)
end

"Variational partitioned Runge-Kutta integrator."
struct IntegratorVPRKpSymmetric{DT, TT, PT <: ParametersVPRKpSymmetric{DT,TT},
                                        ST <: NonlinearSolver{DT},
                                        IT <: InitialGuessPODE{DT,TT}, D, S} <: AbstractIntegratorVPRKwProjection{DT,TT}

    params::PT
    solver::ST
    iguess::IT
    cache::IntegratorCacheVPRK{DT,D,S}
end

function IntegratorVPRKpSymmetric(equation::ET, tableau::TableauVPRK{TT}, Δt::TT) where {DT, TT, ET <: IODE{DT,TT}}
    D = equation.d
    M = equation.n
    S = tableau.s

    # create params
    params = ParametersVPRKpSymmetric(equation, tableau, Δt)

    # create solver
    solver = create_nonlinear_solver(DT, D*(S+2), params)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create cache
    cache = IntegratorCacheVPRK{DT,D,S}(true)

    # create integrator
    IntegratorVPRKpSymmetric{DT, TT, typeof(params), typeof(solver), typeof(iguess), D, S}(params, solver, iguess, cache)
end

function compute_projection_vprk!(x::Vector{ST},
                q::SolutionVector{ST}, p::SolutionVector{ST}, λ::SolutionVector{ST},
                V::Vector{Vector{ST}}, U::Vector{Vector{ST}}, G::Vector{Vector{ST}},
                params::ParametersVPRKpSymmetric{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    local t₀::TT = params.t̅
    local t₁::TT = params.t̅ + params.Δt

    # copy x to λ and q̅
    for k in 1:D
        q[k] = x[D*(S+0)+k]
        λ[k] = x[D*(S+1)+k]
    end

    # compute U=λ
    U[1] .= λ
    U[2] .= λ

    # compute G=g(q,λ)
    params.equ.g(t₀, params.q̅, λ, G[1])
    params.equ.g(t₁, q, λ, G[2])

    # compute p=ϑ(q)⋅p
    params.equ.ϑ(t₁, q, λ, p)
end

"Compute stages of variational partitioned Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST},
                params::ParametersVPRKpSymmetric{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    cache = IntegratorCacheVPRK{ST, D, S}(true)

    quote
        compute_stages!(x, $cache.q̃, $cache.p̃, $cache.λ, $cache.Q, $cache.V, $cache.U, $cache.P, $cache.F, $cache.G, params)

        # compute b = - [P-AF-U]
        compute_rhs_vprk!(b, $cache.P, $cache.F, $cache.G, params)

        # compute b = - [p-bF-G]
        compute_rhs_vprk_projection_p!(b, $cache.p̃, $cache.F, $cache.G, D*(S+0), params)

        # compute b = - [q-bV-U]
        compute_rhs_vprk_projection_q!(b, $cache.q̃, $cache.V, $cache.U, D*(S+1), params)

        compute_rhs_vprk_correction!(b, $cache.V, params)
    end
end

function initial_guess!(int::IntegratorVPRKpSymmetric{DT,TT}, sol::AtomicSolutionPODE{DT,TT}) where {DT,TT}
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                              sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                              int.cache.q̃, int.cache.ṽ,
                              tableau(int).q.c[i])

        for k in eachdim(int)
            int.solver.x[ndims(int)*(i-1)+k] = int.cache.ṽ[k]
        end
    end

    evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                          sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                          int.cache.q̃, one(TT))
    for k in eachdim(int)
        int.solver.x[ndims(int)*(nstages(int)+0)+k] = int.cache.q̃[k]
    end
    for k in eachdim(int)
        int.solver.x[ndims(int)*(nstages(int)+1)+k] = 0
    end
end

"Integrate ODE with variational partitioned Runge-Kutta integrator."
function integrate_step!(int::IntegratorVPRKpSymmetric{DT,TT}, sol::AtomicSolutionPODE{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, sol)

    # compute initial guess
    initial_guess!(int, sol)

    # reset solution
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages and projection vector fields
    compute_stages!(int.solver.x,
                    int.cache.q̃, int.cache.p̃, int.cache.λ,
                    int.cache.Q, int.cache.V, int.cache.U,
                    int.cache.P, int.cache.F, int.cache.G, int.params)

    # compute unprojected solution
    update_solution!(int, sol)

    # add projection to solution
    project_solution!(int, sol, int.params.R)

    # copy solution to initial guess
    update!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f)
end


text:
function getInterpolator(A::AbstractArray{T,3}, grid::RegularGridPositions) where T

  tmp = ( size(A,1)==1 ? NoInterp() : BSpline(Linear()),
          size(A,2)==1 ? NoInterp() : BSpline(Linear()),
          size(A,3)==1 ? NoInterp() : BSpline(Linear()) )

  itp = extrapolate(interpolate(A, tmp), 0.0)
  sitp = scale(itp, range(grid,1), range(grid,2), range(grid,3))
  return sitp
end

function interpolate(A::AbstractArray{T,3}, origin::RegularGridPositions,
                        target::RegularGridPositions) where T

  sitp = getInterpolator(A, origin)
  N = target.shape
  AInterp = zeros(eltype(A), N[1], N[2], N[3])
  rx, ry, rz = range(target,1), range(target,2), range(target,3)

  return _interpolate_inner(AInterp,N,sitp,rx,ry,rz)
end

function _interpolate_inner(AInterp,N,sitp,rx,ry,rz)
  for nz=1:N[3]
    for ny=1:N[2]
      for nx=1:N[1]
        AInterp[nx,ny,nz] = sitp(rx[nx],ry[ny],rz[nz])
      end
    end
  end

  return AInterp
end


text:
######################################################
# Common objects and methods used for all algorithms
######################################################

"""
    Workspace

Abstract Workspace, which is a struct that stores states and information needed for optimization process.
Implementation depends on specific algorithm.
"""
abstract type Workspace end
function reset!(w::Workspace, x0 = nothing)
    if x0 !== nothing
        w.x0 .= x0
    end
    return w
end

"""
    ConvergenceState

A struct that summarizes the convergence state of a solution. The fields in this struct are:
 - `Δx`: the infinity norm of the change in the solution `x`.
 - `Δf`: the change in the objective value `f`.
 - `relΔf`: the ratio of change in the objective value `f`.
 - `kkt_residual`: the Karush-Kuhn-Tucker (KKT) residual of the solution. If [`ScaledKKTCriteria`](@ref) is used instead of [`KKTCriteria`](@ref), the `kkt_residual` will be divided by a factor.
 - `ipopt_residual`: the modified KKT residual used in the IPOPT solver.
 - `infeas`: maximum infeasibility amount. This is 0 if the solution is feasible.
 - `x_converged`: true if `Δx` is less than the `x` tolerance in [`MMAOptions`](@ref).
 - `fabs_converged`: true if `Δf` is less than the `f` tolerance in [`MMAOptions`](@ref).
 - `frel_converged`: true if `relΔf` is less than the `f` tolerance in [`MMAOptions`](@ref).
 - `kkt_converged`: true if the `kkt_residual` is less than the KKT tolerance in [`MMAOptions`](@ref).
 - `ipopt_converged`: true if the `ipopt_residual` is less than the KKT tolerance in [`MMAOptions`](@ref).
 - `infeas_converged`: true if `infeas` is less than the infeasibility tolerance in [`MMAOptions`](@ref).
 - `f_increased`: true if the objective value of the current solution is higher than that of the previous solution.
 - `converged`: true if the solution satisfies the convergence criteria of choice. See [`ConvergenceCriteria`](@ref) for more on the different convergence criteria available.
"""
@with_kw mutable struct ConvergenceState{T}
    Δx::T = Inf
    Δf::T = Inf
    relΔf::T = Inf
    kkt_residual::T = Inf
    ipopt_residual::T = Inf
    infeas::T = Inf
    x_converged::Bool = false
    fabs_converged::Bool = false
    frel_converged::Bool = false
    kkt_converged::Bool = false
    ipopt_converged::Bool = false
    infeas_converged::Bool = false
    f_increased::Bool = true
    converged::Bool = false
end
ConvergenceState(::Type{T}) where {T} = ConvergenceState{T}()

"""
    Tolerance

A struct specifying the different tolerances used to assess the convergence of the algorithms. The following are the fields of `Tolerance`:
 - `x`: the tolerance for `Δx` in [`ConvergenceState`](@ref). `x_converged` will be true if `Δx` is less than the `x` tolerance in `Tolerance`. This is used to assess convergence when the [`GenericCriteria`](@ref) is used as the convergence criteria.
 - `fabs`: the tolerance for `Δf` in [`ConvergenceState`](@ref). `f_converged` will be true if `Δf` is less than the `fabs` tolerance in `Tolerance`. This is used to assess convergence when the [`GenericCriteria`](@ref) is used as the convergence criteria.
 - `frel`: the tolerance for `relΔf` in [`ConvergenceState`](@ref). `f_converged` will be true if `relΔf` is less than the `frel` tolerance in `Tolerance`. This is used to assess convergence when the [`GenericCriteria`](@ref) is used as the convergence criteria.
 - `kkt`: the KKT tolerance. `kkt_converged` in [`ConvergenceState`](@ref) will be true if the `kkt_residual` is less than the KKT tolerance. And `ipopt_converged` in [`ConvergenceState`](@ref) will be true if `ipopt_residual` is less than the KKT tolerance. This is used to assess convergence when the [`KKTCriteria`](@ref), the [`ScaledKKTCriteria`](@ref) or [`IpoptCriteria`](@ref) criteria is used as the convergence criteria.
 - `infeas`: the maximum infeasibility tolerance. `infeas_converged` in [`ConvergenceState`](@ref) will be true if the maximum infeasibility is less than the infeasibility tolerance. This is used to assess convergence regardless of the convergence criteria used.

For more on convergence criteria, see [`GenericCriteria`](@ref), [`KKTCriteria`](@ref), [`ScaledKKTCriteria`](@ref) and [`IpoptCriteria`](@ref).
"""
struct Tolerance{Tx, Tf, Tkkt, Tinfeas}
    x::Tx
    fabs::Tf
    frel::Tf
    kkt::Tkkt
    infeas::Tinfeas
end
function Tolerance(;
    x = 0.0, f = 1e-5, fabs = f, frel = f,
    kkt = 1e-5, infeas = 1e-5,
)
    Tolerance(x, fabs, frel, kkt, infeas)
end
function (tol::Tolerance{<:Function, <:Function, <:Function})(i)
    return Tolerance(tol.x(i), tol.fabs(i), tol.frel(i), tol.kkt(i), tol.infeas)
end

function Base.:*(t::Tolerance, m::Real)
    return Tolerance(t.x * m, t.fabs * m, t.frel * m, t.kkt * m, t.infeas * m)
end

"""
    ConvergenceCriteria

This an abstract type with 4 subtypes:
1. [`GenericCriteria`](@ref)
2. [`KKTCriteria`](@ref)
3. [`ScaledKKTCriteria`](@ref)
4. [`IpoptCriteria`](@ref)
"""
abstract type ConvergenceCriteria end

"""
    GenericCriteria

This is a generic convergence criteria that uses:
1. The maximum change in the solution, `Δx`,
2. The change in the objective value, `Δf`, and
3. The change percentage in the objective value, `Δf`, and
4. The maximum infeasibility `infeas`.
to assess convergence. More details are given in [`assess_convergence!`](@ref).
"""
struct GenericCriteria <: ConvergenceCriteria end

"""
    KKTCriteria

This convergence criteria uses the Karush-Kuhn-Tucker residual and maximum infeasibility to assess convergence. More details are given in [`assess_convergence!`](@ref).
"""
struct KKTCriteria <: ConvergenceCriteria end

"""
    IpoptCriteria

This convergence criteria uses a scaled version of the Karush-Kuhn-Tucker (KKT) residual and maximum infeasibility to assess convergence. This scaled KKT residual is used in the IPOPT nonlinear programming solver as explained in [this paper](http://cepac.cheme.cmu.edu/pasilectures/biegler/ipopt.pdf). More details are given in [`assess_convergence!`](@ref).
"""
struct IpoptCriteria <: ConvergenceCriteria end

"""
    ScaledKKTCriteria

This convergence criteria uses another scaled version of the Karush-Kuhn-Tucker (KKT) residual and maximum infeasibility to assess convergence. In particular if the objective was scaled by a factor `m`, the KKT residual will be scaled down by a factor `max(m, 1/m)`. This scaling was found to make the convergence criteria less sensitive to scale compared to using the traditional KKT residual. More details are given in [`assess_convergence!`](@ref).
"""
struct ScaledKKTCriteria <: ConvergenceCriteria end

"""
    Solution

A struct that stores all the information about a solution. The following are the fields of `Solution`:
 - `x`: the current primal solution
 - `prevx`: the previous primal solution
 - `λ`: the current dual solution
 - `f`: the objective value of `x`
 - `prevf`: the obejctive value of `prevx`
 - `∇f`: the gradient of the objective at `x`
 - `g`: the constraint function value at `x`
 - `∇g`: the Jacobian of the constraint functions at `x`
 - `convstate`: the convergence state of the solution
"""
@params mutable struct Solution
    prevx
    x
    λ
    prevf
    f
    ∇f
    g
    ∇g
    convstate::ConvergenceState
end

"""
    Solution(dualmodel, λ)

Construct an empty solution for the dual model `dualmodel` given a sample dual solution `λ`.
"""
function Solution(dualmodel, λ)
    prevx = copy(getxk(dualmodel))
    x = copy(prevx)
    λ = copy(λ)
    prevf = Inf
    fg = getfk(dualmodel)
    ∇fg = get∇fk(dualmodel)
    f = fg[1]
    g = fg[2:end]
    ∇f = ∇fg[1, :]
    ∇g = ∇fg[2:end, :]
    convstate = ConvergenceState()
    return Solution(prevx, x, λ, prevf, f, ∇f, g, ∇g, convstate)
end

"""
    AbstractResult

An abstract type that stores optimization result.
"""
abstract type AbstractResult end

"""
    GenericResult

A summary result struct returned by [`optimize`](@ref), including following fields:
 - `optimizer`: the optimization algorithm used
 - `initial_x`: the initial primal solution
 - `minimizer`: the best solution found
 - `minimum`: the optimal value
 - `iter`: number of inner iterations run
 - `maxiter_reached`: true if the algorithm stopped due to reaching the maximum number of iterations
 - `tol`: an instance of [`Tolerance`](@ref) that specifies the convergence tolerances
 - `convstate`: an instance of [`ConvergenceCriteria`](@ref) that summarizes the convergence state of the best solution found
 - `fcalls`: the number of times the objective and constraint functions were called during the optimization
"""
@params mutable struct GenericResult <: AbstractResult
    optimizer
    initial_x
    minimizer
    minimum::Real
    iter::Int
    maxiter_reached::Bool
    tol::Tolerance
    convstate::Union{ConvergenceState, Nothing}
    fcalls::Int
end

abstract type AbstractModel end

"""
```
optimize(
    model::AbstractModel,
    optimizer::AbstractOptimizer = MMA02(),
    x0::AbstractVector;
    options,
    convcriteria::ConvergenceCriteria = KKTCriteria(),
    plot_trace::Bool = false,
    callback::Function = plot_trace ? LazyPlottingCallback() : NoCallback(),
    kwargs...,
)
```

Optimizes `model` using the algorithm `optimizer`, e.g. an instance of [`MMA87`](@ref) or [`MMA02`](@ref). `x0` is the initial solution. The keyword arguments are:
 - `options`: used to set the optimization options. It is an instance of [`MMAOptions`](@ref) for [`MMA87`](@ref) and [`MMA02`](@ref).
 - `convcriteria`: an instance of [`ConvergenceCriteria`](@ref) that specifies the convergence criteria of the MMA algorithm.
 - `plot_trace`: a Boolean that if true specifies the callback to be an instance of [`PlottingCallback`](@ref) and plots a live trace of the last 50 solutions.
 - `callback`: a function that is called on `solution` in every iteration of the algorithm. This can be used to store information about the optimization process.

 The details of the MMA optimization algorithms can be found in the original [1987 MMA paper](https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620240207) and the [2002 paper](https://epubs.siam.org/doi/abs/10.1137/S1052623499362822).
"""
function optimize(model::AbstractModel, optimizer::AbstractOptimizer, x0, args...; kwargs...)
    _model, _x0, unflatten = tovecmodel(model, x0)
    r = optimize(_model, optimizer, _x0, args...; kwargs...)
    return @set r.minimizer = unflatten(r.minimizer)
end


text:
function train_net_tracked!(
    net,
    x_data,
    y_data,
    w,
    batch_size,
    n_batches,
    opt;
    show_progress::Bool = true)

    isempty(x_data) && return nothing
    input_size = length(first(x_data))
    output_size = length(first(y_data))

    loss_hist = Vector{Float64}(undef, n_batches)
    X_tot = Matrix{Float32}(undef, input_size, length(x_data))
    Y_tot = Matrix{Float32}(undef, output_size, length(y_data))
    W_tot = copy(w)
    fillmat!(X_tot, x_data, 1:length(w))
    fillmat!(Y_tot, y_data, 1:length(w))

    X = Matrix{Float32}(undef, input_size, batch_size)
    Y = Matrix{Float32}(undef, output_size, batch_size)
    W = Vector{Float32}(undef, batch_size)
    sample_idxs = Vector{Int}(undef, batch_size)
    idxs = 1:length(w)

    Loss = NetLoss(net, X, Y, W)
    Loss_tot = NetLoss(net, X_tot, Y_tot, W_tot)
    p = Flux.params(net)

    prog = Progress(n_batches; enabled=show_progress)
    for i in 1:n_batches
        loss_hist[i] = Loss_tot()

        rand!(sample_idxs, idxs)
        fillmat!(X::Matrix{Float32}, x_data, sample_idxs)
        fillmat!(Y::Matrix{Float32}, y_data, sample_idxs)
        copyto!(W::Vector{Float32}, @view w[sample_idxs])

        gs = gradient(Loss, p)

        Flux.update!(opt, p::Flux.Params, gs)
        next!(prog)
    end

    return loss_hist
end

function train_net_tracked!(net,Mv::AdvantageMemory,bs,n_b,opt;show_progress::Bool = true)
    return train_net_tracked!(net, Mv.I, Mv.r, Mv.t, bs, n_b, opt; show_progress=show_progress)
end

function train_net_tracked!(net,Mπ::StrategyMemory,bs,n_b,opt;show_progress::Bool = true)
    return train_net_tracked!(net, Mπ.I, Mπ.σ, Mπ.t, bs, n_b, opt; show_progress=show_progress)
end

mutable struct LossCache
    txx::Float32
    tx::Vector{Float32}
    t::Float32
end

LossCache(r::Vector{Float32}) = LossCache(0.0f0, zero(r), 0.0f0)

"""
Lowest possible WMSE
"""
function lower_limit_loss(X::Vector{INFO},Y,W) where INFO

    d = Dict{INFO,LossCache}()

    for (k,t,r) in zip(X,W,Y)
        lc = get!(d, k) do # txxsum, txsum, tsum
            LossCache(r)
        end
        lc.txx += t*sum(abs2, r)
        lc.tx .+= t .* r
        lc.t += t
    end

    l = 0.0f0
    for lc in values(d)
        l += lc.txx - sum(abs2, lc.tx)/lc.t
    end
    return l / length(X)
end

lower_limit_loss(Mπ::StrategyMemory) = lower_limit_loss(Mπ.I, Mπ.σ, Mπ.t)
lower_limit_loss(Mv::AdvantageMemory) = lower_limit_loss(Mv.I, Mv.r, Mv.t)

function optimality_distance(net, x_data, y_data, w)
    isempty(x_data) && return 0.0
    input_size = length(first(x_data))
    output_size = length(first(y_data))

    X_tot = Matrix{Float32}(undef, input_size, length(x_data))
    Y_tot = Matrix{Float32}(undef, output_size, length(y_data))
    W_tot = copy(w)

    fillmat!(X_tot, x_data, 1:length(w))
    fillmat!(Y_tot, y_data, 1:length(w))

    l = wmse(net(X_tot), Y_tot, W_tot)
    l_min = lower_limit_loss(x_data, y_data, w)

    if iszero(l_min)
        return l
    else
        return (l - l_min) / l_min
    end
end

optimality_distance(net, Mπ::StrategyMemory)  = optimality_distance(net, Mπ.I, Mπ.σ, Mπ.t)
optimality_distance(net, Mv::AdvantageMemory) = optimality_distance(net, Mv.I, Mv.r, Mv.t)


text:
@recipe function f(x::Vector{RiskParameters{T}}) where T <: EpidemicModel
  xguide --> "Iteration"
  yguide --> "Value"
  legend --> :none
  label --> permutedims(_parameters(x[1]))
  convert(Array{Float64, 2}, x)
end

@recipe function f(iter,
                   x::Vector{RiskParameters{T}}) where T <: EpidemicModel
  xguide --> "Iteration"
  yguide --> "Value"
  legend --> :none
  label --> permutedims(_parameters(x[1]))
  iter, convert(Array{Float64, 2}, x[iter])
end


text:
# This file was generated, do not modify it. # hide
ŷ = predict(tree, rows=test)
@show ŷ[1]

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
module Enrich

    include("enrichr.jl")

    export addGeneList,enrich

end # module

max_stars_repo_name: Fredrikp-ume/Enrich

max_issues_repo_name: Fredrikp-ume/Enrich

max_forks_repo_name: Fredrikp-ume/Enrich


text:
include("gage.jl")
handle = UInt32(393227)

Gage.initialize()
Gage.getsysteminfo(handle)


text:
# Patterns

const OLD_CONVERSION_PATTERN     = "%d{%d-%b %H:%M:%S}:%p:%t:%m%n"
const BASIC_CONVERSION_PATTERN   = "%r %p %c: %m%n"
const DEFAULT_CONVERSION_PATTERN = "%m%n"
const SIMPLE_CONVERSION_PATTERN  = "%p - %m%n"
const TTCC_CONVERSION_PATTERN    = "%r [%t] %-5p %c - %m%n"

""" Pattern format layout

A flexible layout configurable with pattern string. The conversion pattern is closely related to the conversion pattern of the `printf` function in C. A conversion pattern is composed of literal text and format control expressions called *conversion patterns*.

Conversion Pattern |Description
:------------------|:------------
c{precision}|Outputs the name of the logger that published the logging event.
C{precision}|Outputs the fully qualified module name of the caller issuing the
            |logging request. This conversion specifier can be optionally
            |followed by *precision specifier*, that follows the same rules
            |as the logger name converter.
d{pattern}  |Outputs the date of the logging event. The date conversion specifier
            |may be followed by a set of braces containing a date and time
            |pattern string in `strftime` format.
D{pattern}  |Outputs the date of the logging event. The date conversion specifier
            |may be followed by a set of braces containing a date and time
            |pattern string in `Dates.DateFormat` format.
F           |Outputs the file name where the logging request was issued.
K{key}      |Outputs the entries in a `MapMessage`, if one is present in the event.
            |The **K** conversion character can be followed by the key for the map
            |placed between braces, as in **%K{clientNumber}** where `clientNumber`
            |is the key.
l           |Outputs location information of the caller which generated the logging event.
L           |Outputs the line number from where the logging request was issued.
m           |Outputs the application supplied message associated with the logging event.
M           |Outputs the method name where the logging request was issued.
n           |Outputs the platform dependent line separator character or characters.
p           |Outputs the level of the logging event.
r           |Outputs the number of milliseconds elapsed since logger was configured
            |until the creation of the logging event.
t           |Outputs the process identifier that generated the logging event.
u{type}     |Includes either a random or a time-based UUID. `type` parameter specifies
            |type of the created UUID as string: `TIME` for type 1 and `RANDOM` for type 4.
%           |The sequence **%%** outputs a single percent sign.
"""
type PatternLayout <: StringLayout
    pattern::AbstractString
    header::AbstractString
    footer::AbstractString
    #TODO: RegEx replacement of a procuded message
    # replace::Regex
    # replacement::AbstractString

    strftime_format::AbstractString
    date_format::Dates.DateFormat
    conversion::Regex
    initialized::UInt64
    needbacktrace::Bool
    rng::MersenneTwister
    matched::Base.RegexMatchIterator

    function PatternLayout(pattern::AbstractString=DEFAULT_CONVERSION_PATTERN)
        lyt = new()
        lyt.strftime_format = "%Y-%m-%d %H:%M:%S"
        lyt.date_format = Dates.DateFormat("yyyy-mm-ddTHH:MM:SS.sss")
        lyt.conversion = r"%(-?\d+)?(\.\d+)?(c|C|d|D|F|K|l|L|m|M|n|p|r|t|u|%)?({\d+}|{\D+}|{\w+})?"
        lyt.initialized = time_ns()
        lyt.pattern = pattern
        lyt.rng = MersenneTwister(lyt.initialized)

        # match conversions patterns
        lyt.matched = eachmatch(lyt.conversion, lyt.pattern)

        # Check if backtrace is needed
        lyt.needbacktrace = false
        for m in lyt.matched
            lyt.needbacktrace = m.captures[3][1] in Set(Any['l', 'L', 'M', 'F'])
            lyt.needbacktrace && break
        end

        return lyt
    end
end

function PatternLayout(conf::Dict)
    pattern = get(conf, "pattern", DEFAULT_CONVERSION_PATTERN)
    header = get(conf, "header", "")
    footer = get(conf, "footer", "")
    return PatternLayout(pattern)
end

show(io::IO, lyt::PatternLayout) = print(io, "Pattern($(lyt.pattern))")

# Interface implementation

header(lyt::PatternLayout) = UInt8[]

footer(lyt::PatternLayout) = UInt8[]

serialize(lyt::PatternLayout, evnt::Event) = formatpattern(lyt, evnt)

string(lyt::PatternLayout, evnt::Event) = bytestring(serialize(lyt, evnt))

function formatpattern(lyt::PatternLayout, evnt::Event)
    logstring = UInt8[]

    # Find proper frame
    bt = nothing
    if lyt.needbacktrace
        bts = getbacktrace()
        idx = find(bt->bt.func == :log, bts)
        bt = bts[idx[end]+2]
    end

    # process conversion params
    s = 1
    for m in lyt.matched
        append!(logstring, lyt.pattern[s:(m.offset-1)].data)

        # minimum width
        minw = m.captures[1] != nothing ? try parse(Int, m.captures[1]); catch 0 end : 0
        # maximum width
        maxw = m.captures[2] != nothing ? try parse(Int, m.captures[2][2:end]); catch 0 end : 0
        # formating symbol
        sym = m.captures[3][1]
        # formating symbol parameters
        symparam =  m.captures[4] != nothing ? m.captures[4][2:end-1] : ""

        # process formating symbols
        if sym == 'm' # message
            append!(logstring, (evnt |> message |> formatted).data)
        elseif sym == 'n' # newline
            append!(logstring, LOG4JL_LINE_SEPARATOR)
        elseif sym == '%' # %
            push!(logstring, 0x25)
        else
            output = if sym == 'c' # category name (or logger name)
                evnt |> logger
            elseif sym == 'C' # module
                evnt |> fqmn
            elseif sym == 'd' # date: strftime
                tformat = !isempty(symparam) ? symparam : lyt.strftime_format
                Libc.strftime(tformat,timestamp(evnt))
            elseif sym == 'D' # date: Dates.DateFormat
                tformat = !isempty(symparam) ? Dates.DateFormat(symparam) : lyt.date_format
                Dates.format(timestamp(evnt) |> Dates.unix2datetime, tformat)
            elseif sym == 'F' # file
                bt != nothing ? basename(string(bt.file)) : "NA"
            elseif sym == 'l' # func[file:line]
                bt != nothing ? "$(bt.func)[$(basename(string(bt.file))):$(bt.line)]" : "NA"
            elseif sym == 'L' # line
                bt != nothing ? string(bt.line) : "NA"
            elseif sym == 'M' # function
                bt != nothing ? string(bt.func) : "NA"
            elseif sym == 'p' # level
                rpad(evnt |> level |> string, 5, " ")
            elseif sym == 'r' # time elapsed (milliseconds)
                round(Int, (time_ns()-lyt.initialized)/10e6) |> string
            elseif sym == 't' # thread or PID
                getpid() |> string
            elseif sym == 'u' # UUID
                (symparam == "RANDOM" ? Base.Random.uuid4(lyt.rng) : Base.Random.uuid1(lyt.rng)) |> string
            else
                ""
            end

            # adjust output
            lout = length(output)
            if lout > maxw && maxw != 0
                output = output[(lout-maxw+1):end]
                lout = maxw
            end
            if lout < abs(minw) && minw != 0
                output = minw > 0 ? lpad(output, minw, ' ') : rpad(output, -minw, ' ')
            end
            append!(logstring, output.data)
        end
        s = m.offset+length(m.match)
    end
    if s < length(lyt.pattern)
        append!(logstring, lyt.pattern[s:end].data)
    end
    return logstring
end


text:
# BGZF Index
# ==========
#
# An index type for BGZFStream.
#
# The details of the internal is specified in
# https://samtools.github.io/hts-specs/SAMv1.pdf.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# Virtual offset
const VOffset = UInt64

# Data chunk
const Chunk = Tuple{VOffset,VOffset}

# Binning index
const BinIndex = Dict{UInt32,Vector{Chunk}}

# Linear index
const LinearIndex = Vector{VOffset}

# Metadata providing a summary of the number of mappend/unmapped reads.
struct PseudoBin
    # file range of unmapped reads
    unmapped::Chunk

    # number of mapped read segments
    n_mapped::Int64

    # number of unmapped read segments
    n_unmapped::Int64
end

# Index for BGZFStream; used in BAI and Tabix index.
struct BGZFIndex
    # indexes of contigs (chromosomes)
    data::Vector{Tuple{BinIndex,LinearIndex,Nullable{PseudoBin}}}
end

# 16Kbp
const LinearWindowSize = 16 * 1024

# Find chunks overlapping with `(seqid, interval)` in `index`.
function overlapchunks(index::BGZFIndex, seqid::Integer, interval::UnitRange{<:Integer})
    if !(1 ≤ seqid ≤ endof(index.data))
        throw(ArgumentError("sequence id $(seqid) is out of range"))
    end

    if isempty(interval)
        return Chunk[]
    end

    binindex, linindex, pbin = index.data[seqid]
    bins = reg2bins(interval, 14, 5)
    ret = Chunk[]
    idx = cld(first(interval), LinearWindowSize)
    if endof(linindex) ≥ idx
        # `linindex` may be empty for contigs with no records
        offset = linindex[idx]
        for bin in bins
            if haskey(binindex, bin)
                for chunk in binindex[bin]
                    if chunk[2] > offset
                        push!(ret, chunk)
                    end
                end
            end
        end
    end

    # tidy up the list of chunks
    sort!(ret)
    reduce!(ret)

    return ret
end

# Merge chunks so as to minimize the number of seek operations.
function reduce!(chunks)
    @assert issorted(chunks)
    # NOTE: the maximum size of a BGZF block is 64KiB
    merge_threshold = 64 * 1024 * 2
    i = 1
    while i < endof(chunks)
        chunk = chunks[i]
        next = chunks[i+1]
        if chunk[2] > next[1] || next[1] - chunk[2] ≤ merge_threshold
            # merge overlapping or close chunks
            chunks[i] = (chunk[1], max(chunk[2], next[2]))
            deleteat!(chunks, i + 1)
            continue
        end
        i += 1
    end
    return chunks
end

# Read `n_refs` BAI/Tabix-compatible indexes from `input`.
function read_bgzfindex(input, n_refs)
    indexes = Tuple{BinIndex,LinearIndex,Nullable{PseudoBin}}[]
    for _ in 1:n_refs
        # load a binning index (and a pseudo bin)
        n_bins = read(input, Int32)
        binindex = BinIndex()
        pbin = Nullable{PseudoBin}()
        for _ in 1:n_bins
            bin = read(input, UInt32)
            n_chunks = read(input, Int32)
            if bin == 37450
                # pseudo-bin
                @assert n_chunks == 2
                chunk_beg = read(input, UInt64)
                chunk_end = read(input, UInt64)
                n_mapped = read(input, UInt64)
                n_unmapped = read(input, UInt64)
                pbin = Nullable(PseudoBin(
                    (chunk_beg, chunk_end),
                    n_mapped,
                    n_unmapped))
            else
                chunks = Chunk[]
                for i in 1:n_chunks
                    chunk_beg = read(input, UInt64)
                    chunk_end = read(input, UInt64)
                    push!(chunks, (chunk_beg, chunk_end))
                end
                binindex[bin] = chunks
            end
        end

        # load a linear index
        n_intvs = read(input, Int32)
        linindex = LinearIndex()
        for _ in 1:n_intvs
            push!(linindex, read(input, UInt64))
        end

        push!(indexes, (binindex, linindex, pbin))
    end
    return BGZFIndex(indexes)
end


text:
############################################################################################
### get number of parameters
############################################################################################
"""
    n_par(sem_fit::SemFit)
    n_par(model::AbstractSemSingle)
    n_par(model::SemEnsemble)
    n_par(identifier::Dict)

Return the number of parameters.
"""
function n_par end

n_par(fit::SemFit) = n_par(fit.model)

n_par(model::AbstractSemSingle) = n_par(model.imply)

n_par(model::SemEnsemble) = n_par(model.identifier)

n_par(identifier::Dict) = length(identifier)


text:
# Abstract type for Polygonal Meshes
abstract type AbstractPolygonalMesh end
abstract type AbstractCell{dim,V,F} end

struct Node{dim,T}
    x::Vec{dim, T}
end

Node(x::NTuple{dim,T}) where {dim,T} = Node(Vec{dim,T}(x))

"""
get_coords(node::Node) = node.x
get coordinates of a node
"""
@inline get_coordinates(node::Node) = node.x

struct Cell{dim, N, M}
    nodes::NTuple{N,Int}
    faces::NTuple{M,Int}
end

#Common cell types
#Common cell types
const TriangleCell = Cell{2,3,3}
@inline get_cell_name(::TriangleCell) = "Triangle"
@inline reference_edge_nodes(::Type{TriangleCell}) = ((2,3),(3,1),(1,2))

const RectangleCell = Cell{2,4,4}
@inline get_cell_name(::RectangleCell) = "Rectangle"
@inline reference_edge_nodes(::Type{RectangleCell}) = ((1,2),(2,3),(3,4),(4,1))

@inline getnfaces(cell::Cell{dim,N,M}) where {dim,N,M} = M
function topology_elements(cell::Cell{2},element::Int)
    if element == 0
        return cell.nodes
    elseif element == 1
        return cell.faces
    else
        throw("Topology element of order $element not available for cell type")
    end
end

struct PolygonalMesh{dim,N,M,K,T} <: AbstractPolygonalMesh
    cells::Vector{Cell{dim,N,M}}
    nodes::Vector{Node{dim,T}}
    faces::Matrix{Int}                   # v1 v2 cell1 cell2/0
    facesets::Dict{String,Set{Int}}
    nodesets::Dict{String,Set{Int}}
end

function face_orientation(mesh::PolygonalMesh{2,3,3}, cell_idx::Int, face_idx::Int)
    k = reference_edge_nodes(RefTetrahedron, Val{2})[face_idx]
    return mesh.cells[cell_idx].nodes[k[2]] > mesh.cells[cell_idx].nodes[k[1]]
end

function get_vertices_matrix(mesh::PolygonalMesh{dim,N,M,K,T}) where {dim,N,M,K,T}
    nodes_m = Matrix{T}(undef,length(mesh.nodes),dim)
    for (k,node) in enumerate(mesh.nodes)
        nodes_m[k,:] = node.x
    end
    nodes_m
end
function getcells_matrix(mesh::PolygonalMesh{dim,N,M,K,T}) where {dim,N,M,K,T}
    cells_m = Matrix{Int}(undef, getncells(mesh), n_faces_per_cell(mesh))
    for k = 1:getncells(mesh)
        @. cells_m[k,:] = mesh.cells[k].nodes
    end
    cells_m
end
@inline n_faces_per_cell(mesh::PolygonalMesh{dim,N,M}) where {dim,N,M} = M
@inline n_nodes_per_cell(mesh::PolygonalMesh{dim,N,M}) where {dim,N,M} = N
@inline getnfaces(mesh::PolygonalMesh) = size(mesh.faces,1)
@inline getnnodes(mesh::PolygonalMesh) = length(mesh.nodes)
@inline getncells(mesh::PolygonalMesh) = length(mesh.cells)
@inline getfaceset(mesh::PolygonalMesh, set::String) = mesh.facesets[set]
@inline getnodeset(mesh::PolygonalMesh, set::String) = mesh.nodesets[set]
@inline getnodesets(mesh::PolygonalMesh) = mesh.nodesets
"""
    getcoordinates(cell, mesh::PolygonalMesh)

Return a vector with the coordinates of the vertices of cell number `cell`.
"""
@inline function get_coordinates(cell::Cell, mesh::PolygonalMesh{dim,N,M,K,T}) where {dim,N,M,K,T}
    coords = Vector{Vec{dim,T}}(undef, N)
    for (i,j) in enumerate(cell.nodes)
        coords[i] = mesh.nodes[j].x
    end
    return coords
end

@inline function get_coordinates!(coords::Vector{Vec{dim,T}}, cell::Cell, mesh::PolygonalMesh{dim,N,M,K,T}) where {dim,N,M,K,T}
    @assert length(coords) == N
    for (i,j) in enumerate(cell.nodes)
        coords[i] = mesh.nodes[j].x
    end
    return coords
end

@inline function get_cell_coordinates(cell_idx::Int, mesh::PolygonalMesh{dim,N,M,K,T}) where {dim,N,M,K,T}
    coords = Vector{Vec{dim,T}}(undef, N)
    for (i,j) in enumerate(mesh.cells[cell_idx].nodes)
        coords[i] = mesh.nodes[j].x
    end
    return coords
end

@inline function get_cell_coordinates!(coords::Vector{Vec{dim,T}}, cell_idx::Int, mesh::PolygonalMesh{dim,N,M,K,T}) where {dim,N,M,K,T}
    @assert length(coords) == N
    for (i,j) in enumerate(mesh.cells[cell_idx].nodes)
        coords[i] = mesh.nodes[j].x
    end
    return coords
end
@inline function get_coordinates(face::Int, mesh::PolygonalMesh{dim,N,M,K,T}) where {dim,N,M,K,T}
    return [node.x for node in nodes(face,mesh)]::Vector{Vec{dim,T}}
end
@inline getcells(mesh::PolygonalMesh) = mesh.cells
@inline getnodes(mesh::PolygonalMesh) = mesh.nodes
@inline ndims(mesh::PolygonalMesh{dim}) where {dim} = dim
@inline nodes(ele::Cell, mesh::PolygonalMesh) = [mesh.nodes[node] for node in ele.nodes]
@inline nodes(face::Int, mesh::PolygonalMesh{dim,N,M,L}) where {dim,N,M,L} = [mesh.nodes[node] for node in mesh.faces[face,1:L]]
@inline cells(face::Int, mesh::PolygonalMesh{dim,N,M,L}) where {dim,N,M,L} = [mesh.cells[ele] for ele in mesh.faces[face,L+1:end]]
@inline getcell(cell::Int,face::Int, mesh::PolygonalMesh{dim,N,M,L}) where {dim,N,M,L} = mesh.faces[face,L+cell]
@inline getnode(node::Int,face::Int, mesh::PolygonalMesh{dim,N,M,L}) where {dim,N,M,L} = mesh.nodes[mesh.faces[face,node]]

function cell_diameter(mesh::PolygonalMesh{dim,N,M,NN,T}, idx::Int) where {dim,N,M,NN,T}
    K = getcells(mesh)[idx]
    h = zero(T)
     for k in K.faces
        mσ = norm(get_coordinates(getnode(2, k, mesh)) - get_coordinates(getnode(1, k, mesh)))
        h = max(h, mσ)
    end
    h
end

_check_setname(dict, name) = haskey(dict, name) && throw(ArgumentError("there already exists a set with the name: $name"))
_warn_emptyset(set) = length(set) == 0 && @warn("no entities added to set")

function addfaceset!(mesh::PolygonalMesh, name::String, faceid::Set{Int})
    _check_setname(mesh.facesets, name)
    faceset = Set(faceid)
    _warn_emptyset(faceset)
    mesh.facesets[name] = faceset
    mesh
end

#Check if cell has its vertices ordered counter-clockwise
# Works for convex 2D polygons
function check_orientation(mesh::PolygonalMesh{2,N,M}, cell_idx::Int) where {N,M}
    return sum(verts[j][1]*verts[mod1(j+1,N)][2]-verts[mod1(j+1,N)][1]*verts[j][2] for j ∈ 1:N) > 0
end

function cell_volume(mesh::PolygonalMesh{2,N,M}, cell_idx::Int) where {N,M}
    cell = mesh.cells[cell_idx]
    verts = get_coordinates(cell,  mesh)
    return 0.5*abs(sum(verts[j][1]*verts[mod1(j+1,N)][2]-verts[mod1(j+1,N)][1]*verts[j][2] for j ∈ 1:N))
end

function cell_centroid(mesh::PolygonalMesh{2,N,M}, cell_idx::Int) where {N,M}
    cell = mesh.cells[cell_idx]
    verts = get_coordinates(cell,  mesh)
    Ve = cell_volume(mesh, cell_idx)
    xc = 1/(6*Ve)*sum((verts[j][1]*verts[mod1(j+1,N)][2]-verts[mod1(j+1,N)][1]*verts[j][2])*(verts[j][1]+verts[mod1(j+1,N)][1]) for j ∈ 1:N)
    yc = 1/(6*Ve)*sum((verts[j][1]*verts[mod1(j+1,N)][2]-verts[mod1(j+1,N)][1]*verts[j][2])*(verts[j][2]+verts[mod1(j+1,N)][2]) for j ∈ 1:N)
    return Vec{2}((xc,yc))
end


text:
using TrainView

maxframes = 400

max_curvature = 1e-3

## Define the reference camera
    cameraposition_reference = [0.0,0.6,-2.165];
    ψθφ_ref = [0.0,-19.1,-7.28] .* (pi/180.0);

    camera_reference = VideoCamera(cameraposition_reference;
        sensor_height = 720,
        sensor_width = 1280,
        focal_length = 5.8e-3,
        pixelspermeter = 1 / 5.5e-6,
        ψθφ = ψθφ_ref
    );

    ## Define the reference camera
        cameraposition_reference = [0.0,0.6,-2.165];
        ψθφ_ref = [0.0,-19.1,-7.28] .* (pi/180.0);

        camera_reference = VideoCamera(cameraposition_reference;
            sensor_height = 720,
            sensor_width = 1280,
            focal_length = 5.8e-3,
            pixelspermeter = 1 / 5.5e-6,
            ψθφ = ψθφ_ref
        );

## Define disturbances of the whole journey

# horizontal and vertical curvatures
    βs = [max_curvature * cos(8n * 2pi / maxframes)  for n in 1:maxframes]
    αs = [max_curvature * cos(4n * 2pi / maxframes)  for n in 1:maxframes]

# displacements relative to the track
    choose_distortions = [:Y,:Z,:ψ,:θ,:φ];
    max_distortions = [[0.05,0.06]; pi/180 .* [2.0,2.0,0.2]]

    Nsmooth = Int(round(maxframes / 20))
    Ncut_high = Int(round(maxframes / 35))
    Ncut_low = Int(round(maxframes / 130))

    distortions = map(max_distortions) do dsmax
        ds = 2 .* rand(maxframes+Nsmooth) .- 1.0
        fftds = fft(ds)
        fftds[Ncut_high:end-Ncut_high] .= 0.0
        fftds[1:Ncut_low] .= 0.0
        fftds[end-Ncut_low:end] .= 0.0
        ds = real.(ifft(fftds))
        ds = (dsmax / maximum(ds)) .* ds;
    end;

    distortions = permutedims(hcat(distortions...));

    # Define physical track
    trackprop = TrackProperties(track_gauge = 1.435, track_width = 0.11)
    Xspace = LinRange(5.0,20.0,60);

    using Plots
    h = 300;
    w = h * camera.opticalproperties.sensor_width / camera.opticalproperties.sensor_height
    pyplot(size = (w,h))

    # Set the limits of the image
    left_track  = [[x, -trackprop.track_gauge/2.0, 0.0 ] for x = Xspace];
    right_track = [[x, trackprop.track_gauge/2.0, 0.0] for x = Xspace];

    Luvs = camera_image(camera_reference, left_track);
    Ruvs = camera_image(camera_reference, right_track);

# Create a gif
    Lus = [uv[1] for uv in Luvs];
    Lvs = [uv[2] for uv in Luvs];
    Rus = [uv[1] for uv in Ruvs];
    Rvs = [uv[2] for uv in Ruvs];
    ulims = 1.5 .* (min(minimum(Rus),minimum(Lus)),max(maximum(Rus),maximum(Lus)))
    vlims = 1.5 .* (min(minimum(Rvs),minimum(Lvs)),max(maximum(Rvs),maximum(Lvs)))
    # vlims = (vlims[1] * (1 + sign(vlims[1]) * 0.2), vlims[2] * (1 + sign(vlims[2]) * 0.1))

# ns = Int.(round.(LinRange(1,maxframes,140)))[1:35]
#     for i in eachindex(ns)
#         left_track  = [[x, -trackprop.track_gauge/2.0 + βs[ns[i]]*x^2, αs[ns[i]]*x^2] for x = Xspace]
#         right_track = [[x, trackprop.track_gauge/2.0 + βs[ns[i]]*x^2, αs[ns[i]]*x^2] for x = Xspace]
#
#         tracks = TracksAhead(left_track,right_track; trackproperties = trackprop)
#         camera = VideoCamera(camera_reference, Dict(choose_distortions .=> distortions[:,ns[i]]))
#
#         plot(tracks,camera, sleeper_displace = -0.09 * n, xlims = ulims, ylims = vlims, xguide = "", camera_grid = false, yguide = "",axis = false)
#         savefig("simulated-curve-tracks-$i")
#     end

    anim = @animate for n in Int.(round.(LinRange(1,maxframes,140)))
        left_track  = [[x, -trackprop.track_gauge/2.0 + βs[n]*x^2, αs[n]*x^2] for x = Xspace]
        right_track = [[x, trackprop.track_gauge/2.0 + βs[n]*x^2, αs[n]*x^2] for x = Xspace]

        tracks = TracksAhead(left_track,right_track; trackproperties = trackprop)
        camera = VideoCamera(camera_reference, Dict(choose_distortions .=> distortions[:,n]))

        plot(tracks,camera, sleeper_displace = -0.09 * n, xlims = ulims, ylims = vlims, xguide = "", yguide = "",axis = false)
    end

    # gif(anim, "../../images/videos/only-curvature.gif", fps = 7)
    gif(anim, "../../images/videos/simulate-short-trainview.gif", fps = 7)
    # gif(anim, "../../images/videos/simulate-trainview.gif", fps = 7)

# plot scatter from image recognition

    max_pixel_error = 1
    numberofpoints = 60

    # sample more closer to the camera
    len = 60
    mesh_size = LinRange(1e-5,0.2,len);
    Xspace = 1.0 .+ cumsum(cumsum(mesh_size))[1:end-3]
    filter!(x -> x < 40.0, Xspace)

    n = 10
    left_track  = [[x, -trackprop.track_gauge/2.0 + βs[n]*x^2, αs[n]*x^2] for x = Xspace]
    right_track = [[x, trackprop.track_gauge/2.0 + βs[n]*x^2, αs[n]*x^2] for x = Xspace]

    tracks = TracksAhead(left_track,right_track; trackproperties = trackprop)
    camera = VideoCamera(camera_reference, Dict(choose_distortions .=> distortions[:,n]))

    plot(tracks,camera, xlims = ulims, ylims = vlims, xguide = "", yguide = "")

    left_uvs = camera_image(camera, left_track);
    right_uvs = camera_image(camera, right_track);

    # sample only a subset evening across v
    Lvs = [uv[2] for uv in left_uvs];
    Rvs = [uv[2] for uv in right_uvs];

    Lvmax = maximum(Lvs); Lvmin = minimum(Lvs);
    Rvmax = maximum(Rvs); Rvmin = minimum(Rvs);

    vs = reverse(sort(rand(LinRange(Lvmin,Lvmax,5*numberofpoints),numberofpoints)));
    is = [findmin(abs.(v .- Lvs))[2] for v in vs];
    left_uvs = left_uvs[is]

    vs = reverse(sort(rand(LinRange(Rvmin,Rvmax,5*numberofpoints),numberofpoints)));
    is = [findmin(abs.(v .- Rvs))[2] for v in vs];
    right_uvs = right_uvs[is]

    # add error from image recognition
    left_uvs  = [uv +  max_pixel_error .* (rand(2) .- 0.5) for uv in left_uvs]
    right_uvs = [uv +  max_pixel_error .* (rand(2) .- 0.5) for uv in right_uvs]

    Lus = [uv[1] for uv in left_uvs]; Lvs = [uv[2] for uv in left_uvs];
    Rus = [uv[1] for uv in right_uvs]; Rvs = [uv[2] for uv in right_uvs];

    scatter!(Lus,Lvs,colour=:pink, lab="")
    scatter!(Rus,Rvs,colour=:pink, lab="")

    savefig("../../images/pink-image-sampling.pdf")


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

    Standard parameter updates requiring no knowledge about the underlying
    MCMC backend.
        ✔   Metropolis-Hastings updates with random walk proposals
        ✗   Conjugate updates
        ✗   Metropolis-adjusted Langevin algorithm
        ✗   Hamiltionan Monte Carlo

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

#                          HIGHER LEVEL METHODS
#-------------------------------------------------------------------------------

#=
                    these methods MAY be overwritten:
                                                                              =#

"""
    state(ws::GlobalWorkspace, updt::MCMCParamUpdate)

Return a substate of a gloabl state consisting of coordinates that the `updt` is
concerned with.
"""
state(ws::GlobalWorkspace, updt::MCMCParamUpdate) = subidx(state(ws), updt)

state(ws::GlobalWorkspace, updt::MCMCUpdate) = nothing

#                    CUSTOM UPDATES AND LOW LEVEL METHODS
#-------------------------------------------------------------------------------

#=
            These methods MUST be overwritten for each custom
            MCMCParamUpdate
                                                                              =#

"""
    log_transition_density(updt::MCMCParamUpdate, θ, θ°)

Evaluates the log-transition density for an update θ → θ°
"""
function log_transition_density(updt::MCMCParamUpdate, θ, θ°)
    error("log_transition_density not implemented")
end

"""
    proposal!(updt::MCMCParamUpdate, global_ws, ws::LocalWorkspace, step)

Sample a proposal value for the MCMC's local `state` variable.
"""
function proposal!(updt::MCMCParamUpdate, global_ws, ws::LocalWorkspace, step)
    error("proposal! not implemented")
end

"""
    set_parameters!(
        ::Proposal,
        updt::MCMCParamUpdate,
        global_ws::GlobalWorkspace,
        ws::LocalWorkspace,
    )

Set the parameters of the proposal Law to the current value held by the proposal
`state`.
"""
function set_parameters!(
        ::Proposal,
        updt::MCMCParamUpdate,
        global_ws::GlobalWorkspace,
        ws::LocalWorkspace,
    )
    error("set_parameters! for Proposal() is not implemented")
end

"""
    set_parameters!(
        ::Previous,
        updt::MCMCParamUpdate,
        global_ws::GlobalWorkspace,
        ws::LocalWorkspace,
    )

Set the parameters of the currently accepted Law to the current value held by
the accepted `state`.
"""
function set_parameters!(
        ::Previous,
        updt::MCMCParamUpdate,
        global_ws::GlobalWorkspace,
        ws::LocalWorkspace,
    )
    error("set_parameters! for Previous() is not implemented")
end

#=
                    these methods MAY be overwritten:
                                                                              =#

"""
    log_prior(updt::MCMCParamUpdate, θ)

Evaluate log-prior at θ.
"""
log_prior(updt::MCMCParamUpdate, θ) = logpdf(updt.prior, θ)

"""
    coords(updt::MCMCParamUpdate)

Return a list of coordinates that `updt` is concerned with
"""
coords(updt::MCMCParamUpdate) = updt.coords

invcoords(updt::MCMCParamUpdate) = updt.invcoords

"""
    subidx(θ, updt::MCMCParamUpdate)

Return a view of subset of `θ` that `updt` is concerned with
"""
subidx(θ, updt::MCMCParamUpdate) = view(θ, coords(updt))

"""
    compute_gradients_and_momenta!(
        updt::MCMCParamUpdate, ws::LocalWorkspace, ::Any
    )

To be overwritten for methods that are gradient or momentum based.
"""
function compute_gradients_and_momenta!(
        updt::MCMCParamUpdate, ws::LocalWorkspace, ::Any
    )
    nothing
end

#=
                                RandomWalkUpdate
                                                                              =#

"""
    struct RandomWalkUpdate{TRW,TA,K,TP} <: MCMCParamUpdate
        rw::TRW
        adpt::TA
        coords::K
        prior::TP
    end

Definition of an MCMC update that uses Metropolis-Hastings algorithm with random
walk proposals. `rw` is the random walk sampler, `adpt` is the struct
responsible for adaptation of hyperparameters, `coords` lists the
coordinate indeces of the global `state` that the given update operates on and
`prior` is a prior for this update step (i.e. for a given subset of parameter
vector).

    function RandomWalkUpdate(
        rw::TRW,
        idx_of_global::K;
        prior::TP=ImproperPrior(),
        adpt::TA=NoAdaptation()
    ) where {TRW<:RandomWalk,K,TA,TP}

Base constructor.
"""
struct RandomWalkUpdate{TRW,TA,K,TP} <: MCMCParamUpdate
    rw::TRW
    adpt::TA
    coords::K
    invcoords::Dict{Int64,Int64}
    prior::TP

    function RandomWalkUpdate(
            rw::TRW,
            idx_of_global::K;
            prior::TP=ImproperPrior(),
            adpt::TA=NoAdaptation()
        ) where {TRW<:RandomWalk,K,TA,TP}
        invcoords = Dict{Int64,Int64}()
        for (i,coord) in enumerate(idx_of_global)
            invcoords[coord] = i
        end

        new{TRW,TA,K,TP}(rw, adpt, idx_of_global, invcoords, prior)
    end
end

function readjust!(updt::RandomWalkUpdate, adpt, mcmc_iter)
    readjust!(updt.rw, adpt, mcmc_iter)
end

log_transition_density(updt::RandomWalkUpdate, θ, θ°) = logpdf(updt.rw, θ, θ°)

function proposal!(updt::RandomWalkUpdate, global_ws, ws::LocalWorkspace, step)
    rand!(updt.rw, state(ws), state°(ws))
    while (logpdf(updt.prior, ws.sub_ws°.state) === -Inf)
        rand!(updt.rw, state(ws), state°(ws))
    end
end

function set_parameters!(
        ::Proposal,
        updt::RandomWalkUpdate,
        global_ws::GlobalWorkspace,
        ws::LocalWorkspace,
    )
    set_parameters!(global_ws.P°, coords(updt), state°(ws))
end

function set_parameters!(
        ::Previous,
        updt::RandomWalkUpdate,
        global_ws::GlobalWorkspace,
        ws::LocalWorkspace,
    )
    set_parameters!(global_ws.P, coords(updt), subidx(state(global_ws), updt))
end

#TODO
struct MALAUpdate <: MCMCGradientBasedUpdate
end

#TODO
struct HamiltonianMCUpdate <: MCMCGradientBasedUpdate
end

#TODO
#struct ConjugateGaussianUpdate <: MCMCConjugateUpdate
#end