text:
module sNMF_Exact

include("../Misc.jl")
using .Misc

using Distributions, SpecialFunctions, Einsum, Combinatorics
import Base.Iterators: product

export generate, log_marginal

function generate(I::Ƶ, R::Ƶ, T::Ƶ; γ::ℜ=0.1) where {ℜ<:Real, Ƶ<:Int}
    a::ℜ = I*I*γ
    α_R, α_IR = fill(a/R,R), fill(2.0*a/(I*R),I,R)

    θ_R = rand(Dirichlet(α_R))
    θ_IR = reshape([θ_ir for r=1:R for θ_ir=rand(Dirichlet(α_IR[:,r]))],I,R)

    θ_I₁I₂ = Array{ℜ}(undef,I,I)
    @einsum θ_I₁I₂[i₁,i₂] = θ_R[r] * θ_IR[i₁,r] * θ_IR[i₂,r]
    θ_flat = reshape(θ_I₁I₂, I*I)

    X::Array{ℜ} = reshape(rand(Multinomial(T,θ_flat)),I,I)

    return X, (θ_R, θ_IR)
end

function log_marginal(X::Array{ℜ,2}, R::Ƶ; μ::ℜ=1.0, γ::ℜ=0.1) where {ℜ<:Real, Ƶ<:Int}
    I::Ƶ = size(X,1)
    X_flat = reshape(Int.(X),I*I)

    a::ℜ, b::ℜ = I*I*γ, γ/μ
    α_R, α_IR = fill(a/R,R), fill(2.0*a/(I*R),I,R)

    s, S_R, S_IR, S⁺::ℜ = zeros(ℜ,R), zeros(ℜ,R), zeros(ℜ,I,R), sum(X)
    log_PX, log_PS = -Inf, -Inf
    i₁::Ƶ, i₂::Ƶ = 0, 0

    log_C = a*log(b) - (a + S⁺)*log(b + 1.0)
    log_C += sum(lgamma.(2.0 .* α_R)) - sum(lgamma.(α_IR)) - sum(lgamma.(α_R))

    if R == 1
        S_R .= a + S⁺
        S_IR[:,1] = sum(X,dims=1)[1,:] .+ sum(X,dims=2)[:,1] .+ α_IR[:,1]
        log_PS = log_C - sum(lgamma.(X .+ 1.0))
        log_PS += sum(lgamma.(S_R)) + sum(lgamma.(S_IR)) - sum(lgamma.(2.0 .* S_R))
        return log_PS
    end

    for S_div = product(map(X_i₁i₂ -> combinations(1:(X_i₁i₂+R-1),R-1), X_flat)...)
        S_R .= α_R
        S_IR .= α_IR
        log_PS = log_C
        for (i₁i₂, part) = enumerate(S_div)
            i₁, i₂ = mod(i₁i₂-1,I)+1, div(i₁i₂-1,I)+1
            s[1] = part[1]-1
            for r=2:R-1
                s[r] = part[r] - part[r-1] - 1
            end
            s[R] = X[i₁,i₂] + R - 1 - part[R-1]

            S_R .+= s
            S_IR[i₁,:] .+= s
            S_IR[i₂,:] .+= s
            log_PS -= sum(lgamma.(s .+ 1.0))
        end
        log_PS += sum(lgamma.(S_R)) + sum(lgamma.(S_IR)) - sum(lgamma.(2.0 .* S_R))
        log_PX = logsumexp([log_PX,log_PS])
    end

    return log_PX
end
end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text: const BUILTIN_PKGS = Dict(v=>string(k) for (k, v) in Pkg.Types.stdlib())

ext: jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Tcl_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Tcl")
JLLWrappers.@generate_main_file("Tcl", UUID("9fc9b8dd-f768-5d0d-bb0b-e65961e00eb6"))
end  # module Tcl_jll


text:
# Dense Bitmapped Tries
# =====================

abstract DenseBitmappedTrie{T} <: BitmappedTrie{T}

# Why is the shift value of a DenseLeaf 5 instead of 0, and why does
# the shift value of a DenseNode start at 10?
#
# The PersistentVector implements a "tail" optimization, where it
# inserts appended elements into a tail array until that array is 32
# elements long, only then inserting it into the actual bitmapped
# vector trie. This significantly increases the performance of
# operations that touch the very end of the vector (last, append, pop,
# etc.) because you don't have to traverse the trie.
#
# However, it adds a small amount of complexity to the implementation;
# when you query the trie, you now get back a length-32 array instead
# of the actual element. This is why the DenseLeaf has a shift value
# of 5: it leaves an "extra" 5 bits for the PersistentVector to use to
# index into the array returned from the trie. (This also means that a
# DenseNode has to start at shiftby*2.)

immutable DenseNode{T} <: DenseBitmappedTrie{T}
    arr::Vector{DenseBitmappedTrie{T}}
    shift::Int
    length::Int
    maxlength::Int
end

immutable DenseLeaf{T} <: DenseBitmappedTrie{T}
    arr::Vector{T}

    DenseLeaf(arr::Vector) = new(arr)
    DenseLeaf() = new(T[])
end

arrayof(    node::DenseNode) = node.arr
shift(      node::DenseNode) = node.shift
maxlength(  node::DenseNode) = node.maxlength
Base.length(node::DenseNode) = node.length

arrayof(    leaf::DenseLeaf) = leaf.arr
shift(          ::DenseLeaf) = shiftby
maxlength(  leaf::DenseLeaf) = trielen
Base.length(leaf::DenseLeaf) = length(arrayof(leaf))

function promoted{T}(node::DenseBitmappedTrie{T})
    DenseNode{T}(DenseBitmappedTrie{T}[node],
                 shift(node) + shiftby,
                 length(node),
                 maxlength(node) * trielen)
end

function demoted{T}(node::DenseNode{T})
    if shift(node) == shiftby * 2
        DenseLeaf{T}(T[])
    else
        DenseNode{T}(DenseBitmappedTrie{T}[],
                     shift(node) - shiftby,
                     0,
                     round(Int, maxlength(node) / trielen))
    end
end

function witharr{T}(node::DenseNode{T}, arr::Array, lenshift::Int=0)
    DenseNode{T}(arr, shift(node), length(node) + lenshift, maxlength(node))
end
witharr{T}(leaf::DenseLeaf{T}, arr::Array) = DenseLeaf{T}(arr)

function append(leaf::DenseLeaf, el)
    if length(leaf) < maxlength(leaf)
        newarr = copy_to_len(arrayof(leaf), 1 + length(leaf))
        newarr[end] = el
        witharr(leaf, newarr)
    else
        append(promoted(leaf), el)
    end
end
function append{T}(node::DenseNode{T}, el)
    if length(node) == 0
        child = append(demoted(node), el)
        witharr(node, DenseBitmappedTrie{T}[child], 1)
    elseif length(node) < maxlength(node)
        if length(arrayof(node)[end]) == maxlength(arrayof(node)[end])
            newarr = copy_to_len(arrayof(node), 1 + length(arrayof(node)))
            newarr[end] = append(demoted(node), el)
            witharr(node, newarr, 1)
        else
            newarr = arrayof(node)[:]
            newarr[end] = append(newarr[end], el)
            witharr(node, newarr, 1)
        end
    else
        append(promoted(node), el)
    end
end
push(leaf::DenseLeaf, el) = append(leaf, el)
push(node::DenseNode, el) = append(node, el)

Base.getindex(leaf::DenseLeaf, i::Int) = arrayof(leaf)[mask(leaf, i)]
Base.getindex(node::DenseNode, i::Int) = arrayof(node)[mask(node, i)][i]

function assoc{T}(leaf::DenseLeaf{T}, i::Int, el)
    newarr = arrayof(leaf)[:]
    newarr[mask(leaf, i)] = el
    DenseLeaf{T}(newarr)
end
function assoc(node::DenseNode, i::Int, el)
    newarr = arrayof(node)[:]
    idx = mask(node, i)
    newarr[idx] = assoc(newarr[idx], i, el)
    witharr(node, newarr)
end

peek(bt::DenseBitmappedTrie) = bt[end]

# Pop is usually destructive, but that doesn't make sense for an immutable
# structure, so `pop` is defined to return a Trie without its last
# element. Use `peek` to access the last element.
#
pop(leaf::DenseLeaf) = witharr(leaf, arrayof(leaf)[1:end-1])
function pop(node::DenseNode)
    newarr = arrayof(node)[:]
    newarr[end] = pop(newarr[end])
    witharr(node, newarr, -1)
end

# Persistent Vectors
# ==================

immutable PersistentVector{T} <: AbstractArray{T,1}
    trie::DenseBitmappedTrie{Vector{T}}
    tail::Vector{T}
    length::Int

    PersistentVector(trie::DenseBitmappedTrie, tail, length::Int) =
        new(trie, tail, length)
    PersistentVector() = new(DenseLeaf{Vector{T}}(), T[], 0)
    function PersistentVector(arr)
        if length(arr) <= trielen
            new(DenseLeaf{Vector{T}}(), arr, length(arr))
        else
            append(new(DenseLeaf{Vector{T}}(), T[], 0), arr)
        end
    end
end
PersistentVector() = PersistentVector{Any}()
PersistentVector(itr) = PersistentVector{eltype(itr)}(itr)

mask(i::Int) = ((i - 1) & (trielen - 1)) + 1

function boundscheck!(v::PersistentVector, i::Int)
    0 < i <= v.length || error(BoundsError(), " :: Index $i out of bounds ($(v.length))")
end

Base.size(   v::PersistentVector) = v.length
Base.length( v::PersistentVector) = v.length
Base.isempty(v::PersistentVector) = length(v) == 0
Base.endof(  v::PersistentVector) = length(v)

Base.isequal(v1::PersistentVector, v2::PersistentVector) =
    isequal(v1.tail, v2.tail) && isequal(v1.trie, v2.trie)
==(v1::PersistentVector, v2::PersistentVector) =
    v1.tail == v2.tail && v1.trie == v2.trie

function Base.getindex(v::PersistentVector, i::Int)
    boundscheck!(v, i)
    if i > v.length - length(v.tail)
        v.tail[mask(i)]
    else
        v.trie[i][mask(i)]
    end
end

peek(v::PersistentVector) = v[end]

function push{T}(v::PersistentVector{T}, el)
    if length(v.tail) < trielen
        newtail = copy_to_len(v.tail, 1 + length(v.tail))
        newtail[end] = el
        PersistentVector{T}(v.trie, newtail, 1 + v.length)
    else
        # T[el] will give an error when T is an tuple type in v0.3
        # workaround:
        arr = Array(T, 1)
        arr[1] = convert(T, el)
        PersistentVector{T}(append(v.trie, v.tail), arr, 1 + v.length)
    end
end
append{T}(v::PersistentVector{T}, itr) = foldl(push, v, itr)

function assoc{T}(v::PersistentVector{T}, i::Int, el)
    boundscheck!(v, i)
    if i > v.length - length(v.tail)
        newtail = v.tail[1:end]
        newtail[mask(i)] = el
        PersistentVector{T}(v.trie, newtail, v.length)
    else
        newnode = v.trie[i][1:end]
        newnode[mask(i)] = el
        PersistentVector{T}(assoc(v.trie, i, newnode), v.tail, v.length)
    end
end

function pop{T}(v::PersistentVector{T})
    if isempty(v.tail)
        newtail = peek(v.trie)[1:end-1]
        PersistentVector{T}(pop(v.trie), newtail, v.length - 1)
    else
        newtail = v.tail[1:end-1]
        PersistentVector{T}(v.trie, newtail, v.length - 1)
    end
end

immutable ItrState{T}
    index::Int
    leaf::Vector{T}
end

Base.start{T}(v::PersistentVector{T}) = ItrState(1, v.length <= 32 ? v.tail : v.trie[1])
Base.done{T}(v::PersistentVector{T}, state::ItrState{T}) = state.index > v.length

function Base.next{T}(v::PersistentVector{T}, state::ItrState{T})
    i, leaf = state.index, state.leaf
    m = mask(i)
    value = leaf[m]
    i += 1
    if m == 32
        leaf = i > v.length - length(v.tail) ? v.tail : v.trie[i]
    end
    return value, ItrState(i, leaf)
end

function Base.map{T}(f::Function, pv::PersistentVector{T})
    if length(pv) == 0 return PersistentVector{T}() end

    first = f(pv[1])
    v = PersistentVector{typeof(first)}([first])
    for i = 2:length(pv)
        v = append(v, f(pv[i]))
    end
    v
end

function Base.filter{T}(f::Function, pv::PersistentVector{T})
    v = PersistentVector{T}()
    for el in pv
        if f(el)
            v = append(v, el)
        end
    end
    v
end

function Base.hash{T}(pv::PersistentVector{T})
    h = hash(length(pv))
    for el in pv
        h = Base.hash(el, h)
    end
    @compat UInt(h)
end

function print_elements(io, pv, range)
    for i=range
        print(io, "$(pv[i]), ")
    end
end

function print_vec(io::IO, t, head::AbstractString)
    print(io, "$head[")
    if length(t) < 50
        print_elements(io, t, 1:length(t)-1)
    else
        print_elements(io, t, 1:25)
        print(io, "..., ")
        print_elements(io, t, length(t)-25:length(t)-1)
    end
    if length(t) >= 1
        print(io, "$(t[end])]")
    else
        print(io, "]")
    end
end

@compat Base.show{T}(io::IO, ::MIME"text/plain", pv::PersistentVector{T}) =
    print_vec(io, pv, "Persistent{$T}")


text:
@testset "Momentum Gradient Descent" begin
    # TODO: check the skips and exceptions, maybe it's enough to increase number of iterations?
    skip = ("Rosenbrock", "Extended Powell", "Extended Rosenbrock",
            "Trigonometric", "Penalty Function I", "Beale","Paraboloid Random Matrix")
    run_optim_tests(MomentumGradientDescent(),
                    skip = skip,
                    convergence_exceptions = (("Large Polynomial",1),  ("Himmelblau",1),
                                              ("Fletcher-Powell", 1),("Fletcher-Powell", 2),
                                              ("Powell", 1)),
                    minimum_exceptions = (("Large Polynomial", 1), ("Large Polynomial", 2)),
                    iteration_exceptions = (("Paraboloid Diagonal", 10000),
                                            ("Powell", 10000)),
                    f_increase_exceptions = ("Exponential", "Polynomial",
                                             "Paraboloid Random Matrix", "Hosaki"),
                    show_name = debug_printing)
end


text:
"""
    vgg_block(input_filters, output_filters, depth, batchnorm)

A VGG block of convolution layers ([reference](https://arxiv.org/abs/1409.1556v6)).

# Arguments
- `input_filters`: number of input feature maps
- `output_filters`: number of output feature maps
- `depth`: number of convolution/convolution + batch norm layers
- `batchnorm`: set to `true` to include batch normalization after each convolution
"""
function vgg_block(input_filters, output_filters, depth, batchnorm)
    k = (3, 3)
    p = (1, 1)
    layers = []
    for _ in 1:depth
        push!(layers,
              Conv(k, input_filters => output_filters, batchnorm ? identity : relu; pad=p))
        if batchnorm
            push!(layers, BatchNorm(output_filters, relu))
        end
        input_filters = output_filters
    end
    return Chain(layers...)
end

"""
    vgg_convolutional_layers(config, batchnorm, inchannels)

Create VGG convolution layers ([reference](https://arxiv.org/abs/1409.1556v6)).

# Arguments
- `config`: vector of tuples `(output_channels, num_convolutions)` for each block (see [`Metalhead.vgg_block`](#))
- `batchnorm`: set to `true` to include batch normalization after each convolution
- `inchannels`: number of input channels
"""
function vgg_convolutional_layers(config, batchnorm, inchannels)
    layers = []
    input_filters = inchannels
    for c in config
        push!(layers, vgg_block(input_filters, c..., batchnorm))
        push!(layers, MaxPool((2, 2); stride=2))
        input_filters, _ = c
    end
    return Chain(layers...)
end

"""
    vgg_classifier_layers(imsize, nclasses, fcsize, dropout)

Create VGG classifier (fully connected) layers ([reference](https://arxiv.org/abs/1409.1556v6)).

# Arguments
- `imsize`: tuple `(width, height, channels)` indicating the size after the convolution layers (see [`Metalhead.vgg_convolutional_layers`](#))
- `nclasses`: number of output classes
- `fcsize`: input and output size of the intermediate fully connected layer
- `dropout`: the dropout level between each fully connected layer
"""
function vgg_classifier_layers(imsize, nclasses, fcsize, dropout)
    return Chain(FlattenLayer(),
                 Dense(Int(prod(imsize)), fcsize, relu),
                 Dropout(dropout),
                 Dense(fcsize, fcsize, relu),
                 Dropout(dropout),
                 Dense(fcsize, nclasses))
end

"""
    vgg(imsize; config, inchannels, batchnorm = false, nclasses, fcsize, dropout)

Create a VGG model ([reference](https://arxiv.org/abs/1409.1556v6)).

# Arguments
- `imsize`: input image width and height as a tuple
- `config`: the configuration for the convolution layers
- `inchannels`: number of input channels
- `batchnorm`: set to `true` to use batch normalization after each convolution
- `nclasses`: number of output classes
- `fcsize`: intermediate fully connected layer size
            (see [`Metalhead.vgg_classifier_layers`](#))
- `dropout`: dropout level between fully connected layers
"""
function vgg(imsize; config, inchannels, batchnorm=false, nclasses, fcsize, dropout)
    conv = vgg_convolutional_layers(config, batchnorm, inchannels)
    class = vgg_classifier_layers((7, 7, 512), nclasses, fcsize, dropout)
    return Chain(Chain(conv), class)
end

const VGG_CONV_CONFIG = Dict(:A => [(64, 1), (128, 1), (256, 2), (512, 2), (512, 2)],
                             :B => [(64, 2), (128, 2), (256, 2), (512, 2), (512, 2)],
                             :D => [(64, 2), (128, 2), (256, 3), (512, 3), (512, 3)],
                             :E => [(64, 2), (128, 2), (256, 4), (512, 4), (512, 4)])

const VGG_CONFIG = Dict(11 => :A, 13 => :B, 16 => :D, 19 => :E)

function vgg(name::Symbol; kwargs...)
    assert_name_present_in(name,
                           (:vgg11, :vgg11_bn, :vgg13, :vgg13_bn, :vgg16, :vgg16_bn, :vgg19,
                            :vgg19_bn))
    model = if name == :vgg11
        vgg((224, 224); config=VGG_CONV_CONFIG[VGG_CONFIG[11]], inchannels=3,
            batchnorm=false, nclasses=1000,
            fcsize=4096, dropout=0.5f0)
    elseif name == :vgg11_bn
        vgg((224, 224); config=VGG_CONV_CONFIG[VGG_CONFIG[11]], inchannels=3,
            batchnorm=true, nclasses=1000,
            fcsize=4096, dropout=0.5f0)
    elseif name == :vgg13
        vgg((224, 224); config=VGG_CONV_CONFIG[VGG_CONFIG[13]], inchannels=3,
            batchnorm=false, nclasses=1000,
            fcsize=4096, dropout=0.5f0)
    elseif name == :vgg13_bn
        vgg((224, 224); config=VGG_CONV_CONFIG[VGG_CONFIG[13]], inchannels=3,
            batchnorm=true, nclasses=1000,
            fcsize=4096, dropout=0.5f0)
    elseif name == :vgg16
        vgg((224, 224); config=VGG_CONV_CONFIG[VGG_CONFIG[16]], inchannels=3,
            batchnorm=false, nclasses=1000,
            fcsize=4096, dropout=0.5f0)
    elseif name == :vgg16_bn
        vgg((224, 224); config=VGG_CONV_CONFIG[VGG_CONFIG[16]], inchannels=3,
            batchnorm=true, nclasses=1000,
            fcsize=4096, dropout=0.5f0)
    elseif name == :vgg19
        vgg((224, 224); config=VGG_CONV_CONFIG[VGG_CONFIG[19]], inchannels=3,
            batchnorm=false, nclasses=1000,
            fcsize=4096, dropout=0.5f0)
    elseif name == :vgg19_bn
        vgg((224, 224); config=VGG_CONV_CONFIG[VGG_CONFIG[19]], inchannels=3,
            batchnorm=true, nclasses=1000,
            fcsize=4096, dropout=0.5f0)
    end
    return initialize_model(name, model; kwargs...)
end


text:
"""
$(SIGNATURES)

Find active blocks between an opening token (`otoken`) and a closing token
`ctoken`. These can be nested (e.g. braces). Return the list of such blocks.
If `deactivate` is `true`, all the tokens within the block will be marked as
inactive (for further, separate processing).
"""
function find_ocblocks(tokens::Vector{Token}, ocproto::OCProto;
                       inmath=false)::Tuple{Vector{OCBlock}, Vector{Token}}

    ntokens       = length(tokens)
    active_tokens = ones(Bool, length(tokens))
    ocblocks      = OCBlock[]
    nestable      = ocproto.nest

    # go over active tokens check if there's an opening token
    # if so look for the closing one and push
    for (i, τ) ∈ enumerate(tokens)
        # only consider active and opening tokens
        (active_tokens[i] & (τ.name == ocproto.otok)) || continue
        # if nestable, need to keep track of the balance
        if nestable
            # inbalance ≥ 0, 0 if all opening tokens are closed
            inbalance = 1 # we've seen an opening token
            j = i # index for the closing token
            while !iszero(inbalance) & (j < ntokens)
                j += 1
                inbalance += ocbalance(tokens[j], ocproto)
            end
            if inbalance > 0
                throw(OCBlockError("I found at least one opening  token " *
                                   "'$(ocproto.otok)' that is not closed properly.",
                                   context(τ)))
            end
        else
            # seek forward to find the first closing token
            j = findfirst(cτ -> (cτ.name ∈ ocproto.ctok), tokens[i+1:end])
            # error if no closing token is found
            if isnothing(j)
                throw(OCBlockError("I found the opening token '$(τ.name)' " *
                                   "but not the corresponding closing token.",
                                   context(τ)))
            end
            j += i
        end
        push!(ocblocks, OCBlock(ocproto.name, τ => tokens[j]))

        # remove processed tokens and tokens within blocks except if
        # it's a brace block in a math environment.
        span = ifelse((ocproto.name == :LXB) & inmath, [i, j], i:j)
        active_tokens[span] .= false
    end
    return ocblocks, tokens[active_tokens]
end

"""
$(SIGNATURES)

Helper function to update the inbalance counter when looking for the closing
token of a block with nesting. Adds 1 if the token corresponds to an opening
token, removes 1 if it's a closing token and 0 otherwise.
"""
function ocbalance(τ::Token, ocp::OCProto)::Int
    (τ.name == ocp.otok) && return 1
    (τ.name ∈  ocp.ctok) && return -1
    return 0
end

"""
$(SIGNATURES)

Convenience function to find all ocblocks e.g. such as `MD_OCBLOCKS`. Returns a
vector of vectors of ocblocks.
"""
function find_all_ocblocks(tokens::Vector{Token}, ocplist::Vector{OCProto};
                           inmath=false)
    ocbs_all = Vector{OCBlock}()
    for ocp ∈ ocplist
        ocbs, tokens = find_ocblocks(tokens, ocp; inmath=inmath)
        append!(ocbs_all, ocbs)
    end
    return ocbs_all, tokens
end

"""
$SIGNATURES

Deactivate blocks inside other blocks when the outer block is bound to be
reprocessed. There's effectively two cases:
1. (solved here) when a block is inside a larger 'escape' block.
2. (solved in `check_and_merge_indented_blocks!`) when an indented code block
    is contained inside a larger block.
"""
function deactivate_inner_blocks!(blocks::Vector{OCBlock}, nin=MD_OCB_NO_INNER)
    ranges = Vector{Pair{Int,Int}}()
    isempty(blocks) && return ranges
    # see #444 it's important to ensure the blocks are sorted and they now
    # may not be given that we're finding them in 2 passes.
    sort!(blocks, by=(β->from(β)))
    # CASE 1: block inside a larger escape block.
    #   this can happen if there is a code block in an escape block
    #   (see e.g. #151) or if there's indentation in a math block.
    i       = 1
    nb      = length(blocks)
    active  = ones(Bool, nb)
    heads   = from.(blocks)
    while i <= nb
        cblock = blocks[i]
        # is the block active and is it a "no-inner"?
        if active[i] && cblock.name ∈ nin
            # find all blocks within the span of this block, deactivate them
            chead = heads[i]
            ctail = to(cblock)
            # keep track of the range
            push!(ranges, (chead => ctail))
            # look at all blocks starting after the current one that may
            # be within its span (note that, at worst, we have a few thousand
            # blocks here so that it doesn't need to be super optimised...)
            mask = filter(j -> active[j] && chead < heads[j] < ctail, i+1:nb)
            active[mask] .= false
        end
        i += 1
    end
    deleteat!(blocks, map(!, active))
    return ranges
end

"""
$SIGNATURES

Given ranges found by `deactivate_inner_blocks!` deactivate double brace blocks
that are within those ranges (they will be reprocessed later).
"""
function deactivate_inner_dbb!(dbb, ranges)
    isempty(dbb) && return nothing
    isempty(ranges) && return nothing
    active = ones(Bool, length(dbb))
    for (i, d) in enumerate(dbb)
        hd, td = from(d), to(d)
        r = findfirst(r -> r.first < hd && td < r.second, ranges)
        isnothing(r) || (active[i] = false)
    end
    deleteat!(dbb, map(!, active))
    return nothing
end

"""
$(SIGNATURES)

Merge vectors of blocks by order of appearance of the blocks.
"""
function merge_blocks(lvb::Vector{<:AbstractBlock}...)
    blocks = vcat(lvb...)
    sort!(blocks, by=(β->from(β)))
    return blocks
end

"""
$SIGNATURES

Take a list of token and return those corresponding to special characters or
html entities wrapped in `HTML_SPCH` types (will be left alone by the markdown
conversion and be inserted as is in the HTML).
"""
function find_special_chars(tokens::Vector{Token})
    spch = Vector{HTML_SPCH}()
    isempty(tokens) && return spch
    for τ in tokens
        τ.name == :CHAR_ASTERISK    && push!(spch, HTML_SPCH(τ.ss, "*"))
        τ.name == :CHAR_UNDERSCORE  && push!(spch, HTML_SPCH(τ.ss, "_"))
        τ.name == :CHAR_ATSIGN      && push!(spch, HTML_SPCH(τ.ss, "@"))
        τ.name == :CHAR_BACKSPACE   && push!(spch, HTML_SPCH(τ.ss, "\"))
        τ.name == :CHAR_BACKTICK    && push!(spch, HTML_SPCH(τ.ss, "`"))
        τ.name == :CHAR_LINEBREAK   && push!(spch, HTML_SPCH(τ.ss, "<br/>"))
        τ.name == :CHAR_HTML_ENTITY &&
            validate_html_entity(τ.ss) && push!(spch, HTML_SPCH(τ.ss))
    end
    return spch
end


text:
@testset "NaNMath" begin

    X = ones(3,3)
    @test nanmapslices(mean, X, dims=1) == mapslices(mean, X, dims=1)
    @test nanmapslices(mean, X, dims=(1,2)) == mapslices(mean, X, dims=(1,2))
    X[1,2] = NaN
    @test nanmapslices(mean, X, dims=1) == nanmapslices(mean, X, dims=2)' == [1.0 1.0 1.0]
    @test nanmapslices(mean, X, dims=(1,2)) == ones(1,1)

end


text:
export BEARLearner

mutable struct BEARLearner{
    BA1<:NeuralNetworkApproximator,
    BA2<:NeuralNetworkApproximator,
    BC1<:NeuralNetworkApproximator,
    BC2<:NeuralNetworkApproximator,
    V<:NeuralNetworkApproximator,
    L<:NeuralNetworkApproximator,
    R<:AbstractRNG,
} <: AbstractLearner
    policy::BA1
    target_policy::BA2
    qnetwork1::BC1
    qnetwork2::BC2
    target_qnetwork1::BC1
    target_qnetwork2::BC2
    vae::V
    log_α::L
    γ::Float32
    τ::Float32
    λ::Float32
    ε::Float32
    p::Int
    max_log_α::Float32
    min_log_α::Float32
    sample_num::Int
    kernel_type::Symbol
    mmd_σ::Float32
    batch_size::Int
    update_freq::Int
    update_step::Int
    rng::R
    # Logging
    actor_loss::Float32
    critic_loss::Float32
    mmd_loss
end

"""
    BEARLearner(;kwargs...)

See [Stabilizing Off-Policy Q-Learning via Bootstrapping Error Reduction](https://arxiv.org/abs/1906.00949). This implmentation refers to the [official python code](https://github.com/aviralkumar2907/BEAR).

# Keyword arguments
- `policy`, used to get latent action.
- `target_policy`, similar to `policy`, but used to estimate the target.
- `qnetwork1`, used to get Q-values.
- `qnetwork2`, used to get Q-values.
- `target_qnetwork1`, used to estimate the target Q-values.
- `target_qnetwork2`, used to estimate the target Q-values.
- `vae`, used for sampling action to calculate MMD loss. This
can be implemented using a `VAE` in a `NeuralNetworkApproximator`.
- `log_α`, lagrange multiplier implemented by a `NeuralNetworkApproximator`.
- `γ::Float32 = 0.99f0`, reward discount rate.
- `τ::Float32 = 0.005f0`, the speed at which the target network is updated.
- `λ::Float32 = 0.75f0`, used for Clipped Double Q-learning.
- `ε::Float32 = 0.05f0`, threshold of MMD loss.
- `p::Int = 10`, the number of state-action pairs used when calculating the Q value.
- `max_log_α::Float32 = 10.0f0`, maximum value of `log_α`.
- `min_log_α::Float32 = 10.0f0`, minimum value of `log_α`.
- `sample_num::Int = 10`, the number of sample action to calculate MMD loss.
- `kernel_type::Symbol = :laplacian`, the method of calculating MMD loss. Possible values: :laplacian/:gaussian.
- `mmd_σ::Float32 = 10.0f0`, the parameter used for calculating MMD loss.
- `batch_size::Int = 32`
- `update_freq::Int = 50`, the frequency of updating the `approximator`.
- `update_step::Int = 0`
- `rng = Random.GLOBAL_RNG`
"""
function BEARLearner(;
    policy,
    target_policy,
    qnetwork1,
    qnetwork2,
    target_qnetwork1,
    target_qnetwork2,
    vae,
    log_α,
    γ = 0.99f0,
    τ = 0.005f0,
    λ = 0.75f0,
    ε = 0.05f0,
    p = 10,
    max_log_α = 10.0f0,
    min_log_α = -5.0f0,
    sample_num = 10,
    kernel_type = :laplacian,
    mmd_σ = 10.0f0,
    batch_size = 32,
    update_freq = 50,
    update_step = 0,
    rng = Random.GLOBAL_RNG,
)
    copyto!(policy, target_policy)  # force sync
    copyto!(qnetwork1, target_qnetwork1)  # force sync
    copyto!(qnetwork2, target_qnetwork2)  # force sync
    BEARLearner(
        policy,
        target_policy,
        qnetwork1,
        qnetwork2,
        target_qnetwork1,
        target_qnetwork2,
        vae,
        log_α,
        γ,
        τ,
        λ,
        ε,
        p,
        max_log_α,
        min_log_α,
        sample_num,
        kernel_type,
        mmd_σ,
        batch_size,
        update_freq,
        update_step,
        rng,
        0.0f0,
        0.0f0,
        0.0f0,
    )
end

function (l::BEARLearner)(env)
    s = send_to_device(device(l.policy), state(env))
    s = Flux.unsqueeze(s, ndims(s) + 1)
    s = repeat(s, outer=(1, 1, l.p))
    action = l.policy(l.rng, s; is_sampling=true)
    q_value = l.qnetwork1(vcat(s, action))
    idx = argmax(q_value)
    action[idx]
end

function RLBase.update!(l::BEARLearner, batch::NamedTuple{SARTS})
    D = device(l.qnetwork1)
    s, a, r, t, s′ = (send_to_device(D, batch[x]) for x in SARTS)
    γ, τ, λ = l.γ, l.τ, l.λ

    update_vae!(l, s, a)

    repeat_s′ = repeat(s′, outer=(1, 1, l.p))
    repeat_action′ = l.target_policy(l.rng, repeat_s′, is_sampling=true)

    q′_input = vcat(repeat_s′, repeat_action′)

    q′ = maximum(λ .* min.(l.target_qnetwork1(q′_input), l.target_qnetwork2(q′_input)) + (1 - λ) .* max.(l.target_qnetwork1(q′_input), l.target_qnetwork2(q′_input)), dims=3)

    y = r .+ γ .* (1 .- t) .* vec(q′)

    # Train Q Networks
    a = reshape(a, :, l.batch_size)
    q_input = vcat(s, a)

    q_grad_1 = gradient(Flux.params(l.qnetwork1)) do
        q1 = l.qnetwork1(q_input) |> vec
        loss = mse(q1, y)
        ignore() do
            l.critic_loss = loss
        end
        loss
    end
    update!(l.qnetwork1, q_grad_1)

    q_grad_2 = gradient(Flux.params(l.qnetwork2)) do
        q2 = l.qnetwork2(q_input) |> vec
        loss = mse(q2, y)
        ignore() do
            l.critic_loss += loss
        end
        loss
    end
    update!(l.qnetwork2, q_grad_2)

    repeat_s = repeat(s, outer=(1, 1, l.p))
    repeat_a = repeat(a, outer=(1, 1, l.p))
    repeat_q1 = mean(l.target_qnetwork1(vcat(repeat_s, repeat_a)), dims=(1, 3))
    repeat_q2 = mean(l.target_qnetwork2(vcat(repeat_s, repeat_a)), dims=(1, 3))
    q = vec(min.(repeat_q1, repeat_q2))

    alpha = exp(l.log_α.model[1])

    # Train Policy
    p_grad = gradient(Flux.params(l.policy)) do
        raw_sample_action = decode(l.vae.model, repeat(s, outer=(1, 1, l.sample_num)); is_normalize=false)  # action_dim * batch_size * sample_num
        raw_actor_action = l.policy(repeat(s, outer=(1, 1, l.sample_num)); is_sampling=true) # action_dim * batch_size * sample_num

        mmd_loss = maximum_mean_discrepancy_loss(raw_sample_action, raw_actor_action, l.kernel_type, l.mmd_σ)

        actor_loss = mean(-q .+ alpha .* mmd_loss)
        ignore() do
            l.actor_loss = actor_loss
            l.mmd_loss = mmd_loss
        end
        actor_loss
    end
    update!(l.policy, p_grad)

    # Update lagrange multiplier
    l_grad = gradient(Flux.params(l.log_α)) do
        mean(-q .+ alpha .* (l.mmd_loss .- l.ε))
    end
    update!(l.log_α, l_grad)

    clamp!(l.log_α.model, l.min_log_α, l.max_log_α)

    # polyak averaging
    for (dest, src) in zip(
        Flux.params([l.target_policy, l.target_qnetwork1, l.target_qnetwork2]),
        Flux.params([l.policy, l.qnetwork1, l.qnetwork2]),
    )
        dest .= (1 - τ) .* dest .+ τ .* src
    end
end

function update_vae!(l::BEARLearner, s, a)
    a = reshape(a, :, l.batch_size)
    vae_grad = gradient(Flux.params(l.vae)) do
        recon_loss, kl_loss = vae_loss(l.vae.model, s, a)
        0.5f0 * kl_loss + recon_loss
    end
    update!(l.vae, vae_grad)
end


text:
module FashionMNIST_Tests
using Test
using ColorTypes
using FixedPointNumbers
using MLDatasets
using DataDeps

@testset "Constants" begin
    @test FashionMNIST.TRAINIMAGES == "train-images-idx3-ubyte.gz"
    @test FashionMNIST.TRAINLABELS == "train-labels-idx1-ubyte.gz"
    @test FashionMNIST.TESTIMAGES  == "t10k-images-idx3-ubyte.gz"
    @test FashionMNIST.TESTLABELS  == "t10k-labels-idx1-ubyte.gz"
    @test FashionMNIST.classnames() isa Vector{String}
    @test length(FashionMNIST.classnames()) == 10
    @test length(unique(FashionMNIST.classnames())) == 10

    @test DataDeps.registry["FashionMNIST"] isa DataDeps.DataDep
end

@testset "convert2images" begin
    @test FashionMNIST.convert2image === MNIST.convert2image
end

# NOT executed on CI. only executed locally.
# This involves dataset download etc.
if parse(Bool, get(ENV, "CI", "false"))
    @info "CI detected: skipping dataset download"
else
    data_dir = withenv("DATADEPS_ALWAY_ACCEPT"=>"true") do
        datadep"FashionMNIST"
    end

    const _TRAINIMAGES = joinpath(data_dir, FashionMNIST.TRAINIMAGES)
    const _TRAINLABELS = joinpath(data_dir, FashionMNIST.TRAINLABELS)
    const _TESTIMAGES = joinpath(data_dir, FashionMNIST.TESTIMAGES)
    const _TESTLABELS = joinpath(data_dir, FashionMNIST.TESTLABELS)

    @testset "Images" begin
        # Sanity check that the first trainimage is not the first testimage
        @test FashionMNIST.traintensor(1) != FashionMNIST.testtensor(1)
        # Make sure other integer types work as indicies
        @test FashionMNIST.traintensor(0xBAE) == FashionMNIST.traintensor(2990)

        @testset "Test that traintensor are the train images" begin
            for i = rand(1:60_000, 10)
                @test FashionMNIST.traintensor(i) == reinterpret(N0f8, FashionMNIST.Reader.readimages(_TRAINIMAGES, i))
                @test FashionMNIST.traintensor(Float64, i) == FashionMNIST.Reader.readimages(_TRAINIMAGES, i) ./ 255.0
                @test FashionMNIST.traintensor(UInt8, i) == FashionMNIST.Reader.readimages(_TRAINIMAGES, i)
            end
        end
        @testset "Test that testtensor are the test images" begin
            for i = rand(1:10_000, 10)
                @test FashionMNIST.testtensor(i) == reinterpret(N0f8, FashionMNIST.Reader.readimages(_TESTIMAGES, i))
                @test FashionMNIST.testtensor(Float64, i) == FashionMNIST.Reader.readimages(_TESTIMAGES, i) ./ 255.0
                @test FashionMNIST.testtensor(UInt8, i) == FashionMNIST.Reader.readimages(_TESTIMAGES, i)
            end
        end

        # These tests check if the functions return internaly
        # consistent results for different parameters (e.g. index
        # as int or as vector). That means no matter how you
        # specify an index, you will always get the same result
        # for a specific index.
        for (image_fun, T, nimages) in (
                (FashionMNIST.traintensor, Float32, 60_000),
                (FashionMNIST.traintensor, Float64, 60_000),
                (FashionMNIST.traintensor, N0f8,    60_000),
                (FashionMNIST.traintensor, Int,     60_000),
                (FashionMNIST.traintensor, UInt8,   60_000),
                (FashionMNIST.testtensor,  Float32, 10_000),
                (FashionMNIST.testtensor,  Float64, 10_000),
                (FashionMNIST.testtensor,  N0f8,    10_000),
                (FashionMNIST.testtensor,  Int,     10_000),
                (FashionMNIST.testtensor,  UInt8,   10_000)
            )
            @testset "$image_fun with T=$T" begin
                # whole image set
                A = @inferred image_fun(T)
                @test typeof(A) <: Union{Array{T,3},Base.ReinterpretArray{T,3}}
                @test size(A) == (28,28,nimages)

                @test_throws AssertionError image_fun(T,-1)
                @test_throws AssertionError image_fun(T,0)
                @test_throws AssertionError image_fun(T,nimages+1)

                @testset "load single images" begin
                    # Sample a few random images to compare
                    for i = rand(1:nimages, 10)
                        A_i = @inferred image_fun(T,i)
                        @test typeof(A_i) <: Union{Array{T,2},Base.ReinterpretArray{T,2}}
                        @test size(A_i) == (28,28)
                        @test A_i == A[:,:,i]
                    end
                end

                @testset "load multiple images" begin
                    A_5_10 = @inferred image_fun(T,5:10)
                    @test typeof(A_5_10) <: Union{Array{T,3},Base.ReinterpretArray{T,3}}
                    @test size(A_5_10) == (28,28,6)
                    for i = 1:6
                        @test A_5_10[:,:,i] == A[:,:,i+4]
                    end

                    # also test edge cases `1`, `nimages`
                    indices = [10,3,9,1,nimages]
                    A_vec   = image_fun(T,indices)
                    A_vec_f = image_fun(T,Vector{Int32}(indices))
                    @test typeof(A_vec) <: Union{Array{T,3},Base.ReinterpretArray{T,3}}
                    @test typeof(A_vec_f) <: Union{Array{T,3},Base.ReinterpretArray{T,3}}
                    @test size(A_vec)   == (28,28,5)
                    @test size(A_vec_f) == (28,28,5)
                    for i in 1:5
                        @test A_vec[:,:,i] == A[:,:,indices[i]]
                        @test A_vec[:,:,i] == A_vec_f[:,:,i]
                    end
                end
            end
        end
    end

    @testset "Labels" begin
        # Sanity check that the first trainlabel is not also
        # the first testlabel
        @test FashionMNIST.trainlabels(2) != FashionMNIST.testlabels(2)

        # Check a few hand picked examples. I looked at both the
        # pictures and the native output to make sure these
        # values are correspond to the image at the same index.
        @test FashionMNIST.trainlabels(1) === 9
        @test FashionMNIST.trainlabels(5) === 0
        @test FashionMNIST.trainlabels(60_000) === 5
        @test FashionMNIST.testlabels(1) === 9
        @test FashionMNIST.testlabels(10_000) === 5

        # These tests check if the functions return internaly
        # consistent results for different parameters (e.g. index
        # as int or as vector). That means no matter how you
        # specify an index, you will always get the same result
        # for a specific index.
        # -- However, technically these tests do not check if
        #    these are the actual FashionMNIST labels of that index!
        for (label_fun, nlabels) in
                    ((FashionMNIST.trainlabels, 60_000),
                    (FashionMNIST.testlabels,  10_000))
            @testset "$label_fun" begin
                # whole label set
                A = @inferred label_fun()
                @test typeof(A) <: Vector{Int64}
                @test size(A) == (nlabels,)

                @testset "load single label" begin
                    # Sample a few random labels to compare
                    for i = rand(1:nlabels, 10)
                        A_i = @inferred label_fun(i)
                        @test typeof(A_i) <: Int64
                        @test A_i == A[i]
                    end
                end

                @testset "load multiple labels" begin
                    A_5_10 = @inferred label_fun(5:10)
                    @test typeof(A_5_10) <: Vector{Int64}
                    @test size(A_5_10) == (6,)
                    for i = 1:6
                        @test A_5_10[i] == A[i+4]
                    end

                    # also test edge cases `1`, `nlabels`
                    indices = [10,3,9,1,nlabels]
                    A_vec   = @inferred label_fun(indices)
                    A_vec_f = @inferred label_fun(Vector{Int32}(indices))
                    @test typeof(A_vec)   <: Vector{Int64}
                    @test typeof(A_vec_f) <: Vector{Int64}
                    @test size(A_vec)   == (5,)
                    @test size(A_vec_f) == (5,)
                    for i in 1:5
                        @test A_vec[i] == A[indices[i]]
                        @test A_vec[i] == A_vec_f[i]
                    end
                end
            end
        end
    end

    # Check against the already tested tensor and labels functions
    @testset "Data" begin
        for (data_fun, feature_fun, label_fun, nobs) in
                ((FashionMNIST.traindata, FashionMNIST.traintensor, FashionMNIST.trainlabels, 60_000),
                (FashionMNIST.testdata,  FashionMNIST.testtensor,  FashionMNIST.testlabels,  10_000))
            @testset "check $data_fun against $feature_fun and $label_fun" begin
                data, labels = @inferred data_fun()
                @test data == @inferred feature_fun()
                @test labels == @inferred label_fun()

                for i = rand(1:nobs, 10)
                    d_i, l_i = @inferred data_fun(i)
                    @test d_i == @inferred feature_fun(i)
                    @test l_i == @inferred label_fun(i)
                end

                data, labels = @inferred data_fun(5:10)
                @test data == @inferred feature_fun(5:10)
                @test labels == @inferred label_fun(5:10)

                data, labels = @inferred data_fun(Int, 5:10)
                @test data == @inferred feature_fun(Int, 5:10)
                @test labels == @inferred label_fun(5:10)

                indices = [10,3,9,1,nobs]
                data, labels = @inferred data_fun(indices)
                @test data == @inferred feature_fun(indices)
                @test labels == @inferred label_fun(indices)
            end
        end
    end
end

end


text:
# Here we're going to extend the functionality of Stefan Karpinski's Cards.jl

using Cards, Combinatorics, Random, CSV, DataFrames

const ranks = ['1', '2', '3', '4', '5', '6', '7', '8', '9',
        'T', 'J', 'Q', 'K', 'A']
const hand_rankings = convert(Array{String},
    CSV.read(joinpath(dirname(pathof(Poker)), "..", "data", "hand_ranks.csv"),
    DataFrame, header = false, types=repeat([String], 7462)))
const fresh_deck_cards = [2♣,   2♢,  2♡,  2♠,
                          3♣,   3♢,  3♡,  3♠,
                          4♣,   4♢,  4♡,  4♠,
                          5♣,   5♢,  5♡,  5♠,
                          6♣,   6♢,  6♡,  6♠,
                          7♣,   7♢,  7♡,  7♠,
                          8♣,   8♢,  8♡,  8♠,
                          9♣,   9♢,  9♡,  9♠,
                          10♣, 10♢, 10♡, 10♠,
                          J♣,   J♢,  J♡,  J♠,
                          Q♣,   Q♢,  Q♡,  Q♠,
                          K♣,   K♢,  K♡,  K♠,
                          A♣,   A♢,  A♡,  A♠]
const fresh_deck = Hand(fresh_deck_cards)

"""
    +(x, y)

Adds two hands together, combining their cards into a single hand. Probably more
akin to a union, but duplicate cards are not supported, which helps catch bugs
and ensures hands are exactly the length of the two inputs hands added.

"""
function Base.:+(x::Hand, y::Hand)
    new_hand = Array{Card}(undef, 0)
    for card in x
        push!(new_hand, card)
    end
    for card in y
        push!(new_hand, card)
    end
    return Hand(new_hand)
end

"""
    -(x, y)

Subtracts the cards from one hand from another. Subtraction more closely follows
set theory's subtraction, as cards that are attempted to be subtracted that
aren't in the original hand are ignored. This should probably at least warn the
user for most poker-related use cases.

"""
function Base.:-(x::Hand, y::Hand)
    new_hand = Array{Card}(undef, 0)
    for card in x
        if !(card in y)
            push!(new_hand, card)
        end
    end
    return Hand(new_hand)
end

"""
    deal(deck, n)

Returns both a new hand of n cards, drawn without replacement from the deck,
and the deck without the cards that were drawn

"""
function deal(deck::Hand, n::Integer)
    new_hand = Array{Card}(undef, n)
    order = randperm(length(deck))
    for i = 1:n
        new_hand[i] = deck[order[i]]
    end
    deck -= Hand(new_hand)
    return Hand(new_hand), deck
end

"""
    cards(hand)

Converts a Hand (from Cards.jl) to an array of cards

"""
function cards(hand::Hand)
    cards_to_return = Array{Card}(undef, length(hand))
    for i = 1:length(hand)
        cards_to_return[i] = hand[i]
    end
    return cards_to_return
end

"""
    isless(x, y)

Compares if a card is less than another card based on its position in a fresh
deck. This isn't actually used in Poker.jl right now, but is probably a handy
functionality for any extension of Cards.jl

"""
Base.:isless(x::Card, y::Card) = findfirst(isequal(x), fresh_deck_cards) <
    findfirst(isequal(y), fresh_deck_cards)

"""
    isless(x, y)

Compares if the suit of a card is ranked lower than that of another card. This
isn't actually used in Poker.jl right now, but is probably a handy
functionality for any extension of Cards.jl
"""
Base.:isless(x::Suit, y::Suit) = return x.i < y.i

function isflush(hand::Hand)
    if length(hand) != 5
        return false
    else
        return Cards.suit(hand[1]) == Cards.suit(hand[2]) ==
            Cards.suit(hand[3]) == Cards.suit(hand[4]) == Cards.suit(hand[5])
    end
end

"""
    to_string(hand)

Compares if the suit of a card is ranked lower than that of another card. This
isn't actually used in Poker.jl right now, but is probably a handy
functionality for any extension of Cards.jl
"""
to_string(hand::Hand) = join(map(x -> ranks[x],
    Cards.rank.(reverse(sort(cards(hand))))))

"""
    isless(x, y)

Compares two hands to one another, using Texas Hold'em rules. Ties evaluate as
false. If larger hands are provided than five cards, the maximal combination of
five cards is used for the comparison.

"""
function Base.:isless(x::Hand, y::Hand)
    hand1 = length(x) > 5 ? maximum(Hand.(combinations(x, 5))) : x
    hand2 = length(y) > 5 ? maximum(Hand.(combinations(y, 5))) : y
    ranking1 = isflush(hand1) ? first(findall(isequal(to_string(hand1)),
        hand_rankings)) : last(findall(isequal(to_string(hand1)),
        hand_rankings))
    ranking2 = isflush(hand2) ? first(findall(isequal(to_string(hand2)),
        hand_rankings)) : last(findall(isequal(to_string(hand2)),
        hand_rankings))
    return ranking1 > ranking2
end

"""
    <(x, y)

Syntactic sugar for isless(x, y)

"""
Base.:<(x::Hand, y::Hand) = isless(x, y)

"""
    <(x, y)

Syntactic sugar for isless(y, x)

"""
Base.:>(x::Hand, y::Hand) = isless(y, x)


text:
using Trixi, OrdinaryDiffEq

dg = DGMulti(polydeg=4, element_type = Quad(), approximation_type = Polynomial(),
             surface_integral = SurfaceIntegralWeakForm(FluxLaxFriedrichs()),
             volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))

equations = CompressibleEulerEquations2D(1.4)

"""
A compressible version of the double shear layer initial condition. Adapted from
Brown and Minion (1995). See Section 3, equations (27)-(28) for the original
incompressible version.

- David L. Brown and Michael L. Minion (1995)
  Performance of Under-resolved Two-Dimensional Incompressible Flow Simulations.
  [DOI: 10.1006/jcph.1995.1205](https://doi.org/10.1006/jcph.1995.1205)
"""
function initial_condition_BM_vortex(x, t, equations::CompressibleEulerEquations2D)
  pbar = 9.0 / equations.gamma
  delta = 0.05
  epsilon = 30
  H = (x[2] < 0) ? tanh(epsilon * (x[2] + 0.25)) :  tanh(epsilon * (0.25 - x[2]))
  rho = 1.0
  v1 = H
  v2 = delta * cos(2.0 * pi * x[1])
  p = pbar
  return prim2cons(SVector(rho, v1, v2, p), equations)
end
initial_condition = initial_condition_BM_vortex

mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16),
                   coordinates_min=(-0.5, -0.5), coordinates_max=(0.5, 0.5),
                   periodicity=true)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg)

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval=10)
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)

###############################################################################
# run the simulation

tol = 1.0e-8
sol = solve(ode, RDPK3SpFSAL49(), abstol=tol, reltol=tol,
            save_everystep=false, callback=callbacks);

summary_callback() # print the timer summary


text:
@testset "louder" begin
    m = SquareRooter(4)
    control = Warn(_->true, f="42")

    c = IC.louder(control, by=0)
    IC.train!(m, 1)
    state = @test_logs IC.update!(c, m, -1, 1)
    IC.train!(m, 2)
    state = @test_logs (:warn, r"42") IC.update!(c, m, 0, 2, state)
    @test_logs IC.takedown(c, 1, state)
    @test_logs (:warn, r"A ") IC.takedown(c, 2, state)

    c = IC.louder(control, by=1)
    IC.train!(m, 1)
    state = @test_logs IC.update!(c, m, -2, 1)
    IC.train!(m, 2)
    state = @test_logs (:warn, r"42") IC.update!(c, m, -1, 2, state)
    @test_logs IC.takedown(c, 0, state)
    @test_logs (:warn, r"A ") IC.takedown(c, 1, state)
end

@testset "with_state_do" begin
    m = SquareRooter(4)
    test_controls = [Step(2), InvalidValue(), GL(), Callback()]
    _info = fill((:info, r""), 2*length(test_controls))

    @test_logs(_info...,
               for c in test_controls
               d = IC.with_state_do(c)
               state = IC.update!(c, m, 1, 1)
               @test state == IC.update!(d, m, 1, 1)
               state_c = IC.update!(c, m, 0, 2, state)
               state_d =  IC.update!(d, m, 0, 2, state)
               @test state_c == state_c
               @test IC.done(c, state_c) == IC.done(d, state_d)
               end)

    # integration test:
    m = SquareRooter(4)
    v = Any[]
    f(state) = push!(v, state.new_iterations)
    IC.train!(m, IC.with_state_do(Step(2), f=f), NumberLimit(10))
    @test v == 2:2:20
end

@testset "skip" begin
    m = SquareRooter(4)
    test_controls = [Step(2), InvalidValue(), GL(), Callback()]
    _info = fill((:info, r""), 2)

    @test_logs(_info...,
               for c in test_controls
               s = IC.skip(c, predicate=2)
               @test !s.predicate(1)
               @test s.predicate(2)
               @test !s.predicate(3)
               atomic_state = IC.update!(c, m, 1, 1)
               state = IC.update!(s, m, 1, 1)
               @test state == (n = 1, )
               state =  IC.update!(s, m, 1, 2, state)
               @test state == (atomic_state = atomic_state, n = 2)
               state =  IC.update!(s, m, 1, 3, state)
               @test state == (atomic_state = atomic_state, n = 3)
               atomic_state = IC.update!(c, m, 1, 2, atomic_state)
               state =  IC.update!(s, m, 1, 4, state)
               @test state == (atomic_state = atomic_state, n = 4)
               @test IC.done(c, atomic_state) == IC.done(s, state)
               @test IC.takedown(c, 0, atomic_state) == IC.takedown(s, 0, state)
               end)
end


text:
atypes = (ULogarithmic, Logarithmic, CLogarithmic)

# logarithmic + logarithmic

for (i,A) in enumerate(atypes)
	for (j,B) in enumerate(atypes)
		C = i<j ? B : A
		@eval promote_rule(::Type{$A}, ::Type{$B}) = $C
		@eval promote_rule(::Type{$A}, ::Type{$B{T}}) where {T} = $C{T}
		@eval promote_rule(::Type{$A{T}}, ::Type{$B}) where {T} = $C{T}
		@eval promote_rule(::Type{$A{S}}, ::Type{$B{T}}) where {S,T} =
			$C{promote_type(S,T)}
    end
end

# logarithmic + real

# generated for type-stability
@generated promote_rule(::Type{ULogarithmic}, ::Type{R}) where {R<:Real} =
	try
		:($(typeof(ULogarithmic(one(R)))))
	catch
		:(Union{})
	end

promote_rule(::Type{A}, ::Type{R}) where {A<:AnyLog, R<:Real} =
	promote_type(A, promote_type(ULogarithmic, R))

# logarithmic + complex
# only CLogarithmic+Complex is promoted, so that Logarithmic+Complex is promoted to Complex{Logarithmic}

# generated for type-stability
@generated promote_rule(::Type{CLogarithmic}, ::Type{C}) where {C<:Complex} =
	try
		:($(typeof(CLogarithmic(one(C)))))
	catch
		:(Union{})
	end

promote_rule(::Type{CLogarithmic{T}}, ::Type{C}) where {T, C<:Complex} =
	promote_type(CLogarithmic{T}, promote_type(CLogarithmic, C))


text:
using AbstractPlotting, WGLMakie, JSServe

n = 5
line = Node(rand(n))
function sample()
    line[] = rand(n)
    line[] = line[]
end

reading = @async while true
    cond = Condition()
    Timer(x->notify(cond), 1/30)
    t = @async sample()
    wait(cond)
    wait(t)
    sleep(0)
end

handler(session, request) = lines(line)

app = JSServe.Application(handler, "0.0.0.0", 8000)


text:
function handle_macro(@nospecialize(x), state) end
function handle_macro(x::EXPR, state)
    typof(x) !== MacroCall && return
    if typof(x.args[1]) === MacroName
        state(x.args[1])
        if _points_to_Base_macro(x.args[1], "deprecate", state) && length(x.args) == 3
            if bindingof(x.args[2]) !== nothing
                return
            elseif CSTParser.is_func_call(x.args[2])
                # add deprecated method
                # add deprecated function binding and args in new scope
                CSTParser.setbinding!(x.args[2], x)
                CSTParser.mark_sig_args!(x.args[2])
                s0 = state.scope # store previous scope
                state.scope = Scope(s0, Dict(), nothing, false)
                setscope!(x, state.scope) # tag new scope to generating expression
                state(x.args[2])
                state(x.args[3])
                state.scope = s0
            elseif isidentifier(x.args[2])
                CSTParser.setbinding!(x.args[2], x)
            end
        elseif _points_to_Base_macro(x.args[1], "enum", state)
            for i = 2:length(x.args)
                if !(typof(x.args[i]) === PUNCTUATION)
                    if bindingof(x.args[i]) !== nothing
                        break
                    end
                    if i == 3 && typof(x.args[i]) === CSTParser.Begin
                        for j in 1:length(x.args[i].args[2].args)
                            CSTParser.setbinding!(x.args[i].args[2].args[j], x)
                        end
                        break
                    end
                    CSTParser.setbinding!(x.args[i], x)
                end
            end
        elseif _points_to_Base_macro(x.args[1], "goto", state)
            if length(x.args) == 2 && typof(x.args[2]) === CSTParser.IDENTIFIER
                setref!(x.args[2], NoReference)
            end
        elseif _points_to_Base_macro(x.args[1], "label", state)
            if length(x.args) == 2 && typof(x.args[2]) === CSTParser.IDENTIFIER
                CSTParser.setbinding!(x.args[2])
            end
        # elseif _points_to_Base_macro(x.args[1], "nospecialize", state)
        elseif length(x.args[1].args) == 2 && isidentifier(x.args[1].args[2]) && valof(x.args[1].args[2]) == "nospecialize"
            for i = 2:length(x.args)
                if !(typof(x.args[i]) === PUNCTUATION)
                    if bindingof(x.args[i]) !== nothing
                        break
                    end
                    CSTParser.setbinding!(x.args[i], x)
                end
            end
        end
    end
end

function _points_to_Base_macro(x::EXPR, name, state)
    length(x.args) == 2 && isidentifier(x.args[2]) && valof(x.args[2]) == name && refof(x.args[2]) == getsymbolserver(state.server)["Base"].vals[string("@", name)]
end

function _points_to_arbitrary_macro(x::EXPR, module_name, name, state)
    length(x.args) == 2 && isidentifier(x.args[2]) && valof(x.args[2]) == name && refof(x.args[2]) == getsymbolserver(state.server)[module_name].vals[string("@", name)]
end


text:
module FDTDSolver

import LoopVectorization: @turbo
import Parameters: @with_kw
import RecipesBase: @recipe
import FFTW: fft,fftfreq
import Optim: optimize, ParticleSwarm
import LinearAlgebra: norm
include("EnumeTypes.jl")
include("Grid.jl")
include("Waveforms.jl")
include("Fields.jl")
include("FullAnistropicEngine/FullAnistropicEngine.jl")

end


text:
__precompile__()

module StaticFFTs

using Primes, StaticArrays

export fft

include("construct_fft.jl")
include("interface.jl")

end # module

max_issues_repo_name: chriselrod/StaticFFTs

max_forks_repo_name: chriselrod/StaticFFTs


text:
export SWCutoff, ShiftCutoff, SplineCutoff

# this file is include from Potentials.jl
# i.e. it is part of JuLIP.Potentials
# it implements three types of cutoff potentials.

######################## Stillinger-Weber type cutoff

"""
`cutsw(r, Rc, Lc)`

Implementation of the C^∞ Stillinger-Weber type cut-off potential
    1.0 / ( 1.0 + exp( lc / (Rc-r) ) )

`d_cutinf` implements the first derivative
"""
@inline cutsw(r, Rc, Lc) =
    1.0 ./ ( 1.0 + exp( Lc ./ ( max(Rc-r, 0.0) + 1e-2 ) ) )

"derivative of `cutsw`"
@inline function cutsw_d(r, Rc, Lc)
    t = 1 ./ ( max(Rc-r, 0.0) + 1e-2 )    # a numerically stable (Rc-r)^{-1}
    e = 1.0 ./ (1.0 + exp(Lc * t))         # compute exponential only once
    return - Lc * (1.0 - e) .* e .* t.^2
end

# doc below
@pot type SWCutoff <: PairPotential
    Rc::Float64
    Lc::Float64
    e0::Float64
end

"""
`type SWCutoff`: SW type cut-off potential with C^∞ regularity
     1.0 / ( 1.0 + exp( Lc / (Rc-r) ) )

This is not very optimised: one could speed up `evaluate_d` significantly
by avoiding multiple evaluations.

### Parameters

* `Rc` : cut-off radius
* `Lc` : scale
"""
SWCutoff

evaluate(p::SWCutoff, r) = p.e0 * cutsw(r, p.Rc, p.Lc)
evaluate_d(p::SWCutoff, r) = p.e0 * cutsw_d(r, p.Rc, p.Lc)
cutoff(p::SWCutoff) = p.Rc

# TODO: hack to make TB work (reconsider this)
evaluate(p::SWCutoff, r, R) = evaluate(p, r)
evaluate_d(p::SWCutoff, r, R) = evaluate_d(p, r)

# simplified constructor to ensure compatibility
SWCutoff(Rc, Lc) = SWCutoff(Rc, Lc, 1.0)

# kw-constructor (original SW parameters, pair term)
SWCutoff(; Rc=1.8, Lc=1.0, e0=1.0) = SWCutoff(Rc, Lc, e0)

######################## Shift-Cutoff: one should not use this!
# TODO: why is this commented out? probably uncomment and check.

# """
# `ShiftCutoff` : takes the pair-potential and shifts and truncates it
#     f_cut(r) = (f(r) - f(rcut)) .* (r <= rcut)
#
# Note this is not constructed by multiplying the cutoff potential with the
# actual potential, but it is constructed by
# ```julia
#     ShiftCutoff(Rc, pp)
# ```
#
# **WARNING:** this creates potentials that are not even C^1, which can
# cause all sorts of hell.
# """
# @pot type ShiftCutoff{T <: PairPotential} <: PairPotential
#     pp::T
#     Rc::Float64
#     Jc::Float64
# end
# ShiftCutoff(pp, Rc) = ShiftCutoff(pp, Rc, pp(Rc))
# evaluate(p::ShiftCutoff, r) = (p.pp(r) - p.Jc) .* (r .<= p.Rc)
# evaluate_d(p::ShiftCutoff, r) = (@D p.pp(r)) .* (r .<= p.Rc)
# cutoff(p::SWCutoff) = p.Rc

######################## Quintic Spline cut-off:

function fcut(r, r0, r1)
    s = 1 - (r-r0) / (r1-r0)
    return ((s .>= 1) + (0 .<= s .< 1) .* (6 * s.^5 - 15 * s.^4 + 10 * s.^3) )
end

"Derivative of `fcut`; see documentation of `fcut`."
function fcut_d(r, r0, r1)
    s = 1-(r-r0) / (r1-r0)
    return ( - (30*s.^4 - 60 * s.^3 + 30 * s.^2) / (r1-r0)
             .* (0 .< s .< 1) )
end

@pot type SplineCutoff <: PairPotential
   r0::Float64
   r1::Float64
end

"""
`SplineCutoff` : Piecewise quintic C^{2,1} regular polynomial.

Parameters:

* r0 : inner cut-off radius
* r1 : outer cut-off radius.
"""
SplineCutoff

evaluate(p::SplineCutoff, r) = fcut(r, p.r0, p.r1)
evaluate_d(p::SplineCutoff, r) = fcut_d(r, p.r0, p.r1)
cutoff(p::SplineCutoff) = p.r1
Base.string(p::SplineCutoff) = "SplineCutoff(r0=$(p.r0), r1=$(p.r1))"


text:
# ------------------------------------------------------------------------------------------
# ## Comenzando con Jupyter Notebooks
#
# Un documento de **Jupyter** está compuesto por celdas. Las celdas se delimitan por cajas
# con bordes grises o verdes, dependiendo si se están editando o no.
#
# Las celdas se clasifican en dos tipos: celdas de código y de texto tipo *markdown*.
#
# ### Editando celdas tipo *Markdown*
# Para ver el texto plano de esta celda, solamente debes hacer doble click sobre el texto.
# Automáticamente podrás editar el texto y caerás en cuenta lo sencillo es darle formato.
#
# Aquí encontrarás un [tutorial](https://www.markdowntutorial.com) y [manual de referencias
# rápido](https://en.support.wordpress.com/markdown-quick-reference/) con la información
# básica que vas a necesitar para darle un formato pulcro a tus archivos.
#
# ### Corriendo una celda de código
# Para ejecutar una celda de código, seleccióna la celda y (1) presiona `Shift` + `Enter` o
# (2) haz click en el botón de correr que se encuentra en la barra de herramientas superior:
# ------------------------------------------------------------------------------------------

1 + 1
2 + 2

# ------------------------------------------------------------------------------------------
# Si eres buen observador podrás notar que solamente la última línea de cada celda se
# imprime a pantalla al ejecutarse. Se puede suprimir esta impresión a pantalla con un
# semicolon.
# ------------------------------------------------------------------------------------------

1 + 1
2 + 2;

# ------------------------------------------------------------------------------------------
# ### ¿Cómo buscar información sobre las funciones?
#
# ¡Igualito que en la REPL!
#
# ------------------------------------------------------------------------------------------

?println

# ------------------------------------------------------------------------------------------
# ### ¿Cómo usar comandos del shell?
#
# ¡Lo mismo!
#
# ------------------------------------------------------------------------------------------

;ls

;pwd

# ------------------------------------------------------------------------------------------
# Como ya podrás imaginarte, todos los comandos que hemos usados en la REPL se podrán usar
# acá....
#
# Para una referencia completa sobre lo que se puede hacer en esta interfaz, como siempre,
# ver la [documentación oficial](https://jupyter-
# notebook.readthedocs.io/en/stable/examples/Notebook/Notebook%20Basics.html#)
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
### A Pluto.jl notebook ###
# v0.12.17

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 6d5ac8a4-47dc-11eb-3098-abac2e02369b
begin
using Plots
using PlutoUI
end

# ╔═╡ 7525b260-47dc-11eb-172e-a721065970f4
begin
	r= @bind range html"<input type=range min=5 max=15>"
	p= @bind power Slider(-2:5)

	md"""
	$(r) $(p)

	power= $(range)
	range= $(power)
	"""
end

# ╔═╡ 72e5bc98-47dc-11eb-1473-675b6bb0958c
begin
	f(x)=x.^power
	x=-range:range
	plot(x,f(x))
end

# ╔═╡ 13611b20-47e4-11eb-05c0-bb4c4d7be8d3
@bind dims html"""
<canvas width="200" height="200" style="position: relative"></canvas>

<script>
// 🐸 `currentScript` is the current script tag - we use it to select elements 🐸 //
const canvas = currentScript.closest('pluto-output').querySelector("canvas")
const ctx = canvas.getContext("2d")

var startX = 80
var startY = 40

function onmove(e){
	// 🐸 We send the value back to Julia 🐸 //
	canvas.value = [e.layerX - startX, e.layerY - startY]
	canvas.dispatchEvent(new CustomEvent("input"))

	ctx.fillStyle = '#ffecec'
	ctx.fillRect(0, 0, 200, 200)
	ctx.fillStyle = '#3f3d6d'
	ctx.fillRect(startX, startY, ...canvas.value)
}

canvas.onmousedown = e => {
	startX = e.layerX
	startY = e.layerY
	canvas.onmousemove = onmove
}

canvas.onmouseup = e => {
	canvas.onmousemove = null
}

// Fire a fake mousemoveevent to show something
onmove({layerX: 130, layerY: 160})

</script>
"""

# ╔═╡ Cell order:
# ╠═6d5ac8a4-47dc-11eb-3098-abac2e02369b
# ╠═72e5bc98-47dc-11eb-1473-675b6bb0958c
# ╠═7525b260-47dc-11eb-172e-a721065970f4
# ╠═13611b20-47e4-11eb-05c0-bb4c4d7be8d3

max_stars_repo_name: punkmath/Misc

max_issues_repo_name: punkmath/Misc

max_forks_repo_name: punkmath/Misc


text:
#=
Some notes to myself because I always forget it:
-> Estimate normalizing marginal likelihood term in P(s_1:t | e_1:t) via mean(w_1:t)
-> calculate w_1:t iteratively via w_t-k:t
-> in SIS we calculate normalizing W_t ∝ W_t-1 * w_t-k:t
-> in PF we always calculate normalizing W_t ∝ w_t-k:t, because W_t-1 is constant from resampling step.
-> in adaptive Res. PF, we still need normalizing W_t ∝ W_t-1 * w_t-k:t. In case resampled, W_t is constant, if not need SIS calculation.
-> Sometimes people define unnormalized w_t via w_t-1 * g(..) -> but this refers to w_1:t for each step! Not what we need in PF, as we calculate normalizing term iteratively.
=#

############################################################################################
"""
$(TYPEDEF)
Super type for parameter weighting techniques.
"""
abstract type ParameterWeighting end

"""
$(SIGNATURES)
Weight function that is dispatched on BaytesCore.ParameterWeighting types.

# Examples
```julia
```

"""
function weight!() end

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

Container for weights and normalized particles weights at current time step t.

# Fields
$(TYPEDFIELDS)
"""
struct ParameterWeights
    "Used for log likelihood calculation."
    ℓweights::Vector{Float64}
    "exp(ℓweightsₙ) used for resampling particles and sampling trajectories ~ kept in log space for numerical stability."
    ℓweightsₙ::Vector{Float64}
    "A buffer vector with same size as ℓweights and ℓweightsₙ."
    buffer::Vector{Float64}
    function ParameterWeights(Nparameter::I) where {I<:Integer}
        ArgCheck.@argcheck Nparameter > 0 "Number of particles need to be positive."
        return new(
            fill(log(1.0 / Nparameter), Nparameter),
            fill(log(1.0 / Nparameter), Nparameter),
            fill(1.0 / Nparameter, Nparameter),
        )
    end
end

############################################################################################
#!NOTE: Cannot document functor
# Returns log weights and normalized log weights at time > 1.
function (weights::ParameterWeights)(ℓevidenceₜ::AbstractVector{T}) where {T<:Real}
    weights.ℓweights .= ℓevidenceₜ #Incremental weight
    normalize!(weights)
    return nothing
end

############################################################################################
"""
$(SIGNATURES)
Set weights back to equal.

# Examples
```julia
```

"""
function update!(weights::ParameterWeights)
    weights.ℓweights .= weights.ℓweightsₙ .= log(1.0 / length(weights.buffer))
    return nothing
end

############################################################################################
"""
$(SIGNATURES)
Inplace-Normalize parameter weights, accounting for ℓweightsₙ at previous iteration.

# Examples
```julia
```

"""
function normalize!(weights::ParameterWeights)
    # Wₜ ∝ Wₜ₋₁×wₜ #!NOTE: Relevant if no resampling step performed
    weights.ℓweightsₙ .= (weights.ℓweightsₙ .+ weights.ℓweights)
    weights.ℓweightsₙ .-= logsumexp(weights.ℓweightsₙ)
    return nothing
end

############################################################################################
"""
$(SIGNATURES)
Stable version of computing weighted incremental weight average.

# Examples
```julia
```

"""
function weightedincrement(weights::ParameterWeights)
    return logsumexp(weights.ℓweightsₙ[iter] + weights.ℓweights[iter] for iter in eachindex(weights.ℓweightsₙ))
end

############################################################################################
"""
$(SIGNATURES)
Stable version of computing effectice sample size of particle filter via normalized log weights.

# Examples
```julia
```

"""
function computeESS(weightsₙ::Vector{T}) where {T<:AbstractFloat}
    return min(T(length(weightsₙ)), 1.0 / sum(abs2, weightsₙ))
end
function computeESS(weights::ParameterWeights)
    weights.buffer .= exp.(weights.ℓweightsₙ)
    return min(eltype(weights.buffer)(length(weights.buffer)), 1.0 / sum(abs2, weights.buffer))
end

############################################################################################
"""
$(SIGNATURES)
Draw proposal trajectory index.

# Examples
```julia
```

"""
function draw!(_rng::Random.AbstractRNG, weights::ParameterWeights)
    ## Draw new path
    weights.buffer .= exp.(weights.ℓweightsₙ)
    path = randcat(_rng, weights.buffer)
    ## Return trajectory
    return path
end

############################################################################################
# Export
export
    ParameterWeighting,
    weight!,
    ParameterWeights,
    draw!,
    normalize!,
    weightedincrement,
    computeESS


text:
using RHEOS
using SpecialFunctions
using Test

const tol = (eps(RHEOS.RheoFloat))^(0.125)

include("definitions.jl")
include("base.jl")
include("datagen.jl")
include("IO.jl")
include("processing.jl")

include(joinpath(@__DIR__, "models", "elements.jl"))
include(joinpath(@__DIR__, "models", "maxwell.jl"))
include(joinpath(@__DIR__, "models", "kelvinvoigt.jl"))
include(joinpath(@__DIR__, "models", "zener.jl"))
include(joinpath(@__DIR__, "models", "poynting-thomson.jl"))


text:
println("Testing custom complexities.")
using SymbolicRegression, Test

x1, x2, x3 = Node("x1"), Node("x2"), Node("x3")

# First, test regular complexities:
function make_options(; kw...)
    return Options(; binary_operators=(+, -, *, /, ^), unary_operators=(cos, sin), kw...)
end
options = make_options()
tree = sin((x1 + x2 + x3)^2.3)
@test compute_complexity(tree, options) == 8

options = make_options(; complexity_of_operators=[sin => 3])
@test compute_complexity(tree, options) == 10
options = make_options(; complexity_of_operators=[sin => 3, (+) => 2])
@test compute_complexity(tree, options) == 12

# Real numbers:
options = make_options(; complexity_of_operators=[sin => 3, (+) => 2, (^) => 3.2])
@test compute_complexity(tree, options) == round(Int, 12 + (3.2 - 1))

# Now, test other things, like variables and constants:
options = make_options(;
    complexity_of_operators=[sin => 3, (+) => 2], complexity_of_variables=2
)
@test compute_complexity(tree, options) == 12 + 3 * 1
options = make_options(;
    complexity_of_operators=[sin => 3, (+) => 2],
    complexity_of_variables=2,
    complexity_of_constants=2,
)
@test compute_complexity(tree, options) == 12 + 3 * 1 + 1
options = make_options(;
    complexity_of_operators=[sin => 3, (+) => 2],
    complexity_of_variables=2,
    complexity_of_constants=2.6,
)
@test compute_complexity(tree, options) == 12 + 3 * 1 + 1 + 1

println("Passed.")


text:
"""
    MolecularProperty(molar_mass, p_crit, T_crit, V_crit, acentric_factor = 0.0)

Type that defines the static properties of a molecular species.
"""
struct MolecularProperty{R}
    mw::R
    p_c::R
    T_c::R
    V_c::R
    ω::R
    function MolecularProperty(mw, p_c, T_c, V_c, ω = 0.0)
        @assert p_c >= 0
        @assert mw > 0
        @assert T_c > 0
        @assert V_c > 0
        @assert -1.0 <= ω "Expected -1 ≤ ω, was $ω"
        new{typeof(mw)}(mw, p_c, T_c, V_c, ω)
    end
end

"""
    MolecularProperty("Name")

Convenience constructor that looks up molecular properties from a table in `tabulated_properties`.

The properties are taken from the wonderful MIT-licensed [CoolProp](http://coolprop.org). Please note that
the equations of state included in this module may not be approprioate for all the available fluids, especially for mixtures!

See list of species [at CoolProp website](http://coolprop.org/fluid_properties/PurePseudoPure.html#list-of-fluids).

"""
function MolecularProperty(s::String)
    t = tabulated_properties();
    @assert haskey(t, s)
    return t[s]
end
"""Get the critical pressure for a species."""
@inline critical_pressure(M::MolecularProperty{R}) where R = M.p_c::R
"""Get the critical temperature for a species."""
@inline critical_temperature(M::MolecularProperty{R}) where R = M.T_c::R
"""Get the critical volume for a species."""
@inline critical_volume(M::MolecularProperty{R}) where R = M.V_c::R
"""Get the molar weight for a species."""
@inline molar_weight(M::MolecularProperty{R}) where R = M.mw::R
"""Get the acentric factorfor a species."""
@inline acentric_factor(M::MolecularProperty{R}) where R = M.ω::R

"""
    MultiComponentMixture(properties::NTuple{N, MolecularProperty}; A_ij = nothing, names = ["C1", "C2", ...], name = "UnnamedMixture")

Create a multicomponent mixture with an optional binary interaction coefficient matrix `A_ij`.
"""
struct MultiComponentMixture{R, N}
    name::String
    component_names::Vector{String}
    properties::NTuple{N, MolecularProperty{R}}
    binary_interaction::Union{Matrix{R}, Nothing}
    function MultiComponentMixture(properties; A_ij = nothing, names = ["C$d" for d in 1:length(properties)], name = "UnnamedMixture")
        properties = Tuple(properties)
        realtype = typeof(properties[1].mw)
        n = length(properties)
        if !isnothing(A_ij)
            @assert size(A_ij) == (n, n)
        end
        new{realtype, n}(name, names, properties, A_ij)
    end
end

"""
    number_of_components(mixture)

Return number of components in the `MultiComponentMixture`.
"""
number_of_components(m::MultiComponentMixture{R, N}) where {R, N} = N
function molecular_properties(mixture::MultiComponentMixture{R, N}) where {R, N}
    return mixture.properties::NTuple{N, MolecularProperty{R}}
end
@inline molecular_property(mixture::MultiComponentMixture{R}, i) where R = @inbounds molecular_properties(mixture)[i]::MolecularProperty{R}
@inline reduced_pressure(mixture::MultiComponentMixture, cond, i) = cond.p/critical_pressure(molecular_property(mixture, i))
@inline reduced_temperature(mixture::MultiComponentMixture, cond, i) = cond.T/critical_temperature(molecular_property(mixture, i))
@inline reduced_pT(arg...) = (reduced_pressure(arg...), reduced_temperature(arg...))


text:
# Rank-based fitness assignment
# sp - selective linear presure in [1.0, 2.0]
function ranklinear(sp::Float64)
  @assert 1.0 <= sp <= 2.0 "Selective pressure has to be in range [1.0, 2.0]."
  function rank(fitness::Vector{Float64}, N::Int)
      λ = float(length(fitness))
      idx = sortperm(fitness)
      ranks = zeros(Integer(λ))
      for i in 1:Integer(λ)
          ranks[i] = ( 2.0- sp + 2.0*(sp - 1.0)*(idx[i] - 1.0) / (λ - 1.0) ) /λ
      end
      return sample(1:Integer(λ),weights(ranks),N,replace = false)
  end
  return rank
end

# Tournament selection
function tournament(groupSize::Integer)
  groupSize <= 0 && error("Group size needs to be positive")
  function tournamentN(fitness::Vector{Float64},N::Integer)
    selection = Array(Int, N)
    nF = length(fitness)

    for i in 1:N
      contenders = sample(1:nF,groupSize,replace = true)
      winner = first(contenders)
      wf = fitness[winner]
      for k in 2:groupSize
        contender = contenders[k]
        if fitness[contender] > wf
          winner = contender
          wf = fitness[contender]
        end
      end
      selection[i] = winner
    end
    return selection
  end
  return tournamentN
end

# Roulette wheel (proportionate selection) selection
function roulette(fitness::Vector{Float64}, N::Integer)
  n = length(fitness)
  wv = weights(fitness)
  sample(1:n,wv,N,replace = true)
end

# Stochastic universal sampling (SUS)
function sus(fitness::Vector{Float64}, N::Integer)
  n = length(fitness)
  @assert N <= n "Number of selections (N) must be less than numeber of elements"
  wv = weights(fitness)
  sample(1:n,wv,N,replace = false)
end


text:
# ---
# title: 1095. Find in Mountain Array
# id: problem1095
# author: Indigo
# date: 2022-03-08
# difficulty: Hard
# categories: Binary Search
# link: <https://leetcode.com/problems/find-in-mountain-array/description/>
# hidden: true
# ---
#
# _(This problem is an   **interactive problem**.)_
#
# You may recall that an array `A` is a _mountain array_ if and only if:
#
#   * `A.length >= 3`
#   * There exists some `i` with `0 < i < A.length - 1` such that:
#     * `A[0] < A[1] < ... A[i-1] < A[i]`
#     * `A[i] > A[i+1] > ... > A[A.length - 1]`
#
# Given a mountain array `mountainArr`, return the **minimum**  `index` such
# that `mountainArr.get(index) == target`.  If such an `index` doesn't exist,
# return `-1`.
#
# **You can 't access the mountain array directly.**  You may only access the
# array using a `MountainArray` interface:
#
#   * `MountainArray.get(k)` returns the element of the array at index `k` (0-indexed).
#   * `MountainArray.length()` returns the length of the array.
#
# Submissions making more than `100` calls to `MountainArray.get` will be judged
# _Wrong Answer_.  Also, any solutions that attempt to circumvent the judge will
# result in disqualification.
#
#
#
# **Example 1:**
#
#
#
#     Input: array = [1,2,3,4,5,3,1], target = 3
#     Output: 2
#     Explanation: 3 exists in the array, at index=2 and index=5. Return the minimum index, which is 2.
#
# **Example 2:**
#
#
#
#     Input: array = [0,1,2,4,2,1], target = 3
#     Output: -1
#     Explanation: 3 does not exist in the array, so we return -1.
#
#
#
#
# **Constraints:**
#
#   * `3 <= mountain_arr.length() <= 10000`
#   * `0 <= target <= 10^9`
#   * `0 <= mountain_arr.get(index) <= 10^9`
#
#
## @lc code=start
using LeetCode

struct MountainArray
    vec::Vector{Int}
end

Base.length(ma::MountainArray) = length(ma.vec)
Base.getindex(ma::MountainArray, i) = ma.vec[i]

function find_in_mountain_array(ma::MountainArray, target::Int)
    len = length(ma)
    i, j = 1, len
    while i != j
        mid = (i + j) >> 1
        ma[mid] > ma[mid + 1] ? (j = mid) : (i = mid + 1)
    end
    idx1 = searchsortedfirst(1:i, i; by=i -> ma[i] >= target)
    ma[idx1] == target && return idx1
    idx2 = searchsortedfirst((i + 1):len, i; by=i -> ma[i] <= target) + i
    ma[idx2] == target && return idx2
    return -1
end
## @lc code=end


text:
## Curve Manipulations
struct RateCombination{T,U,V} <: AbstractYieldCurve
    r1::T
    r2::U
    op::V
end
__ratetype(::Type{RateCombination{T,U,V}}) where {T,U,V}= __ratetype(T)

rate(rc::RateCombination, time) = rc.op(rate(rc.r1, time), rate(rc.r2, time))
function discount(rc::RateCombination, time)
    a1 = discount(rc.r1, time)^(-1 / time) - 1
    a2 = discount(rc.r2, time)^(-1 / time) - 1
    return 1 / (1 + rc.op(a1, a2))^time
end

Base.zero(rc::RateCombination, time) = zero(rc,time,Periodic(1))
function Base.zero(rc::RateCombination, time, cf::C) where {C<:CompoundingFrequency}
    d = discount(rc,time)
    i = Periodic(1/d^(1/time)-1,1)
    return convert(cf, i) # c.zero is a curve of continuous rates represented as floats. explicitly wrap in continuous before converting
end

"""
    Yields.AbstractYieldCurve + Yields.AbstractYieldCurve

The addition of two yields will create a `RateCombination`. For `rate`, `discount`, and `accumulation` purposes the spot rates of the two curves will be added together.
"""
function Base.:+(a::AbstractYieldCurve, b::AbstractYieldCurve)
    return RateCombination(a, b, +)
end

function Base.:+(a::Constant, b::Constant)
    a_kind = rate(a).compounding
    rate_new_basis = rate(convert(a_kind, rate(b)))
    return Constant(
        Rate(
            rate(a.rate) + rate_new_basis,
            a_kind
        )
    )
end

function Base.:+(a::T, b) where {T<:AbstractYieldCurve}
    return a + Constant(b)
end

function Base.:+(a, b::T) where {T<:AbstractYieldCurve}
    return Constant(a) + b
end

"""
    Yields.AbstractYieldCurve - Yields.AbstractYieldCurve

The subtraction of two yields will create a `RateCombination`. For `rate`, `discount`, and `accumulation` purposes the spot rates of the second curves will be subtracted from the first.
"""
function Base.:-(a::AbstractYieldCurve, b::AbstractYieldCurve)
    return RateCombination(a, b, -)
end

function Base.:-(a::Constant, b::Constant)
    a_kind = rate(a).compounding
    rate_new_basis = rate(convert(a_kind, rate(b)))
    return Constant(
        Rate(
            rate(a.rate) - rate_new_basis,
            a_kind
        )
    )
end

function Base.:-(a::T, b) where {T<:AbstractYieldCurve}
    return a - Constant(b)
end

function Base.:-(a, b::T) where {T<:AbstractYieldCurve}
    return Constant(a) - b
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module Types

export VersionInterval, VersionSet, Requires, Available, Fixed, merge_requires!, satisfies
import Base: show, isempty, in, intersect, ==, hash, copy, deepcopy_internal

immutable VersionInterval
    lower::VersionNumber
    upper::VersionNumber
end
VersionInterval(lower::VersionNumber) = VersionInterval(lower,typemax(VersionNumber))
VersionInterval() = VersionInterval(typemin(VersionNumber))

show(io::IO, i::VersionInterval) = print(io, "[$(i.lower),$(i.upper))")
isempty(i::VersionInterval) = i.upper <= i.lower
in(v::VersionNumber, i::VersionInterval) = i.lower <= v < i.upper
intersect(a::VersionInterval, b::VersionInterval) = VersionInterval(max(a.lower,b.lower), min(a.upper,b.upper))
==(a::VersionInterval, b::VersionInterval) = a.lower == b.lower && a.upper == b.upper
hash(i::VersionInterval, h::UInt) = hash((i.lower, i.upper), h + (0x0f870a92db508386 % UInt))

immutable VersionSet
    intervals::Vector{VersionInterval}
end
function VersionSet(versions::Vector{VersionNumber})
    intervals = VersionInterval[]
    if isempty(versions)
        push!(intervals, VersionInterval())
    else
        isodd(length(versions)) && push!(versions, typemax(VersionNumber))
        while !isempty(versions)
            push!(intervals, VersionInterval(shift!(versions), shift!(versions)))
        end
    end
    VersionSet(intervals)
end
VersionSet(versions::VersionNumber...) = VersionSet(VersionNumber[versions...])

show(io::IO, s::VersionSet) = join(io, s.intervals, " ∪ ")
isempty(s::VersionSet) = all(isempty, s.intervals)
in(v::VersionNumber, s::VersionSet) = any(i->in(v,i), s.intervals)
function intersect(A::VersionSet, B::VersionSet)
    ivals = vec([ intersect(a,b) for a in A.intervals, b in B.intervals ])
    ivals = copy(ivals) # temporary bandaid for issue #4592
    filter!(i->!isempty(i), ivals)
    sort!(ivals, by=i->i.lower)
    VersionSet(ivals)
end
==(A::VersionSet, B::VersionSet) = A.intervals == B.intervals
hash(s::VersionSet, h::UInt) = hash(s.intervals, h + (0x2fd2ca6efa023f44 % UInt))
deepcopy_internal(vs::VersionSet, ::ObjectIdDict) = VersionSet(copy(vs.intervals))

typealias Requires Dict{String,VersionSet}

function merge_requires!(A::Requires, B::Requires)
    for (pkg,vers) in B
        A[pkg] = haskey(A,pkg) ? intersect(A[pkg],vers) : vers
    end
    return A
end

satisfies(pkg::AbstractString, ver::VersionNumber, reqs::Requires) =
    !haskey(reqs, pkg) || in(ver, reqs[pkg])

immutable Available
    sha1::String
    requires::Requires
end

==(a::Available, b::Available) = a.sha1 == b.sha1 && a.requires == b.requires
hash(a::Available, h::UInt) = hash((a.sha1, a.requires), h + (0xbc8ae0de9d11d972 % UInt))
copy(a::Available) = Available(a.sha1, copy(a.requires))

show(io::IO, a::Available) = isempty(a.requires) ?
    print(io, "Available(", repr(a.sha1), ")") :
    print(io, "Available(", repr(a.sha1), ",", a.requires, ")")

immutable Fixed
    version::VersionNumber
    requires::Requires
end
Fixed(v::VersionNumber) = Fixed(v,Requires())

==(a::Fixed, b::Fixed) = a.version == b.version && a.requires == b.requires
hash(f::Fixed, h::UInt) = hash((f.version, f.requires), h + (0x68628b809fd417ca % UInt))

show(io::IO, f::Fixed) = isempty(f.requires) ?
    print(io, "Fixed(", repr(f.version), ")") :
    print(io, "Fixed(", repr(f.version), ",", f.requires, ")")

# TODO: Available & Fixed are almost the same – merge them?
# Free could include the same information too, it just isn't
# required by anything that processes these things.

end # module


text:
"""
    SubDataFrame{<:AbstractDataFrame, <:AbstractIndex, <:AbstractVector{Int}} <: AbstractDataFrame

A view of an `AbstractDataFrame`. It is returned by a call to the `view` function
on an `AbstractDataFrame` if a collections of rows and columns are specified.

A `SubDataFrame` is an `AbstractDataFrame`, so expect that most
DataFrame functions should work. Such methods include `describe`,
`summary`, `nrow`, `size`, `by`, `stack`, and `join`.

If the selection of columns in a parent data frame is passed as `:` (a colon)
then `SubDataFrame` will always have all columns from the parent,
even if they are added or removed after its creation.

# Examples
```jldoctest
julia> df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
                      b = repeat([2, 1], outer=[4]),
                      c = 1:8)
8×3 DataFrame
 Row │ a      b      c
     │ Int64  Int64  Int64
─────┼─────────────────────
   1 │     1      2      1
   2 │     2      1      2
   3 │     3      2      3
   4 │     4      1      4
   5 │     1      2      5
   6 │     2      1      6
   7 │     3      2      7
   8 │     4      1      8

julia> sdf1 = view(df, :, 2:3) # column subsetting
8×2 SubDataFrame
 Row │ b      c
     │ Int64  Int64
─────┼──────────────
   1 │     2      1
   2 │     1      2
   3 │     2      3
   4 │     1      4
   5 │     2      5
   6 │     1      6
   7 │     2      7
   8 │     1      8

julia> sdf2 = @view df[end:-1:1, [1, 3]]  # row and column subsetting
8×2 SubDataFrame
 Row │ a      c
     │ Int64  Int64
─────┼──────────────
   1 │     4      8
   2 │     3      7
   3 │     2      6
   4 │     1      5
   5 │     4      4
   6 │     3      3
   7 │     2      2
   8 │     1      1

julia> sdf3 = groupby(df, :a)[1]  # indexing a GroupedDataFrame returns a SubDataFrame
2×3 SubDataFrame
 Row │ a      b      c
     │ Int64  Int64  Int64
─────┼─────────────────────
   1 │     1      2      1
   2 │     1      2      5
```
"""
struct SubDataFrame{D<:AbstractDataFrame, S<:AbstractIndex, T<:AbstractVector{Int}} <: AbstractDataFrame
    parent::D
    colindex::S
    rows::T # maps from subdf row indexes to parent row indexes
end

Base.@propagate_inbounds function SubDataFrame(parent::DataFrame, rows::AbstractVector{Int}, cols)
    @boundscheck if !checkindex(Bool, axes(parent, 1), rows)
        throw(BoundsError(parent, (rows, cols)))
    end
    SubDataFrame(parent, SubIndex(index(parent), cols), rows)
end
Base.@propagate_inbounds SubDataFrame(parent::DataFrame, ::Colon, cols) =
    SubDataFrame(parent, axes(parent, 1), cols)
@inline SubDataFrame(parent::DataFrame, row::Integer, cols) =
    throw(ArgumentError("invalid row index: $row of type $(typeof(row))"))

Base.@propagate_inbounds function SubDataFrame(parent::DataFrame, rows::AbstractVector{<:Integer}, cols)
    if any(x -> x isa Bool, rows)
        throw(ArgumentError("invalid row index of type `Bool`"))
    end
    return SubDataFrame(parent, convert(Vector{Int}, rows), cols)
end

Base.@propagate_inbounds function SubDataFrame(parent::DataFrame, rows::AbstractVector{Bool}, cols)
    if length(rows) != nrow(parent)
        throw(ArgumentError("invalid length of `AbstractVector{Bool}` row index " *
                            "(got $(length(rows)), expected $(nrow(parent)))"))
    end
    return SubDataFrame(parent, _findall(rows), cols)
end

Base.@propagate_inbounds function SubDataFrame(parent::DataFrame, rows::AbstractVector, cols)
    if !all(x -> (x isa Integer) && !(x isa Bool), rows)
        throw(ArgumentError("only `Integer` indices are accepted in `rows`"))
    end
    return SubDataFrame(parent, convert(Vector{Int}, rows), cols)
end

Base.@propagate_inbounds SubDataFrame(sdf::SubDataFrame, rowind, cols) =
    SubDataFrame(parent(sdf), rows(sdf)[rowind], parentcols(index(sdf), cols))
Base.@propagate_inbounds SubDataFrame(sdf::SubDataFrame, rowind::Bool, cols) =
    throw(ArgumentError("invalid row index of type Bool"))
Base.@propagate_inbounds SubDataFrame(sdf::SubDataFrame, rowind, ::Colon) =
    if index(sdf) isa Index # sdf was created using : as row selector
        SubDataFrame(parent(sdf), rows(sdf)[rowind], :)
    else
        SubDataFrame(parent(sdf), rows(sdf)[rowind], parentcols(index(sdf), :))
    end
Base.@propagate_inbounds SubDataFrame(sdf::SubDataFrame, rowind::Bool, ::Colon) =
    throw(ArgumentError("invalid row index of type Bool"))
Base.@propagate_inbounds SubDataFrame(sdf::SubDataFrame, ::Colon, cols) =
    SubDataFrame(parent(sdf), rows(sdf), parentcols(index(sdf), cols))
@inline SubDataFrame(sdf::SubDataFrame, ::Colon, ::Colon) = sdf

rows(sdf::SubDataFrame) = getfield(sdf, :rows)
Base.parent(sdf::SubDataFrame) = getfield(sdf, :parent)
Base.parentindices(sdf::SubDataFrame) = (rows(sdf), parentcols(index(sdf)))

Base.@propagate_inbounds Base.view(adf::AbstractDataFrame, rowinds, colind::ColumnIndex) =
    view(adf[!, colind], rowinds)
Base.@propagate_inbounds Base.view(adf::AbstractDataFrame, ::typeof(!), colind::ColumnIndex) =
    view(adf[!, colind], :)
@inline Base.view(adf::AbstractDataFrame, rowinds, colind::Bool) =
    throw(ArgumentError("invalid column index $colind of type `Bool`"))
Base.@propagate_inbounds Base.view(adf::AbstractDataFrame, rowinds,
                                   colinds::MultiColumnIndex) =
    SubDataFrame(adf, rowinds, colinds)
Base.@propagate_inbounds Base.view(adf::AbstractDataFrame, rowinds::typeof(!),
                                   colinds::MultiColumnIndex) =
    SubDataFrame(adf, :, colinds)
Base.@propagate_inbounds Base.view(adf::AbstractDataFrame, rowinds::Not,
                                   colinds::MultiColumnIndex) =
    SubDataFrame(adf, axes(adf, 1)[rowinds], colinds)

##############################################################################
##
## AbstractDataFrame interface
##
##############################################################################

index(sdf::SubDataFrame) = getfield(sdf, :colindex)

nrow(sdf::SubDataFrame) = ncol(sdf) > 0 ? length(rows(sdf))::Int : 0
ncol(sdf::SubDataFrame) = length(index(sdf))

Base.@propagate_inbounds Base.getindex(sdf::SubDataFrame, rowind::Integer, colind::ColumnIndex) =
    parent(sdf)[rows(sdf)[rowind], parentcols(index(sdf), colind)]
Base.@propagate_inbounds Base.getindex(sdf::SubDataFrame, rowind::Bool, colind::ColumnIndex) =
    throw(ArgumentError("invalid row index of type Bool"))

Base.@propagate_inbounds Base.getindex(sdf::SubDataFrame, rowinds::Union{AbstractVector, Not},
                                       colind::ColumnIndex) =
    parent(sdf)[rows(sdf)[rowinds], parentcols(index(sdf), colind)]
Base.@propagate_inbounds Base.getindex(sdf::SubDataFrame, ::Colon, colind::ColumnIndex) =
    parent(sdf)[rows(sdf), parentcols(index(sdf), colind)]
Base.@propagate_inbounds Base.getindex(sdf::SubDataFrame, ::typeof(!), colind::ColumnIndex) =
    view(parent(sdf), rows(sdf), parentcols(index(sdf), colind))

Base.@propagate_inbounds Base.getindex(sdf::SubDataFrame, rowinds::Union{AbstractVector, Not},
                                       colinds::MultiColumnIndex) =
    parent(sdf)[rows(sdf)[rowinds], parentcols(index(sdf), colinds)]
Base.@propagate_inbounds Base.getindex(sdf::SubDataFrame, ::Colon,
                                       colinds::MultiColumnIndex) =
    parent(sdf)[rows(sdf), parentcols(index(sdf), colinds)]
Base.@propagate_inbounds Base.getindex(df::SubDataFrame, row_ind::typeof(!),
                                       col_inds::MultiColumnIndex) =
    select(df, col_inds, copycols=false)

Base.@propagate_inbounds function Base.setindex!(sdf::SubDataFrame, val::Any, idx::CartesianIndex{2})
    return setindex!(sdf, val, idx[1], idx[2])
end

Base.@propagate_inbounds function Base.setindex!(sdf::SubDataFrame, val::Any, ::Colon, colinds::Any)
    if colinds isa SymbolOrString && columnindex(sdf, colinds) == 0
        if !is_column_insertion_allowed(sdf)
            throw(ArgumentError("creating new columns in a SubDataFrame that subsets " *
                                "columns of its parent data frame is disallowed"))
        end
        if !(val isa AbstractVector && nrow(sdf) == length(val))
            throw(ArgumentError("Assigned value must be a vector with length " *
                                "equal to number of rows in the SubDataFrame"))
        end
        T = eltype(val)
        newcol = similar(val, Union{T, Missing}, nrow(parent(sdf)))
        fill!(newcol, missing)
        newcol[rows(sdf)] = val
        parent(sdf)[!, colinds] = newcol
    else
        parent(sdf)[rows(sdf), parentcols(index(sdf), colinds)] = val
    end
    return sdf
end

# TODO: in the future, when refactoring source code
# (presumably when we would first define all the types that the package provides
#  and then define methods for them)
# consider merging SubDataFrame and DataFrame setindex! methods

function Base.setindex!(sdf::SubDataFrame, v::AbstractVector,
                        ::typeof(!), col_ind::ColumnIndex)
    if col_ind isa Union{Signed, Unsigned} && !(1 <= col_ind <= ncol(sdf))
        throw(ArgumentError("Cannot assign to non-existent column: $col_ind"))
    end
    if col_ind isa SymbolOrString && columnindex(sdf, col_ind) == 0
        if !is_column_insertion_allowed(sdf)
            throw(ArgumentError("creating new columns in a SubDataFrame that subsets " *
                                "columns of its parent data frame is disallowed"))
        end
        sdf[:, col_ind] = v
    else
        pdf = parent(sdf)
        p_col_ind = parentcols(index(sdf), col_ind)
        old_col = pdf[!, p_col_ind]
        T = eltype(old_col)
        S = eltype(v)
        newcol = similar(old_col, promote_type(T, S), length(old_col))
        newcol .= old_col
        newcol[rows(sdf)] = v
        pdf[!, p_col_ind] = newcol
    end
    return sdf
end

for T in MULTICOLUMNINDEX_TUPLE
    @eval function Base.setindex!(sdf::SubDataFrame,
                                  new_df::AbstractDataFrame,
                                  row_inds::typeof(!),
                                  col_inds::$T)
        idxs = index(sdf)[col_inds]
        if view(_names(sdf), idxs) != _names(new_df)
            throw(ArgumentError("Column names in source and target data frames do not match"))
        end
        for (j, col) in enumerate(idxs)
            sdf[!, col] = new_df[!, j]
        end
        return sdf
    end

    @eval function Base.setindex!(sdf::SubDataFrame,
                                  mx::AbstractMatrix,
                                  row_inds::typeof(!),
                                  col_inds::$T)
        idxs = index(sdf)[col_inds]
        if size(mx, 2) != length(idxs)
            throw(DimensionMismatch("number of selected columns ($(length(idxs))) " *
                                    "and number of columns in " *
                                    "matrix ($(size(mx, 2))) do not match"))
        end
        for (j, col) in enumerate(idxs)
            sdf[!, col] = view(mx, :, j)
        end
        return sdf
    end
end

Base.@propagate_inbounds function Base.setindex!(sdf::SubDataFrame, val::Any, rowinds::Any, colinds::Any)
    parent(sdf)[rows(sdf)[rowinds], parentcols(index(sdf), colinds)] = val
    return sdf
end

Base.@propagate_inbounds Base.setindex!(sdf::SubDataFrame, val::Any, rowinds::Bool, colinds::Any) =
    throw(ArgumentError("invalid row index of type Bool"))

Base.setproperty!(df::SubDataFrame, col_ind::Symbol, v::AbstractVector) =
    (df[!, col_ind] = v)
Base.setproperty!(df::SubDataFrame, col_ind::AbstractString, v::AbstractVector) =
    (df[!, col_ind] = v)
Base.setproperty!(::SubDataFrame, col_ind::Symbol, v::Any) =
    throw(ArgumentError("It is only allowed to pass a vector as a column of a SubDataFrame. " *
                        "Instead use `df[!, col_ind] .= v` if you want to use broadcasting."))
Base.setproperty!(::SubDataFrame, col_ind::AbstractString, v::Any) =
    throw(ArgumentError("It is only allowed to pass a vector as a column of a SubDataFrame. " *
                        "Instead use `df[!, col_ind] .= v` if you want to use broadcasting."))

##############################################################################
##
## Miscellaneous
##
##############################################################################

Base.copy(sdf::SubDataFrame) = parent(sdf)[rows(sdf), parentcols(index(sdf), :)]

Base.delete!(df::SubDataFrame, ind) =
    throw(ArgumentError("SubDataFrame does not support deleting rows"))

function DataFrame(sdf::SubDataFrame; copycols::Bool=true)
    if copycols
        sdf[:, :]
    else
        DataFrame(collect(eachcol(sdf)), _names(sdf), copycols=false)
    end
end

Base.convert(::Type{DataFrame}, sdf::SubDataFrame) = DataFrame(sdf)

# this function tests if it is allowed to add columns to passed SubDataFrame
# currently it is only allowed when SubDataFrame was created with : as column selector
# which results in using Index as its index (as opposed to other columns selectors
# which result in SubIndex)
function is_column_insertion_allowed(df::AbstractDataFrame)
    if df isa DataFrame
        return true
    elseif df isa SubDataFrame
        return getfield(df, :colindex) isa Index
    end
    throw(ArgumentError("Unsupported data frame type"))
end

function _replace_columns!(sdf::SubDataFrame, newdf::DataFrame)
    colsmatch = _names(sdf) == _names(newdf)

    if !(colsmatch || is_column_insertion_allowed(sdf))
        throw(ArgumentError("changing the sequence of column names in a SubDataFrame " *
                            "that subsets columns of its parent data frame is disallowed"))
    end

    for colname in _names(newdf)
        # This will allocate a fresh column in parent(sdf) for each colname
        sdf[!, colname] = newdf[!, colname]
    end

    # If columns did not match this means that we have either:
    # 1. inserted some columns into pdf
    # or
    # 2. requested to reorder the existing columns
    # and that operation was allowed.
    # Therefore we need to update the parent of sdf in place to make sure
    # it holds only the required target columns in a correct order.
    if !colsmatch
        pdf = parent(sdf)
        @assert pdf isa DataFrame
        select!(pdf, _names(newdf))
    end

    return sdf
end