text:
@testset "Pairwise Dispatch" begin
    N = 100
    points = Vortex.Point.(rand(ComplexF64, N), rand(N))
    sheet  = Vortex.Sheet(rand(ComplexF64, N), accumulate(+, rand(N)), rand())
    blobs = sheet.blobs

    z = rand(ComplexF64)

    wp = induce_velocity(z, points, 0.0)
    wb = induce_velocity(z, blobs, 0.0)
    ws = induce_velocity(z, sheet, 0.0)

    @test wb == ws

    w = induce_velocity(z, (points, blobs, sheet), 0.0)
    @test w == (wp + wb + ws)

    w = induce_velocity(z, ((points, blobs), sheet), 0.0)
    @test w == (wp + wb + ws)

    zs = getfield.(points, :z)
    wps = induce_velocity(zs, points, 0.0)
    wbs = induce_velocity(zs, blobs, 0.0)
    wss = induce_velocity(zs, sheet, 0.0)

    @test wss == wbs
    @test wbs == induce_velocity(points, blobs, 0.0)
    @test wss == induce_velocity(points, sheet, 0.0)

    @test wbs .+ wss == induce_velocity(points, (blobs, sheet), 0.0)

    reset_velocity!(wss)
    induce_velocity!(wss, points, sheet, 0.0)

    @test wss == wbs

    reset_velocity!(wps)
    induce_velocity!(wps, points, (blobs, sheet), 0.0)
    @test wps == wss .+ wbs

    # Distribute to a group of sources
    sys = (sheet, blobs)
    ws = allocate_velocity(sys)
    induce_velocity!(ws, sys, points, 0.0)

    @test ws[1] == induce_velocity(sheet, points, 0.0)
    @test ws[2] == induce_velocity(blobs, points, 0.0)

    # Self-induce velocities
    reset_velocity!(ws)
    self_induce_velocity!(ws, sys, 0.0)
    @test ws == self_induce_velocity(sys, 0.0)

    wb = allocate_velocity(blobs)
    induce_velocity!(wb, blobs, blobs, 0.0)
    induce_velocity!(wb, blobs, sheet, 0.0)

    @test ws[1] ≈ wb

    reset_velocity!(wb)
    induce_velocity!(wb, sheet, sheet, 0.0)
    induce_velocity!(wb, sheet, blobs, 0.0)

    @test ws[2] ≈ wb

    points = Vortex.Point.(rand(ComplexF64, 2N), rand(2N))
    reset_velocity!(ws, (points, blobs))
    @test length.(ws) == (2N, N)

    reset_velocity!(ws, (points, points))
    reset_velocity!(wb, points)
    self_induce_velocity!(wb, points, 0.0)
    mutually_induce_velocity!(ws[1], ws[2], points, points, 0.0)
    @test ws[1] ≈ ws[2]
    @test ws[1] ≈ wb

    reset_velocity!(wb, points)
    induce_velocity!(wb, points, points, 0.0)
    @test ws[1] ≈ wb

    @testset "defaults" begin
        # Minimum implementation of a vortex point source
        @eval struct NewPoint <: Elements.Element
            z::ComplexF64
            Γ::Float64
        end

        Elements.position(p::NewPoint) = p.z
        Elements.circulation(p::NewPoint) = p.Γ
        Elements.impulse(p::NewPoint) = -im*p.z*p.Γ

        Elements.kind(::NewPoint) = Elements.Singleton
        Elements.kind(::Type{NewPoint}) = Elements.Singleton

        function Motions.induce_velocity(z::ComplexF64, p::NewPoint, t)
            p.Γ*Vortex.Points.cauchy_kernel(z - p.z)
        end

        zs = rand(ComplexF64, 100)
        Γs = rand(100)

        new_points = NewPoint.(zs, Γs)
        points = Vortex.Point.(zs, Γs)

        wn = allocate_velocity(new_points)
        wp = allocate_velocity(points)

        self_induce_velocity!(wn, new_points, 0.0)
        self_induce_velocity!(wp, points, 0.0)

        @test wn ≈ wp
    end

    @testset "Property Macro" begin
        import PotentialFlow.Properties: @property

        @eval Elements (@property begin
                      signature = max_blob_radius(s::Source)
                      stype = Int
                      reduce = max
                      end)

        @eval Elements (@property begin
                      signature = induce_count(t::Target, tc::Target, s::Source, sc::Source)
                      preallocator = allocate_count
                      stype = Int
                      end)

        @eval Vortex.Points begin
            Elements.induce_count(::ComplexF64, tcount, ::Point, scount) = scount
            Elements.max_blob_radius(::Point) = 0
        end
        @eval Vortex.Blobs begin
            Elements.induce_count(::ComplexF64, tcount, ::Blob, scount) = tcount
            Elements.max_blob_radius(b::Blob) = b.δ
        end

        N = rand(1:100)
        points = Vortex.Point.(rand(ComplexF64, N), rand(N))
        pcounts = fill(1, length(points))

        blobs  = Vortex.Blob.(rand(ComplexF64, N), rand(N), rand())
        bcounts = fill(2, length(blobs))

        sources  = Source.Blob.(rand(ComplexF64, N), rand(N), rand())
        bcounts = fill(2, length(blobs))

        targets = rand(ComplexF64, rand(1:100))
        tcounts = fill(3, length(targets))

        @test Elements.max_blob_radius(points) == 0
        @test Elements.max_blob_radius(blobs) == blobs[1].δ
        @test Elements.max_blob_radius((points, blobs)) == blobs[1].δ

        @test Elements.induce_count(targets[1], tcounts[1], points, pcounts) == length(points)
        @test Elements.induce_count(targets[1], tcounts[1], blobs, bcounts) == 3length(blobs)
        @test Elements.induce_count(targets[1], tcounts[1],
                                  (points, blobs), (pcounts, bcounts)) == 3length(blobs) + length(points)

        out = Elements.induce_count(targets, tcounts, (points, blobs), (pcounts, bcounts))
        @test length(out) == length(targets)

        out₂ = Elements.allocate_count(targets)
        Elements.induce_count!(out₂, targets, tcounts, (points, blobs), (pcounts, bcounts))

        @test out == out₂
        @test all(out .== 3length(blobs) + length(points))
    end
end


text:
_preprocess_args(p, args, s) = args #needs to modify still_to_process
preprocessArgs!(p) = p
is_seriestype_supported(st) = true
finalize_subplot!(plt, att) = nothing

function _process_userrecipes(plt, plotattributes::AbstractDict{Symbol,Any}, args)
    still_to_process = RecipesBase.RecipeData[]
    args = _preprocess_args(plotattributes, args, still_to_process)

    # for plotting recipes, swap out the args and update the parameter dictionary
    # we are keeping a stack of series that still need to be processed.
    # each pass through the loop, we pop one off and apply the recipe.
    # the recipe will return a list a Series objects... the ones that are
    # finished (no more args) get added to the kw_list, the ones that are not
    # are placed on top of the stack and are then processed further.
    kw_list = Dict{Symbol,Any}[]
    while !isempty(still_to_process)
        # grab the first in line to be processed and either add it to the kw_list or
        # pass it through apply_recipe to generate a list of RecipeData objects (data + attributes)
        # for further processing.
        next_series = popfirst!(still_to_process)
        # recipedata should be of type RecipeData.  if it's not then the inputs must not have been fully processed by recipes
        if !(typeof(next_series) <: RecipesBase.RecipeData)
            error("Inputs couldn't be processed... expected RecipeData but got: $next_series")
        end
        if isempty(next_series.args)
            _process_userrecipe(plt, kw_list, next_series)
        else
            rd_list = RecipesBase.apply_recipe(next_series.plotattributes, next_series.args...)
            prepend!(still_to_process,rd_list)
        end
    end

    # don't allow something else to handle it
    plotattributes[:smooth] = false
    kw_list
end

# plot recipes

# Grab the first in line to be processed and pass it through apply_recipe
# to generate a list of RecipeData objects (data + attributes).
# If we applied a "plot recipe" without error, then add the returned datalist's KWs,
# otherwise we just add the original KW.
function _process_plotrecipe(plt, kw::AbstractDict{Symbol,Any}, kw_list::Vector{Dict{Symbol,Any}}, still_to_process::Vector{Dict{Symbol,Any}})
    if !isa(get(kw, :seriestype, nothing), Symbol)
        # seriestype was never set, or it's not a Symbol, so it can't be a plot recipe
        push!(kw_list, kw)
        return
    end
    try
        st = kw[:seriestype]
        st = kw[:seriestype] = get(_typeAliases, st, st) #TODO this requires access to the const_typeAliases
        datalist = RecipesBase.apply_recipe(kw, Val{st}, plt)
        for data in datalist
            preprocessArgs!(data.plotattributes)
            if data.plotattributes[:seriestype] == st
                error("Plot recipe $st returned the same seriestype: $(data.plotattributes)")
            end
            push!(still_to_process, data.plotattributes)
        end
    catch err
        if isa(err, MethodError)
            push!(kw_list, kw)
        else
            rethrow()
        end
    end
    return
end

# -------------------------------------------------------------------------------

# this method recursively applies series recipes when the seriestype is not supported
# natively by the backend
function _process_seriesrecipe(plt, plotattributes::AbstractDict{Symbol,Any})
    #println("process $(typeof(plotattributes))")
    # replace seriestype aliases
    st = Symbol(plotattributes[:seriestype])
    st = plotattributes[:seriestype] = get(_typeAliases, st, st) #TODO here again

    # shapes shouldn't have fillrange set
    if plotattributes[:seriestype] == :shape
        plotattributes[:fillrange] = nothing
    end

    # if it's natively supported, finalize processing and pass along to the backend, otherwise recurse
    if is_seriestype_supported(st)
        finalize_subplot!(plt, plotattributes)

    else
        # get a sub list of series for this seriestype
        datalist = RecipesBase.apply_recipe(plotattributes, Val{st}, plotattributes[:x], plotattributes[:y], plotattributes[:z])

        # assuming there was no error, recursively apply the series recipes
        for data in datalist
            if isa(data, RecipesBase.RecipeData)
                preprocessArgs!(data.plotattributes)
                if data.plotattributes[:seriestype] == st
                    error("The seriestype didn't change in series recipe $st.  This will cause a StackOverflow.")
                end
                _process_seriesrecipe(plt, data.plotattributes)
            else
                @warn("Unhandled recipe: $(data)")
                break
            end
        end
    end
    nothing
end

ext: jl


text:
struct ClusteringResult{C<:AbstractMatrix{<:AbstractFloat},D<:Real}
    centers::C                 # cluster centers (d x k)
    assignments::Vector{Int}   # assignments (n)
    costs::Vector{D}           # cost of the assignments (n)
    counts::Vector{Int}        # number of points assigned to each cluster (k)
    sd::C                      # standard deviation of each cluster (k) in each dimension (d)
    iter::Int                  # number of iterations required to acquire result
end

const _default_θn = Float64(0.15)
const _default_θe = Float64(1)
const _default_θc = Float64(0.5)
const _default_L = Integer(2)
const _default_iter = Int(10)
const _default_threshold = Float64(0.98)

"""
    base(X, k, [...]) -> ClusteringResult

    ISODATA clustering of the ``d×n`` data matrix `X` (each column of `X`
    is a ``d``-dimensional data point) into maximum of `k` clusters.

    # Arguments
    - `θn`:     minimum number of cluster members (as percentage)
                    *determines discarding clusters
    - `θe`:     maximum standard deviation allowed for a cluster
                    *determines splitting clusters
    - `θc`:     minimum distance between clusters,
                    *determines lumping clusters
    - `L`:      number of clusters to merge during lumping
    - `iter`:   number of iterations to run
"""
function base(
    X::AbstractMatrix{<:Real},
    k::Integer;
    θn::Float64 = _default_θn,
    θe::Float64 = _default_θe,
    θc::Float64 = _default_θc,
    L::Integer = 2,
    iter::Integer = _default_iter,
)
    d, n = size(X)

    # assign initial cluster centers
    # centers = kmppcenters(X, k)
    centers = randomcenters(X, k)
    assignments = Vector{Int}(undef, n)

    dist = pairwise(Euclidean(), centers, X, dims = 2)
    D = eltype(dist)

    counts = Vector{Int}(undef, k)
    costs = Vector{D}(undef, n)
    sd = Array{D}(undef, d, k)

    for i = 1:iter
        # compute number of clusters
        NROWS = size(centers, 2)
        #println(NROWS)

        # assign cluster members
        counts = Vector{Int}(undef, NROWS)
        dist = pairwise(Euclidean(), centers, X, dims = 2)
        update_assignments!(dist, assignments, costs, counts)
        update_centers!(X, assignments, centers, counts)

        dist = pairwise(Euclidean(), centers, X, dims = 2)
        AVEDIST = Vector{Float64}(undef, NROWS)
        AD = Float64(0)
        AD = update_avedist!(AVEDIST, AD, assignments, costs, counts)

        # compute std dev of each component in each cluster
        sd = Array{D}(undef, d, NROWS)
        update_sd!(X, assignments, centers, counts, sd)

        # discard small clusters
        centers = discard_clusters!(X, centers, counts, θn)

        # split or lump
        if size(centers, 2) == NROWS
            if mod(i, 2) != 0 || i == iter
                centers = lump_clusters!(centers, counts, L, θc)
            elseif mod(i, 2) == 0
                centers = split_clusters!(X, centers, counts, sd, k, AVEDIST, AD, θn, θe)
            end
        end
    end

    return ClusteringResult(centers, assignments, costs, counts, sd, iter)
end

"""
    base_conv(X, k, [...]) -> ClusteringResult

    ISODATA clustering of the ``d×n`` data matrix `X` (each column of `X`
    is a ``d``-dimensional data point) into maximum of `k` clusters. Rather
    than taking a parameter for the number of iterations, this function will
    iterate until convergence.

    # Arguments
    - `θn`:     minimum number of cluster members (as percentage)
                    *determines discarding clusters
    - `θe`:     maximum standard deviation allowed for a cluster
                    *determines splitting clusters
    - `θc`:     minimum distance between clusters
                    *determines lumping clusters
    - `L`:      number of clusters to merge during lumping
    - `T`:      maximum number of samples whose cluster assignments
                    is allowed to be unchanged between iterations.
                    *terminates when this threshold is reached
    - `iseeds`: function to determine initial seeds. accepted values:
                    randomcenters, kmppcenters, isodata_s
"""
function base_conv(
    X::AbstractMatrix{<:Real},
    k::Integer;
    θn::Float64 = _default_θn,
    θe::Float64 = _default_θe,
    θc::Float64 = _default_θc,
    L::Integer = 2,
    T::Float64 = _default_threshold,
    iseeds = randomcenters
)
    d, n = size(X)

    # assign initial cluster centers
    centers = iseeds(X, k)
    assignments = Vector{Int}(undef, n)

    dist = pairwise(Euclidean(), centers, X, dims = 2)
    D = eltype(dist)

    costs = Vector{D}(undef, n)
    counts = Vector{Int}(undef, k)
    sd = Array{D}(undef, d, k)
    iter = 0

    while true
        iter += 1
        if iter > 100
            return ClusteringResult(centers, assignments, costs, counts, sd, iter)
        end

        # compute number of clusters
        NROWS = size(centers, 2)

        # assign cluster members
        counts = Vector{Int}(undef, NROWS)
        dist = pairwise(Euclidean(), centers, X, dims = 2)

        prev_assignments = deepcopy(assignments)
        update_assignments!(dist, assignments, costs, counts)
        update_centers!(X, assignments, centers, counts)

        # check convergence threshold
        updates = assignments .== prev_assignments
        if (sum(updates)/length(updates)) > T
            return ClusteringResult(centers, assignments, costs, counts, sd, iter)
        end

        dist = pairwise(Euclidean(), centers, X, dims = 2)
        AVEDIST = Vector{Float64}(undef, NROWS)
        AD = Float64(0)
        AD = update_avedist!(AVEDIST, AD, assignments, costs, counts)

        # compute std dev of each component in each cluster
        sd = Array{D}(undef, d, NROWS)
        update_sd!(X, assignments, centers, counts, sd)

        # discard small clusters
        centers, deleted = discard_clusters!(X, centers, counts, θn)

        # split or lump
        if !deleted && size(centers, 2) == NROWS
            if mod(iter, 2) != 0
                centers = lump_clusters!(centers, counts, L, θc)
            elseif mod(iter, 2) == 0
                centers = split_clusters!(X, centers, counts, sd, k, AVEDIST, AD, θn, θe)
            end
        end
    end

    return ClusteringResult(centers, assignments, costs, counts, sd, iter)
end

"""
    update_assignments!(dist, assignments, costs, counts) -> None

    Update `assignments` in place using `dist`. Also updates `costs`
    and `counts`.

    # Arguments
    -        `dist`: array of the distance of each point to each cluster
    - `assignments`: vector of length `n` with the assigned cluster
    -       `costs`: vector with the cost to assign data point `i` to
                        the nearest cluster center
    -      `counts`: array of the number of patterns assigned to each cluster
"""
function update_assignments!(
    dist::Matrix{<:Real},       # in
    assignments::Vector{Int},   # out
    costs::Vector{<:Real},      # out
    counts::Vector{Int},
)        # out
    k, n = size(dist)

    fill!(counts, 0)
    fill!(costs, 0)

    for j = 1:n
        c, a = findmin(view(dist, :, j))
        assignments[j] = a
        costs[j] += c
        counts[a] += 1
    end
end

"""
    update_centers!(X, assignments, centers, counts) -> None

    Calculate new means and update `centers` in place using `assignments`.

    # Arguments
    -           `X`: array of the input data
    - `assignments`: vector of length `n` with the assigned cluster
    -     `centers`: array of cluster centers
    -      `counts`: array of the number of patterns assigned to each cluster
"""
function update_centers!(
    X::AbstractMatrix{<:Real},                  # in
    assignments::Vector{Int},                   # in
    centers::AbstractMatrix{<:AbstractFloat},   # out
    counts::Vector{Int},                        # in
)
    d, n = size(X)
    k = size(centers, 2)

    fill!(centers, 0)

    # for each data point
    for j = 1:n
        cj = assignments[j]
        # for each dimension
        for i = 1:d
            centers[i, cj] += X[i, j]
        end
    end

    # for each cluster
    for j = 1:k
        # for each dimension
        for i = 1:d
            centers[i, j] /= counts[j]
        end
    end
end

"""
    update_avedist!(AVEDIST, AD, assignments, costs, counts) -> Float64

    Calculate the average distance (AVEDIST) of each cluster to each of it's member
    patterns and the overall average distance (AD) of clusters to their members.

    # Arguments
    -     `AVEDIST`: vector with the average distance of a cluster center to each of its members
    -          `AD`: the overall average distance of cluster centers to their members
    - `assignments`: vector of length `n` with the assigned cluster
    -       `costs`: vector with the cost to assign data point `i` to
                        the nearest cluster center
    -      `counts`: array of the number of patterns assigned to each cluster
"""
function update_avedist!(
    AVEDIST::Vector{Float64},   # in/out
    AD::Float64,                # out
    assignments::Vector{Int},   # in
    costs::Vector{<:Real},      # in
    counts::Vector{Int},
)        # in
    n = size(costs, 1)
    k = size(counts, 1)

    for j = 1:k
        cost = 0
        for i = 1:n
            if assignments[i] == j
                cost += costs[i]
            end
        end
        AVEDIST[j] = cost / counts[j]
        AD += AVEDIST[j]
    end

    AD /= k
    return AD
end

"""
    update_sd!(X, assignments, centers, counts, sd) -> None

    Calculate new means and update `centers` in place using `assignments`.

    # Arguments
    -           `X`: array of the distance of each point to each cluster
    - `assignments`: vector of length `n` with the assigned cluster
    -     `centers`: array of cluster centers
    -      `counts`: array of the number of patterns assigned to each cluster
    -          `sd`: array of the standard deviation of each cluster
                        in each dimension
"""
function update_sd!(
    X::AbstractMatrix{<:Real},
    assignments::Vector{Int},
    centers::AbstractMatrix{<:AbstractFloat},
    counts::Vector{Int},
    sd::AbstractMatrix{<:AbstractFloat},
)
    d, n = size(X)
    k = size(centers, 2)

    fill!(sd, 0)

    # for each data point
    for j = 1:n
        cj = assignments[j]
        # for each dimension
        for i = 1:d
            sd[i, cj] += (X[i, j] - centers[cj])^2
        end
    end

    # for each cluster
    for j = 1:k
        # for each dimension
        for i = 1:d
            sd[i, j] = sqrt(sd[i, j] / counts[j])
        end
    end
end


text:
using TextGrids
using Test

@testset "TextGrids.jl" begin
    # Write your tests here.
end


text:
# Quotes
# module Quotes # not a module anymore
type Quote
  value::Float64
  is_valid::Bool

  Quote(value::Float64, is_valid::Bool = true) = new(value, is_valid)
end

# end


text:
const symbol_numbers = Dict{Symbol, Int}()

"""Return a new, unique symbol like _z3_"""
function make_symbol(s::Symbol)  # default is :z

    i = get(symbol_numbers, s, 0)
    symbol_numbers[s] = i + 1

    if i == 0
        return Symbol("_", s)
    else
        return Symbol("_", s, i)
    end
end

make_symbol(c::Char) = make_symbol(Symbol(c))

let current_symbol = 'a'
    global function make_symbol()
        current_sym = current_symbol

        if current_sym < 'z'
            current_symbol += 1
        else
            current_symbol = 'a'
        end

        return make_symbol(current_sym)

    end
end

make_symbols(n::Integer) = [make_symbol() for i in 1:n]

make_symbols(v::Vector{Symbol}) = make_symbols(length(v))

# The following function is not used
"""Check if a symbol like `:a` has been uniqued to `:_a_1_`"""
function isuniqued(s::Symbol)
    ss = string(s)
    contains(ss, "_") && isdigit(ss[end-1])
end

# Types for representing a flattened AST:

# Combine Assignment and FunctionAssignment ?

struct Assignment
    lhs
    op
    args
end

struct FunctionAssignment
    f  # function name
    args  # input arguments
    return_arguments
    intermediate  # tuple of intermediate variables
end

# Close to single assignment form
mutable struct FlatAST
    top  # topmost variable(s)
    input_variables::Vector{Symbol}
    intermediate::Vector{Symbol}  # generated vars
    code # ::Vector{Assignment}
    variables::Vector{Symbol}  # cleaned version of input_variables
end

function Base.show(io::IO, flatAST::FlatAST)
    println(io, "top: ", flatAST.top)
    println(io, "input vars: ", flatAST.input_variables)
    println(io, "intermediate vars: ", flatAST.intermediate)
    println(io, "code: ", flatAST.code)
end

FlatAST() = FlatAST([], [], [], [], [])

export FlatAST

##

set_top!(flatAST::FlatAST, vars) = flatAST.top = vars  # also returns vars

add_variable!(flatAST::FlatAST, var) = push!(flatAST.input_variables, var)

add_intermediate!(flatAST::FlatAST, var::Symbol) = push!(flatAST.intermediate, var)
add_intermediate!(flatAST::FlatAST, vars::Vector{Symbol}) = append!(flatAST.intermediate, vars)

add_code!(flatAST::FlatAST, code) = push!(flatAST.code, code)

export flatten

function flatten(ex, var = [])
    ex = MacroTools.striplines(ex)
    flatAST = FlatAST()
    if !isempty(var)
        for i in var push!(flatAST.input_variables, i) end
    end
    top = flatten!(flatAST, ex, var)

    return top, flatAST
end

"""`flatten!` recursively converts a Julia expression into a "flat" (one-dimensional)
structure, stored in a FlatAST object. This is close to SSA (single-assignment form,
https://en.wikipedia.org/wiki/Static_single_assignment_form).

Variables that are found are considered `input_variables`.
Generated variables introduced at intermediate nodes are stored in
`intermediate`.
Returns the variable at the top of the current piece of the tree."""

# TODO: Parameters

# numbers:
function flatten!(flatAST::FlatAST, ex::ModelingToolkit.Constant, var)
    return ex.value  # nothing to do the AST; return the number
end

function flatten!(flatAST::FlatAST, ex, var)
    return ex  # nothing to do to the AST; return the number
end

# symbols:
function flatten!(flatAST::FlatAST, ex::Variable, var)  # symbols are leaves
        if isempty(var)
            add_variable!(flatAST, Symbol(ex))  # add the discovered symbol as an input variable
        end
    return Symbol(ex)
end

function flatten!(flatAST::FlatAST, ex::Symbol, var)
    if isempty(var)
        add_variable!(flatAST, ex)  # add the discovered symbol as an input variable
    end
   return ex
end

function flatten!(flatAST::FlatAST, ex::Expr, var = [])
    local top

    if ex.head == :$    # constants written as $a
        top = process_constant!(flatAST, ex)

    elseif ex.head == :call  # function calls
        top = process_call!(flatAST, ex, var)

    elseif ex.head == :(=)  # assignments
        top = process_assignment!(flatAST, ex)

    elseif ex.head == :block
        top = process_block!(flatAST, ex)

    elseif ex.head == :tuple
        top = process_tuple!(flatAST, ex)

    elseif ex.head == :return
        top = process_return!(flatAST, ex)

    else
        error("Currently unable to process expressions with ex.head=$(ex.head)")

    end

    set_top!(flatAST, top)
end

function flatten!(flatAST::FlatAST, ex::Operation, var)
    # top = process_operation!(flatAST, ex, var)
    # set_top!(flatAST, top)

    if ex.op isa Variable
        return flatten!(flatAST, ex.op, var)
    else
       top = process_operation!(flatAST, ex, var)
       set_top!(flatAST, top)
    end
end

function process_constant!(flatAST::FlatAST, ex)
    return esc(ex.args[1])  # interpolate the value of the external constant
end

"""A block represents a linear sequence of Julia statements.
They are processed in order.
"""

function process_block!(flatAST::FlatAST, ex)
    local top

    for arg in ex.args
        isa(arg, LineNumberNode) && continue
        top = flatten!(flatAST, arg)
    end

    return top  # last variable assigned
end

# function process_iterated_function!(flatAST::FlatAST, ex)

function process_tuple!(flatAST::FlatAST, ex)
    # println("Entering process_tuple")
    # @show flatAST
    # @show ex

    top_args = [flatten!(flatAST, arg) for arg in ex.args]

    # top_args = []  # the arguments returned for each element of the tuple
    # for arg in ex.args
    #     top = flatten!(flatAST, arg)
    #     push!(top_args, top)
    # end

    return top_args

end

"""An assigment is of the form a = f(...).
The name a is currently retained.
TODO: It should later be made unique.
"""
function process_assignment!(flatAST::FlatAST, ex)
    # println("process_assignment!:")
    #  @show ex
    #  @show ex.args[1], ex.args[2]

    top = flatten!(flatAST, ex.args[2])
    # @show top

    var = ex.args[1]
    # @show var

    # TODO: Replace the following by multiple dispatch
    if isa(var, Expr) && var.head == :tuple
        vars = [var.args...]

    elseif isa(var, Tuple)
        vars = [var...]

    elseif isa(var, Vector)
        vars = var

    else
        vars = [var]
    end

    add_intermediate!(flatAST, vars)

    top_level_code = Assignment(vars, :(), top)  # empty operation
    add_code!(flatAST, top_level_code)

    # j@show flatAST

    return var

end

"""Processes something of the form `(f↑4)(x)` (write as `\\uparrow<TAB>`)
by rewriting it to the equivalent set of iterated functions"""
function process_iterated_function!(flatAST::FlatAST, ex)
    total_function_call = ex.args[1]
    args = ex.args[2:end]

    # @show args

    function_name = total_function_call.args[2]
    power = total_function_call.args[3]  # assumed integer

    new_expr = :($function_name($(args...)))

    # @show new_expr

    for i in 2:power
        new_expr = :($function_name($new_expr))
    end

    # @show new_expr

    flatten!(flatAST, new_expr)  # replace the current expression with the new one
end

"""A call is something like +(x, y).
A new variable is introduced for the result; its name can be specified
    using the new_var optional argument. If none is given, then a new, generated
    name is used.
"""
function process_call!(flatAST::FlatAST, ex, var = [], new_var=nothing)

    #println("Entering process_call!")
    #@show ex
    #@show flatAST
    #@show new_var
    op = ex.args[1]
    #@show op

    if isa(op, Expr)
        if op.head == :line
            return

        elseif op.head == :call && op.args[1]==:↑   # iterated function like f ↑ 4
            return process_iterated_function!(flatAST, ex)
        end

    end

    # rewrite +(a,b,c) as +(a,+(b,c)) by recursive splitting
    # TODO: Use @match here!

    if op in (:+, :*) && length(ex.args) > 3
        return flatten!(flatAST, :( ($op)($(ex.args[2]), ($op)($(ex.args[3:end]...) )) ), var)
    end

    top_args = []
    for arg in ex.args[2:end]

        isa(arg, LineNumberNode) && continue
        top = flatten!(flatAST, arg, var)

        if isa(top, Vector)  # TODO: make top always a Vector?
            append!(top_args, top)

        else
            push!(top_args, top)
        end
    end

    top_level_code = quote end

    #@show op

    if op ∈ keys(reverse_operations)  # standard operator
        if new_var == nothing
            new_var = make_symbol()
        end

        add_intermediate!(flatAST, new_var)

        top_level_code = Assignment(new_var, op, top_args)

    else
        if haskey(registered_functions, op)

            f = registered_functions[op]

            # make enough new variables for all the returned arguments:
            return_args = make_symbols(f.return_arguments)

            intermediate = make_symbols(f.intermediate) #registered_functions[op].intermediate  # make_symbol(:z_tuple)

            add_intermediate!(flatAST, return_args)
            add_intermediate!(flatAST, intermediate)

            top_level_code = FunctionAssignment(op, top_args, return_args, intermediate)

            new_var = return_args

        else

            throw(ArgumentError("Function $op not available. Use @function to define it."))
        end
    end

    add_code!(flatAST, top_level_code)

    return new_var

end

function process_operation!(flatAST::FlatAST, ex, var, new_var=nothing)

    # println("\n\n--")
    # @show flatAST
    # @show ex
    # @show var
    # @show new_var

    op = ex.op

    if op in (+, *) && length(ex.args) > 2
        return flatten!(flatAST, Expression( (op)((ex.args[1]), (op)((ex.args[2:end]...) )) ), var)
    end

    top_args = []
    for arg in ex.args[1:end]

        isa(arg, LineNumberNode) && continue
        top = flatten!(flatAST, arg, var)

        if isa(top, Vector)
            append!(top_args, top)

        else
            push!(top_args, top)
        end
    end

    top_level_code = quote end

    #@show op

    if Symbol(op) ∈ keys(reverse_operations)  # standard operator
        if new_var == nothing
            new_var = make_symbol()
        end

        add_intermediate!(flatAST, new_var)

        top_level_code = Assignment(new_var, Symbol(op), top_args)

    else
        if haskey(registered_functions, Symbol(op))

            f = registered_functions[Symbol(op)]

            # make enough new variables for all the returned arguments:
            return_args = make_symbols(f.return_arguments)

            intermediate = make_symbols(f.intermediate) #registered_functions[op].intermediate  # make_symbol(:z_tuple)

            add_intermediate!(flatAST, return_args)
            add_intermediate!(flatAST, intermediate)

            top_level_code = FunctionAssignment(Symbol(op), top_args, return_args, intermediate)

            new_var = return_args

        else

            throw(ArgumentError("Function $op not available. Use @function to define it."))
        end
    end

    add_code!(flatAST, top_level_code)
    return new_var

end


text:
module EMLUserTray
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("UserTray/User_trayBasic.jl")
	include("UserTray/User_tray.jl")
end


text:
abstract ClassicalSimulator{T} <: Simulator{T}
abstract Classical2dSimulator{T} <: ClassicalSimulator{T}
abstract Classical3dSimulator{T} <: ClassicalSimulator{T}

function initializeCount!(simulator::ClassicalSimulator)
	simulator.countiteration = 0
end

function getSpaceCount(simulator::ClassicalSimulator)
	return getCount(simulator)
end

function getTrotterCount(simulator::ClassicalSimulator)
	return 0
end

function getFirstTerm(simulator::ClassicalSimulator)
	return getTemperature(simulator)
end


text:
"""
    log_histogram(logger, name, (bins,weights); step=step(logger))

Logs a histogram under the tag `name` on given `step`. The histogram must be
passed as a tuple holding the `N+1` bin edges and the height of the `N` bins.

You can also pass the raw data, and a binning algorithm from `StatsBase.jl` will
be used to bin the data.
"""
function log_histogram(logger::TBLogger, name::AbstractString, (bins,weights)::Tuple{AbstractVector, AbstractArray};
                       step=nothing)
    weights = collect(vec(weights))
    summ = SummaryCollection(histogram_summary(name, collect(bins), weights))
    write_event(logger.file, make_event(logger, summ, step=step))
end

"""
    log_histogram(logger, name, data::Vector; step=step(logger))

Bins the values found in `data` and logs them as an histogram under the tag
`name`.
"""
function log_histogram(logger::TBLogger, name::AbstractString, data::AbstractArray;
                       step=nothing)
    data = vec(data)
    hvals = fit(Histogram, data)
    summ = SummaryCollection(histogram_summary(name, collect(hvals.edges[1]), hvals.weights))
    write_event(logger.file, make_event(logger, summ, step=step))
end

"""
    log_vector(logger, name, data::Vector; step=step(logger))

Logs the vector found in `data` as an histogram under the name `name`.
"""
function log_vector(logger::TBLogger, name::AbstractString, data::AbstractVector; step=nothing)
    summ = SummaryCollection(histogram_summary(name, collect(0:length(data)),data))
    write_event(logger.file, make_event(logger, summ, step=step))
end

function histogram_summary(name::AbstractString, edges::AbstractVector{<:Number}, hist_vals::AbstractVector{<:Number})
    @assert length(edges) == length(hist_vals)+1

    hp = HistogramProto(min=minimum(edges), max=maximum(edges),
                        bucket_limit=edges[2:end],
                        bucket=hist_vals)
    Summary_Value(tag=name, histo=hp)
end


text:
# This file is a part of ROOTFramework.jl, licensed under the MIT License (MIT).

using Cxx

cxxinclude("TString.h")

import Base: string, show

export TString, TStringPtr
const TString = cxxt"TString"
const TStringPtr = pcpp"TString"

TString() = @cxx TString()
TString(str::AbstractString) = @cxx TString(pointer(str))

string(tstr::TString) = unsafe_string((@cxx tstr->Data()), (@cxx tstr->Length()))
show(io::IO, x::TString) = print(io, "TString(\"$(escape_string(string(x)))\")")


text:
using Random, Test, TestSetExtensions
using LatexLH, StructLH, ModelObjectsLH, ModelParams, CollegeEntry

ce = CollegeEntry;

# The code uses the fact that reshape undoes vec. This is tested here.
function reshape_test()
    @testset ExtendedTestSet "reshape" begin
        rng = MersenneTwister(123);
        sizeV = (4,3,2);
        x = rand(rng, sizeV...);
        v = vec(x);
        @test isapprox(x, reshape(v, sizeV...))
    end
end

function best_available_test()
    @testset ExtendedTestSet "best available" begin
        nc = 4;
        nl = 3;
        avail_clM = falses(nc, nl);
        avail_clM[1,2] = true;
        @test ce.best_available(avail_clM) == 1
        avail_clM[nc, nl] = true;
        @test ce.best_available(avail_clM) == nc
        avail_clM[1 : (nc-1),:] .= true;
        @test ce.best_available(avail_clM) == nc
	end
end

function entry_decisions_test(switches :: AbstractEntrySwitches{F1},
    admissionS, prefShocks :: Bool; takeSubset :: Bool = false) where F1

    @testset ExtendedTestSet "Entry Decisions $(typeof(switches)), $(typeof(admissionS))" begin
        rng = MersenneTwister(12);
        # println(switches);
        @test StructLH.describe(switches) isa Array{String};
        objId = ObjectId(:entryOneStep);
        st = ce.make_test_symbol_table();
        entryS = init_entry_decision(objId, switches, st);
        nc = n_colleges(switches);
        admProbFct = ce.make_test_admprob_fct_logistic(nc);
        # println(entryS, "  Pref shocks: $prefShocks");

        if takeSubset
            idxV = 1 : 3 : n_types(switches);
            subset_types!(entryS, idxV);
            @test n_types(entryS) == length(idxV)
            @test validate_es(entryS.switches)
        end

        J = n_types(switches);
        nl = n_locations(switches)
        vWork_jV, vCollege_jcM = CollegeEntry.values_for_test(rng, J, nc, nl);
        hsGpaPctV = collect(range(0.1, 0.9, length = J));
        rank_jV = vcat(2 : 2 : J, 1 : 2 : J);

        er = entry_decisions(entryS, admissionS, admProbFct,
            vWork_jV, vCollege_jcM, hsGpaPctV, rank_jV; prefShocks = prefShocks);
        @test validate_er(er; validateFracLocal = prefShocks)
        entryProb_jcM = entry_probs_jc(er);
        eValM = expected_values_jl(er);
        # Need interior entry probs. Otherwise adjust values.
        # Unless there are no pref shocks - then probs are always 0 and 1
        if all(entry_probs_j(er) .< 0.2)  &&  prefShocks
            @warn "Entry probs too low:  $(entry_probs_j(er))"
            @test false
        end
        @test any(entry_probs_j(er) .< 0.8)

        # Switching off preference shocks should not affect entry
        if prefShocks
            er2 = entry_decisions(entryS, admissionS, admProbFct,
                vWork_jV, vCollege_jcM, hsGpaPctV, rank_jV; prefShocks = false);
            @test isapprox(enrollment_cl(er2), enrollment_cl(er), rtol = 0.005);
            @test isapprox(frac_local_c(er2), frac_local_c(er), atol = 0.005);
        end

        # # Computing fracLocal across colleges and across types should give the same answer
        # fracLocal = frac_local(er);
        # fracLocal2 = sum(frac_local_c(er) .* enrollment_c(er)) / sum(enrollment_c(er));
        # @test isapprox(fracLocal, fracLocal2)
        # entryMass_jV = entry_probs_j(er, :all) .* type_mass_j(er);
        # fracLocal3 = sum(frac_local_j(er) .* entryMass_jV) / sum(entryMass_jV);
        # @test isapprox(fracLocal3, fracLocal)

        # Entry results methods
        if n_locations(er) == 1
            @test frac_local(er) == 1.0
            @test all(frac_local_j(er) .== 1.0)
            @test all(frac_local_c(er) .== 1.0)
        elseif prefShocks
            @test 0.0 < frac_local(er) < 1.0
            @test all(0.0 .< frac_local_j(er) .< 1.0)
            # Frac local can be 1 when a college is local only
            @test all(0.0 .< frac_local_c(er) .<= 1.0)
        end

        # Check implied properties
        totalMass = sum(type_mass_jl(er));

        # Solving one student at a time should give the same answer IF no colleges are full.
        for j = 1 : J
            if nl == 1
                entryProb_cV, eVal, entryProbBest_clM =
                    ce.entry_decisions_one_student(
                        entryS, admissionS, admProbFct,
                        vWork_jV[j], vCollege_jcM[j,:], hsGpaPctV[j],
                        fill(false, nc);
                        prefShocks = prefShocks);
                @test sum(entryProb_cV) .<= 1.0
                # Not having full colleges should increase value
                @test eVal >= eValM[j]
                # And it should increase entry - unless we have no pref shocks
                if prefShocks
                    @test sum(entryProb_cV) > sum(entryProb_jcM[j,:] .- 1e-5)
                end
                if !any(colleges_full(er))
                    @test isapprox(entryProb_cV, entryProb_jcM[j,:])
                    @test isapprox(eVal, eValM[j])
                end

            else
                # Solver for student in one location
                for l = 1 : nl
                    entryProb_clM, eVal, entryProbBest_clM = ce.entry_decisions_one_student(
                        entryS, admissionS, admProbFct,
                        vWork_jV[j], vCollege_jcM[j,:],
                        hsGpaPctV[j], fill(false, nc, nl), l;
                        prefShocks = prefShocks);
                    @test all(sum(entryProb_clM, dims = 2) .<= 1.0)
                    # Not having full colleges should increase value
                    @test eVal >= eValM[j, l]
                end
            end
        end
    end
end

# function subset_types_test(admissionS)
#     @testset "Subset types" begin
#         J = 8; nc = 3; nl = 4;

#         vWork_jV = collect(range(-0.5, 1.5, length = J));
#         vCollege_jcM = range(-0.2, 1.2, length = J) * range(1.0, 2.0, length = nc)';
#         hsGpaPctV = collect(range(0.1, 0.9, length = J));
#         rank_jV = vcat(2 : 2 : J, 1 : 2 : J);

#         # solve for all
#         e = CollegeEntry.make_test_entry_sequ_multiloc(J, nc, nl, Float64(1e8));
#         er = entry_decisions(e, admissionS,
#             vWork_jV, vCollege_jcM, hsGpaPctV, rank_jV);
#         @test validate_er(er)
#         entryProb_jcM = entry_probs_jc(er);
#         eValM = expected_values_jl(er);

#         # solve for subset; should yield same solution without capacity constraints
#         typeV = collect(2 : 2 : J);
#         e2 = copy_for_subset_of_types(e, typeV);
#         er2 = entry_decisions(e2, admissionS,
#             vWork_jV[typeV], vCollege_jcM[typeV], hsGpaPctV[typeV],
#             1 : length(typeV));

#         @test isapprox(entryProb_jcM[typeV,:], entry_probs_jc(e2))
#         @test isapprox(eValM[typeV,:], expected_values_jl(e2))
# 	end
# end

function sim_entry_one_test(switches :: AbstractEntrySwitches{F1},
    admissionS) where F1

    @testset ExtendedTestSet "Sim entry $(typeof(switches)), $(typeof(admissionS))" begin
        rng = MersenneTwister(32);
        objId = ObjectId(:entryOneStep);
        st = ce.make_test_symbol_table();
        entryS = init_entry_decision(objId, switches, st);
        nc = n_colleges(switches);
        admProbFct = CollegeEntry.make_test_admprob_fct_logistic(nc);
        # println(entryS);

        nSim = Int(1e5);
        J = n_types(switches);
        nl = n_locations(switches)
        vWork_jV, vCollege_jcM = CollegeEntry.values_for_test(rng, J, nc, nl);
        vWork = vWork_jV[1];
        vCollege_cV = vCollege_jcM[1,:];
        endowPct = 0.7;
        full_clM = rand(rng, Bool, nc, nl);
        full_clM[1,1] = false;
        l = nl;

        prob_clM, eVal, entryProbBest_clM = ce.entry_decisions_one_student(
            entryS, admissionS, admProbFct, vWork, vCollege_cV,
            endowPct, full_clM, l);
        # Check that bounded away from 0, 1
        @test sum(prob_clM) < 0.85
        @test sum(prob_clM) > 0.15

        @test CollegeEntry.check_prob_array(prob_clM)
        # If all corners, nothing to be tested
        @test all(prob_clM .< 0.9)

        probSim_clM, eValSim = CollegeEntry.sim_one_student(
            entryS, admissionS, admProbFct, vWork, vCollege_cV,
            endowPct, full_clM, l, nSim, rng);

        probGap = maximum(abs.(probSim_clM .- prob_clM));
        println("Max prob_cl gap:  ", probGap);
        if probGap > 1e-2
            show_matrix(probSim_clM; header = "Simulated");
            show_matrix(prob_clM; header = "Solved");
        end
        @test isapprox(probSim_clM, prob_clM, atol = 1e-2)
        @test isapprox(eValSim, eVal, rtol = 2e-2)
    end
end

# Simultaneous entry (no capacity constraints) with prob entry a function of a single indicator
# function admission_onevar_test(prefShocks)
#     J = 8;
#     nc = 4;
#     @testset "Admissions based on single indicator" begin
#         admissionS = ce.make_test_admissions_onevar(nc);
#         # No capacity constraints
#         switches = make_test_entry_sequ_multiloc(J, nc, nc + 1, 100 * J * nc;
#             localOnlyV = [1]);
#         for prefShocks ∈ [true, false]
#             entry_decisions_test(switches, admissionS, prefShocks);
#         end
#         sim_entry_one_test(switches, admissionS);
#     end
# end

@testset ExtendedTestSet "Entry decisions" begin
    J = 8;
    nc = 4;
    for admissionS in (
        ce.make_test_admissions_cutoff(nc),
        ce.make_test_admissions_onevar(nc)
        )

        if admissionS isa AdmissionsOneVar
            # Currently only works without capacity constraints +++++
            capacities = (:high, );
            capFactor = 100.0;
        else
            capacities = (:high, :low);
            capFactor = 0.3;
        end

        # The last case ensures that some colleges are full
        for switches in test_entry_switches(J, nc; capacities = capacities);
            for prefShocks ∈ [true, false]
                entry_decisions_test(switches, admissionS, prefShocks);
            end
            sim_entry_one_test(switches, admissionS);
        end

        # Test with subsetting
        J = 21;
        switches =
            make_test_entry_sequ_multiloc(J, nc, nc + 1, capFactor * J * nc);
        entry_decisions_test(switches, admissionS, true; takeSubset = true);
    end

    reshape_test()
    # subset_types_test(admissionS)
end

# -----------


text:
## Visualisation of Geometric Benchmarks
using Compose
using Gadfly

function make_point_pairs(lines)
  b = Array(Any, length(lines)-1)
  for i = 1:length(lines) - 1
    j = i + 1
    b[i] = [lines[i][:,1] lines[j][:,1]]
  end
  b
end

pair(x) = x[1],x[2]

function to_lines(points)
  [line([pair(points[:,i]),pair(points[:,i+1])]) for i = 1:(size(points,2)-1)]
end

function make_compose_lines(point_pairs)
  [Compose.line([pair(o[:,1]), pair(o[:,2])]) for o in point_pairs]
end

function get_lines(lines...)
  all_lines = apply(vcat,lines)
  x = map(l->(context(units=UnitBox(0, 0, 10, 10)),
              l,
              linewidth(.5mm),
              stroke(rand_color()),
              fill(nothing)),
          all_lines)
  x
end

function draw_lines(lines...)
  all_lines = apply(vcat,lines)
  x = map(l->(context(units=UnitBox(0, 0, 10, 10)),
              l,
              linewidth(.5mm),
              stroke(rand_color()),
              fill(nothing)),
          all_lines)
  apply(compose,vcat(context(), x))
end
rand_color() = Compose.RGB(rand(),rand(),rand())


text:
using LinearAlgebra: diagind

struct ConstantFriction{T} <: ElectronicFrictionProvider
    ndofs::Int
    γ::T
end

function friction!(model::ConstantFriction, F::AbstractMatrix, ::AbstractMatrix)
    F[diagind(F)] .= model.γ
end

NQCModels.ndofs(model::ConstantFriction) = model.ndofs


text:
const PlainDS = Array{Face{3,Int},1}

struct FaceRing{T}
    v::Int
    t::T # topology
end

"""
    EdgeRing(v,t)

Construct an edge ring iterator at vertex `v` from a given topology `t`.
"""
struct EdgeRing{T}
    v::Int
    t::T # topology
end

"""
    VertexRing(v,t)

Construct a vertex ring iterator at vertex `v` from a given topology `t`.
"""
struct VertexRing{T}
    v::Int
    start::Union{Int,Nothing}
    t::T # topology
end

VertexRing(v::Int,t) = VertexRing(v,nothing,t)
### Perhaps one may also consider to make a trait

### With this pairiterator I now can initialize EdgeRing as iterator from VertexRing.

struct PairIterator
    iter
end

function Base.iterate(iter::PairIterator)
    i1,state1 = iterate(iter.iter)
    i2, state2 = iterate(iter.iter,state1)
    return Face(i1,i2),(state2,i2)
end

function Base.iterate(iter::PairIterator,state)

    if state==nothing
        return nothing
    end

    state1,i1 = state
    step = iterate(iter.iter,state1)

    if step==nothing
        i2,state2 = iterate(iter.iter)
        return Face(i1,i2), nothing
    else
        i2, state2 = step
        return Face(i1,i2),(state2,i2)
    end
end

function Base.collect(iter::Union{FaceRing,VertexRing})
    collection = Int[]
    for i in iter
        push!(collection,i)
    end
    return collection
end

function Base.collect(iter::EdgeRing)
    collection = Face{2,Int}[]
    for i in iter
        push!(collection,i)
    end
    return collection
end

function Base.collect(iter::PairIterator)
    collection = Face{2,Int}[]
    for i in iter
        push!(collection,i)
    end
    return collection
end

### There is a room for unstructured circulators iterators for performance reasons.

function find_other_triangle_edge(v1::Integer,v2::Integer,skip::Integer,t::PlainDS)
    for i in 1:length(t)
        if in(v1,t[i]) & in(v2,t[i]) & !(i==skip)
            return i
        end
    end
    return -1
end

function find_triangle_vertex(v::Integer,t::PlainDS)

    for i in 1:length(t)
        if in(v,t[i])
            return i
        end
    end
    return length(t) + 1 # -1
end


text:
function to_image(image::AbstractMatrix{<: AbstractFloat}, colormap::AbstractVector{<: Colorant}, colorrange)
    return interpolated_getindex.((to_value(colormap),), image, (to_value(colorrange),))
end

"""
    resample(A::AbstractVector, len::Integer)
Resample a vector with linear interpolation to have length `len`
"""
function resample(A::AbstractVector, len::Integer)
    length(A) == len && return A
    return interpolated_getindex.((A,), range(0.0, stop=1.0, length=len))
end

"""
    resampled_colors(attributes::Attributes, levels::Integer)

Resample the color attribute from `attributes`. Resamples `:colormap` if present,
or repeats `:color`.
"""
function resampled_colors(attributes, levels::Integer)
    cols = if haskey(attributes, :color)
        c = get_attribute(attributes, :color)
        c isa AbstractVector ? resample(c, levels) : repeated(c, levels)
    else
        c = get_attribute(attributes, :colormap)
        resample(c, levels)
    end
end

"""
Like `get!(f, dict, key)` but also calls `f` and replaces `key` when the corresponding
value is nothing
"""
function replace_automatic!(f, dict, key)
    haskey(dict, key) || return (dict[key] = f())
    val = dict[key]
    to_value(val) == automatic && return (dict[key] = f())
    val
end

is_unitrange(x) = (false, 0:0)
is_unitrange(x::AbstractRange) = (true, x)
function is_unitrange(x::AbstractVector)
    length(x) < 2 && return false, 0:0
    diff = x[2] - x[1]
    length(x) < 3 && return true, x[1]:x[2]
    last = x[3]
    for elem in drop(x, 3)
        diff2 = elem - last
        diff2 != diff && return false, 0:0
    end
    return true, range(first(x), diff, length(x))
end

function ngrid(x::AbstractVector, y::AbstractVector)
    xgrid = [Float32(x[i]) for i = 1:length(x), j = 1:length(y)]
    ygrid = [Float32(y[j]) for i = 1:length(x), j = 1:length(y)]
    xgrid, ygrid
end

function nan_extrema(array)
    mini, maxi = (Inf, -Inf)
    for elem in array
        isnan(elem) && continue
        mini = min(mini, elem)
        maxi = max(maxi, elem)
    end
    Vec2f0(mini, maxi)
end

function extract_expr(extract_func, dictlike, args)
    if args.head != :tuple
        error("Usage: args need to be a tuple. Found: $args")
    end
    expr = Expr(:block)
    for elem in args.args
        push!(expr.args, :($(esc(elem)) = $(extract_func)($(esc(dictlike)), $(QuoteNode(elem)))))
    end
    push!(expr.args, esc(args))
    expr
end

"""
usage @exctract scene (a, b, c, d)
"""
macro extract(scene, args)
    extract_expr(getindex, scene, args)
end

"""
    @get_attribute scene (a, b, c, d)

This will extract attribute `a`, `b`, `c`, `d` from `scene` and apply the correct attribute
conversions + will extract the value if it's a signal.
It will make those attributes available as variables and return them as a tuple.
So the above is equal to:
will become:
```example
begin
    a = get_attribute(scene, :a)
    b = get_attribute(scene, :b)
    c = get_attribute(scene, :c)
    (a, b, c)
end
```
"""
macro get_attribute(scene, args)
    extract_expr(get_attribute, scene, args)
end

@inline getindex_value(x::Union{Dict,Attributes,AbstractPlot}, key::Symbol) = to_value(x[key])
@inline getindex_value(x, key::Symbol) = to_value(getfield(x, key))

"""
usage @extractvalue scene (a, b, c, d)
will become:
```example
begin
    a = to_value(scene[:a])
    b = to_value(scene[:b])
    c = to_value(scene[:c])
    (a, b, c)
end
```
"""
macro extractvalue(scene, args)
    extract_expr(getindex_value, scene, args)
end

bs_length(x::NativeFont) = 1 # these are our rules, and for what we do, Vecs are usually scalars
bs_length(x::VecTypes) = 1 # these are our rules, and for what we do, Vecs are usually scalars
bs_length(x::AbstractArray) = length(x)
bs_length(x::AbstractString) = length(x)
bs_length(x) = 1

bs_getindex(x::NativeFont, i) = x # these are our rules, and for what we do, Vecs are usually scalars
bs_getindex(x::VecTypes, i) = x # these are our rules, and for what we do, Vecs are usually scalars
bs_getindex(x::AbstractArray, i) = x[i]
bs_getindex(x::AbstractString, i) = x[i]
bs_getindex(x, i) = x

"""
Like broadcast but for foreach. Doesn't care about shape and treats Tuples && StaticVectors as scalars.
"""
function broadcast_foreach(f, args...)
    lengths = bs_length.(args)
    maxlen = maximum(lengths)

    # all non scalars should have same length
    if any(x -> !(x in (0, 1, maxlen)), lengths)
        error("All non scalars need same length, Found lengths for each argument: $lengths, $(typeof.(args))")
    end

    # skip if there's a zero length element (like an empty annotations collection, etc)
    # this differs from standard broadcasting logic in which all non-scalar shapes have to match
    0 in lengths && return

    for i in 1:maxlen
        f(bs_getindex.(args, i)...)
    end
    return
end

"""
    from_dict(::Type{T}, dict)
Creates the type `T` from the fields in dict.
Automatically converts to the correct node types.
"""
function from_dict(::Type{T}, dict) where T
    T(map(fieldnames(T)) do name
        convert(fieldtype(T, name), dict[name])
    end...)
end

same_length_array(array, value::NativeFont) = repeated(value, length(array))
same_length_array(array, value) = repeated(value, length(array))
function same_length_array(arr, value::Vector)
    if length(arr) != length(value)
        error("Array lengths do not match. Found: $(length(arr)) of $(eltype(arr)) but $(length(value)) $(eltype(value))")
    end
    value
end
same_length_array(arr, value, key) = same_length_array(arr, convert_attribute(value, key))

function to_ndim(T::Type{<: VecTypes{N,ET}}, vec::VecTypes{N2}, fillval) where {N,ET,N2}
    T(ntuple(Val(N)) do i
        i > N2 && return ET(fillval)
        @inbounds return vec[i]
    end)
end

dim3(x) = ntuple(i -> x, Val(3))
dim3(x::NTuple{3,Any}) = x

dim2(x) = ntuple(i -> x, Val(2))
dim2(x::NTuple{2,Any}) = x

lerp(a::T, b::T, val::AbstractFloat) where {T} = (a .+ (val * (b .- a)))

function merged_get!(defaults::Function, key, scene, input::Vector{Any})
    return merged_get!(defaults, key, scene, Attributes(input))
end

function merged_get!(defaults::Function, key, scene::SceneLike, input::Attributes)
    d = defaults()
    if haskey(theme(scene), key)
        # we need to merge theme(scene) with the defaults, because it might be an incomplete theme
        # TODO have a mark that says "theme uncomplete" and only then get the defaults
        d = merge!(to_value(theme(scene, key)), d)
    end
    return merge!(input, d)
end

struct Key{K} end
macro key_str(arg)
    :(Key{$(QuoteNode(Symbol(arg)))})
end
Base.broadcastable(x::Key) = (x,)

to_vector(x::AbstractVector, len, T) = convert(Vector{T}, x)
function to_vector(x::AbstractArray, len, T)
    if length(x) in size(x) # assert that just one dim != 1
        to_vector(vec(x), len, T)
    else
        error("Can't convert to a Vector. Please supply a range/vector/interval")
    end
end
function to_vector(x::ClosedInterval, len, T)
    a, b = T.(extrema(x))
    range(a, stop=b, length=len)
end

"""
Returns (N1, N2) with `N1 x N2 == n`. N2 might become 1
"""
function close2square(n::Real)
    # a cannot be greater than the square root of n
    # b cannot be smaller than the square root of n
    # we get the maximum allowed value of a
    amax = floor(Int, sqrt(n));
    if 0 == rem(n, amax)
        # special case where n is a square number
        return (amax, div(n, amax))
    end
    # Get its prime factors of n
    primeFactors  = factor(n);
    # Start with a factor 1 in the list of candidates for a
    candidates = [1]
    for (f, _) in primeFactors
        # Add new candidates which are obtained by multiplying
        # existing candidates with the new prime factor f
        # Set union ensures that duplicate candidates are removed
        candidates = union(candidates, f .* candidates)
        # throw out candidates which are larger than amax
        filter!(x -> x <= amax, candidates)
    end
    # Take the largest factor in the list d
    (candidates[end], div(n, candidates[end]))
end

"""
A colorsampler maps numnber values from a certain range to values of a colormap
```
x = ColorSampler(colormap, (0.0, 1.0))
x[0.5] # returns color at half point of colormap
```
"""
struct ColorSampler{Data <: AbstractArray}
    colormap::Data
    color_range::Tuple{Float64,Float64}
end

function Base.getindex(cs::ColorSampler, value::Number)
    return interpolated_getindex(cs.colormap, value, cs.color_range)
end

# This function was copied from GR.jl,
# written by Josef Heinen.
"""
    peaks([n=49])

Return a nonlinear function on a grid.  Useful for test cases.
"""
function peaks(n=49)
    x = LinRange(-3, 3, n)
    y = LinRange(-3, 3, n)
    3 * (1 .- x').^2 .* exp.(-(x'.^2) .- (y .+ 1).^2) .- 10 * (x' / 5 .- x'.^3 .- y.^5) .* exp.(-x'.^2 .- y.^2) .- 1 / 3 * exp.(-(x' .+ 1).^2 .- y.^2)
end

get_dim(x, ind, dim, size) = get_dim(LinRange(extrema(x)..., size[dim]), ind, dim, size)
get_dim(x::AbstractVector, ind, dim, size) = x[Tuple(ind)[dim]]
get_dim(x::AbstractMatrix, ind, dim, size) = x[ind]

"""
    surface_normals(x, y, z)
Normals for a surface defined on the grid xy
"""
function surface_normals(x, y, z)
    function normal(i)
        i1, imax = CartesianIndex(1, 1), CartesianIndex(size(z))
        ci(x, y) = min(max(i + CartesianIndex(x, y), i1), imax)
        of = (ci(-1, -1), ci(1, -1), ci(-1, 1), ci(1, 1))
        function offsets(off)
            s = size(z)
            return Vec3f0(get_dim(x, off, 1, s), get_dim(y, off, 2, s), z[off])
        end
        return normalize(mapreduce(offsets, +, init=Vec3f0(0), of))
    end
    return vec(map(normal, CartesianIndices(z)))
end

function attribute_names(PlotType)
    # TODO, have all plot types store their attribute names
    return keys(default_theme(nothing, PlotType))
end

"""
    attributes_from(PlotType, plot)

Gets the attributes from plot, that are valid for PlotType
"""
function attributes_from(PlotType, plot)
    result = Attributes()
    for key in attribute_names(PlotType)
        if haskey(plot, key)
            result[key] = plot[key]
        end
    end
    return result
end


text:
# Tutorial on conventional 2D SLAM example
# This tutorial shows how to use some of the commonly used factor types
# This tutorial follows from the ContinuousScalar example from IncrementalInference
using GraffSDK
using GraffSDK.DataHelpers
using ProgressMeter
using UUIDs
using IncrementalInference
using RoME
using Arena.Amphitheatre
using Colors
##

# 1a. Create a Configuration
config = loadGraffConfig();
#Create a hexagonal sessions
config.sessionId = "HexVisDemo"*string(uuid4())[1:6]

@info "Session backlog (queue length) = $(getSessionBacklog())"

# 2. Confirm that the robot already exists, create if it doesn't.
println(" - Creating or retrieving robot '$(config.robotId)'...")
robot = nothing
if isRobotExisting()
    println(" -- Robot '$(config.robotId)' already exists, retrieving it...")
    robot = getRobot();
else
    # Create a new one programatically - can also do this via the UI.
    println(" -- Robot '$(config.robotId)' doesn't exist, creating it...")
    newRobot = RobotRequest(config.robotId, "My New Bot", "Description of my neat robot", "Active");
    robot = addRobot(newRobot);
end
println(robot)

# 3. Create or retrieve the session.
# Get sessions, if it already exists, add to it.
println(" - Creating or retrieving data session '$(config.sessionId)' for robot...")
session = nothing
if isSessionExisting()
    println(" -- Session '$(config.sessionId)' already exists for robot '$(config.robotId)', retrieving it...")
    session = getSession()
else
    # Create a new one
    println(" -- Session '$(config.sessionId)' doesn't exist for robot '$(config.robotId)', creating it...")
    newSessionRequest = SessionDetailsRequest(config.sessionId, "A test dataset demonstrating data ingestion for a wheeled vehicle driving in a hexagon.", "Pose2")
    session = addSession(newSessionRequest)
end
println(session)

# 4. Drive around in a hexagon
addVariable(:x0, Pose2, ["POSE"])
# Add at a fixed location PriorPose2 to pin :x0 to a starting location (10,10, pi/4)
addFactor([:x0], IIF.Prior( MvNormal([0.0; 0.0; 0.0], Matrix(Diagonal([0.1;0.1;0.05].^2)) )))

println(" - Adding hexagonal driving pattern to session...")
@showprogress for i in 0:5
    psym = Symbol("x$i")
    nsym = Symbol("x$(i+1)")
    println(" - Pose $nsym: Adding new odometry measurement...")
    addVariable(nsym, Pose2, ["POSE"])

    pp = Pose2Pose2(MvNormal([1.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    addFactor([psym;nsym], pp)
end

# 5. Now lets add a couple landmarks
# Ref: https://github.com/dehann/RoME.jl/blob/master/examples/Slam2dExample.jl#L35
addVariable(:l1, Point2, ["LANDMARK"])
p2br2 = Pose2Point2BearingRange(Normal(0,0.1),Normal(2.0,0.2))
addFactor([:x0; :l1], p2br2)
addFactor([:x6; :l1], p2br2)

# Landmarks generally require more work once they're created, e.g. creating factors,
# so they are not set to ready by default. Once you've completed all the factor links and want to solve,
# call putReady to tell the solver it can use the new nodes. This is added to the end of the processing queue.
putReady(true)

# 5. Create AbstractAmphitheatre container to hold different visualizers
visdatasets = AbstractAmphitheatre[]
graffVis = BasicGraffPose(config)

push!(visdatasets, graffVis)

vis, vistask = visualize(visdatasets)

# 6. now create a local fg hexslam
# start with an empty factor graph object
fg = initfg()

# Add the first pose :x0
addVariable!(fg, :x0, Pose2, labels=["POSE"])

# Add at a fixed location PriorPose2 to pin :x0 to a starting location (10,10, pi/4)
addFactor!(fg, [:x0], IIF.Prior( MvNormal([4; 0; -pi], Matrix(Diagonal([0.1;0.1;0.05].^2)) )))

#
romeVis = BasicFactorGraphPose("DemoRobot","LocalHexVisDemo"*string(uuid4())[1:6], fg, meanmax=:mean, poseProp = plDrawProp(0.3, 0.1, RGBA(0,1,1,0.5)))

push!(visdatasets, romeVis)

# Drive around in a hexagon
for i in 0:5
  psym = Symbol("x$i")
  nsym = Symbol("x$(i+1)")
  addVariable!(fg, nsym, Pose2, labels=["POSE"])
  pp = Pose2Pose2(MvNormal([1.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
  addFactor!(fg, [psym;nsym], pp )
end

# Add landmarks with Bearing range measurements
addVariable!(fg, :l1, Point2, labels=["LANDMARK"])
    p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(2.0,0.2))
    addFactor!(fg, [:x0; :l1], p2br)

# Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.1),Normal(2.0,0.2))
    addFactor!(fg, [:x6; :l1], p2br2)

# solve
batchsolve = @async batchSolve!(fg)

##
@info "To stop call stopAmphiVis!()"
stopAmphiVis!()


text:
using OrdinaryDiffEq
using Distributions
using DiffEqUncertainty
using Plots

function my_centralmoment(n::Int, g::Function, args...; kwargs...)
    if n < 2 return Float64[] end
    # Compute the expectations of g, g^2, ..., g^n
    sol = expectation(x -> [g(x)^i for i in 1:n], args...; nout = n, kwargs...)
    exp_set = sol[:]
    mu_g = popfirst!(exp_set)
    # Combine according to binomial expansion
    const_term(n) = (-1)^(n-1) * (n-1) * mu_g^n
    binom_term(n, k, mu, exp_gi) = binomial(n, k) * (-mu)^(n - k) * exp_gi
    binom_sum = function (exp_vals)
        m = length(exp_vals) + 1
        sum([binom_term(m, k + 1, mu_g, v) for (k,v) in enumerate(exp_vals)]) + const_term(m)
    end
    return [zero(exp_set[1]), [binom_sum(exp_set[1:i]) for i in 1:length(exp_set)]...]
end

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    @inbounds begin
        du[1] = -β*c*I/N*S
        du[2] = β*c*I/N*S - γ*I
        du[3] = γ*I
    end
    nothing
end;

δt = 1.0
tmax = 40.0
tspan = (0.0,tmax);
t = 0:δt:tmax;

u0 = [990.0,10.0,0.0]
p = [0.05,10,0.25]
prob_ode = ODEProblem(sir_ode!,u0,tspan,p);

p_dist = [Uniform(0.01,0.1),
     Uniform(5,20.0),
     Uniform(0.1,1.0)];

g(sol) = sol(t)
g(sol,x,i) = sol(x)[i];

n_samples = 1000
sol_ode_mean_mc = expectation(g, prob_ode, u0, p_dist, MonteCarlo(), Tsit5(); trajectories = n_samples)
sol_ode_mean_mc = Array(sol_ode_mean_mc)'

sol_ode_std_mc = [[sqrt(my_centralmoment(2, (sol) -> g(sol,x,i), prob_ode, u0, p_dist, MonteCarlo(), Tsit5(); trajectories = n_samples)[2]) for x in t] for i in 1:3]
sol_ode_std_mc = hcat(sol_ode_std_mc...)

sol_ode_mean_k = expectation(g, prob_ode, u0, p_dist, Koopman(), Tsit5())
sol_ode_mean_k = Array(sol_ode_mean_k)'

sol_ode_std_k = [[sqrt(my_centralmoment(2, (sol) -> g(sol,x,i), prob_ode, u0, p_dist, Koopman(), Tsit5())[2]) for x in t] for i in 1:3]
sol_ode_std_k = hcat(sol_ode_std_k...)

l = @layout [a b]
pl1 = plot(t,
     sol_ode_mean_mc,
     ribbon=sol_ode_std_mc,
     fillalpha=0.15,
     label=["S" "I" "R"],
     xlabel="Time",
     ylabel="Number",
     title="Monte Carlo")
pl2 = plot(t,
     sol_ode_mean_k,
     ribbon=sol_ode_std_k,
     fillalpha=0.15,
     label=["S" "I" "R"],
     xlabel="Time",
     ylabel="Number",
     title="Koopman")
plot(pl1,pl2,layout=l)


text:
# bayes_meta_analysis.jl
# functions for performing a Bayesian meta analysis
# Pkg.checkout("Plots","dev")
using Distributions  #,  Plots, StatPlots
# Normalize values for each experiment to allow comparison

function standardize_scores(means,stds,nobs)
  # Normalizing scores to [0,1] interval
  # Dividing by week 0 score

  # if no data, return all zeros
  for i in 1:length(nobs)
    if nobs[i] == 0
      nobs[i] = 1
    end
  end

  stdz_mn = means./means[1]
  # scale SD
  stdz_sd = stds./means[1]
  # standard error of scaled mean
  stdz_se = stdz_sd./sqrt.(nobs)
  return stdz_mn, stdz_se
end

### Cannot divide by std for inference!
#function standardize_scores_sddiv(means,stds,nobs)
  # Normalizing scores to [0,1] interval
  # Dividing by week 0 score
#  stdz_mn = (means[1] .- means) ./ stds
  # scale SD
  #stdz_sd = stds/means[1]
  # standard error of scaled mean
#  stdz_se = stdz_sd./sqrt(nobs)
#  return stdz_mn, stdz_se
#end

"""
        update_mean(m1,m0,s21,s22,n1,n0)

Compute updated (combined) posterior summary statistics given
two sets of posterior statistics assuming both Gaussian.

 - For uninformative prior, set m0 = s20 = n0 = 0.0
 - see Gelman et al. (2014), eqn. (3.7), p. 68

Input:

        m0 = posterior mean for first sample (or prior mean)

        m1 = posterior mean for second sample (or likelihood mean)

        s0 = posterior standard error for sample 1 mean (NOT SD of sample)

        s1 = posterior standard error for sample 1 mean (NOT SD of sample)

        n0 = number of obs. in sample 1

        n1 = number of obs. in sample 2

N.B.: s = std(x)/sqrt(n)

Returns:

        m2 = posterior mean of combined sample

        s2 = posterior standard error for combined sample

        n2 = n0 + n1
"""
function update_mean(m1,m0,s1,s0,n1,n0)
  s21 = s1^2
  s20 = s0^2
  n2 = n1 + n0
  m2 = (n1/n2)*m1 + (n0/n2)*m0
  v2 = n2 - 1.
  vs22 = (n1-1.)*s21 + (n0-1.)*s20 +
              ((n1*n0)/n2)*(m1 - m0)^2.
  s2 = sqrt(vs22/v2)
  return m2, s2, n2
end

"""
        marginal_posterior_mu(m,s,n, M)

  Compute M marginal Student-t posterior draws.

Input:

          posterior mean, m, standard error, s, and number of obs. n

          M = number of draws (optional: default M = 100,000)

Returns:

          M draws from Student-t posterior density.
"""
function marginal_posterior_mu(m,s,n;M=100000)
   df = n-1
   s2 = 1.0./rand(Gamma(df/2, 2/df),M)
   tv = m .+ s.*sqrt.(s2).*rand(Normal(),M)
   return tv
end

"""
draw from posteriors for treatment and nontreatment
mn = mean, s2 = variance, n = no. of obs.
"""
function post_treat_notreat(mntreat,s2treat,ntreat, mnnotreat, s2notreat,nnotreat; M = 10000)
  # compute posteriors for treatment and nontreatment
  # group scores for each week/period
  wks = length(mntreat)
  drawst = zeros(M,wks)
  drawsnt = zeros(M,wks)
  for i in 1:wks
    drawst[:,i] = marginal_posterior_mu(mntreat[i],s2treat[i],ntreat[i], M = M)
    drawsnt[:,i] = marginal_posterior_mu(mnnotreat[i],s2notreat[i],nnotreat[i], M = M)
  end
  return drawst, drawsnt
end

"""
Compute posterior difference in treatment vs. nontreatment means
"""
function treat_pcbo_difference(drawst,drawsnt)
  difdraws = zeros(length(drawst[:,1]),length(drawst[1,:]))
  for i in 1:length(drawst[1,:])
    difdraws[:,i] = drawst[:,i] - drawsnt[:,i]
  end
  return difdraws
end

"""
Trajectory plot with 0.95 CIs given posterior means
"""
function traj_plot(means;overlay = false,nwks=3,ci=true)
  if (nwks == 3)
    z = [1 2 3]'   # for up to week 8 = 3 data points
  else
    z = [1 2 3 4]'   # for up to week 12 = 4 data points
  end
  mn = means
  lz = log.(z)
  Xloglin = hcat(ones(nwks),lz)
  b = Xloglin \ mn[1:nwks]
  s2 = sum((mn[1:nwks] - Xloglin*b).^2)/(nwks-1.)
  println("coeffs = ",b)
  println("s^2 (eqn. variance) = ",s2)

  # compute R^2
  tss = sum((mn[1:nwks]' - mean(mn[1:nwks])).^2)
  R2 = 1.0 - s2*2.0/tss
  println("Rsquared = ",R2)
  zs = 1.0:0.1:4.2
  score = b[1] + b[2].*log.(zs')

  #plot(1.0:0.1:4.2,score',size=(400,200),label="")
  mns = mn[1:nwks]
  #scatter!(z',mns',label="")

  s = sqrt(s2)
  tcrit5 = 4.3  # 5% crit. value from Student-t with 2 d.f.
  up95 = score .+ s*tcrit5
  lo95 = score .- s*tcrit5
  #  zz = 0.0:0.37:12.2
  zz = 1.0:0.1:4.2
  if (overlay == false)
    labels = [0,4,8,12]
#    zz = 0.0:0.37:12.2
 zz = 1.0:0.1:4.2
    plot(zz,score',size=(400,200),label="")
  else
    plot!(zz,score',size=(400,200),label="")
  end
  wks = [0 4 8 12]
  #scatter!(wks[1:nwks]',mns',label="")
  scatter!(z',mns',label="")

  if ci == true
    plot!(zz,up95',size=(400,200),label="",color=:green)
    plot!(zz,lo95',size=(400,200),label="",color=:green)
  end
#  hline!([1.0], color=:red,label="")
end

function pvals_by_week(m)
# p-value computation:
  df = csv_plot(m)
  x=df[:,1]   # scale horizontal axis to weeks
  wks = convert_to_12weeks(x)
  b,s,r2 = traj_reg_stats(m)
  score = df[:,2]   # predictive means
  bpval = zeros(length(score))
  for i in 1:length(score)
    #tvals = marg_post_mu(score[i],s^2,nwks)
    #bpval[i] = bayespval(tvals)
    t = score[i]/s
    bpval[i] = cdf(TDist(3),t)

  end
  zz = zeros(7,3)
  j = 1
  for i = 1:length(score)
    if (wks[i] == 0.0) || (wks[i] == 2.0) || (wks[i] == 4.0) || (wks[i] == 6.0) || (wks[i] == 8.0) || (wks[i] == 10.0) || (wks[i] == 12.0)
      println([ round(score[i],3),wks[i],bpval[i] ])
      zz[j,1], zz[j,2], zz[j,3] =  round(score[i],3), wks[i], bpval[i]
      j += 1
    end
  end
  return zz
end

function traj_cis(means;nwks=3)
  if (nwks == 3)
    z = [1 2 3]'   # for up to week 8 = 3 data points
  else
    z = [1 2 3 4]'   # for up to week 12 = 4 data points
  end

  mn = means
  lz = log.(z)
  Xloglin = hcat(ones(nwks),lz)
  b = Xloglin \ mn[1:nwks]
  s2 = sum((mn[1:nwks] - Xloglin*b).^2)/(nwks-1.)
  zs =  [1 2 3 4]'
  score = b[1] + b[2].*log.(zs')

  s = sqrt(s2)
  tcrit5 = 4.3  # 5% crit. value from Student-t with 2 d.f.
  up95 = score .+ s*tcrit5
  lo95 = score .- s*tcrit5
  return lo95, up95
end

function cond_post_mu(m,s2; M = 10000)
  # need to either evaluate pdf of t, or
  # draw from t and plot pseudo-sample
  ts = rand(Normal(m,sqrt(s2)),M)
  return ts
end

function traj_reg_stats(means; nwks=3)
  if (nwks == 3)
    z = [1 2 3]'   # for up to week 8 = 3 data points
  else
    z = [1 2 3 4]'   # for up to week 12 = 4 data points
  end

  mn = means
  lz = log.(z)
  Xloglin = hcat(ones(nwks),lz)
  b = Xloglin \ mn[1:nwks]
  s2 = sum((mn[1:nwks] - Xloglin*b).^2)/(nwks-1.)
  println("coeffs = ",b)
  println("s^2 (eqn. variance) = ",s2)

  # compute R^2
  tss = sum((mn[1:nwks]' - mean(mn[1:nwks])).^2)
  R2 = 1 - s2*2/tss
  println("Rsquared = ",R2)
  s = sqrt(s2)
  return b, s, R2
end


text:
import SpecialFunctions: loggamma
using Distributions
using Optim

"""
    ν(x) ∝ exp(p*β*x - p*logГ(x)) ≡ exp(γ*x - p*logГ(x))
"""
struct GammaShapeLikelihood{ T <: Real } <: ContinuousUnivariateDistribution
    p :: T
    γ :: T # p * β
end

Distributions.@distr_support GammaShapeLikelihood 0 Inf

Distributions.support(dist::GammaShapeLikelihood) = Distributions.RealInterval(minimum(dist), maximum(dist))

Base.show(io::IO, distribution::GammaShapeLikelihood{T}) where T = print(io, "GammaShapeLikelihood{$T}(π = $(distribution.p), γ = $(distribution.γ))")

Distributions.logpdf(distribution::GammaShapeLikelihood, x::Real) = distribution.γ * x - distribution.p * loggamma(x)

prod_analytical_rule(::Type{ <: GammaShapeLikelihood }, ::Type{ <: GammaShapeLikelihood }) = ProdAnalyticalRuleAvailable()

function prod(::ProdAnalytical, left::GammaShapeLikelihood, right::GammaShapeLikelihood)
    return GammaShapeLikelihood(left.p + right.p, left.γ + right.γ)
end

# TODO over form constraints
# function approximate_prod_expectations(approximation::GaussLaguerreQuadrature, left::GammaDistributionsFamily, right::GammaShapeLikelihood)
#     b = rate(left)

#     """
#     q(x)    ∝ v(x)*v(x)
#             ∝ exp(γ*x - p*ln(Г(x))) * exp((a-1)*ln(x) - b*x)
#             = exp(-x) * exp((γ-b+1)*x + (a-1)*ln(x) - p*ln(Г(x)))
#     """
#     f = let p = right.p, a = shape(left), γ = right.γ, b = b
#         x -> exp((γ - b + 1) * x - p * loggamma(x) + (a - 1) * log(x))
#     end

#     logf = let p = right.p, a = shape(left), γ = right.γ, b = b
#         x -> (γ - b + 1) * x - p * loggamma(x) + (a - 1) * log(x)
#     end

#     # calculate log-normalization constant
#     logC = log_approximate(approximation, logf)

#     # mean function without explicitly calculating the normalization constant
#     mf = let logf = logf, logC = logC
#         x -> x * exp(logf(x) - logC)
#     end

#     # calculate mean
#     m = approximate(approximation, mf)

#     # variance function without explicitly calculating the normalization constant
#     vf = let logf = logf, logC = logC, m = m
#         x -> (x - m) ^ 2 * exp(logf(x) - logC)
#     end

#     # calculate variance
#     v = approximate(approximation, vf)

#     return m, v
# end

# function approximate_prod_expectations(approximation::ImportanceSamplingApproximation, left::GammaDistributionsFamily, right::GammaShapeLikelihood)

#     f = let p = right.p, γ = right.γ
#         x -> exp(γ * x - p * loggamma(x))
#     end

#     m, v = approximate_meancov(approximation, f, GammaShapeScale(shape(left), scale(left)))

#     return m, v
# end

# function prod(::ProdAnalytical, left::GammaShapeLikelihood, right::GammaDistributionsFamily)
#     return prod(ProdAnalytical(), right, left)
# end

# function prod(::ProdAnalytical, left::GammaDistributionsFamily, right::GammaShapeLikelihood)
#     m, v = approximate_prod_expectations(right.approximation, left, right)

#     a = m ^ 2 / v
#     b = m / v

#     return GammaShapeRate(a, b)
# end

# Expectation maximisation
# TODO: same do over form constraints

# function prod(::ProdExpectationMaximisation, left::GammaShapeLikelihood, right::GammaShapeLikelihood)
#     @assert left.approximation == right.approximation "Different approximation types for $(left) and $(right) messages"
#     return GammaShapeLikelihood(left.p + right.p, left.γ + right.γ, left.approximation)
# end

# function prod(::ProdExpectationMaximisation, left::GammaShapeLikelihood, right::GammaDistributionsFamily)
#     return prod(ProdExpectationMaximisation(), right, left)
# end

# function prod(::ProdExpectationMaximisation, left::GammaDistributionsFamily, right::GammaShapeLikelihood)

#     a, b = shape(left), rate(left)
#     γ, p = right.γ, right.p

#     f(x) = (a-1)*log(x[1]) - b*x[1] + γ*x[1] - p*loggamma(x[1]) - loggamma(a) + a*log(b)

#     x_0 = [ mean(left) ]
#     res = optimize(x -> -f(x), [ 0.0 ], [ Inf ], x_0, Fminbox(GradientDescent()))

#     â = Optim.minimizer(res)[1]

#     return PointMass(â)
# end


text:
##############################################################################
# Model and alternative constructors.
"""`Model(...)`

The `Model()` function constructs an ERM model. The typical invocation is
`Model(U, V, Loss(), Reg())`, where `U` and `V` specify raw inputs and targets,
respectively, and `Loss()` specifies some type of training loss (default: `SquareLoss()`)
and `Reg()` specifies some type of regularizer (default: `L2Reg()`).
For more details, see the description of ERM models in the usage notes.
"""

using Random

mutable struct Model
    D::Data
    loss::Loss
    regularizer::Regularizer
    solver::Solver
    S::DataSource
    X
    Y
    regweights
    istrained::Bool
    verbose::Bool
    xydataisinvalid::Bool
    disinvalid::Bool
    Uestnumcols
    Vestnumcols
    embedallwarning::Bool
end

function Model(U, V; loss=SquareLoss(), reg=L2Reg(),
               Unames = nothing, Vnames = nothing,
               embedall = true, verbose=true,
               Uestnumcols = 0, Vestnumcols=0, kwargs...)
    S = makeFrameSource(U, V, Unames, Vnames; kwargs...)
    M =  Model(NoData(),
               loss,
               reg,
               DefaultSolver(),
               S,
               nothing, nothing, #X,Y
               nothing, #regweights
               false, #istrained
               verbose, #verbose
               true, # xydataisinvalid
               true,  # disinvalid
               Uestnumcols, Vestnumcols,
               false)
    setdata(M)
    if embedall
        if M.verbose
            println("Model: applying default embedding")
        end
        defaultembedding(M; kwargs...)
        M.embedallwarning = true
    end
    return M
end
##############################################################################
function getfoldrows(n, nfolds)
    # split into groups
    groups = [convert(Int64, round(x*n/nfolds)) for x in 1:nfolds]
    unshift!(groups,0)
    # i'th fold has indices groups[i]+1:groups[i+1]
    p = randperm(n)
    foldrows = Any[]
    nonfoldrows = Any[]
    for i=1:nfolds
        push!(foldrows, sort(p[groups[i]+1:groups[i+1]]))
        push!(nonfoldrows, sort([ p[1:groups[i]] ; p[groups[i+1]+1:end]]))
    end
    return foldrows, nonfoldrows
end

function splitrows(n, trainfrac::Array{Int64,1})
    trainrows = trainfrac
    allrows = Set(1:n)
    trainset = Set(copy(trainrows))
    testrows = sort(collect(setdiff(allrows, trainset)))
    return trainrows, testrows
end

function splitrows(n, trainfrac::Number; splitmethod=0)
    if splitmethod == 0
        ntrain = convert(Int64, round(trainfrac*n))
        p = randperm(n)
        trainrows = sort(p[1:ntrain])
        testrows = sort(p[ntrain+1:n])
    else
        # pick by Bernoulli
        testrows = Int64[]
        trainrows = Int64[]
        for i=1:n
            if rand()>trainfrac
                push!(testrows,i)
            else
                push!(trainrows,i)
            end
        end
    end
    return trainrows, testrows
end

function setdata(M::Model)
    M.X, M.Y  = getXY(M.S; Uestnumcols = M.Uestnumcols, Vestnumcols = M.Vestnumcols, verbose=M.verbose)
    M.xydataisinvalid = false
    setregweights(M)
    return
end

function setregweights(M::Model)
    if M.xydataisinvalid
        return
    end
    X = M.X
    d = size(X,2)
    R = ones(d)
    for i=1:d
        if Statistics.var(X[:,i]) == 0
            if norm(X[:,i]) != 0
                R[i] = 0
            end
        end
    end
    M.regweights = R
end

# function defaultembedding(M::Model; stand=true)
#     addfeatureV(M, 1, stand=stand)
#     addfeatureU(M, etype="one")
#     d = size(getU(M),2)
#     for i=1:d
#         addfeatureU(M, i, stand=stand)
#     end
# end

function defaultembedding(M::Model; stand=true, kwargs...)
    addfeatureV(M, 1, stand=stand)
    addfeatureU(M, etype="all", stand=stand, addones=true, kwargs...)
end

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

function SplitData(X, Y, trainfrac)
    trainrows, testrows = splitrows(size(X,1), trainfrac)
    Xtrain = X[trainrows,:]
    Ytrain = Y[trainrows,:]
    Xtest = X[testrows,:]
    Ytest = Y[testrows,:]
    return SplitData(Xtrain, Ytrain, Xtest, Ytest, trainrows, testrows, trainfrac, NoResults())
end

function FoldedData(X, Y, nfolds)
    foldrows, nonfoldrows = getfoldrows(size(X,1), nfolds)
    return FoldedData(X, Y, nfolds, foldrows, nonfoldrows, NoResults())
end

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

function warndata(M)
    if length(M.S.Xmaps)==0
        println("Warning: Model has no X data. Use addfeatureU or Model(...,embedall=true)")
    end
    if length(M.S.Ymaps)==0
        println("Warning: Model has no Y data. Use addfeatureV or Model(...,embedall=true)")
    end
end
function splittraintestx(M, trainfrac)
    setdata(M)
    if M.verbose
        println("Model: splitting data")
    end
    M.D = SplitData(M.X, M.Y, usetrainfrac(M, trainfrac))
    M.disinvalid = false
end

function usetrainfrac(M::Model, trainfrac)
    if trainfrac==nothing
        if isa(M.D, SplitData)
            return M.D.trainfrac
        end
        return 0.8
    end
    return trainfrac
end

function splittraintest(M::Model; trainfrac=nothing, resplit=false, force=false)
    if resplit || force
        M.disinvalid = true
    end
    if !M.disinvalid && isa(M.D, SplitData)
        if trainfrac != nothing && trainfrac != M.D.trainfrac
            M.disinvalid = true
        end
    end
    if M.disinvalid ||  !isa(M.D, SplitData)
        splittraintestx(M, trainfrac)
    end
end

function splitfolds(M::Model, nfolds; resplit=false, force=false)
    if resplit || force
        M.disinvalid = true
    end
    if !M.disinvalid && isa(M.D, FoldedData)
        if M.D.nfolds != nfolds
            M.disinvalid = true
        end
    end

    if M.disinvalid || !isa(M.D, FoldedData)
        setdata(M)
        M.D = FoldedData(M.X, M.Y, nfolds)
        M.disinvalid = false
    end
end

##############################################################################
# fit

function trainx(M::Model, lambda, Xtrain, Xtest, Ytrain, Ytest; theta_guess = nothing)
    assignsolver(M)
    if M.verbose
        println("Model: calling solver: ", M.solver)
        for i=1:length(M.regweights)
            if M.regweights[i] == 0
                println("Model: Not regularizing constant feature X[:,$(i)]")
            end
        end
    end
    theta = solve(M.solver, M.loss, M.regularizer, M.regweights, Xtrain, Ytrain, lambda;
                  theta_guess = theta_guess)
    trainloss = loss(M.loss, predict(M, Xtrain, theta), Ytrain)
    testloss = loss(M.loss,  predict(M, Xtest,  theta), Ytest)
    return PointResults(theta, lambda, trainloss, testloss)
end

function trainfoldsx(M::Model, lambda, nfolds)
    results = Array{PointResults}(nfolds)
    for i=1:nfolds
        results[i] =  trainx(M, lambda, Xtrain(M,i), Xtest(M,i), Ytrain(M,i), Ytest(M,i))
    end
    return FoldResults(results)
end

function trainpathx(M::Model, lambda::Array; quiet=true, kwargs...)
    m = length(lambda)
    results = Array{PointResults}(m)
    for i=1:m
        if i>1
            tg = results[i-1].theta
        else
            tg = nothing
        end
        if M.verbose
            println("lambda = ", lambda[i])
        end
        results[i] = trainx(M, lambda[i], Xtrain(M), Xtest(M), Ytrain(M), Ytest(M);
                            theta_guess = tg)
    end
    imin =  findmin([x.testloss for x in results])[2]
    return RegPathResults(results, imin)
end

"""`trainfolds(M [,lambda=1e-10, nfolds=5])`

The `trainfolds` function carries out n-fold cross validation on
a model `M`. Specify regularization weight through optional argument
`lambda`, and the number of folds through `nfolds`. Default:
`lambda=1e-10`, and `nfolds=5`.
"""
function trainfolds(M::Model; lambda=1e-10, nfolds=5,
                    resplit=false, features=nothing, kwargs...)
    warndata(M)
    splitfolds(M, nfolds, resplit)
    M.D.results = trainfoldsx(M, lambda, nfolds)
    M.istrained = true
    if M.verbose
        status(M)
    end
end

"""`trainpath(M [,lambda=logspace(-5, 5, 100), trainfrac=0.8])`

The `trainpath` function trains a model `M` over a set of regularization weights.
Specify these weights by invoking the optional argument `lambda`, and set a train-test
ratio by using the optional argument `trainfrac`.

Defaults:
 -  `lambda=logspace(-5,5,100)` so training occurs over `lambda` between 1e-5 and
1e5.
 - `trainfrac=0.8`, so training occurs with a 80-20 train-test split.

Example: `trainpath(M, lambda=logspace(-1, 1, 100))` trains over `lambda` between
0.1 and 10.

Example `trainpath(M, trainfrac=0.75)` trains w/ 75-25 train-test split.
"""
function trainpath(M::Model; lambda=logspace(-5,5,100), trainfrac=0.8,
                   resplit=false, features=nothing, kwargs...)
    warndata(M)
    splittraintest(M; trainfrac=trainfrac, resplit=resplit)
    M.D.results = trainpathx(M, lambda)
    M.istrained = true
    if M.verbose
        status(M)
    end
end

"""`train(M [, lambda=1e-10, trainfrac=nothing])`
This function trains a model `M`. The usual invocation is
`train(M)`. Users may choose to specify a different choice of
regularization weight `lambda`. For example to specify
a weight of `lambda = 0.01`, one invokes
`train(M, lambda=0.001)`, and to specify a different train split,
one invokes `train(M, trainfrac=0.75)`, which means that
75 percent of the data will be used for training and only 25 percent will
be used for test. The default parameters are
`lambda = 1e-10` and `trainfrac=nothing`, which will result in a
80-20 train-test split."""
function train(M::Model; lambda=1e-10, trainfrac=nothing,
               resplit=false, kwargs...)
    warndata(M)
    splittraintest(M; trainfrac=trainfrac, resplit=resplit)
    M.D.results = trainx(M, lambda, Xtrain(M), Xtest(M), Ytrain(M), Ytest(M); kwargs...)
    M.istrained = true
    if M.verbose
        status(M)
    end
end

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

function assignsolver(M::Model, force=false)
    if force || isa(M.solver, DefaultSolver)
        M.solver = getsolver(M.loss, M.regularizer)
    end
end

function setloss(M::Model, l)
    M.loss = l
    assignsolver(M, true)
end
function setreg(M::Model, r)
    M.regularizer = r
    assignsolver(M, true)
end
function setsolver(M::Model, s)
    if s == "default"
        assignsolver(M, true)
        return
    end
    M.solver = s
end

##############################################################################
# querying

getU(M::Model) = getU(M.S)
getV(M::Model) = getV(M.S)
#getXY(M::Model) = getXY(M.S)

Xtest(M::Model) = Xtest(M.D)
Xtrain(M::Model) =  Xtrain(M.D)
Xtrain(D::SplitData) = D.Xtrain
Xtest(D::SplitData) = D.Xtest
Utrain(M::Model) = getU(M.S)[M.D.trainrows,:]
Utest(M::Model)  = getU(M.S)[M.D.testrows,:]
Vtrain(M::Model) = getV(M.S)[M.D.trainrows,:]
Vtest(M::Model)  = getV(M.S)[M.D.testrows,:]

Xtest(M::Model, fold) = Xtest(M.D, fold)
Xtrain(M::Model, fold) = Xtrain(M.D, fold)
Xtrain(D::FoldedData, fold) = D.X[D.nonfoldrows[fold],:]
Ytrain(D::FoldedData, fold) = D.Y[D.nonfoldrows[fold],:]
Xtest(D::FoldedData, fold) = D.X[D.foldrows[fold],:]
Ytest(D::FoldedData, fold) = D.Y[D.foldrows[fold],:]
Ytest(M::Model, fold) = Ytest(M.D, fold)
Ytrain(M::Model, fold) = Ytrain(M.D, fold)

Ytest(M::Model) = Ytest(M.D)
Ytrain(M::Model) = Ytrain(M.D)
Ytrain(D::SplitData) = D.Ytrain
Ytest(D::SplitData) = D.Ytest

import Base.split
split(M::Model; kwargs...) = splittraintest(M; force=true, kwargs...)

function warnembeddings(M)
    if !M.verbose
        return
    end
    if !M.embedallwarning
        return
    end
    println("Warning: You are adding features to a model which was created with embedall=true")
end

function addfeatureU(M::Model; kwargs...)
    warnembeddings(M)
    M.xydataisinvalid = true
    M.disinvalid = true
    addfeatureU(M.S, nothing; kwargs...)
end

function addfeatureV(M::Model; kwargs...)
    warnembeddings(M)
    M.xydataisinvalid = true
    M.disinvalid = true
    addfeatureV(M.S, nothing; kwargs...)
end

function addfeatureU(M::Model, col; kwargs...)
    warnembeddings(M)
    M.xydataisinvalid = true
    M.disinvalid = true
    addfeatureU(M.S, col; kwargs...)
end

function addfeatureV(M::Model, col; kwargs...)
    warnembeddings(M)
    M.xydataisinvalid = true
    M.disinvalid = true
    addfeatureV(M.S, col; kwargs...)
end

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

function status(io::IO, R::RegPathResults)
    println(io, "----------------------------------------")
    println(io, "Optimal results along regpath")
    println(io, "  optimal lambda: ",  lambdaopt(R))
    println(io, "  optimal test loss: ", testloss(R))
end

function status(io::IO, R::PointResults)
    println(io, "----------------------------------------")
    println(io, "Results for single train/test")
    println(io, "  training  loss: ", trainloss(R))
    println(io, "  test loss: ", testloss(R))
end

"""
Prints and returns the status of the model.
"""
function status(io::IO, M::Model)
    status(io, M.D.results)
    println(io, "  training samples: ", length(Ytrain(M)))
    println(io, "  test samples: ", length(Ytest(M)))
    println(io, "  columns in X: ", size(M.X,2))
    println(io, "----------------------------------------")
end

"""`status(M)` prints the status of the model after the most recent action performed on it."""
status(M::Model; kwargs...)  = status(stdout, M; kwargs...)


text:
module NavNodeTests

using Test

using Compat
import Documenter: Documents, Builder
import Documenter.Documents: NavNode

mutable struct FakeDocumentInternal
    pages   :: Dict{String, Void}
    navlist :: Vector{NavNode}
    FakeDocumentInternal() = new(Dict(), [])
end
mutable struct FakeDocument
    internal :: FakeDocumentInternal
    FakeDocument() = new(FakeDocumentInternal())
end

@testset "NavNode" begin
    @test fieldtype(FakeDocumentInternal, :navlist) == fieldtype(Documents.Internal, :navlist)

    pages = [
        "page1.md",
        "Page2" => "page2.md",
        "Section" => [
            "page3.md",
            "Page4" => "page4.md",
            "Subsection" => [
                "page5.md",
            ],
        ],
        "page6.md",
    ]
    doc = FakeDocument()
    doc.internal.pages = Dict(map(i -> "page$i.md" => nothing, 1:8))
    navtree = Builder.walk_navpages(pages, nothing, doc)
    navlist = doc.internal.navlist

    @test length(navlist) == 6
    for (i,navnode) in enumerate(navlist)
        @test navnode.page == "page$i.md"
    end

    @test isa(navtree, Vector{NavNode})
    @test length(navtree) == 4
    @test navtree[1] === navlist[1]
    @test navtree[2] === navlist[2]
    @test navtree[4] === navlist[6]

    section = navtree[3]
    @test section.title_override == "Section"
    @test section.page === nothing
    @test length(section.children) == 3

    navpath = Documents.navpath(navlist[5])
    @test length(navpath) == 3
    @test navpath[1] === section
    @test navpath[3] === navlist[5]
end

end