text:
###########################################################################
### ProposalDensities wrapping distributions from the Distributions package
###########################################################################
type SingleSymmetricDistributionWrapper{D<:Distributions.Distribution} <: SymmetricDensity
    distribution::D
end

type SingleASymmetricDistributionWrapper{D<:Distributions.Distribution} <: ASymmetricDensity
    distribution::D
end

type CompoundSymmetricDistributionWrapper{D<:Distributions.UnivariateDistribution} <: SymmetricDensity
    distributions::Vector{D}
end

type CompoundASymmetricDistributionWrapper{D<:Distributions.UnivariateDistribution} <: ASymmetricDensity
    distributions::Vector{D}
end

typealias SingleDistributionWrapper Union{SingleSymmetricDistributionWrapper, SingleASymmetricDistributionWrapper}
typealias CompoundDistributionWrapper Union{CompoundSymmetricDistributionWrapper, CompoundASymmetricDistributionWrapper}
typealias DistributionWrapper Union{SingleDistributionWrapper,CompoundDistributionWrapper}

densityname(d::DistributionWrapper) = "DistributionWrapper"

typealias NormalDensity SingleSymmetricDistributionWrapper{Distributions.MvNormal}
typealias LogNormalDensity SingleASymmetricDistributionWrapper{Distributions.MvLogNormal}
typealias UniformDensity CompoundSymmetricDistributionWrapper{Distributions.Uniform}
typealias LaplaceDensity CompoundSymmetricDistributionWrapper{Distributions.Laplace}
typealias TriangularDensity CompoundSymmetricDistributionWrapper{Distributions.SymTriangularDist}
typealias BactrianDensity CompoundSymmetricDistributionWrapper{Bactrian}

_density(::Type{Val{:normal}},x::AbstractVector,s::AbstractArray) = NormalDensity(_distribution(Val{:normal},x,s))
_density(::Type{Val{:normal}},s::AbstractArray) = NormalDensity(_distribution(Val{:normal},s))
_density(::Type{Val{:lognormal}},x::AbstractVector,s::AbstractArray) = LogNormalDensity(_distribution(Val{:lognormal},x,s))
_density(::Type{Val{:uniform}},x::AbstractVector,s::AbstractVector) = UniformDensity(_distributions(Val{:uniform},x,s))
_density(::Type{Val{:laplace}},x::AbstractVector,s::AbstractVector) = LaplaceDensity(_distributions(Val{:laplace},x,s))
_density(::Type{Val{:triangular}},x::AbstractVector,s::AbstractVector) = TriangularDensity(_distributions(Val{:triangular},x,s))
_density(::Type{Val{:bactrian}},x::AbstractVector,s::AbstractVector,subtype::Symbol,mixing::AbstractFloat) = BactrianDensity(_distributions(Val{:bactrian},x,s,subtype,mixing))

###Propose points by calling rand! on the corresponding distribution
@inline propose!(d::SingleDistributionWrapper,v::AbstractArray) = (rand!(d.distribution,v) ; v)

@inline function propose!(d::CompoundDistributionWrapper,v::AbstractArray)
    nrows = size(v,1)
    for j=1:size(v,2)
        @simd for i=1:nrows
            @inbounds v[i,j] = rand(d.distributions[i])
        end
    end
    v
end

###Calculate the logprobability by calling logpdf! on the corresponding distribution
@inline logprobability!(r::AbstractVector,d::SingleDistributionWrapper,v::AbstractVector) = (@inbounds r[1] = logpdf(d.distribution,v) ; r)
@inline logprobability!(r::AbstractVector,d::SingleDistributionWrapper,v::AbstractMatrix) = (logpdf!(r,d.distribution,v) ; r)
@inline logprobability(d::SingleDistributionWrapper,v::AbstractArray) = logprobability!(zeros(eltype(v),size(v,2)),d,v)

@inline function logprobability!(r::AbstractVector,d::CompoundDistributionWrapper,v::AbstractArray)
    nrows = size(v,1)
    for j=1:size(v,2)
        @inbounds r[j] = zero(eltype(v))
        @simd for i=1:nrows
            @inbounds r[j] += logpdf(d.distributions[i],v[i,j])
        end
    end
    r
end

@inline logprobability(d::CompoundDistributionWrapper,v::AbstractArray) = logprobability!(zeros(eltype(v),size(v,2)),d,v)

#update functions for single distributions
condition!(d::SingleDistributionWrapper,x::AbstractVector) = (d.distribution = recenter(d.distribution,x) ; d)
scale!(d::SingleDistributionWrapper,s::AbstractFloat) = (d.distribution = rescale(d.distribution,s) ; d)
update!(d::SingleDistributionWrapper,x::AbstractVector,s::AbstractArray) = (d.distribution = update(d.distribution,x,s) ; d)
update!(d::NormalDensity,s::AbstractArray) = (d.distribution = update(d.distribution,s) ; d)

#update functions for compound distributions
condition!(d::CompoundDistributionWrapper,x::AbstractVector) = (map!(recenter,d.distributions,d.distributions,x) ; d)
scale!(d::CompoundDistributionWrapper,s::AbstractFloat) = (map!((d1)->rescale(d1,s),d.distributions,d.distributions) ; d)
update!(d::CompoundDistributionWrapper,x::AbstractVector,s::AbstractVector) = (map!(update,d.distributions,d.distributions,x,s) ; d)

###Base.show for DistributionWrappers
function show(io::IO,d::DistributionWrapper)
    println(io,"$(densityname(d)) with distribution:")
    show(io,d.distribution)
    println(io)
    nothing
end


text:
"""
    $SIGNATURES

Create a `Graphs.incdict` object and populate with nodes (states) and edges (transitions)
according to the contents of parameters passed in.

Notes:
- Current implementation repeats duplicate transitions as new edges.
"""
function histGraphStateMachineTransitions(stateVisits, allStates::Vector{Symbol})

  g = Graphs.incdict(Graphs.ExVertex,is_directed=true)
  lookup = Dict{Symbol, Int}()

  # add all required states as nodes to the visualization graph
  fid = 0
  for state in allStates
    fid += 1
    lbl = string(state)
    exvert = Graphs.ExVertex(fid, lbl)
    exvert.attributes["label"] = lbl
    Graphs.add_vertex!(g, exvert)
    lookup[state] = fid
  end

  # add all edges to graph
  count = 0
  for (from, tos) in stateVisits
    for to in tos
      count += 1
      exvf = g.vertices[lookup[from]]
      exvt = g.vertices[lookup[to]]
      # add the edge fom one to the next state
      edge = Graphs.make_edge(g, exvf, exvt)
      Graphs.add_edge!(g, edge)
    end
  end

  return g, lookup
end

function renderStateMachineFrame(vg,
                                 frame::Int;
                                 title::String="",
                                 viewerapp::String="eog",
                                 fext::String="png",
                                 engine::String="dot",
                                 show::Bool=true,
                                 folder::String="fsm_animation",
                                 timest::String="",
                                 rmfirst::Bool=false  )
  #
  folderpath = "/tmp/$folder/"
  if rmfirst
    @warn "removing contents of $(folderpath)"
    Base.rm(folderpath, recursive=true, force=true)
  end
  mkpath(folderpath)
  dotfile = "$folderpath/csm_$frame.dot"
  filepath= "$folderpath/csm_$frame.$(fext)"

  # dot file handle
  fid = open(dotfile,"w")
  write(fid,Graphs.to_dot(vg))
  close(fid)

  # build the dot file somewhat manually
  fid = open("$folderpath/dotscript.sh","w")
  str = "head -n `wc -l $dotfile | awk '{print \$1-1}'` $dotfile > $folderpath/tmpdot.dot"
  println(fid, str)
  println(fid, "echo \"graph [label=\\\"$title, #$step, $(timest)\\\", labelloc=t];\" >> $folderpath/tmpdot.dot")
  println(fid, "echo \"}\" >> $folderpath/tmpdot.dot")
  close(fid)
  run(`chmod u+x $folderpath/dotscript.sh`)
  run(`sh $folderpath/dotscript.sh`)
  Base.rm(dotfile)
  Base.rm("$folderpath/dotscript.sh")
  run(`mv $folderpath/tmpdot.dot $dotfile`)

  # compile output and maybe show to user
  run(`$(engine) $(dotfile) -T$(fext) -o $(filepath)`)
  show ? (@async run(`$viewerapp $filepath`)) : nothing
  return filepath
end

function setVisGraphOnState!(vg, vertid; xlabel::String="", appendxlabel::String="")
  #
  vg.vertices[vertid].attributes["fillcolor"] = "red"
  vg.vertices[vertid].attributes["style"] = "filled"
  if length(xlabel) > 0
    vg.vertices[vertid].attributes["xlabel"] = xlabel
  end
  if haskey(vg.vertices[vertid].attributes, "xlabel")
    vg.vertices[vertid].attributes["xlabel"] = vg.vertices[vertid].attributes["xlabel"]*appendxlabel
  elseif length(appendxlabel) > 0
    vg.vertices[vertid].attributes["xlabel"] = appendxlabel
  end
end

function clearVisGraphAttributes!(vg)
  for (vid,vert) in vg.vertices
    haskey(vert.attributes, "fillcolor") ? delete!(vert.attributes, "fillcolor") : nothing
    haskey(vert.attributes, "style") ? delete!(vert.attributes, "style") : nothing
    haskey(vert.attributes, "xlabel") ? delete!(vert.attributes, "xlabel") : nothing
  end
  nothing
end

function drawStateTransitionStep(hist,
                                 step::Int,
                                 vg,
                                 lookup::Dict{Symbol,Int};
                                 title::String="",
                                 viewerapp::String="eog",
                                 fext::String="png",
                                 engine::String="dot",
                                 show::Bool=true,
                                 folder::String="",
                                 frame::Int=step  )
  #

  lbl = getStateLabel(hist[step][3])
  vertid = lookup[lbl]
  vert = vg.vertices[vertid]

  fillcolorbefore = haskey(vert.attributes, "fillcolor") ? deepcopy(vg.vertices[vertid].attributes["fillcolor"]) : nothing
  stylebefore = haskey(vert.attributes, "style") ? deepcopy(vg.vertices[vertid].attributes["style"]) : nothing
  xlabelbefore = haskey(vg.vertices[vertid].attributes,"xlabel") ? deepcopy(vg.vertices[vertid].attributes["xlabel"]) : nothing
  # delete!(vert.attributes, "fillcolor")
  # delete!(vert.attributes, "style")

  # identify and set the node
  xlabel = length(title) > 0 ? (xlabelbefore != nothing ? xlabelbefore*"," : "")*title : ""
  setVisGraphOnState!(vg, vertid, xlabel=xlabel)

  # render state machine frame
  filepath = renderStateMachineFrame(vg,
                                     frame,
                                     title=title,
                                     viewerapp=viewerapp,
                                     fext=fext,
                                     engine=engine,
                                     show=show,
                                     folder=folder,
                                     timest=string(split(string(hist[step][1]),'T')[end]),
                                     rmfirst=false)
  #

  # clean up the vg structure
  fillcolorbefore == nothing ? delete!(vert.attributes, "fillcolor") : (vert.attributes["fillcolor"]=fillcolorbefore)
  stylebefore == nothing ? delete!(vert.attributes, "style") : (vert.attributes["style"]=stylebefore)
  xlabelbefore == nothing ? delete!(vert.attributes, "xlabel") : (vert.attributes["xlabel"]=xlabelbefore)

  return filepath
end

function drawStateMachineHistory(hist; show::Bool=false, folder::String="")

  stateVisits, allStates = histStateMachineTransitions(hist)

  vg, lookup = histGraphStateMachineTransitions(stateVisits, allStates)

  for i in 1:length(hist)
    drawStateTransitionStep(hist, i, vg, lookup, folder=folder, show=show)
  end

  return nothing
end

"""
    $SIGNATURES

Draw simultaneously separate time synchronized frames from each of the desired
state machines.  These images can be produced into synchronous side-by-side videos
which allows for easier debugging and comparison of concurrently running state
machines.
"""
function animateStateMachineHistoryByTime(hist::Vector{Tuple{DateTime, Int, <: Function, T}};
                                          frames::Int=100,
                                          folder="animatestate",
                                          title::String="",
                                          show::Bool=false,
                                          startT=hist[1][1],
                                          stopT=hist[end][1],
                                          rmfirst::Bool=true  ) where T
  #
  stateVisits, allStates = histStateMachineTransitions(hist)

  vg, lookup = histGraphStateMachineTransitions(stateVisits, allStates)

  totT = stopT - startT
  totT = Millisecond(round(Int, 1.05*totT.value))
  # totT *= 1.05

  step = 1
  len = length(hist)
  @showprogress "exporting state machine images, $title " for i in 1:frames
    aniT = Millisecond(round(Int, i/frames*totT.value)) + startT
    # aniT = i/frames*totT + startT
    if hist[step][1] < aniT && step < len
      step += 1
    end
    drawStateTransitionStep(hist, step, vg, lookup, title=title, folder=folder, show=show, frame=i)
  end

  nothing
end

function animateStateMachineHistoryByTimeCompound(hists::Dict{Symbol, Vector{Tuple{DateTime, Int, <: Function, T}}},
                                                  startT,
                                                  stopT;
                                                  frames::Int=100,
                                                  folder="animatestate",
                                                  title::String="",
                                                  show::Bool=false,
                                                  clearstale::Bool=true,
                                                  rmfirst::Bool=true  ) where T
  #
  # Dict{Symbol, Vector{Symbol}}
  stateVisits = Dict{Symbol, Vector{Symbol}}()
  allStates = Vector{Symbol}()
  for (csym,hist) in hists
    stateVisits, allStates = histStateMachineTransitions(hist,allStates=allStates, stateVisits=stateVisits  )
  end

  #
  vg, lookup = histGraphStateMachineTransitions(stateVisits, allStates)

  # total draw time and step initialization
  totT = stopT - startT
  totT = Millisecond(round(Int, 1.05*totT.value))
  histsteps = ones(Int, length(hists))

  # clear any stale state
  clearstale ? clearVisGraphAttributes!(vg) : nothing

  # loop across time
  @showprogress "exporting state machine images, $title " for i in 1:frames
    # calc frame time
    aniT = Millisecond(round(Int, i/frames*totT.value)) + startT

    # loop over all state machines
    histidx = 0
    for (csym, hist) in hists
      histidx += 1
      step = histsteps[histidx]
      len = length(hist)
      if hist[step][1] < aniT && step < len
        histsteps[histidx] += 1
      end
      # redefine after +1
      step = histsteps[histidx]

      # modify vg for each history
      lbl = getStateLabel(hist[step][3])
      vertid = lookup[lbl]
      setVisGraphOnState!(vg, vertid, appendxlabel=string(csym)*",")
    end

    # finally render one frame
    renderStateMachineFrame(vg,
                            i,
                            title=title,
                            show=false,
                            folder=folder,
                            timest=string(split(string(aniT),' ')[1]),
                            rmfirst=false  )
    #
    clearVisGraphAttributes!(vg)
  end

end

ext: jl


text:
"""
    IndependentMOKernel(k::Kernel)

Kernel for multiple independent outputs with kernel `k` each.

# Definition

For inputs ``x, x'`` and output dimensions ``p_x, p_{x'}'``, the kernel ``\\widetilde{k}``
for independent outputs with kernel ``k`` each is defined as
```math
\\widetilde{k}\\big((x, p_x), (x', p_{x'})\\big) = \\begin{cases}
    k(x, x') & \\text{if } p_x = p_{x'}, \\\\
    0 & \\text{otherwise}.
\\end{cases}
```
Mathematically, it is equivalent to a matrix-valued kernel defined as
```math
\\widetilde{K}(x, x') = \\mathrm{diag}\\big(k(x, x'), \\ldots, k(x, x')\\big) \\in \\mathbb{R}^{m \\times m},
```
where ``m`` is the number of outputs.
"""
struct IndependentMOKernel{Tkernel<:Kernel} <: MOKernel
    kernel::Tkernel
end

function (κ::IndependentMOKernel)((x, px)::Tuple{Any,Int}, (y, py)::Tuple{Any,Int})
    return κ.kernel(x, y) * (px == py)
end

function _kernelmatrix_kron_helper(::MOInputIsotopicByFeatures, Kfeatures, B)
    return kron(Kfeatures, B)
end

function _kernelmatrix_kron_helper(::MOInputIsotopicByOutputs, Kfeatures, B)
    return kron(B, Kfeatures)
end

function kernelmatrix(
    k::IndependentMOKernel, x::MOI, y::MOI
) where {MOI<:IsotopicMOInputsUnion}
    @assert x.out_dim == y.out_dim
    Kfeatures = kernelmatrix(k.kernel, x.x, y.x)
    mtype = eltype(Kfeatures)
    return _kernelmatrix_kron_helper(x, Kfeatures, Eye{mtype}(x.out_dim))
end

if VERSION >= v"1.6"
    function _kernelmatrix_kron_helper!(K, ::MOInputIsotopicByFeatures, Kfeatures, B)
        return kron!(K, Kfeatures, B)
    end

    function _kernelmatrix_kron_helper!(K, ::MOInputIsotopicByOutputs, Kfeatures, B)
        return kron!(K, B, Kfeatures)
    end

    function kernelmatrix!(
        K::AbstractMatrix, k::IndependentMOKernel, x::MOI, y::MOI
    ) where {MOI<:IsotopicMOInputsUnion}
        @assert x.out_dim == y.out_dim
        Ktmp = kernelmatrix(k.kernel, x.x, y.x)
        mtype = eltype(Ktmp)
        return _kernelmatrix_kron_helper!(
            K, x, Ktmp, Matrix{mtype}(I, x.out_dim, x.out_dim)
        )
    end
end

function Base.show(io::IO, k::IndependentMOKernel)
    return print(io, string("Independent Multi-Output Kernel\n\t", string(k.kernel)))
end


text:
using ROOTHistograms
using Base.Test

include("histroot.jl")


text:
# TODO:
# This is a punt for now. It does not handle any sort of custom resolving tags,
# only matching default implicits.

const DEFAULT_SCALAR_TAG = "tag:yaml.org,2002:str"
const DEFAULT_SEQUENCE_TAG = "tag:yaml.org,2002:seq"
const DEFAULT_MAPPING_TAG = "tag:yaml.org,2002:map"

const default_implicit_resolvers =
    [
         ("tag:yaml.org,2002:bool",
          r"^(?:true|True|TRUE|false|False|FALSE)$"x),

         ("tag:yaml.org,2002:int",
          r"^(?:[-+]?0b[0-1_]+
            |[-+]? [0-9]+
            |0o [0-7]+
            |0x [0-9a-fA-F]+)$"x),

         ("tag:yaml.org,2002:float",
          r"^(?:[-+]? ( \. [0-9]+ | [0-9]+ ( \. [0-9]* )? ) ( [eE] [-+]? [0-9]+ )?
            |[-+]? (?: \.inf | \.Inf | \.INF )
            |\.nan | \.NaN | \.NAN)$"x),

         ("tag:yaml.org,2002:merge",
          r"^(?:<<)$"),

         ("tag:yaml.org,2002:null",
          r"^(?:~|null|Null|NULL|)$"x),

         ("tag:yaml.org,2002:timestamp",
          r"^(\d{4})-    (?# year)
             (\d\d?)-    (?# month)
             (\d\d?)     (?# day)
            (?:
              (?:[Tt]|[ \t]+)
              (\d\d?):      (?# hour)
              (\d\d):       (?# minute)
              (\d\d)        (?# second)
              (?:\.(\d*))?  (?# fraction)
              (?:
                [ \t]*(Z|([+\-])(\d\d?)
                  (?:
                      :(\d\d)
                  )?)
              )?
            )?$"x),

         ("tag:yaml.org,2002:value",
          r"^(?:=)$"),

         ("tag:yaml.org,2002:yaml",
          r"^(?:!|&|\*)$")
    ]

mutable struct Resolver
    implicit_resolvers::Vector

    function Resolver()
        new(copy(default_implicit_resolvers))
    end
end

function resolve(resolver::Resolver, ::Type{ScalarNode}, value, implicit)
    if implicit[1]
        for (tag, pat) in resolver.implicit_resolvers
            if occursin(pat, value)
                return tag
            end
        end
    end

    DEFAULT_SCALAR_TAG
end

function resolve(resolver::Resolver, ::Type{SequenceNode}, value, implicit)
    DEFAULT_SEQUENCE_TAG
end

function resolve(resolver::Resolver, ::Type{MappingNode}, value, implicit)
    DEFAULT_MAPPING_TAG
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    NeighborhoodSearch(object, neighborhood)

A method for searching neighbors in spatial `object` inside `neighborhood`.
"""
struct NeighborhoodSearch{O,N,K} <: NeighborSearchMethod
  # input fields
  object::O
  neigh::N

  # state fields
  kdtree::K
end

function NeighborhoodSearch(object::O, neigh::N) where {O,N}
  kdtree = if neigh isa BallNeighborhood
    KDTree(coordinates(object), metric(neigh))
  else
    nothing
  end
  NeighborhoodSearch{O,N,typeof(kdtree)}(object, neigh, kdtree)
end

# search method for any neighborhood
function search(xₒ::AbstractVector, method::NeighborhoodSearch; mask=nothing)
  object = method.object
  neigh  = method.neigh
  N = ncoords(object)
  T = coordtype(object)
  n = nelms(object)

  inds = mask ≠ nothing ? view(1:n, mask) : 1:n

  x = MVector{N,T}(undef)

  neighbors = Vector{Int}()
  @inbounds for ind in inds
    coordinates!(x, object, ind)
    if isneighbor(neigh, xₒ, x)
      push!(neighbors, ind)
    end
  end

  neighbors
end

# search method for ball neighborhood
function search(xₒ::AbstractVector, method::NeighborhoodSearch{O,N,K};
                mask=nothing) where {O,N<:BallNeighborhood,K}
  inds = inrange(method.kdtree, xₒ, radius(method.neigh))
  if mask ≠ nothing
    neighbors = Vector{Int}()
    @inbounds for ind in inds
      if mask[ind]
        push!(neighbors, ind)
      end
    end
    neighbors
  else
    inds
  end
end


text:
module DECAES

using LinearAlgebra, SpecialFunctions, Statistics, Random
import ArgParse, BangBang, Dierckx, DocStringExtensions, ForwardDiff, Logging, MAT, NIfTI, NLopt, Parameters, PolynomialRoots, ProgressMeter, SIMD, StaticArrays, TupleTools, UnPack, UnsafeArrays
using ArgParse: @add_arg_table!, ArgParseSettings, add_arg_group!, add_arg_table!, parse_args
using BangBang: setindex!!, setproperty!!, setproperties!!
using Base.MathConstants: φ
using DocStringExtensions: FIELDS, SIGNATURES, TYPEDFIELDS, TYPEDSIGNATURES
using ForwardDiff: DiffResults
using Logging: ConsoleLogger, with_logger
using Parameters: @with_kw, @with_kw_noshow
using ProgressMeter: Progress, BarGlyphs, finish!, next!
using StaticArrays: FieldVector, SA, SArray, SVector, SMatrix, SizedVector, MVector
using SIMD: FloatingTypes, Vec, shufflevector
using UnPack: @unpack, @pack!
using UnsafeArrays: @uviews, uviews, uview

include("NNLS.jl")
using .NNLS

include("ParXRec.jl")
using .ParXRec

include("NormalHermiteSplines.jl")
using .NormalHermiteSplines

include("types.jl")
include("utils.jl")
include("splines.jl")
include("lsqnonneg.jl")
include("lsqnonneg_reg.jl")
include("EPGdecaycurve.jl")
include("T2mapSEcorr.jl")
include("T2partSEcorr.jl")
include("main.jl")

# Exported symbols
export MAT, NIfTI, ParXRec, load_image
export T2mapSEcorr, T2mapOptions, T2partSEcorr, T2partOptions
export EPGdecaycurve, EPGdecaycurve!, EPGdecaycurve_work
export lsqnonneg, lsqnonneg_chi2, lsqnonneg_gcv, lsqnonneg_lcurve, lcurve_corner
export main, julia_main

# include("precompile.jl")
# _precompile_()

end


text:
#Julia Winston Plotting

using Winston
# fix the random seed so we have reproducible results
srand(111)
# generate a plot
pl = plot(cumsum(rand(100) .- 0.5), "g", cumsum(rand(100) .- 0.5), "b")
# display the plot
display(pl)


text:
"""
    AnovaData

Container for the complete results of an ANOVA test.
"""
mutable struct AnovaData{N}
    effects::Vector{AnovaEffect}
    errors::Vector{<:AnovaEffect}
    ncrossedfactorlevels::Vector{Int}
    npercrossedcell::Int
    crossedcellmeans::Array{Float64, N}
end

import Base.show
function show(io::IO, ad::AnovaData)
    colnames = ["Effect", "SS", "DF", "MS", "F", "p", "ω²"]
    rownames = [e.name for e ∈ ad.effects]
    nrows = length(ad.effects)

    compactshow(x) = sprint(show, x, context=:compact=>true)

    ss = [e.ss |> compactshow for e ∈ ad.effects]
    df = [e.df |> Int |> compactshow for e ∈ ad.effects]
    ms = [typeof(e) ∈ [AnovaFactor, AnovaResult] ? e.ms |> compactshow  : "" for e ∈ ad.effects]
    f = [e isa AnovaResult ? e.f |> compactshow : "" for e ∈ ad.effects]
    p = [e isa AnovaResult ? e.p |> compactshow : "" for e ∈ ad.effects]
    ω² = [e isa AnovaResult ? e.ω² |> compactshow : "" for e ∈ ad.effects]

    ndec = [length.(last.(split.(values, "."))) for values ∈ [ss, ms, f, p, ω²]]
    maxndec = maximum.(ndec)
    rpadding = [maxndec[i] .- ndec[i] for i ∈ 1:5]

    ss = [rpad(ss[i], rpadding[1][i] + length(ss[i])) for i ∈ 1:nrows]
    ms = [rpad(ms[i], rpadding[2][i] + length(ms[i])) for i ∈ 1:nrows]
    f = [rpad(f[i], rpadding[3][i] + length(f[i])) for i ∈ 1:nrows]
    p = [rpad(p[i], rpadding[4][i] + length(p[i])) for i ∈ 1:nrows]
    ω² = [rpad(ω²[i], rpadding[5][i] + length(ω²[i])) for i ∈ 1:nrows]

    colwidths = [length.(values) |> maximum for values ∈ [[colnames[1]; rownames],
                                                          [colnames[2]; ss],
                                                          [colnames[3]; df],
                                                          [colnames[4]; ms],
                                                          [colnames[5]; f],
                                                          [colnames[6]; p],
                                                          [colnames[7]; ω²]]]
    colwidths[2:end] .+= 2

    headerrow = join(lpad.(colnames, colwidths))
    separator = repeat("-", sum(colwidths))
    rows = [lpad(rownames[i], colwidths[1]) *
            lpad(ss[i], colwidths[2]) *
            lpad(df[i], colwidths[3]) *
            lpad(ms[i], colwidths[4]) *
            lpad(f[i], colwidths[5]) *
            lpad(p[i], colwidths[6]) *
            lpad(ω²[i], colwidths[7]) for i ∈ 1:nrows]

    println(io)
    println(io, "Analysis of Variance Results")
    println(io)
    println(io, headerrow)
    println(io, separator)
    foreach(i -> println(io, rows[i]), 1:length(rows))
end


text:
export @symbolic_wrap, @wrapped

# Turn A{X} <: B{Int, X} into
#
# B{Int, X} where X
function set_where(subt, supert)
    if !(supert isa Expr && supert.head == :curly)
        return supert
    end

    Ss = []
    if subt isa Expr && subt.head == :curly
        Ss = subt.args[2:end]
    end

    Ts = intersect(supert.args[2:end], Ss)
    Expr(:where, supert, Ts...)
end

getname(x::Symbol) = x

function getname(x::Expr)
    @assert x.head == :curly
    return x.args[1]
end

macro symbolic_wrap(expr)
    @assert expr isa Expr && expr.head == :struct
    @assert expr.args[2].head == :(<:)
    sig = expr.args[2]
    T = getname(sig.args[1])
    supertype = set_where(sig.args[1], sig.args[2])

    quote
        $expr

        Symbolics.has_symwrapper(::Type{<:$supertype}) = true
        Symbolics.wrapper_type(::Type{<:$supertype}) = $T
        Symbolics.wraps_type(::Type{$T}) = $supertype
        Symbolics.iswrapped(::$T) = true
    end |> esc
end

iswrapped(x) = false
unwrap(x) = x

function wrap(x)
    T = SymbolicUtils.symtype(x)
    Symbolics.has_symwrapper(T) ?
        Symbolics.wrapper_type(T)(x) : x
end

function wrapper_type end
function wraps_type end

has_symwrapper(::Type) = false

function wrap_func_expr(mod, expr)
    @assert expr.head == :function || (expr.head == :(=) &&
                                       expr.args[1] isa Expr &&
                                       expr.args[1].head == :call)

    def = splitdef(expr)

    sig = expr.args[1]
    body = def[:body]

    fname = def[:name]
    args = get(def, :args, [])
    kwargs = get(def, :kwargs, [])

    impl_name = Symbol(fname,"_", hash(string(args)*string(kwargs)))

    function kwargname(kwarg)
        if kwarg isa Expr && kwarg.head == :kw
            kwarg.args[1]
        elseif kwarg isa Expr && kwarg.head == :(...)
            kwarg.args[1]
        else
            kwarg
        end
    end

    function argname(arg)
        if arg isa Expr && (arg.head == :(::) || arg.head == :(...))
            arg.args[1]
        elseif arg isa Expr
            error("$arg not supported as an argument")
        else
            arg
        end
    end

    names = vcat(argname.(args), kwargname.(kwargs))

    function type_options(arg)
        if arg isa Expr && arg.head == :(::)
            T = Base.eval(mod, arg.args[2])
            has_symwrapper(T) ? (T, :(SymbolicUtils.Symbolic{<:$T}), wrapper_type(T))  : (T,)
        elseif arg isa Expr && arg.head == :(...)
            Ts = type_options(arg.args[1])
            map(x->Vararg{x},Ts)
        else
            (Any,)
        end
    end

    types = map(type_options, args)

    impl = :(function $impl_name($(names...))
        $body
    end)
    # TODO: maybe don't drop first lol
    methods = map(Iterators.drop(Iterators.product(types...), 1)) do Ts
        method_args = map(names, Ts) do n, T
            :($n::$T)
        end

        fbody = :(if any($iswrapped, ($(names...),))
                      $wrap($impl_name($([:($unwrap($arg)) for arg in names]...)))
                  else
                      $impl_name($(names...))
                  end)

        if isempty(kwargs)
            :(function $fname($(method_args...))
                  $fbody
              end)
        else
            :(function $fname($(method_args...); $(kwargs...))
                  $fbody
              end)
        end
    end

    quote
        $impl
        $(methods...)
    end |> esc
end

macro wrapped(expr)
    wrap_func_expr(__module__, expr)
end


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# Auxiliary, non-exported functions are declared here.

""" meshgrid(x, y, z)

    Create a z-coordinate matrix from `x`, `y` coordinates and a function `f`,
    such that `z[row,col] = f(x[col], y[row)`"""
function meshgrid(x,y,f)
    z = zeros(length(y),length(x))
    for (yi,yy) in enumerate(y)
        for (xi,xx) in enumerate(x)
            z[yi,xi] = f(xx,yy)
        end
    end
    return z
end

### Debug mode
function debug(msg, f="")
    if config[:debug]
        s = split(msg, "\n", keepempty=false)
        if !isempty(s)
            if !isempty(f)
                printstyled("┌ Gaston in function $f\n", color=:yellow)
            else
                printstyled("┌ Gaston\n", color=:yellow)
            end
            for ss in s
                println("│ $ss")
            end
            printstyled("└\n", color=:yellow)
        end
    end
end

# create x,y coordinates for a histogram, from a sample vector, using a number
# of bins
function hist(s,bins)
    # When adding an element s to a bin, we use an iequality m < s <= M.
    # In order to account for elements s==m, we need to special-case
    # the computation for the first bin
    ms = minimum(s)
    Ms = maximum(s)
    bins = max(bins, 1)
    if Ms == ms
        # compute a "natural" scale
        g = (10.0^floor(log10(abs(ms)+eps()))) / 2
        ms, Ms = ms - g, ms + g
    end
    delta = (Ms-ms)/bins
    x = ms:delta:Ms
    y = zeros(bins)
    # this is special-cased because we want to include the minimum in the
    # first bin
    y[1] = sum(ms .<= s .<= x[2])
    for i in 2:length(x)-2
        y[i] = sum(x[i] .< s .<= x[i+1])
    end
    # this is special-cased because there is no guarantee that x[end] == Ms
    # (because of how ranges work)
    if length(y) > 1 y[end] = sum(x[end-1] .< s .<= Ms) end

    if bins != 1
        # We want the left bin to start at ms and the right bin to end at Ms
        x = (ms+delta/2):delta:Ms
    else
        # add two empty bins on the sides to provide a scale to gnuplot
        x = (ms-delta/2):delta:(ms+delta/2)
        y = [0.0, y[1], 0.0]
    end
    return x,y
end

function Base.display(x::Figure)
    debug("Entering display()")
    isempty(x) && return nothing
    if config[:term] in term_text
        show(stdout, "text/plain", x)
        return nothing
    end
    if config[:mode] == "null"
        return nothing
    end
    llplot(x)
end

Base.show(io::IO, ::MIME"text/plain", x::Figure) = show(io, x)

function Base.show(io::IO, x::Figure)
    debug("showable = $(config[:showable])\nmode = $(config[:mode])", "show():text/plain")

    isempty(x) && return nothing
    config[:mode] == "null" && return nothing
    gnuplot_state.gnuplot_available || return nothing
    config[:term] in term_text || return nothing

    tmpfile = tempname()
    save(term = config[:term], saveopts = config[:termopts], output=tmpfile)
    while !isfile(tmpfile) end  # avoid race condition with read in next line
    write(io, read(tmpfile))
    rm(tmpfile, force=true)
    nothing
end

function Base.show(io::IO, ::MIME"image/svg+xml", x::Figure)
    debug("showable = $(config[:showable])\nmode = $(config[:mode])", "show():image/svg+xml")

    isempty(x) && return nothing
    config[:mode] == "null" && return nothing
    gnuplot_state.gnuplot_available || return nothing

    tmpfile = tempname()
    save(term="svg", saveopts = config[:termopts], output=tmpfile)
    while !isfile(tmpfile) end  # avoid race condition with read in next line
    write(io, read(tmpfile))
    rm(tmpfile, force=true)
    nothing
end

function Base.show(io::IO, ::MIME"image/png", x::Figure)
    debug("showable = $(config[:showable])\nmode = $(config[:mode])", "show():image/png")

    isempty(x) && return nothing
    config[:mode] == "null" && return nothing
    gnuplot_state.gnuplot_available || return nothing

    tmpfile = tempname()
    save(term="png", saveopts = config[:termopts], output=tmpfile)
    while !isfile(tmpfile) end  # avoid race condition with read in next line
    write(io, read(tmpfile))
    rm(tmpfile, force=true)
    nothing
end

# build a string with plot commands according to configuration
function plotstring(sp::SubPlot)
    curves = sp.curves
    file = sp.datafile
    # We have to insert "," between plot commands. One easy way to do this
    # is create the first plot command, then the rest
    # We also need to keep track of the current index (starts at zero)
    p = Vector{String}()
    for (i, curve) in enumerate(curves)
        push!(p, "'$file' i $(i-1) $(curve.conf)")
    end
    cmd = "plot "
    sp.dims == 3 && (cmd = "splot ")
    return cmd*join(p, ", ")
end

# Build a "set term" string appropriate for the terminal type
function termstring(handle ; term=config[:term], termopts=config[:termopts])
    global gnuplot_state

    # build termstring
    ts = "set term $term "
    term ∈ term_window && (ts = ts*string(handle)*" ")
    ts = ts*termopts
    return ts
end

# Calculating palettes is expensive, so store them in a cache. The cache is
# pre-populated with gnuplot's gray palette
Palette_cache = Dict{Symbol, String}(:gray => "set palette gray")

# Calculating linetypes from a colorscheme is expensive, so we use a cache.
Linetypes_cache = Dict{Symbol, String}()

# Convert a symbol to string, converting all '_' to spaces and surrounding it with ' '
function symtostr(s)
    s isa Symbol || return s
    s = string(s)
    return "'$(join(split(s,"_")," "))'"
end

# parse plot configuration
function parse(kwargs)
    # string that will be returned
    curveconf = String[]
    # the `curveconf` argument overrides all others
    if :curveconf in keys(kwargs)
        curveconf = [kwargs[:curveconf]]
    else
        K = keys(kwargs)
        ## These keys are processed first and in order.
        # Keys that take arguments.
        for kw in (:every, :e, :skip, :using, :u, :smooth, :bins)
            if kw in K
                push!(curveconf, " $kw $(kwargs[kw]) ")
                break
            end
        end
        # Keys that don't take arguments.
        if :noautoscale in keys(kwargs) && kwargs[:noautoscale] in (:on, :true, "on", "true")
            push!(curveconf, " noautoscale ")
        end
        # with or plotstyle
        for kw in (:with, :w, :plotstyle)
            if kw in K
                push!(curveconf, " with $(kwargs[kw]) ")
                break
            end
        end
        # pointtype
        for kw in (:pointtype, :pt, :marker)
            if kw in K
                val = kwargs[kw]
                if val isa String
                    push!(curveconf, " pointtype $(pointtypes(val)) ")
                elseif val isa Char
                    push!(curveconf, " pointtype '$val' ")
                else
                    push!(curveconf, " pointtype $val ")
                end
                break
            end
        end
        # pointsize
        for kw in (:pointsize, :ps, :markersize, :ms)
            if kw in K
                push!(curveconf, " pointsize $(kwargs[kw]) ")
            end
        end
        # other curveconf elements with arguments
        for kw in (:linecolor, :lc, :ls, :linestyle, :lt, :linetype, :lw,
                   :linewidth, :fill, :fs, :fillcolor, :fc, :dashtype, :dt,
                   :pointinterval, :pi, :pointnumber, :pn)
            if kw in K
                val = symtostr(kwargs[kw])
                push!(curveconf, " $kw $val ")
            end
        end
        # legend
        for kw in (:legend, :leg, :title, :t)
            if kw in K
                val = symtostr(kwargs[kw])
                push!(curveconf, " title $val ")
                break
            end
        end
        # elements with no arguments
        for kw in (:nohidden3d, :nocontours, :nosurface, :noautoscale)
            if kw in K
                if kwargs[kw] in (true, :on, :true, "on", "true")
                    push!(curveconf, " $kw ")
                end
            end
        end
    end
    return join(curveconf, " ")
end

function parse(axes::Axes)
    axesconf = String[]
    K = keys(axes)
    for k in K
        # axesconf is a string with gnuplot commands
        if k == :axesconf
            push!(axesconf, axes[:axesconf])
            continue
        end
        # palette; code inspired by @gcalderone's Gnuplot.jl
        if k in (:pal, :palette)
            val = axes[k]
            val isa Vector || (val = [val])
            for v in val
                if v isa Symbol
                    if haskey(Palette_cache, v)
                        push!(axesconf, Palette_cache[v])
                        continue
                    end
                    cm = colorschemes[v]
                    colors = String[]
                    for i in range(0, 1, length=length(cm))
                        c = get(cm, i)
                        push!(colors, "$i $(c.r) $(c.g) $(c.b)")
                    end
                    s = "set palette defined ("*join(colors, ", ")*")\nset palette maxcolors $(length(cm))"
                    push!(Palette_cache, (v => s))
                    push!(axesconf, s)
                else
                    push!(axesconf, "set palette $v")
                end
            end
            continue
        end
        # linetype definitions; code inspired by @gcalderone's Gnuplot.jl
        if k in (:lt, :linetype)
            val = axes[k]
            val isa Vector || (val = [val])
            for v in val
                if v isa Symbol
                    if haskey(Linetypes_cache, v)
                        push!(axesconf, Linetypes_cache[v])
                        continue
                    end
                    cm = colorschemes[v]
                    linetypes = String[]
                    for i in 1:length(cm)
                        c = cm[i]
                        s = join(string.( round.(Int, 255 .*(c.r, c.g, c.b)), base=16, pad=2))
                        push!(linetypes, "set lt $i lc rgb '#$s'")
                    end
                    s = join(linetypes,"\n")*"\nset linetype cycle $(length(cm))"
                    push!(Linetypes_cache, (v => s))
                    push!(axesconf, s)
                else
                    push!(axesconf, "set linetype $v")
                end
            end
            continue
        end
        # tics
        if k in (:xtics, :ytics, :ztics, :tics)
            val = axes[k]
            val isa Vector || (val = [val])
            for v in val
                if v isa AbstractRange
                    push!(axesconf,"set $k $(first(v)),$(step(v)),$(last(v))")
                elseif v isa Tuple
                    tics = v[1]
                    labs = v[2]
                    tics isa AbstractRange && (tics = collect(tics))
                    println(tics)
                    s = """("$(labs[1])" $(tics[1])"""
                    for i in 2:length(tics)
                        s *= """, "$(labs[i])" $(tics[i])"""
                    end
                    s *= ")"
                    push!(axesconf,"set $k $s")
                elseif v in (:off, :false, false, "false")
                    push!(axesconf,"unset $k")
                else
                    push!(axesconf,"set $k $v")
                end
            end
            continue
        end
        # axis type
        if k == :axis
            val = symtostr(axes[k])
            if val == "semilogx"
                push!(axesconf, "set logscale x")
            elseif val == "semilogy"
                push!(axesconf, "set logscale y")
            elseif val == "semilogz"
                push!(axesconf, "set logscale z")
            elseif val == "loglog"
                push!(axesconf, "set logscale xyz")
            else
                push!(axesconf, "set $val")
            end
            continue
        end
        # range
        if k in (:xrange, :yrange, :zrange, :cbrange)
            val = axes[k]
            val isa Vector || (val = [val])
            for v in val
                if v isa Vector || v isa Tuple
                    push!(axesconf, "set $k [$(ifelse(isinf(v[1]),*,v[1])):$(ifelse(isinf(v[2]),*,v[2]))]")
                else
                    push!(axesconf, "set $k $v")
                end
            end
            continue
        end
        # view
        if k == :view
            val = axes[k]
            val isa Vector || (val = [val])
            for v in val
                if v isa Tuple
                    push!(axesconf, "set view $(join(v, ", "))")
                else
                    push!(axesconf, "set view $v")
                end
            end
            continue
        end
        # dashtypes
        if k in (:dt, :dashtype)
            val = axes[k]
            val isa Vector || (val = [val])
            for v in val
                push!(axesconf, "set dashtype $val")
            end
            continue
        end
        ### handle remaining keys
        val = axes[k]
        val isa Vector || (val = [val])
        for v in val
            # on/off arguments
            if v in (true, :on, :true, "on", "true")
                push!(axesconf, "set $k")
            elseif v in (false, :off, :false, "false")
                push!(axesconf, "unset $k")
            else
                push!(axesconf, "set $k $(symtostr(v))")
            end
        end
    end
    return join(axesconf, "\n")
end

# Define pointtype synonyms
const pt_syns = Dict("dot"      => 0,
                     "⋅"        => 0,
                     "+"        => 1,
                     "plus"     => 1,
                     "x"        => 2,
                     "*"        => 3,
                     "star"     => 3,
                     "esquare"  => 4,
                     "fsquare"  => 5,
                     "ecircle"  => 6,
                     "fcircle"  => 7,
                     "etrianup" => 8,
                     "ftrianup" => 9,
                     "etriandn" => 10,
                     "ftriandn" => 11,
                     "edmd"     => 12,
                     "fdmd"     => 13
                    )
const pt_syns_aqua = Dict("dot"      => 0,
                          "⋅"        => 0,
                          "+"        => 1,
                          "plus"     => 1,
                          "x"        => 2,
                          "*"        => 3,
                          "star"     => 3,
                          "esquare"  => 4,
                          "edmd"     => 5,
                          "etrianup" => 6
                         )
function pointtypes(pt)
    if config[:term] == "aqua"
        if pt in keys(pt_syns_aqua)
            return pt_syns_aqua[pt]
        end
    else
        if pt in keys(pt_syns)
            return pt_syns[pt]
        end
    end
    return "'$pt'"
end

# write commands to gnuplot's pipe
function gnuplot_send(s)
    debug(s, "gnuplot_send")
    if gnuplot_state.gnuplot_available
        w = write(P.gstdin, s*"\n")
        # check that data was accepted by the pipe
        if !(w > 0)
            @warn "Something went wrong writing to gnuplot STDIN."
            return
        end
        flush(P.gstdin)
    end
end


text:
## Do inference in parallel

using Distributed

## add number of processors desired

addprocs(6)

## Load relevant code into scope for each worker

@everywhere begin
    using Pkg;Pkg.activate(".")
end
@everywhere using JLD2,LinearAlgebra,TransformVariables,LogDensityProblems,MCMCChains,Dates
@everywhere import KenyaCoVSD

## Load data

@load("data/linelist_data_with_pos_neg_20feb_to_27apr_c_age.jld2")
@load("data/serological_data_with_20feb_to_31dec_age.jld2")
@load("data/deaths_data20210521.jld2")
@load("data/N_kenya.jld2")

## Parameter transformations and initial KE guess for HMC

trans = as((ct_min1 = as(Real, 0.4, 0.999),
        R₀ = as(Real, 0.0, 7.0),
        ϵ = as(Real, 0.7, 0.999),
        χ₁=as(Real, 0.1, 30.0),χ₂=as(Real, 0.1, 30.0),
        p_test₁=as(Real, 5e-6, 5e-3),p_test₂ = as(Real, 5e-6, 5e-3),
        p_choose1=as(Real, 0.1, 0.999),
        P_eff=as(Real, 0.5, 0.999),
        schooleffect = as(Real, 0.5, 0.999),
        extra_transmissibility = as(Real, 0.5, 2.5),
        influx_exposed_new_variant = as(Real, 0.0, 10e3),
        p_test_boost = as(Real, 1.0, 3.0),
        χ_boost= as(Real, 1.0, 3.0),
        E₀=as(Real, 0.0, 10e3)))

## Transformation of variables
D = Diagonal(0.1*ones(TransformVariables.dimension(trans)))

## Define scenarios
scenarios = [(σ = 0.0,ω = 1/180,scenario_name = "_no_waning_immunity"),
              (σ = 0.5,ω = 1/180,scenario_name = "_med_susceptibility"),
              (σ = 1.0,ω = 1/180,scenario_name = "_full_susceptibility"),
              (σ = 0.16,ω = 1/90,scenario_name = "_fast_reversion"),
              (σ = 0.16,ω = 1/365,scenario_name = "_slow_reversion"),
              (σ = 1.0,ω = 1/(5*365),scenario_name = "_slow_reversion_to_full_susceptibility")]

##Fit each scenario
@distributed for i = 1:length(scenarios)
    scenario = scenarios[i]
    model = KenyaCoVSD.CoVAreaModel("Nairobi",KenyaCoVSD.ll_twogroup_newvariant_infboost,KenyaCoVSD.priors_twogroup_newvariant_cities;
                case_data = linelist_data_with_pos_neg,
                sero_data = serological_data,
                death_data = deaths_data,
                pop_data = N_kenya)

    KenyaCoVSD.inferparameters!(model,2000,trans,0.05,D;serowaningrate = 0.0,σ = scenario.σ,ω = scenario.ω)

    @save("sensitivity_analysis/Nairobi_model$(scenario.scenario_name).jld2",model)
end


text:
using GeneralizedGenerated
using TupleVectors: chainvec

export testparams
EmptyNTtype = NamedTuple{(),Tuple{}} where T<:Tuple

testparams(d::AbstractMeasure) = testvalue(d)

@inline function testparams(mc::ModelClosure; cfg = NamedTuple(), ctx=NamedTuple())
    gg_call(mc, testparams, cfg, ctx, DropReturn())
end

###############################################################################
# ctx::NamedTuple

@inline function tilde(::typeof(testparams), lens::typeof(identity), xname, x, d, cfg, ctx::NamedTuple)
    xnew = testparams(d)
    ctx′ = merge(ctx, NamedTuple{(dynamic(xname),)}((xnew,)))
    (xnew, ctx′, ctx′)
end

@inline function tilde(::typeof(testparams), lens, xname, x, d, cfg, ctx::NamedTuple)
    xnew = set(x, Lens!!(lens), testparams(d))
    ctx′ = merge(ctx, NamedTuple{(dynamic(xname),)}((xnew,)))
    (xnew, ctx′, ctx′)
end


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

# This implements the functions needed by the BKL and RRR methods
# for the choice of the next move

module DeltaE

using Random
using ExtractMacro
using ..Common
using ..Interface
@static VERSION < v"1.1" && using Future: copy!

include("ArraySets.jl")
using .ArraySets

import .ArraySets: check_consistency

include("DynamicSamplers.jl")
using .DynamicSamplers

export DeltaECache, gen_ΔEcache, check_consistency, compute_staged!, apply_staged!,
       compute_reverse_probabilities!, rand_move, rand_skip, apply_move!,
       get_class_f, get_z, gen_EOcache

# findk(ΔElist, dE) = findfirst(ΔElist, abs(dE))

function binsearchexpr(imin::Int, imax::Int, i::Int, ΔElist::Symbol, dE::Symbol)
    if imax - imin < 10
        ex = :(return 0)
        for j = imax:-1:imin
            ex = quote
                if $dE == $ΔElist[$j]
                    return $j
                else
                    $ex
                end
            end
        end
        return ex
    end
    quote
        cE = $ΔElist[$i]
        cE == $dE && return $i
        if cE < $dE
            $(binsearchexpr(i + 1, imax, (i + 1 + imax) ÷ 2, ΔElist, dE))
        else
            $(binsearchexpr(imin, i - 1, (imin + i - 1) ÷ 2, ΔElist, dE))
        end
    end
end

@generated function findk(ΔElist::NTuple{L,ET}, dE::ET) where {L,ET}
    pre = :(dE = abs(dE))
    ex = binsearchexpr(1, L, (1 + L) ÷ 2, :ΔElist, :dE)
    return quote
        $pre
        @inbounds $ex
    end
end

mutable struct DeltaECache{ET,L}
    N::Int
    ΔElist::NTuple{L,ET}
    ft::Vec
    T::Vec
    z::Float64
    T′::Vec
    z′::Float64
    ascache::Vector{ArraySet}
    pos::IVec
    staged::Vector{NTuple{3,Int}}
    function DeltaECache{ET,L}(X::DiscrGraph, C::Config, ΔElist::NTuple{L,ET}, β::Float64, rrr::Bool) where {ET,L}
        N = getN(X)
        @assert C.N == N
        ascache = [ArraySet(N) for k = 1:2L]
        pos = zeros(Int, N)
        for i = 1:N
            ΔE = delta_energy(X, C, i)
            aki = findk(ΔElist, ΔE)
            #@assert aki > 0 (ΔE, ΔElist)
            upi = ΔE > 0 || (ΔE == 0 && C.s[i] == 1)
            ki = aki + L * upi
            pos[i] = ki
            push!(ascache[ki], i)
            #check_consistency(ascache[ki])
        end
        staged = empty!(Array{NTuple{3,Int}}(undef, N))

        ft = Float64[exp(-β * ΔE) for ΔE in allΔE(X)]
        T = zeros(2L)
        z = 0.0
        for k = 1:2L
            f = get_class_f(ft, k, L)
            x = length(ascache[k]) * f
            z += x
            T[k] = x
        end
        T′ = rrr ? similar(T) : T
        z′ = z
        new{ET,L}(N, ΔElist, ft, T, z, T′, z′, ascache, pos, staged)
    end
end

# @generated function DeltaECache{ET}(X::DiscrGraph{ET}, C::Config, β::Float64, rrr::Bool = true)
#     ΔElist = allΔE(X)
#     Expr(:call, Expr(:curly, :DeltaECache, ET, length(ΔElist)), :X, :C, ΔElist, :β, :rrr)
# end

function DeltaECache(X::DiscrGraph{ET}, C::Config, β::Float64, rrr::Bool = true) where {ET}
    ΔElist = allΔE(X)
    DeltaECache{ET,length(ΔElist)}(X, C, ΔElist, β, rrr)
end

gen_ΔEcache(X::DiscrGraph, C::Config, β::Float64, rrr::Bool = true) = DeltaECache(X, C, β, rrr)

get_z(ΔEcache::DeltaECache) = ΔEcache.z

function check_consistency(ΔEcache::DeltaECache{ET,L}) where {ET,L}
    @extract ΔEcache : ΔElist ascache pos
    for as in ascache
        check_consistency(as)
    end
    for (i,k) in enumerate(pos)
        @assert 1 ≤ k ≤ 2L
        as = ascache[k]
        p = as.pos[i]
        @assert 1 ≤ p ≤ length(as)
        for k1 = 1:2L
            k1 == k && continue
            as1 = ascache[k1]
            @assert as1.pos[i] == 0
        end
    end
end

get_class_f(ΔEcache::DeltaECache{ET,L}, k::Integer) where {ET,L} = get_class_f(ΔEcache.ft, k, L)
get_class_f(ft::Vec, k::Integer, L::Integer) = k > L ? ft[k - L] : 1.0

function rand_skip(ΔEcache::DeltaECache)
    @extract ΔEcache : z N
    return floor(Int, Base.log1p(-rand()) / Base.log1p(-z / N))
end

function rand_move(ΔEcache::DeltaECache{ET,L}) where {ET,L}
    @extract ΔEcache: ΔElist ascache ft T z
    r = rand() * z
    k = 0
    cT = 0.0
    for outer k = 1:2L
        cT += T[k]
        r < cT && break
    end
    r < cT || while T[k] == 0
        k -= 1
    end

    if k ≤ L
        ΔE = -ΔElist[k]
    else
        ΔE = ΔElist[k - L]
    end
    move = rand(ascache[k])

    return move, ΔE
end

function apply_staged!(ΔEcache::DeltaECache)
    @extract ΔEcache : ΔElist ascache pos staged

    for (j,k0,k1) in staged
        #@assert k0 ≠ k1
        as0 = ascache[k0]
        as1 = ascache[k1]
        delete!(as0, j)
        push!(as1, j)
        pos[j] = k1
    end
    ΔEcache.T, ΔEcache.T′, ΔEcache.z = ΔEcache.T′, ΔEcache.T, ΔEcache.z′
    return ΔEcache
end

function compute_reverse_probabilities!(ΔEcache::DeltaECache{ET,L}) where {ET,L}
    @extract ΔEcache : staged ft T T′ z

    z′ = z
    T′ ≢ T && copy!(T′, T)
    for (_,k0,k1) in staged
        f0 = get_class_f(ft, k0, L)
        f1 = get_class_f(ft, k1, L)

        T′[k0] -= f0
        T′[k1] += f1
        z′ += f1 - f0
    end
    ΔEcache.z′ = z′

    return z′
end

function compute_staged!(X::DiscrGraph{ET}, C::Config, i::Int, ΔEcache::DeltaECache{ET,L}) where {ET,L}
    @extract C : N s
    @extract ΔEcache : ΔElist pos staged

    spinflip!(X, C, i)
    empty!(staged)

    @inbounds begin
        for j in neighbors(X, i)
            k0 = pos[j]
            upj0 = k0 > L

            dE1 = delta_energy(X, C, j)
            ak1 = findk(ΔElist, dE1)
            upj1 = dE1 > 0 || (dE1 == 0 && s[j] == 1)

            k1 = ak1 + L * upj1

            k0 == k1 && continue

            push!(staged, (j,k0,k1))
        end
        k0 = pos[i]
        k1 = k0 - L * (2 * (k0 > L) - 1)
        push!(staged, (i,k0,k1))
    end

    spinflip!(X, C, i)
end

function apply_move!(X::Union{DiscrGraph{ET},DoubleGraph{DiscrGraph{ET}}}, C::Config,
                     move::Int, ΔEcache::DeltaECache{ET,L}) where {ET,L}
    ## equivalent to:
    #
    # compute_staged!(X, C, move, ΔEcache)
    # compute_reverse_probabilities!(ΔEcache)
    # spinflip!(X, C, move)
    # apply_staged!(ΔEcache)

    @extract C : s
    @extract ΔEcache : ΔElist ascache pos ft T z

    spinflip!(X, C, move)

    X0 = inner_graph(X)

    z′ = z
    @inbounds begin
        for j in neighbors(X0, move)
            k0 = pos[j]
            upj0 = k0 > L

            dE1 = delta_energy(X0, C, j)
            ak1 = findk(ΔElist, dE1)
            upj1 = dE1 > 0 || (dE1 == 0 && s[j] == 1)

            k1 = ak1 + L * upj1

            k0 == k1 && continue

            f0 = get_class_f(ft, k0, L)
            f1 = get_class_f(ft, k1, L)

            T[k0] -= f0
            T[k1] += f1
            z′ += f1 - f0

            as0 = ascache[k0]
            as1 = ascache[k1]
            delete!(as0, j)
            push!(as1, j)
            pos[j] = k1
        end
        k0 = pos[move]
        k1 = k0 - L * (2 * (k0 > L) - 1)

        f0 = get_class_f(ft, k0, L)
        f1 = get_class_f(ft, k1, L)

        T[k0] -= f0
        T[k1] += f1
        z′ += f1 - f0

        as0 = ascache[k0]
        as1 = ascache[k1]
        delete!(as0, move)
        push!(as1, move)
        pos[move] = k1
    end

    c = z / z′
    ΔEcache.z = z′
    return c
end

prior(x) = x > 0 ? exp(-x) : 1.0

mutable struct DeltaECacheCont{ET}
    dynsmp::DynamicSampler
    ΔEs::Vector{ET}
    β::Float64
    staged::Vector{Tuple{Int,ET,Float64}}
    function DeltaECacheCont{ET}(X::AbstractGraph, C::Config, β::Float64) where {ET}
        N = getN(X)
        @assert C.N == N
        ΔEs = [delta_energy(X, C, i) for i = 1:N]
        dynsmp = DynamicSampler(prior(β * ΔE) for ΔE in ΔEs)
        staged = empty!(Array{Tuple{Int,ET,Float64}}(undef, N))
        return new{ET}(dynsmp, ΔEs, β, staged)
    end
end

# the rrr argument here is just for consistency but it's unused
gen_ΔEcache(X::AbstractGraph{ET}, C::Config, β::Float64, rrr::Bool = true) where {ET} = DeltaECacheCont{ET}(X, C, β)

get_z(ΔEcache::DeltaECacheCont) = ΔEcache.dynsmp.z

function rand_skip(ΔEcache::DeltaECacheCont)
    @extract ΔEcache : dynsmp
    @extract dynsmp : z N
    b = clamp(z / N, floatmin(Float64), 1.0)
    return floor(Int, log1p(-rand()) / log1p(-b))
end

function rand_move(ΔEcache::DeltaECacheCont)
    @extract ΔEcache: dynsmp ΔEs

    move = rand(dynsmp)
    ΔE = ΔEs[move]
    return move, ΔE
end

function apply_staged!(ΔEcache::DeltaECacheCont)
    @extract ΔEcache : dynsmp ΔEs staged

    @inbounds for (j,ΔE,p) in staged
        ΔEs[j] = ΔE
        dynsmp[j] = p
    end
    return ΔEcache
end

function compute_reverse_probabilities!(ΔEcache::DeltaECacheCont)
    @extract ΔEcache : dynsmp staged
    @extract dynsmp : z N

    @inbounds for (j,_,p) in staged
        z += p - dynsmp[j]
    end
    z = clamp(z, floatmin(Float64), N)

    return z
end

function compute_staged!(X::AbstractGraph{ET}, C::Config, i::Int, ΔEcache::DeltaECacheCont{ET}) where {ET}
    @extract C : N s
    @extract ΔEcache : β staged

    spinflip!(X, C, i)
    empty!(staged)

    ΔE = delta_energy(X, C, i)
    p = prior(β * ΔE)
    push!(staged, (i,ΔE,p))
    for j in neighbors(X, i)
        ΔE = delta_energy(X, C, j)
        p = prior(β * ΔE)
        push!(staged, (j,ΔE,p))
    end

    spinflip!(X, C, i)
end

get_inner(X, ::Type{Val{true}}) = inner_graph(X)
get_inner(X, ::Type{Val{false}}) = X

function apply_move!(X::AbstractGraph{ET}, C::Config,
                     move::Int, ΔEcache::DeltaECacheCont{ET},
                     inner = Val{true}) where {ET}
    ## equivalent to:
    #
    # compute_staged!(X, C, move, ΔEcache)
    # compute_reverse_probabilities!(ΔEcache)
    # spinflip!(X, C, move)
    # apply_staged!(ΔEcache)

    @extract C : s
    @extract ΔEcache : dynsmp ΔEs β

    spinflip!(X, C, move)

    X0 = get_inner(X, inner)

    z = dynsmp.z
    @inbounds begin
        ΔE = delta_energy(X0, C, move)
        ΔEs[move] = ΔE
        dynsmp[move] = prior(β * ΔE)

        for j in neighbors(X0, move)
            ΔE = delta_energy(X0, C, j)
            ΔEs[j] = ΔE
            dynsmp[j] = prior(β * ΔE)
        end
    end
    z′ = dynsmp.z
    c = z / z′
    return c
end

function findks(ΔElist, ΔE, L, has_zero)
    ak = findk(ΔElist, ΔE)
    #@assert aki > 0 (ΔE, ΔElist)
    if ΔE ≥ 0
        k = ak + L - has_zero
    else
        k = L + 1 - ak
    end
    return k
end

mutable struct EOCache{ET,L}
    N::Int
    ΔElist::NTuple{L,ET}
    has_zero::Bool
    fτ::Vec
    z::Float64
    ascache::Vector{ArraySet}
    pos::IVec
    function EOCache{ET,L}(X::DiscrGraph, C::Config, ΔElist::NTuple{L,ET}, τ::Float64) where {ET,L}
        N = getN(X)
        @assert C.N == N
        has_zero = zero(ET) ∈ ΔElist
        K = 2L - has_zero
        ascache = [ArraySet(N) for k = 1:K]
        pos = zeros(Int, N)
        for i = 1:N
            ΔE = delta_energy(X, C, i)
            ki = findks(ΔElist, ΔE, L, has_zero)
            pos[i] = ki
            push!(ascache[ki], i)
            #check_consistency(ascache[ki])
        end

        fτ = cumsum([j^(-τ) for j = 1:N])
        z = fτ[end]

        new{ET,L}(N, ΔElist, has_zero, fτ, z, ascache, pos)
    end
end

# @generated function EOCache{ET}(X::DiscrGraph{ET}, C::Config, τ::Float64)
#     ΔElist = allΔE(X)
#     Expr(:call, Expr(:curly, :EOCache, ET, length(ΔElist)), :X, :C, ΔElist, :τ)
# end

function EOCache(X::DiscrGraph{ET}, C::Config, τ::Float64) where {ET}
    ΔElist = allΔE(X)
    EOCache{ET,length(ΔElist)}(X, C, ΔElist, τ)
end

gen_EOcache(X::DiscrGraph, C::Config, τ::Float64) = EOCache(X, C, τ)

function check_consistency(ΔEcache::EOCache{ET,L}) where {ET,L}
    @extract ΔEcache : ΔElist has_zero ascache pos
    for as in ascache
        check_consistency(as)
    end
    K = 2L - has_zero
    for (i,k) in enumerate(pos)
        @assert 1 ≤ k ≤ K
        as = ascache[k]
        p = as.pos[i]
        @assert 1 ≤ p ≤ length(as)
        for k1 = 1:K
            k1 == k && continue
            as1 = ascache[k1]
            @assert as1.pos[i] == 0
        end
    end
end

function rand_move(ΔEcache::EOCache{ET,L}) where {ET,L}
    @extract ΔEcache: N ΔElist has_zero ascache pos fτ z
    K = 2L - has_zero
    r = (1 - rand()) * z

    i = searchsortedfirst(fτ, r)
    # @show i, map(length, ascache)
    @assert 1 ≤ i ≤ N

    # i = 0
    # c = 0.0
    # for i = 1:N
    #     c += fτ[i]
    #     r < c && break
    # end
    # r < c || while fτ[i] == 0
    #     i -= 1
    # end

    k = 0
    t = 0
    while i > t
        k += 1
        t += length(ascache[k])
    end

    if k ≤ L
        ΔE = -ΔElist[L + 1 - k]
    else
        ΔE = ΔElist[k - L + has_zero]
    end
    move = rand(ascache[k])

    return move, ΔE
end

function apply_move!(X::DiscrGraph{ET}, C::Config,
                     move::Int, ΔEcache::EOCache{ET,L}) where {ET,L}
    @extract C : s
    @extract ΔEcache : ΔElist has_zero ascache pos fτ z

    spinflip!(X, C, move)

    @inbounds begin
        for j in neighbors(X, move)
            k0 = pos[j]
            dE1 = delta_energy(X, C, j)
            k1 = findks(ΔElist, dE1, L, has_zero)

            k0 == k1 && continue

            as0 = ascache[k0]
            as1 = ascache[k1]
            delete!(as0, j)
            push!(as1, j)
            pos[j] = k1
        end
        k0 = pos[move]
        dE1 = delta_energy(X, C, move)
        k1 = findks(ΔElist, dE1, L, has_zero)

        if k0 ≠ k1
            as0 = ascache[k0]
            as1 = ascache[k1]
            delete!(as0, move)
            push!(as1, move)
            pos[move] = k1
        end
    end
end

# WARNING: very sub-optimal...
mutable struct EOCacheCont{ET}
    N::Int
    ΔEs::Vector{ET}
    rank::IVec
    fτ::Vec
    z::Float64
    function EOCacheCont{ET}(X::AbstractGraph, C::Config, τ::Float64) where {ET}
        N = getN(X)
        @assert C.N == N
        ΔEs = [delta_energy(X, C, i) for i = 1:N]
        rank = sortperm(ΔEs)

        fτ = cumsum([j^(-τ) for j = 1:N])
        z = fτ[end]

        new{ET}(N, ΔEs, rank, fτ, z)
    end
end

gen_EOcache(X::AbstractGraph{ET}, C::Config, τ::Float64) where {ET} = EOCacheCont{ET}(X, C, τ)

function rand_move(ΔEcache::EOCacheCont{ET}) where {ET}
    @extract ΔEcache: N ΔEs rank fτ z
    r = (1 - rand()) * z

    i = searchsortedfirst(fτ, r)
    # @show i, map(length, ascache)
    @assert 1 ≤ i ≤ N

    move = rank[i]
    ΔE = ΔEs[move]

    return move, ΔE
end

function apply_move!(X::AbstractGraph{ET}, C::Config, move::Int, ΔEcache::EOCacheCont{ET}) where {ET}
    @extract C : s
    @extract ΔEcache : ΔEs rank fτ z

    spinflip!(X, C, move)

    @inbounds begin
        ΔE = delta_energy(X, C, move)
        ΔEs[move] = ΔE

        for j in neighbors(X, move)
            ΔE = delta_energy(X, C, j)
            ΔEs[j] = ΔE
        end

        sortperm!(rank, ΔEs, initialized=true)
        rankshuffle!(rank, ΔEs)
    end
end

function rankshuffle!(rank, ΔEs)
    N = length(rank)
    @assert length(ΔEs) == N
    # @assert isperm(rank)
    # @assert issorted(ΔEs[rank])
    i0, i1 = 1, 2
    while i0 < N
        ΔE = ΔEs[rank[i0]]
        while i1 ≤ N && ΔEs[rank[i1]] == ΔE
            i1 += 1
        end
        if (i1-1) > i0
            rank[i0:(i1-1)] = shuffle!(rank[i0:i1-1])
        end
        i0 = i1
        i1 += 1
    end
    if (i1-1) > i0
        rank[i0:(i1-1)] = shuffle!(rank[i0:i1-1])
    end
    # @assert isperm(rank)
    # @assert issorted(ΔEs[rank])
    return rank
end

end # module


text:
module MolarWeight

    # Atomic weight data collection
    include("Data.jl")

    # Internal functions
    include("getaw.jl")
    include("askmform.jl")

    # Exported functions
    include("molarweght.jl")
    export molarweight

end


text:
module GridWorld
    include("common.jl")
    include("aspect_ratio.jl")
    include("color.jl")
    include("margins.jl")
    include("offset.jl")
    include("label.jl")
    include("table.jl")
    include("cell.jl")
    include("svg.jl")
    include("grid.jl")
    include("table_cell.jl")

    export Grid, Margins, Color, getwidth, getheight, getlinecolor, getnumrows,
    getnumcols, getcellsize, Label, Cell, cell, Table, getrowheights,
    getcolwidths, getrowheight, getcolheight, tablecell, TableCell, Piece,
    AbstractGridState, AbstractGridGame, getpieceidx, getpieces, getpiecesdict,
    getstate, getpiece, statetype, SVG, Drawable, Offset, draw, savesvg

    export func
    """
    func(x)

    Returns double the number `x` plus `1`.
    """
    func(x) = 2x + 1

end


text:
export ContinuousSpace, MultiContinuousSpace, DiscreteSpace, MultiDiscreteSpace

import StatsBase:sample

include("continuous_space.jl")
include("discrete_space.jl")
include("multi_continuous_space.jl")
include("multi_discrete_space.jl")


text:
const test_legacy = false

using BenchmarkTools

using Strs
using Strs:V6_COMPAT

using REPL

@api test StrAPI, StrBase, StrRegex

const ver = "v0.$(VERSION.minor)"
const git = "https://github.com/JuliaString/"

const pkgdir    = dirname(pathof(Strs))
const RC        = REPL.REPLCompletions
_uninit(T, len) = T(undef, len)

@static test_legacy && (using LegacyStrings)

const STRS_SETUP = true


text:
module BrainWaveExperimental

export calcDelayedPeakPoints, calcPeaksFromExpectedLatencies, calcPeaksInflFromExpectedLatencies, calcPeaksInflFromExpectedLatencies2, calcPeaksInflFromExpectedLatencies3, calcPeaksTroughsInflFromExpectedLatencies,
findABRWaves, findInflectionsDown, findInflectionsUp,
selectLargestPeakInWindow, selectLargestTroughInWindow, selectStrongestInflectionInWindow

using BrainWave, DataFrames, DocStringExtensions

include("findABRPeaks.jl")
include("utils.jl")

end #module


text:
using OrdinaryDiffEq, DiffEqCallbacks
using DiffEqSensitivity, Zygote, Test
import ForwardDiff
import FiniteDiff

abstol = 1e-6
reltol = 1e-6
savingtimes = 0.1

function test_discrete_callback(cb, tstops, g)
  function fiip(du, u, p, t)
    #du[1] = dx = p[1]*u[1]
    du[:] .= p[1]*u
  end

  p = Float64[0.8123198]
  u0 = Float64[1.0]

  prob = ODEProblem(fiip, u0, (0.0, 1.0), p)

  @show g(solve(prob, Tsit5(), callback=cb, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes))

  du01, dp1 = Zygote.gradient(
    (u0, p) -> g(solve(prob, Tsit5(), u0=u0, p=p, callback=cb, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes,
      sensealg=ForwardDiffSensitivity(;convert_tspan=true))),
    u0, p)

  dstuff1 = ForwardDiff.gradient(
    (θ) -> g(solve(prob, Tsit5(), u0=θ[1:1], p=θ[2:2], callback=cb, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes)),
    [u0; p])

  dstuff2 = FiniteDiff.finite_difference_gradient(
    (θ) -> g(solve(prob, Tsit5(), u0=θ[1:1], p=θ[2:2], callback=cb, tstops=tstops, abstol=abstol, reltol=reltol, saveat=savingtimes)),
    [u0; p])

  @show du01 dp1 dstuff1 dstuff2
  @test du01 ≈ dstuff1[1:1] atol=1e-6
  @test dp1 ≈ dstuff1[2:2] atol=1e-6
  @test du01 ≈ dstuff2[1:1] atol=1e-6
  @test dp1 ≈ dstuff2[2:2] atol=1e-6
end

@testset "ForwardDiffSensitivity: Discrete callbacks" begin
  g(u) = sum(u.^2)
  @testset "reset to initial condition" begin
    affecttimes = range(0.0, 1.0, length=6)[2:end]
    u0 = [1.0]
    condition(u, t, integrator) = t ∈ affecttimes
    affect!(integrator) = (integrator.u .= u0; @show "triggered!")
    cb = DiscreteCallback(condition, affect!, save_positions=(false,false))
    test_discrete_callback(cb, affecttimes, g)
  end
end


text:
# install Python dependencies
ENV["PYTHON"] = ""
using Pkg
Pkg.activate(".")
Pkg.add(["PyCall", "Conda"])
Pkg.build("PyCall")
using Conda
Conda.add("pandas")
Conda.add("pytables")

@testset "HDF5" begin
    import Pandas
    fname = joinpath(testpath, "test.hdf")
    write_table!(fname, "/data", df)
    @test filesize(fname) > 0
    @test list_tables(fname) == ["/data"]
    df_recovered = DataFrame(read_table(fname, "/data"); copycols=false)
    @test compare_df_ignore_order(df,df_recovered)
    fname = joinpath(testpath, "test2.hdf")
    write_table!(fname, "my_data", nt)
    @test filesize(fname) > 0
    nt_recovered = DataFrame(read_table(fname, "my_data"); copycols=false)
    @test compare_df_ignore_order(DataFrame(nt), nt_recovered)
end


text:
require("DataFrames")
using DataFrames

##############################################################################
#
# Introduction to Basic Types and Constructors
#
##############################################################################

#
# Missingness is indicated using a new type: NAtype
#

NAtype
typeof(NAtype)
super(NAtype)
dump(NAtype)

#
# There is a single value that has the type NAtype: it can
# created using either NAtype() or the shorter alias NA.
#

NAtype()
NA
typeof(NA)
dump(NA)

#
# As we'll see later, NA poisons all other scalar objects and
# recycles during interactions with non-scalar objects to poison
# them entrywise
#
# For now, we'll just demonstrate a single example of that poisoning
# behavior
#

NA + 1

#
# NA's represent missingness for scalar data. To express Arrays
# with missingness, we've define a new set of DataArray objects,
# a new parametric composite type that wraps a standard Julia array
# and complements this Julia array with an additional bit vector
# specifying which entries of the primary vector are missing.
#
# These two core components are stored in a DataArray's data and isna
# fields respectively. To save space, the isna field is always a BitArray.
#

DataArray
typeof(DataArray)
typeof(DataArray{Int64})
super(DataArray)
super(super(DataArray))
dump(DataArray)
DataArray.names

#
# The fundamental constructor for a DataArray specifies the data
# values and the missingness metadata
#

dv = DataArray([1, 2, 3], falses(3))

#
# A simpler type of constructor assumes that no data is missing and
# only specifies the vector of non-missing values.
#

dv = DataArray([1, 2, 3])

#
# Yet another constructor, which is conveniently brief, converts
# a Julian Range1 object into a standard Julia Vector and wraps
# this in a DataArray.
#

dv = DataArray(1:3)

#
# To order to provide a very concise constructor, we've hacked in
# a very cute constructor for DataVector's using a trick involving
# their getindex() method.
#

dv = DataVector[1, 2, NA, 4]

#
# Because the missingness of a DataVector is stored in metadata rather than
# using existing values like NaN's, it is easy to create DataVector's
# that store arbitrary Julia types like ComplexPair's and Bool's.
#

dv = DataArray([1 + 2im, 3 - 1im])
dv = DataArray([true, false])

#
# Sometimes you want to convert all of the values inside a DataVector to
# a new type while preserving missingness. Where Julia has specialized
# conversion functions, DataArray's mimic them. Use the data* prefix to find
# them in the REPL
#

dv = DataArray([1.1, 2.1])

databool(dv)
dataint(dv)
datafloat(dv)

#
# Julia also provides convenience functions for generating initialized
# vectors with standard values.
#

datazeros(5)
datazeros(Int64, 5)
dataones(5)
dataones(Int64, 5)
datafalses(5)
datatrues(5)

#
# If you know the type of the DataArray you want to create, but not the values
# you can create DataArray's of length N that are NA everywhere
#

dv = DataArray(Int64, 5)
dv = DataArray(Float64, 5)
dv = DataArray(ComplexPair, 5)

#
# While DataArray's are a very powerful tool for dealing with missing data,
# they only bring us part of the way towards representing real-world data
# sets in Julia. The missing piece is a tabular data structure like those
# found in relational databases or spreadsheets.
#
# To represent these kinds of data sets, we provide the DataFrame type. The
# DataFrame type is a new composite type with two fields:
#
# * columns: A Julia Vector{Any}, each element of which is a single
#            column. The typical column is a DataVector.
#
# * colindex: An Index object that allows one to access entries in the
#             columns. An Index can be constructed from an array of
#             ByteStrings. The details of the Index type are described
#             later.
#
# The DataFrame protocol requires that a DataFrame df behave as follows:
#
# * df is a table with M rows and N columns.
#
# * Every column of df has its own type, which is why a DataFrame is not
#   a matrix of DataVector's.
#
# * All columns of df are guaranteed to be of length M.
#
# * All columns of df are guaranteed to be capable of storing NA if it
#   is ever inserted. (NB: There is ongoing debate about whether a column
#	should always be a DataVector or whether it should only be converted to
#   a DataVector if an NA is introduced by assignment operations.)
#

df_columns = {datazeros(5), datafalses(5)}
df_colindex = Index(["A", "B"])

df = DataFrame(df_columns, df_colindex)

#
# In practice, many other constructors are more convenient to use. The
# simplest is to provide only the columns, which will produce default
# names for all the columns.
#

df = DataFrame(df_columns)

#
# One often would like to construct DataFrame's from columns which may
# not yet be DataVector's. This is possible using the same type of
# constructor. All columns that are not yet DataVector's will be converted
# to DataVector's. (THIS LAST CLAIM IS NOT YET TRUE AND TOM DOES NOT WANT IT.)
#

df = DataFrame({ones(5), zeros(5)})

#
# Often one wishes to convert an existing matrix into a DataFrame.
# This is also possible.
#

df = DataFrame(ones(5, 3))

#
# Like DataVector's, it is possible to create empty DataFrame's in which
# all of the default values are NA.
#
# In the simplest version, we specify a type, the number of rows
# and the number of columns.
#

df = DataFrame(Int64, 10, 5)

#
# Alternatively, one can specify a Vector of types. This implicitly
# defines the number of columns, but one must still specify the number
# of rows.
#

df = DataFrame({Int64, Float64}, 4)

#
# When you know what the names of the columns will be, but not the values
# it is possible to specify them as well.
#
# SHOULD THIS BE DataFrame(types, nrow, names) INSTEAD?
#

DataFrame({Int64, Float64}, ["A", "B"], 10)

# DataFrame({Int64, Float64}, Index(["A", "B"]), 10) NEED TO MAKE THIS WORK

#
# A more uniquely Julian way of creating DataFrame's exploits Julia's ability to
# quote Expressions in order to produce behavior like R's delayed evaluation
# strategy.
#

df = DataFrame(quote
				A = rand(5)
				B = datatrues(5)
			   end)

##############################################################################
#
# Indexing and Assignment Elements
#
##############################################################################

#
# Because a DataVector is a 1-dimensional Array, indexing into it is trivial and
# behaves exactly like indexing into a standard Julia vector.
#

dv = dataones(5)
dv[1]
dv[5]
dv[end]
dv[1:3]
dv[[true, true, false, false, false]]

dv[1] = 3
dv[5] = 5.3
dv[end] = 2.1
dv[1:3] = [3.2, 3.2, 3.1]
# dv[[true, true, false, false, false]] = dataones(2) # SHOULD WE MAKE THIS WORK?

#
# In contrast, a DataFrame is a random-access data structure that can be
# indexing into and assigned to in many different ways. We walk through
# them one-by-one below.
#

#
# Simple numeric indexing
#

df = DataFrame(Int64, 5, 3)
df[1, 3]
df[1]

#
# Range numeric indexing:
#

df = DataFrame(Int64, 5, 3)

df[1, :]
df[:, 3]
df[1:2, 3]
df[1, 1:3]
df[:, :]

#
# Column name indexing
#

df["x1"]
df[1, "x1"]
df[1:3, "x1"]

df[["x1", "x2"]]
df[1, ["x1", "x2"]]
df[1:3, ["x1", "x2"]]

##############################################################################
#
# Unary Operators for NA, DataVector's and DataFrame's
#
# +
# -
# !
#
# MISSING: The transpose unary operator
#
##############################################################################

+NA
-NA
!NA

+dataones(5)
-dataones(5)
!datafalses(5)

##############################################################################
#
# Binary Operators
#
# * Arithmetic Operators
#   +, .+, -, .-, *, .*, /, ./, ^, .^
#
# * Bit Operators
#   MISSING
#
# * Comparison Operators
#   ==, !=, <, <=, >, >=, .==, .!=, .<, .<=, .>, .>=
#
##############################################################################

#
# The standard arithmetic operators work on DataVector's when
# they interact with Number's, NA's or other DataVector's.
#

dv = dataones(5)
dv[1] = NA
df = DataFrame(quote
				 a = 1:5
			   end)

#
# NA's with NA's
#

NA + NA
NA .+ NA
# And so on for -, .- , *, .*, /, ./, ^, .^

#
# NA's with scalars and scalars with NA's
#

1 + NA
1 .+ NA
NA + 1
NA .+ 1
# And so on for -, .- , *, .*, /, ./, ^, .^

#
# NA's with DataVector's
#

dv + NA
dv .+ NA
NA + dv
NA .+ dv
# And so on for -, .- , *, .*, /, ./, ^, .^

#
# DataVector's with scalars
#

dv + 1
dv .+ 1
# And so on for -, .- , .*, ./, .^

#
# Scalars with DataVector's
#

1 + dv
1 .+ dv
# And so on for -, .- , *, .*, /, ./, ^, .^

# HOW MUCH SHOULD WE HAVE OPERATIONS W/ DATAFRAMES?

NA + df
df + NA
1 + df
df + 1
# dv + df # SHOULD THIS EXIST?
# df + dv # SHOULD THIS EXIST?
df + df
# And so on for -, .- , .*, ./, .^

#
# The standard bit operators work on DataVector's
#

# MISSING

#
# The standard comparison operators work on DataVector's
#

NA .< NA
NA .< "a"
NA .< 1
NA .== dv

dv .< NA
dv .< "a"
dv .< 1
dv .== dv

df .< NA
df .< "a"
df .< 1
# df .== dv # SHOULD THIS EXIST?
df .== df

##############################################################################
#
# Elementwise Functions
#
# abs, sign, acos, acosh, asin, asinh, atan, atan2, atanh, sin, sinh
# cos, cosh, tan, tanh, ceil, floor, round, trunc, signif, exp, log
# log10, log1p, log2, exponent, sqrt
#
##############################################################################

#
# Standard functions that apply to scalar values of type Number
# return NA when applied to NA.
#

abs(NA)

#
# Standard functions are broadcast to the elements of DataVector's and
# DataFrame's for elementwise application.
#

dv = dataones(5)
df = DataFrame({dv})

abs(dv)
abs(df)

##############################################################################
#
# Pairwise Functions
#
# diff
#
##############################################################################

#
# Functions that operate on pairs of entries of a Vector work on DataVector's
# and insert NA where it would be produced by other operator rules
#

diff(dv)

##############################################################################
#
# Cumulative Functions
#
# cumprod, cumsum, cumsum_kbn
#
# MISSING: cummin, cummax
#
##############################################################################

#
# Functions that operate cumulatively on the entries of a Vector work on
# DataVector's and insert NA where it would be produced by other operator rules
#

cumprod(dv)
cumsum(dv)
cumsum_kbn(dv)

##############################################################################
#
# Aggregative Functions
#
# min, max, prod, sum, mean, median, std, var, fft, norm
#
##############################################################################

min(dv)

#
# To broadcast these to individual columns, use the col*s versions:
# colmins, colmaxs, colprods, colsums, colmeans, colmedians,
# colstds, colvars, colffts, colnorms
#

colmins(df)

##############################################################################
#
# Loading Standard Data Sets
#
##############################################################################

#
# The DataFrames package is easiest to explore if you also install the
# RDatasets package, which provides access to 570 classic data sets.
#

require("RDatasets")

iris = RDatasets.data("datasets", "iris")
dia = RDatasets.data("ggplot2", "diamonds")

##############################################################################
#
# Split-Apply-Combine
#
##############################################################################

#
# The basic mechanism for spliting data is the groupby() function,
# which will produce a GroupedDataFrame() object that is easiest to
# interact with by iterating over its entries
#

for df in groupby(iris, "Species")
	println("A DataFrame with $(nrow(df)) rows")
end

#
# The | (pipe) operator for GroupedDataFrame's allows you to run
# simple functions on the columns of the induced DataFrame's.
#
# You pass a simple function by producing a symbol with its name.
#

groupby(iris, "Species") | :mean

#
# Another simple way to split-and-apply (without clear combining) is to
# use the map() function
#

map(df -> mean(df[1]), groupby(iris, "Species"))

##############################################################################
#
# Reshaping
#
##############################################################################

#
# If you are looking for the equivalent of the R "Reshape" packages melt()
# and cast() functions, you can use stack() and unstack()
#
# Note that these functions have exactly the oppposite syntax as melt()
# and cast()
#

stack(iris, ["Petal.Length", "Petal.Width"])


text:
using Base.Broadcast: broadcastable
using BenchmarkTools
using ChainRulesCore
using LinearAlgebra
using LinearAlgebra.BLAS: ger!, gemv!, gemv, scal!
using StaticArrays
using SparseArrays
using Test

@testset "ChainRulesCore" begin
    @testset "differentials" begin
        include("differentials/abstract_zero.jl")
        include("differentials/thunks.jl")
        include("differentials/composite.jl")
        include("differentials/notimplemented.jl")
    end

    include("accumulation.jl")
    include("projection.jl")

    include("rules.jl")
    include("rule_definition_tools.jl")
    include("config.jl")

    include("deprecated.jl")
end


text:
function olMap(extrema::Array{<:Tuple}, colors::NamedTuple; port::Int = 3141)
    (xmin, xmax), (ymin, ymax) = extrema

    return js"""
        function (container){

            ////////////////////////////////////////////////// tile layer
            var projection = new $ol.proj.Projection({
                code: 'raster',

                extent: $([ymin,-xmax,ymax,-xmin]),
                units: 'pixels',
            })

            $ol.proj.addProjection(projection);
            var tiles = new $ol.layer.Tile({

                source: new $ol.source.XYZ({
                    tileUrlFunction: function (zxy) {

                        [z,x,y] = zxy
                        return 'http://localhost:$port/'+[z,y,x].join("/")+'.png?'+'channel='+document.getElementById('channel').value+'&scale='+document.getElementById('scale').value+'&seed='+Math.random()
                    },

                    projection: 'raster',
                    wrapX: false,
                })
            })

            ////////////////////////////////////////////////// polygon layer
            var polygons = new $ol.source.Vector({wrapX:false})
            var gates = new $ol.layer.Vector({
                source: polygons,

                style: new $ol.style.Style({
                    stroke: new $ol.style.Stroke({ color: 'rgb(0, 136, 170)', width: 3 }),
                }),
            })

            //////////////////////////////////////////////////////// initialise map
            var map = new $ol.Map({
                layers: [ tiles, gates ],
                target: 'map',

                view: new $ol.View({
                    projection: 'raster',
                    center: [0, 0],
                    zoom: 0,
                }),
            })

            var overlay = new $ol.Overlay({ element: document.getElementById('js-overlay') })
            map.addOverlay(overlay) // delete button pop-up

            document.getElementById("map").tiles = tiles.getSource()
            document.getElementById("map").gates = gates

            ///////////////////////////////////////////////////////// polygon interactions
            var selectedFeature
            var interaction
            var features

            document.querySelectorAll('#polygon-interactions button').forEach(
                button => button.addEventListener('click', function(event) {
                    var features

                    map.removeInteraction(interaction)
                    switch (event.target.id) {

                        case 'polygon':
                            interaction = new $ol.interaction.Draw({
                                type: 'Polygon', source: polygons })

                            map.addInteraction(interaction)
                            polygons.on('addfeature', function(event) {

                                console.assert( event.feature.getGeometry().getType()=="Polygon" )
                                event.feature.setId( uuidv4() )

                                ////////////////////////////////////////////// send gate to julia
                                fetch( 'http://localhost:$port/gate', {
                                    headers: { 'Content-Type': 'application/json' },
                                    method: 'POST',

                                    body: JSON.stringify({
                                        id: event.feature.getId(),
                                        coordinates: event.feature.getGeometry().getCoordinates()[0]
                                    })

                                //////////////////////////////////////////// update violins/stats
                                }).then( response => response.text() ).then( data => {
                                    var data = $d3.csvParse(data,$d3.autoType)

                                    violins({

                                        binCentres: Array.from( $d3.group(data,d=>d.channelValue).keys() ),
                                        density: Object.keys(data[0]).filter(x=>x!="channelValue").map( function(channel) {

                                            return { id: event.feature.getId(),
                                                name: channel, values: data.map(x=>x[channel]).map(
                                                y=>y/$d3.max(data.map(x=>x[channel])) ) }
                                        })
                                    })

                                }).catch( error => {
                                    console.error('Error:',error)
                                })
                            })
                            break

                        case 'modify':
                            interaction = new $ol.interaction.Modify({
                                features: new $ol.Collection(polygons.getFeatures()) })

                            map.addInteraction(interaction)
                            break

                        case 'delete':
                            interaction = new $ol.interaction.Select({
                                condition: $ol.events.condition.click, layers: [gates] })

                            map.addInteraction(interaction)
                            interaction.on('select', function(event) {

                                selectedFeature = event.selected[0]
                                selectedFeature ? overlay.setPosition(selectedFeature.getGeometry().getLastCoordinate()): overlay.setPosition(undefined)
                            })
                            break

                        default:
                            break
                    }

                    var snap = new $ol.interaction.Snap({source: polygons })
                    map.addInteraction(snap)
                })
            )

            document.getElementById('js-remove').addEventListener('click', function() {
                polygons.removeFeature(selectedFeature)

                var id = CSS.escape(selectedFeature.getId())
                document.querySelectorAll("#"+id).forEach( violin=> violin.remove())

                overlay.setPosition(undefined)
                interaction.getFeatures().clear()
            })

            var sidebar = new $ol.control.Sidebar({ element: 'sidebar', position: 'left' })
            map.addControl(sidebar)
        }
    """
end

max_stars_repo_name: gszep/flow-explorer

max_issues_repo_name: gszep/flow-explorer

max_forks_repo_name: gszep/flow-explorer


text:
# standalone generator
using Knet, JLD
include("model.jl")

function tokenize(vocabulary, prediction)
    index = findmax(prediction)[2]
    return filter(x->vocabulary[x] == index, collect(keys(vocabulary)))[1]
end

function gentoken(parameters, states, inputseq, vocabulary)
    result = Any[]
    hlayers = length(states) / 4
    hlayers = convert(Int, hlayers)

    # go through the encoder
    enstates = states[1:2hlayers]
    final_hidden = nothing
    for i=1:length(inputseq)
        input = oftype(parameters[1], inputseq[i])
        x = input * parameters[end-1]
        final_hidden = forward(parameters[1:2hlayers], enstates, x)
    end

    # initialize the hidden state of the decoder, decstates[1] is the initial hidden of the decoder lstm
    decstates = states[2hlayers+1:4hlayers]
    decstates[1] = final_hidden

    # token skeleton for decoder input, batchsize considered to be 1
    token = oftype(parameters[1], zeros(1, length(vocabulary)))

    input = copy(token)
    input[1] = oftype(input[1], 1.0) # first token is <s>
    it = tokenize(vocabulary, input)
    while (it != "</s>")
        x = input * parameters[end]
        hidden = forward(parameters[2hlayers+1:4hlayers], decstates, x)
        ypred = hidden * parameters[end-2] .+ parameters[end-3]
        ynorm = logp(ypred, 2)

        it = tokenize(vocabulary, ynorm)
        (it != "</s>") && println("predicted token is $it") # Don't want to see </s>
        input = copy(token)
        input[vocabulary[it]] = oftype(input[1], 1.0)
    end

end


text:
module AgentBasedModels

using DataFrames: AbstractAggregate, getiterator
using CUDA: findfirst, atomic_add!
using Base: add_with_overflow
using Random
using LinearAlgebra
using Distributions
using CUDA
using DataFrames
using CSV
import MacroTools: postwalk, @capture, inexpr, prettify, gensym_ids, flatten, unblock
#using WriteVTK

export Agent, @agent, add
export Model, compile
export Community, CommunityInTime, saveCSV, loadCommunityFromCSV, loadCommunityInTimeFromCSV
export SimulationFree, SimulationGrid
export SimulationSpace, FlatBoundary, Periodic, NonPeriodic, Bound
export MediumFlat, δMedium_
export configurator_
export extrude
export compactHexagonal

#Constants
include("./constants/constants.jl")
include("./constants/abstractTypes.jl")

#Agent
include("./agent/agentStructure.jl")
include("./agent/constructAgent.jl")

#Model
include("./model/structProgram.jl")

#Simulation Space
include("./simulationSpace/abstractTypes.jl")
include("./simulationSpace/medium/mediumStructures.jl")
include("./simulationSpace/simulationFree.jl")
include("./simulationSpace/simulationGrid.jl")

#Model
include("./model/model.jl")
include("./model/agentCode/addParameters.jl")
include("./model/agentCode/checkBounds.jl")
include("./model/agentCode/cleanInteraction.jl")
include("./model/agentCode/eventDeath.jl")
include("./model/agentCode/eventDivision.jl")
include("./model/agentCode/updateGlobal.jl")
include("./model/agentCode/updateLocal.jl")
include("./model/agentCode/updateLocalInteraction.jl")
include("./model/compile.jl")

#Community
include("./community/community.jl")
include("./community/baseModuleExtensions.jl")
include("./community/constructors/latices/hexagonal.jl")
include("./community/constructors/extrude.jl")
include("./community/IO/save.jl")
include("./community/IO/load.jl")

#Random
include("./model/random/distribution.jl")
include("./model/random/randomAdapt.jl")

#Integrators
include("./model/agentCode/integrator/euler.jl")
include("./model/agentCode/integrator/heun.jl")
include("./model/agentCode/integrator/integrators.jl")

#Integrators Medium
include("./model/agentCode/integratorMedium/ftcs.jl")
include("./model/agentCode/integratorMedium/lax.jl")
include("./model/agentCode/integratorMedium/leapfrog.jl")
include("./model/agentCode/integratorMedium/integratorsMedium.jl")
include("./model/agentCode/updateMediumInteraction.jl")

#Saving
include("./model/agentCode/saving/saveRAM.jl")
include("./model/agentCode/saving/saveCSV.jl")
include("./model/agentCode/saving/saving.jl")

#Auxiliar function
include("./auxiliar/checkDeclared.jl")
include("./auxiliar/clean.jl")
include("./auxiliar/substitution.jl")
include("./auxiliar/vectorize.jl")
include("./auxiliar/arguments.jl")
include("./auxiliar/wrapping.jl")
include("./auxiliar/emptyquote.jl")
include("./auxiliar/symbols.jl")
include("./auxiliar/updates.jl")
include("./auxiliar/substitution2.jl")
include("./auxiliar/extract.jl")

#Cuda specific functions
include("./cuda/cudaAdapt.jl")
include("./cuda/cudaConfigurator.jl")

end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia