text:
using Test
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCell
include("useful_routines.jl")

function displacement(alpha, x)
    u1 = alpha * x[2] * sin(pi * x[1])
    u2 = alpha * (x[1]^3 + cos(pi * x[2]))
    return [u1, u2]
end

function stress_field(lambda, mu, alpha, x)
    s11 =
        (lambda + 2mu) * alpha * pi * x[2] * cos(pi * x[1]) -
        lambda * alpha * pi * sin(pi * x[2])
    s22 =
        -(lambda + 2mu) * alpha * pi * sin(pi * x[2]) +
        lambda * alpha * pi * x[2] * cos(pi * x[1])
    s12 = alpha * mu * (3x[1]^2 + sin(pi * x[1]))
    return [s11, s22, s12]
end

function body_force(lambda, mu, alpha, x)
    b1 = alpha * (lambda + 2mu) * pi^2 * x[2] * sin(pi * x[1])
    b2 =
        -alpha * (6mu * x[1] + (lambda + mu) * pi * cos(pi * x[1])) +
        alpha * (lambda + 2mu) * pi^2 * cos(pi * x[2])
    return [b1, b2]
end

function onboundary(x, L, W)
    return x[2] ≈ 0.0 || x[1] ≈ L || x[2] ≈ W || x[1] ≈ 0.0
end

function solve_elasticity(distancefunc, nelmts, polyorder, numqp, penaltyfactor)
    L = 1.0
    W = 1.0
    lambda, mu = 1.0, 2.0
    dx = 1.0 / nelmts
    penalty = penaltyfactor / dx * (lambda + mu)
    alpha = 0.1
    stiffness = CutCell.HookeStiffness(lambda, mu, lambda, mu)

    basis = TensorProductBasis(2, polyorder)
    mesh = CutCell.Mesh([0.0, 0.0], [L, W], [nelmts, nelmts], basis)
    levelset = InterpolatingPolynomial(1, basis)
    levelsetcoeffs = CutCell.levelset_coefficients(distancefunc, mesh)

    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    facequads = CutCell.FaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition =
        CutCell.coherent_interface_condition(basis, interfacequads, stiffness, cutmesh, penalty)

    displacementbc = CutCell.DisplacementCondition(
        x -> displacement(alpha, x),
        basis,
        facequads,
        stiffness,
        cutmesh,
        x -> onboundary(x, L, W),
        penalty,
    )

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
    CutCell.assemble_interface_condition!(sysmatrix, interfacecondition, cutmesh)
    CutCell.assemble_body_force_linear_form!(
        sysrhs,
        x -> body_force(lambda, mu, alpha, x),
        basis,
        cellquads,
        cutmesh,
    )
    CutCell.assemble_penalty_displacement_bc!(sysmatrix, sysrhs, displacementbc, cutmesh)

    matrix = CutCell.make_sparse(sysmatrix, cutmesh)
    rhs = CutCell.rhs(sysrhs, cutmesh)

    sol = matrix \ rhs

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_stress_mass_matrix!(sysmatrix, basis, cellquads, cutmesh)
    CutCell.assemble_stress_linear_form!(sysrhs, basis, cellquads, stiffness, sol, cutmesh)

    matrix = CutCell.make_sparse_stress_operator(sysmatrix, cutmesh)
    rhs = CutCell.stress_rhs(sysrhs, cutmesh)

    stressvec = matrix \ rhs
    stress = reshape(stressvec, 3, :)
    err = mesh_L2_error(
        stress,
        x -> stress_field(lambda, mu, alpha, x),
        basis,
        cellquads,
        cutmesh,
    )

    return err
end

function mean(v)
    return sum(v) / length(v)
end

function convergence_rate(v, dx)
    diff(log.(v)) ./ diff(log.(dx))
end

x0 = [0.5, 0.0]
normal = [1.0, 0.0]
polyorder = 1
numqp = required_quadrature_order(polyorder) + 2
penaltyfactor = 1e2
powers = [3,4,5]
nelmts = [2^p + 1 for p in powers]
dx = 1.0 ./ nelmts

err = [
    solve_elasticity(
        x -> plane_distance_function(x, normal, x0),
        ne,
        polyorder,
        numqp,
        penaltyfactor,
    ) for ne in nelmts
]
serr = [[er[i] for er in err] for i = 1:3]

rates = [convergence_rate(v, dx) for v in serr]
@test all([allapprox(rates[i],repeat([1.5],length(rates[i])),0.05) for i = 1:3])

# powers = [3,4,5]
# nelmts = [2^p + 1 for p in powers]
# dx = 1.0 ./ nelmts
# polyorder = 2
# numqp = required_quadrature_order(polyorder) + 2
# err = [
#     solve_elasticity(
#         x -> plane_distance_function(x, normal, x0),
#         ne,
#         polyorder,
#         numqp,
#         penaltyfactor,
#     ) for ne in nelmts
# ]
# serr = [[er[i] for er in err] for i = 1:3]
# rates = [convergence_rate(v, dx) for v in serr]
# @test all([allapprox(rates[i],repeat([2.0],length(rates[i])),0.05) for i = 1:3])

powers = [3,4,5]
nelmts = [2^p + 1 for p in powers]
dx = 1.0 ./ nelmts
xc = [1.0, 0.5]
radius = 0.45
polyorder = 2
numqp = required_quadrature_order(polyorder) + 2
err = [
    solve_elasticity(
        x -> circle_distance_function(x, xc, radius),
        ne,
        polyorder,
        numqp,
        penaltyfactor,
    ) for ne in nelmts
]
serr = [[er[i] for er in err] for i = 1:3]
rates = [convergence_rate(v, dx) for v in serr]
@test all([all(rates[i] .> 1.95) for i = 1:3])

# powers = [2,3,4]
# nelmts = [2^p + 1 for p in powers]
# dx = 1.0 ./ nelmts
# polyorder = 3
# numqp = required_quadrature_order(polyorder) + 2
# err = [
#     solve_elasticity(
#         x -> plane_distance_function(x, normal, x0),
#         ne,
#         polyorder,
#         numqp,
#         penaltyfactor,
#     ) for ne in nelmts
# ]
# serr = [[er[i] for er in err] for i = 1:3]
# rates = [convergence_rate(v, dx) for v in serr]
# @test all([all(rates[i] .> 3.0) for i = 1:3])

powers = [3,4,5]
nelmts = [2^p + 1 for p in powers]
dx = 1.0 ./ nelmts
xc = [1.0, 0.5]
radius = 0.45
polyorder = 3
numqp = required_quadrature_order(polyorder) + 2
err = [
    solve_elasticity(
        x -> circle_distance_function(x, xc, radius),
        ne,
        polyorder,
        numqp,
        penaltyfactor,
    ) for ne in nelmts
]
serr = [[er[i] for er in err] for i = 1:3]
rates = [convergence_rate(v, dx) for v in serr]
@test all([all(rates[i] .> 2.95) for i = 1:3])


text:
module ModuleNoise
using MultipleViewGeometry.ModuleTypes, MultipleViewGeometry.ModuleMathAliases, MultipleViewGeometry.ModuleOperators
using StaticArrays
export perturb
include("perturb.jl")
end

ext: jl


text:
module SolverLogging
using Printf
using Formatting
using Crayons
using Logging

include("utils.jl")
include("conditional_crayon.jl")
include("logger.jl")
# include("default_logger.jl")
include("setentry.jl")

# Set entries in default logger
"""
    @log [logger] args...

Logs the data with `logger`, which is the default logger if not provided.
Can be any of the following forms:

    @log "name" 1.0    # literal value
    @log "name" 2a     # an expression
    @log "name" name   # a variable
    @log name          # shortcut for previous

If the specified entry is active at the current verbosity level,
"""
macro log(args...)
    return log_expr(args...)
end

# Use the symbol as the string name
# e.g. `@log α`` logs the value of α in the "α" field
function log_expr(logger, ex::Symbol)
    name = string(ex)
    _log_expr(logger, name, ex)
end

# Pass-through to actual function
log_expr(logger, name::String, ex) = _log_expr(logger, name, ex)
log_expr(logger, name::String, ex, op::QuoteNode) = _log_expr(logger, name, ex, op)

function _log_expr(logger, name::String, expr, op::QuoteNode=QuoteNode(:replace))
    quote
        let lg = $(esc(logger))
            # @debug "Calling @log on $name"
            if isenabled(lg)
                # espec = lg.fmt[$name]
                # if lg.opts.curlevel <= espec.level
                _log!(lg, $name, $(esc(expr)), $op)
                # end
            end
        end
    end
end

export
    setentry,
    Logger,
    @log,
    printheader,
    printrow,
    printlog,
    setlevel!,
    ConditionalCrayon

end # module


text:
"""
	nodes(t)
	leaves(t)
	internals(t)

Iterator over all nodes / leaves / internal nodes of a tree.

# Note
`length` cannot be called on `internals(t)` as the latter is based on `Iterators.filter`.
  A way to get the number of internal nodes of a tree is for example by calling
  `length(nodes(t)) - length(leaves(t))`.
"""
nodes, leaves, internals
nodes(t::Tree) = values(t.lnodes)
leaves(t::Tree) = values(t.lleaves)
internals(t::Tree) = Iterators.filter(x->!x.isleaf, values(t.lnodes))

nodes(f::Function, t::Tree) = filter(f, values(t.lnodes))

#=
Post order traversal iterators
=#
abstract type POTIterator end

Base.IteratorSize(::Type{POTIterator}) = Iterators.HasLength()
Base.IteratorEltype(::Type{POTIterator}) = HasEltype()

Base.iterate(itr::POTIterator) = firststate(itr, itr.root)

"""
	struct POT{T<:TreeNodeData} <: POTIterator
		root::TreeNode{T}
	end
"""
struct POT{T<:TreeNodeData} <: POTIterator
	root::TreeNode{T}
end
Base.eltype(::Type{POT{T}}) where T = TreeNode{T}
POT(t::Tree) = POT(t.root)
function Base.length(iter::POT)
	return count(x->true, iter.root)
end

"""
	struct POTleaves{T<:TreeNodeData} <: POTIterator
		root::TreeNode{T}
	end
"""
struct POTleaves{T<:TreeNodeData} <: POTIterator
	root::TreeNode{T}
end
Base.eltype(::Type{POTleaves{T}}) where T = TreeNode{T}
POTleaves(t::Tree) = POTleaves(t.root)
function Base.length(iter::POTleaves)
	return count(isleaf, iter.root)
end

struct POTState{T<:TreeNodeData}
	n::TreeNode{T}
	i::Int # Position of n in list of siblings -- `n.anc.child[i]==n`
	direction::Symbol
end

"""
- `state.n.isleaf`: go to sibling and down or ancestor and up (stop if root)
- Otherwise: go to deepest child and up.
"""
function Base.iterate(itr::POTIterator, state::POTState)
	if state.direction == :down
		return go_down(itr, state)
	elseif state.direction == :up
		return go_up(itr, state)
	elseif state.direction == :stop
		return nothing
	end
end

function go_down(itr::POTIterator, state::POTState{T}) where T
	if state.n.isleaf # Go back to ancestor or sibling anyway
		if state.n.isroot || state.n == itr.root
			return (state.n, POTState{T}(n, 0, :stop))
		elseif state.i < length(state.n.anc.child) # Go to sibling
			return (state.n, POTState{T}(state.n.anc.child[state.i+1], state.i+1, :down))
		else # Go back to ancestor
			return (state.n, POTState{T}(state.n.anc, get_sibling_number(state.n.anc), :up))
		end
	end
	return firststate(itr, state.n) # Go to deepest child of `n`
end
function go_up(itr::POT{T}, state::POTState{T}) where T
	if state.n.isroot || state.n == itr.root
		return (state.n, POTState{T}(state.n, 0, :stop))
	elseif state.i < length(state.n.anc.child) # Go to sibling
		return (state.n, POTState{T}(state.n.anc.child[state.i+1], state.i+1, :down))
	else # Go back to ancestor
		return (state.n, POTState{T}(state.n.anc, get_sibling_number(state.n.anc), :up))
	end
end
function go_up(itr::POTleaves{T}, state::POTState{T}) where T
	if state.n.isleaf
		if state.n.isroot || state.n == itr.root
			return (state.n, POTState{T}(state.n, 0, :stop))
		elseif state.i < length(state.n.anc.child) # Go to sibling
			return (state.n, POTState{T}(state.n.anc.child[state.i+1], state.i+1, :down))
		else # Go back to ancestor
			return (state.n, POTState{T}(state.n.anc, get_sibling_number(state.n.anc), :up))
		end
	else
		if state.n.isroot || state.n == itr.root
			return nothing
		elseif state.i < length(state.n.anc.child) # Go to sibling
			iterate(itr, POTState(state.n.anc.child[state.i+1], state.i+1, :down))
		else # Go back to ancestor
			iterate(itr, POTState(state.n.anc, get_sibling_number(state.n.anc), :up))
		end
	end
end

"""
Go to deepest child of `a`.
"""
function firststate(itr::POTIterator, a::TreeNode{T}) where T
	if a.isleaf
		return iterate(itr, POTState{T}(a, 1, :up))
	end
	firststate(itr, a.child[1])
end

function get_sibling_number(n::TreeNode)
	if n.isroot
		return 0
	end
	for (i,c) in enumerate(n.anc.child)
		if n == c
			return i
		end
	end
	@error "Could not find $(n.label) in children of $(n.anc.label)."
end

max_stars_repo_name: PierreBarrat/TreeTools

max_issues_repo_name: PierreBarrat/TreeTools


text:
type MultiCubicSpline{T<:Real}
    x::Vector{T}
    fx::Matrix{T}
    fdp::Matrix{T}
    npts::Ref{Cint}
    dirty::Bool
end

function CubicSpline(x::Vector{Float64}, fx::Matrix{Float64})
    N, M = size(fx)
    npts = Ref{Cint}(N)
    fdp = zeros(fx)

    for i in 1:M
        _dodp!(x, view(fx, :, i), view(fdp, :,i), npts)
    end

    MultiCubicSpline(x, fx, fdp, npts, false)
end

function update!(cs::MultiCubicSpline, fx::Matrix{Float64})
    cs.fdp .= 0.0
    cs.fx .= fx .+ 0.0
    M = size(fx)[2]

    for i in 1:M
        _dodp!(cs.x, cs.fx[:,i], cs.fdp[:,i], cs.npts)
    end

    cs.dirty = false
end

function interp_level(cs::MultiCubicSpline{Float64},
                      x::Float64, i::Int)
    y = zeros(1)
    yp = zeros(1)
    ydp = zeros(1)

    _interp!(x, cs.x, cs.fx[:,i], cs.fdp[:,i], cs.npts,
             1, 0, 0, y, yp, ydp)
    y[1]
end

function(itp::MultiCubicSpline)(x::Real, i::Int, update::Bool=true)
    if update && itp.dirty
        update!(itp, itp.fx)
    end
    interp_level(itp, x, i)
end


text:
function optimize_acquisition(dim::Int64,
                              #max_feval::Int64,
                              max_time::Float64,
                              x_hint::AbstractVector,
                              y_opt::Real,
                              X::Matrix,
                              K, a, k;
                              verbose::Bool=true,
                              prng=Random.GLOBAL_RNG)
    α(x) = -expected_improvement(x, y_opt, X, K, a, k)
    dims = length(x_hint)

    res = CMAEvolutionStrategy.minimize(
        α,
        x_hint,
        0.5*√(dim);
        lower=zeros(dim),
        upper=ones(dim),
        seed=prng.seed[1],
        xtol=1e-4,
        ftol=1e-5,
        #maxfevals=max_feval,
        maxtime=max_time,
        multi_threading=true,
        verbosity=verbose ? 3 : 0
    )

    x_sol = CMAEvolutionStrategy.xbest(res)
    y_sol = α(x_sol)
    x_sol = clamp.(x_sol, 0, 1)
    return x_sol, y_sol
end

function optimize_mean(dim::Int64,
                       max_feval::Real,
                       X::Matrix,
                       K, a, k)
    data_μs = [ gp_predict(X[:,i], X, K, a, k)[1] for i = 1:size(X,2) ]
    opt_idx    = argmax(data_μs)
    x_opt      = X[:,opt_idx]
    y_opt      = data_μs[opt_idx]
    x_opt, y_opt, opt_idx
end


text:
using Flux, Flux.Tracker
using Distributions

struct VP
  b0
  b1
  log_sig

  VP() = new(param(randn(2)), param(randn(2)), param(randn(2)))
end

function lpdf_logx(dist, logx)
  return logpdf(dist, exp(logx)) + logx
end

function loglike(y, x, b0, b1, sig, N)
  return sum(logpdf.(Normal.(b0 .+ b1 .* x, sig), y)) * (N / length(y))
end

function log_p(b0, b1, log_sig)
  return logpdf(Normal(0, 1), b0) + logpdf(Normal(0, 1), b1) + lpdf_logx(Gamma(1, 1), log_sig)
end

lpdf_vp(p, x) = logpdf(Normal(p[1], exp(p[2])), x)
function log_q(b0, b1, log_sig, vp)
  return lpdf_vp(vp.b0, b0) + lpdf_vp(vp.b1, b1) + lpdf_vp(vp.log_sig, log_sig)
end

rsample(p) = randn() * exp(p[2]) + p[1]
function elbo(y, x, vp, N)
  b0 = rsample(vp.b0)
  b1 = rsample(vp.b1)
  log_sig = rsample(vp.log_sig)
  sig = exp(log_sig)

  return loglike(y, x, b0, b1, sig, N) + log_p(b0, b1, log_sig) - log_q(b0, b1, log_sig, vp)
end


text:
using Ramnet.Partitioners:
partition,
RandomPartitioner,
random_partitioning,
LinearPartitioner

@testset "Random partitioning (and generic methods)" begin
    input_v        = Bool[1, 0, 0, 1, 0]
    input_m        = Bool[0 0 1 1 1; 1 0 0 0 0]

    shorter_v      = Bool[1, 0, 1]
    longer_v       = Bool[1, 1, 1, 1, 1, 0, 0, 1]

    shorter_m      = Bool[1 0 1; 1 0 1]
    longer_m       = Bool[1 1 1 1 1 0 0 1; 1 1 1 1 1 0 0 1]

    partition_size = 2

    partitioner = RandomPartitioner(length(input_v), partition_size)

    expected_len = ceil(Int, length(input_v) / partition_size)

    @test eltype(partition(partitioner, input_v)) == typeof(input_v)
    @test eltype(partition(partitioner, input_m)) == typeof(input_m)

    @test length(partitioner) == expected_len
    @test length(partition(partitioner, input_v)) == expected_len
    @test length(partition(partitioner, input_m)) == expected_len

    @test length(collect(partition(partitioner, input_v))) == expected_len
    @test length(collect(partition(partitioner, input_m))) == expected_len
    @test length(collect(random_partitioning(input_v, partition_size))) == expected_len
    @test length(collect(random_partitioning(input_m, partition_size))) == expected_len

    @test_throws DimensionMismatch partition(partitioner, shorter_v)
    @test_throws DimensionMismatch partition(partitioner, longer_v)
    @test_throws DimensionMismatch partition(partitioner, shorter_m)
    @test_throws DimensionMismatch partition(partitioner, longer_m)
end

@testset "Linear partitioning" begin
    input_v = Int[1, 2, 3, 4, 5]
    input_m = Int[1 2 3 4 5; 6 7 8 9 10]

    partition_size = 3
    partitioner = LinearPartitioner(length(input_v), partition_size)

    @test collect(partition(partitioner, input_v)) == [
        [1, 2, 3],
        [4, 5]
    ]

    @test collect(partition(partitioner, input_m)) == [
        [1 2 3; 6 7 8],
        [4 5; 9 10]
    ]
end

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet


text:
using BinaryOptimization

# hamiltonian cycle
function hamiltonian_cycle_qubo(n::Int, vars::Variables, A1::Real=1, A2::Real=A1)
    @assert vars.indexcapacities == [n, n]
    hobo = HOBO(vars)

    for t=1:n
        hobo_tmp = HOBO(vars)
        for j=1:n
            hobo_tmp += vars[t, j]
        end
        hobo += A1 * (1. - hobo_tmp)^2
    end

    for j=1:n
        hobo_tmp = HOBO(vars)
        for t=1:n
            hobo_tmp += vars[t, j]
        end
        hobo += A2 * (1. - hobo_tmp)^2
    end
    hobo
end

function hamiltonian_cycle_qubo(n::Int, name::String="Hamiltonian Cycle Comlete Graph")
    hamiltonian_cycle_qubo(n, Variables(name, [n, n], ["timestep", "vertex"]))
end

# tsp
function tsp_qubo(w::Matrix{T}, vars::Variables, B::Real=1., A1::Real=2*B*maximum(w), A2::Real=2*B*maximum(w)) where T<:Real
    n = size(w, 1)
    @assert size(w, 1) == size(w, 2)
    @assert vars.indexcapacities == [n, n]
    @assert A1 >= 0 && A2>=0 && B >= 0

    hobo_hamiltonian = hamiltonian_cycle_qubo(n, vars, A1, A2)
    hobo_tsp = HOBO(vars)
    #remove according to graph
    for u = 1:n, v = 1:n, t = 1:n
        if u != v
            hobo_tsp += w[u, v] * vars[t, u] * vars[t % n + 1, v]
        end
    end
    hobo_hamiltonian + B*hobo_tsp
end

function tsp_qubo(w::Matrix; name::String="TSP", B::Real=1., A::Real=2*B*maximum(w))
    tsp_qubo(w, Variables(name, collect(size(w)), ["timestep", "city"]), B, A, A)
end

max_stars_repo_name: iitis/ec-qaoa-code

max_issues_repo_name: iitis/ec-qaoa-code

max_forks_repo_name: iitis/ec-qaoa-code


text:
## Extending docs format support --------------------------------------------------------

parsedocs(ds::Docile.Interface.Docs{:md}) = Markdown.parse(data(ds))

## Common -------------------------------------------------------------------------------
const MDHTAGS = ["#", "##", "###", "####", "#####", "######"]
const MDSTYLETAGS = ["", "*", "**"]
const MD_SUBHEADER_OPTIONS = [:skip, :simple, :category]

file"docs/config.md"
type Config
    category_order         :: Vector{Symbol}
    metadata_order         :: Vector{Symbol}
    include_internal       :: Bool
    mathjax                :: Bool
    mdstyle_header         :: ASCIIString
    mdstyle_objname        :: ASCIIString
    mdstyle_meta           :: ASCIIString
    mdstyle_subheader      :: ASCIIString
    mdstyle_index_mod      :: ASCIIString
    md_subheader           :: Symbol
    md_index_modprefix     :: ByteString
    md_permalink           :: Bool

    const defaults = Dict{Symbol, Any}([
        (:category_order         , [:module, :function, :method, :type,
                                    :typealias, :macro, :global, :comment]),
        (:metadata_order         , [:source]),
        (:include_internal       , true),
        (:mathjax                , false),
        (:mdstyle_header         , "#"),
        (:mdstyle_objname        , "####"),
        (:mdstyle_meta           , "*"),
        (:mdstyle_subheader      , "##"),
        (:mdstyle_index_mod      , "##"),
        (:md_subheader           , :simple),
        (:md_index_modprefix     , "MODULE: "),
        (:md_permalink           , true)
        ])

    """
    Returns a default Config. If any args... are given these will overwrite the defaults.

    ```
    using Lexicon
    config = Config(md_permalink = false, mathjax = true)

    ```
    """
    function Config(; args...)
        return update_config!(new(), merge(defaults, Dict(args)))
    end
end

function update_config!(config::Config, args::Dict)
    for (k, v) in args
        try
            k in fieldnames(Config) ? setfield!(config, k, v) : warn("Invalid setting: '$(k) = $(v)'.")
        catch err # e.g. TypeError
            warn("Invalid setting: '$(k) = $(v)'. Error: $err")
        end
    end

    for k in [:mdstyle_header, :mdstyle_objname, :mdstyle_meta, :mdstyle_subheader, :mdstyle_index_mod]
        getfield(config, k) in vcat(MDHTAGS, MDSTYLETAGS) ||
                error("""Invalid mdstyle : config-item `$k -> $(getfield(config, k))`.
                      Valid values: [$(join(vcat(MDHTAGS, MDSTYLETAGS), ", "))].""")

        config.md_subheader in MD_SUBHEADER_OPTIONS ||
                    error("""Invalid md_subheader : config-item `$k -> $(getfield(config, k))`.
                          Valid values: $MD_SUBHEADER_OPTIONS.""")
    end
    return config
end

type Entries
    sourcepaths :: Vector{UTF8String}
    modulenames :: Vector{UTF8String}
    include_internal::Bool
    exported::Dict{Symbol, Vector{@compat(Tuple{Module, Any, AbstractEntry, AbstractString})}}
    internal::Dict{Symbol, Vector{@compat(Tuple{Module, Any, AbstractEntry, AbstractString})}}

    Entries(config::Config) = new(UTF8String[], UTF8String[],
                                    config.include_internal,
                                    Dict([(c, []) for c in config.category_order]),
                                    Dict([(c, []) for c in config.category_order]))
end

function has_items(entries::Dict)
    return sum([length(x) for x in values(entries)]) > 0
end

function push!(ents::Entries, modulename::Module, obj, ent::AbstractEntry,
                                anchorname::AbstractString,  cat::Symbol)
    if isexported(modulename, obj)
        push!(ents.exported[cat], (modulename, obj, ent, anchorname))
    elseif ents.include_internal
        push!(ents.internal[cat], (modulename, obj, ent, anchorname))
    end
end

type Index
    entries::Vector{Entries}
end
Index() = Index(Vector{Entries}[])

function update!(index::Index, ents::Entries)
    push!(index.entries, ents)
end

function mainsetup(io::IO, mime::MIME, doc::Metadata, ents::Entries,
                              filepath::AbstractString, config::Config)
    # Root may be a file or directory. Get the dir.
    rootdir = isfile(root(doc)) ? dirname(root(doc)) : root(doc)
    for file in manual(doc)
        writemime(io, mime, readall(joinpath(rootdir, file)))
    end
    # TODO: this index should be renamed to not confuse with the `Index` type
    idx = Dict{Symbol, Any}()
    for (obj, entry) in entries(doc)
        addentry!(idx, obj, entry)
    end

    if !isempty(idx)
        ents = prepare_entries(idx, ents, doc, config)
        if (has_items(ents.exported) || has_items(ents.internal))
            push!(ents.sourcepaths, abspath(filepath))
            push!(ents.modulenames, string(modulename(doc)))
            process_entries(io, mime, "Exported", ents.exported, config)
            process_entries(io, mime, "Internal", ents.internal, config)
        end
    end
    return ents
end

function prepare_entries(idx::Dict{Symbol, Any}, ents::Entries, doc::Metadata, config::Config)
    pageanchors = Dict{Symbol, Dict{AbstractString, Int}}([])
    for k in config.category_order
        haskey(idx, k) || continue
        k in keys(pageanchors) || (pageanchors[k] = Dict())
        basenames = pageanchors[k]
        for (s, obj) in idx[k]
            ent = entries(doc)[obj]
            if k == :comment
                basename = "comment"
            else
                # adjust basename to be a valid html id
                basename = generate_html_id(string(k == :macro ? macroname(ent) : name(obj)))
            end
            basename in keys(basenames)                  ?
                    anchornum = basenames[basename] += 1 :
                    anchornum = basenames[basename] = 1
            string(k, "_", generate_html_id(s))
            push!(ents, modulename(doc), obj, entries(doc)[obj],
                    "$(string(k))__$(basename).$(anchornum)", k)
        end
    end
    return ents
end

file"docs/save.md"
function save(file::AbstractString, modulename::Module, config::Config; args...)
    config = update_config!(deepcopy(config), Dict(args))
    mime = MIME("text/$(strip(last(splitext(file)), '.'))")
    index_entries = save(file, mime, documentation(modulename), config)
    return index_entries
end
save(file::AbstractString, modulename::Module; args...) = save(file, modulename, Config(); args...)

"""
Saves an *API-Index* to `file`.
"""
function save(file::AbstractString, index::Index, config::Config; args...)
    config = update_config!(deepcopy(config), Dict(args))
    save(file, MIME("text/$(strip(last(splitext(file)), '.'))"), index, config)
end
save(file::AbstractString, index::Index; args...) = save(file, index, Config(); args...)

# returns true if the Entry has metadata which should be included in the output
function has_output_metadata(entry::Entry, config::Config)
    output_metadata = false
    for m in keys(metadata(entry))
        if  m in config.metadata_order
            output_metadata = true
            break
        end
    end
    return output_metadata
end

# Convert's a string to a valid html id
function generate_html_id(s::AbstractString)
    # http://www.w3.org/TR/html4/types.html#type-id
    # ID tokens must begin with a letter ([A-Za-z]) and may be followed by any number of letters,
    # digits ([0-9]), hyphens ("-"), underscores ("_"), colons (":"), and periods (".").
    valid_chars = Set(Char['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M',
                           'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z',
                           'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm',
                           'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z',
                           '0', '1', '2', '3', '4', '5', '6', '7', '8', '9',
                           '-', '_', ':', '.'])

    replace_chars = Set(Char[' ', '"', '#', '$', '%', '&', '\'', '(', ')', '*', '+', ',',
                             '/', ';', '<', '=', '>', '?', '@', '[', ']', '^', '`', '{',
                             '¦', '}', '~'])
    skip_chars = Set(Char['!'])

    io = IOBuffer()
    for c in s
        if c in skip_chars
            continue
        elseif c in valid_chars
            write(io, lowercase(string(c)))
        elseif c in replace_chars
            write(io, "_")
        else
            write(io, string(@compat(Int(c))))
        end
    end
    # Note: In our case no need to check for begins with letter or is empty
    #  we prepend always the category
    return takebuf_string(io)
end

# Dispatch container for metadata display.
type Meta{keyword}
    content
end

# Cleanup object signatures. Remove method location links.

writeobj(any, entry)       = replace(string(any), ",", ", ")
writeobj(m::Method, entry) = replace(first(split(string(m), " at ")), ",", ", ")

function writeobj(f::Function, entry::Entry{:macro})
    replace(string("@", metadata(entry)[:signature]), ",", ", ")
end

function addentry!{category}(idx, obj, entry::Entry{category})
    section, pair = get!(idx, category, @compat(Tuple{AbstractString, Any})[]), (writeobj(obj, entry), obj)
    insert!(section, searchsortedlast(section, pair, by = x -> first(x)) + 1, pair)
end

# from base/methodshow.jl
function url(m::Meta{:source})
    line, file = m.content
    try
        d = dirname(file)
        u = Pkg.Git.readchomp(`config remote.origin.url`, dir=d)
        u = match(Pkg.Git.GITHUB_REGEX,u).captures[1]
        root = cd(d) do # dir=d confuses --show-toplevel, apparently
            Pkg.Git.readchomp(`rev-parse --show-toplevel`)
        end
        root = @windows? replace(root, "/", "\\") : root
        if startswith(file, root)
            commit = Pkg.Git.readchomp(`rev-parse HEAD`, dir=d)
            return "https://github.com/$u/tree/$commit/"*file[length(root)+2:end]*"#L$line"
        else
            return Base.fileurl(file)
        end
    catch
        return Base.fileurl(file)
    end
end

## Format-specific rendering ------------------------------------------------------------

include("render/plain.jl")
include("render/html.jl")
include("render/md.jl")


text:
input = [1,12,0,20,8,16]

struct Seen
    before::Int
    earlier::Int
end
Seen(index) = Seen(index, 0)

function remember(last, what, round)
    if !haskey(last, what)
        last[what] = Seen(round)
    else
        last[what] = Seen(round, last[what].before)
    end
end

function ex(last::Dict{Int,Seen}, which, round, current)
    while round <= which
        p = last[current]
        if p.earlier == 0 && p.before == round - 1
            current = 0
        else
            @assert p.earlier != 0
            current = p.before - p.earlier
        end
        remember(last, current, round)
        round += 1
    end
    current
end

function ex(nums, which)
    which <= length(nums) && return nums[which]

    d = Dict{Int,Seen}()
    for (idx, elm) in enumerate(nums)
        remember(d, elm, idx)
    end

    ex(d, which, length(nums) + 1, nums[end])
end

const YEAR = 2020
const HUGE = 30_000_000
println("The $(YEAR)th number spoken is: ", ex(input, YEAR))
println("The $(HUGE) number is: ", ex(input, HUGE))

using Test

@testset "Rambunctious Recitation" begin
    @testset "example 1" begin
        @test ex([1,3,2], YEAR) == 1
        @test ex([2,1,3], YEAR) == 10
        @test ex([1,2,3], YEAR) == 27
        @test ex([2,3,1], YEAR) == 78
        @test ex([3,2,1], YEAR) == 438
        @test ex([3,1,2], YEAR) == 1836
    end

    #@testset "example 2" begin
    #    @test ex([0,3,6], HUGE) == 175594
    #    @test ex([1,3,2], HUGE) == 2578
    #    @test ex([2,1,3], HUGE) == 3544142
    #    @test ex([1,2,3], HUGE) == 261214
    #    @test ex([2,3,1], HUGE) == 6895259
    #    @test ex([3,2,1], HUGE) == 18
    #    @test ex([3,1,2], HUGE) == 362
    #end

    @testset "results" begin
        @test ex(input, YEAR) == 273
        @test ex(input, HUGE) == 47205
    end
end

max_stars_repo_name: CmdQ/AoC


text:
using BrainFlow
using DataFrames
using VegaLite
using Query

BrainFlow.enable_dev_logger(BrainFlow.BOARD_CONTROLLER)

params = BrainFlowInputParams()
board_shim = BrainFlow.BoardShim(BrainFlow.GFORCE_PRO_BOARD, params)

BrainFlow.prepare_session(board_shim)
BrainFlow.start_stream(board_shim)
sleep(5)
BrainFlow.stop_stream(board_shim)
data = transpose(BrainFlow.get_board_data(board_shim))
BrainFlow.release_session(board_shim)

emg_channels = BrainFlow.get_emg_channels(BrainFlow.GFORCE_PRO_BOARD)
emg_data = data[:, emg_channels]

emg_names = Symbol.(["emg$x" for x in 1:8])
df = DataFrame(emg_data)
DataFrames.rename!(df, emg_names)

sampling_rate = BrainFlow.get_sampling_rate(BrainFlow.GFORCE_PRO_BOARD)
nrows = size(df)[1]
df.time = (1:nrows)/sampling_rate

df2 = stack(df, 1:8)
rename!(df2, [:time, :channel, :value])

df2 |> @filter(_.time > 0.2) |>
@vlplot(
    :line,
    x=:time,
    y=:value,
    row=:channel,
    width=400,
    height=25
)


text:
using Deconvolution
using Test, StableRNGs, FFTW

rng = StableRNG(42) # Fixed random seed

##### Wiener deconvolution

x = range(0, stop = 10, length = 15)
s = sinpi.(x) .- 1.5cos.(x)
n = rand(rng, length(s))
blurring_ft  = exp.(-0.001 * (range(-div(length(x), 2), length = length(x)) .^ 2) .^ (5 // 6))
blurred_s_ft = fftshift(blurring_ft) .* fft(s) .+ fft(n)
blurred_s    = real(ifft(blurred_s_ft))
blurring     = ifft(fftshift(blurring_ft))
@testset "No blurring with noise scalar parameter" begin
    @test wiener(s + n, s, 0.5) ≈
        [-1.0688616250172964, -0.32510758040709303, -0.4928689108095881, 1.644396339118468,
        1.866916504419297, 0.9245465020787589, 1.5461081346912358, 0.12677498594640266,
        -1.4157929359946475, -0.0247233772043046, -0.9345328207346199, 0.03628973013108201,
        2.414101393218902, 0.568743737331393, 1.5290604254094078]
end

@testset "Blurring with noise scalar parameter" begin
    @test wiener(blurred_s, s, 0.5, blurring) ≈
        [-1.0652058986210413, -0.3304166624246508, -0.4875855887501838, 1.6438301222753278,
        1.8625054116642619, 0.9303156760177144, 1.5413569689140112, 0.12849334093024348,
        -1.4141988631895932, -0.025180391629450346, -0.9343071603438413,
        0.03626923211984258, 2.4171155004554263, 0.5643293761857023, 1.5290637671248288]
end

@testset "No blurring with noise array" begin
    @test wiener(blurred_s, s, ifftshift(n)) ≈
        [-1.5060672134546482, -0.7416509948811038, -0.8868133826324668, 1.2108272657591024,
         1.414853150573584, 0.5580080918195482, 1.0351041993041612, -0.2750339384343115,
         -1.73737743125768, -0.5146485193142672, -1.3081993526690698, -0.37649294337342376,
         1.9491773714002598, 0.19435316342110073, 1.1019019675089878]
end

@testset "Blurring with noise array" begin
    @test wiener(blurred_s, s, ifftshift(n), blurring) ≈
        [-1.524937548601649, -0.7346444798363926, -0.8990054207462668, 1.2192518520141087,
        1.4207985817464395, 0.5541712011813319, 1.0415798812308625, -0.2742145640600878,
        -1.7491892004352156, -0.5071574425862486, -1.3165036261429737, -0.38315131571856315,
        1.9700562217921265, 0.18020610144859178, 1.1205661509755478]
end

##### Richardson-Lucy deconvolution

@testset "Richardson-Lucy deconvolution tests" begin
    k=-0.001
    blurring_ft  = exp.(-k * (range(-div(length(x), 2), length = length(x)) .^ 2) .^ (5 // 6))
    blurred_s_ft = fftshift(blurring_ft) .* fft(s)
    blurred_s    = real(ifft(blurred_s_ft))
    blurring     = ifft(fftshift(blurring_ft))

    estimated = lucy(blurred_s, blurring, iterations=1)

    @test sum(abs.(s .- estimated)) < sum(abs.(s .- blurred_s))

    k =-0.0015 # should improve also when blur params are not perfectly estimated
    blurring_ft  = exp.(-k * (range(-div(length(x), 2), length = length(x)) .^ 2) .^ (5 // 6))
    blurring     = ifft(fftshift(blurring_ft))

    estimated = lucy(blurred_s, blurring, iterations=1)

    @test sum(abs.(s .- estimated)) < sum(abs.(s .- blurred_s))
end


text:
"""
A partition divides a region into several subregions. The main algorithm of
every partition is deciding in which subregion a given point x lies, implemented
by the `partition_index` function.

The intent of a partition is to be combined with a function set on each subregion,
effectively creating a piecewise function.

The pieces of a partition have an index ranging from 1 to length(partition).
"""
abstract type Partition end

"""
A PartitionPiece represents a piece of a partition - it essentially combines
a partition and an index into that partition. It is constructed by indexing
a partition.
"""
struct PartitionPiece{P <: Partition,I}
    part    ::  P
    idx     ::  I
end

partition(piece::PartitionPiece) = piece.part
index(piece::PartitionPiece) = piece.idx

left(piece::PartitionPiece) = left(partition(piece), index(piece))
right(piece::PartitionPiece) = right(partition(piece), index(piece))

Base.in(x, piece::PartitionPiece) = in_partition(partition(piece), index(piece), x)

function checkbounds(part::Partition, i::Int)
    1 <= i <= length(part) || throw(BoundsError())
end

function getindex(part::Partition, i)
    checkbounds(part, i)
    PartitionPiece(part, i)
end

eachindex(part::Partition) = 1:length(part)

"""
A PiecewiseInterval is a subdivision of an interval into a number of subintervals.
The subintervals are determined by a set of breakpoints ``t_i``, ``i=0,...,n``. The main
interval is ``[t_0,t_n]``. The subintervals always contain the right endpoint, i.e.
the first subinterval is ``[t_0,t_1]``, the second is the halfopen interval
``(t_1,t_2]`` and so on until the n-th interval ``(t_{n-1},t_n]``.
"""
struct PiecewiseInterval{T} <: Partition
    points  ::  Array{T,1}
end

# n partitions of equal size on the interval [a,b]
PiecewiseInterval(a, b, n::Int) = PiecewiseInterval(collect(linspace(a,b,n+1)))

length(part::PiecewiseInterval) = length(part.points)-1

support(part::PiecewiseInterval) = Interval(part.points[1],part.points[end])

# Not sure this is a good idea, but let's respect the convention that the
    # interval does not contain the left endpoint, unless it is the first interval.
support(part::PiecewiseInterval, i::Int) = i==1 ? Interval(part.points[1],part.points[2]) : Interval{:open,:closed}(part.points[i],part.points[i+1])
left(part::PiecewiseInterval) = part.points[1]

left(part::PiecewiseInterval, i::Int) = part.points[i]

right(part::PiecewiseInterval) = part.points[end]

right(part::PiecewiseInterval, i::Int) = part.points[i+1]

"Return the index of the region in which the given point x lies."
function partition_index(part::PiecewiseInterval, x)
    @assert in(x,support(part))

    idx = 1
    while x > part.points[idx+1]
        idx += 1
    end
    idx
end

in_partition(part::PiecewiseInterval, i, x) = in(x,support(part,i))

# TODO: applying partition_index repeatedly, say for all points in a grid, could
# be done more efficiently if we can exploit ordering of the points in the grid.
# Perhaps partition_index applied to an iterable set of points could return an
# iterable set of indices?

split_interval(part::PiecewiseInterval, i::Int, x) = PiecewiseInterval([part.points[1:i]...,x,part.points[i+1:end]...])

split_interval(part::PiecewiseInterval, x) = split(part, partition_index(part, x), x)


text:
function kaplan_meier(sort_ev::SortedEvents)
    fs, ls, events = sort_ev.fs, sort_ev.ls, sort_ev.events
    ds = ls-fs+1
    ns = length(events)- fs +1
    surv = cumprod(1.-ds./ns)
    return getfield.(events[fs],[:time]), surv
end

kaplan_meier(events::Array) = kaplan_meier(SortedEvents(events))

if isdefined(:NullableArray)
    kaplan_meier(events::NullableArray, args...) = kaplan_meier(events.values[!events.isnull], args...)
end
kaplan_meier(events::AbstractArray, args...) = kaplan_meier(convert(Array,events), args...)


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

include("../scripts-part1/part1.jl")
H, xc_g = diffusion_3D(; do_visu=false)

# Reference test using ReferenceTests.jl
"Compare all dict entries"
comp(d1, d2) = keys(d1) == keys(d2) && all([ isapprox(v1, v2; atol = 1e-5) for (v1,v2) in zip(values(d1), values(d2))])
inds = Int.(ceil.(LinRange(1, length(xc_g), 12)))
d = Dict(:X=> xc_g[inds], :H=>H[inds, inds, 15])

@testset "Ref-test diffusion 3D dual-time" begin
    @test_reference "reftest-files/test_1.bson" d by=comp
end


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

# print a figure to a file
function printfigure(;handle::Handle    = gnuplot_state.current,
                     term::String       = config[:print][:print_term],
                     args...)

    # process arguments
    PA = PlotArgs()
    for k in keys(args)
        # substitute synonyms
        key = k
        for s in syn
            key ∈ s && (key = s[1]; break)
        end
        # substitute print arguments
        key == :font && (key = :printfont)
        key == :size && (key = :printsize)
        key == :linewidth && (key = :printlinewidth)
        key == :background && (key = :printbackground)
        # store arguments
        for f in fieldnames(PlotArgs)
            f == key && (setfield!(PA, f, string(args[k])); break)
        end
    end

    h = findfigure(handle)
    h == 0 && throw(DomainError(h, "requested figure does not exist."))
    isempty(PA.output) && throw(DomainError("Please specify an output filename."))

    # set figure's print parameters
    term == "pdf" && (term = "pdfcairo")
    term == "png" && (term = "pngcairo")
    term == "eps" && (term = "epscairo")
    valid_file_term(term)

    PA.printfont == "" && (printfont = TerminalDefaults[term][:font])
    PA.printsize == "" && (printsize = TerminalDefaults[term][:size])
    PA.printlinewidth == "" && (printlinewidth = "1")
    PA.printbackground == "" && (printbackground = TerminalDefaults[term][:background])

    fig = gnuplot_state.figs[h]
    pc = PrintConf(print_term       = term,
                   print_font       = PA.printfont,
                   print_size       = PA.printsize,
                   print_linewidth  = PA.printlinewidth,
                   print_background = PA.printbackground)
    fig.print = pc
    llplot(fig,print=true)

    return nothing
end


text:
__precompile__()
module NamedTuples

abstract NamedTuple <: Associative

Base.keys( t::NamedTuple ) = fieldnames( t )
Base.values( t::NamedTuple ) = [ getfield( t, i ) for i in fieldnames( t )]
Base.length( t::NamedTuple ) = length( fieldnames( t ))
# Iteration
Base.start( t::NamedTuple ) = 1
Base.done( t::NamedTuple, iter ) = iter>length( fieldnames( t ))
Base.next( t::NamedTuple, iter ) = ( getfield( t, iter ), iter + 1 )
Base.endof( t::NamedTuple ) = length( t )
Base.last( t::NamedTuple ) = t[end]
function Base.show( io::IO, t::NamedTuple )
    print(io, "(")
    first = true
    for (k,v) in zip(keys(t),values(t))
        !first && print(io, ", ")
        print(io, k, " => "); show(io, v)
        first = false
    end
    print(io, ")")
end
# Make this indexable so that it works like a Tuple
Base.getindex( t::NamedTuple, i::Int ) = getfield( t, i )
Base.getindex( t::NamedTuple, i::AbstractVector) = slice( t, i )
# We also support indexing by symbol
Base.getindex( t::NamedTuple, i::Symbol ) = getfield( t, i )
Base.getindex( t::NamedTuple, i::Symbol, default ) = get( t, i, default)
# This is a linear lookup...
Base.get( t::NamedTuple, i::Symbol, default ) = i in keys(t) ? t[i] : default
# Deep compare

import Base: ==

function ==( lhs::NamedTuple, rhs::NamedTuple)
    ( lhs === rhs ) && return true
    ( typeof( lhs ) != typeof( rhs )) && return false
    for i in 1:length( lhs )
        if ( lhs[i] != rhs[i])
            return false
        end
    end
    return true
end
# Deep hash

function Base.hash( nt::NamedTuple, hs::UInt64)
    h = 17
    for v in values(nt)
        h = h * 23 + hash( v, hs )
    end
    return h
end

# Helper type, for transforming parse tree to NameTuple definition
immutable ParseNode{T} end

function trans( ::Type{ParseNode{:Symbol}}, expr::Expr)
    (expr.args[1],nothing,nothing)
end

# Only escape items that need to be escaped.
escape( e::Expr ) = esc( e )
escape( e::Symbol ) = esc( e )
escape( e ) = e

function trans( ::Type{ParseNode{:(=>)}}, expr::Expr)
    (sym, typ ) = trans( expr.args[1])
    return (sym, typ, escape( expr.args[2] ))
end

# Allow unary
function trans( ::Type{ParseNode{:(::)}}, expr::Expr)
    if( length( expr.args ) > 1 )
        return ( expr.args[1], expr.args[2], nothing)
    else
        return ( nothing, expr.args[1], nothing)
    end
end

function trans( expr::Expr )
    trans( ParseNode{expr.head}, expr)
end

function trans( sym::Symbol )
    return (sym, nothing, nothing)
end

function trans( ::Type{ParseNode{:quote}}, expr::Expr )
    return trans( expr.args[1] )
end

# Literal nodes
function trans{T}( lit::T )
    return (nothing, nothing, escape(lit) )
end

function trans{T}( ::Type{ParseNode{T}}, expr::Expr)
    return (nothing, nothing, escape(expr) )
end

#
# Create a NameTuple in the context of this module
# this is only done if the tuple has not already been
# constructed.
function create_tuple( fields::Vector{Symbol})
    escaped_fieldnames = [replace(string(i), "_", "__") for i in fields]
    name = Symbol( string( "_NT_", join( escaped_fieldnames, "_")) )
    if !isdefined(NamedTuples, name)
        len = length( fields )
        types = [Symbol("T$n") for n in 1:len]
        tfields = [ Expr(:(::), Symbol( fields[n] ), Symbol( "T$n") ) for n in 1:len ]
        def = Expr(:type, false, Expr( :(<:), Expr( :curly, name, types... ), :NamedTuple ), Expr(:block, tfields...) )
        eval(NamedTuples, def)
    end
    return name
end

#
# Given a symbol list create the NamedTuple
#
@doc doc"Given a symbol vector create the `NamedTuple`" ->
function make_tuple( syms::Vector{Symbol} )
    return create_tuple( syms )
end

#
# Given an expression vector create the NamedTuple
#
@doc doc"Given an expression vector create the `NamedTuple`" ->
function make_tuple( exprs::Vector)
    len    = length( exprs )
    fields = Array(Symbol, len)
    values = Array(Any, len)
    typs   = Array(Any, len)

    # Are we directly constructing the type, if so all values must be
    # supplied by the caller, we use this state to ensure this
    construct = false
    # handle the case where this is defining a datatype
    for i in 1:len
        expr = exprs[i]
        ( sym, typ, val ) = trans( expr )
        if( construct == true && val == nothing || ( i > 1 && construct == false && val != nothing ))
            error( "Invalid tuple, all values must be specified during construction @ ($expr)")
        end
        construct  = val != nothing
        fields[i]  = sym != nothing?sym:Symbol( "_$(i)_")
        typs[i] = typ
        # On construction ensure that the types are consitent with the declared types, if applicable
        values[i]  = ( typ != nothing && construct)? Expr( :call, :convert, typ, val ) : val
    end

    name = create_tuple( fields )

    # Either call the constructor with the supplied values or return the type
    if( !construct )
        return Expr( :curly, name, typs... )
    else
        return Expr( :call, name, values ... )
    end
end

@doc doc"""
Syntax

    @NT( a, b )                 -> Defines a tuple with a and b as members
    @NT( a::Int64, b::Float64 ) -> Defines a tuple with the specific arg types as members
    @NT( a => 1, b => "hello")  -> Defines and constructs a tuple with the specifed members and values
    @NT( a, b )( 1, "hello")    -> Is equivalent to the above definition
    @NT( a::Int64 )( 2.0 )      -> Calls `convert( Int64, 2.0 )` on construction and sets `a`

NamedTuples may be used anywhere you would use a regular Tuple, this includes method definition and return values.

    module Test
    using NamedTuples

    function foo( y )
        a = 1
        x = 3
        return  @NT( a => 1, b => "world", c => "hello", d=>a/x, y => a/y  )
    end
    function bar( nt::@NT( a::Int64, c::ASCIIString ))
        return repeat( nt.c, nt.a )
    end

    end

    Test.foo( 1 ) # Returns a NamedTuple of 5 elements
    Test.bar( @NT( a=> 2, c=>"hello")) # Returns `hellohello`
""" ->
macro NT( expr... )
    return make_tuple( collect( expr ))
end

# Helper function for 0.4 compat
if VERSION < v"0.5.0"
    getfieldname( t, i ) = fieldnames(t)[i]
else
    getfieldname( t, i ) = fieldname( t, i )
end

@doc doc"""
Create a slice of an existing NamedTuple using a UnitRange. Construct a new NamedTuple with
the result.
This copies the underlying data.
""" ->
function Base.slice( t::NamedTuple, rng::AbstractVector )
    names = unique( Symbol[ isa(i,Symbol) ? i : getfieldname(typeof(t),i) for i in rng ] )
    name = create_tuple( names )
    getfield(NamedTuples,name)([ getfield( t, i ) for i in names ]...)
end

@doc doc"""
Merge two NamedTuples favoring the lhs
Order is preserved lhs names come first.
This copies the underlying data.
""" ->
function Base.merge( lhs::NamedTuple, rhs::NamedTuple )
    nms = unique( vcat( fieldnames( lhs ), fieldnames( rhs )) )
    name = create_tuple( nms )
    # FIXME should handle the type only case
    vals = [ haskey( lhs, nm ) ? lhs[nm] : rhs[nm] for nm in nms ]
    getfield(NamedTuples,name)(vals...)
end

@doc doc"""
Create a new NamedTuple with the new value set on it, either overwriting
the old value or appending a new value.
This copies the underlying data.
""" ->
function setindex{V}( t::NamedTuple, key::Symbol, val::V)
    nt = getfield( NamedTuples,create_tuple( [key] ))( val )
    return merge( t, nt )
end

@doc doc"""
Create a new NamedTuple with the specified element removed.
""" ->
function delete( t::NamedTuple, key::Symbol )
    nms = filter( x->x!=key, fieldnames( t ) )
    name = create_tuple( nms )
    vals = [ getindex( t, nm ) for nm in nms ]
    return getfield(NamedTuples, name)(vals...)
end

export @NT, NamedTuple, setindex, delete

end # module


text:
using DataFrames, Random, VegaLite, SimpleWeightedGraphs, LinearAlgebra
using VegaGraphs
using Test

@testset "graphplot()" begin
    Random.seed!(7)
    incidence = rand([0,1],20,20);
    m         = incidence'*incidence
    m[diagind(m)] .= 0.0
    g = SimpleWeightedGraph(m)
    field = rand([100,200,300],20)
    p     = VegaGraphs.graphplot(g,tooltip=true,
        node_label=false,node_colorfield=field,
        node_colorfieldtype="o",node_sizefield=field)
    @test typeof(p) == VegaLite.VLSpec
end


text:
# runtime.jl

using SPECTrecon

include("helper.jl")
include("rotate2.jl")
include("fft_convolve.jl")
include("SPECTplan.jl")
include("project.jl")
include("backproject.jl")


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

mutable struct VirtualMachineScaleSetUpdatePublicIPAddressConfiguration <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    properties::Any # spec type: Union{ Nothing, VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties } # spec name: properties

    function VirtualMachineScaleSetUpdatePublicIPAddressConfiguration(;name=nothing, properties=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetUpdatePublicIPAddressConfiguration, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(VirtualMachineScaleSetUpdatePublicIPAddressConfiguration, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type VirtualMachineScaleSetUpdatePublicIPAddressConfiguration

const _property_map_VirtualMachineScaleSetUpdatePublicIPAddressConfiguration = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("properties")=>Symbol("properties"))
const _property_types_VirtualMachineScaleSetUpdatePublicIPAddressConfiguration = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("properties")=>"VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties")
Base.propertynames(::Type{ VirtualMachineScaleSetUpdatePublicIPAddressConfiguration }) = collect(keys(_property_map_VirtualMachineScaleSetUpdatePublicIPAddressConfiguration))
Swagger.property_type(::Type{ VirtualMachineScaleSetUpdatePublicIPAddressConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetUpdatePublicIPAddressConfiguration[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetUpdatePublicIPAddressConfiguration }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetUpdatePublicIPAddressConfiguration[property_name]

function check_required(o::VirtualMachineScaleSetUpdatePublicIPAddressConfiguration)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetUpdatePublicIPAddressConfiguration }, name::Symbol, val)
end


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

initial_condition = initial_condition_convergence_test
source_terms = source_terms_convergence_test
boundary_conditions = boundary_condition_periodic

###############################################################################
# Get the DG approximation space

solver = DGSEM(polydeg=6, surface_flux=FluxRotated(flux_hll))

###############################################################################
# Get the curved quad mesh from a file (downloads the file if not available locally)

default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh")
isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh",
                                       default_mesh_file)
mesh_file = default_mesh_file

mesh = UnstructuredMesh2D(mesh_file, periodicity=true)

###############################################################################
# create the semi discretization object

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    source_terms=source_terms,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

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

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)

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

sol = solve(ode, SSPRK43(), save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
using Documenter, NCBITaxonomy

makedocs(
    sitename="NCBITaxonomy",
    authors="Timothée Poisot",
    modules=[NCBITaxonomy],
    pages=[
        "Index" => "index.md",
        "Finding taxa" => "namefinding.md",
        "Navigating lineages" => "lineages.md",
        "Portal data use-case" => "portal.md",
        "Phylo.jl use-case" => "phylo.md"
        ]
)

deploydocs(
    deps=Deps.pip("pygments", "python-markdown-math"),
    repo="github.com/EcoJulia/NCBITaxonomy.jl.git",
    devbranch="main",
    push_preview=true
)


text:
using LocalPackageServer, LocalRegistry, Pkg, Test, Inflate
using LocalPackageServer: fetch, GitStorageServer, Config

const TEST_GITCONFIG = Dict(
    "user.name" => "LocalRegistryTests",
    "user.email" => "localregistrytests@example.com"
)

# Hook into the LocalRegistry testing infrastructure.
include(joinpath(dirname(dirname(pathof(LocalRegistry))), "test", "utils.jl"))

mktempdir(@__DIR__) do test_dir
    # Start by creating an empty registry.
    upstream_registry_dir = joinpath(test_dir, "upstream_registry.git")
    mkpath(upstream_registry_dir)
    run(`git -C $(upstream_registry_dir) init --bare`)
    registry_dir = joinpath(test_dir, "TestRegistry")
    packages_dir = joinpath(test_dir, "packages")
    upstream_registry_url = string("file://", upstream_registry_dir)
    registry_uuid = "ed6ca2f6-392d-11ea-3224-d3daf7fee369"
    create_registry(registry_dir, upstream_registry_url, push = true,
                    gitconfig = TEST_GITCONFIG,
                    uuid = registry_uuid)

    # Next configure a git backed StorageServer, serving the newly
    # created registry.
    config_dict = Dict{String, Any}()
    config_dict["host"] = "127.0.0.1"
    config_dict["port"] = 8080
    config_dict["local_registry"] = upstream_registry_url
    config_dict["cache_dir"] = joinpath(test_dir, "cache")
    config_dict["git_clones_dir"] = joinpath(test_dir, "data")
    config_dict["min_time_between_registry_updates"] = 0
    config_dict["gitconfig"] = TEST_GITCONFIG
    config = Config(config_dict)
    mkpath(config.cache_dir)

    # We won't actually run this server, but test its backend
    # functions, specifically through `fetch`.
    path = fetch(config, "/registries")
    @test isfile(path)
    registry_git = gitcmd(registry_dir, TEST_GITCONFIG)
    hash = readchomp(`$(registry_git) rev-parse --verify HEAD:`)
    initial_registry_resource = readchomp(path)
    @test initial_registry_resource == "/registry/$(registry_uuid)/$(hash)"
    path = fetch(config, initial_registry_resource)
    @test isfile(path)
    @test startswith(inflate_gzip(path), "Registry.toml")

    # Create a test package and register it.
    prepare_package(packages_dir, "FirstTest1.toml")
    first_test_dir = joinpath(packages_dir, "FirstTest")
    first_test_url = "file://$(first_test_dir)"
    register(first_test_dir, registry = registry_dir,
             repo = first_test_url,
             gitconfig = TEST_GITCONFIG, push = true)
    first_test_uuid = "d7508571-2240-4c50-b21c-240e414cc6d2"

    # This should now give a different hash for the registry.
    path = fetch(config, "/registries")
    @test isfile(path)
    @test readchomp(path) != initial_registry_resource
    path = fetch(config, readchomp(path))
    @test isfile(path)

    # Verify that the package resource is available.
    git = gitcmd(first_test_dir, TEST_GITCONFIG)
    hash = readchomp(`$git rev-parse --verify HEAD:`)
    path = fetch(config, "/package/$(first_test_uuid)/$(hash)")
    @test isfile(path)
    dir = joinpath(test_dir, "data", "packages", first_test_uuid)
    @test isdir(dir)
    git = gitcmd(dir, TEST_GITCONFIG)
    @test readchomp(`$git rev-parse --verify HEAD:`) == hash

    # Test the server metadata.
    meta = LocalPackageServer.collect_meta(config)
    @test meta["julia_version"] == string(VERSION)
    @test haskey(meta, "pkgserver_version")
    @test meta["packages_cached"] == 1
    @test meta["artifacts_cached"] == 0
end


text:
using GeometricIntegrators
using Weave

GeometricIntegrators.set_config(:verbosity, 0)

Weave.set_chunk_defaults!(:echo => false, :results => "raw")

include("../src/massless-charged-particle.jl")

weave("../weave/massless-charged-particle-spark-glvprk.jmd",
         out_path = "src/massless-charged-particle",
         doctype = "github",
         mod = MasslessChargedParticleSPARK)

weave("../weave/massless-charged-particle-spark-glspark.jmd",
         out_path = "src/massless-charged-particle",
         doctype = "github",
         mod = MasslessChargedParticleSPARK)

weave("../weave/massless-charged-particle-spark-lobspark.jmd",
         out_path = "src/massless-charged-particle",
         doctype = "github",
         mod = MasslessChargedParticleSPARK)

weave("../weave/massless-charged-particle-spark-pinternal.jmd",
         out_path = "src/massless-charged-particle",
         doctype = "github",
         mod = MasslessChargedParticleSPARK)

weave("../weave/massless-charged-particle-spark-plobatto.jmd",
         out_path = "src/massless-charged-particle",
         doctype = "github",
         mod = MasslessChargedParticleSPARK)

weave("../weave/massless-charged-particle-spark-pmidpoint.jmd",
         out_path = "src/massless-charged-particle",
         doctype = "github",
         mod = MasslessChargedParticleSPARK)

weave("../weave/massless-charged-particle-spark-psymmetric.jmd",
         out_path = "src/massless-charged-particle",
         doctype = "github",
         mod = MasslessChargedParticleSPARK)


text:
module SolveDSGE
using LinearAlgebra

include("types.jl")
include("extra_functions.jl")
include("vec_to_vech.jl")
include("vech_to_vec.jl")
include("derivative.jl")
include("hessian.jl")
include("newton.jl")
include("doubling.jl")
include("dlyap.jl")
include("check_model.jl")
include("solve_re.jl")
include("solve_disc.jl")
include("solve_commit.jl")
include("solve_quasi.jl")
include("solve_timeless.jl")
include("impulses.jl")

export derivative,
       hessian,
       newton,
       solve_re,
       check_model_form,
       convert_second_order

export Blanchard_Kahn_Form,
       Blanchard_Kahn_Soln,
       Klein_Form,
       Klein_Soln,
       Binder_Pesaran_Form,
       Binder_Pesaran_Soln,
       Sims_Form,
       Sims_Soln,
       Gomme_Klein_Form,
       Gomme_Klein_Soln,
       Lombardo_Sutherland_Form,
       Lombardo_Sutherland_Soln

export doubling,
       dlyap,
       solve_disc,
       solve_commit,
       solve_quasi,
       solve_timeless

export State_Space_Objective,
       Structural_Objective

export State_Space_Form,
       Generalized_State_Space_Form,
       Structural_Form,
       Generalized_Structural_Form,
       State_Space_Soln,
       Structural_Soln

export Perturbable_Soln,
       Second_Order_State_Space_Soln

export impulses

end


text:
# Unit test functions
function unit_test_dynamicarray_empty()
    d = DynamicArray{Float64}()
    @test d.capacity == 0

    return nothing
end

function unit_test_dynamicarray_size()
    d = DynamicArray{Float64}()
    @test size(d) == 0

    return nothing
end

function unit_test_dynamicarray_isempty()
    d = DynamicArray{Float64}()
    @test isempty(d)

    return nothing
end

function unit_test_dynamicarray_push()
    d = DynamicArray{Float64}()
    val = 1.0
    push!(d, val)
    @test d[1] == val

    return nothing
end

function unit_test_dynamicarray_pop()
    d = DynamicArray{Float64}()
    val = 1.0
    push!(d, val)
    @test pop!(d) == val

    return nothing
end

function unit_test_dynamicarray_getindex()
    d = DynamicArray{Float64}()
    val = 2.0
    push!(d, val)
    @test d[1] == val

    return nothing
end

function unit_test_dynamicarray_setindex()
    d = DynamicArray{Float64}()
    val1 = 2.0
    val2 = 3.0
    push!(d, val1)
    d[1] = val2
    @test d[1] == val2

    return nothing
end

function unit_test_dynamicarray_resize()
    d = DynamicArray{Float64}()
    capacity = 8
    resize!(d, capacity)
    @test d.capacity == capacity

    return nothing
end

function unit_test_dynamicarray_out_of_bounds()
    d = DynamicArray{Float64}()
    val = 2.0
    push!(d, val)
    @test_throws BoundsError d[2]

    return nothing
end

# Run the unit tests
@testset "Dynamic array" begin
    unit_test_dynamicarray_empty()
    unit_test_dynamicarray_size()
    unit_test_dynamicarray_isempty()
    unit_test_dynamicarray_push()
    unit_test_dynamicarray_pop()
    unit_test_dynamicarray_getindex()
    unit_test_dynamicarray_setindex()
    unit_test_dynamicarray_resize()
    unit_test_dynamicarray_out_of_bounds()
end


text:
using FEASTSolver
using LinearAlgebra
using NonlinearEigenproblems: nep_gallery, compute_Mder

function info(Λ, X, residuals, c, r)
    inside(x) = in_contour(x, c, r)
    in_eig = Λ[inside.(Λ)]
    in_res = residuals[inside.(Λ)]
    print("number inside : ")
    println(size(Λ[inside.(Λ)])[1])
    if sum(inside.(Λ)) > 0
        in_res_conv = in_res[in_res .<= 1e-5]
        in_eig_conv = in_eig[in_res .<= 1e-5]
        print("number inside converged : ")
        println(size(in_eig_conv)[1])
        print("max res inside: ")
        println(maximum(residuals[inside.(Λ)]))
        if size(in_res_conv, 1) > 0
            print("max res inside non spurious: ")
            println(maximum(in_res_conv))
        end
    end
    # display(in_eig_conv)
end

nep = nep_gallery("nlevp_native_hadeler", 100, 200)
T(x) = compute_Mder(nep, x)

n = size(nep, 1)

#
# function hadeler(n=200, b₀=100)
#     M0 = zeros(ComplexF64, n, n)
#     M1 = zeros(ComplexF64, n, n)
#     M2 = zeros(ComplexF64, n, n)
#
#     M0 .= b₀ * Matrix(I, n, n)
#     for i=1:n
#         for j=1:n
#             M1[i,j] = (n+1-max(i,j))*i*j
#             M2[i,j] = (n*(i == j)) + 1/(i+j)
#         end
#     end
#
#     function prob(B0, B1, B2)
#         f(z) =  (exp(z)-1).*B1 .+ z^2 .* B2 .- B0
#     end
#
#     return prob(M0, M1, M2)
# end
#
# C = complex(-3.5, 0.0)
# R = 0.05
C = complex(-30.0, 0.0)
R = 10.0
# T = hadeler()
# e, v, res = nlfeast!(T, rand(ComplexF64,n,10), 2^3, 0, c=C, r=R)
e, v, res = nlfeast!(T, (rand(ComplexF64,n,15)), 2^3, 30, c=C, r=R, debug=true, ϵ=10e-16, spurious=1e-2)
# e, v, res = nlfeast_moments!(T, (rand(ComplexF64,1000,20)), 2^6, 10,store=false, moments=2, c=C, r=R, debug=true, ϵ=10e-16, spurious=1e-6)
# e, v, res = @profile nlfeast!(T, rand(ComplexF64,64,30), 2^4, 20, c=complex(1.0,1.0), r=0.5, debug=true, ϵ=10e-16)
# Profile.print()
# e, v, res = nlfeast_it!(T, rand(ComplexF64,64,30), 2^3, 0, c=complex(1.0,1.0), r=0.5)
# e, v, res = @timev nlfeast_it!(T, rand(ComplexF64,64,30), 2^3, 20, c=complex(1.0,1.0), r=0.5, debug=true, ϵ=10e-14)
# display(e)
info(e,v, res, C, R)

inside(x) = in_contour(x, C, R)
in_eig = e[inside.(e)]
in_res = res[inside.(e)]
# in_res_conv = in_res[in_res .<= 1e-3]
# in_eig_conv = in_eig[in_res .<= 1e-3]
# print("\nmax res inside: ")
# println(maximum(res[inside.(e)]))
# print("max res inside non spurious: ")
# println(maximum(in_res_conv))
# print("number inside : ")
# println(size(e[inside.(e)])[1])
# print("number inside converged : ")
# println(size(in_eig_conv)[1])
display(in_eig)
display(in_res)


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

@testset "MIT,Duration" begin
    # mit2yp conversions
    @test mit2yp(MIT{Quarterly}(5)) == (1, 2)
    @test mit2yp(MIT{Quarterly}(4)) == (1, 1)
    @test mit2yp(MIT{Quarterly}(3)) == (0, 4)
    @test mit2yp(MIT{Quarterly}(2)) == (0, 3)
    @test mit2yp(MIT{Quarterly}(1)) == (0, 2)
    @test mit2yp(MIT{Quarterly}(0)) == (0, 1)
    @test mit2yp(MIT{Quarterly}(-1)) == (-1, 4)
    @test mit2yp(MIT{Quarterly}(-2)) == (-1, 3)
    @test mit2yp(MIT{Quarterly}(-3)) == (-1, 2)
    @test mit2yp(MIT{Quarterly}(-4)) == (-1, 1)
    @test mit2yp(MIT{Quarterly}(-5)) == (-2, 4)
    @test mit2yp(MIT{Quarterly}(-6)) == (-2, 3)
    # subtractions
    @test typeof(qq(2020, 1) - qq(2019, 2)) == Duration{Quarterly}
    @test typeof(qq(2020, 1) - 2) == MIT{Quarterly}
    @test typeof(qq(2020, 1) - Duration{Quarterly}(2)) == MIT{Quarterly}
    @test typeof(Duration{Quarterly}(5) - 2) == Duration{Quarterly}
    @test typeof(Duration{Quarterly}(5) - Duration{Quarterly}(2)) == Duration{Quarterly}
    @test_throws ArgumentError qq(2020, 1) - mm(2019, 2)
    @test_throws ArgumentError qq(2020, 1) - Duration{Monthly}(5)
    @test_throws ArgumentError Duration{Quarterly}(8) - Duration{Monthly}(5)
    # equality
    @test qq(2020, 1) == qq(2020, 1)
    @test qq(2020, 1) != qq(2020, 2)
    @test qq(2020, 1) != mm(2020, 1)
    @test 5 == qq(2020, 1) - (qq(2020, 1) - 5)
    @test Duration{Quarterly}(5) == qq(2020, 1) - (qq(2020, 1) - 5)
    @test 5 == qq(2020, 1) - (qq(2020, 1) - 5)
    @test Duration{Monthly}(5) != qq(2020, 1) - (qq(2020, 1) - 5)
    @test Duration{Quarterly}(5) != MIT{Quarterly}(5)
    @test 5 == MIT{Quarterly}(5)
    # order
    @test qq(2000, 1) < qq(2000, 2)
    @test qq(2000, 1) <= qq(2000, 1)
    @test_throws ArgumentError qq(2000, 1) < mm(2000, 1)
    @test qq(0, 1) == 0
    @test mm(0, 1) == 0
    @test qq(0, 1) != mm(0, 1)
    @test_throws ArgumentError qq(0, 1) <= mm(0, 1)
    @test_throws ArgumentError qq(0, 1) < mm(0, 1)
    @test Duration{Quarterly}(5) < Duration{Quarterly}(6)
    @test !(Duration{Quarterly}(5) < Duration{Quarterly}(5))
    @test Duration{Quarterly}(5) <= Duration{Quarterly}(5)
    @test Duration{Quarterly}(5) == 5
    @test Duration{Monthly}(5) == 5
    @test !(Duration{Quarterly}(5) == Duration{Monthly}(5))
    @test_throws ArgumentError Duration{Quarterly}(5) < Duration{Monthly}(5)
    @test_throws ArgumentError Duration{Quarterly}(5) <= Duration{Monthly}(5)
    @test (MIT{Quarterly}(5) == 5) && (Duration{Quarterly}(5) == 5) && (MIT{Quarterly}(5) != Duration{Quarterly}(5))
    @test_throws ArgumentError MIT{Quarterly}(5) < Duration{Quarterly}(5)
    @test_throws ArgumentError MIT{Quarterly}(5) <= Duration{Quarterly}(5)
    @test_throws ArgumentError Duration{Quarterly}(5) < MIT{Quarterly}(5)
    @test_throws ArgumentError Duration{Quarterly}(5) <= MIT{Quarterly}(5)
    # addition
    @test qq(2020, 1) + 4 == qq(2021, 1)
    @test_throws ArgumentError qq(2020, 1) + qq(1, 0)
    @test_throws ArgumentError qq(2020, 1) + mm(1, 0)
    @test qq(2020, 1) + Duration{Quarterly}(4) == qq(2021, 1)
    @test Duration{Quarterly}(5) + Duration{Quarterly}(2) == 7
    @test Duration{Quarterly}(5) + 2 == 7
    @test Duration{Quarterly}(5) + 2 isa Duration{Quarterly}
    @test 2 + Duration{Quarterly}(5) == 7
    @test 2 + Duration{Quarterly}(5) isa Duration{Quarterly}
    @test_throws ArgumentError Duration{Quarterly}(5) + Duration{Monthly}(2)
    @test Duration{Quarterly}(5) + Duration{Quarterly}(2) isa Duration{Quarterly}
    @test_throws ArgumentError 20Q1 + Duration{Monthly}(2)
    # conversions to float (for plotting)
    @test 2000Q1 + 1 == 2000Q2
    @test 2000Q1 + 1.0 == 2001.0
    @test 2000Q1 + 1.2 == 2001.2
    @test 1.2 + 5U == 6.2
    @test 5U + 1.2 == 6.2
    # promotions
    @test_throws ArgumentError promote(1, 1Q1)
    @test_throws ArgumentError promote(1Q1, 1)
    @test_throws ArgumentError promote(1Q1 - 1Q1, 1)
    @test_throws ArgumentError promote(1, 1Q1 - 1Q1)
    @test promote(1.1, 1Q1) === (1.1, 1.0)
    @test promote(1Q1, 1.2) === (1.0, 1.2)
end

@testset "Range" begin
    rng = 2020Q1:2020Q4
    @test rng isa UnitRange{MIT{Quarterly}}
    @test isempty(2020Q1:2019Q1)
    @test length(rng) isa Int
    @test length(rng) == 4
    @test step(rng) isa Int
    @test step(rng) == 1
    for (i, m) in enumerate(rng)
    @test m isa MIT{Quarterly}
        @test first(rng) <= m <= last(rng)
        @test rng[i] == m
    end
    @test_throws ArgumentError 2020Q1:2020M12
    @test union(3U:5U, 4U:6U) === 3U:6U
    @test_throws ArgumentError union(3U:5U, 4Q1:6Q1)
end

@testset "FPConst" begin
    @test 8U === MIT{Unit}(8)
    @test 2000Y === yy(2000)
    @test 1999Q1 === qq(1999, 1)
    @test 1999Q2 === qq(1999, 2)
    @test 1999Q3 === qq(1999, 3)
    @test 1999Q4 === qq(1999, 4)
    @test 1988M12 === mm(1988, 12)
    @test 1988M11 === mm(1988, 11)
    @test 1988M10 === mm(1988, 10)
    @test 1988M9 === mm(1988, 9)
    @test 1988M8 === mm(1988, 8)
    @test 1988M7 === mm(1988, 7)
    @test 1988M6 === mm(1988, 6)
    @test 1988M5 === mm(1988, 5)
    @test 1988M4 === mm(1988, 4)
    @test 1988M3 === mm(1988, 3)
    @test 1988M2 === mm(1988, 2)
    @test 1988M1 === mm(1988, 1)
end

@testset "MITops" begin
    @test 5U < 8U
    @test 5U <= 8U
    @test 5U <= 5U
    @test 5U >= 5U
    @test 5U == 5U
    @test 8U >= 5U
    @test 8U > 5U
    @test 2001Q1 >= 2000Q3
    @test_throws ArgumentError 1M1 > 2Q1
    @test_throws ArgumentError 1M1 <= 2Q1
    @test 1M1 != 2Q1

    @test 2001Y + 5 == 2006Y
    @test 6 + 2001Q3 == 2003Q1
    @test 2003Q1 - 2001Q3 == 6
    @test 2003Q1 - 6 == 2001Q3
    @test_throws ArgumentError 6 - 2003Q1
    @test_throws ArgumentError 2003Q1 + 2003Q1
    @test_throws ArgumentError 2003Q1 + 2003Y
end

@testset "MIT.show" begin
    let io = IOBuffer()
        show(io, 2020Q1)
        show(io, pp(20, 3; N=6))
        show(io, 5U)
        println(io, 2020M1 - 2019M1)
        show(io, 3U - 2U)
        show(io, 2000M12 - 2000M1)
        println(io, Duration{Yearly}(7))
        show(io, Q1)
        show(io, 1U)
        println(io, M1, M12, ".")
        foo = readlines(seek(io, 0))
        @test foo == ["2020Q120P35U12", "1117", "1Q11U1M11M12."]
    end
end

@testset "frequencyof" begin
    @test frequencyof(qq(2000, 1)) == Quarterly
    @test frequencyof(mm(2000, 1)) == Monthly
    @test frequencyof(yy(2000)) == Yearly
    @test frequencyof(1U) == Unit
    @test frequencyof(qq(2001, 1):qq(2002, 1)) == Quarterly
    @test_throws ArgumentError frequencyof(1)
    @test_throws ArgumentError frequencyof(Int)
    @test frequencyof(qq(2000, 1) - qq(2000, 1)) == Quarterly
    @test frequencyof(mm(2000, 1) - mm(2000, 1)) == Monthly
    @test frequencyof(yy(2000, 1) - yy(2000, 1)) == Yearly
    @test frequencyof(5U - 3U) == Unit
    @test frequencyof(TSeries(yy(2000), zeros(5))) == Yearly
end

@testset "mm, qq, yy" begin
    @test mm(2020, 1) == MIT{Monthly}(2020 * 12)
    @test qq(2020, 1) == MIT{Quarterly}(2020 * 4)
    @test yy(2020) == MIT{Yearly}(2020)
end

@testset "year, period" begin
    @test_throws ArgumentError year(1U)
    let val = pp(2020, 2; N=6)
        @test year(val) == 2020
        @test period(val) == 2
        @test frequencyof(val) === YPFrequency{6}
    end
    @test year(mm(2020, 12)) == 2020
    @test period(mm(2020, 12)) == 12
end


text:
using Test
using BigBed
using Distributions
using GenomicFeatures

import Random
import ColorTypes: RGB
import FixedPointNumbers: N0f8

using BioCore

using FormatSpecimens

# Generate an array of n random Interval{Int} object.
# With sequence names samples from seqnames, and intervals drawn to lie in [1, maxpos].
function random_intervals(seqnames, maxpos::Int, n::Int)
    seq_dist = Categorical(length(seqnames))
    strand_dist = Categorical(2)
    length_dist = Normal(1000, 1000)
    intervals = Vector{Interval{Int}}(undef, n)
    for i in 1:n
        intlen = maxpos
        while intlen >= maxpos || intlen <= 0
            intlen = ceil(Int, rand(length_dist))
        end
        first = rand(1:maxpos-intlen)
        last = first + intlen - 1
        strand = rand(strand_dist) == 1 ? STRAND_POS : STRAND_NEG
        intervals[i] = Interval{Int}(seqnames[rand(seq_dist)], first, last, strand, i)
    end
    return intervals
end

@testset "BigBed" begin
    @testset "empty" begin
        buffer = IOBuffer()
        data = buffer.data
        writer = BigBed.Writer(buffer, [("chr1", 1000)])
        close(writer)
        reader = BigBed.Reader(IOBuffer(data))
        @test length(collect(reader)) == 0
    end

    @testset "small" begin
        buffer = IOBuffer()
        data = buffer.data
        writer = BigBed.Writer(buffer, [("chr1", 1000)])
        write(writer, ("chr1", 50, 100, "name1"))
        close(writer)
        reader = BigBed.Reader(IOBuffer(data))
        records = collect(reader)
        @test length(records) == 1
        @test BigBed.chrom(records[1]) == "chr1"
        @test BigBed.chromstart(records[1]) === 50
        @test BigBed.chromend(records[1]) === 100
        @test BigBed.name(records[1]) == "name1"
        @test !BigBed.hasscore(records[1])
        @test BigBed.optionals(records[1]) == ["name1"]

        buffer = IOBuffer()
        data = buffer.data
        writer = BigBed.Writer(buffer, [("chr1", 1000)])
        write(writer, ("chr1", 1, 100, "some name", 100, '+', 10, 90, RGB(0.5, 0.1, 0.2), 2, [4, 10], [10, 20]))
        close(writer)
        reader = BigBed.Reader(IOBuffer(data))
        records = collect(reader)
        @test length(records) == 1
        @test BigBed.haschrom(records[1]) === BioCore.hasseqname(records[1]) === true
        @test BigBed.chrom(records[1]) == BioCore.seqname(records[1]) == "chr1"
        @test BigBed.haschromstart(records[1]) === BioCore.hasleftposition(records[1]) === true
        @test BigBed.chromstart(records[1]) === BioCore.leftposition(records[1]) === 1
        @test BigBed.haschromend(records[1]) === BioCore.hasrightposition(records[1]) === true
        @test BigBed.chromend(records[1]) === BioCore.rightposition(records[1]) === 100
        @test BigBed.hasname(records[1])
        @test BigBed.name(records[1]) == "some name"
        @test BigBed.hasscore(records[1])
        @test BigBed.score(records[1]) === 100
        @test BigBed.hasstrand(records[1])
        @test BigBed.strand(records[1]) === STRAND_POS
        @test BigBed.hasthickstart(records[1])
        @test BigBed.thickstart(records[1]) === 10
        @test BigBed.hasthickend(records[1])
        @test BigBed.thickend(records[1]) === 90
        @test BigBed.hasitemrgb(records[1])
        @test BigBed.itemrgb(records[1]) === convert(RGB{N0f8}, RGB(0.5, 0.1, 0.2))
        @test BigBed.hasblockcount(records[1])
        @test BigBed.blockcount(records[1]) === 2
        @test BigBed.hasblocksizes(records[1])
        @test BigBed.blocksizes(records[1]) == [4, 10]
        @test BigBed.hasblockstarts(records[1])
        @test BigBed.blockstarts(records[1]) == [10, 20]
        @test BigBed.optionals(records[1]) == ["some name", "100", "+", "9", "90", "128,26,51", "2", "4,10,", "9,19,"]
    end

    @testset "large" begin
        buffer = IOBuffer()
        data = buffer.data
        binsize = 32
        writer = BigBed.Writer(buffer, [("chr1", 100_000), ("chr2", 90_000)], binsize=binsize)
        for i in 1:10_000
            write(writer, ("chr1", (i-1)*10+1, i*10, string("name", i)))
        end
        n = 0
        p = 1
        while p ≤ 90_000
            sz = min(rand(1:100), 90_000 - p)
            write(writer, ("chr2", p, p+sz, string("name", n + 1)))
            n += 1
            p += sz + 1
        end
        close(writer)
        reader = BigBed.Reader(IOBuffer(data))
        records = collect(reader)
        @test length(records) == 10_000 + n
        records = collect(eachoverlap(reader, Interval("chr1", 50_001, 50_165)))
        @test length(records) == 17
    end

    @testset "round trip" begin
        function test_round_trip(filepath)
            reader = open(BigBed.Reader, filepath)
            buffer = IOBuffer()
            data = buffer.data
            writer = BigBed.Writer(buffer, BigBed.chromlist(reader))
            original = []
            for record in reader
                t = (BigBed.chrom(record), BigBed.chromstart(record), BigBed.chromend(record), BigBed.optionals(record)...)
                write(writer, t)
                push!(original, t)
            end
            close(writer)
            close(reader)

            reader = BigBed.Reader(IOBuffer(data))
            copy = []
            for record in reader
                t = (BigBed.chrom(record), BigBed.chromstart(record), BigBed.chromend(record), BigBed.optionals(record)...)
                push!(copy, t)
            end
            close(reader)

            @test original == copy
        end

        dir = path_of_format("BBI")
        for specimen in list_valid_specimens("BBI")
            if hastag(specimen, "bigbed")
                test_round_trip(joinpath(dir, specimen["filename"]))
            end
        end
    end

    @testset "overlap" begin
        chromlen = 1_000_000
        Random.seed!(1234)
        chroms = ["one", "two", "three", "four", "five"]
        intervals = IntervalCollection([Interval(GenomicFeatures.seqname(i), GenomicFeatures.leftposition(i), GenomicFeatures.rightposition(i)) for i in random_intervals(chroms, chromlen, 10_000)], true)

        buffer = IOBuffer()
        data = buffer.data
        writer = BigBed.Writer(buffer, [(chrom, chromlen) for chrom in chroms])
        for i in intervals
            write(writer, i)
        end
        close(writer)

        reader = BigBed.Reader(IOBuffer(data))
        queries = random_intervals(chroms, chromlen, 1000)
        triplet(x::Interval) = GenomicFeatures.seqname(x), GenomicFeatures.leftposition(x), GenomicFeatures.rightposition(x)
        triplet(x::BigBed.Record) = BigBed.chrom(x), BigBed.chromstart(x), BigBed.chromend(x)
        @test all(triplet.(collect(eachoverlap(intervals, q))) == triplet.(collect(eachoverlap(reader, q))) for q in queries)
        close(reader)
    end
end


text:
@test_skip solve
@test_skip φΩ
@test_skip φΣ
@test_skip rfenergy


text:
export linear_region, linear_regions, estimate_boxsizes, linreg
export boxcounting_dim, capacity_dim, generalized_dim,
information_dim, estimate_boxsizes, kaplanyorke_dim
#####################################################################################
# Functions and methods to deduce linear scaling regions
#####################################################################################
using Statistics
using Statistics: covm, varm
# The following function comes from a version in StatsBase that is now deleted
# StatsBase is copyrighted under the MIT License with
# Copyright (c) 2012-2016: Dahua Lin, Simon Byrne, Andreas Noack, Douglas Bates,
# John Myles White, Simon Kornblith, and other contributors.
"""
    linreg(x, y) -> a, b
Perform a linear regression to find the best coefficients so that the curve:
`y = a + b*x` has the least squared error.
"""
function linreg(x::AbstractVector, y::AbstractVector)
    # Least squares given
    # Y = a + b*X
    # where
    # b = cov(X, Y)/var(X)
    # a = mean(Y) - b*mean(X)
    if size(x) != size(y)
        throw(DimensionMismatch("x has size $(size(x)) and y has size $(size(y)), " *
            "but these must be the same size"))
    end
    mx = Statistics.mean(x)
    my = Statistics.mean(y)
    # don't need to worry about the scaling (n vs n - 1)
    # since they cancel in the ratio
    b = covm(x, mx, y, my)/varm(x, mx)
    a = my - b*mx
    return a, b
end

slope(x, y) = linreg(x, y)[2]

"""
    linear_region(x, y; dxi::Int = 1, tol = 0.2) -> ([ind1, ind2], slope)
Call [`linear_regions`](@ref), identify the largest linear region
and approximate the slope of the entire region using `linreg`.
Return the indices where
the region starts and stops (`x[ind1:ind2]`) as well as the approximated slope.
"""
function linear_region(x::AbstractVector, y::AbstractVector;
    dxi::Int = 1, tol::Real = 0.2)

    # Find biggest linear region:
    reg_ind = max_linear_region(linear_regions(x,y; dxi=dxi, tol=tol)...)
    # least squares fit:
    xfit = view(x, reg_ind[1]:reg_ind[2])
    yfit = view(y, reg_ind[1]:reg_ind[2])
    approx_tang = slope(xfit, yfit)
    return reg_ind, approx_tang
end

"""
    max_linear_region(lrs::Vector{Int}, tangents::Vector{Float64})
Find the biggest linear region and return it.
"""
function max_linear_region(lrs::Vector{Int}, tangents::Vector{Float64})
    dis = 0
    tagind = 0
    for i in 1:length(lrs)-1
        if lrs[i+1] - lrs[i] > dis
            dis = lrs[i+1] - lrs[i]
            tagind = i
        end
    end
    return [lrs[tagind], lrs[tagind+1]]
end

"""
    linear_regions(x, y; dxi::Int = 1, tol = 0.2) -> (lrs, tangents)
Identify regions where the curve `y(x)` is linear, by scanning the
`x`-axis every `dxi` indices (e.g. at `x[1] to x[5], x[5] to x[10], x[10] to x[15]`
and so on if `dxi=5`).

If the slope (calculated using `LsqFit`) of a region of width `dxi` is
approximatelly equal to that of the previous region,
within tolerance `tol`,
then these two regions belong to the same linear region.

Return the indices of `x` that correspond to linear regions, `lrs`,
and the approximated `tangents` at each region. `lrs` is a vector of `Int`.
Notice that `tangents` is _not_ accurate: it is not recomputed at every step,
but only when its error exceeds the tolerance `tol`! Use [`linear_region`](@ref)
to obtain a correct estimate for the slope of the largest linear region.
"""
function linear_regions(x::AbstractVector, y::AbstractVector;
    dxi::Int = 1, tol::Real = 0.2)

    maxit = length(x) ÷ dxi

    tangents = Float64[slope(view(x, 1:max(dxi, 2)), view(y, 1:max(dxi, 2)))]

    prevtang = tangents[1]
    lrs = Int[1] #start of first linear region is always 1
    lastk = 1

    # Start loop over all partitions of `x` into `dxi` intervals:
    for k in 1:maxit-1
        # tang = linreg(view(x, k*dxi:(k+1)*dxi), view(y, k*dxi:(k+1)*dxi))[2]
        tang = slope(view(x, k*dxi:(k+1)*dxi), view(y, k*dxi:(k+1)*dxi))
        if isapprox(tang, prevtang, rtol=tol)
            # Tanget is similar with initial previous one (based on tolerance)
            continue
        else
            # Tangent is not similar.
            # Push new tangent for a new linear region
            push!(tangents, tang)

            # Set the START of a new linear region
            # which is also the END of the previous linear region
            push!(lrs, k*dxi)
            lastk = k
        end

        # Set new previous tangent (only if it was not the same as current)
        prevtang = tang
    end
    push!(lrs, length(x))
    return lrs, tangents
end

#######################################################################################
# Dimensions
#######################################################################################
"""
    estimate_boxsizes(dataset::AbstractDataset; k::Int = 12, z = -1, w = 1)
Return `k` exponentially spaced values: `10 .^ range(lower+w, upper+z, length = k)`.

`lower` is the magnitude of the
minimum pair-wise distance between datapoints while `upper` is the magnitude
of the maximum difference between
greatest and smallest number among each timeseries.

"Magnitude" here stands for order of magnitude, i.e. `round(log10(x))`.
"""
function estimate_boxsizes(data::AbstractDataset{D, T};
    k::Int = 12, z = -1.0, w = 1.0) where {D, T<:Number}

    mi, ma = minmaxima(data)
    upper = round(log10(maximum(ma - mi)))

    mindist = min_pairwise_distance(data)[2]
    lower = ceil(log10(mindist)) # ceil necessary to not use smaller distance.

    if lower ≥ upper
        error(
        "Boxsize estimation failed: `upper` was found ≥ than "*
        "`lower`. Adjust keywords or provide a bigger dataset.")
    end
    if lower + w + 2 ≥ upper + z
        @warn "Boxsizes limits do not differ by 2 orders of magnitude or more. "*
        "Setting `w -= 0.5; z += 0.5`. Please adjust keywords or provide a bigger dataset."
        w -= 0.5; z += 0.5
    end

    return 10.0 .^ range(lower+w, stop = upper+z, length = k)
end

estimate_boxsizes(ts::AbstractMatrix; kwargs...) =
estimate_boxsizes(convert(Dataset, ts); kwargs...)

"""
    generalized_dim(α, dataset [, sizes]) -> D_α
Return the `α` order generalized dimension of the `dataset`, by calculating
the [`genentropy`](@ref) for each `ε ∈ sizes`.

## Description
The returned dimension is approximated by the
(inverse) power law exponent of the scaling of the [`genentropy`](@ref)
versus the box size `ε`, where `ε ∈ sizes`.

Calling this function performs a lot of automated steps:

  1. A vector of box sizes is decided by calling `sizes = estimate_boxsizes(dataset)`,
     if `sizes` is not given.
  2. For each element of `sizes` the appropriate entropy is
     calculated, through `d = genentropy.(α, sizes, dataset)`.
     Let `x = -log.(sizes)`.
  3. The curve `d(x)` is decomposed into linear regions,
     using [`linear_regions`](@ref)`(x, d)`.
  4. The biggest linear region is chosen, and a fit for the slope of that
     region is performed using the function [`linear_region`](@ref).
     This slope is the return value of `generalized_dim`.

By doing these steps one by one yourself, you can adjust the keyword arguments
given to each of these function calls, refining the accuracy of the result.

The following aliases are provided:

  * α = 0 : `boxcounting_dim`, `capacity_dim`
  * α = 1 : `information_dim`
"""
function generalized_dim(α, data::AbstractDataset, sizes = estimate_boxsizes(data))
    dd = genentropy.(α, sizes, Ref(data))
    return linear_region(-log.(sizes), dd)[2]
end
generalized_dim(α, matrix::AbstractMatrix, args...) =
generalized_dim(α, convert(AbstractDataset, matrix), args...)

# Aliases
"capacity_dim(args...) = generalized_dim(0, args...)"
capacity_dim(args...) = generalized_dim(0, args...)
boxcounting_dim = capacity_dim

"information_dim(args...) = generalized_dim(1, args...)"
information_dim(args...) = generalized_dim(1, args...)

"""
    kaplanyorke_dim(lyapunovs::AbstractVector)
Calculate the Kaplan-Yorke dimension, a.k.a. Lyapunov dimension [1].

## Description
The Kaplan-Yorke dimension is simply the point where
`cumsum(lyapunovs)` becomes zero (interpolated). If
the sum of the exponents never becomes negative the function
will return the length of the input vector.

Useful in combination with [`lyapunovs`](@ref).

## References

[1] :  J. Kaplan & J. Yorke,
*Chaotic behavior of multidimensional difference equations*,
Lecture Notes in Mathematics vol. **730**, Springer (1979)
"""
function kaplanyorke_dim(v::AbstractVector)
    issorted(v, rev = true) || throw(ArgumentError(
    "The lyapunov vector must be sorted from most positive to most negative"))

    s = cumsum(v); k = length(v)
    # Find k such that sum(λ_i for i in 1:k) is still possitive
    for i in eachindex(s)
        if s[i] < 0
            k = i-1
            break
        end
    end

    if k == 0
        return zero(v[1])
    elseif k < length(v)
        return k + s[k]/abs(v[k+1])
    else
        return typeof(v[1])(length(v))
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 2dd190ae-3888-11eb-2e68-090f885617d6
using DataStructures

# ╔═╡ 47c147e6-3888-11eb-1574-cd2d6aea7b1c
file = open(f->read(f, String), "day_8.txt")

# ╔═╡ 019df412-388a-11eb-17c1-2bfe6376975b
#first match get command, second get value of comand
regex = r"(\w{3}) \+?(-?\d+)"

# ╔═╡ f875f352-3952-11eb-3d4a-979e6ccb162f
to_int(x) = parse(Int64,x)

# ╔═╡ 913adf24-3955-11eb-2cf3-5bbfeee310db
program =[ (x.captures[1], to_int(x.captures[2])) for x in eachmatch(regex,file)]

# ╔═╡ 90ffacfe-388a-11eb-3224-377942614d33
begin
	mutable struct Register
		pointer::Int64
		value::Int64

	end

	Register() = Register(1,0)
	Register(pointer,value) = Register(pointer,value)
end

# ╔═╡ 2ffb9de0-3953-11eb-02f2-71323a50ca88
begin
	function nop(register, x)
		register.pointer += 1
	end

	function acc(register, x)
		register.value += x
		register.pointer += 1
	end

	function jmp(register, x)
		register.pointer += x
	end
end

# ╔═╡ 0edb82ea-3953-11eb-0547-5158b0364ddc
begin
	commands = Dict()
	commands["nop"] = nop;
	commands["acc"] = acc;
	commands["jmp"] = jmp;

	function execute_line(register)
		line = program[register.pointer]
		commands[line[1]](register, line[2])
	end
end

# ╔═╡ ab476d98-3954-11eb-117e-dbb761d47b74
begin

	function execute_program()
		register= Register()
		visited = Dict() #every line executed

		while(haskey(visited,register.pointer) == false)
			visited[register.pointer] = true
			execute_line(register)
		end

		return (register.value , register.pointer)
	end
end

# ╔═╡ 11487890-3954-11eb-14f9-139bb83e3c6c
execute_program()

# ╔═╡ aab19fbe-3956-11eb-0e2e-4d1994109154
md"# Question 2"

# ╔═╡ ee71662e-3958-11eb-0204-09256c308b5f
begin

	function execute_program2(register,visited,has_switched)

		while(haskey(visited,register.pointer) == false)

			#Got end of program
				if size(program)[1]+1 == register.pointer
					return register.value
			end

			visited[register.pointer] = true

			line = program[register.pointer]
			cmd = commands[line[1]]
			cmd_value = line[2]

			#swap jump or swap
			if(has_switched == false)

				if(cmd == jmp)
					new_pointer = register.pointer + 1
					new_register = Register(new_pointer,register.value)
					result =execute_program2(new_register,copy(visited),true)

					if(result isa Int)
						return result
					end

				elseif (cmd == nop)
					new_pointer = register.pointer + cmd_value
					new_register = Register(new_pointer,register.value)
					result =execute_program2(new_register,copy(visited),true)

					if(result isa Int)
						return result
					end
				end
			end

			cmd(register, cmd_value)
		end

		#program visited one line 2 times
		return false

	end
end

# ╔═╡ 5139d282-395b-11eb-15f9-9b3bd7c3eec3
execute_program2(Register(),Dict(),false)

# ╔═╡ e1945406-395b-11eb-0d7a-192a3e1f408d
size(program)[1]+1

# ╔═╡ Cell order:
# ╠═2dd190ae-3888-11eb-2e68-090f885617d6
# ╠═47c147e6-3888-11eb-1574-cd2d6aea7b1c
# ╠═019df412-388a-11eb-17c1-2bfe6376975b
# ╠═f875f352-3952-11eb-3d4a-979e6ccb162f
# ╠═913adf24-3955-11eb-2cf3-5bbfeee310db
# ╠═90ffacfe-388a-11eb-3224-377942614d33
# ╠═2ffb9de0-3953-11eb-02f2-71323a50ca88
# ╠═0edb82ea-3953-11eb-0547-5158b0364ddc
# ╠═ab476d98-3954-11eb-117e-dbb761d47b74
# ╠═11487890-3954-11eb-14f9-139bb83e3c6c
# ╠═aab19fbe-3956-11eb-0e2e-4d1994109154
# ╠═ee71662e-3958-11eb-0204-09256c308b5f
# ╠═5139d282-395b-11eb-15f9-9b3bd7c3eec3
# ╠═e1945406-395b-11eb-0d7a-192a3e1f408d


text:
using Ballistics
using Test

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


text:
# Algorithm 2.3
# Newton's method for finfing roots

function bisect_approx(f, fp, p0, tolerance, maxIter)
	for i ∈ 1:maxIter
		p = p0 - (f(p0)/fp(p0))
		if abs(p - p0) < tolerance
			return p
		end
		p0 = p
	end
	throw(ArgumentError("Method failed"))
end

println(bisect_approx(x->(cos(x) - x), x->(-sin(x) - 1), π/4, 0.1, 100))


text:
## Let's look at the Monte Hall Problem.
# The first thing we'll do is to build a function to simulate
# a Monty Hall experiment. That is, we will have three doors.
# Behind one is a car. Behind the other two are goats.
# You will pick one door (at random).
# -- Work for yourself -- try picking just one specific door,
# -- say door 1.
# Then Monty Hall will show you one door that has a goat.
# And this is a door that you didn't pick...
# Then we have two strategies, stay and switch. Switch will
# will go to the door that isn't revealed.

"""
`monty_hall`
============

Input
-----
- `switch` a boolean (true/false) that says if we switch or not.

Returns:
- return true if we won the car or return false otherwise
"""
function monty_hall(switch::Bool)
    cardoor = rand(1:3) # this picks a random number from the range 1:3
    yourdoor = rand(1:3) # this picks a random door

    alldoors = Set(1:3)  # this creates a set
    delete!(alldoors, cardoor) # monty can't show you this door
    delete!(alldoors, yourdoor) # monty can't show you this door
    montydoor = rand(alldoors) # monty can show you this door or doors
    delete!(alldoors, montydoor)
    if length(alldoors) == 0
        switchdoor = cardoor
    else
        @assert(length(alldoors) == 1)
        switchdoor = first(alldoors) # gets the first thing in the set.
    end
    if switch
        yourdoor = switchdoor
    end
    return yourdoor == cardoor
end
monty_hall(false)

###
""" Run multiple trials and show the overall average success. """
function monty_trials(n::Int, switch::Bool)
    ncars = 0
    for i=1:n
        ncars += monty_hall(switch)
    end
    return ncars/n
end
monty_trials(10000, false)

## What happens if we DO switch?
monty_trials(10000, true)


text:
import MKL # Github Actions uses Intel CPUs, where MKL is faster than OpenBLAS
import Documenter
import Literate
import WhitenedRBMs

ENV["JULIA_DEBUG"] = "Documenter,Literate,WhitenedRBMs"

#=
We place Literate.jl source .jl files and the generated .md files inside docs/src/literate.
=#
const literate_dir = joinpath(@__DIR__, "src/literate")

#=
Helper function to remove all "*.md" files from a directory.
=#
function clear_md_files(dir::String)
    if isdir(dir)
        for (root, dirs, files) in walkdir(dir), file in files
            if endswith(file, ".md")
                rm(joinpath(root, file))
            end
        end
    end
end

#=
Remove previously Literate.jl generated files. This removes all "*.md" files inside
`literate_dir`. This is a precaution: if we build docs locally and something fails,
and then change the name of a source file (".jl"), we will be left with a lingering
".md" file which will be included in the current docs build. The following line makes
sure this does not happen.
=#
clear_md_files(literate_dir)

#=
Run Literate.jl on the .jl source files within docs/src/literate (recursively).
For each .jl file, this creates a markdown .md file at the same location as
and with the same name as the corresponding .jl file, but with the extension
changed (.jl -> .md).
=#
if isdir(literate_dir)
    for (root, dirs, files) in walkdir(literate_dir), file in files
        if endswith(file, ".jl")
            Literate.markdown(joinpath(root, file), root)
        end
    end
else
    @warn "There was nothing at" literate_dir
end

#=
Build docs.
=#
Documenter.makedocs(
    modules = [WhitenedRBMs],
    sitename = "WhitenedRBMs.jl",
    pages = [
        "Home" => "index.md",
        "Examples" => [
            "MNIST identity" => "literate/MNIST_identity.md",
            "MNIST centered" => "literate/MNIST_center.md",
            "MNIST whitened" => "literate/MNIST_white.md",
        ],
        "Reference" => "reference.md"
    ],
    strict = true
)

#=
After the docs have been compiled, we can remove the *.md files generated by Literate.
=#
clear_md_files(literate_dir)

#=
Deploy docs to Github pages.
=#
Documenter.deploydocs(
    repo = "github.com/cossio/WhitenedRBMs.jl.git",
    devbranch = "master"
)