text:
# ## `projective/bounds.jl`
#
# To make it easier to handle spatial bounds, we define
# some helpers here.
#
# [`makebounds`](#) creates bounds from side lengths or ranges.

"""
    makebounds(sz[, T])
    makebounds(ranges[, T])

Helper for creating spatial bounds.

## Examples

{cell=makebounds}
```julia
using DataAugmentation: makebounds, showbounds
makebounds((100, 100), Float32)
```
{cell=makebounds}
```julia
makebounds((100, 100)) == makebounds((1:100, 1:100))
```

{cell=makebounds}
```julia
bounds = makebounds((100, 100))
showbounds(bounds)
```
"""
function makebounds(sz::NTuple{N, Int}, T = Float32) where N
    return makebounds(Tuple(1:a for a in sz), T)
end

function makebounds(ranges::NTuple{N, R}, T = Float32) where {N, R<:AbstractUnitRange}
    return collect(SVector{N, T}, Iterators.product(((r[begin]-1, r[end]) for r in ranges)...))
end

# [`boundsextrema`](#), [`boundsranges`](#) and [`boundssize`](#) represent
# the inverse of `makebounds`.

function boundsextrema(bounds::AbstractArray{<:SVector{N}}) where N
    mins = Tuple(floor(Int, minimum(getindex.(bounds, i))) for i = 1:N)
    maxs = Tuple(ceil(Int, maximum(getindex.(bounds, i))) for i = 1:N)
    return mins, maxs
end

function boundsranges(bounds)
    mins, maxs = boundsextrema(bounds)
    return UnitRange.(mins .+ 1, maxs)
end

"""
    boundssize(bounds)

`(100, 100) |> makebounds |> boundssize == (100, 100)`
"""
boundssize(bounds) = length.(boundsranges(bounds))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct MachineConfiguration <: SwaggerModel
    vcpu_count::Union{ Nothing, Int32 } # vcpu_count
    mem_size_mib::Union{ Nothing, Int32 } # mem_size_mib
    ht_enabled::Union{ Nothing, Bool } # ht_enabled
    cpu_template::Union{ Nothing, CpuTemplate } # cpu_template

    function MachineConfiguration(;vcpu_count=nothing, mem_size_mib=nothing, ht_enabled=nothing, cpu_template=nothing)
        o = new()
        set_field!(o, :vcpu_count, vcpu_count)
        set_field!(o, :mem_size_mib, mem_size_mib)
        set_field!(o, :ht_enabled, ht_enabled)
        set_field!(o, :cpu_template, cpu_template)
        o
    end
end # type MachineConfiguration

const _name_map_MachineConfiguration = Dict{String,Symbol}(["vcpu_count"=>:vcpu_count, "mem_size_mib"=>:mem_size_mib, "ht_enabled"=>:ht_enabled, "cpu_template"=>:cpu_template])
const _field_map_MachineConfiguration = Dict{Symbol,String}([:vcpu_count=>"vcpu_count", :mem_size_mib=>"mem_size_mib", :ht_enabled=>"ht_enabled", :cpu_template=>"cpu_template"])
Swagger.name_map(::Type{ MachineConfiguration }) = _name_map_MachineConfiguration
Swagger.field_map(::Type{ MachineConfiguration }) = _field_map_MachineConfiguration

function check_required(o::MachineConfiguration)
    true
end

function validate_field(o::MachineConfiguration, name::Symbol, val)
    if name === :vcpu_count
        Swagger.validate_param(name, "MachineConfiguration", :maximum, val, 32, false)
        Swagger.validate_param(name, "MachineConfiguration", :minimum, val, 1, false)
    end
end

ext: jl


text:
#export DistPred

## type DGlmResp                    # distributed response in a glm model
##     d::Distribution
##     l::Link
##     eta::DArray{Float64,1,1}     # linear predictor
##     mu::DArray{Float64,1,1}      # mean response
##     offset::DArray{Float64,1,1}  # offset added to linear predictor (usually 0)
##     wts::DArray{Float64,1,1}     # prior weights
##     y::DArray{Float64,1,1}       # response
##     ## FIXME: Add compatibility checks here
## end

## function DGlmResp(d::Distribution, l::Link, y::DArray{Float64,1,1})
##     wt     = darray((T,d,da)->ones(T,d), Float64, size(y), distdim(y), y.pmap)
##     offset = darray((T,d,da)->zeros(T,d), Float64, size(y), distdim(y), y.pmap)
##     mu     = similar(y)
##     @sync begin
##         for p = y.pmap
##             @spawnat p copy_to(localize(mu), d.mustart(localize(y), localize(wt)))
##         end
##     end
##     dGlmResp(d, l, map_vectorized(link.linkFun, mu), mu, offset, wt, y)
## end

## DGlmResp(d::Distribution, y::DArray{Float64,1,1}) = DGlmResp(d, canonicallink(d), y)

## At_mul_B{T <: Real}(A::DArray{T, 2, 1}, B::DArray{T, 2, 1}) = Ac_mul_B(A, B)

## function Ac_mul_B{T <: Real}(A::DArray{T, 2, 1}, B::DArray{T, 2, 1})
##     if (all(procs(A) != procs(B)) || all(dist(A) != dist(B)))
##         # FIXME: B should be redistributed to match A
##         error("Arrays A and B must be distributed similarly")
##     end
##     if is(A, B)
##         return mapreduce(+, fetch, {@spawnat p BLAS.syrk('T', localize(A)) for p in procs(A)})
##     end
##     mapreduce(+, fetch, {@spawnat p Ac_mul_B(localize(A), localize(B)) for p in procs(A)})
## end

## function Ac_mul_B{T <: Real}(A::DArray{T, 2, 1}, B::DArray{T, 1, 1})
##     if (all(procs(A) != procs(B)) || all(dist(A) != dist(B)))
##         # FIXME: B should be redistributed to match A
##         error("Arrays A and B must be distributed similarly")
##     end
##     mapreduce(+, fetch, {@spawnat p Ac_mul_B(localize(A), localize(B)) for p in procs(A)})
## end

## type DistPred{T} <: LinPred   # predictor with distributed (on rows) X
##     X::DArray{T, 2, 1}        # model matrix
##     beta::Vector{T}           # coefficient vector
##     r::CholeskyDense{T}
##     function DistPred(X, beta)
##         if size(X, 2) != length(beta) error("dimension mismatch") end
##         new(X, beta, chold(X'X))
##     end
## end

## function (\)(A::DArray{Float64,2,1}, B::DArray{Float64,1,1})
##     R   = Cholesky(A'A)
##     LAPACK.potrs!('U', R, A'B)
## end


text:
"""
### function plot_lsim(y, t)

Plot the result computed by `lsim` function.

##### Args

* y: Output vector.
* t: Time vector.

##### Parameters

* title: Title of the figure.

"""

function plot_lsim(y, t; title="")
    # Get the number of outputs in y.
    num_outputs = size(y[1])[1]

    # Compute the grid of the subplots.
    nx = (num_outputs > 3) ? 3 : num_outputs
    ny = (num_outputs > 3) ? num_outputs%3 : 1

    if (plot_system == :pyplot)
        fig = plt.figure(title)

        ax = Array{Any}(num_outputs)

        for i = 1:num_outputs
            ax[i] = fig[:add_subplot](nx, ny, i)
        end

        for i = 1:num_outputs
            ax[i][:plot](t, map(x->x[i], y), linewidth=plw)
            ax[i][:set_title]("Y$i")
            ax[i][:set_xlabel]("Time")
            ax[i][:set_ylabel]("Output")

            # Change figure border width.
            ax[i][:spines]["top"][:set_linewidth](blw)
            ax[i][:spines]["bottom"][:set_linewidth](blw)
            ax[i][:spines]["left"][:set_linewidth](blw)
            ax[i][:spines]["right"][:set_linewidth](blw)
        end

        fig[:tight_layout]()
    end
end

"""
### function plot_lsim(y, t, u::Function)

Plot the result computed by `lsim` function together with the input function
`u`.

##### Args

* y: Output vector.
* t: Time vector.
* u: Input function.

##### Parameters

* title: Title of the figure.

##### Remarks

It is assumed that `u` function always returns a scalar.

"""

function plot_lsim(y, t, u::Function; title="")
    # Get the number of outputs in y.
    num_outputs = size(y[1])[1]

    # Compute the grid of the subplots.
    nx = (num_outputs > 3) ? 3 : num_outputs
    ny = (num_outputs > 3) ? num_outputs%3 : 1

    if (plot_system == :pyplot)
        fig = plt.figure(title)

        ax = Array{Any}(num_outputs)

        for i = 1:num_outputs
            ax[i] = fig[:add_subplot](nx, ny, i)
        end

        for i = 1:num_outputs
            # Plot the input function.
            ax[i][:plot](t, map(t->u(t), t),
                         color="#BBBBBB", linewidth=plw, linestyle="--")

            ax[i][:plot](t, map(x->x[i], y), linewidth=plw)
            ax[i][:set_title]("Y$i")
            ax[i][:set_xlabel]("Time")
            ax[i][:set_ylabel]("Output")

            # Change figure border width.
            ax[i][:spines]["top"][:set_linewidth](blw)
            ax[i][:spines]["bottom"][:set_linewidth](blw)
            ax[i][:spines]["left"][:set_linewidth](blw)
            ax[i][:spines]["right"][:set_linewidth](blw)
        end

        fig[:tight_layout]()
    end
end


text:
# utilities

_has_group(df::DataFrames.AbstractDataFrame, group::Any) = false
_has_group(df::DataFrames.AbstractDataFrame, group::Symbol) = hasproperty(df, group)
function _has_group(df::DataFrames.AbstractDataFrame, group::Vector{Symbol})
    all(x -> hasproperty(df, x), group)
end

_group_name(df::DataFrames.AbstractDataFrame, group::Symbol) = df[1, group]
function _group_name(df::DataFrames.AbstractDataFrame, groups::Vector{Symbol})
    string("(", join([df[1, g] for g in groups], ", "), ")")
end

"""
$(SIGNATURES)
Build a trace of kind `kind`, using the columns of `df` where possible. In
particular for all keyword arguments, if the value of the keyword argument is a
Symbol and matches one of the column names of `df`, replace the value of the
keyword argument with the column of `df`

If `group` is passed and is a Symbol that is one of the column names of `df`,
then call `by(df, group)` and construct one trace per SubDataFrame, passing
all other keyword arguments. This means all keyword arguments are passed
applied to all traces

Also, when using this routine you can pass a function as a value for any
keyword argument. This function will be replaced by calling the function on the
DataFrame. For example, if I were to pass `name=(df) -> "Wage (average =
\$(mean(df[!, :X1])))"` then the `name` attribute on the trace would be replaced by
the  `Wage (average = XX)`, where `XX` is the average of the `X1` column in the
DataFrame.

The ability to pass functions as values for keyword arguments is particularly
useful when using the `group` keyword arugment, as the function will be applied
to each SubDataFrame. In the example above, the name attribute would set a
different mean for each group.
"""
function GenericTrace(df::DataFrames.AbstractDataFrame; group=nothing, kind="scatter", kwargs...)
    d = Dict{Symbol,Any}(kwargs)
    if _has_group(df, group)
        _traces = []
        for dfg in DataFrames.groupby(df, group)
            push!(_traces,  GenericTrace(dfg; kind=kind, name=_group_name(dfg, group), kwargs...))
        end
        return GenericTrace[t for t in _traces]
    else
        if (group !== nothing)
            @warn "Unknown group $(group), skipping"
        end
    end

    for (k, v) in d
        if isa(v, Symbol) && hasproperty(df, v)
            d[k] = df[!, v]
        elseif isa(v, Function)
            d[k] = v(df)
        end
    end
    GenericTrace(kind; d...)
end

"""
$(SIGNATURES)
Pass the provided values of `x` and `y` as keyword arguments for constructing
the trace from `df`. See other method for more information
"""
function GenericTrace(df::DataFrames.AbstractDataFrame, x::Symbol, y::Symbol; kwargs...)
    GenericTrace(df; x=x, y=y, kwargs...)
end

"""
$(SIGNATURES)
Pass the provided value `y` as keyword argument for constructing the trace from
`df`. See other method for more information
"""
function GenericTrace(df::DataFrames.AbstractDataFrame, y::Symbol; kwargs...)
    GenericTrace(df; y=y, kwargs...)
end

"""
$(SIGNATURES)
Construct a plot using the columns of `df` if possible. For each keyword
argument, if the value of the argument is a Symbol and the `df` has a column
whose name matches the value, replace the value with the column of the `df`.

If `group` is passed and is a Symbol that is one of the column names of `df`,
then call `by(df, group)` and construct one trace per SubDataFrame, passing
all other keyword arguments. This means all keyword arguments are passed
applied to all traces
"""
function Plot(df::DataFrames.AbstractDataFrame, l::Layout=Layout();
              style::Style=CURRENT_STYLE[], kwargs...)
    Plot(GenericTrace(df; kwargs...), l, style=style)
end

"""
$(SIGNATURES)
Construct a plot from `df`, passing the provided values of x and y as keyword
arguments. See docstring for other method for more information.
"""
function Plot(d::DataFrames.AbstractDataFrame, x::Symbol, y::Symbol, l::Layout=Layout();
              style::Style=CURRENT_STYLE[], kwargs...)
    Plot(d, l; x=x, y=y, style=style, kwargs...)
end

"""
$(SIGNATURES)
Construct a plot from `df`, passing the provided value y as a keyword argument.
See docstring for other method for more information.
"""
function Plot(d::DataFrames.AbstractDataFrame, y::Symbol, l::Layout=Layout();
              style::Style=CURRENT_STYLE[], kwargs...)
    Plot(d, l; y=y, style=style, kwargs...)
end

for t in _TRACE_TYPES
    str_t = string(t)
    @eval $t(df::DataFrames.AbstractDataFrame; kwargs...) = GenericTrace(df; kind=$(str_t), kwargs...)
end


text: include("blapi.jl")


text:
Mesh(::MetaGraph)::Mesh
Mesh(::MetaDiGraph)::Mesh

MetaGraph(::Mesh)::MetaGraph
MetaDiGraph(::Mesh)::MetaDiGraph


text:
struct EmbeddedMap{T<:Integer}
    P::PlanarMap{T}
    locs::Array{<:Tuple{<:Real,<:Real}}
end

function EmbeddedMap(RNG::Random.AbstractRNG,
                     P::PlanarMap;
                     outeredge=(1,neighbors(P,1)[1]),
                     corners=nothing)
    T = triangulation(RNG,P,outeredge=outeredge,corners=corners)
    outface = corners ≠ nothing ? Face(corners...) : face(T,outeredge...)
    WT = schnyderwood(RNG,T,outface=outface)
    n = length(P)
    map((a,b)->ϕ(a,b,n),EmbeddedMap(P,schnyder_coords(WT)))
end

EmbeddedMap(P::PlanarMap;kwargs...) =
        EmbeddedMap(Random.MersenneTwister(0),P;kwargs...)

length(E::EmbeddedMap) = length(E.P)

import Base.map
function map(f::Function,E::EmbeddedMap)
    EmbeddedMap(E.P,[f(a,b) for (a,b) in E.locs])
end

ϕ(a::Real,b::Real,n::Real) = tuple(([1 1/2; 0 sqrt(3)/2]*[a; b]/n)...)
ϕ(a::Real,b::Real) = ϕ(a,b,1)
ϕ(t::NTuple{2}) = ϕ(t...)

function bbox_size(E::EmbeddedMap)
    My = maximum(y for (x,y) in E.locs)
    my = minimum(y for (x,y) in E.locs)
    Mx = maximum(x for (x,y) in E.locs)
    mx = minimum(x for (x,y) in E.locs)
    max(Mx-mx,My-my)
end

struct GraphPenSet
    pointpen
    diskpen
    labelpen
    edgepen
    facepen
    pointpens
    diskpens
    labelpens
    edgepens
    facepens
end

function GraphPenSet(;kwargs...)
    args = DataStructures.OrderedDict(
            :pointpen=>AsyPlots.Pen(color="DarkGreen"),
            :diskpen=>AsyPlots.Pen(color="Black"),
            :labelpen=>AsyPlots.Pen(fontsize=12,color="White"),
            :edgepen=>AsyPlots.Pen(color="Navy",linewidth=1.2),
            :facepen=>AsyPlots.Pen(color="LightBlue"),
            :pointpens=>nothing,
            :diskpens=>nothing,
            :labelpens=>nothing,
            :edgepens=>nothing,
            :facepens=>nothing)

    merge!(args,Dict(kwargs))
    GraphPenSet(values(args)...)
end

for name in (:point,:disk,:label,:edge,:face)
    namepen = Symbol(string(name)*"pen")
    namepens = Symbol(string(name)*"pens")
    @eval function $namepen(G::GraphPenSet,k)
        if hasmethod(getindex,map(typeof,(G.$namepens,k)))
            G.$namepens[k]
        elseif !isempty(methods(G.$namepens))
            G.$namepens(k)
        else
            G.$namepen
        end
    end
end

function draw(E::EmbeddedMap;
              outeredge=(1,neighbors(E.P,1)[1]),
              drawpoints=true,
              drawedges=true,
              fillfaces=false,
              drawlabels=true,
              graphpenset=GraphPenSet(),
              kwargs...)

    if length(kwargs) > 0
        G = GraphPenSet(;kwargs...)
    else
        G = graphpenset
    end
    s = bbox_size(E) # for sizing the disk that the label goes in
    grlist = AsyPlots.GraphicElement2D[]
    if fillfaces
        append!(grlist,[AsyPlots.Polygon2D([E.locs[i] for i in fc],
                            pen=AsyPlots.NoPen(),
                            fillpen=facepen(G,fc)) for fc in interiorfaces(
                                faces(E.P;outeredge=outeredge))])
    end
    if drawedges
        append!(grlist,[AsyPlots.Path([E.locs[i],E.locs[j]],pen=edgepen(G,(i,j)))
                                                    for (i,j) in edges(E.P)])
    end
    if drawlabels
        append!(grlist,vcat([[AsyPlots.Circle(E.locs[i],
                                              s*labelpen(G,i).fontsize/500,
                                     pen=diskpen(G,i),
                                     fillpen=pointpen(G,i)),
                              AsyPlots.Label(string(i),E.locs[i];
                                     pen=labelpen(G,i))] for i in 1:length(E.P)]...))
    elseif drawpoints
        append!(grlist,[AsyPlots.Point(E.locs[i],
                                  pen=pointpen(G,i)) for i in 1:length(E.P)])
    end
    AsyPlots.Plot(grlist)
end

"""
    draw(P,outeredge,corners)

Draw the planar map P with outer face either:
    - the face to the left of the `outeredge`
    - the face containing the sequence of vertices `corners`

- `outeredge` is an ordered pair
- `corners` is either `nothing` (to indicate deference to
    specification via `outeredge`) or an ordered triple
"""
function draw(P::PlanarMap;
              outeredge=(1,neighbors(P,1)[1]),
              corners=nothing,
              kwargs...)
    draw(EmbeddedMap(P;outeredge=outeredge,corners=corners);kwargs...)
end

function draw(WT::WoodedTriangulation;
              rotation=0,
              linewidth=10/length(WT),
              pointsize=3,
              pointcolor=AsyPlots.NamedColor(0.65,0.65,0.65),
              linecolor="default",
              labelcolor="black",
              labelsize=max(1,floor(log10(length(WT))))*min(12,max(144/length(WT))),
              kwargs...)
    ASYCOLORS =     Dict(zip((:blue, :red, :green,     :gray),
        AsyPlots.NamedColor.(("Navy","Red","DarkGreen","Gray"))))
    edgedict = Dict((e.tail,e.head)=>e.color for e in edges(WT))
    edgepens(t) = linecolor == "default" ?
                        AsyPlots.Pen(linewidth=linewidth,
                                     color=ASYCOLORS[edgedict[t]]) :
                                AsyPlots.Pen(linewidth=linewidth,
                                             color=linecolor)
    function facepens(F::Face)
        AsyPlots.Pen(color=sum(1/3 * AsyPlots.NamedColor(
                            string(edgedict[e])) for e in edges(F)))
    end
    GPS = GraphPenSet(;edgepens=edgepens,
                       facepens=facepens,
                       diskpen=AsyPlots.Pen(linewidth=10/length(WT)),
                       pointpen=AsyPlots.Pen(linewidth=pointsize,
                                             color=pointcolor),
                       labelpen=AsyPlots.Pen(fontsize=labelsize,color=labelcolor))
    draw(EmbeddedMap(WT.P,
         map((t->reim(cis(rotation)*(t[1]+im*t[2])))∘ϕ,schnyder_coords(WT))),
         outeredge=(blueroot(WT),greenroot(WT));
         graphpenset=GPS,kwargs...)
end


text:
@testset "polyagammanegativemultinomial" begin
    # TODO add more basic tests
end


text:
execute = isempty(ARGS) || ARGS[1] == "run"

#using Pkg
#Pkg.add(url="https://github.com/StevenWhitaker/PERK.jl")

using PERK
using Documenter
using Literate

# https://juliadocs.github.io/Documenter.jl/stable/man/syntax/#@example-block
ENV["GKSwstype"] = "100"
ENV["GKS_ENCODING"] = "utf-8"

# generate examples using Literate
lit = joinpath(@__DIR__, "lit")
src = joinpath(@__DIR__, "src")
gen = joinpath(@__DIR__, "src/generated")

base = "StevenWhitaker/PERK.jl"
repo_root_url =
    "https://github.com/$base/blob/main/docs/lit/examples"
nbviewer_root_url =
    "https://nbviewer.org/github/$base/tree/gh-pages/dev/generated/examples"
binder_root_url =
    "https://mybinder.org/v2/gh/$base/gh-pages?filepath=dev/generated/examples"

DocMeta.setdocmeta!(PERK, :DocTestSetup, :(using PERK); recursive=true)

for (root, _, files) in walkdir(lit), file in files
    splitext(file)[2] == ".jl" || continue # process .jl files only
    ipath = joinpath(root, file)
    opath = splitdir(replace(ipath, lit => gen))[1]
    Literate.markdown(ipath, opath, documenter = execute; # run examples
        repo_root_url, nbviewer_root_url, binder_root_url)
    Literate.notebook(ipath, opath; execute = false, # no-run notebooks
        repo_root_url, nbviewer_root_url, binder_root_url)
end

# Documentation structure
ismd(f) = splitext(f)[2] == ".md"
pages(folder) =
    [joinpath("generated/", folder, f) for f in readdir(joinpath(gen, folder)) if ismd(f)]

isci = get(ENV, "CI", nothing) == "true"

format = Documenter.HTML(;
    prettyurls = isci,
    edit_link = "main",
    canonical = "https://StevenWhitaker.github.io/PERK.jl/stable",
#   assets = String[],
)

makedocs(;
    modules = [PERK],
    authors = "Steven Whitaker and contributors",
    sitename = "PERK.jl",
    format,
    pages = [
        "Home" => "index.md",
        "Methods" => "methods.md",
        "Examples" => pages("examples")
    ],
)

if isci
    deploydocs(;
        repo = "github.com/StevenWhitaker/PERK.jl.git",
        devbranch = "main",
        devurl = "dev",
        versions = ["stable" => "v^", "dev" => "dev"],
        forcepush = true,
#       push_preview = true,
    )
end


text:
# This file holds the definition of the functions pertaining to the Transition
#   type.

export  setTransitionName!,
        setTransitionNode!,
        setTransitionIsOut!,
        setTransitionSchedule!,
        setTransitionMaxAttempts!,
        setTransitionFluxLimits!,
        setTransitionHasPriority!,
        addTransitionCondition!,
        clearTransitionConditions!,
        setTransitionConditions!,
        addTransitionAttributeChange!,
        clearTransitionAttributeChanges!,
        setTransitionAttributeChanges!,
        setTransitionProbabilities!

"""
```
setTransitionName!(
    transition::Transition,
    name::String )
```
This function sets the name of the transition `transition` to `name`.

This function returns `true`, indication that the transition's name has been successfully set.
"""
function setTransitionName!( transition::Transition, name::String )::Bool

    transition.name = name
    return true

end  # setTransitionName!( trans, name )

"""
```
setTransitionNode!(
    transition::Transition,
    node::String,
    isTargetNode::Bool = false )
```
This function sets one of the nodes (source/target) of the transition `transitoin` to `node`. If `istargetNode` is `true`, the target node will be set, otherwise the source state will be set. Note that the target node will be ignored if the transition's `isOutTransition` flag is set to `true`.

This function returns `true`, indicating that the node has been successfully set.
"""
function setTransitionNode!( transition::Transition, node::String,
    isTargetNode::Bool = false )::Bool

    tmpNode = lowercase( node ) == "dummy" ? "dummy" : node

    if isTargetNode
        transition.targetNode = tmpNode
    else
        transition.sourceNode = tmpNode
    end  # if isTargetNode

    return true

end  # setTransitionNode!( transition, node, isTargetNode )

"""
```
setTransitionIsOut!(
    transition::Transition,
    isOutTransition::Bool )
```
This function sets the `isOutTransition` flag of the transition `transition` to `isOutTransition`.

This function returns `true`, indicating the flag has been successfully set.
"""
function setTransitionIsOut!( transition::Transition,
    isOutTransition::Bool )::Bool

    transition.isOutTransition = isOutTransition
    return true

end  # setTransitionIsOut!( transition, isOutTransition )

"""
```
setTransitionSchedule!(
    transition::Transition,
    freq::Real,
    offset::Real = 0.0 )
```
This function sets the schedule of the transition `transition`. This transition will be checked every `freq` time units with an offset of `offset` with respect to the start of the simulation.

If the entered period is ⩽ 0, the function issues a warning and will make no changes.

This function returns `true` if the schedule has been successfully set, and `false` if the entered period was negative.
"""
function setTransitionSchedule!( transition::Transition, freq::Real,
    offset::Real = 0.0 )::Bool

    if freq <= 0.0
        @warn "Time between two transition checks must be > 0.0, not making any changes."
        return false
    end  # if freq <= 0.0

    transition.freq = freq
    transition.offset = offset % freq + ( offset < 0.0 ? freq : 0.0 )
    return true

end  # setTransitionSchedule!( transition, freq, offset )

"""
```
setTransitionMaxAttempts!(
    transition::Transition,
    maxAttempts::Integer )
```
This function sets the maximum number of times a personnel member may attempt the transition `transition` to `maxAttempts`.

A value < 0 for `maxAttempts`  means that there is no restriction on number of attempts, and 0 means there are as many attempts as entries in the vector of probabilities.

This function returns `true`, indicating that the maximum number of attempts at the transition is successfully set.
"""
function setTransitionMaxAttempts!( transition::Transition,
    maxAttempts::Integer )::Bool

    transition.maxAttempts = max( maxAttempts, -1 )
    return true

end  # setTransitionMaxAttempts!( transition, maxAttempts )

"""
```
setTransitionFluxLimits!(
    transition::Transition,
    minFlux::Integer,
    maxFlux::Integer )
```
This function sets the limits of the flux of the transition `transition` to the range `minFlux` to `maxFlux`.

A negative value for the maximum flux means that there is no upper bound on the flux.

If the minimum flux is larger than the maximum flux, the function gives a warning and doesn't make any changes.

This function returns `true` if the limits on the transition's flux have been successfully set, and `false` if the entered minimum flux was larger than the maximum flux.
"""
function setTransitionFluxLimits!( transition::Transition, minFlux::Integer, maxFlux::Integer )::Bool

    if ( maxFlux >= 0 ) && ( minFlux > maxFlux )
        @warn "The minimum flux of the transition cannot be higher than the maximum flux, not making any changes."
        return false
    end  # if ( maxFlux >= 0 ) && ...

    transition.minFlux = max( minFlux, 0 )
    transition.maxFlux = max( maxFlux, -1 )
    return true

end  # setTransitionFluxLimits!( transition, minFlux, maxFlux )

"""
```
setTransitionHasPriority!(
    transition::Transition,
    hasPriority::Bool )
```
This function sets the `hasPriority` flag of the transition `transition` to `hasPriority`, where `true` means the transition can overrule the target state's population target.

This function returns `true`, indicating that the flas is successfully set.
"""
function setTransitionHasPriority!( transition::Transition,
    hasPriority::Bool )::Bool

    transition.hasPriority = hasPriority

    if hasPriority && ( transition.priority > 0 )
        transition.priority *= -1
    elseif !hasPriority && ( transition.priority < 0 )
        transition.priority *= -1
    end  # if hasPriority && ( transition.priority > 0 )

    return true

end  # setTransitionHasPriority!( transition, hasPriority )

"""
```
addTransitionCondition!(
    transition::Transition,
    conditions::MPcondition... )
```
This function adds the conditions in  `conditions` as extra conditions to the transition `transition`. This function does NOT check if these conditions are contradictory with each other, with the existing ones, or with the source node.

This function returns `true`, indicating that the condition has been successfully added.
"""
function addTransitionCondition!( transition::Transition,
    conditions::MPcondition... )::Bool

    append!( transition.extraConditions, collect( conditions ) )
    return true

end  # addTransitionCondition!( transition, conditions )

"""
```
clearTransitionConditions!( transition::Transition )
```
This function clear all extra conditions from the transition `transition`.

This function returns `true`, indicating that the list of conditions has been successfully cleared.
"""
function clearTransitionConditions!( transition::Transition )::Bool

    empty!( transition.extraConditions )
    return true

end  # clearTransitionConditions!( transition )

"""
```
setTransitionConditions!(
    transition::Transition,
    conditions::Vector{MPcondition} )
```
This function sets the extra conditions for the transition `transition` to the list of conditions `conditions`. This function does NOT check if these conditions are contradictory with each other or with the source node.
"""
function setTransitionConditions!( transition::Transition,
    conditions::Vector{MPcondition} )::Bool

    transition.extraConditions = conditions
    return true

end  # setTransitionConditions!( transition, conditions )

"""
```
addTransitionAttributeChange!(
    transition::Transition,
    attrVals::NTuple{2,String}... )
```
This function adds the extra attribute changes in the list `attrVals` to the transition `transition`. These changes take priority over the changes imposed by moving from the source to the target node.

If the list contains multiple changes for the attribute, the function issues a warning and doesn't make any changes.

This function returns `true` if the attribute changes have been successfully added, and `false` if there were duplicate entries in the list.
"""
function addTransitionAttributeChange!( transition::Transition,
    attrVals::NTuple{2,String}... )::Bool

    newAttrs = map( attrVal -> attrVal[1], collect( attrVals ) )

    if length( newAttrs ) != length( unique( newAttrs ) )
        @warn "Duplicate entries in the attribute/value list, not making any changes."
        return false
    end  # if length( newAttrs ) != length( unique( newAttrs ) )

    for attrVal in attrVals
        transition.extraChanges[attrVal[1]] = attrVal[2]
    end  # for attrVal in attrVals

    return true

end  # addTransitionAttributeChange!( transition, attrVals )

"""
```
clearTransitionAttributeChanges!( transition::Transition )
```
This function clears all attribute changes from the transition `transition`.

This function returns `true`, indicating that the list of attribute changes has been successfully cleared.
"""
function clearTransitionAttributeChanges!( transition::Transition )::Bool

    empty!( transition.extraChanges )
    return true

end  # clearTransitionAttributeChanges!( trans )

"""
```
setTransitionAttributeChanges!(
    transition::Transition,
    attrVals::Dict{String,String} )
```
This function sets the additional attribute changes made by the transition `transition` to the list in `attrVals`. These changes take priority over the changes imposed by moving from the source to the target node.

This function returns `true`, indicating the list of changes has been successfully set.
"""
function setTransitionAttributeChanges!( transition::Transition,
    attrVals::Dict{String,String} )::Bool

    transition.extraChanges = deepcopy( attrVals )
    return true

end  # setTransitionAttributeChanges!( transition, attrVals )

"""
```
setTransitionAttributeChanges!(
    transition::Transition,
    attrVals::NTuple{2,String}... )
```
This function sets the additional attribute changes made by the transition `transition` to the list in `attrVals`. These changes take priority over the changes imposed by moving from the source to the target node.

If the list contains multiple entries for the same attribute, the function issues a warning and doesn't make any changes.

This function returns `true` if the list of changes has been successfully set, and `false` if that list contained duplicate entries.

"""
function setTransitionAttributeChanges!( transition::Transition, attrVals::NTuple{2,String}... )::Bool

    # Check for duplicates.
    newAttrs = map( attrVal -> attrVal[1], collect( attrVals ) )

    if length( newAttrs ) != length( unique( newAttrs ) )
        @warn "Duplicate entries in the attribute/value list, not making any changes."
        return false
    end  # if length( newAttrs ) != length( unique( newAttrs ) )

    # Convert the val/weight pairs to a dictionary.
    attrValDict = Dict{String,String}()

    for attrVal in attrVals
        attrValDict[attrVal[1]] = attrVal[2]
    end  # for attrVal in attrVals

    return setTransitionAttributeChanges!( transition, attrValDict )

end  # setTransitionAttributeChanges!( transition, attrVals )

"""
```
setTransitionAttributeChanges!(
    transition::Transition,
    attributes::Vector{String},
    values::Vector{String} )
```
This function sets the additional attribute changes made by the transition `transition` to the attributes in `attributes` and the values in `values`. These changes take priority over the changes imposed by moving from the source to the target node.

If the list contains multiple entries for the same attribute, the function issues a warning and doesn't make any changes.

If there is a mismatch between the lengths of the list of attributes and the list of values, the function issues a warning and doesn't make any changes.

This function returns `true` if the list of changes has been successfully set, and `false` if that list contained duplicate entries.
"""
function setTransitionAttributeChanges!( transition::Transition,
    attributes::Vector{String}, values::Vector{String} )::Bool

    if length( attributes ) != length( values )
        @warn "Mismatched lengths of vector of attributes and values, not making any changes."
        return false
    end  # if length( attributes ) != length( values )

    if length( attributes ) != length( unique( attributes ) )
        @warn "Duplicate entries in the attributes list, not making any changes."
        return false
    end  # if length( attributes ) != length( unique( attributes ) )

    attrValDict = Dict{String,String}()

    for ii in eachindex( attributes )
        attrValDict[attributes[ii]] = values[ii]
    end  # for ii in eachindex( attributes )

    return setTransitionAttributeChanges!( transition, attrValDict )

end  # setTransitionAttributeChanges!( transition, attributes, values )

"""
```
setTransitionProbabilities!(
    transition::Transition,
    probabilities::Vector{Float64} )
```
This function sets the execution probabilities of the state transition `trans`
to the valid entries in `probs`. Valid entries are between 0.0 and 1.0
inclusive, and invalid entries are removed. Additionally, all leading 0.0
entries will be removed as well.

This function returns `nothing`. If there are no valid probabilities, the
function will issue a warning to that effect, and won't make any changes.
"""
function setTransitionProbabilities!( transition::Transition, probabilities::Vector{Float64} )::Bool

    # Filter out all invalid probabilities.
    tmpProbs = filter( prob -> 0.0 <= prob <= 1.0, probabilities )

    if isempty( tmpProbs )
        @warn "No valid probabilities entered, not making any changes."
        return false
    end  # if isempty( tmpProbs )

    transition.probabilityList = tmpProbs
    return true

end  # setTransitionProbabilities!( trans, probs )

function Base.show( io::IO, transition::Transition )

    print( io, "Transition '", transition.name, "': '",
        transition.sourceNode, "' to '",
        transition.isOutTransition ? "external" : transition.targetNode, "'" )
    print( io, "\n  Occurs with period ", transition.freq, " (offset ",
        transition.offset, ")" )

    if !isempty( transition.extraConditions )
        print( io, "\n  Extra conditions: ", transition.extraConditions )
    end  # if !isempty( transition.extraConditions )

    if !isempty( transition.extraChanges )
        print( io, "\n  Extra changes: ",
            join( map( attr -> string( attr, " to ",
            transition.extraChanges[attr] ),
            collect( keys( transition.extraChanges ) ) ), ", " ) )
    end  # if !isempty( transition.extraChanges )

    if transition.maxAttempts == -1
        print( io, "\n  Infinite number of" )
    else
        nAttempts = transition.maxAttempts == 0 ?
            length( transition.probabilityList ) : transition.maxAttempts
        print( io, "\n  Max ", nAttempts )
    end  # if transition.maxAttempts == -1

    print( io, " attempts" )

    if transition.minFlux > 0
        print( io, "\n  Min flux: ", transition.minFlux )
    end  # if transition.minFlux > 0

    if transition.maxFlux >= 0
        print( io, "\n  Max flux: ", transition.maxFlux )

        if transition.hasPriority
            print( io, " (overrides state pop. constraint)" )
        end  # if transition.hasPriority
    end  # if transition.maxFlux >= 0

    print( io, "\n  Transition execution probabilities: ",
        join( transition.probabilityList .* 100, "%, " ), '%' )
    print( io, "\n Transition at priority ", - transition.priority )

end  # Base.show( io, transition )


text:
module Ctrl

# Types
import RIOPA: Config, DataSet, DataSetConfig, DataStreamConfig, PayloadGroup
# Functions
import RIOPA: read_config, default_config, get_evolution_function
# Modules
import RIOPA: DataGen, IO

function configure_stream(stream_cfg::Config)
    evolve_func = get_evolution_function(get(stream_cfg, :evolution, nothing))
    groups = map(
        grp_cfg -> PayloadGroup(grp_cfg[:size_range], grp_cfg[:ratio]),
        stream_cfg[:proc_payloads],
    )
    return DataStreamConfig(stream_cfg[:name], evolve_func, groups)
end

function configure_dataset(sub_cfg::Config)
    compute_seconds = get(sub_cfg, :compute_seconds, 0)
    datagen_tag_str = get(sub_cfg, :datagen_backend, nothing)
    io_tag_str = get(sub_cfg, :io_backend, nothing)
    step_factor = get(sub_cfg, :step_conversion_factor, 1)
    streams = map(configure_stream, sub_cfg[:data_streams])
    ds = DataSet(
        DataSetConfig(
            sub_cfg[:name],
            sub_cfg[:basename],
            DataGen.get_tag(datagen_tag_str),
            IO.get_tag(io_tag_str),
            sub_cfg[:nsteps],
            step_factor,
            compute_seconds,
            streams,
        ),
    )
    DataGen.initialize_streams!(ds.cfg.datagen_backend_tag, ds)
    return ds
end

generate_data!(ds::DataSet) = DataGen.generate!(ds.cfg.datagen_backend_tag, ds)

function perform_io_step!(ds::DataSet)
    IO.perform_step(ds.cfg.io_backend_tag, ds)
    ds.curr_step += 1
    ds.timestamp = time()
end

function steps_remain(ds::DataSet)
    return ds.curr_step <= ds.cfg.nsteps
end

struct Controller
    datasets::Vector{DataSet}
end

function (c::Controller)()
    cfgtimes = map(ds -> ds.cfg.compute_seconds, c.datasets)
    enabled = map(steps_remain, c.datasets)
    get_reset_time = (i, enable) -> enable ? cfgtimes[i] : Inf
    times = map(((i, enable),) -> get_reset_time(i, enable), enumerate(enabled))
    while any(enabled)
        ta = time()
        dur, idx = findmin(times)
        ds = c.datasets[idx]
        generate_data!(ds)
        diff = time() - ta
        dur = max(0, dur - diff)
        sleep(dur)
        times .-= dur
        perform_io_step!(ds)
        enabled[idx] = steps_remain(ds)
        times[idx] = get_reset_time(idx, enabled[idx])
    end
end

function run(config::Config)
    ctrl = Controller(map(configure_dataset, config[:datasets]))
    ctrl()
end

run(filename::AbstractString) = run(read_config(filename))

run(::Nothing) = run(default_config())

end


text:
push!(LOAD_PATH,"../src/")

using Documenter, SimpleWebsockets

makedocs(
    sitename = "SimpleWebsockets",
    format = Documenter.HTML(),
    modules = [SimpleWebsockets],
    pages = [
        "Introduction" => "index.md",
        "Server" => [
            "Server Usage" => "WebsocketServer.md",
            "Server Options" => "ServerOptions.md",
        ],
        "Client" => [
            "Client Usage" => "WebsocketClient.md",
            "Client Options" => "ClientOptions.md",
        ],
        "Websocket Connection" => "WebsocketConnection.md",
        "Authorization" => "Authorization.md",
        "Error handling" => "Errors.md",
        "Acknowledgments" => "Acknowledgments.md"
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/citkane/SimpleWebsockets.jl.git"
)

max_stars_repo_name: citkane/SimpleWebsockets

max_issues_repo_name: citkane/SimpleWebsockets

max_forks_repo_name: citkane/SimpleWebsockets


text:
using BayesianQuadrature
using Documenter

makedocs(;
    modules=[BayesianQuadrature],
    authors="Theo Galy-Fajou <theo.galyfajou@gmail.com> and contributors",
    repo="https://github.com/theogf/BayesianQuadrature.jl/blob/{commit}{path}#L{line}",
    sitename="BayesianQuadrature.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://theogf.github.io/BayesianQuadrature.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Theory" => "theory.md",
        "User Guide" => "userguide.md",
        "API" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/theogf/BayesianQuadrature.jl", devbranch="main", push_preview=true
)


text:
#type definitions

mutable struct SampledData
    DF::DataFrame
    VAF::Array{Float64,1}
    counts::Array{Int64,1}
    depth::Array{Int64,1}
end

mutable struct AnalysedData
    DF::DataFrame
    VAF::Array{Float64,1}
end

mutable struct RsqObj
    metric::Float64
    mu::Float64
    pval::Float64
end

mutable struct MetricObj
    metric::Float64
    pval::Float64
end

mutable struct AllMetrics
    rsq::RsqObj
    area::MetricObj
    Dk::MetricObj
    meanD::MetricObj
end

mutable struct Simulation
  input::InputParameters
  output::SimResult
  sampleddata::SampledData
end

mutable struct SimulationM
  input::InputParametersM
  output::SimResultM
  sampleddata::SampledData
end

mutable struct SimulationStemCells
  input::InputParameters
  output::StemCellSimResult
  sampleddata::SampledData
end

###############################################################################

function sampledhist(AF::Array{Float64, 1}, cellnum::Int64; detectionlimit = 0.1, ploidy = 2.0, read_depth = 100.0, cellularity = 1.0)

    AF = AF./ploidy
    AF = AF .* cellularity
    filter!(x -> x > detectionlimit * cellnum, AF)
    samp_percent = read_depth/cellnum
    #depth = rand(Binomial(cellnum,samp_percent), length(AF))
    depth = rand(Poisson(read_depth), length(AF))
    samp_alleles = map((n, p) -> rand(Binomial(n, p)), depth, AF/cellnum)

    VAF = samp_alleles./depth

    #data for histogram
    x = 0.005:0.01:1.005
    y = fit(Histogram, VAF, x, closed=:right)
    DFhist = DataFrame(VAF = x[1:end-1], freq = y.weights)

    SampledData(DFhist, VAF, samp_alleles, depth)
end

function betabinom(p, n, ρ)
    μ = p * n
    shape1 = (μ / n) * ((1 / ρ) - 1)
    shape2 = n * shape1/μ - shape1
    return rand(Binomial(n, rand(Beta(shape1, shape2))))
end

function sampledhist(AF, cellnum, ρ; detectionlimit = 0.1, ploidy = 2.0, read_depth = 100.0, cellularity = 1.0)

    AF = AF./ploidy
    AF = AF .* cellularity
    filter!(x -> x > detectionlimit * cellnum, AF)
    samp_percent = read_depth/cellnum
    #depth = rand(Binomial(cellnum, samp_percent), length(AF))
    depth = rand(Poisson(read_depth), length(AF))
    samp_alleles = map((x, y) -> betabinom(x, y, ρ), AF/cellnum, depth)
    VAF = samp_alleles./depth

    #data for histogram
    x = 0.005:0.01:1.005
    y = fit(Histogram, VAF, x, closed=:right)
    DFhist = DataFrame(VAF = x[1:end-1], freq = y.weights)

    SampledData(DFhist, VAF, samp_alleles, depth)
end

function cumulativedist(sresult; fmin = 0.1, fmax = 0.3)

    VAF = sresult.sampleddata.VAF

    #calculate cumulative sum
    steps = fmax:-0.001:fmin
    cumsum = Array{Int64}(undef, 0)
    v = Array{Float64}(undef, 0)

    for i in steps
        push!(cumsum, sum(VAF .>= i))
        push!(v, i)
    end
    cumsum = cumsum .- cumsum[1]

    DF = DataFrame(cumsum = map(Float64, cumsum), v = v)
    DF[:invf] = 1 ./ DF[:v] .- 1 ./ fmax

    DF[:theory] = Mcdf(DF[:v], fmin, fmax)
    DF[:normalized] = DF[:cumsum] ./ maximum(DF[:cumsum])

    lmfit = fit(LinearModel, @formula(cumsum ~ invf + 0), DF)
    DF[:prediction] = predict(lmfit)

    return AnalysedData(DF, VAF)
end

function Mcdf(f,fmin,fmax)
    (1.0./f .- 1.0/fmax) ./ (1.0/fmin .- 1.0/fmax)
end

function rsq(AD, fmin, fmax, metricp)

    #fit constrained fit
    lmfit = fit(LinearModel, @formula(cumsum ~ invf + 0), AD.DF)
    #calculate R^2 value
    rsqval = 1 .- (sum(residuals(lmfit) .^ 2) ./ sum((AD.DF[:cumsum] .- 0) .^ 2))
    #extract coefficient for mutation rate
    mu = coef(lmfit)[1]
    #get pvalue
    pval = metricp[:pval][searchsortedlast(metricp[:invrsqmetric],1 - rsqval)]

    return RsqObj(rsqval, mu, pval)
end

function kolmogorovD(AD, fmin, fmax, metricp)

    xfiltered = AD.VAF[fmin .< AD.VAF .< fmax]
    n = length(xfiltered)
    cdfs = 1 .- map(x->Mcdf(x, fmin, fmax), sort(xfiltered))
    δp = maximum((1:n) / n - cdfs)
    δn = -minimum((0:n-1) / n - cdfs)
    δ = max(δn, δp)
    pval = metricp[:pval][searchsortedlast(metricp[:Dkmetric], δ)]

    MetricObj(δ, pval)
end

function Draw(DF1, DF2)
    maximum(abs(DF1[:cumsum] .- DF2[:cumsum]))
end

function meanDraw(DF1, DF2)
    #mean(abs(DF1[:cumsum] - DF2[:cumsum]))
    sqrt(sum((DF1[:cumsum] .- DF2[:cumsum]).^2))
end

function kolmogorovDmean(AD, fmin, fmax, metricp)
    D = mean(abs.(convert(Array, AD.DF[:theory]) .- convert(Array, AD.DF[:normalized])))
    pval = metricp[:pval][searchsortedlast(metricp[:meanDmetric], D)]
    return MetricObj(D, pval)
end

function trapz(x::Vector{T}, y::Vector{T}) where T <: Number
    local n = length(x)
    if (length(y) != n)
        error("Vectors 'x', 'y' must be of same length")
    end
    if n == 1; return 0.0; end
    r = 0.0
    for i in 2:n
        r += (x[i] .- x[i-1]) .* (y[i] .+ y[i-1])
    end
    r / 2.0
end

function areametric(AD, fmin, fmax, metricp)
    area = abs.(trapz(convert(Array, AD.DF[:invf]), convert(Array, AD.DF[:normalized])) - trapz(convert(Array, AD.DF[:invf]), convert(Array, AD.DF[:theory])))
    area = area / (1 / fmin - 1 / fmax)
    pval = metricp[:pval][searchsortedlast(metricp[:areametric], area)]
    return MetricObj(area, pval)
end

function areametricraw(AD, DFABC; fmin = 0.12, fmax = 0.8)
    area = abs.(trapz(convert(Array, AD.DF[:v]), convert(Array, AD.DF[:cumsum])) - trapz(convert(Array, AD.DF[:v]), convert(Array, DFABC[:cumsum])))
    area = area / (1 / fmin - 1 / fmax)
    return area
end

"""
    simulate(; <keyword arguments>)

Simulate a stochastic model of tumour growth with a single subclone introduced at a random time and with a random fitness advantage. Output return synthetic sequencing data.
...
## Arguments
- `read_depth = 200.0`: Mean read depth of the target data set
- `detectionlimit = 5/read_depth`: Ability to detect low frequency variants. Assumes 5 reads are needed to call a variant.
- `μ = 10.0`: Mutation rate per division per genome (this will timesed by ploidy for the mutation rate per cell). At each division a Poisson random variable with mean μ is sampled.
- `clonalmutations = 100.0`: Number of clonal mutations present in the first cell.
- `nclones = 1`: Number of subclones introduced
- `Nmax = 10^4`: Maximum population size.
- `ρ = 0.0`: Overdispersion parameter for beta-binomial model of sequencing data. ρ = 0.0 means model is binomial sampling
- `timefunction = exptime`: Function for KMC algorithm timestep. exptime returns an exponentially distributed random variable, if you would rather return the mean of the distribution define a function that returns 1. ie `returnone() = 1`.
- `ploidy = 2`: ploidy of the genome
- `d = 0.0`: Death rate of the thost population in the tumour
- `b = log(2)`: Birth rate of the population. Set to `log(2)` so that tumour doubles with each unit increase in t in the absence of cell death.
- `fixedmu = false`: If set to false number of mutations per division is fixed and not sampled from a poisson distribution.
...
"""
function simulate(; nclones = 1, ploidy = 2, read_depth = 100.0, detectionlimit = 5 ./ read_depth, clonalmutations = 100.0, μ = 10.0, d = 0.0, b = log(2), ρ = 0.0, Nmax = 10^4, s = repeat([1.0], inner = nclones), tevent = collect(1.0:0.5:100.0)[1:nclones], cellularity = 1.0, fixedmu = false, timefunction::Function = exptime, maxclonefreq = 200)

    nclones == length(s) || error("Number of clones is $(nclones), size of selection coefficient array is $(length(s)), these must be the same size ")
    nclones == length(tevent) || error("Number of clones is $(nclones), size of selection coefficient array is $(length(tevent)), these must be the same size ")
    IP = InputParameters(nclones,
    Nmax,
    detectionlimit,
    ploidy,
    read_depth,
    clonalmutations,
    s,
    (μ/2) * ploidy,
    b,
    d,
    tevent,
    ρ,
    cellularity,
    fixedmu,
    timefunction,
    maxclonefreq)

    #get simulation data
    simresult, IP = run1simulation(IP, 0.0, 1.0)

    #get sampled VAFs
    if IP.ρ > 0.0
        sampleddata = sampledhist(simresult.trueVAF, IP.Nmax, IP.ρ, detectionlimit = IP.detectionlimit, ploidy = IP.ploidy, read_depth = IP.read_depth, cellularity = IP.cellularity)
    else
        sampleddata = sampledhist(simresult.trueVAF, IP.Nmax, detectionlimit = IP.detectionlimit, ploidy = IP.ploidy, read_depth = IP.read_depth, cellularity = IP.cellularity)
    end

    return Simulation(IP, simresult, sampleddata)
end

"""
    simulate(minclonesize, maxclonesize; <keyword arguments>)

Return simulation with frequency of subclones >minclones & <maxclonesize.
"""
function simulate(minclonesize, maxclonesize; nclones = 1, ploidy = 2, read_depth = 100.0, detectionlimit = 5 ./read_depth, clonalmutations = 100.0, μ = 10.0, d = 0.0, b = log(2), ρ = 0.0, Nmax = 10^4, cellularity = 1.0, fixedmu = false, tmin = 6.0, tmax = 20.0, smin = 0.0, smax = 25.0, timefunction::Function = exptime, maxclonefreq = 200)

    correctnc = false
    IP, simresult = 0, 0

    while correctnc == false
      freq = sort(rand(Uniform(minclonesize, maxclonesize), nclones), rev = true)
      tevent = sort(rand(Uniform(tmin, log(Nmax)/log(2)), nclones))
      if nclones == 1
          s1 = selection(log(2), minclonesize, log(Nmax)/log(2), tevent[1])
          s2 = selection(log(2), maxclonesize, log(Nmax)/log(2), tevent[1])
          s = rand(Uniform(s1, 3*s2), nclones) #increasing this s2 makes the distribution of clone sizes more uniform
      elseif nclones == 2
          s1, s2 = selection2clone(log(2), freq[1], freq[2], log(Nmax)/log(2), tevent[1], tevent[2])
          sm = maximum([s1, s2])
          s = rand(Uniform(s1[1], 3*sm[1]), nclones)
      else
          s = rand(Uniform(smin, smax), nclones)
      end

      IP = InputParameters(nclones,
      Nmax,
      detectionlimit,
      ploidy,
      read_depth,
      clonalmutations,
      s,
      (μ/2) * ploidy, #divide by 2 because both daughter cells can accumulate mutations so this ensures the *per division* mutation rate is accurate
      b,
      d,
      tevent,
      ρ,
      cellularity,
      fixedmu,
      timefunction,
      maxclonefreq)

      simresult, IP = run1simulation(IP, minclonesize, maxclonesize)

      nctemp = IP.numclones
      if nctemp == nclones
          correctnc = true
      end
    end
    #println(i)
    #get sampled VAFs
    if IP.ρ > 0.0
        sampleddata = sampledhist(simresult.trueVAF, IP.Nmax, IP.ρ, detectionlimit = IP.detectionlimit, ploidy = IP.ploidy, read_depth = IP.read_depth, cellularity = IP.cellularity)
    else
        sampleddata = sampledhist(simresult.trueVAF, IP.Nmax, detectionlimit = IP.detectionlimit, ploidy = IP.ploidy, read_depth = IP.read_depth, cellularity = IP.cellularity)
    end

    return Simulation(IP, simresult, sampleddata)
end

"""
    simulate(minclonesize, maxclonesize, independentclones::Bool; <keyword arguments>)

Return simulation with frequency of subclones >minclones & <maxclonesize and specify whether subclones are independent or not (ie nested or not). Only applicable to >1 subclone.
"""
function simulate(minclonesize, maxclonesize, independentclones::Bool; nclones = 1, ploidy = 2, read_depth = 100.0, detectionlimit = 5 ./read_depth, clonalmutations = 100.0, μ = 10.0, d = 0.0, b = log(2), ρ = 0.0, Nmax = 10^4, cellularity = 1.0, fixedmu = false, tmin = 3.0, tmax = 20.0, smin = 0.0, smax = 25.0, timefunction::Function = exptime, maxclonefreq = 200)

  ct = 1
  x = 0.0

  while (ct >= 1)
    x = simulate(minclonesize, maxclonesize; nclones = nclones, ploidy = ploidy, read_depth = read_depth, detectionlimit = detectionlimit, clonalmutations = clonalmutations, μ = μ, d = d, b = b, ρ = ρ, Nmax = Nmax, cellularity = cellularity, fixedmu = fixedmu, tmin = tmin, tmax = tmax, smin = smin, smax = smax, timefunction = timefunction, maxclonefreq = maxclonefreq)
    ct = sum(x.output.clonetype)
  end

  return x
end

"""
    simulate(minclonesize, maxclonesize, mindiff; <keyword arguments>)

Return simulation with frequency of subclones >minclones & <maxclonesize and subclone frequency are at least mindiff apart.
"""
function simulate(minclonesize, maxclonesize, mindiff::Float64; nclones = 1, ploidy = 2, read_depth = 100.0, detectionlimit = 5 ./read_depth, clonalmutations = 100.0, μ = 10.0, d = 0.0, b = log(2), ρ = 0.0, Nmax = 10^4, cellularity = 1.0, fixedmu = false, tmin = 3.0, tmax = 20.0, smin = 0.0, smax = 25.0, timefunction::Function = exptime, maxclonefreq = maxclonefreq)

  nclones == 2 || error("nclones must be = 2 for this method as it is a function to simulate until we arrive at 2 clones that are greater than mindiff apart")

  ct = 1
  x = 0.0
  freqdiff = false

  while (freqdiff == false)
    x = simulate(minclonesize, maxclonesize; nclones = nclones, ploidy = ploidy, read_depth = read_depth, detectionlimit = detectionlimit, clonalmutations = clonalmutations, μ = μ, d = d, b = b, ρ = ρ, Nmax = Nmax, cellularity = cellularity, fixedmu = fixedmu, tmin = tmin, tmax = tmax, smin = smin, smax = smax, timefunction = timefunction, maxclonefreq = maxclonefreq)
    if abs(x.output.clonefreq[2] - x.output.clonefreq[1]) > mindiff
      freqdiff = true
    end
  end

  return x
end

"""
    simulate(minclonesize, maxclonesize, mindiff; <keyword arguments>)

Return simulation with frequency of subclones >minclones & <maxclonesize and subclone frequency are at least mindiff and have at least minmutations mutations.
"""
function simulate(minclonesize, maxclonesize, mindiff::Float64, minmutations::Int64; nclones = 1, ploidy = 2, read_depth = 100.0, detectionlimit = 5 ./read_depth, clonalmutations = 100.0, μ = 10.0, d = 0.0, b = log(2), ρ = 0.0, Nmax = 10^4, cellularity = 1.0, fixedmu = false, tmin = 3.0, tmax = 20.0, smin = 0.0, smax = 25.0, timefunction::Function = exptime, maxclonefreq = maxclonefreq)

  nclones == 2 || error("nclones must be = 2 for this method as it is a function to simulate until we arrive at 2 clones that are greater than mindiff apart")

  ct = 1
  x = 0.0
  savesim = false

  while (savesim == false)
    x = simulate(minclonesize, maxclonesize, mindiff; nclones = nclones, ploidy = ploidy, read_depth = read_depth, detectionlimit = detectionlimit, clonalmutations = clonalmutations, μ = μ, d = d, b = b, ρ = ρ, Nmax = Nmax, cellularity = cellularity, fixedmu = fixedmu, tmin = tmin, tmax = tmax, smin = smin, smax = smax, timefunction = timefunction, maxclonefreq = maxclonefreq)
    if sum(x.output.subclonemutations .> minmutations) == 2
      savesim = true
    end
  end

  return x
end

function getmetrics(AD, metricp; fmin = 0.1, fmax = 0.3)

    rsq1 = rsq(AD, fmin, fmax, metricp)
    Dk = kolmogorovD(AD, fmin, fmax, metricp)
    meanD = kolmogorovDmean(AD, fmin, fmax, metricp)
    area = areametric(AD, fmin, fmax, metricp)

    return(AllMetrics(rsq1, area, Dk, meanD))
end

"""
    simulatestemcells(; <keyword arguments>)

Simulate a stochastic model of tumour growth with a stem cell architecture. Output return synthetic sequencing data.
...
## Arguments
- `α = 0.1`: Symmetric division rate
- `maxdivisions = 5`: Maximum number of divisions of differentiated cells.
- `read_depth = 200.0`: Mean read depth of the target data set
- `detectionlimit = 5/read_depth`: Ability to detect low frequency variants. Assumes 5 reads are needed to call a variant.
- `μ = 10.0`: Mutation rate per division. At each division a Poisson random variable with mean μ is sampled.
- `clonalmutations = 100.0`: Number of clonal mutations present in the first cell.
- `Nmax = 10^4`: Maximum population size.
- `ρ = 0.0`: Overdispersion parameter for beta-binomial model of sequencing data. ρ = 0.0 means model is binomial sampling
- `timefunction = exptime`: Function for KMC algorithm timestep. exptime returns an exponentially distributed random variable, if you would rather return the mean of the distribution define a function that returns 1. ie `returnone() = 1`.
- `ploidy = 2`: ploidy of the genome
- `d = 0.0`: Death rate of the thost population in the tumour
- `b = log(2)`: Birth rate of the population. Set to `log(2)` so that tumour doubles with each unit increase in t in the absence of cell death.
- `fixedmu = false`: If set to false number of mutations per division is fixed and not sampled from a poisson distribution.
...
"""
function simulatestemcells(;ploidy = 2, α = 0.1, maxdivisions = 5, read_depth = 100.0, nclones = 0, detectionlimit = 5 ./read_depth, clonalmutations = 100.0, μ = 10.0, d = 0.0, b = log(2), ρ = 0.0, Nmax = 10^4, s = repeat([1.0], inner = nclones), tevent = collect(1.0:0.5:100.0)[1:nclones], cellularity = 1.0, fixedmu = false, timefunction::Function = exptime)

    nclones == length(s) || error("Number of clones is $(nclones), size of selection coefficient array is $(length(s)), these must be the same size ")

    nclones == length(tevent) || error("Number of clones is $(nclones), size of selection coefficient array is $(length(tevent)), these must be the same size ")

    IP = InputParameters(nclones,
    Nmax,
    detectionlimit,
    ploidy,
    read_depth,
    clonalmutations,
    s,
    μ,
    b,
    d,
    tevent,
    ρ,
    cellularity,
    fixedmu,
    timefunction)

    #get simulation data
    simresult = run1simulationstem(Nmax; α = α, maxdivisions = maxdivisions, d = d, μ = μ, clonalmutations = clonalmutations)

    #get sampled VAFs
    if IP.ρ > 0.0
        sampleddata = sampledhist(simresult.trueVAF, IP.Nmax, IP.ρ, detectionlimit = IP.detectionlimit, ploidy = IP.ploidy, read_depth = IP.read_depth, cellularity = IP.cellularity)
    else
        sampleddata = sampledhist(simresult.trueVAF, IP.Nmax, detectionlimit = IP.detectionlimit, ploidy = IP.ploidy, read_depth = IP.read_depth, cellularity = IP.cellularity)
    end

    return SimulationStemCells(IP, simresult, sampleddata)
end

function simulatedifferentmutations(;ploidy = 2, read_depth = 100.0, detectionlimit = 5 ./read_depth, clonalmutations = 100.0, μp = 10.0, μd = 0.001, μneg = 0.0, d = 0.0, b = log(2), ρ = 0.0, Nmax = 10^4, cellularity = 1.0, timefunction::Function = exptime, s = sfunc, fitnessfunc = nonmultiplicativefitness)

    IP = InputParametersM(
    Nmax,
    detectionlimit,
    ploidy,
    read_depth,
    clonalmutations,
    μp,
    μd,
    μneg,
    b,
    d,
    ρ,
    cellularity,
    s,
    timefunction,
    fitnessfunc)

    #get simulation data
    simresult, IP = run1simulation(IP)

    #get sampled VAFs
    if IP.ρ > 0.0
        sampleddata = sampledhist(simresult.trueVAF, IP.Nmax, IP.ρ, detectionlimit = IP.detectionlimit, ploidy = IP.ploidy, read_depth = IP.read_depth, cellularity = IP.cellularity)
    else
        sampleddata = sampledhist(simresult.trueVAF, IP.Nmax, detectionlimit = IP.detectionlimit, ploidy = IP.ploidy, read_depth = IP.read_depth, cellularity = IP.cellularity)
    end

    return SimulationM(IP, simresult, sampleddata)
end

function simulatedifferentmutationsmoran(tmax; ploidy = 2, read_depth = 100.0, detectionlimit = 5 ./read_depth, clonalmutations = 100.0, μp = 10.0, μd = 0.001, μneg = 0.0, d = 0.0, b = log(2), ρ = 0.0, Nmax = 10^4, cellularity = 1.0, timefunction::Function = exptime, s = sfunc, fitnessfunc = nonmultiplicativefitness)

    IP = InputParametersM(
    Nmax,
    detectionlimit,
    ploidy,
    read_depth,
    clonalmutations,
    μp,
    μd,
    μneg,
    b,
    d,
    ρ,
    cellularity,
    s,
    timefunction,
    nonmultiplicativefitness)

    #get simulation data
    simresult, IP = run1simulation(IP, tmax)

    #get sampled VAFs
    if IP.ρ > 0.0
        sampleddata = sampledhist(simresult.trueVAF, IP.Nmax, IP.ρ, detectionlimit = IP.detectionlimit, ploidy = IP.ploidy, read_depth = IP.read_depth, cellularity = IP.cellularity)
    else
        sampleddata = sampledhist(simresult.trueVAF, IP.Nmax, detectionlimit = IP.detectionlimit, ploidy = IP.ploidy, read_depth = IP.read_depth, cellularity = IP.cellularity)
    end

    return SimulationM(IP, simresult, sampleddata)
end

function getsummary(inandout; sname = "", fmin = 0.1, fmax = 0.3)

  simresult = inandout.output
  IP = inandout.input
  sampleddata = inandout.sampleddata
  #get cumulativedistributions
  AD = cumulativedist(inandout, fmin = fmin, fmax = fmax)
  allmetrics = getmetrics(AD, metricp, fmin = fmin, fmax = fmax)
  fitout = rsq(AD, fmin, fmax, metricp)

  DF = DataFrame(
  sname = sname,
  Nmax = IP.Nmax,
  numclones = IP.numclones,
  Nevent = join(simresult.cloneN, ","),
  teventin = join(IP.tevent, ","),
  s = join(IP.selection, ","),
  clonefreq = join(simresult.clonefreq, ","),
  teventtrue = join(simresult.clonetime, ","),
  subclonemutations = join(simresult.subclonemutations, ","),
  clonetype = join(simresult.clonetype, ","),
  br = join(simresult.birthrates, ","),
  dr = join(simresult.deathrates, ","),
  mu = IP.μ,
  muout = fitout.mu,
  b = IP.b,
  d = IP.d,
  tend = simresult.tend,
  clonalmutations = IP.clonalmutations,
  area = allmetrics.area.metric,
  area_pval = allmetrics.area.pval,
  Dk = allmetrics.Dk.metric,
  Dk_pval = allmetrics.Dk.pval,
  meanD = allmetrics.meanD.metric,
  meanD_pval = allmetrics.meanD.pval,
  rsq = allmetrics.rsq.metric,
  rsq_pval = allmetrics.rsq.pval,
  fmin = fmin,
  fmax = fmax,
  detectionlimit = IP.detectionlimit,
  ploidy = IP.ploidy,
  read_depth = IP.read_depth,
  num_muts = AD.DF[:cumsum][end]
  )

  return DF

end

function getsummary(inandout::SimulationStemCells; sname = "", fmin = 0.1, fmax = 0.3)

  simresult = inandout.output
  IP = inandout.input
  sampleddata = inandout.sampleddata
  #get cumulativedistributions
  AD = cumulativedist(inandout, fmin = fmin, fmax = fmax)
  allmetrics = getmetrics(AD, metricp, fmin = fmin, fmax = fmax)
  fitout = rsq(AD, fmin, fmax, metricp)

  DF = DataFrame(
  sname = sname,
  Nmax = IP.Nmax,
  numclones = IP.numclones,
  mu = IP.μ,
  muout = fitout.mu,
  clonalmutations = IP.clonalmutations,
  area = allmetrics.area.metric,
  area_pval = allmetrics.area.pval,
  Dk = allmetrics.Dk.metric,
  Dk_pval = allmetrics.Dk.pval,
  meanD = allmetrics.meanD.metric,
  meanD_pval = allmetrics.meanD.pval,
  rsq = allmetrics.rsq.metric,
  rsq_pval = allmetrics.rsq.pval,
  fmin = fmin,
  fmax = fmax,
  detectionlimit = IP.detectionlimit,
  ploidy = IP.ploidy,
  read_depth = IP.read_depth,
  num_muts = AD.DF[:cumsum][end]
  )

  return DF

end


text:
module ErlangTerm

export serialize, deserialize

const NEW_FLOAT = UInt8(70)
const SMALL_INTEGER = UInt8(97)
const INTEGER = UInt8(98)
const ATOM = UInt8(100)
const SMALL_TUPLE = UInt8(104)
const LARGE_TUPLE = UInt8(105)
const NIL = UInt8(106)
const STRING = UInt8(107)
const LIST = UInt8(108)
const BINARY = UInt8(109)
const SMALL_BIG = UInt8(110)
const LARGE_BIG = UInt8(111)
const MAP = UInt8(116)
const ATOM_UTF8 = UInt8(118)
const SMALL_ATOM_UTF8 = UInt8(119)

const VERSION = UInt8(131)

"""
    deserialize(io::IO)
    deserialize(binary::Array{UInt8})

Deserialize a Julia value encoded in the
[Erlang external term format](http://erlang.org/doc/apps/erts/erl_ext_dist.html)
from a stream or an array of bytes.
"""
deserialize

"""
    serialize(io::IO, value)
    serialize(value)

Serialize a Julia value to the
[Erlang external term format](http://erlang.org/doc/apps/erts/erl_ext_dist.html).
"""
serialize

function deserialize(io::IO)
    version = read(io, UInt8)
    version != VERSION && throw(ArgumentError("Unknown external term format version: $(Int(version))."))

    tag = read(io, UInt8)
    deserialize(io, Val(tag))
end

deserialize(binary::Array{UInt8}) = deserialize(IOBuffer(binary))

function serialize(io::IO, data)
    write(io, VERSION)
    _serialize(io, data)
end

serialize(data) = take!(serialize(IOBuffer(), data))

function _serialize(io, val::Integer)
    if val > 0 && val <= typemax(UInt8)
        write(io, SMALL_INTEGER)
        write(io, UInt8(val))
    elseif abs(val) <= typemax(Int32)
        write(io, INTEGER)
        write(io, hton(Int32(val)))
    else
        arr = digits(UInt8, abs(val), base=0x100)
        n = length(arr)
        if n < 0x100
            write(io, SMALL_BIG)
            write(io, UInt8(n))
        else
            write(io, LARGE_BIG)
            write(io, hton(UInt32(n)))
        end
        write(io, UInt8(sign(val) == 1 ? 0 : 1))
        write(io, arr)
    end
    io
end

function deserialize(io, ::Val{SMALL_INTEGER})
    Int(read(io, UInt8))
end

function deserialize(io, ::Val{INTEGER})
    ntoh(read(io, Int32))
end

function _serialize(io, val::Float64)
    write(io, NEW_FLOAT)
    write(io, hton(val))
    io
end

function deserialize(io, ::Val{NEW_FLOAT})
    ntoh(read(io, Float64))
end

function _serialize(io, val::Symbol)
    str = Array{UInt8}(string(val))
    n = length(str)
    if n > typemax(UInt16)
        throw(ArgumentError("Cannot serialize symbols with more than $(typemax(UInt16)) characters."))
    elseif n > typemax(UInt8)
        write(io, ATOM_UTF8)
        write(io, hton(UInt16(n)))
    else
        write(io, SMALL_ATOM_UTF8)
        write(io, hton(UInt8(n)))
    end
    write(io, str)
    io
end

function deserialize(io, ::Val{ATOM})
    n = Int(ntoh(read(io, UInt16)))
    Symbol(String(read(io, n)))
end

function deserialize(io, ::Val{ATOM_UTF8})
    n = Int(ntoh(read(io, UInt16)))
    Symbol(String(read(io, n)))
end

function deserialize(io, ::Val{SMALL_ATOM_UTF8})
    n = Int(ntoh(read(io, UInt8)))
    Symbol(String(read(io, n)))
end

function _serialize(io, val::String)
    write(io, BINARY)
    str = Array{UInt8}(val)
    n = length(str)
    n > typemax(UInt32) && throw(ArgumentError("Strings longer than $(typemax(UInt32)) cannot be serialized."))
    write(io, hton(UInt32(n)))
    write(io, str)
    io
end

function deserialize(io, ::Val{BINARY})
    n = Int(ntoh(read(io, UInt32)))
    String(read(io, n))
end

function deserialize(io, ::Val{SMALL_BIG})
    n = Int(read(io, UInt8))
    sign = Int32(read(io, UInt8) == 0 ? 1 : -1)
    arr = Array{UInt8,1}(read(io, n))
    _bignum(n, sign, arr)
end

function deserialize(io, ::Val{LARGE_BIG})
    n = Int(ntoh(read(io, UInt32)))
    sign = Int32(read(io, UInt8) == 0 ? 1 : -1)
    arr = Array{UInt8,1}(read(io, n))
    _bignum(n, sign, arr)
end

function _bignum(n, sign, arr)
    itype = n < 4 ? Int32 :
            n < 8 ? Int64 :
            n < 16 ? Int128 : BigInt
    sign * dig2int(itype, arr)
end

dig2int(arr) = sum(((i, x),) -> x * 256^(i - 1), pairs(arr))
dig2int(::Type{Int32}, arr) = Int32(dig2int(arr))
dig2int(::Type{Int64}, arr) = sum(((i, x),) -> x * Int64(256)^(i - 1), pairs(arr))
dig2int(::Type{Int128}, arr) = sum(((i, x),) -> x * Int128(256)^(i - 1), pairs(arr))
dig2int(::Type{BigInt}, arr) = sum(((i, x),) -> x * BigInt(256)^(i - 1), pairs(arr))

function _serialize(io, array::AbstractArray)
    if isempty(array)
        write(io, NIL)
        return io
    end

    write(io, LIST)
    n = length(array)
    n > typemax(UInt32) && throw(ArgumentError("Arrays longer than $(typemax(UInt32)) cannot be serialized."))
    write(io, hton(UInt32(n)))
    for el in array
        _serialize(io, el)
    end
    write(io, NIL)
    io
end

function deserialize(io, ::Val{NIL})
    []
end

function deserialize(io, ::Val{LIST})
    n = Int(ntoh(read(io, UInt32)))
    i = 0
    array = []
    while i < n
        tag = read(io, UInt8)
        push!(array, deserialize(io, Val(tag)))
        i += 1
    end
    tail = read(io, UInt8)
    tail != NIL && throw(ArgumentError("Improper lists are not supported."))
    array
end

function _serialize(io, dict::Dict)
    write(io, MAP)
    n = length(dict)
    n > typemax(UInt32) && throw(ArgumentError("Dicts with more than $(typemax(UInt32)) pairs cannot be serialized."))
    write(io, hton(UInt32(n)))
    for (key, value) in dict
        _serialize(io, key)
        _serialize(io, value)
    end
    io
end

function deserialize(io, ::Val{MAP})
    n = Int(ntoh(read(io, UInt32)))
    i = 0
    dict = Dict()
    while i < n
        keytag = read(io, UInt8)
        key = deserialize(io, Val(keytag))
        valuetag = read(io, UInt8)
        value = deserialize(io, Val(valuetag))
        push!(dict, key=>value)
        i += 1
    end
    dict
end

function _serialize(io, val::Tuple)
    n = length(val)
    if n > typemax(UInt32)
        throw(ArgumentError("Tuples longer than $(typemax(UInt32)) cannot be serialized."))
    elseif n > typemax(UInt8)
        write(io, LARGE_TUPLE)
        write(io, hton(UInt32(n)))
    else
        write(io, SMALL_TUPLE)
        write(io, hton(UInt8(n)))
    end
    for el in val
        _serialize(io, el)
    end
    io
end

function deserialize(io, ::Val{SMALL_TUPLE})
    n = Int(ntoh(read(io, UInt8)))
    i = 0
    array = []
    while i < n
        tag = read(io, UInt8)
        push!(array, deserialize(io, Val(tag)))
        i += 1
    end
    Tuple(array)
end

function deserialize(io, ::Val{LARGE_TUPLE})
    n = Int(ntoh(read(io, UInt32)))
    i = 0
    array = []
    while i < n
        tag = read(io, UInt8)
        push!(array, deserialize(io, Val(tag)))
        i += 1
    end
    Tuple(array)
end

function deserialize(io, ::Val{STRING})
    n = Int(ntoh(read(io, UInt16)))
    i = 0
    array = []
    while i < n
        push!(array, ntoh(read(io, UInt8)))
        i += 1
    end
    array
end

function _serialize(io, ::Nothing)
    _serialize(io, :nil)
    io
end

function _serialize(io, unknown)
    _serialize(io, repr(unknown))
end

end # module


text:
selector(items) = @rpc select(items)

input() = @rpc input()

info(s) = @msg info(string(s))

"""
  progress(x = [0..1])

Set Atom's progress bar to the given value.
"""
progress(x::Void = nothing) = @msg progress(x)

progress(x::Real) =
  @msg progress(x < 0.01 ? nothing :
                x > 1 ? 1 :
                x)

macro !(ex)
  quote
    result = $(esc(ex))
    display(span(c(render(Inline(), $(Expr(:quote, ex))),
                   " = ",
                   render(Inline(), result))))
    result
  end
end

# Blink stuff

type Shell <: AtomShell.Shell end

AtomShell.active(::Shell) = true

AtomShell.raw_window(::Shell, opts) =
  @rpc createWindow(merge(AtomShell.window_defaults, opts))

AtomShell.dot(::Shell, win::Integer, code; callback = true) =
  (callback ? rpc : msg)(:withWin, win, Blink.jsstring(code))

AtomShell.active(::Shell, win::Integer) = @rpc winActive(win)

plotsize() = @rpc plotsize()

ploturl(url::String) = @msg ploturl(url)

function blinkplot()
  p = Page()
  ploturl(Blink.localurl(p))
  return wait(p)
end


text:
using RemarkableFloats
using Test

a64 = Float64(0.5)
a32 = Float32(0.5)
a16 = Float16(0.5)

@test getmark(a64) === Marker0
@test getmark(a32) === Marker0
@test getmark(a16) === Marker0

@test getmarkvalue(a64) === 0
@test getmarkvalue(a32) === 0
@test getmarkvalue(a16) === 0

@test has0marks(a64)
@test has0marks(a32)
@test has0marks(a16)

b64 = setmark(a64, Marker1)
b32 = setmark(a32, Marker2)
b16 = setmark(a16, Marker3)

@test getmark(b64) === Marker1
@test getmark(b32) === Marker2
@test getmark(b16) === Marker3

@test getmarkvalue(b64) === 1
@test getmarkvalue(b32) === 2
@test getmarkvalue(b16) === 3

@test has1mark(b64)
@test has2marks(b32)
@test has3marks(b16)

@test has1or2marks(b64)
@test has1or2marks(b32)
@test has2or3marks(b16)
@test has1or3marks(b16)

a64 = Float64(0.5)
a32 = Float32(0.5)
a16 = Float16(0.5)

@setmark!(a64, Marker2)
@incrmark!(a32)

@test has2marks(a64)
@test has1mark(a32)

@decrmark!(a64)
@unmark!(a32)

@test has1mark(a64)
@test has0marks(a32)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct DiskUpdateProperties <: SwaggerModel
    osType::Any # spec type: Union{ Nothing, String } # spec name: osType
    diskSizeGB::Any # spec type: Union{ Nothing, Int32 } # spec name: diskSizeGB
    encryptionSettingsCollection::Any # spec type: Union{ Nothing, EncryptionSettingsCollection } # spec name: encryptionSettingsCollection
    diskIOPSReadWrite::Any # spec type: Union{ Nothing, Int64 } # spec name: diskIOPSReadWrite
    diskMBpsReadWrite::Any # spec type: Union{ Nothing, Int64 } # spec name: diskMBpsReadWrite
    diskIOPSReadOnly::Any # spec type: Union{ Nothing, Int64 } # spec name: diskIOPSReadOnly
    diskMBpsReadOnly::Any # spec type: Union{ Nothing, Int64 } # spec name: diskMBpsReadOnly
    maxShares::Any # spec type: Union{ Nothing, Int32 } # spec name: maxShares
    encryption::Any # spec type: Union{ Nothing, Encryption } # spec name: encryption
    networkAccessPolicy::Any # spec type: Union{ Nothing, NetworkAccessPolicy } # spec name: networkAccessPolicy
    diskAccessId::Any # spec type: Union{ Nothing, String } # spec name: diskAccessId

    function DiskUpdateProperties(;osType=nothing, diskSizeGB=nothing, encryptionSettingsCollection=nothing, diskIOPSReadWrite=nothing, diskMBpsReadWrite=nothing, diskIOPSReadOnly=nothing, diskMBpsReadOnly=nothing, maxShares=nothing, encryption=nothing, networkAccessPolicy=nothing, diskAccessId=nothing)
        o = new()
        validate_property(DiskUpdateProperties, Symbol("osType"), osType)
        setfield!(o, Symbol("osType"), osType)
        validate_property(DiskUpdateProperties, Symbol("diskSizeGB"), diskSizeGB)
        setfield!(o, Symbol("diskSizeGB"), diskSizeGB)
        validate_property(DiskUpdateProperties, Symbol("encryptionSettingsCollection"), encryptionSettingsCollection)
        setfield!(o, Symbol("encryptionSettingsCollection"), encryptionSettingsCollection)
        validate_property(DiskUpdateProperties, Symbol("diskIOPSReadWrite"), diskIOPSReadWrite)
        setfield!(o, Symbol("diskIOPSReadWrite"), diskIOPSReadWrite)
        validate_property(DiskUpdateProperties, Symbol("diskMBpsReadWrite"), diskMBpsReadWrite)
        setfield!(o, Symbol("diskMBpsReadWrite"), diskMBpsReadWrite)
        validate_property(DiskUpdateProperties, Symbol("diskIOPSReadOnly"), diskIOPSReadOnly)
        setfield!(o, Symbol("diskIOPSReadOnly"), diskIOPSReadOnly)
        validate_property(DiskUpdateProperties, Symbol("diskMBpsReadOnly"), diskMBpsReadOnly)
        setfield!(o, Symbol("diskMBpsReadOnly"), diskMBpsReadOnly)
        validate_property(DiskUpdateProperties, Symbol("maxShares"), maxShares)
        setfield!(o, Symbol("maxShares"), maxShares)
        validate_property(DiskUpdateProperties, Symbol("encryption"), encryption)
        setfield!(o, Symbol("encryption"), encryption)
        validate_property(DiskUpdateProperties, Symbol("networkAccessPolicy"), networkAccessPolicy)
        setfield!(o, Symbol("networkAccessPolicy"), networkAccessPolicy)
        validate_property(DiskUpdateProperties, Symbol("diskAccessId"), diskAccessId)
        setfield!(o, Symbol("diskAccessId"), diskAccessId)
        o
    end
end # type DiskUpdateProperties

const _property_map_DiskUpdateProperties = Dict{Symbol,Symbol}(Symbol("osType")=>Symbol("osType"), Symbol("diskSizeGB")=>Symbol("diskSizeGB"), Symbol("encryptionSettingsCollection")=>Symbol("encryptionSettingsCollection"), Symbol("diskIOPSReadWrite")=>Symbol("diskIOPSReadWrite"), Symbol("diskMBpsReadWrite")=>Symbol("diskMBpsReadWrite"), Symbol("diskIOPSReadOnly")=>Symbol("diskIOPSReadOnly"), Symbol("diskMBpsReadOnly")=>Symbol("diskMBpsReadOnly"), Symbol("maxShares")=>Symbol("maxShares"), Symbol("encryption")=>Symbol("encryption"), Symbol("networkAccessPolicy")=>Symbol("networkAccessPolicy"), Symbol("diskAccessId")=>Symbol("diskAccessId"))
const _property_types_DiskUpdateProperties = Dict{Symbol,String}(Symbol("osType")=>"String", Symbol("diskSizeGB")=>"Int32", Symbol("encryptionSettingsCollection")=>"EncryptionSettingsCollection", Symbol("diskIOPSReadWrite")=>"Int64", Symbol("diskMBpsReadWrite")=>"Int64", Symbol("diskIOPSReadOnly")=>"Int64", Symbol("diskMBpsReadOnly")=>"Int64", Symbol("maxShares")=>"Int32", Symbol("encryption")=>"Encryption", Symbol("networkAccessPolicy")=>"NetworkAccessPolicy", Symbol("diskAccessId")=>"String")
Base.propertynames(::Type{ DiskUpdateProperties }) = collect(keys(_property_map_DiskUpdateProperties))
Swagger.property_type(::Type{ DiskUpdateProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DiskUpdateProperties[name]))}
Swagger.field_name(::Type{ DiskUpdateProperties }, property_name::Symbol) =  _property_map_DiskUpdateProperties[property_name]

const _allowed_DiskUpdateProperties_osType = ["Windows", "Linux"]

function check_required(o::DiskUpdateProperties)
    true
end

function validate_property(::Type{ DiskUpdateProperties }, name::Symbol, val)
    if name === Symbol("osType")
        Swagger.validate_param(name, "DiskUpdateProperties", :enum, val, _allowed_DiskUpdateProperties_osType)
    end
end


text:
abstract type AbstractFDI <: AbstractEnv end
include("DelayFDI.jl")
include("LPFFDI.jl")


text:
include("./LogData.jl")

struct Simulation

  psf :: String
  pdb :: String
  dcd :: String

  natoms :: Int64
  atoms :: Vector{Atom}

  nbonds :: Int64
  bonds :: Vector{Bond}

  nanbles :: Int64
  angles :: Vector{Angle}

  ndihedrals :: Int64
  dihedrals :: Vector{Dihedral}

  nimpropers :: Int64
  improper :: Vector{Improper}

  nframes :: Int64
  dcdaxis :: Bool

  vmd :: String

  FortranDCD :: FortranFile

  logfile :: String
  log :: LogData

  parfiles :: Vector{String}

end

function Base.show( io :: IO, simulation :: Simulation )
  println(" Namd simulation with ", simulation.natoms," atoms. ")
  if simulation.logfile != nothing
    println("   LOG file: ", simulation.logfile )
  end
  if simulation.psf != nothing
    println("   PSF file: ", simulation.psf )
  end
  if simulation.psf != nothing
    println("   DCD file: ", simulation.dcd )
  end
end


text:
import ACE.Transforms: GetVal

@doc raw"""
`Ylm1pBasis`

One-particle basis component of the form
```math
\phi_{lm}({\bm r}) = Y_l^m(\hat{\br r})
```
Fundamental building block of ACE basis sets of the form
```math
   R_{nl}^{\mu_i, \mu_j}(r_{ij}) Y_l^m(\hat{\bm r})
```
This type basically just translates the `SHBasis` into a valid one-particle
basis.
"""
function Ylm1pBasis(L::Integer; varsym = :rr, lsym = :l, msym = :m,
                                  label = "Y$lsym$msym")
   SH = SphericalHarmonics.SHBasis(L)
   NTLM = NamedTuple{(lsym, msym)}
   spec = [ NTLM(tuple(SphericalHarmonics.idx2lm(i)...))
            for i = 1:length(SH) ]
   degrees = [ b[lsym] for b in spec ]
   return B1pComponent(SH, GetVal{varsym}(), spec, degrees, label)
end


text:
"""
    ElectronOutput <: AbstractInteractOutput

    ElectronOutput(init; ruleset, tspan, kw...)

A html output using Interact.jl and an Electron window through Blink.jl
ElectronOutput automatically generates sliders to control simulations
in realtime. args and kwargs are passed to [`InteractOutput`](@ref).

# Arguments

- `init`: initialisation `AbstractArray` or `NamedTuple` of `AbstractArray`.

# Keywords

$INTERACT_OUTPUT_KEYWORDS
"""
mutable struct ElectronOutput{T,F,I<:InteractOutput{T,F}} <: AbstractInteractOutput{T,F}
    interface::I
    window::Blink.AtomShell.Window
end
function ElectronOutput(; kwargs...)
    interface = InteractOutput(; kwargs...)
    window = _newelectronwindow(interface)
    ElectronOutput(interface, window)
end

interface(o::ElectronOutput) = o.interface

# Forward output methods to InteractOutput: ElectronOutput is just a wrapper.
for f in (:length, :size, :firstindex, :lastindex)
    @eval $f(o::ElectronOutput) = $f(interface(o))
end

for f in (:push!, :append!)
    @eval $f(o::ElectronOutput, x) = $f(interface(o), x)
end

for f in (:getindex, :view, :setindex!)
    @eval $f(o::ElectronOutput, i::Union{Colon,AbstractArray}) = $f(interface(o), i)
    @eval $f(o::ElectronOutput, i::Int) = $f(interface(o), i)
end

# DynamicGrids.jl interface
for f in (
    :frames, :init, :aux, :mask, :extent, :graphicconfig, :imageconfig,
    :stoppedframe, :tspan, :fps, :isasync, :isstored, :store, :isshowable,
    :minval, :maxval, :renderer
)
   @eval DG.$f(o::ElectronOutput) = DG.$f(interface(o))
end

for f in (
    :setfps!, :setrunning!, :settspan!, :settimestamp!, :setstoppedframe!,
    :initialise!, :finalise!, :initialisegraphics, :finalisegraphics, :maybesleep,
)
   @eval DG.$f(o::ElectronOutput, args...) = DG.$f(interface(o), args...)
end

DG.storeframe!(o::ElectronOutput, data::DG.AbstractSimData) =
    DG.storeframe!(interface(o), data)
DG.showframe(o::ElectronOutput, data::DG.AbstractSimData, args...) =
    DG.showframe(interface(o), data, args...)

# Running checks depend on the blink window still being open
DG.isrunning(o::ElectronOutput) = _isalive(o) && DG.isrunning(interface(o))

_isalive(o::ElectronOutput) = o.window.content.sock.state == WebSockets.ReadyState(1)

# Display

function Base.display(o::ElectronOutput)
    if !_isalive(o)
        o.window = _newelectronwindow(interface(o))
    end
    return nothing
end

function _newelectronwindow(interface)
    window = Blink.AtomShell.Window()
    body!(window, interface.page)
    return window
end