text:
function range2ind(ts::AbstractRange, t)
  indx = round(Int, (t-first(ts))/step(ts))
  indx += one(indx)
  indx = minimum((maximum((indx,one(indx))),length(ts)))
  return indx
end

function range2ind(ts::AbstractVector, t)
  r = searchsorted(ts, t)
  r1 = minimum((first(r),length(ts)))
  r2 = maximum((last(r),one(last(r))))
  return abs(ts[r1] - t) < abs(ts[r2] - t) ? r1 : r2
end

unpackx(a) = @view a[1:end-1]
mypack(a::SArray,c::Number) = SVector(a..., c)
mypack(a,c::Number) = [a; c]

# guided drift
function (G::GuidingDriftCache)(du,u,p,t)
  @unpack k, message = G
  @unpack f, g, constant_diffusity = k
  @unpack ktilde, ts, soldis, sol, filter = message

  x = unpackx(u)
  dx = unpackx(du)
  d = length(x)

  # find cursor
  @inbounds cur_time = range2ind(ts, t)

  if isapprox(t, ts[cur_time]; atol = 1000eps(typeof(t)), rtol = 1000eps(t))
    # non-interpolating version
    # take care for multivariate case here if P isa Matrix, ν  isa Vector, c isa Scalar
    # ν, P, c or F, H, c parametrization = μ,Σ,c
    μ = soldis[cur_time][1]
    Σ = soldis[cur_time][2]
  else
    μ = @view sol(t)[1:d]
    Σ = reshape(@view(sol(t)[d+1:d+d*d]), d, d)
  end

  if !(filter isa InformationFilter)
    r = Σ\(μ - x)
    du[end] = dot(f(x,p,t) - ktilde.f(x,ktilde.p,t), r)
    if !constant_diffusity
      du[end] -= 0.5*tr((outer_(g(x,p,t)) - outer_(ktilde.g(x,ktilde.p,t)))*(inv(Σ) - outer_(r)))
    end
    dx[:] .= vec(f(x,p,t) + (outer_(g(x,p,t))*r)) # evolution guided by observations
  else
    du[end] = ..
    dx[:] .= ..
  end

  return nothing
end

function (G::GuidingDriftCache)(u,p,t)
  @unpack k, message = G
  @unpack f, g, constant_diffusity = k
  @unpack ktilde, ts, soldis, sol = message

  x = unpackx(u)
  d = length(x)

  # find cursor
  @inbounds cur_time = range2ind(ts, t)

  if isapprox(t, ts[cur_time]; atol = 1000eps(typeof(t)), rtol = 1000eps(t))
    # non-interpolating version
    # take care for multivariate case here if P isa Matrix, ν  isa Vector, c isa Scalar
    # ν, P, c or F, H, c parametrization = μ,Σ,c
    μ = soldis[cur_time][1]
    Σ = soldis[cur_time][2]
  else
    μ = @view sol(t)[1:d]
    Σ = reshape(@view(sol(t)[d+1:d+d*d]), d, d)
  end

  if !(filter isa InformationFilter)
    r = Σ\(μ .- x)

    dl = dot(f(x,p,t) -  ktilde.f(x,ktilde.p,t), r)
    if !constant_diffusity
      dl -= 0.5*tr((outer_(g(x,p,t)) - outer_(ktilde.g(x,ktilde.p,t)))*(inv(Σ) - outer_(r)))
    end
    dx = vec(f(x,p,t) + outer_(g(x,p,t))*r) # evolution guided by observations
  else
    dl = ..
    dx = ..
  end

  return mypack(dx, dl)
end

# guided diffusion
function (G::GuidingDiffusionVectorCache)(du,u,p,t)
  @unpack g = G

  x = @view u[1:end-1]
  du[1:end-1] .= g(x,p,t)
  return nothing
end

function (G::GuidingDiffusionVectorCache)(u,p,t)
  @unpack g = G

  x = @view u[1:end-1]
  dx = g(x,p,t)
  return mypack(dx, zero(eltype(u)))
end

function (G::GuidingDiffusionCache)(du,u,p,t)
  @unpack g = G

  x = @view u[1:end-1]
  du[1:end-1,1:end-1] .= g(x,p,t)
  return nothing
end

function (G::GuidingDiffusionCache)(u,p,t)
  @unpack g, padded_size = G

  x = @view u[1:end-1]
  dx = g(x,p,t)
  return PaddedView(zero(eltype(dx)), dx, padded_size)
end

function construct_forwardguiding_Problem(k::Union{SDEKernel,SDEKernel!}, message, u0, Z, inplace, alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm})
  @unpack f, g, trange, p, noise_rate_prototype = k

  guided_f = GuidingDriftCache(k,message)

  if noise_rate_prototype===nothing
    # diagonal noise in StochasticDiffEq
    guided_g = GuidingDiffusionVectorCache(g)
    _noise_rate_prototype = nothing
  else
    # non-diagonal noise in StochasticDiffEq
    padded_size = size(noise_rate_prototype) .+ 1
    _noise_rate_prototype = PaddedView(zero(eltype(noise_rate_prototype)), noise_rate_prototype, padded_size)
    guided_g = GuidingDiffusionCache(g,padded_size)
  end

  if Z===nothing
    noise = Z
  else
    if Z.curW isa AbstractVector
      if Z isa DiffEqNoiseProcess.NoiseGrid
        noise = AugmentedNoiseGrid(Z.t,Z.W,Z.Z)
      else
        error("Please pass a NoiseGrid to forwardguiding.")
      end
    else
      # scalar case
      noise = Z
    end
  end

  prob = SDEProblem{inplace}(guided_f, guided_g, u0, get_tspan(trange), p, noise=noise,
                             noise_rate_prototype=_noise_rate_prototype)

  return prob
end

function construct_forwardguiding_Problem(k::Union{SDEKernel,SDEKernel!}, message, (u0,ll0), Z, inplace, alg::AbstractInternalSolver; P=nothing)

  if P===nothing
    if inplace
      k isa SDEKernel && error("SDEKernel cannot be used with the inplace=true option. SDEKernel! has to be used instead.")
      P = GuidedSDE!(k, message)
    else
      P = GuidedSDE(k, message)
    end
  end

  # du [index, time, space, ll]
  u, dz, du, Z = construct_sample_Problem(k, u0, Z, alg, ll0)

  return u, dz, du, Z, P
end

function forwardguiding(k::SDEKernel, message, (x0, ll0), alg=EM(false), Z=nothing;
  dt=get_dt(k.trange),
  inplace=true, kwargs...)

  # check that message.ts is sorted
  !issorted(message.ts) && error("Something went wrong. Message.ts is not sorted! Please report this.")

  sol, ll = _forwardguiding(k::SDEKernel, message, (x0, ll0), alg, Z;
    dt=dt, inplace=inplace, kwargs...)
end

function _forwardguiding(k::SDEKernel, message, (x0, ll0), alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}, Z;
    dt=dt, isadaptive=StochasticDiffEq.isadaptive(alg),
    numtraj=nothing, ensemblealg=EnsembleThreads(), output_func=(sol,i) -> (sol[end],false),
    inplace=inplace, kwargs...)

  @unpack f, g, trange, p = k

  u0 = mypack(x0,ll0)

  # construct guiding SDE problem
  prob = construct_forwardguiding_Problem(k, message, u0, Z, inplace, alg)

  if numtraj==nothing
    if !isadaptive
      sol = solve(prob, alg, tstops=message.ts, adaptive=isadaptive; kwargs...)
    else
      sol = solve(prob, alg, dt=dt, adaptive=isadaptive; kwargs...)
    end
    if Z===nothing
      W = [x[1:length(x0)] for x in sol.W.W]
    else
      W = Z.W
    end
    return (sol[end][1:length(x0)],sol[end][end]), (sol.t, [x[1:length(x0)] for x in sol.u], W)
  else
    ensembleprob = EnsembleProblem(prob, output_func = output_func)
    if !isadaptive
      sol = solve(ensembleprob, alg, ensemblealg=ensemblealg,
          tstops=message.ts, trajectories=numtraj, adaptive=isadaptive; kwargs...)
    else
      sol = solve(ensembleprob, alg, ensemblealg=ensemblealg,
          dt=dt, adaptive=isadaptive, trajectories=numtraj; kwargs...)
    end
    return ([x[1:length(x0)] for x in sol.u], last.(sol.u)), nothing
  end
end

function _forwardguiding(k::SDEKernel, message, (x0, ll0), alg::AbstractInternalSolver, Z; P=nothing, save=true, inplace=inplace, kwargs...)

  @unpack f, g, trange, p = k

  u0 = mypack(x0,ll0)

  u, dz, du, Z, P = construct_forwardguiding_Problem(k, message, (x0, ll0), Z, inplace, alg; P=P)

  if save
    uu = typeof(u)[]
  else
    uu = nothing
  end
  uu, uT = solve!(alg, uu, u, Z, P)
  if save
    rest = (getindex.(uu,2), getindex.(uu,3), getindex.(Z,3))
  else
    rest = nothing
  end
  return (uT[3],uT[end]), rest
end

function tangent!(du, u, dz, P::GuidedSDE!)
  @unpack kernel, message = P
  @unpack ktilde, ts, soldis = message
  @unpack f, gstep!, ws, p, noise_rate_prototype = kernel
  x, t = u[3], u[2]
  cur_time = u[1]

  μ = soldis[cur_time][1]
  Σ = soldis[cur_time][2]
  r = Σ\(μ .- x)

  f(du[3], u[3], p, u[2])
  du[3] += outer_(g(x,p,t))*r
  du[3] .*= dz[2]
  gstep!(du[3], ws, u[3], p, u[2], dz[3], noise_rate_prototype)

  dl = dot(f(x,p,t) - ktilde.f(x,ktilde.p,t), r)
  if !constant_diffusity
    dl -= 0.5*tr((outer_(g(x,p,t)) - outer_(ktilde.g(x,ktilde.p,t)))*(inv(Σ) - outer_(r)))
  end
  (dz[1], dz[2], du[3], dl)
end

function dZ!(u, dz, Z, P::GuidedSDE)
  i = u[1]
  (Z[i+1][1] - Z[i][1], Z[i+1][2] - Z[i][2],  Z[i+1][3] -  Z[i][3])
end

function exponential_map!(u, du, P::GuidedSDE)
  (u[1] + du[1], u[2] + du[2], u[3] + du[3], u[4] + du[4])
end

function tangent!(du, u, dz, P::GuidedSDE)
  @unpack kernel, message = P
  @unpack ktilde, ts, soldis = message
  @unpack f, g, p, noise_rate_prototype, constant_diffusity = kernel
  x, t = u[3], u[2]
  cur_time = u[1]

  μ = soldis[cur_time][1]
  Σ = soldis[cur_time][2]

  r = Σ\(μ - x)

  du3 = (f(x,p,t) + outer_(g(x,p,t))*r)*dz[2]
  if noise_rate_prototype===nothing
    du3 += g(x,p,t).*dz[3]
  else
    du3 += g(x,p,t)*dz[3]
  end

  dl = dot(f(x,p,t) - ktilde.f(x,ktilde.p,t), r)
  if !constant_diffusity
    dl -= 0.5*tr((outer_(g(x,p,t)) - outer_(ktilde.g(x,ktilde.p,t)))*(inv(Σ) - outer_(r)))
  end
  (dz[1], dz[2], du3, dl)
end


text:
using Pkg

reqs = ["Conda", "PyCall"]

# By kristoffer.carlsson
# https://discourse.julialang.org/t/how-to-use-pkg-dependencies-instead-of-pkg-installed/36416/10
isinstalled(pkg::String) = any(x -> x.name == pkg && x.is_direct_dep, values(Pkg.dependencies()))

function pkgSetup()
	toInstall = filter!(isinstalled, reqs)
	println(toInstall)
	map(x -> Pkg.add(x), toInstall)
end

using Conda
function condaSetup()
	Conda.add("magic-wormhole",channel="conda-forge")
end

using PyCall
function pySetup()
	os = pyimport("os")
end

function doSetup()
	pkgSetup()
	condaSetup()
	return pySetup()
end

osPyModule = doSetup()

function test()
	osPyModule.system("conda run wormhole --help")
end

function runCommand(str)
	osPyModule.system("conda activate base & $str")
end

function receive(word)
	runCommand("wormhole receive $word")
end

function sendText(txt)
	osPyModule.system("conda run wormhole send --text $txt")
	#osPyModule.system("conda activate base & wormhole send --text $txt")
end

#Todo macro for ease of use (no quotes)

ext: jl


text:
using HttpParser
using HttpCommon
using Base.Test
using Compat
import Compat: String

FIREFOX_REQ = tuple("GET /favicon.ico HTTP/1.1\r\n",
         "Host: 0.0.0.0=5000\r\n",
         "User-Agent: Mozilla/5.0 (X11; U; Linux i686; en-US; rv:1.9) Gecko/2008061015 Firefox/3.0\r\n",
         "Accept: text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8\r\n",
         "Accept-Language: en-us,en;q=0.5\r\n",
         "Accept-Encoding: gzip,deflate\r\n",
         "Accept-Charset: ISO-8859-1,utf-8;q=0.7,*;q=0.7\r\n",
         "Keep-Alive: 300\r\n",
         "Connection: keep-alive\r\n",
         "\r\n")

MALFORMED = tuple("GET /MALFORMED HTTP/1.1\r\n",
         "aaaaaaaaaaaaa:++++++++++\r\n",
         "\r\n")

TWO_CHUNKS_MULT_ZERO_END = tuple("POST /two_chunks_mult_zero_end HTTP/1.1\r\n",
         "Transfer-Encoding: chunked\r\n",
         "\r\n",
         "5\r\nhello\r\n",
         "6\r\n world\r\n",
         "000\r\n",
         "\r\n")

WEBSOCK = tuple("DELETE /chat HTTP/1.1\r\n",
        "Host: server.example.com\r\n",
        "Upgrade: websocket\r\n",
        "Connection: Upgrade\r\n",
        "Sec-WebSocket-Key: dGhlIHNhbXBsZSBub25jZQ==\r\n",
        "Origin: http://example.com\r\n",
        "Sec-WebSocket-Protocol: chat, superchat\r\n",
        "Sec-WebSocket-Version: 13\r\n",
        "\r\n",)

r = Request("", "", Dict{AbstractString,AbstractString}(), "")

function on_message_begin(parser)
    # Clear the resource when the message starts
    r.resource = ""
    return 0
end

function on_url(parser, at, len)
    # Concatenate the resource for each on_url callback
    r.resource = string(r.resource, unsafe_string(convert(Ptr{UInt8}, at), Int(len)))
    return 0
end

function on_status_complete(parser)
    return 0
end

function on_header_field(parser, at, len)
    header = unsafe_string(convert(Ptr{UInt8}, at), Int(len))
    # set the current header
    r.headers["current_header"] = header
    return 0
end

function on_header_value(parser, at, len)
    s = unsafe_string(convert(Ptr{UInt8}, at), Int(len))
    # once we know we have the header value, that will be the value for current header
    r.headers[r.headers["current_header"]] = s
    # reset current_header
    r.headers["current_header"] = ""
    return 0
end

function on_headers_complete(parser)
    p = unsafe_load(parser)
    # get first two bits of p.type_and_flags

    # The parser type are the bottom two bits
    # 0x03 = 00000011
    ptype = p.type_and_flags & 0x03
    # flags = p.type_and_flags >>> 3
    if ptype == 0
        r.method = http_method_str(convert(Int, p.method))
    end
    if ptype == 1
        r.headers["status_code"] = string(convert(Int, p.status_code))
    end
    r.headers["http_major"] = string(convert(Int, p.http_major))
    r.headers["http_minor"] = string(convert(Int, p.http_minor))
    r.headers["Keep-Alive"] = string(http_should_keep_alive(parser))
    return 0
end

function on_body(parser, at, len)
    append!(r.data, unsafe_wrap(Array, convert(Ptr{UInt8}, at), (len,)))
    return 0
end

function on_message_complete(parser)
    return 0
end

function on_chunk_header(parser)
    return 0
end

function on_chunk_complete(parser)
    return 0
end

c_message_begin_cb = cfunction(on_message_begin, HttpParser.HTTP_CB...)
c_url_cb = cfunction(on_url, HttpParser.HTTP_DATA_CB...)
c_status_complete_cb = cfunction(on_status_complete, HttpParser.HTTP_CB...)
c_header_field_cb = cfunction(on_header_field, HttpParser.HTTP_DATA_CB...)
c_header_value_cb = cfunction(on_header_value, HttpParser.HTTP_DATA_CB...)
c_headers_complete_cb = cfunction(on_headers_complete, HttpParser.HTTP_CB...)
c_body_cb = cfunction(on_body, HttpParser.HTTP_DATA_CB...)
c_message_complete_cb = cfunction(on_message_complete, HttpParser.HTTP_CB...)
c_body_cb = cfunction(on_body, HttpParser.HTTP_DATA_CB...)
c_message_complete_cb = cfunction(on_message_complete, HttpParser.HTTP_CB...)
c_chunk_header_cb = cfunction(on_chunk_header, HttpParser.HTTP_CB...)
c_chunk_complete_cb = cfunction(on_chunk_complete, HttpParser.HTTP_CB...)

function init(test::Tuple)
    # reset request
    # Moved this up for testing purposes
    r.method = ""
    r.resource = ""
    r.headers = Dict{AbstractString, AbstractString}()
    r.data = ""
    parser = Parser()
    http_parser_init(parser)
    settings = ParserSettings(c_message_begin_cb, c_url_cb,
                              c_status_complete_cb, c_header_field_cb,
                              c_header_value_cb, c_headers_complete_cb,
                              c_body_cb, c_message_complete_cb,
                              c_chunk_header_cb, c_chunk_complete_cb)

    for i=1:length(test)
        size = http_parser_execute(parser, settings, test[i])
    end

    # errno = parser.errno_and_upgrade & 0xf3
    # upgrade = parser.errno_and_upgrade >>> 7
end

init(FIREFOX_REQ)
@test r.method == "GET"
@test r.resource == "/favicon.ico"
@test r.headers["Host"] == "0.0.0.0=5000"
@test r.headers["User-Agent"] == "Mozilla/5.0 (X11; U; Linux i686; en-US; rv:1.9) Gecko/2008061015 Firefox/3.0"
@test r.headers["Accept"] == "text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8"
@test r.headers["Accept-Language"] == "en-us,en;q=0.5"
@test r.headers["Accept-Encoding"] == "gzip,deflate"
@test r.headers["Accept-Charset"] == "ISO-8859-1,utf-8;q=0.7,*;q=0.7"
@test r.headers["Keep-Alive"] == "1"
@test r.headers["Connection"] == "keep-alive"
@test isempty(r.data)
init(MALFORMED)
@test r.method == "GET"
@test r.resource == "/MALFORMED"
init(TWO_CHUNKS_MULT_ZERO_END)
@test r.method == "POST"
@test r.resource == "/two_chunks_mult_zero_end"
@test String(r.data) == "hello world"
# @test r.data == "hello\r\n5 world\r\n6"
init(WEBSOCK)
@test r.method == "DELETE"
@test r.resource == "/chat"

# URL Parser
const ParsedUrl = Dict{Symbol,AbstractString}

@test parse_url("//some_path") == (ParsedUrl(
    :UF_PATH=>"//some_path"
), 0)

@test parse_url("HTTP://www.example.com/") == (ParsedUrl(
    :UF_HOST=>"www.example.com",
    :UF_PATH=>"/",
    :UF_SCHEMA=>"HTTP"
), 0)

@test parse_url("HTTP://www.example.com") == (ParsedUrl(
    :UF_HOST=>"www.example.com",
    :UF_SCHEMA=>"HTTP"
), 0)

@test parse_url("http://user:aaa@www.example.com/") == (ParsedUrl(
    :UF_HOST=>"www.example.com",
    :UF_SCHEMA=>"http",
    :UF_USERINFO=>"user:aaa",
    :UF_PATH=>"/",
), 0)

@test parse_url("http://x.com/path?that%27s#all,%20folks") == (ParsedUrl(
    :UF_HOST=>"x.com",
    :UF_QUERY=>"that%27s",
    :UF_FRAGMENT=>"all,%20folks",
    :UF_SCHEMA=>"http",
    :UF_PATH=>"/path",
), 0)

@test parse_url("http://user:pass@example.com:8000/foo/bar?baz=quux#frag") == (ParsedUrl(
    :UF_PORT=>"8000",
    :UF_HOST=>"example.com",
    :UF_QUERY=>"baz=quux",
    :UF_FRAGMENT=>"frag",
    :UF_USERINFO=>"user:pass",
    :UF_PATH=>"/foo/bar",
    :UF_SCHEMA=>"http"
), 8000)

@test HttpParser.version() >= v"2.6"

p = Parser()
HttpParser.pause(p)
err = HttpParser.errno(p)
@test err == 0x1f
@test HttpParser.errno_name(err) == "HPE_PAUSED"
@test HttpParser.errno_description(err) == "parser is paused"
HttpParser.resume(p)
@test HttpParser.errno(p) == 0x00

ex = HttpParser.HttpParserError(0x1f)
@test ex.errno == 0x1f

info("All assertions passed!")


text:
function view_merafile(filename::String;
                        verbose::Bool=verbose_mode)

         return view_merafile(filename=filename,
                              verbose=verbose)
end

function view_merafile(;
                       filename::String="",
                       verbose::Bool=verbose_mode)

    viewoutput = Dict()

    if verbose
        println("Mera-file Contains:")
        println()
    end

    file = h5open(filename, "r")
    stored_datatypes = names(file)
    for datatype in stored_datatypes
        output = readheader(file, datatype, verbose)
        viewoutput[datatype] = output
    end

    s = filesize(filename)
    svalue, sunit = humanize(Float64(s), 3, "memory")
    viewoutput["FileSize"] = (svalue, sunit)

    if verbose
        println("Total file size: ", svalue, " ", sunit)
        println("-----------------------------------")
        println()
    end

    close(file)
    return viewoutput
end

function readheader(file, datatype::String, verbose::Bool)
    output = Dict()

    if known_datatype(datatype)
        group = file[datatype]
        path = group["info/path"]
        path = read(path)
        output["path"] = path

        merafile_version = group["merafile_version"]
        merafile_version = read(merafile_version)
        output["merafile_version"] = merafile_version

        simcode = group["info/simcode"]
        simcode = read(simcode)
        output["SimCode"] = simcode

        column_names = group["data/names"]
        column_names = Symbol.(read(column_names) )
        output["Variables"] = column_names

        total_fields = group["total_fields"]
        total_fields = read(total_fields)
        output["TotalFields"] = total_fields

        ctype = group["compression_type"]
        ctype = read(ctype)
        cn = group["compression"]
        cn = read(cn)
        output["CompressionType"] = ctype
        output["Compression"] = cn

        mem = group["memory"]
        mem = read(mem)
        mem1 = parse(Float64, mem[1])
        output["Memory"] = mem

        comment = group["comment"]
        comment = read(comment)
        output["Comment"] = comment

        if verbose
            println(datatype)
            println("-----------")
            println("Directory: ", path )
            println("merafile_version: ", merafile_version)
            println("Simulation code: ", simcode )
            println("Data variables: ", column_names )
            println("Total fields: ", total_fields)
            println("Compression type: ", ctype, " with: ", cn )
            println("Uncompressed memory size: ", round(mem1, digits=3)," ", mem[2] )
            println("Comment: ", comment)
            println()
            println()
        end
    else

        group = file[datatype]
        output[datatype] = typeof(read(group))
        if verbose
            println("Further content: ", datatype)
            println("-----------")
            println()
            println()
        end
    end

    return output
end

function known_datatype(datatype::String)
    knowntypes = ["hydro", "particles", "clumps"]
    if in(datatype, knowntypes)
        return true
    else
        return false
    end
end


text:
export sample_psdd

using StatsFuns
using LogicCircuits

"Samples an element from a Binomial distribution with p=0.5."
function sample_row(n::Int)::Int
    P = Float64[binompdf(n, 0.5, k) for k ∈ 1:n]
    c = 0.0
    r = rand()
    for (i, p) ∈ enumerate(P)
        c += p
        if r < c return i end
    end
    return n
end

"Samples a single combination uniformly following Robert Floyd's algorithm."
function sample_comb(n::Int, k::Int)::BitSet
    S = BitSet()
    for i ∈ n-k+1:n
        r = rand(1:i)
        if r ∉ S push!(S, r)
        else push!(S, i) end
    end
    return S
end

@inline function random_weights(n::Int)::Vector{Float64}
    W = rand(Float64, n)
    return log.(W /= sum(W))
end

@inline function get_lit(l::Int32, V::Vtree, L::Dict{Int32, StructProbLiteralNode})::StructProbLiteralNode
    if !haskey(L, l)
        node = StructProbLiteralNode(l, V)
        L[l] = node
        return node
    end
    return L[l]
end

"Options for SamplePSDD."
@enum SamplingOpts begin
    plain    = 0  # No compressions or merges.
    compress = 1  # Only compressions are allowed.
    merge    = 2  # Only merges are allowed.
    full     = 3  # Compressions and merges are allowed.
end

"""
Samples a partial partition.
"""
function sample_partition(ϕ::Bdd, Sc::BitSet, p::Real, k::Integer, ⊤_k::Integer, exact::Bool)::Dict{Bdd, Vector{Bdd}}
    X = intersect!(scope(ϕ), Sc)
    idem = isempty(X) || is_⊤(ϕ)
    O = shuffle!(idem ? collect(Sc) : X)
    E = Dict{Bdd, Vector{Bdd}}()
    if idem
        E[ϕ] = Vector{Bdd}()
        sample_idem_primes!(O, p, E[ϕ], ⊤_k)
    else sample_primes!(ϕ, O, E, k, exact) end
    return E
end

"""
Samples primes for the ⊤ case.
"""
function sample_idem_primes!(O::Vector{Int}, p::Real, P::Vector{Bdd}, k::Integer)
    e_count, Sc_len = 1, length(O)
    Q = Tuple{Int, Vector{Int}}[(1, Vector{Int}())]
    while !isempty(Q)
        i, V = popfirst!(Q)
        if e_count >= k || i > Sc_len
            push!(P, and(V))
            continue
        end
        x = O[i]
        c = rand() > p
        if i > 1 && c
            push!(Q, (i+1, V))
        else
            push!(Q, (i+1, push!(copy(V), x)))
            push!(Q, (i+1, push!(copy(V), -x)))
            e_count += 1
        end
    end
    nothing
end

"""
Samples primes for partial partition.
"""
function sample_primes!(ϕ::Bdd, O::Vector{Int}, E::Dict{Bdd, Vector{Bdd}}, k::Integer,
        exact::Bool)
    e_count, Sc_len = 1, length(O)
    Q = Tuple{Bdd, Int, Vector{Int}}[(ϕ, 1, Vector{Int}())]
    while !isempty(Q)
        ψ, i, V = popfirst!(Q)
        if (e_count >= k && !exact) || (i > Sc_len) || is_⊤(ψ)
            if !haskey(E, ψ) E[ψ] = Bdd[and(V)]
            else push!(E[ψ], and(V)) end
            continue
        end
        x = O[i]
        if x ∉ ψ push!(Q, (ψ, i+1, V))
        else
            α, β = ψ|x, ψ|-x
            if !is_⊥(α) push!(Q, (α, i+1, push!(copy(V), x))) end
            if !is_⊥(β) push!(Q, (β, i+1, push!(copy(V), -x))) end
            e_count += 1
        end
    end
    nothing
end

"""
Samples a PSDD from a BDD `ϕ` and vtree `V` with at most `k` elements in each disjunction node.
"""
@inline function sample_psdd(ϕ::Bdd, V::Vtree, k::Integer, D::DataFrame; opts::SamplingOpts = full,
        randomize_weights::Bool = false, pseudocount::Real = 1.0, fact_on_⊤::Bool = false,
        ⊤_k::Integer = k, p_mr::Real = 0.5, always_compress::Bool = false,
        always_merge::Bool = false, merge_branch::Real = 0.0, maxiter::Integer = 0)::StructProbCircuit
    memo = Dict{Tuple{Vtree, Bdd}, StructSumNode}()
    C = sample_psdd_r(ϕ, V, k, Dict{Int32, StructProbLiteralNode}(), randomize_weights, opts,
                      fact_on_⊤, ⊤_k, p_mr, always_compress, always_merge, memo, merge_branch > 0.0,
                      merge_branch, false, false)
    if maxiter > 0
        # Optionally grow the circuit by Strudel.
        loss(x) = heuristic_loss(x, D)
        C = struct_learn(C; primitives = [split_step], kwargs = Dict(split_step => (loss = loss,)),
                         maxiter, verbose = false)
    end
    !randomize_weights && estimate_parameters!(C, D; pseudocount)
    return C
end

function sample_psdd_r(ϕ::Bdd, V::Vtree, k::Integer, leaves::Dict{Int32, StructProbLiteralNode},
        randomize_weights::Bool, opts::SamplingOpts, fact_on_⊤::Bool, ⊤_k::Integer, p_mr::Real,
        always_compress::Bool, always_merge::Bool, repeats::Dict{Tuple{Vtree, Bdd}, StructSumNode},
        merge_this::Bool, merge_branch_pr::Float64, exact::Bool, anc_exact::Bool)::StructProbCircuit
    merge_branch = merge_branch_pr > 0.0
    if merge_branch
        r_p = (V, ϕ)
        if merge_this && (merge_branch_pr > rand()) && haskey(repeats, r_p) return repeats[r_p] end
    end
    if isleaf(V)
        v = convert(Int32, V.var)
        if v ∈ ϕ
            if is_lit(ϕ) return get_lit(to_lit(ϕ), V, leaves) end
            if is_⊤(ϕ|v) return get_lit(v, V, leaves) end
            if is_⊤(ϕ|-v) return get_lit(-v, V, leaves) end
            S = StructSumNode([get_lit(v, V, leaves), get_lit(-v, V, leaves)], V)
            if merge_branch repeats[r_p] = S end
            if randomize_weights S.log_probs = random_weights(2) end
            return S
        end
        S = StructSumNode([get_lit(v, V, leaves), get_lit(-v, V, leaves)], V)
        if merge_branch repeats[r_p] = S end
        if randomize_weights S.log_probs = random_weights(2) end
        return S
    elseif fact_on_⊤ && (is_⊤(ϕ) || isempty(intersect!(scopeset(ϕ), variables(V))))
        # When ϕ ≡ ⊤ and we want to simplify the circuit, fully factorize.
        left = sample_psdd_r(⊤, V.left, k, leaves, randomize_weights, opts, fact_on_⊤, ⊤_k, p_mr,
                             always_compress, always_merge, repeats, merge_branch, merge_branch_pr,
                             false, false)
        right = sample_psdd_r(⊤, V.right, k, leaves, randomize_weights, opts, fact_on_⊤, ⊤_k, p_mr,
                              always_compress, always_merge, repeats, merge_branch,
                              merge_branch_pr, false, false)
        S = StructSumNode([StructMulNode(left, right, V)], V)
        if merge_branch repeats[r_p] = S end
        return S
    elseif length(V.variables) == 2 # when |Sc(a,b)|=2, we have only two elements: (a,ϕ|a) and (¬a,ϕ|¬a).
        prime_var = convert(Int, V.left.var)
        C = Vector{StructProbCircuit}()
        # Left element.
        left_sub_ϕ = ϕ|prime_var
        if !is_⊥(left_sub_ϕ)
            left_prime = sample_psdd_r(bdd_var(prime_var), V.left, k, leaves,
                                       randomize_weights, opts, fact_on_⊤, ⊤_k, p_mr,
                                       always_compress, always_merge, repeats, merge_branch,
                                       merge_branch_pr, true, true)
            left_sub = sample_psdd_r(left_sub_ϕ, V.right, k, leaves, randomize_weights, opts,
                                     fact_on_⊤, ⊤_k, p_mr, always_compress, always_merge, repeats,
                                     merge_branch, merge_branch_pr, true, true)
            push!(C, StructMulNode(left_prime, left_sub, V))
        end
        right_sub_ϕ = ϕ|-prime_var
        if !is_⊥(right_sub_ϕ)
            right_prime = sample_psdd_r(bdd_var(-prime_var), V.left, k, leaves,
                                        randomize_weights, opts, fact_on_⊤, ⊤_k, p_mr,
                                        always_compress, always_merge, repeats, merge_branch,
                                        merge_branch_pr, true, true)
            right_sub = sample_psdd_r(right_sub_ϕ, V.right, k, leaves, randomize_weights, opts,
                                      fact_on_⊤, ⊤_k, p_mr, always_compress, always_merge, repeats,
                                      merge_branch, merge_branch_pr, true, true)
            push!(C, StructMulNode(right_prime, right_sub, V))
        end
        S = StructSumNode(C, V)
        if merge_branch repeats[r_p] = S end
        if randomize_weights S.log_probs = random_weights(length(C)) end
        return S
    end
    L, R = variables(V.left), variables(V.right)
    force_exact = (exact || anc_exact) && !is_⊤(ϕ)
    E = sample_partition(ϕ, L, p_mr, k, ⊤_k, force_exact)
    C = Vector{StructProbCircuit}()
    for (s, P) ∈ E
        # Single element with sub s.
        if length(P) == 1
            # Create (p, s) and move on.
            p = first(P)
            prime_node = sample_psdd_r(p, V.left, k, leaves, randomize_weights, opts, fact_on_⊤,
                                       ⊤_k, p_mr, always_compress, always_merge, repeats,
                                       merge_branch, merge_branch_pr, true, true)
            sub_node = sample_psdd_r(s, V.right, k, leaves, randomize_weights, opts, fact_on_⊤,
                                     ⊤_k, p_mr, always_compress, always_merge, repeats,
                                     merge_branch && isempty(intersect!(scopeset(s), R)),
                                     merge_branch_pr, false, anc_exact)
            push!(C, StructMulNode(prime_node, sub_node, V))
            continue
        end
        # Else, there are length(P) elements with sub s.
        if opts == plain # no compressions nor merges.
            for p ∈ P
                prime_node = sample_psdd_r(p, V.left, k, leaves, randomize_weights, opts, fact_on_⊤,
                                           ⊤_k, p_mr, always_compress, always_merge, repeats,
                                           merge_branch, merge_branch_pr, true, true)
                sub_node = sample_psdd_r(s, V.right, k, leaves, randomize_weights, opts, fact_on_⊤,
                                         ⊤_k, p_mr, always_compress, always_merge, repeats,
                                         merge_branch && isempty(intersect!(scopeset(s), R)),
                                         merge_branch_pr, false, anc_exact)
                push!(C, StructMulNode(prime_node, sub_node, V))
            end
            continue
        end
        n = length(P)
        do_compress = Int(opts) & Int(compress) > 0
        if do_compress
            if always_compress
                K = Tuple{Bdd, Vector{Int}}[(reduce(∨, P), Int[1])]
            else
                c = sample_row(n)
                if c == n # Compress everyone.
                    K = Tuple{Bdd, Vector{Int}}[(reduce(∨, P), Int[1])]
                elseif c == 1 # Compress no one.
                    K = Tuple{Bdd, Vector{Int}}[(P[i], Int[i]) for i ∈ 1:n]
                else # Compress c randomly selected elements.
                    K = Vector{Tuple{Bdd, Vector{Int}}}(undef, n-c+1)
                    comb = sample_comb(n, c)
                    j, ψ = 2, ⊥
                    @inbounds for i ∈ 1:n
                        if i ∈ comb
                            ψ = ψ ∨ P[i]
                        else
                            K[j] = (P[i], Int[j])
                            j += 1
                        end
                    end
                    K[1] = (ψ, Int[1])
                end
            end
        else
            K = Vector{Tuple{Bdd, Vector{Int}}}(undef, n)
            @inbounds for i ∈ 1:n K[i] = (P[i], Int[i]) end
        end
        m = length(K)
        if (Int(opts) & Int(merge)) > 0 # then merge.
            c = sample_row(m)
            if always_merge || c == m
                M = Vector{Int}(undef, m)
                @inbounds for i ∈ 2:m
                    M[i] = i
                    K[i] = (K[i][1], Int[])
                end
                M[1] = 1
                K[1] = (K[1][1], M)
            elseif c != 1 # if c == 1, no need to merge.
                comb = sample_comb(m, c)
                M = Vector{Int}(undef, c)
                j = 1
                for i ∈ 1:m
                    if i ∈ comb
                        M[j] = i
                        if j == 1 K[i] = (K[i][1], M)
                        else K[i] = (K[i][1], Int[]) end
                        j += 1
                    end
                end
            end
        end
        # Create subs.
        subs = Vector{StructProbCircuit}(undef, m)
        for (i, c) ∈ enumerate(K)
            M = c[2]
            if !isempty(M)
                sub_node = sample_psdd_r(s, V.right, k, leaves, randomize_weights, opts, fact_on_⊤,
                                         ⊤_k, p_mr, always_compress, always_merge, repeats,
                                         merge_branch && isempty(intersect!(scopeset(s), R)),
                                         merge_branch_pr, false, anc_exact)
                for j ∈ M subs[j] = sub_node end
            end
        end
        for (i, c) ∈ enumerate(K)
            α, M = c[1], c[2]
            prime_node = sample_psdd_r(α, V.left, k, leaves, randomize_weights, opts, fact_on_⊤,
                                       ⊤_k, p_mr, always_compress, always_merge, repeats,
                                       merge_branch, merge_branch_pr, true, true)
            prod_node = StructMulNode(prime_node, subs[i], V)
            push!(C, prod_node)
        end
    end
    S = StructSumNode(C, V)
    if merge_branch repeats[r_p] = S end
    if randomize_weights S.log_probs = random_weights(length(C)) end
    return S
end


text:
abstract type Phone end

struct DeskPhone <: Phone
    book::Dict{String,String}
end

abstract type Camera end

struct kodak
    roll::Vector{Array{Int32,2}}
end

struct CellPhone <: Phone
    book::Dict{String,String}
    roll::Vector{AbstractVector}
end

function dialnumber(phone::CellPhone)
    println("beep beep")
end

function dialnumber(phone::Phone)
    println("tat tat tat tat")
end

function snap(camera, img)
    println("click")
    push!(roll, img)
end

dphone = DeskPhone(Dict(["information" => "411"]))
cphone = CellPhone(Dict(["emergency" => "911"]), [[]])

dialnumber(dphone)
dialnumber(cphone)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using cascadeddecay, QuantumOptics, meanfield

const γ = 1.
const e_dipole = [0,0,1.]
const T = [0:0.05:10]

#atoms = Atom[Atom(position, e_dipole) for position=cube(0.627)]
#atoms = Atom[Atom(position, e_dipole) for position=cube(0.2)]
atoms = Atom[Atom(position, e_dipole) for position=chain(0.4,30)]
system = CascadedDecaySystem(atoms, γ)

# Meanfield timeevolution
N = length(system.atoms)
sz0 = zeros(ComplexF64, N)
sp0 = zeros(ComplexF64, N) .- 0.5

println("meanfield +z")
@time T_mf, sz_mf, sp_mf = meanfield_timeevolution(T, system, sz0, sp0)
@time T_mf, sz_mf, sp_mf = meanfield_timeevolution(T, system, sz0, sp0)

# Meanfield xyz timeevolution
N = length(system.atoms)
sx0 = zeros(Float64, N) .- 1
sy0 = zeros(Float64, N)
sz0 = zeros(Float64, N)

println("meanfield xyz")
@time T_xmf, sx_xmf, sy_xmf, sz_xmf= meanfield.meanfield_timeevolution2(T, system, sx0, sy0, sz0)
@time T_xmf, sx_xmf, sy_xmf, sz_xmf= meanfield.meanfield_timeevolution2(T, system, sx0, sy0, sz0)

# Correlation included timeevolution
N = length(system.atoms)
sz0 = zeros(ComplexF64, N)
sp0 = zeros(ComplexF64, N) .- 0.5
Cpm0 = zeros(ComplexF64, N, N)
Cpz0 = zeros(ComplexF64, N, N)
Cpp0 = zeros(ComplexF64, N, N)
Czz0 = zeros(ComplexF64, N, N)

println("correlation +z")
@time T_cor, sz_cor, sp_cor = correlation_timeevolution(T, system, sz0, sp0, Cpm0, Cpz0, Cpp0, Czz0)
@time T_cor, sz_cor, sp_cor = correlation_timeevolution(T, system, sz0, sp0, Cpm0, Cpz0, Cpp0, Czz0)

# Correlation included timeevolution xyz
N = length(system.atoms)
sx0 = zeros(Float64, N) .- 1
sy0 = zeros(Float64, N)
sz0 = zeros(Float64, N)
C0 = Dict{AbstractString, Matrix{Float64}}()
C0["xx"] = zeros(Float64, N, N) .+ 1
C0["yy"] = zeros(Float64, N, N)
C0["zz"] = zeros(Float64, N, N)
C0["xy"] = zeros(Float64, N, N)
C0["xz"] = zeros(Float64, N, N)
C0["yz"] = zeros(Float64, N, N)

#@time T_xcor, sx_xcor, sy_xcor, sz_xcor, C_xcor = meanfield.correlation_timeevolution2(T, system, sx0, sy0, sz0, C0)
#@time T_xcor, sx_xcor, sy_xcor, sz_xcor, C_xcor = meanfield.correlation_timeevolution2(T, system, sx0, sy0, sz0, C0)

println("correlation xyz")
@time T_xcor, sx_xcor, sy_xcor, sz_xcor, C_xcor = meanfield.correlation_timeevolution3(T, system, sx0, sy0, sz0, C0)
@time T_xcor, sx_xcor, sy_xcor, sz_xcor, C_xcor = meanfield.correlation_timeevolution3(T, system, sx0, sy0, sz0, C0)

println(norm(sz_xcor[end]-sz_cor[end]))


text:
getoutputlength(version::Int32, x::Integer) = version < 2 ? x : padding(x)

function validatedata(data::AbstractVector{UInt8})
    if length(data) < MIN_FILE_LENGTH
        throw(ArgumentError("Data is not in feather format: total length of buffer: $(length(data))"))
    end
    header = data[1:4]
    footer = data[(end-3):end]
    if header ≠ FEATHER_MAGIC_BYTES || footer ≠ FEATHER_MAGIC_BYTES
        throw(ArgumentError(string("Data is not in feather format: header = $header, ",
                                   "footer = $footer.")))
    end
    data
end

loaddata(data::AbstractVector{UInt8}; use_mmap::Bool=true) = validatedata(data)
loaddata(io::IOBuffer; use_mmap::Bool=true) = validatedata(io.data)
loaddata(io::IO; use_mmap::Bool=true) = Base.read(io)
function loaddata(filename::AbstractString; use_mmap::Bool=SHOULD_USE_MMAP)
    isfile(filename) || throw(ArgumentError("'$filename' does not exist."))
    data = SHOULD_USE_MMAP ? Mmap.mmap(filename) : Base.read(filename)
    loaddata(data)
end

function metalength(data::AbstractVector{UInt8})
    Base.read(IOBuffer(data[(length(data)-7):(length(data)-4)]), Int32)
end

function metaposition(data::AbstractVector{UInt8}, metalen::Integer=metalength(data))
    length(data) - (metalen+7)
end

function rootposition(data::AbstractVector{UInt8}, mpos::Integer=metaposition(data))
    Base.read(IOBuffer(data[mpos:(mpos+4)]), Int32)
end

function getctable(data::AbstractVector{UInt8})
    metapos = metaposition(data)
    rootpos = rootposition(data, metapos)
    ctable = FlatBuffers.read(Metadata.CTable, data, metapos + rootpos - 1)
    if ctable.version < FEATHER_VERSION
        @warn("This feather file is old and may not be readable.")
    end
    ctable
end


text:
using BioMedQuery
using DataFrames

"""
    get_semantic_occurrences_df(mesh, umls_concepts...)

Given a mesh descriptor and filtering UMLS concepts, returns a DataFrame with
PMIDs and related mesh descriptor filtered by the concepts.
"""
function get_semantic_occurrences_df(mesh::String, umls_concepts::String...)

    concept_string = "'" * join(umls_concepts,"','") * "'"

    # check for parentheses in mesh and replace with spaces per docs

    articles = BioMedQuery.Processes.pubmed_search_and_parse("", """("$mesh"[mh:noexp])""", 999999999999) #"""("$mesh"[mh:noexp])"""

    # # select pmids from medline which use the mesh descriptor
    # query_string = """ SELECT mh.pmid, mrc.str as descriptor
    #                 FROM medline.mesh_heading mh
    #
    #                 JOIN (select uid from medline.mesh_desc where name = '$mesh') md
    #                 ON md.uid <> mh.desc_uid
    #
    #                 JOIN pubmed_comorbidities.ALL_MESH mrc
    #                 ON mrc.uid = mh.desc_uid
    #
    #                 JOIN umls_meta.MRSTY mrs
    #                 ON mrc.cui = mrs.cui
    #
    #                 JOIN medline.mesh_heading mh2
    #                 ON mh2.pmid = mh.pmid
    #
    #                 join (select uid from medline.mesh_desc where name = '$mesh') md2
    #                 on md2.uid = mh2.desc_uid
    #
    #                 WHERE mrs.sty in ($concept_string); """;
    #
    # articles_df = MySQL.query(db, query_string, DataFrame)
    #
    # unique!(articles_df)

    return articles

end


text:
# This file was generated, do not modify it. # hide
using Gnuplot, Random
Random.seed!(123)
goffset = "set offsets graph .05, graph .05, graph .05, graph .05"
@gp rand(30) rand(30) "w p pt 4 ps 3 lc '#0072B2' lw 2 t 'marker'" goffset
save(term="pngcairo size 600,400", output="plt_ex10.png") # hide


text:
# ---
# title: 1380. Lucky Numbers in a Matrix
# id: problem1380
# author: zhwang
# date: 2022-02-15
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/lucky-numbers-in-a-matrix/description/>
# hidden: true
# ---
#
# Given a `m * n` matrix of **distinct** numbers, return all lucky numbers in
# the matrix in **any** order.
#
# A lucky number is an element of the matrix such that it is the minimum element
# in its row and maximum in its column.
#
#
#
# **Example 1:**
#
#
#
#     Input: matrix = [[3,7,8],[9,11,13],[15,16,17]]
#     Output: [15]
#     Explanation: 15 is the only lucky number since it is the minimum in its row and the maximum in its column
#
#
# **Example 2:**
#
#
#
#     Input: matrix = [[1,10,4,2],[9,3,8,7],[15,16,17,12]]
#     Output: [12]
#     Explanation: 12 is the only lucky number since it is the minimum in its row and the maximum in its column.
#
#
# **Example 3:**
#
#
#
#     Input: matrix = [[7,8],[1,2]]
#     Output: [7]
#
#
#
#
# **Constraints:**
#
#   * `m == mat.length`
#   * `n == mat[i].length`
#   * `1 <= n, m <= 50`
#   * `1 <= matrix[i][j] <= 10^5`.
#   * All elements in the matrix are distinct.
#
#
## @lc code=start
using LeetCode
function lucky_numbers(matrix::Vector{Vector{Int}})::Vector{Int}
    m, n = length(matrix), length(matrix[1])
    row_min = Dict{Int,Tuple{Int,Int}}(i=>(matrix[i][1], 1) for i in 1:m)
    col_max = Dict{Int,Tuple{Int,Int}}(j=>(0,0) for j in 1:n)
    for i in 1:m
        for j in 1:n
            c = matrix[i][j]
            if c < row_min[i][1]
                row_min[i] = (c, j)
            end
            if c > col_max[j][1]
                col_max[j] = (c, i)
            end
        end
    end
    return [row_min[i][1] for i in 1:m if col_max[row_min[i][2]][2] == i]
end
## @lc code=end

## @lc test=start
@testset "1380.lucky-numbers-in-a-matrix.jl" begin
    @test lucky_numbers([[3,7,8],[9,11,13],[15,16,17]]) == [15]
    @test lucky_numbers([[1,10,4,2],[9,3,8,7],[15,16,17,12]]) == [12]
    @test lucky_numbers([[7,8],[1,2]]) == [7]
end
## @lc test=end

## @lc info=start
# link: [solution to 1380](https://leetcode-cn.com/problems/lucky-numbers-in-a-matrix/solution/1380-ju-zhen-zhong-de-xing-yun-shu-julia-ehjx/)
## @lc info=end


text:
# Formulas for representing and working with linear-model-type expressions
# Original by Harlan D. Harris.  Later modifications by John Myles White
# and Douglas M. Bates.

## Formulas are written as expressions and parsed by the Julia parser.
## For example :(y ~ a + b + log(c))
## In Julia the & operator is used for an interaction.  What would be written
## in R as y ~ a + b + a:b is written :(y ~ a + b + a&b) in Julia.
## The equivalent R expression, y ~ a*b, is the same in Julia

## The lhs of a one-sided formula is 'nothing'
## The rhs of a formula can be 1

type Formula
    lhs::Union(Symbol, Expr,Nothing)
    rhs::Union(Symbol, Expr, Integer)
end

type Terms
    terms::Vector
    eterms::Vector                    # evaluation terms
    factors::Matrix{Int8}             # maps terms to evaluation terms
    order::Vector{Int}                # orders of rhs terms
    response::Bool       # indicator of a response, which is eterms[1] if present
    intercept::Bool      # is there an intercept column in the model matrix?
end

type ModelFrame
    df::AbstractDataFrame
    terms::Terms
    msng::BitArray
end

type ModelMatrix{T<:Union(Float32,Float64)}
    m::Matrix{T}
    assign::Vector{Int}
end

Base.size(mm::ModelMatrix) = size(mm.m)
Base.size(mm::ModelMatrix,dim...) = size(mm.m,dim...)

function Formula(ex::Expr)
    aa = ex.args
    if aa[1] != :~
        error("Invalid formula, top-level argument must be '~'.  Check parentheses.")
    end
    if length(aa) == 2 return Formula(nothing, aa[2]) end
    Formula(aa[2], aa[3])
end

function Base.show(io::IO, f::Formula)
    print(io, string("Formula: ", f.lhs == nothing ? "" : f.lhs, " ~ ", f.rhs))
end

## Return, as a vector of symbols, the names of all the variables in
## an expression or a formula
function allvars(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    [[allvars(a) for a in ex.args[2:]]...]
end
allvars(f::Formula) = unique(vcat(allvars(f.rhs), allvars(f.lhs)))
allvars(sym::Symbol) = [sym]
allvars(v) = Array(Symbol,0)

const specials = Set(:+,:-,:*,:/,:&,:|,:^) # special operators in formulas

function dospecials(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    a1 = ex.args[1]
    if !(a1 in specials) return ex end
    excp = copy(ex)
    excp.args = vcat(a1,map(dospecials, ex.args[2:]))
    if a1 != :* return excp end
    aa = excp.args
    a2 = aa[2]
    a3 = aa[3]
    if length(aa) > 3
        excp.args = vcat(a1, aa[3:])
        a3 = dospecials(excp)
    end
    :($a2 + $a3 + $a2 & $a3)
end
dospecials(a) = a

const associative = Set(:+,:*,:&)       # associative special operators

## If the expression is a call to the function s return its arguments
## Otherwise return the expression
function ex_or_args(ex::Expr,s::Symbol)
    if ex.head != :call error("Non-call expression encountered") end
    excp = copy(ex)
    a1 = ex.args[1]
    a2 = map(condense, ex.args[2:])
    if a1 == s return a2 end
    excp.args = vcat(a1, a2)
    excp
end
ex_or_args(a,s::Symbol) = a

## Condense calls like :(+(a,+(b,c))) to :(+(a,b,c))
## Also need to work out how to distribute & over +
function condense(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    a1 = ex.args[1]
    if !(a1 in associative) return ex end
    excp = copy(ex)
    excp.args = vcat(a1, map(x->ex_or_args(x,a1), ex.args[2:])...)
    excp
end
condense(a) = a

getterms(ex::Expr) = (ex.head == :call && ex.args[1] == :+) ? ex.args[2:] : ex
getterms(a) = a

ord(ex::Expr) = (ex.head == :call && ex.args[1] == :&) ? length(ex.args)-1 : 1
ord(a) = 1

const nonevaluation = Set(:&,:|)        # operators constructed from other evaluations
## evaluation terms - the (filtered) arguments for :& and :|, otherwise the term itself
function evt(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    if !(ex.args[1] in nonevaluation) return ex end
    filter(x->!isa(x,Number), vcat(map(getterms, ex.args[2:])...))
end
evt(a) = {a}

function Terms(f::Formula)
    rhs = condense(dospecials(f.rhs))
    tt = getterms(rhs)
    if !isa(tt,AbstractArray) tt = [tt] end
    tt = tt[!(tt .== 1)]             # drop any explicit 1's
    noint = (tt .== 0) | (tt .== -1) # should also handle :(-(expr,1))
    tt = tt[!noint]
    oo = int(map(ord, tt))           # orders of interaction terms
    if !issorted(oo)                 # sort terms by increasing order
        pp = sortperm(oo)
        tt = tt[pp]
        oo = oo[pp]
    end
    etrms = map(evt, tt)
    haslhs = f.lhs != nothing
    if haslhs
        unshift!(etrms, {f.lhs})
        unshift!(oo, 1)
    end
    ev = unique(vcat(etrms...))
    facs = int8(hcat(map(x->(s=Set(x...);map(t->int8(t in s), ev)),etrms)...))
    Terms(tt, ev, facs, oo, haslhs, !any(noint))
end

## Default NA handler.  Others can be added as keyword arguments
function na_omit(df::DataFrame)
    cc = complete_cases(df)
    df[cc,:], cc
end

## Trim the pool field of da to only those levels that occur in the refs
function dropUnusedLevels!(da::PooledDataArray)
    rr = da.refs
    uu = unique(rr)
    if length(uu) == length(da.pool) return da end
    T = eltype(rr)
    su = sort!(uu)
    dict = Dict(su, one(T):convert(T,length(uu)))
    da.refs = map(x->dict[x], rr)
    da.pool = da.pool[uu]
    da
end
dropUnusedLevels!(x) = x

function ModelFrame(f::Formula, d::AbstractDataFrame)
    trms = Terms(f)
    df,msng = na_omit(DataFrame(map(x->with(d,x),trms.eterms)))
    names!(df, map(string, trms.eterms))
    for c in df dropUnusedLevels!(c[2]) end
    ModelFrame(df, trms, msng)
end
ModelFrame(ex::Expr, d::AbstractDataFrame) = ModelFrame(Formula(ex), d)

function model_response(mf::ModelFrame)
    mf.terms.response || error("Model formula one-sided")
    array(mf.df[bool(mf.terms.factors[:,1])][:,1])
end

function contr_treatment(n::Integer, contrasts::Bool, sparse::Bool, base::Integer)
    if n < 2 error("not enought degrees of freedom to define contrasts") end
    contr = sparse ? speye(n) : eye(n) .== 1.
    if !contrasts return contr end
    if !(1 <= base <= n) error("base = $base is not allowed for n = $n") end
    contr[:,vcat(1:(base-1),(base+1):end)]
end
contr_treatment(n::Integer,contrasts::Bool,sparse::Bool) = contr_treatment(n,contrasts,sparse,1)
contr_treatment(n::Integer,contrasts::Bool) = contr_treatment(n,contrasts,false,1)
contr_treatment(n::Integer) = contr_treatment(n,true,false,1)
cols(v::PooledDataVector) = contr_treatment(length(v.pool))[v.refs,:]
cols(v::DataVector) = reshape(float64(v.data), (length(v),1))

function isfe(ex::Expr)                 # true for fixed-effects terms
    if ex.head != :call error("Non-call expression encountered") end
    ex.args[1] != :|
end
isfe(a) = true

## Expand the columns in an interaction term
function expandcols(trm::Vector)
    if length(trm) == 1 return float64(trm[1]) end
    if length(trm) == 2
        a = float64(trm[1])
        b = float64(trm[2])
        nca = size(a,2)
        ncb = size(b,2)
        return hcat([a[:,i].*b[:,j] for i in 1:nca, j in 1:ncb]...)
    end
    error("code for 3rd and higher order interactions not yet written")
end

nc(trm::Vector) = *([size(x,2) for x in trm]...)

function ModelMatrix(mf::ModelFrame)
    trms = mf.terms
    aa = {{ones(size(mf.df,1),int(trms.intercept))}}
    asgn = zeros(Int, (int(trms.intercept)))
    fetrms = bool(map(isfe, trms.terms))
    if trms.response unshift!(fetrms,false) end
    ff = trms.factors[:,fetrms]
    ## need to be cautious here to avoid evaluating cols for a factor with many levels
    ## if the factor doesn't occur in the fetrms
    rows = vec(bool(sum(ff,[2])))
    ff = ff[rows,:]
    cc = [cols(x[2]) for x in mf.df[:,rows]]
    for j in 1:size(ff,2)
        trm = cc[bool(ff[:,j])]
        push!(aa, trm)
        asgn = vcat(asgn, fill(j, nc(trm)))
    end
    ModelMatrix{Float64}(hcat([expandcols(t) for t in aa]...), asgn)
end

function coefnames(fr::ModelFrame)
    if fr.terms.intercept
        vnames = UTF8String["(Intercept)"]
    else
        vnames = UTF8String[]
    end
    # Need to only include active levels
    for term in fr.terms.terms
        if isa(fr.df[term], PooledDataArray)
            for lev in levels(fr.df[term])[2:end]
                push!(vnames, string(term, " - ", lev))
            end
        else
            push!(vnames, string(term))
        end
    end
    return vnames
end

# Expand dummy variables and equations
## function model_matrix(mf::ModelFrame)
##     ex = mf.formula.rhs[1]
##     # BUG: complete_cases doesn't preserve grouped columns
##     df = mf.df#[complete_cases(mf.df),1:ncol(mf.df)]
##     rdf = df[mf.y_indexes]
##     mdf = expand(ex, df)

##     # TODO: Convert to Array{Float64} in a cleaner way
##     rnames = colnames(rdf)
##     mnames = colnames(mdf)
##     r = Array(Float64,nrow(rdf),ncol(rdf))
##     m = Array(Float64,nrow(mdf),ncol(mdf))
##     for i = 1:nrow(rdf)
##       for j = 1:ncol(rdf)
##         r[i,j] = float(rdf[i,j])
##       end
##       for j = 1:ncol(mdf)
##         m[i,j] = float(mdf[i,j])
##       end
##     end

##     include_intercept = true
##     if include_intercept
##       m = hcat(ones(nrow(mdf)), m)
##       unshift!(mnames, "(Intercept)")
##     end

##     ModelMatrix(m, r, mnames, rnames)
##     ## mnames = {}
##     ## rnames = {}
##     ## for c in 1:ncol(rdf)
##     ##   r = hcat(r, float(rdf[c]))
##     ##   push!(rnames, colnames(rdf)[c])
##     ## end
##     ## for c in 1:ncol(mdf)
##     ##   m = hcat(m, mdf[c])
##     ##   push!(mnames, colnames(mdf)[c])
##     ## end
## end

## model_matrix(f::Formula, d::AbstractDataFrame) = model_matrix(ModelFrame(f, d))
## model_matrix(ex::Expr, d::AbstractDataFrame) = model_matrix(ModelFrame(Formula(ex), d))

## # TODO: Make a more general version of these functions
## # TODO: Be able to extract information about each column name
## function interaction_design_matrix(a::AbstractDataFrame, b::AbstractDataFrame)
##    cols = {}
##    col_names = Array(ASCIIString,0)
##    for i in 1:ncol(a)
##        for j in 1:ncol(b)
##           push!(cols, DataArray(a[:,i] .* b[:,j]))
##           push!(col_names, string(colnames(a)[i],"&",colnames(b)[j]))
##        end
##    end
##    DataFrame(cols, col_names)
## end

## function interaction_design_matrix(a::AbstractDataFrame, b::AbstractDataFrame, c::AbstractDataFrame)
##    cols = {}
##    col_names = Array(ASCIIString,0)
##    for i in 1:ncol(a)
##        for j in 1:ncol(b)
##            for k in 1:ncol(b)
##               push!(cols, DataArray(a[:,i] .* b[:,j] .* c[:,k]))
##               push!(col_names, string(colnames(a)[i],"&",colnames(b)[j],"&",colnames(c)[k]))
##            end
##        end
##    end
##    DataFrame(cols, col_names)
## end

## # Temporary: Manually describe the interactions needed for DataFrame Array.
## function all_interactions(dfs::Array{Any,1})
##     d = DataFrame()
##     if length(dfs) == 2
##       combos = ([1,2],)
##     elseif length(dfs) == 3
##       combos = ([1,2], [1,3], [2,3], [1,2,3])
##     else
##       error("interactions with more than 3 terms not implemented (yet)")
##     end
##     for combo in combos
##        if length(combo) == 2
##          a = interaction_design_matrix(dfs[combo[1]],dfs[combo[2]])
##        elseif length(combo) == 3
##          a = interaction_design_matrix(dfs[combo[1]],dfs[combo[2]],dfs[combo[3]])
##        end
##        d = insert!(d, a)
##     end
##     return d
## end

## # string(Expr) now quotes, which we don't want. This hacks around that, stealing
## # from print_to_string
## function formula_string(ex::Expr)
##     s = memio(0, false)
##     Base.show_unquoted(s, ex)
##     takebuf_string(s)
## end

## #
## # The main expression to DataFrame expansion function.
## # Returns a DataFrame.
## #

## function expand(ex::Expr, df::AbstractDataFrame)
##     f = eval(ex.args[1])
##     if method_exists(f, (FormulaExpander, Vector{Any}, DataFrame))
##         # These are specialized expander functions (+, *, &, etc.)
##         f(FormulaExpander(), ex.args[2:end], df)
##     else
##         # Everything else is called recursively:
##         expand(with(df, ex), formula_string(ex), df)
##     end
## end

## function expand(s::Symbol, df::AbstractDataFrame)
##     expand(with(df, s), string(s), df)
## end

## # TODO: make this array{symbol}?
## function expand(args::Array{Any}, df::AbstractDataFrame)
##     [expand(x, df) for x in args]
## end

## function expand(x, name::ByteString, df::AbstractDataFrame)
##     # If this happens to be a column group, then expand each and concatenate
##     if is_group(df, name)
##       preds = get_groups(df)[name]
##       dfs = [expand(symbol(x), df) for x in preds]
##       return cbind(dfs...)
##     end
##     # This is the default for expansion: put it right in to a DataFrame.
##     DataFrame({x}, [name])
## end

## #
## # Methods for expansion of specific data types
## #

## # Expand a PooledDataVector into a matrix of indicators for each dummy variable
## # TODO: account for NAs?
## function expand(poolcol::PooledDataVector, colname::ByteString, df::AbstractDataFrame)
##     newcol = {DataArray([convert(Float64,x)::Float64 for x in (poolcol.refs .== i)]) for i in 2:length(poolcol.pool)}
##     newcolname = [string(colname, ":", x) for x in poolcol.pool[2:length(poolcol.pool)]]
##     DataFrame(newcol, convert(Vector{ByteString}, newcolname))
## end

## #
## # Methods for Formula expansion
## #
## type FormulaExpander; end # This is an indictor type.

## function +(::FormulaExpander, args::Vector{Any}, df::AbstractDataFrame)
##     d = DataFrame()
##     for a in args
##         d = insert!(d, expand(a, df))
##     end
##     d
## end
## function (&)(::FormulaExpander, args::Vector{Any}, df::AbstractDataFrame)
##     interaction_design_matrix(expand(args[1], df), expand(args[2], df))
## end
## function *(::FormulaExpander, args::Vector{Any}, df::AbstractDataFrame)
##     d = +(FormulaExpander(), args, df)
##     d = insert!(d, all_interactions(expand(args, df)))
##     d
## end


text:
module PythonBufferInterface

import ArrayInterface

export pybuf_info, pybuf_buffer, pybuf_ptr, pybuf_readonly, pybuf_itemsize, pybuf_format, pybuf_ndim, pybuf_shape, pybuf_stride, pybuf_strides, pybuf_suboffset, pybuf_suboffsets, pybuf_len

const TYPE_TO_TYPESTR = Dict(
    Cchar => "b",
    Cuchar => "B",
    Cshort => "h",
    Cushort => "H",
    Cint => "i",
    Cuint => "I",
    Clong => "l",
    Culong => "L",
    Clonglong => "q",
    Culonglong => "Q",
    Float16 => "e",
    Cfloat => "f",
    Cdouble => "d",
    Complex{Float16} => "Ze",
    Complex{Cfloat} => "Zf",
    Complex{Cdouble} => "Zd",
    Bool => "?",
    Ptr{Cvoid} => "P",
)

### Interface

"""
    pybuf_buffer(x)

If `x` does not satisfy the Python buffer interface, return `nothing`. Otherwise return an
object `b` which can be queried for information using functions such as `pybuf_ptr`,
`pybuf_ndim`, etc.
"""
pybuf_buffer(x) = nothing
pybuf_buffer(x::AbstractArray) = is_pybuf_eltype(eltype(x)) && ArrayInterface.defines_strides(typeof(x)) && ArrayInterface.device(typeof(x)) == ArrayInterface.CPUPointer() ? x : nothing

"""
    pybuf_ptr(b)

Pointer to the start of the data in buffer `b`.
"""
pybuf_ptr(b) = Ptr{Cvoid}(Base.unsafe_convert(Ptr{eltype(b)}, b))

"""
    pybuf_readonly(b)

True if the buffer `b` is read only.
"""
pybuf_readonly(b) = !ArrayInterface.can_setindex(typeof(b))

"""
    pybuf_itemsize(b)

The size in bytes of a single item in the buffer `b`.
"""
pybuf_itemsize(b) = Int(sizeof(eltype(b)))

"""
    pybuf_format(b)

The format string describing an element of the buffer `b`.
"""
pybuf_format(b) = pybuf_format(eltype(b))

"""
    pybuf_ndim(b)

The number of dimensions in buffer `b`.
"""
pybuf_ndim(b) = Int(ndims(b))

"""
    pybuf_shape(b, i)

The size of the buffer `b` along dimension `i`.
"""
pybuf_shape(b, i) = Int(size(b, i))

"""
    pybuf_stride(b, i)

The distance in bytes between consecutive elements of buffer `b` along dimension `i`.
"""
pybuf_stride(b, i) = Int(stride(b, i) * Base.aligned_sizeof(eltype(b)))

"""
    pybuf_suboffset(b, i)

The suboffsets in buffer `b` along dimension `i`.
"""
pybuf_suboffset(b, i) = -1

### Utils

"""
    is_pybuf_eltype(T::Type)

Return true if items of type `T` are compatible with the buffer interface.

That is, `T` must be a bits type or a struct of these.
"""
@generated is_pybuf_eltype(::Type{T}) where {T} = (
    isconcretetype(T) &&
    Base.allocatedinline(T) &&
    isa(T, DataType) &&
    (
        isbitstype(T) ||
        (
            isstructtype(T) &&
            all(i->is_pybuf_eltype(fieldtype(T, i)), 1:fieldcount(T))
        )
    )
)

"""
    pybuf_format(T::Type)

The format string describing items of type `T`.
"""
@generated pybuf_format(::Type{T}) where {T} =
    if !is_pybuf_eltype(T)
        error("invalid Python buffer element type")
    elseif haskey(TYPE_TO_TYPESTR, T)
        TYPE_TO_TYPESTR[T]
    elseif isstructtype(T)
        n = fieldcount(T)
        flds = []
        for i = 1:n
            nm = fieldname(T, i)
            tp = fieldtype(T, i)
            push!(flds, string(pybuf_format(tp), nm isa Symbol ? ":$nm:" : ""))
            off0 = fieldoffset(T, i) + sizeof(tp)
            off1 = i == n ? sizeof(T) : fieldoffset(T, i + 1)
            d = off1 - off0
            @assert d ≥ 0
            d > 0 && push!(flds, "$(d)x")
        end
        string("T{", join(flds, " "), "}")
    elseif isbitstype(T)
        "$(sizeof(T))x"
    else
        @assert false
    end

"""
    pybuf_shape(b)

The size of the buffer `b` along each dimension.
"""
pybuf_shape(b) = ntuple(i->pybuf_shape(b, i)::Int, pybuf_ndim(b)::Int)

"""
    pybuf_strides(b)

The distance in bytes between consecutive elements of buffer `b` along each dimension.
"""
pybuf_strides(b) = ntuple(i->pybuf_stride(b, i)::Int, pybuf_ndim(b)::Int)

"""
    pybuf_suboffsets(b)

The suboffsets in buffer `b` along each dimension.
"""
pybuf_suboffsets(b) = ntuple(i->pybuf_suboffset(b, i)::Int, pybuf_ndim(b)::Int)

"""
    pybuf_len(b)

The size of the buffer `b` in bytes if it were stored contiguously.
"""
pybuf_len(b) = (prod(pybuf_shape(b)) * pybuf_itemsize(b))::Int

"""
    pybuf_info(x)

If `x` satisfies the buffer protocol, return a named tuple describing it. Otherwise return
`nothing`.

The named tuple has fields `buffer`, `ptr`, `len`, `readonly`, `itemsize`, `format`,
`ndim`, `shape`, `strides`, `suboffsets` corresponding to the outputs of the `pybuf_*`
functions.

The information in the tuple remains valid as long as the `buffer` is not GC'd.
"""
pybuf_info(x) = if (b = pybuf_buffer(x)) !== nothing
    (
        buffer = b,
        ptr = pybuf_ptr(b)::Ptr{Cvoid},
        len = pybuf_len(b)::Int,
        readonly = pybuf_readonly(b)::Bool,
        itemsize = pybuf_itemsize(b)::Int,
        format = pybuf_format(b)::String,
        ndim = pybuf_ndim(b)::Int,
        shape = pybuf_shape(b)::Tuple{Vararg{Int}},
        strides = pybuf_strides(b)::Tuple{Vararg{Int}},
        suboffsets = pybuf_suboffsets(b)::Tuple{Vararg{Int}},
    )
end

end # module


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

module DirectGaussianSimulation

using GeoStatsBase

using Variography
using LinearAlgebra

import GeoStatsBase: preprocess, solve_single

export DirectGaussSim

"""
    DirectGaussSim(var₁=>param₁, var₂=>param₂, ...)

Direct Gaussian simulation (a.k.a. LU simulation).

## Parameters

* `variogram` - theoretical variogram (default to GaussianVariogram())

### References

Alabert 1987. *The practice of fast conditional simulations through the
LU decomposition of the covariance matrix.*
"""
@simsolver DirectGaussSim begin
  @param variogram = GaussianVariogram()
  @param mean = nothing
end

function preprocess(problem::SimulationProblem, solver::DirectGaussSim)
  # retrieve problem info
  pdata = data(problem)
  pdomain = domain(problem)

  # result of preprocessing
  preproc = Dict{Symbol,Tuple}()

  for (var, V) in variables(problem)
    # get user parameters
    if var ∈ keys(solver.params)
      varparams = solver.params[var]
    else
      varparams = DirectGaussSimParam()
    end

    # determine variogram model
    γ = varparams.variogram

    # check stationarity
    @assert isstationary(γ) "variogram model must be stationary"

    # retrieve data locations in domain and data values
    datalocs = Vector{Int}()
    z₁ = Vector{V}()
    for (loc, dataloc) in datamap(problem, var)
      push!(datalocs, loc)
      push!(z₁, pdata[dataloc,var])
    end

    # retrieve simulation locations
    simlocs = [l for l in 1:npoints(pdomain) if l ∉ datalocs]

    # covariance between simulation locations
    C₂₂ = sill(γ) .- pairwise(γ, pdomain, simlocs)

    if isempty(datalocs)
      d₂  = zero(V)
      L₂₂ = cholesky(Symmetric(C₂₂)).L
    else
      # covariance beween data locations
      C₁₁ = sill(γ) .- pairwise(γ, pdomain, datalocs)
      C₁₂ = sill(γ) .- pairwise(γ, pdomain, datalocs, simlocs)

      L₁₁ = cholesky(Symmetric(C₁₁)).L
      B₁₂ = L₁₁ \ C₁₂
      A₂₁ = B₁₂'

      d₂ = A₂₁ * (L₁₁ \ z₁)
      L₂₂ = cholesky(Symmetric(C₂₂ - A₂₁*B₁₂)).L
    end

    if !isnothing(varparams.mean) && !isempty(datalocs)
      @warn "mean can only be specified in unconditional simulation"
    end

    # mean for unconditional simulation
    μ = isnothing(varparams.mean) ? zero(V) : varparams.mean

    preproc[var] = (z₁, d₂, L₂₂, μ, datalocs, simlocs)
  end

  preproc
end

function solve_single(problem::SimulationProblem, var::Symbol,
                      solver::DirectGaussSim, preproc)
  # retrieve problem info
  pdomain = domain(problem)

  # unpack preprocessed parameters
  z₁, d₂, L₂₂, μ, datalocs, simlocs = preproc[var]

  # allocate memory for result
  realization = Vector{eltype(z₁)}(undef, npoints(pdomain))

  # set hard data
  realization[datalocs] = z₁

  # simulate the rest
  w₂ = randn(size(L₂₂, 2))
  y₂ = d₂ .+ L₂₂*w₂
  realization[simlocs] = y₂

  # adjust mean in case of unconditional simulation
  isempty(datalocs) && (realization .+= μ)

  realization
end

end


text:
using MLStyle
export ErasedType, ERArrow, ERNom, ERTuple, ERAny
export tyErase

@data ErasedType begin
    ERArrow(ErasedType, ErasedType)
    ERNom(Symbol)
    ERTuple(Vector{ErasedType})
    ERType(ErasedType)
    ERAny
end

function tyErase(hm::HMT)
    @match hm begin
        Arrow(a, r) => ERArrow(tyErase(a), tyErase(r))
        Forall(_, p) => tyErase(p)
        Var() => error("type variable should be solved")
        Bound() => ERAny
        Tup(xs) => ERTuple(ErasedType[tyErase(x) for x in xs])
        Nom(x) => ERNom(x)
        App(Nom(:Type), x) => ERType(tyErase(x))
        App() => ERAny
        Implicit(t) => tyErase(t)
    end
end


text:
## MDP definition for the merging problem

const MergingState = Vector{Float64} #SVector{5, Float64} # se, ve, ae, so, vo

"""
    MergingMDP <: MDP{MergingState, Int64}
A discrete state and action MDP representing a merging scenario with only two agents:
the ego vehicle on the merge lane and another traffic participant on the main lane.
"""
@with_kw mutable struct MergingMDP <: MDP{MergingState, Int64}
    env::MergingEnvironment = MergingEnvironment()
    cardef::VehicleDef = CARDEF
    dt::Float64 = 0.2
    min_acc::Float64 = -4.0
    pos_res::Float64 = 2.0*M2L
    vel_res::Float64 = 2.0
    acc_res::Float64 = 1.0
    jerk_levels::SVector{5, Float64} = SVector(-0.2, -0.1, 0, 0.1, 0.2)
    human_actions::Vector{Float64} = [0.0]# [-2.0, 0.0, 2.0]
    ego_lane::Lane = env.roadway[LaneTag(MERGE_LANE_ID, 1)]
    other_lane::Lane = env.roadway[LaneTag(MAIN_LANE_ID, 1)]
    ego_grid::LinRange{Float64} = LinRange(0.0, env.merge_lane_length + env.after_merge_length, 40)
    human_grid::LinRange{Float64} = LinRange(0.0, env.merge_lane_length + env.after_merge_length, 40)
    velocity_grid::LinRange{Float64} = LinRange(0.0, env.main_lane_vmax, 10)
    acceleration_grid::LinRange{Float64} = LinRange(-4.0, 2, 4)
    state_grid::RectangleGrid{5} = RectangleGrid(ego_grid, velocity_grid, acceleration_grid, human_grid, velocity_grid)
    goal_reward::Float64 = 1.0
    collision_cost::Float64 = 0.0
    discount_factor::Float64 = 1.0
end

const AbsentPos = -1.0

POMDPs.discount(mdp::MergingMDP) = mdp.discount_factor
POMDPs.actions(mdp::MergingMDP) = 1:7
POMDPs.n_actions(mdp::MergingMDP) = 7
POMDPs.actionindex(mdp::MergingMDP, a::Int64) = a

function POMDPs.states(mdp::MergingMDP)
    prod = Iterators.product(mdp.ego_grid,
                             mdp.velocity_grid,
                             mdp.acceleration_grid,
                             vcat(mdp.human_grid, AbsentPos),
                             mdp.velocity_grid)
    vecs = imap(x->Float64[x...], prod)
end
POMDPs.n_states(mdp::MergingMDP) = reduce(*, length.([mdp.ego_grid, mdp.velocity_grid, mdp.acceleration_grid, vcat(mdp.human_grid, AbsentPos), mdp.velocity_grid]))

function POMDPs.stateindex(mdp::MergingMDP, s::MergingState)
    se, ve, ae, so, vo = s
    sei = findfirst(x->isapprox(x, se), mdp.ego_grid)
    vei = findfirst(x->isapprox(x, ve), mdp.velocity_grid)
    aei = findfirst(x->isapprox(x, ae), mdp.acceleration_grid)
    soi = findfirst(x->isapprox(x, so), mdp.human_grid)
    if so == AbsentPos
        soi = length(mdp.human_grid) + 1
    end
    voi = findfirst(x->isapprox(x, vo), mdp.velocity_grid)
    return LinearIndices((length(mdp.ego_grid), length(mdp.velocity_grid), length(mdp.acceleration_grid),
                          length(mdp.human_grid)+1, length(mdp.velocity_grid)))[sei, vei, aei, soi, voi]
end

function POMDPs.reward(mdp::MergingMDP, s::MergingState, a::Int64, sp::MergingState)
    if first(sp) >= last(mdp.ego_grid)
        return mdp.goal_reward
    end
    if collisioncheck(mdp, sp)
        return mdp.collision_cost
    end
    return 0.0
end

POMDPs.isterminal(mdp::MergingMDP, s::MergingState) = first(s) >= last(mdp.ego_grid) || collisioncheck(mdp, s)

function POMDPs.transition(mdp::MergingMDP, s::MergingState, a::Int64)
    if first(s) >= last(mdp.ego_grid)
        return SparseCat([s], [1.0])
    end
    se, ve, ae = ego_transition(mdp, s, a)
    sitps = MergingState[]
    witps = Float64[]
    if s[4] == AbsentPos
        so = AbsentPos
        vo = first(mdp.velocity_grid)
        sitp, witp = interpolants(mdp.state_grid, SVector(se, ve, ae, so, vo))
        sitp = ind2x.(Ref(mdp.state_grid), sitp)
        for s in sitp
            s[4] = AbsentPos
        end
        sitps = vcat(sitps, sitp)
        witps = vcat(witps, witp)
    else
        for ah in mdp.human_actions
            so, vo = other_transition(mdp, s, ah)
            sitp, witp = interpolants(mdp.state_grid, SVector(se, ve, ae, so, vo))
            sitp = ind2x.(Ref(mdp.state_grid), sitp)
            if so == AbsentPos
                for s in sitp
                    s[4] = AbsentPos
                end
            end
            sitps = vcat(sitps, sitp)
            witps = vcat(witps, witp)
        end
    end
    sitps, witps = remove_doublons(sitps, witps)
    normalize!(witps, 1)
    return SparseCat(sitps, witps)
end

function remove_doublons(s::Vector{S}, w::Vector{Float64}) where S
    ns = S[]
    ws = Float64[]
    for (i, s) in enumerate(s)
        si = findfirst(isequal(s), ns)
        if si != nothing
            ws[si] += w[i]
        else
            push!(ns, s)
            push!(ws, w[i])
        end
    end
    return ns, ws
end

function ego_transition(mdp::MergingMDP, s::MergingState, a::Int64)
    se, v, acc, so, vo = s
    acc_ = acc
    if a == 1 # full brake
        acc_ = mdp.min_acc
    elseif a == 7 # release throttle
        acc_ = 0.0
    else
        acc_ += mdp.jerk_levels[a - 1]
    end
    v_ = v + acc_ * mdp.dt
    s_ = se + v_*mdp.dt
    return s_, v_, acc_
end

function other_transition(mdp::MergingMDP, s::MergingState, ah::Float64)
    se, v, acc, so, vo = s
    vo_ = vo + ah* mdp.dt
    so_ = so + vo_*mdp.dt
    if so_ >= last(mdp.human_grid)
        so_ = AbsentPos
    end
    return so_, vo_
end

## helpers
function collisioncheck(mdp::MergingMDP, s::MergingState)
    ego_posG = get_posG(Frenet(mdp.ego_lane, s[1]), mdp.env.roadway)
    other_posG = get_posG(Frenet(mdp.other_lane, s[4]), mdp.env.roadway)
    ego_p = polygon(ego_posG, mdp.cardef)
    other_p = polygon(other_posG, mdp.cardef)
    return overlap(ego_p, other_p)
end

function getegostate(mdp::MergingMDP, s::MergingState)
    ego_posF = Frenet(mdp.ego_lane, s[1])
    return VehicleState(ego_posF, mdp.env.roadway, s[2])
end

function getotherstate(mdp::MergingMDP, s::MergingState)
    other_posF = Frenet(mdp.other_lane, s[4])
    return VehicleState(other_posF, mdp.env.roadway, s[5])
end

## rendering

@with_kw struct MergingViz
    mdp::MergingMDP = MergingMDP()
    s::MergingState =  [mdp.ego_grid[1], mdp.velocity_grid[1], mdp.acceleration_grid[1], mdp.human_grid[1], mdp.velocity_grid[1]]
    a::Union{Nothing, Int64} = nothing
    action_values::Union{Nothing, Vector{Float64}} = nothing
end

function AutoViz.render!(rm::RenderModel, v::MergingViz)
    render!(rm , mdp.env.roadway)
    se, ve, acc, so, vo = v.s
    text_ego = @sprintf("v: %2.2f", ve)
    act = v.a
    if act != nothing
        as = "a = "
        if act == 1
            as *= "full brake"
        elseif act == 7
            as *= "release"
        else
            as *=  @sprintf("%1.2f", mdp.jerk_levels[act - 1])
        end
        text_ego *= ", " * as
    end
    av = v.action_values
    if av != nothing
        text_ego *= "Probas " * reduce(*, @sprintf("%1.2f, ", pa) for pa in av)
    end
    ego_posG = get_posG(Frenet(mdp.ego_lane, se), mdp.env.roadway)
    ego_acar = ArrowCar(ego_posG.x, ego_posG.y, ego_posG.θ,
                        color=COLOR_CAR_EGO, text=text_ego,
                        length=mdp.cardef.length, width=mdp.cardef.width)
    other_posG = get_posG(Frenet(mdp.other_lane, so), mdp.env.roadway)
    if so == AbsentPos
        other_color = RGB(0.5,0.5,0.5)
    else
        other_color = COLOR_CAR_OTHER
    end
    other_acar = ArrowCar(other_posG.x, other_posG.y, other_posG.θ,
                          color=other_color, text="v: $vo",
                          length=mdp.cardef.length, width=mdp.cardef.width)
    render!(rm, ego_acar)
    render!(rm, other_acar)
end

function AutoViz.getcenter(s::MergingState)
    se, v, acc, so, vo = s
    ego_posG = get_posG(Frenet(mdp.ego_lane, se), mdp.env.roadway)
    other_posG = get_posG(Frenet(mdp.other_lane, so), mdp.env.roadway)
    return 0.5*SVector(ego_posG.x + other_posG.x, ego_posG.y + other_posG.y)
end


text:
"""
    name_aliases()

Returns a `Dict` containing:
 - a key (`String`), which we consider to be our core variable name
 - values (`Tuple` of `String`s), which we consider to be aliases of our core variable name
"""
function name_aliases()
    dict = Dict(
        "zc" => ("z_half",),
        "zf" => ("z",),
        "α0_c" => ("alpha0_half",),
        "α0_f" => ("alpha0",),
        "p0_c" => ("p0_half",),
        "p0_f" => ("p0",),
        "ρ0_c" => ("rho0_half",),
        "ρ0_f" => ("rho0",),
        "updraft_area" => ("updraft_fraction",),
        "updraft_thetal" => ("updraft_thetali",),
        "thetal_mean" => ("thetali_mean", "theta_mean"),
        "total_flux_h" => ("resolved_z_flux_thetali", "resolved_z_flux_theta"),
        "total_flux_qt" => ("resolved_z_flux_qt", "qt_flux_z"),
        "u_mean" => ("u_translational_mean",),
        "v_mean" => ("v_translational_mean",),
        "tke_mean" => ("tke_nd_mean",),
        "total_flux_s" => ("s_flux_z",),
    )
    return dict
end

"""
    get_nc_data(ds::NCDatasets.Dataset, var::String)

Returns the data for variable `var`, trying first its aliases
defined in `name_aliases`, in the `ds::NCDatasets.Dataset`.
"""
function get_nc_data(ds, var::String)
    dict = name_aliases()
    key_options = haskey(dict, var) ? (dict[var]..., var) : (var,)

    for key in key_options
        if haskey(ds, key)
            return ds[key]
        else
            for group_option in ["profiles", "reference", "timeseries"]
                haskey(ds.group, group_option) || continue
                if haskey(ds.group[group_option], key)
                    return ds.group[group_option][key]
                end
            end
        end
    end
    return nothing
end


text:
"""
TagData(data1, data2; batch=128, ftype=Float32, dense=false) creates a
data generator that can be used with a Tagger model.  The source data1
and target data2 should be sequence generators, i.e. next(data1)
should deliver a vector of Ints that represent the next sequence.
This division of labor allows different file formats to be supported.
maxtoken(data1) should give the largest integer produced by data1.

The following transformations are performed by an TagData generator:

* sequences are minibatched according to the batch argument.
* sequences in a minibatch padded to all be the same length.
* the end-of-sequence is indicated using a "nothing" token.

"""
type TagData; x; y; mask;
    function TagData(xgen, ygen; batchsize=128, ftype=Float32, dense=false, o...)
        new(TagBatch(xgen, batchsize, ftype, dense),
            TagBatch(ygen, batchsize, ftype, dense),
            zeros(Cuchar, batchsize))
    end
end

type TagBatch; sgen; state; sent; word; done;
    function TagBatch(sgen, batchsize, ftype, dense)
        word = (dense ? zeros : sponehot)(ftype, maxtoken(sgen), batchsize)
        sent = Array(Any, batchsize)
        new(sgen, nothing, sent, word, false)
    end
end

import Base: start, done, next

# the TagData state is just nword, which is the number of words
# completed in the current batch.

function start(d::TagData)
    d.x.done = d.y.done = false
    d.x.state = start(d.x.sgen)
    d.y.state = start(d.y.sgen)
    nextbatch(d.x)
    nextbatch(d.y)
    return 0
end

# We stop when there is not enough data to fill a batch

done(d::TagData, state)=(d.x.done && d.y.done)

# TagData.next returns the next token

function next(d::TagData, nword)
    nword += 1
    maxlen = maximum(map(length, d.x.sent))
    if nword > maxlen
        nextbatch(d.x)
        nextbatch(d.y)
        return (nothing, 0)
    end
    for s=1:length(d.x.sent)
        n=length(d.x.sent[s])
        xword = (nword <= n ? d.x.sent[s][nword] : 0)
        yword = (nword <= n ? d.y.sent[s][nword] : 0)
        # @show s,n,xword,yword
        setrow!(d.x.word, xword, s)
        setrow!(d.y.word, yword, s)
        d.mask[s] = (xword == 0 ? 0 : 1)
    end
    return ((d.x.word, d.y.word, d.mask), nword)
end

function nextbatch(b::TagBatch)
    for i=1:length(b.sent)
        done(b.sgen, b.state) && (b.done=true; return)
        (b.sent[i], b.state) = next(b.sgen, b.state)
    end
end

maxtoken(s::TagData,i)=maxtoken(i==1 ? s.x.sgen : i==2 ? s.y.sgen : error())
maxtoken(x::Vector{Vector{Int}}) = maximum(map(maximum,x))

sponehot(ftype, m, n)=sparse(ones(Int32,n), Int32[1:n;], one(ftype), m, n)

# # TODO: these assume one hot columns, make them more general.
# setrow!(x::SparseMatrixCSC,i,j)=(i>0 ? (x.rowval[j] = i; x.nzval[j] = 1) : (x.rowval[j]=1; x.nzval[j]=0))
# setrow!(x::Array,i,j)=(x[:,j]=0; i>0 && (x[i,j]=1))

# function readvocab(file) # TODO: test with cmd e.g. `zcat foo.gz`
#     d = Dict{Any,Int}()
#     open(file) do f
#         for l in eachline(f)
#             for w in split(l)
#                 get!(d, w, 1+length(d))
#             end
#         end
#     end
#     return d
# end


text:
module TestIteration
    using Base.Test
    using DataFrames

    dv = @data([1, 2, NA])
    dm = DataArray([1 2; 3 4])
    dt = DataArray(zeros(2, 2, 2))

    df = DataFrame(A = 1:2, B = 2:3)

    for el in dv
        @test ndims(el) == 0
    end

    for el in dm
        @test ndims(el) == 0
    end

    for el in dt
        @test ndims(el) == 0
    end

    for row in eachrow(df)
        @test isa(row, DataFrameRow)
        @test row[:B]-row[:A] == 1
    end

    for col in eachcol(df)
        @test isa(col, (Symbol, AbstractDataVector))
    end

    @test isequal(map(x -> minimum(array(x)), eachrow(df)), {1,2})
    @test isequal(map(minimum, eachcol(df)), DataFrame(A = [1], B = [2]))

    row = DataFrameRow(df, 1)

    row[:A] = 100
    @test df[1, :A] == 100

    row[1] = 101
    @test df[1, :A] == 101

    df = DataFrame(A = 1:4, B = ["M", "F", "F", "M"])

    s1 = sub(df, 1:3)
    s1[2,:A] = 4
    @test df[2, :A] == 4
    @test sub(s1, 1:2) == sub(df, 1:2)

    s2 = sub(df, 1:2:3)
    s2[2, :B] = "M"
    @test df[3, :B] == "M"
    @test sub(s2, 1:1:2) == sub(df, [1,3])

    # @test_fail for x in df; end # Raises an error
end


text:
using WallTimeProgress
using Test

import WallTimeProgress: next_period, _with_underscores, increment!, reset!

@testset "utilities" begin
    @test _with_underscores(199) == "199"
    @test _with_underscores(0) == "0"
    @test _with_underscores(1099) == "1_099"
    @test _with_underscores(1700000) == "1_700_000"
end

@testset "next period" begin
    @test next_period(1, 10) == 10
    @test next_period(9, 10) == 10
    @test next_period(0, 10) == 10
    @test next_period(10, 10) == 20
end

@testset "WallTimeTracker" begin
    io = IOBuffer()
    wtt = WallTimeTracker(10; output = io)
    increment!(wtt)
    increment!(wtt, 5)
    increment!(wtt, 4)
    increment!(wtt, 9)
    increment!(wtt, 9)
    output = split(String(take!(io)), '\n')
    println(output)
    @test occursin(r"processed 10 records, [-e.0-9]+ s/record", output[1])
    @test occursin(r"processed 28 records, [-e.0-9]+ s/record", output[2])
    @test count(wtt) == 28
    @test wtt.next_report == 30
end


text:
function angle(vec1::Union{Array{F,1},Array{Int,1}},vec2::Union{Array{F,1},Array{Int,1}}) where {F <: AbstractFloat}
    acos(dot(vec1, vec2) / (norm(vec1) * norm(vec2))) * 180 / pi
    #TODO: fix DomainError with 1.0000000000000002 problem
end

function vecs_to_triangle(veci::Union{Array{F,1},BitArray{1}},
                          vecj::Union{Array{F,1},BitArray{1}},
                          veck::Union{Array{F,1},BitArray{1}}) where {F <: AbstractFloat}
    seg_ij = vecj - veci
    seg_ik = veck - veci
    seg_jk = veck - vecj

    angle1 = angle(seg_ij, seg_ik)
    angle2 = angle(-seg_ij, seg_jk)
    angle3 = 180.0 - angle1 - angle2
    return (angle1, angle2, angle3)
end

function brute_force_triangles(pts::Array{F,2}) where {F <: AbstractFloat}

    n,d = size(pts)
    # n choose 3 triangles
    triangles = Array{Tuple{Tuple{Int,Int,Int},Tuple{F,F,F}},1}(undef,binomial(n,3))

    index = 1
    for i =1:n
        for j =1:i-1
            for k = 1:j-1
                x = ((k,j,i),vecs_to_triangle(pts[k,:],pts[j,:],pts[i,:]))
                triangles[index] = x
                index += 1
            end
        end
    end

    return triangles
end

function brute_force_triangles(pts::Array{BitArray{1},1}) where {F <: AbstractFloat}

    n = length(pts)
    # n choose 3 triangles
    #TODO: fix F type error
    triangles = Array{Tuple{Tuple{Int,Int,Int},Tuple{Float64,Float64,Float64}},1}(undef,binomial(n,3))

    index = 1
    for i =1:n
        for j =1:i-1
            for k = 1:j-1
                x = ((k,j,i),vecs_to_triangle(pts[k],pts[j],pts[i]))
                triangles[index] = x
                index += 1
            end
        end
    end

    return triangles

end


text:
module FighterJets

export FighterJet

using Base: Bool, Float64
#FigherJet is avery fast vehicle with powerful weapons. It
mutable struct FighterJet

    "power status: true - on, false - off"
    power::Bool

    "current direction in radians"
    direction::Float64

    "current position coordinate (x, y)"
    position::Tuple{Float64, Float64}
end

# Import generic functions from
include("../../Vehicles/src/Vehicles.jl")
using .Vehicles: power_on!, power_off!, turn!, move!, position!, Vehicles

# Implementation of Vehicle interface
function power_on!(fj::FighterJet)
    fj.power =true
    println("Powered on: ", fj)
    nothing
end

function power_off!(fj::FighterJet)
    fj.power = false
    println("Powered off: ", fj)
    nothing
end

function turn!(fj::FighterJet, direction)
    fj.direction = direction
    println("Changed direction to ", direction, ": ", fj)
    nothing
end

function move!(fj::FighterJet, distance)
    x, y = fj.position
    dx =  round(distance * cos(fj.direction), digits = 2)
    dy =  round(distance * sin(fj.direction), digits = 2)
    fj.position = (x + dx, y + dy)
    println("Moved (", dx, ",", dy, "): ", fj)
    nothing
end

function position(fj::FighterJet)
    fj.position
end

end # module


text:
#
# Sensitivity analysis related types and functions
#

struct Sensitivity <: Wrapper
    pyo::PyObject

    function Sensitivity(pyo::PyObject)
        new(pyo)
    end
end

struct SAResult  <: Wrapper
    pyo::PyObject

    function SAResult(pyo::PyObject)
        new(pyo)
    end
end

function sa(m::Model, response; policy = Dict(), method = "sobol", nsamples = 1000, kwargs...)
    my_sa = py"my_sa"
    pyo = my_sa(m.pyo, response; policy = policy, method = method, nsamples = nsamples, kwargs...)
    return SAResult(pyo."sa_result")   # access without conversion
end

function plot(sar::SAResult; kwargs...)
    return rhodium.SAResult.plot(sar.pyo; kwargs...)
end

function plot_sobol(sar::SAResult;
                    radSc = 2.0, scaling = 1, widthSc = 0.5, STthick = 1, varNameMult = 1.3,
                    colors = nothing, groups = nothing, gpNameMult = 1.5, threshold = "conf")

    fig = rhodium.SAResult.plot_sobol(sar.pyo; radSc = radSc, scaling = scaling, widthSc = widthSc,
                               STthick = STthick, varNameMult = varNameMult,
                               colors = colors, groups = groups, gpNameMult = gpNameMult,
                               threshold = threshold)
    return fig
end

function oat(m::Model, response; policy = Dict(), nsamples::Int = 100, kwargs...)
    fig = rhodium.oat(m.pyo, response; policy = policy, nsamples = nsamples, kwargs...)
    return fig
end


text:
using Distributions, StatsBase
export Hollywood

struct Hollywood
    α::Real
    θ::Real
    ν::DiscreteUnivariateDistribution
    K::Real
    function Hollywood(α::Real, θ::Real, ν::DiscreteUnivariateDistribution, K::Real)
        # @show θ, -K*α
        @assert ((0 < α < 1.0) & (θ > -α)) | ((α < 0) & (θ == (- K * α))) "Check parameters satisfy required constraints"
        new(α, θ, ν, K)
    end
end

Hollywood(α::Real, θ::Real, ν::DiscreteUnivariateDistribution) = Hollywood(α, θ, ν, Inf)
Hollywood(α::Real, ν::DiscreteUnivariateDistribution, K::Int) = Hollywood(α, (-K)*α, ν, K)

function StatsBase.sample!(
    out::InteractionSequence,
    model::Hollywood
)
    # Init
    t = 0
    V = 0
    counts = Int[]
    prob_new(t, V, θ, α) = (θ + α * V) / (t - 1 + θ)
    for i in 1:length(out)
        m = rand(model.ν)
        I_tmp = zeros(Int, m)
        for j in 1:m
            t += 1
            if t==1
                I_tmp[j] = 1
                V += 1
                push!(counts, 1)
            elseif rand() < prob_new(t, V, model.θ, model.α)
                I_tmp[j] = V+1
                V+=1
                push!(counts, 1)
            else
                μ = cumsum(counts) / sum(counts)
                pushfirst!(μ, 0.0)
                # @show μ
                val, p = rand_multivariate_bernoulli(μ)
                I_tmp[j] = val
                counts[val] += 1
            end
            out[i] = I_tmp
        end
    end
end

function StatsBase.sample(
    model::Hollywood,
    n::Int
    )
    out = fill(Int[], n)
    sample!(out, model)
    return out
end


text:
using .SampleChainsDynamicHMC
using Random

export sample

using ..Soss

function sample(rng::AbstractRNG,
    ::Type{DynamicHMCChain},
    m::ConditionalModel,
    nsamples::Int=1000,
    nchains::Int=4)

    ℓ(x) = logdensity(m, x)
    tr = xform(m)

    chains = initialize!(rng, nchains, DynamicHMCChain, ℓ, tr)
    drawsamples!(chains, nsamples - 1)
    return chains
end

function sample(
    ::Type{DynamicHMCChain},
    m::ConditionalModel,
    nsamples::Int=1000,
    nchains::Int=4)

    sample(Random.GLOBAL_RNG, DynamicHMCChain, m, nsamples, nchains)
end


text:
struct AnnulusView{T,N,M<:AbstractArray{T,N},IT} <: AbstractArray{T,N}
    parent::M
    rmin::Float64
    rmax::Float64
    indices::IT
    fill::T
end

"""
    AnnulusView(cube::AbstractArray{T,3};
                inner=0, outer=last(size(parent))/2 + 0.5,
                fill=0)

Cut out an annulus with inner radius `inner` and outer radius `outer`. Values that fall outside of this region will be replaced with `fill`. This does not copy any data, it is merely a view into the data.
"""
function AnnulusView(parent::AbstractArray{T,3};
                     inner=0,
                     outer=(size(parent, 3) + 1) / 2,
                     fill=zero(T)) where T
    # check inputs
    0 ≤ inner < outer || error("Invalid annulus region [$inner, $outer]")
    time_axis = axes(parent, 1)
    space_indices = CartesianIndices((axes(parent, 2), axes(parent, 3)))
    space_axis = filter(idx -> inside_annulus(inner, outer, center(parent)[2:3], idx), space_indices)
    return AnnulusView(parent, Float64(inner), Float64(outer), (time_axis, space_axis), T(fill))
end

Base.parent(view::AnnulusView) = view.parent
Base.size(view::AnnulusView) = size(parent(view))
Base.copy(view::AnnulusView) = AnnulusView(copy(parent(view)), view.rmin, view.rmax, view.indices, view.fill)

inside_annulus(view::AnnulusView, args...) = inside_annulus(view.rmin, view.rmax, center(view)[2:3], args...)

@propagate_inbounds function Base.getindex(view::AnnulusView{T,N}, idx::Vararg{<:Integer,N}) where {T,N}
    @boundscheck checkbounds(parent(view), idx...)
    ifelse(inside_annulus(view, idx...), convert(T, parent(view)[idx...]), view.fill)
end

@propagate_inbounds function Base.setindex!(view::AnnulusView{T,N}, val, idx::Vararg{<:Integer,N}) where {T,N}
    @boundscheck checkbounds(parent(view), idx...)
    if inside_annulus(view, idx...)
        parent(view)[idx...] = val
    else
        view.fill
    end
end

"""
    (::AnnulusView)(asview=false)

Return the pixels that fall within the annulus as a matrix. This matrix is equivalent to unrolling each frame and then spatially filtering the pixels outside the annulus. If `asview` is true, the returned values will be a view of the parent array instead of a copy.

# Examples
```jldoctest
julia> ann = AnnulusView(ones(10, 101, 101); inner=5, outer=20);

julia> X = ann();

julia> size(X)
(10, 1188)
```
"""
function (view::AnnulusView)(asview=false)
    if asview
        @view parent(view)[view.indices...]
    else
        parent(view)[view.indices...]
    end
end

"""
    copyto!(::AnnulusView, mat::AbstractMatrix)

Copy the pixels from `mat` into the pixels in the annulus. `mat` should have the same size as the matrix output from [`AnnulusView`](@ref)

# Examples
```jldoctest
julia> ann = AnnulusView(ones(10, 101, 101); inner=5, outer=20);

julia> X = ann();

julia> new_ann = copyto!(ann, -X);

julia> new_ann() ≈ -X
true
```
"""
function Base.copyto!(view::AnnulusView, mat::AbstractMatrix)
    inverse!(view, view.parent, mat)
    return view
end

"""
    inverse!(::AnnulusView, out, mat)

In-place version of [`inverse`](@ref) that fills `out` in-place.
"""
function inverse!(view::AnnulusView, out, mat)
    out[view.indices...] = mat
    return out
end

"""
    inverse(::AnnulusView, mat::AbstractMatrix)

Generate a cube similar to the view with the pixels from `mat`. `mat` should have the same size as the matrix output from [`AnnulusView`](@ref)

# Examples
```jldoctest
julia> ann = AnnulusView(ones(10, 101, 101); inner=5, outer=20);

julia> X = ann();

julia> out = inverse(ann, -X);

julia> out ≈ -ann
true
```
"""
function inverse(view::AnnulusView, mat)
    out = fill!(similar(parent(view)), view.fill)
    inverse!(view, out, mat)
end


text:
struct BGGReview # review scraped from game page
    id::Int
    name::String
    username::String
    rating::Float32
    comment::String
end

function Base.show(io::IO, r::BGGReview)
    return print(io, "$(r.rating)/10 for $(r.name) from $(r.username)")
end

struct BGGGameInfo
    id::Int
    name::String
    yearpublished::Int
    minplayers::Int
    maxplayers::Int
    playingtime::Int
    minplaytime::Int
    maxplaytime::Int
    minage::Int
    usersrated::Int
    average::Float32
    bayesaverage::Float32
    stddev::Float32
    median::Float32
    owned::Int
    trading::Int
    wanting::Int
    wishing::Int
    numcomments::Int
    numweights::Int
    averageweight::Float32
    mechanics::Vector{String}
    families::Vector{String}
    suggested_numplayers::Dict{String,Tuple{Int64,Int64,Int64}}
end

"""
    gamereviews(id::Integer; kwargs...)

Scrape all reviews for a game

# Keyword arguments
  - `pagesize::Integer`:
    Number of reviews per page request. Defaults to API maximum 100.
  - `waittime::Real`:
    Wait time between page requests in seconds. Defaults to `2.0f0`.
"""
function gamereviews(id::Integer; pagesize::Integer=100, waittime=2.0f0)
    doc = get_xml("$XMLAPI2/thing?id=$(id)&stats=1")
    name = findfirst("/items/item/name[@type='primary']", doc)["value"]
    usersrated = parse(
        Int, findfirst("/items/item/statistics/ratings/usersrated", doc)["value"]
    )
    maxpage = cld(usersrated, pagesize) # using `pagesize` reviews per page
    maxpage == 0 && (maxpage = 1)
    # The actual number of reviews is usually slightly higher than `usersrated`.
    # Count number of reviews on last page to get the real count.
    doc = get_xml(
        "$XMLAPI2/thing?id=$(id)&ratingcomments=1&page=$(maxpage)&pagesize=$(pagesize)"
    )
    review_count =
        (maxpage - 1) * pagesize + length(findall("/items/item/comments/comment", doc))

    if review_count == 0
        return Vector{BGGReview}()
    end

    reviews = Vector{BGGReview}(undef, review_count) # pre-allocate
    @showprogress 1 "Scraping $review_count reviews for $name..." for page in 1:maxpage
        doc = get_xml(
            "$XMLAPI2/thing?id=$(id)&ratingcomments=1&page=$(page)&pagesize=$(pagesize)"
        )
        pagereviews = findall("/items/item/comments/comment", doc)
        for (i, r) in enumerate(pagereviews)
            reviews[(page - 1) * pagesize + i] = BGGReview(
                id, name, r["username"], parse(Float32, r["rating"]), r["value"]
            )
        end
        sleep(waittime)
    end
    return reviews
end

"""
    gameinfo(id::Integer)
    gameinfo(ids::AbstractArray{<:Integer})

Get game information and summary of reviews.
"""
function gameinfo(id::Integer; waittime=0.0f0)
    doc = get_xml("$XMLAPI2/thing?id=$(id)&stats=1")
    game = findfirst("/items/item", doc)

    # Primary name of game
    name = findfirst("name[@type='primary']", game)["value"]

    # Parse mechanics and families
    boardgamemechanic = [
        n["value"] for n in findall("link[@type='boardgamemechanic']", game)
    ]
    boardgamefamily = [n["value"] for n in findall("link[@type='boardgamefamily']", game)]

    # Find first index after list of names
    yearpublished = parse(Int, findfirst("yearpublished", game)["value"])
    minplayers = parse(Int, findfirst("minplayers", game)["value"])
    maxplayers = parse(Int, findfirst("maxplayers", game)["value"])

    # Vote counts for player count recommendations
    suggested_numplayers = Dict{String,Tuple{Int64,Int64,Int64}}()
    for res in findall("poll[@name='suggested_numplayers']/results", game)
        numplayers = res["numplayers"]
        suggested_numplayers[numplayers] = tuple(
            parse.(Int, r["numvotes"] for r in elements(res))...
        )
    end

    playingtime = parse(Int, findfirst("playingtime", game)["value"])
    minplaytime = parse(Int, findfirst("minplaytime", game)["value"])
    maxplaytime = parse(Int, findfirst("maxplaytime", game)["value"])
    minage = parse(Int, findfirst("minage", game)["value"])

    # Parse ratings summary
    ratings = findfirst("statistics/ratings", game)
    usersrated = parse(Int, findfirst("usersrated", ratings)["value"])
    average = parse(Float32, findfirst("average", ratings)["value"])
    bayesaverage = parse(Float32, findfirst("bayesaverage", ratings)["value"])
    stddev = parse(Float32, findfirst("stddev", ratings)["value"])
    median = parse(Float32, findfirst("median", ratings)["value"])
    owned = parse(Int, findfirst("owned", ratings)["value"])
    trading = parse(Int, findfirst("trading", ratings)["value"])
    wanting = parse(Int, findfirst("wanting", ratings)["value"])
    wishing = parse(Int, findfirst("wishing", ratings)["value"])
    numcomments = parse(Int, findfirst("numcomments", ratings)["value"])
    numweights = parse(Int, findfirst("numweights", ratings)["value"])
    averageweight = parse(Float32, findfirst("averageweight", ratings)["value"])

    sleep(waittime)
    return BGGGameInfo(
        id,
        name,
        yearpublished,
        minplayers,
        maxplayers,
        playingtime,
        minplaytime,
        maxplaytime,
        minage,
        usersrated,
        average,
        bayesaverage,
        stddev,
        median,
        owned,
        trading,
        wanting,
        wishing,
        numcomments,
        numweights,
        averageweight,
        boardgamemechanic,
        boardgamefamily,
        suggested_numplayers,
    )
end

function gameinfo(ids::AbstractArray{<:Integer}; waittime=2.0f0) # add waittime by default
    gamecount = length(ids)
    games = Vector{BGGGameInfo}(undef, gamecount) # pre-allocate
    @showprogress 1 "Scraping $gamecount games..." for (i, n) in enumerate(ids)
        games[i] = gameinfo(n; waittime=waittime)
    end
    return games
end


text:
# Autogenerated wrapper script for argp_standalone_jll for armv6l-linux-musleabihf
export argp_h, libargp

JLLWrappers.@generate_wrapper_header("argp_standalone")
JLLWrappers.@declare_file_product(argp_h)
JLLWrappers.@declare_file_product(libargp)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_file_product(
        argp_h,
        "include/argp.h",
    )

    JLLWrappers.@init_file_product(
        libargp,
        "lib/libargp.a",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
mutable struct BinarySearchNode{T} <: BinaryNode{T}
    key
    data::T
    leftchild::BinarySearchNode{T}
    rightchild::BinarySearchNode{T}
    function BinarySearchNode{T}() where T
        n = new()
        n.leftchild = n.rightchild = n
        n
    end
    function BinarySearchNode{T}(key, data) where T
        n = new()
        n.key = key
        n.data = data
        n.leftchild = n.rightchild = n
        n
    end
    function BinarySearchNode{T}(key, data, leftchild, rightchild) where T
        new(key, data, leftchild, rightchild)
    end
end

function BinarySearchNode(key, data, lc::BinarySearchNode{T}, rc::BinarySearchNode{T}) where T
    BinarySearchNode{T}(key, data, lc, rc)
end
BinarySearchNode(key, data::T) where T = BinarySearchNode{T}(key, data)

getvalue(n::BinarySearchNode) = n.data
setvalue!(n::BinarySearchNode, item) = (n.data = item)

getkey(n::BinarySearchNode) = n.key
setkey!(n::BinarySearchNode, key) = (n.key = key)

leftchild(n::BinarySearchNode) = n.leftchild
setleftchild!(n::BinarySearchNode, node) = (n.leftchild = node)

rightchild(n::BinarySearchNode) = n.rightchild
setrightchild!(n::BinarySearchNode, node) = (n.rightchild = node)

isleaf(n::BinarySearchNode) = (n === n.leftchild === n.rightchild)