text:
function SymbolicWedderburn.decompose(
    k::MP.AbstractPolynomialLike,
    hom::SymbolicWedderburn.InducedActionHomomorphism,
)
    # correct only if features(hom) == monomials

    indcs = [hom[mono] for mono in MP.monomials(k)]
    coeffs = MP.coefficients(k)

    return indcs, coeffs
end

function SymbolicWedderburn.ExtensionHomomorphism(
    action::SymbolicWedderburn.Action,
    basis::MB.MonomialBasis,
)
    monos = collect(basis.monomials)
    mono_to_index = Dict(monos[i] => i for i in eachindex(monos))
    return SymbolicWedderburn.ExtensionHomomorphism(action, monos, mono_to_index)
end

struct VariablePermutation <: SymbolicWedderburn.ByPermutations end
_map_idx(f, v::AbstractVector) = map(f, eachindex(v))
_tuple_map_idx(f, ::Tuple{}, i) = tuple()
_tuple_map_idx(f, v::Tuple, i) = (f(i), _tuple_map_idx(f, Base.tail(v), i + 1)...)
_map_idx(f, v::Tuple) = _tuple_map_idx(f, v, 1)
function SymbolicWedderburn.action(::VariablePermutation, p, mono::MP.AbstractMonomial)
    v = MP.variables(mono)
    MP.substitute(MP.Eval(), mono, v => _map_idx(i -> v[i^p], v))
end
abstract type OnMonomials <: SymbolicWedderburn.ByLinearTransformation end

function SymbolicWedderburn.action(
    a::Union{VariablePermutation,OnMonomials},
    el,
    term::MP.AbstractTerm,
)
    return MP.coefficient(term) * SymbolicWedderburn.action(a, el, MP.monomial(term))
end
function SymbolicWedderburn.action(
    a::Union{VariablePermutation,OnMonomials},
    el,
    poly::MP.AbstractPolynomial,
)
    return sum([SymbolicWedderburn.action(a, el, term) for term in MP.terms(poly)])
end

# TODO Move it to MultivariateBases
function MP.polynomialtype(
    ::Type{<:MB.AbstractPolynomialVectorBasis{PT}},
    T::Type,
) where {PT}
    C = MP.coefficienttype(PT)
    U = MA.promote_operation(*, C, T)
    V = MA.promote_operation(+, U, U)
    return MP.polynomialtype(PT, V)
end

struct Ideal{C,GT,AT<:SymbolicWedderburn.Action} <: SumOfSquares.Certificate.AbstractIdealCertificate
    pattern::Pattern{GT,AT}
    certificate::C
end
SumOfSquares.Certificate.get(certificate::Ideal, attr::SumOfSquares.Certificate.Cone) = SumOfSquares.Certificate.get(certificate.certificate, attr)
function SumOfSquares.matrix_cone_type(::Type{<:Ideal{C}}) where {C}
    return SumOfSquares.matrix_cone_type(C)
end
SumOfSquares.Certificate.get(::Type{<:Ideal}, ::SumOfSquares.Certificate.GramBasisType) = Vector{Vector{MB.FixedPolynomialBasis}}
SumOfSquares.Certificate.zero_basis_type(::Type{<:Ideal}) = MB.MonomialBasis
SumOfSquares.Certificate.zero_basis(::Ideal) = MB.MonomialBasis
function SumOfSquares.Certificate.get(certificate::Ideal, attr::SumOfSquares.Certificate.ReducedPolynomial, poly, domain)
    return SumOfSquares.Certificate.get(certificate.certificate, attr, poly, domain)
end

function matrix_reps(cert, R, basis, ::Type{T}, form) where {T}
    polys = R * basis.monomials
    return map(SymbolicWedderburn.gens(cert.pattern.group)) do g
        S = Matrix{T}(undef, length(polys), length(polys))
        for i in eachindex(polys)
            p = polys[i]
            q = SymbolicWedderburn.action(cert.pattern.action, g, p)
            coefs = MP.coefficients(q, basis.monomials)
            S[:, i] = _linsolve(R, coefs, form)
        end
        return S
    end
end

function SumOfSquares.Certificate.get(cert::Ideal, attr::SumOfSquares.Certificate.GramBasis, poly)
    basis = SumOfSquares.Certificate.get(cert.certificate, attr, poly)
    T = SumOfSquares._complex(Float64, SumOfSquares.matrix_cone_type(typeof(cert)))
    summands = SymbolicWedderburn.symmetry_adapted_basis(T, cert.pattern.group, basis, cert.pattern.action)
    # We have a new basis `b = vcat(R * basis.monomials for R in summands)``.
    # SymbolicWedderburn guarantees that the invariant subspace spanned by the
    # polynomials of the vector `R * basis.monomials` is invariant under the
    # action of the group. That is, the matrix representation `ρ(g)` induced by
    # this basis is block-diagonal with one block for each summand.
    # That block is the matrix `S` computed below.
    # So an invariant solution `b'*Q*b` satisfies `Diagonal(S' for S in ...) * Q * Diagonal(S for S in ...) = Q`.
    # Or in equivalently: `Q * Diagonal(S for S in ...) = Diagonal(inv(S') for S in ...) * Q`.
    return map(summands) do summand
        R = SymbolicWedderburn.basis(summand)
        #@show typeof(R)
        m = SymbolicWedderburn.multiplicity(summand)
        N = size(R, 1)
        d = SymbolicWedderburn.degree(summand)
        display(R)
        @show R * basis.monomials
        @show m
        @show d
        #display.(matrix_reps(cert, Matrix(T(1) * LinearAlgebra.I, 4, 4), basis, T, _RowEchelonMatrix()))
        S = matrix_reps(cert, R, basis, T, _RowEchelonMatrix())
        decomose_semisimple = d > 1
        if decomose_semisimple
            # If it's not orthogonal, how can we conclude that we can still use the semisimple summands block-decomposition ?
            # In Example 1.7.2 of Sagan's book, he uses Corollary 1.6.6 which requires that `X` and `Y` are irreducible
            # (where `X` and `Y` are here the `S` corresponding to two different summands).
            # Here, given semisimple representations `X` and `Y`, they are not irreducible if `m > 1`.
            # Furthermore, as they are not orthogonal, we have something like `T * X = inv(Y') * T`, so how can we know that
            # `X` and `inv(Y')` are not equivalent (to exclude the case 1. of Corollary 1.6.6) ?
            if !all(is_orthogonal, S)
                R = orthogonalize(R)
                display(R)
                S = matrix_reps(cert, R, basis, T, _OrthogonalMatrix())
                display.(S)
                for i in 1:size(R, 1)
                    R[i, :] = LinearAlgebra.normalize(R[i, :])
                end
                display(R)
                S = matrix_reps(cert, R, basis, T, _OrthogonalMatrix())
                display.(S)
                if !all(is_orthogonal, S)
                    error("The matrix representation induced from the action on the polynomial basis is not orthogonal.")
                    # We would like to just throw this warning and just not decompose the semisimple summand but
                    # as explained in the comment above, it's not even clear that the diagonalization induced by the simple summands is correct.
                    #@warn("The matrix representation induced from the action on the polynomial basis is not orthogonal. The $(m * d)-dimensional semisimple summand can be decomposed onto $m simple summands of degree $d so that the $(m * d) x $(m * d) diagonal block is reduced to $d identical copied of a single $m x $m diagonal block. However, as the action is not orthogonal, this decomposition will not happen.")
                    #decomose_semisimple = false
                end
            end
        end
        F = convert(Matrix{T}, R)
        if decomose_semisimple
            if m > 1
                U = ordered_block_diag(S, d)
            else
                U = Matrix{T}(LinearAlgebra.I, N, N)
            end
            if U === nothing
                error("Could not simultaneously block-diagonalize into $m identical $dx$d blocks")
            end
            # From Example 1.7.3 of
            # Sagan, The symmetric group, Springer Science & Business Media, 2001
            # we know that there exists `C` such that `Q = kron(C, I)` if we use
            # `(U[1:d] * F)' * basis.monomials`, `(U[d+1:2d] * F)' * basis.monomials`, ...
            # where `C` are some complex numbers as they are eigenvalues (see Corollary 1.6.8).
            # As `Q` is symmetric, we know the eigenvalues are real so we can take `C` real as well.
            # Moreover, `Q = kron(C, I)` is not block diagonal but we can get a block-diagonal
            # `Q = kron(I, Q)` by permuting the rows and columns:
            # `(U[1:d:(1+d*(m-1))] * F)' * basis.monomials`, `(U[2:d:(2+d*(m-1))] * F)' * basis.monomials`, ...
            map(1:d) do i
                MB.FixedPolynomialBasis(
                    (transpose(U[:, i:d:(i + d * (m - 1))]) * F) * basis.monomials,
                )
            end
        else
            F = convert(Matrix{T}, R)
            [MB.FixedPolynomialBasis(F * basis.monomials)]
        end
    end
end


text:
using .RasterDataSources

using .RasterDataSources: RasterDataSource

const RDS = RasterDataSources

"""
    Raster(T::Type{<:RasterDataSource}, [layer]; kw...) => Raster

Load a `RasterDataSource` as an `Raster`. `T` and `layers` are are passed to
`RasterDataSources.getraster`, while `kw` args are for both `getraster` and `Raster`.

# Keywords

- `month`: an `Int` between `1` and `12`, usually for `Climate` datasets
- `date`: a `DateTime` object, usually for `Weather` datasets.
- `res`: a `String` resolution, for datasets with multiple resolutions.

Other `Raster` keywords are passed to the `Raster` constructor.

See the docs for
[`RasterDatasources.getraster`](http://docs.ecojulia.org/RasterDataSources.jl/stable/#getraster)
for more specific details about data sources, layers and keyword arguments.
"""
function Raster(T::Type{<:RasterDataSource}, layer; crs=_source_crs(T), kw...)
    rds_kw, gd_kw = _filterkw(kw)
    filename = getraster(T, layer; rds_kw...)
    Raster(filename; name=RDS.layerkeys(T, layer), crs, gd_kw...)
end

"""
    RasterStack(T::Type{<:RasterDataSource}, [layers::Union{Symbol,AbstractArray,Tuple}]; kw...) => RasterStack

Load a `RasterDataSource` as an `RasterStack`. `T` and `layers` are passed to
`RasterDataSources.getraster`, while `kw` args are for both `getraster` and `RasterStack`.

# Keywords

- `month`: an `Int` between `1` and `12`, usually for `Climate` datasets.
- `date`: a `DateTime` object, usually for `Weather` datasets.
- `res`: a `String` resolution, for datasets with multiple resolutions.

Other `RasterStack` keywords are passed to the `RasterStack` constructor.

See the docs for
[`RasterDatasources.getraster`](http://docs.ecojulia.org/RasterDataSources.jl/stable/#getraster)
for more specific details about data sources, layers and keyword arguments.
"""
RasterStack(T::Type{<:RasterDataSource}; kw...) = RasterStack(T, RDS.layers(T); kw...)
RasterStack(T::Type{<:RasterDataSource}, layer::Symbol; kw...) = RasterStack(T, (layer,); kw...)
function RasterStack(T::Type{<:RasterDataSource}, layers::Tuple; crs=_source_crs(T), kw...)
    rds_kw, gd_kw = _filterkw(kw)
    filenames = map(l -> getraster(T, l; rds_kw...), layers)
    RasterStack(filenames; keys=RDS.layerkeys(T, layers), crs, gd_kw...)
end

"""
    RasterSeries(T::Type{<:RasterDataSource}, [layers::Union{Symbol,AbstractArray,Tuple}]; kw...) => AbstractRasterSeries

Load a `RasterDataSource` as an `AbstractRasterSeries`. `T`, `args` are are passed to
`RasterDataSource.getraster`, while `kw` args are for both `getraster` and `RasterSeries`.

# Keywords

- `month`: a `Vector` or range of `Int` between `1` and `12`, usually for `Climate` datasets.
- `date`: a `Vector` of `DateTime` objects, usually for `Weather` datasets.
- `res`: a `String` resolution, for datasets with multiple resolutions.

Other `RasterSeries` keywords are passed to the `RasterSeries` constructor.

See the docs for
[`RasterDatasources.getraster`](http://docs.ecojulia.org/RasterDataSources.jl/stable/#getraster)
for more specific details about data sources, layers and keyword arguments.
"""
RasterSeries(T::Type{<:RasterDataSource}; kw...) = RasterSeries(T, RDS.layers(T); kw...)
# DateTime time-series
function RasterSeries(T::Type{<:RasterDataSource}, layers;
    resize=_mayberesize(T), crs=_source_crs(T), mappedcrs=nothing, kw...
)
    monthdim = if haskey(values(kw), :month) values(kw)[:month] isa AbstractArray
        Dim{:month}(values(kw)[:month]; lookup=Sampled(; sampling=Intervals(Start())))
    else
        nothing
    end
    datedim = if haskey(values(kw), :date)
        dates = if values(kw)[:date] isa Tuple
            dates = RasterDataSources.date_sequence(T, values(kw)[:date])
            Ti(dates; lookup=Sampled(; sampling=Intervals(Start())))
        elseif values(kw)[:date] isa AbstractArray
            Ti(dates; lookup=Sampled(; sampling=Intervals(Start())))
        else
            nothing
        end
    else
        nothing
    end
    if isnothing(monthdim) && isnothing(datedim)
        throw(ArgumentError("A RasterSeries can only be constructed from a data source with `date` or `month` keywords that are AbstractArray. For other sources, use RasterStack or Raster directly"))
    end

    filenames = getraster(T, layers; kw...)
    can_duplicate = RDS.has_constant_dims(T) && RDS.has_constant_metadata(T)

    if filenames isa AbstractVector{<:AbstractVector}
        series = [RasterSeries(inner_fns, monthdim; resize, crs, mappedcrs, duplicate_first=can_duplicate) for inner_fns in filenames]
        return RasterSeries(series, datedim)
    elseif filenames isa AbstractVector
        dim = isnothing(datedim) ? monthdim : datedim
        return RasterSeries(filenames, dim; resize, crs, mappedcrs, duplicate_first=can_duplicate)
    end
end

_mayberesize(T) = RDS.has_matching_layer_size(T) ? nothing : crop

_source_crs(T) = nothing
_source_crs(T::Type{AWAP}) = crs=EPSG(4326)
_source_crs(T::Type{ALWB}) = crs=EPSG(4326)

function _filterkw(kw)
    rds = []; gd = []
    for p in kw
        dest = first(p) in (:date, :month, :res) ? rds : gd
        push!(dest, p)
    end
    rds, gd
end

ext: jl


text:
# Sets up to test random_neutral_walk() which is in src/Evolve.jl.

#include("../test/miller_examples.jl")

# tries to find a small numactive 3-bit multiplier circuit starting with
#   saved 3-bit multiplier circuits
g3m = goal_for_three_bit_multiplier()
funcs = [AND, XOR ]
c3m1 = build_chromosome( # First saved 3-bit multiplier circuit
       (1,2,3,4,5,6),
       ((XOR,Integer[3, 6]),(AND,Integer[6, 3]),(AND,Integer[5, 2]),(AND,Integer[3, 9]),(AND,Integer[4, 7]),
        (AND,Integer[1, 5]),(AND,Integer[4, 2]),(AND,Integer[1, 4]),(AND,Integer[5, 3]),(AND,Integer[8, 14]),
        (AND,Integer[6, 2]),(XOR,Integer[12, 16]),(XOR,Integer[4, 1]),(XOR,Integer[19, 10]),(XOR,Integer[11, 9]),
        (AND,Integer[6, 20]),(AND,Integer[11, 16]),(AND,Integer[22, 21]),(AND,Integer[13, 18]),(XOR,Integer[24, 13]),
        (AND,Integer[15, 16]),(XOR,Integer[18, 13]),(AND,Integer[10, 26]),(AND,Integer[17, 12]),
        (XOR,Integer[29, 30]),(AND,Integer[8, 8]),(AND,Integer[28, 31]),(XOR,Integer[14, 33]),(XOR,Integer[12, 16]),
        (XOR,Integer[25, 27]),(XOR,Integer[31, 28]),(XOR,Integer[15, 17]),(AND,Integer[32, 32]),
        (XOR,Integer[21, 22]),(XOR,Integer[34, 36])),
       (36,37,38,39,40,41))
c3m2 = build_chromosome( # Second saved 3-bit multiplier circuit
         (1,2,3,4,5,6),
         ((AND,Integer[3, 3]),(AND,Integer[7, 5]),(AND,Integer[2, 4]),(AND,Integer[1, 4]),
          (AND,Integer[8, 9]),(AND,Integer[5, 2]),(AND,Integer[1, 5]),(AND,Integer[2, 6]),
          (XOR,Integer[12, 11]),(AND,Integer[9, 13]),(AND,Integer[14, 8]),(AND,Integer[7, 5]),
          (AND,Integer[6, 1]),(AND,Integer[19, 12]),(XOR,Integer[4, 15]),(XOR,Integer[9, 13]),
          (XOR,Integer[15, 10]),(XOR,Integer[17, 4]),(XOR,Integer[11, 22]),(AND,Integer[6, 7]),
          (XOR,Integer[20, 23]),(XOR,Integer[27, 12]),(XOR,Integer[12, 19]),(AND,Integer[26, 27]),
          (AND,Integer[7, 24]),(AND,Integer[25, 30]),(XOR,Integer[30, 20]),(AND,Integer[16, 12]),
          (XOR,Integer[32, 21]),(XOR,Integer[14, 18]),(XOR,Integer[34, 32]),(AND,Integer[26, 26]),
          (XOR,Integer[25, 33]),(AND,Integer[28, 35]),(XOR,Integer[29, 31])),
         (36,37,38,39,40,41))
rnw1 = random_neutral_walk( c3m1, [g3m], funcs, 20000, num_mutations=2 )
rnw2 = random_neutral_walk( c3m1, [g3m], funcs, 20000, num_mutations=2, reduce_numactive_reps=500 )


text:
using CSV, GLM, Plots; pyplot()

data = CSV.read("../data/polynomialData.csv")
data.X2 = abs2.(data.X)
model1 = lm(@formula(Y ~ X + X2), data)
model2 = lm(@formula(Y ~ X*X), data)
println(model1)
println(model2)

yHat(x) = coef(model1)[1] + coef(model1)[2]*x + coef(model1)[3]*x^2

xGrid = -5:0.1:5
scatter(data.X, data.Y, c=:blue, msw=0)
plot!(xGrid, yHat.(xGrid),
    xlabel="X", ylabel="Y", c=:red, legend=:none)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
import Base: map, map!
import Core.Compiler: PhiNode, PiNode, ssamap, userefs
import MacroTools: walk

walk(x::PhiNode, inner, outer) =
  outer(PhiNode(x.edges,
                [inner(isassigned(x.values, i) ? x.values[i] : undef)
                 for i in 1:length(x.values)]))

walk(x::PiNode, inner, outer) = outer(PiNode(inner(x.val), x.typ))

xcall(mod::Module, f::Symbol, args...) = Expr(:call, GlobalRef(mod, f), args...)
xcall(f::Symbol, args...) = xcall(Base, f, args...)
xcall(f, args...) = Expr(:call, f, args...)

map(f, br::Branch) = Branch(br, condition = f(br.condition), args = f.(br.args))

function map(f, b::BasicBlock)
  stmts = map(x -> Statement(x, expr = f(x.expr)), b.stmts)
  branches = map(br -> map(f, br), b.branches)
  BasicBlock(stmts, b.args, b.argtypes, branches)
end

function map!(f, b::BasicBlock)
  map!(x -> Statement(x, expr = f(x.expr)), b.stmts, b.stmts)
  map!(br -> Branch(br, condition = f(br.condition), args = f.(br.args)), b.branches, b.branches)
  return b
end

function map!(f, b::Block)
  map!(f, basicblock(b))
end

function map(f, ir::IR)
  IR(ir.defs, map.(f, ir.blocks), ir.lines, ir.meta)
end

function map!(f, ir::IR)
  for b in blocks(ir)
    map!(f, b)
  end
  return ir
end

walk(st::Statement, inner, outer) = Statement(st, expr = inner(st.expr))
walk(bb::BasicBlock, inner, outer) = map(inner, bb)
walk(bb::Branch, inner, outer) = map(inner, bb)
walk(b::Block, inner, outer) = walk(basicblock(b), inner, outer)

walk(ir::IR, inner, outer) = outer(map(inner, ir))

prewalk!(f, ir::Union{IR,Block}) = map!(x -> prewalk(f, x), ir)
postwalk!(f, ir::Union{IR,Block}) = map!(x -> postwalk(f, x), ir)

varmap(f, x) = prewalk(x -> x isa Variable ? f(x) : x, x)

exprtype(x::GlobalRef) = isconst(x.mod, x.name) ? Typeof(getfield(x.mod, x.name)) : Any

exprtype(ir::IR, x::GlobalRef) = exprtype(x)
exprtype(ir::IR, x::QuoteNode) = Typeof(x.value)
exprtype(ir::IR, x::Expr) = error(x)
exprtype(ir::IR, x) = Typeof(x)

function exprtype(ir::IR, x::Variable)
  b, i = get(ir.defs, x.id, (-1, -1))
  b == -1 && error("No such variable $x")
  if i > 0
    widenconst(ir[x].type)
  else
    widenconst(ir.blocks[b].argtypes[-i])
  end
end

function exprline(ir::IR, x::Variable)
  b, i = get(ir.defs, x.id, (-1, -1))
  i > 0 || return
  get(ir.lines, ir[x].line, nothing)
end


text:
"""
    isleaf(node::Node)
Test whether a node is a leaf or not.
"""
isleaf(node::Node) = node.children === nothing

"""
    isroot(node::Node)
Return `true` if `node` is the root node (meaning, it has no parent).
"""
isroot(node::Node) = node.parent === nothing

"""
    children(node::Node)

Return the immediate children of `node`.
"""
function children(node::Node)
    isleaf(node) ? Vector{Node}() : collect(values(node.children))
end

"""
    ordered_children(node)

Return the children as an array, ordered first by "+"
"""
function ordered_children(node)
    links_chnodes = Array{Node,1}()
    for (name, chnode) in node.children
        if chnode.MTG.link == "+"
            pushfirst!(links_chnodes, chnode)
        else
            push!(links_chnodes, chnode)
        end
    end
    return links_chnodes
end

"""
Add a new child to a parent node, and add the parent node as the parent.

See also [`insert_child!`](@ref) for a more user-friendly approach.
"""
function addchild!(parent::Node, id::Int, MTG::M, attributes) where {M<:AbstractNodeMTG}
    child = Node(id, parent, MTG, attributes)
    addchild!(parent, child)
end

function addchild!(parent::Node, id::Int, MTG::M) where {M<:AbstractNodeMTG}
    child = Node(id, parent, MTG)
    addchild!(parent, child)
end

function addchild!(parent::Node, child::Node; force = false)

    if child.parent === missing || force == true
        child.parent = parent
    elseif child.parent != parent && force == false
        error("The node already has a parent. Hint: use `force=true` if needed.")
    end

    if parent.children === nothing
        parent.children = Dict{Int,Node}(child.id => child)
    else
        push!(parent.children, child.id => child)
    end
end

"""
Find the root node of a tree, given any node in the tree.
"""
function get_root(node::Node)
    if isroot(node)
        return (node)
    else
        get_root(node.parent)
    end
end

"""
    siblings(node::Node)

Return the siblings of `node` as a vector of nodes (or `nothing` if non-existant).
"""
function siblings(node::Node)
    # If there is no parent, no siblings, return nothing:
    node.parent !== nothing || return nothing

    all_siblings = children(node.parent)

    all_siblings[findall(x -> x != node, all_siblings)]
end

function nextsibling(node::Node)
    # If there is no parent, no siblings, return nothing:
    node.parent !== nothing || return nothing
    # If the siblings field is not empty, return its value:
    node.siblings !== nothing && return node.siblings
    # Else, compute the siblings:
    all_siblings = children(node.parent)
    # Get the index of the current node in the siblings:
    node_index = findfirst(x -> x == node, all_siblings)
    if node_index < length(all_siblings)
        all_siblings[node_index+1]
    else
        nothing
    end
end

"""
    max_id(mtg)

Returns the maximum id of the mtg
"""
function max_id(mtg)
    maxid = [0]

    function update_maxname(id, maxid)
        id > maxid[1] ? maxid[1] = id : nothing
    end

    traverse!(get_root(mtg), x -> update_maxname(x.id, maxid))

    return maxid[1]
end

"""
    new_id(mtg)
    new_id(mtg, max_id)

Make a new unique identifier by incrementing on the maximum node id.
Hint: prefer using `max_id = max_id(mtg)` and then `new_id(mtg, max_is)` for performance
if you do it repeatidely.
"""
function new_id(max_id::Int)
    max_id + 1
end

function new_id(mtg::Node)
    new_id(max_id(mtg))
end

"""
    get_attributes(node)

Get all attributes names available on the node and its children.
"""
get_attributes(node) = unique(vcat(traverse(node, node -> collect(keys(node.attributes)))...))

"""
    names(node)

Get all attributes names available on the node and its children. This is an alias for
[`get_attributes`](@ref).
"""
names(node::Node) = get_attributes(node)

"""
    list_nodes(node)

List all nodes IDs in the subtree of `node`.
"""
list_nodes(node) = traverse(node, node -> node.id)


text:
# Utility
@inline function _map(f::F, a::Tuple{Vararg{Any, N}}) where {F, N}
    ntuple(Val(N)) do i
        f(a[i])
    end
end

@inline function _map(f::F, a::Tuple{Vararg{Any, N}}, b::Tuple{Vararg{Any, N}}) where {F, N}
    ntuple(Val(N)) do i
        f(a[i], b[i])
    end
end


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

"""
    pairwise(γ, X)

Evaluate variogram `γ` between all n² pairs of columns in a
m-by-n matrix `X` efficiently.
"""
function pairwise(γ::Variogram, X::AbstractMatrix)
  m, n = size(X)
  R = result_type(γ, X, X)
  Γ = Array{R}(undef, n, n)
  pairwise!(Γ, γ, X)

  Γ
end

"""
    pairwise!(Γ, γ, X)

Non-allocating pairwise evaluation.
"""
function pairwise!(Γ, γ::Variogram, X::AbstractMatrix)
  m, n = size(X)
  @inbounds for j=1:n
    xj = view(X, :, j)
    for i=j+1:n
      xi = view(X, :, i)
      Γ[i,j] = γ(xi, xj)
    end
    Γ[j,j] = γ(xj, xj)
    for i=1:j-1
      Γ[i,j] = Γ[j,i] # leverage the symmetry
    end
  end
end

"""
    pairwise(γ, domain, locations)

Evaluate variogram `γ` between all `locations` in `domain`.
"""
function pairwise(γ::Variogram, domain::AbstractDomain{T,N},
                  locations::AbstractVector{Int}) where {T<:Real,N}
  xi = MVector{N,T}(undef)
  xj = MVector{N,T}(undef)
  n = length(locations)
  R = result_type(γ, xi, xj)
  Γ = Matrix{R}(undef, n, n)
  @inbounds for j=1:n
    coordinates!(xj, domain, locations[j])
    for i=j+1:n
      coordinates!(xi, domain, locations[i])
      Γ[i,j] = γ(xi, xj)
    end
    Γ[j,j] = γ(xj, xj)
    for i=1:j-1
      Γ[i,j] = Γ[j,i] # leverage the symmetry
    end
  end

  Γ
end

"""
    pairwise(γ, domain, locations₁, locations₂)

Evaluate variogram `γ` between `locations₁` and `locations₂` in `domain`.
"""
function pairwise(γ::Variogram, domain::AbstractDomain{T,N},
                  locations₁::AbstractVector{Int},
                  locations₂::AbstractVector{Int}) where {T<:Real,N}
  xi = MVector{N,T}(undef)
  xj = MVector{N,T}(undef)
  m = length(locations₁)
  n = length(locations₂)
  R = result_type(γ, xi, xj)
  Γ = Array{R}(undef, m, n)
  @inbounds for j=1:n
    coordinates!(xj, domain, locations₂[j])
    for i=1:m
      coordinates!(xi, domain, locations₁[i])
      Γ[i,j] = γ(xi, xj)
    end
  end

  Γ
end


text:
module HamiltonianABC

using Distributions

export
    logdensity,
    simulate!,
    propose,
    mcmc,
    RWMH

"""
    logdensity(model, θ)

Return the log density for `model` at parameter `θ`.

A *model* is an object that can behave as a probability distribution over some
space, possibly by simulation (see [`simulate`](@ref)).
"""
function logdensity end

"""
    simulate!(model)

Re-generate the simulated shocks for `model`. For Monte Carlo calculation of
integrals, `model` should save appropriate shocks, independent of
parameters. Evaluations of [`logdensity`](@ref) then transform these using the
parameters.
"""
function simulate! end

"""
Random walk Metropolis-Hastings.

In this library, it is provided purely for testing and pedagogical purposes. See
[`NUTS`](@ref) for a more practical alternative.
"""
struct RWMH{S <: AbstractMatrix}
    "Covariance matrix for proposals."
    Σ::S
end

"""
    propose(transition, model, θ)

Generate a Metropolis-Hastings proposal using the given `transition` algorithm
and `model`, starting from parameters `θ`.

Return `θ′, a` where `θ′` is the proposed new parameter and `a` is the
acceptance probability.

The method `logdensity(model, θ)` needs to be supported by the arguments.
"""
function propose(transition::RWMH, model, θ)
    θ′ = rand(MvNormal(θ, transition.Σ))
    a = min(1, exp(logdensity(model, θ′) - logdensity(model, θ)))
    θ′, a
end

"""
    mcmc(transition, model, θ, N)

Runs a Markov Chain Monte Carlo algorithm by calling [`propose`](@ref)
repeatedly `N` times. After each iteration, re-simulate the shocks for model
using `simulate!`.

Return the vector of parameters, and the average acceptance rate.

The following interface needs to be supported by the arguments: -
`logposterior(model, θ)`, `simulate!(model)`.
"""
function mcmc{T}(transition, model, θ::T, N)
    chain = Vector{T}(N)
    as = Vector{Float64}(N)
    for i in 1:N
        θ′, a = propose(transition, model, θ)
        if rand() ≤ a
            θ = θ′
        end
        as[i] = a
        chain[i] = θ
        simulate!(model)
    end
    chain, mean(as)
end

end # module


text:
"""
    get_format_filenames(format, folder)

Returns an array of filenames in folder (string) with suffix format (string)
"""
get_format_filenames(format::String, folder::String, prefix::String) =[string(folder,f) for f in  readdir(folder) if (endswith(f, format) && startswith(f,prefix))]

"""
    get_spatiotemporal_var_atts()

Return the variable attributes for the spatiotemporal variable (x,y,y)
"""
function get_spatiotemporal_var_atts()
    x_atts = Dict("longname" => "x co-ordinates of ice grid points (h grid)",  "units" => "m")
    y_atts = Dict("longname" => "y co-ordinates of ice grid points (h grid)",  "units" => "m")
    time_atts = Dict("longname" => "Time", "units" => "years");
    return x_atts, y_atts, time_atts
end

"""
    get_spatial_dimensions()

Return one-dimensional arrays of the spatial variables
"""
function get_spatial_dimensions(fname)
    format = return_extension(fname)
    if format == "mat"
        vars = matread(fname)
    elseif format == "jld2"
        vars = load(fname)
    end
    X = vars["x"]
    X =  X[:,1]
    Y = vars["y"]
    Y = Y[1,:]
    return X, Y
end

"""
    get_time_output(filenames)

Return the times associated with filenames
"""

function get_times(filenames)
    t = zeros(length(filenames))
    for i = 1:length(filenames)
        format = return_extension(filenames[i])
        vars = get_output_as_dict(filenames[i],format)
        t[i] = vars["t"]
    end
    return t
end
"""
    get_output_as_dict(filename,format)

Read the data in filename according to different format specification
"""
function get_output_as_dict(filename,format)
    if format == "mat"
        vars= matread(filename)
    elseif format == "jld2"
        vars = load(filename)
    end
return vars
end

"""
    return_extension(file)

Return the extension of the input file
"""
return_extension(file) =  file[(findlast(isequal('.'),file)+1):end];

"""
make_ncfile(folder,format)

Wrapper script to zip the output files in "folder" with type "format" to an nc file with name nc_name_full (including path)
"""
function make_ncfile(format, folder, nc_name, prefix)
    #check that the input format is
    filenames = get_format_filenames(format, folder, prefix)
    if ~isempty(filenames)
        make_ncfile_from_filenames(filenames, format, nc_name)
    else
        println("attempted to zip the outputs to nc format, but did not find any files...")
    end
    return nothing
end

"""
    make_ncfile_from_filenames(filenames, format)

Output an nc file from filenames, which have format "format" to a file with name nc_name.
nc_name must contain the path as well!
"""
function make_ncfile_from_filenames(filenames, format, nc_name_full)
    #get the spatial and temporal variables as arrays of size N x 1
    x, y = get_spatial_dimensions(filenames[1])
    t    = get_times(filenames)

    #setup attributes for spatiotemporal variables
    x_atts, y_atts, time_atts = get_spatiotemporal_var_atts()

    #setup Dimensions
    t_ncdim = NcDim("TIME", t, atts = time_atts)
    x_ncdim = NcDim("x", x, atts = x_atts)
    y_ncdim = NcDim("y", y, atts = y_atts)

    #create the nc file variables
    output_dict = Dict() #dictionary for the output values
    output_vars = Dict() #dictionary for the output variables

    #get the keys
    if format == "mat"
        filekeys = keys(matread(filenames[1]))
    elseif format == "jld2"
        filekeys = keys(load(filenames[1]))
        #println(typeof(filekeys))
    end

    for key in filekeys
        if ~(key in ["x", "y", "t"]) #if this isn't a spatial dimensions
            var_out = zeros(length(x), length(y), length(t))

            #check the size of this variable in the first file
            sz = size(get_output_as_dict(filenames[1],format)[key])
            if sz == (length(x), length(y))
                #if the size is OK, create a variable and add to the array
                ncvar_key = NcVar(key, Array([x_ncdim, y_ncdim, t_ncdim]))
                output_vars[key] = ncvar_key

                #get the data in an array
                for i = 1:length(filenames)
                    var_out[:,:,i] = get_output_as_dict(filenames[i],format)[key]
                end

                #add this to dictionary
                output_dict[key] = var_out
            else
                @warn string("found an output variable (", key, ") who's spatial dimensions do not match the co-ordinates. Skipping this variable from the nc output...")
            end
        end
    end

    #make the nc file
    isfile(nc_name_full) && rm(nc_name_full)
    NetCDF.create(nc_name_full, [val for val in values(output_vars)]) do nc
        # Writing data to the file is done using putvar
        for key in keys(output_dict) #for every variable of the correct size
            NetCDF.putvar(nc,key,output_dict[key])
        end
    end
    return nothing
end

"""
    zip_output(simulation)

Zip all of the output files from simulation.
"""
function zip_output(simulation)
    @unpack output_params = simulation
    if output_params.zip_format == "nc"
        nc_name_full = string(output_params.output_path, output_params.prefix, ".nc")
        make_ncfile(output_params.output_format, output_params.output_path, nc_name_full, output_params.prefix)
    end
    return nothing
end


text:
module breeder
importall common, selection, representation, mutation, variation, crossover, species
using core, utility, individual
export Breeder, BreederDefinition

"""
The base type used by all breeders.
"""
abstract Breeder

"""
The base type used by all breeder definitions.
"""
abstract BreederDefinition

# Include each of the breeders.
include("breeder/flat.jl")
end


text:
using PoissonBox
using Test
using IterativeSolvers
import PoissonBox as PB
using Random: MersenneTwister

@testset "laplace eqn operator" begin
    forward = PoissonBox.forward
    @testset "1d" begin
        @test forward(1:5) ≈ [1,0,0,0,5]
        @test forward([0,0,1,0,0]) ≈ [0,0.5,-1,0.5,0]
    end
    @testset "2d" begin
        @test forward(
            [0 0 0 0 0
                 0 0 0 0 0
                 0 0 1 0 0
                 0 0 0 0 0
                 0 0 0 0 0
                ]) ≈ [
                 0 0    0    0    0
                 0 0    0.25 0    0
                 0 0.25 -1   0.25 0
                 0 0    0.25 0    0
                 0 0    0    0    0
                ]
    end
    @testset "3d" begin
        arr = fill(0, (5,5,5))
        arr[3,3,3] = 1
        expected = fill(0.0,(5,5,5))
        expected[3,3,3] = -1
        expected[2,3,3] = 1/6
        expected[4,3,3] = 1/6
        expected[3,2,3] = 1/6
        expected[3,4,3] = 1/6
        expected[3,3,2] = 1/6
        expected[3,3,4] = 1/6
        @test forward(arr) ≈ expected
    end
    @testset "3d bdry" begin
        arr = randn(3,4,5)
        out = forward(arr)
        expected = copy(arr)
        expected[2:2, 2:3, 2:4] .= 0
        out[2:2, 2:3, 2:4] .= 0
        @test expected ≈ out
    end
end

@testset "solve" begin
    solve = PoissonBox.solve
    y = Float64[1, 0, 0, 0, 5]
    @test solve(y) ≈ 1:5
    for shape in [(5,),
                  (4,3),
                  (5,4),
                  (4,3,4),
                 ]
        y = rand(MersenneTwister(0), shape...)
        x = solve(y)
        @test PB.forward(x) ≈ y rtol=1e-3
    end
end

@testset "solve low level" begin
    op, v = PB.poisson_problem([1, 0, 0, 0, 5])
    sol = bicgstabl(op, v, 2)
    u = PB.reshape_solution(op, sol)
    @test u ≈ 1:5
end


text:
logbesseli(ν, x) = log(besselix(ν * one(real(x)), x)) + abs(real(x))

# hypergeometric functions with matrix arguments
function pFq end

pFq(::Tuple{}, ::Tuple{}, x) = exp(tr(x)) # ₀F₀
# TODO: Implement ₀F₁
# see P Koev, A Edelman. The efficient evaluation of the hypergeometric function of a
# matrix argument. Math. Comp. 75 (2006), 833-846
# pFq(a::NTuple{p}, b::NTuple{q}, x)

logpFq(a, b, z) = log(pFq(a, b, z))


text:
using MLStyle
using PrettyPrinting

function jsonlit(ex)
    (∉)  = !isa
    function mk_kv(kv)
        (key, value) = @match kv begin
            :($k : $v) => (string(k), jsonlit(v))
            k          => (string(k), k)
        end
        :($key => $value)
    end

    @match ex begin
        :{ $(kvs...)} =>
            begin
                kvs = (mk_kv(kv) for kv in kvs if kv ∉ LineNumberNode)
                :(Dict($(kvs...)))
            end
        :[$(elts...)] => :[$(map(jsonlit, elts)...)]
        a => a
    end
end

macro jsonlit(ex)
    jsonlit(ex) |> esc
end

key = 1
value = 2
pprint(@jsonlit {
    key   : value,
    "key" :  "value",
    key   : value,
    (key+value),
    value : {
        key,
        value,
        key: {
            value
        }
    }
})

Dict(
    "key" => 2,
    "key + value" => 3,
    "value" => Dict(
                "key" => Dict("value" => 2),
                "value" => 2
               )
)

max_stars_repo_name: thautwarm/Slides

max_issues_repo_name: thautwarm/Slides

max_forks_repo_name: thautwarm/Slides


text:
using Pkg

pkg"add https://github.com/kirstvh/BioCCP"
pkg"precompile"

max_issues_repo_name: kirstvh/BioCCP

max_forks_repo_name: kirstvh/BioCCP


text:
"""
    TestIterator{T,IS<:Base.IteratorSize,IE<:Base.IteratorEltype}

A configurable iterator for testing purposes.

    TestIterator(data, itersize, itereltype)
    TestIterator(data)

The iterator wraps another iterator `data`, such as an array, that must have at least as
many features implemented as the test iterator and have a `FiniteDifferences.to_vec`
overload. By default, the iterator it has the same features as `data`.

The optional methods `eltype`, `length`, and `size` are automatically defined and forwarded
to `data` if the type arguments indicate that they should be defined.
"""
struct TestIterator{T,IS,IE}
    data::T
end
function TestIterator(data, itersize::Base.IteratorSize, itereltype::Base.IteratorEltype)
    return TestIterator{typeof(data),typeof(itersize),typeof(itereltype)}(data)
end
TestIterator(data) = TestIterator(data, Base.IteratorSize(data), Base.IteratorEltype(data))

Base.iterate(iter::TestIterator) = iterate(iter.data)
Base.iterate(iter::TestIterator, state) = iterate(iter.data, state)

Base.IteratorSize(::Type{<:TestIterator{<:Any,IS}}) where {IS} = IS()

Base.IteratorEltype(::Type{<:TestIterator{<:Any,<:Any,IE}}) where {IE} = IE()

Base.eltype(::Type{<:TestIterator{T,<:Any,Base.HasEltype}}) where {T} = eltype(T)

Base.length(iter::TestIterator{<:Any,Base.HasLength}) = length(iter.data)
Base.length(iter::TestIterator{<:Any,<:Base.HasShape}) = length(iter.data)

Base.size(iter::TestIterator{<:Any,<:Base.HasShape}) = size(iter.data)

Base.:(==)(iter1::T, iter2::T) where {T<:TestIterator} = iter1.data == iter2.data

Base.isequal(iter1::T, iter2::T) where {T<:TestIterator} = isequal(iter1.data, iter2.data)

function Base.hash(iter::TestIterator{<:Any,IT,IS}) where {IT,IS}
    return mapreduce(hash, hash, (iter.data, IT, IS))
end

# To make it a valid differential: needs at very least `zero` and `+`
Base.zero(::Type{<:TestIterator}) = ZeroTangent()
function Base.:+(iter1::TestIterator{T,IS,IE}, iter2::TestIterator{T,IS,IE}) where {T,IS,IE}
    return TestIterator{T,IS,IE}(map(+, iter1.data, iter2.data))
end

# For testing purposes:

function rand_tangent(rng::AbstractRNG, x::TestIterator{<:Any,IS,IE}) where {IS,IE}
    ∂data = rand_tangent(rng, x.data)
    return TestIterator{typeof(∂data),IS,IE}(∂data)
end

function FiniteDifferences.to_vec(iter::TestIterator)
    iter_vec, back = to_vec(iter.data)
    function TestIterator_from_vec(v)
        return TestIterator(back(v), Base.IteratorSize(iter), Base.IteratorEltype(iter))
    end
    return iter_vec, TestIterator_from_vec
end


text:
"""
	gmtsimplify(cmd0::String="", arg1=nothing, kwargs...)

Line reduction using the Douglas-Peucker algorithm.

Full option list at [`gmtsimplify`]($(GMTdoc)gmtsimplify.html)

Parameters
----------

- **T** | **tol** | **tolerance** :: [Type => Str | Number]    `Arg = tolerance[unit]`

    Specifies the maximum mismatch tolerance in the user units. If the data is not Cartesian then append the distance unit.
    ($(GMTdoc)gmtsimplify.html#t)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_swap_xy)
"""
function gmtsimplify(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("gmtsimplify", cmd0, arg1)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode

	cmd, = parse_common_opts(d, "", [:V_params :b :d :e :f :g :h :i :o :yx])
	cmd  = add_opt(d, cmd, 'T', [:T :tol :tolerance])

	common_grd(d, cmd0, cmd, "gmtsimplify ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
gmtsimplify(arg1, cmd0::String=""; kw...) = gmtsimplify(cmd0, arg1; kw...)


text:
export NoteSequence, midifile

using MIDI
using MIDI: TimeSignatureEvent, KeySignatureEvent, SetTempoEvent
using Base.Iterators

"""
    SeqNote <: Any

Stores the [`MIDI.Note`](@ref) data along with it's `program` and `instrument`.
"""
mutable struct SeqNote
    pitch::Int
    velocity::Int
    start_time::Int
    end_time::Int
    instrument::Int
    program::Int
end

"""
    Tempo <: Any

Stores the tempo information from a [`MIDI.SetTempoEvent`](@ref).

`qpm` is the tempo in quarter notes per minute.
"""
mutable struct Tempo
    time::Int
    qpm::Float64
end

"""
    PitchBend <: Any

Stores the [`MIDI.PitchBendEvent`](@ref) data along with it's `program` and `instrument`.
"""
mutable struct PitchBend
    time::Int
    bend::Int
    instrument::Int
    program::Int
end

"""
    ControlChange <: Any

Stores the [`MIDI.ControlChangeEvent`](@ref) data along with it's `program` and `instrument`.
"""
mutable struct ControlChange
    time::Int
    controller::Int
    value::Int
    instrument::Int
    program::Int
end

"""
    NoteSequence <: Any

`NoteSequence` contains a symbolic sequence of music represented by data
in the form of notes, control changes, pitch bend, etc.
It closely follows the structure of a MIDI file.

Unlike a MIDI file, each event in a `NoteSequence` also stores the corresponding program
and instrument number.

## Fields
* `tpq::Int=220`: Resolution in ticks per quarter.
* `isquantized::Bool=false`: If true, the sequence is quantized and all the time information
   in the sequence is represented in integer steps. Otherwise, the sequence is unquantized
   and all the time information in the sequence is represented in ticks.
* `sps::Int=-1`: Amount of quantization steps per second (only used if the sequence is quantized).
* `total_time::Int=0`: Total time of the sequence. It is the end time of the final the note.
   It does not count the end time of other events, if they occur after the final note.
* `timesignatures::Vector{TimeSignatureEvent}=[]`: Stores the time signature information
   in the sequence.
* `keysignatures::Vector{KeySignatureEvent}=[]`:  Stores the key signature information
   in the sequence.
* `tempos::Vector{Tempo}=[]`: Stores the tempo (quarter notes per minute) changes in the sequence.
* `notes::Vector{SeqNote}=[]`: Stores the musical notes in the sequence.
* `pitchbends::Vector{PitchBend}=[]`: Stores the `PitchBend` data that's useful for performance
   and re-synthesis.
* `controlchanges::Vector{ControlChange}=[]`: Stores the `ControlChange` data that's useful for
   performance and re-synthesis.
"""
Base.@kwdef mutable struct NoteSequence
    tpq::Int = 220
    isquantized::Bool = false
    sps::Int = -1
    total_time::Int = 0
    timesignatures::Vector{TimeSignatureEvent} = []
    keysignatures::Vector{KeySignatureEvent} = []
    tempos::Vector{Tempo} = []
    notes::Vector{SeqNote} = []
    pitchbends::Vector{PitchBend} = []
    controlchanges::Vector{ControlChange} = []
end

function NoteSequence(tpq::Int; isquantized::Bool = false, sps::Int = -1)
    if !isquantized && sps != -1
        throw(ArgumentError("`sps` must be not be given for an unquantized note sequence"))
    end

    if isquantized && sps <= 0
        throw(ArgumentError("`sps` must be greater than zero for a quantized note sequence"))
    end

    NoteSequence(tpq=tpq, isquantized=isquantized, sps=sps)
end

function Base.getproperty(ns::NoteSequence, sym::Symbol)
    if sym === :steps_per_quarter
        if !ns.isquantized
            steps_per_quarter = -1
        else
            steps_per_quarter = 60.0 / ns.tempos[1].qpm * ns.sps
        end
        return steps_per_quarter
    else
        return getfield(ns, sym)
    end
end

function Base.show(io::IO, ns::NoteSequence)
    print(io, "NoteSequence(tpq=$(ns.tpq), isquantized=$(ns.isquantized), sps=$(ns.sps))\n")
    unit = (ns.isquantized) ? "steps" : "ticks"
    print(io, "  Total time = $(ns.total_time) $unit\n")
    T = length(ns.timesignatures)
    K = length(ns.keysignatures)
    Te = length(ns.tempos)
    print(io, "  $T TimeSignatures, $K KeySignatures, $Te Tempos\n")
    N = length(ns.notes)
    P = length(ns.pitchbends)
    C = length(ns.controlchanges)
    print(io, "  $N Notes, $P PitchBends, $C ControlChanges\n")
end

"""
    NoteSequence(midi::MIDIFile)

Return [`NoteSequence`](@ref) from a `MIDIFile`.
"""
function NoteSequence(midi::MIDIFile)
    midi = deepcopy(midi)
    ns = NoteSequence(Int(midi.tpq), isquantized=false)

    # Load meta info
    # No need to convert midifile to absolute time since tempochanges returns absolute time automatically
    for (time, qpm) in tempochanges(midi)
        push!(ns.tempos, Tempo(time, qpm))
    end

    # Convert midifile to absolute time for the rest of the events
    toabsolutetime!(midi)
    for event in midi.tracks[1].events
        if event isa TimeSignatureEvent
            push!(ns.timesignatures, event)
        elseif event isa KeySignatureEvent
            push!(ns.keysignatures, event)
        end
    end

    instruments = getinstruments(midi, :absolute)
    for (ins_num, ins) in enumerate(instruments)
        for note in ins.notes
            seqnote = SeqNote(note.pitch, note.velocity, note.position,
                              note.position + note.duration,
                              ins_num, ins.program)

            push!(ns.notes, seqnote)

            if ns.total_time == 0 || seqnote.end_time > ns.total_time
                ns.total_time = seqnote.end_time
            end
        end

        for event in ins.pitchbends
            push!(ns.pitchbends, PitchBend(event.dT, event.pitch, ins_num, ins.program))
        end

        for event in ins.controlchanges
            push!(ns.controlchanges, ControlChange(event.dT, event.controller, event.value, ins_num, ins.program))
        end
    end

    return ns
end

"""
    getinstruments(ns::NoteSequence)

Extract [`Instrument`](@ref)s from the [`NoteSequence`](@ref).
"""
function getinstruments(ns::NoteSequence)
    # Create dicts that map (instrument_num, program) to their events
    instrument_notes = DefaultOrderedDict{Tuple{Int64, Int64}, Notes}(Notes)
    instrument_pb = DefaultOrderedDict{Tuple{Int64, Int64}, Vector{PitchBendEvent}}(Vector{PitchBendEvent})
    instrument_cc = DefaultOrderedDict{Tuple{Int64, Int64}, Vector{ControlChangeEvent}}(Vector{ControlChangeEvent})

    for note in ns.notes
        key = (note.instrument, note.program)
        push!(instrument_notes[key], Note(note.pitch, note.velocity, note.start_time, note.end_time - note.start_time))
    end

    for pb in ns.pitchbends
        key = (pb.instrument, pb.program)
        push!(instrument_pb[key], PitchBendEvent(pb.time, pb.pitch))
    end

    for cc in ns.controlchanges
        key = (cc.instrument, cc.program)
        push!(instrument_cc[key], ControlChangeEvent(cc.time, cc.controller, cc.value))
    end

    # Obtain unique keys from each of those dicts
    instrument_info = collect(union(keys.([instrument_notes, instrument_pb, instrument_cc])...))
    sort!(instrument_info)

    instruments = Vector{Instrument}()
    for (ins_num, program) in instrument_info
        ins = Instrument(program=program)
        key = (ins_num, program)
        append!(ins.notes, instrument_notes[key])
        append!(ins.pitchbends, instrument_pb[key])
        append!(ins.controlchanges, instrument_cc[key])
        push!(instruments, ins)
    end

    instruments
end

"""
    midifile(ns::NoteSequence)

Return a MIDIFile from a [`NoteSequence`](@ref)
"""
function midifile(ns::NoteSequence)
    midifile = MIDIFile(1, ns.tpq, MIDITrack[])
    metatrack = MIDITrack()

    isempty(ns.tempos) && push!(ns.tempos, Tempo(0, 120.0))
    for tempo in ns.tempos
        # Convert qpm to microseconds
        μs = Int(6e7 ÷ tempo.qpm)
        push!(metatrack.events, SetTempoEvent(tempo.time, μs))
    end

    # Add default time signature if time signatures are missing
    if isempty(ns.timesignatures)
        push!(ns.timesignatures, TimeSignatureEvent(0, 4, 4, 24, 8))
    end

    append!(metatrack.events, ns.timesignatures)
    append!(metatrack.events, ns.keysignatures)
    push!(midifile.tracks, metatrack)

    instruments = getinstruments(ns)
    instrument_tracks = getmiditracks(instruments)
    append!(midifile.tracks, instrument_tracks)

    midifile
end

"""
    quantizedstep(ticks::Int, tpq::Int, qpm::Float64, sps::Int, cutoff=QUANTIZE_CUTOFF)

Quantizes `ticks` to the nearest steps, given the steps per second (`sps`).
The ticks are converted to seconds using the ticks per quarter (`tpq`)
and quarter notes per minute (`qpm`).
The seconds are then further converted to steps before being quantized.

Notes above the cutoff are rounded up to the closest step
and notes below the cutoff are rounded down to the closest step.
For example,
if 1.0 <= event <= 1.75, it will be quantized to step 1
if 1.75 < event <= 2.0, it will be quantized to step 2
"""
function quantizedstep(ticks::Int, tpq::Int, qpm::Float64, sps::Int, cutoff=QUANTIZE_CUTOFF)
    seconds = ticks * ms_per_tick(tpq, qpm) / 1e3
    steps = seconds * sps

    floor(Int, steps + (1 - cutoff))
end

function quantize_notes!(ns::NoteSequence, sps::Int)
    qpm = ns.tempos[1].qpm

    for note in ns.notes
        note.start_time = quantizedstep(note.start_time, ns.tpq, qpm, sps)
        note.end_time = quantizedstep(note.end_time, ns.tpq, qpm, sps)

        if note.start_time == note.end_time
            note.end_time += 1
        end

        if ns.total_time == -1 || note.end_time > ns.total_time
            ns.total_time = note.end_time
        end

        if note.start_time < 0 || note.end_time < 0
            throw(error("Note has negative time"))
        end
    end

    for cc in ns.controlchanges
        cc.time = quantizedstep(cc.time, ns.tpq, qpm, sps)
        if cc.time < 0
            throw(error("Control change event has negative time"))
        end
    end

    ns.isquantized = true
    ns.sps = sps

    ns
end

"""
    absolutequantize!(ns::NoteSequence, sps::Int)
    absolutequantize(ns::NoteSequence, sps::Int)

Quantize a NoteSequence to absolute time based on the given steps per second (`sps`).
"""
function absolutequantize!(ns::NoteSequence, sps::Int)
    ns.isquantized && throw(error("The NoteSequence is already quantized"))

    if isempty(ns.tempos)
        push!(ns.tempos, Tempo(0, DEFAULT_QPM))
    else
        for tempo in ns.tempos[2:end]
            if tempo.qpm != ns.tempos[1].qpm
                throw(error("The NoteSequence has multiple tempo changes"))
            end
        end
    end

    qpm = ns.tempos[1].qpm

    ns.total_time = quantizedstep(ns.total_time, ns.tpq, qpm, sps)
    quantize_notes!(ns, sps)
end

function absolutequantize(ns::NoteSequence, sps::Int)
    ns = deepcopy(ns)
    absolutequantize!(ns, sps)
end

"""
    relativequantize!(ns::NoteSequence, steps_per_quarter::Int)
    relativequantize(ns::NoteSequence, steps_per_quarter::Int)

Quantize a NoteSequence relative to its tempo. `steps_per_quarter` is the amount of
quantized time steps each quarter note will be divided into.
"""
function relativequantize!(ns::NoteSequence, steps_per_quarter::Int)
    ns.isquantized && throw(error("The NoteSequence is already quantized"))

    if isempty(ns.timesignatures)
        push!(ns.timesignatures, MIDI.TimeSignatureEvent(0, 4, 4, 24, 8))
    else
        timesignatures = sort(ns.timesignatures, by=timesig->timesig.dT)
        if timesignatures[1].dT != 0 &&
                !(timesignatures[1].numerator == 4 && timesignatures[1].denominator == 4)
            throw(error("Found multiple time signature changes due to an implicit time signature"))
        end

        for timesignature in timesignatures[2:end]
            if (timesignature.numerator != ns.timesignatures[1].numerator ||
                    timesignatures.denominator != ns.timesignatures[1].denominator)
                throw(error("Found multiple time signatures in the note sequence."))
            end
        end
    end

    if isempty(ns.tempos)
        push!(ns.tempos, NoteSequences.Tempo(0, DEFAULT_QPM))
    else
        tempos = sort(ns.tempos, by=tempo->tempo.time)
        tempos[1].time != 0 && tempos[1].qpm != DEFAULT_QPM &&
            throw(error("Found multiple tempo changes due to an implicit tempo change."))

        for tempo in tempos
            tempo.qpm != ns.tempos[1].qpm &&
                throw(error("Found multiple tempo changes in the note sequence."))
        end
    end

    qpm = ns.tempos[1].qpm

    # Convert steps per quarter to quantization steps per second
    sps = round(Int, steps_per_quarter * qpm / 60.0)

    ns.total_time =  quantizedstep(ns.total_time, ns.tpq, qpm, sps)
    quantize_notes!(ns, sps)

    ns
end

function relativequantize(ns::NoteSequence, steps_per_quarter::Int)
    ns = deepcopy(ns)
    relativequantize!(ns, steps_per_quarter)
end

"""
    temporalstretch!(ns::NoteSequence, factor::Real)

Apply a constant temporal stretch to a NoteSequence.
`factor` is the amount of stretch to be applied.
If the `factor` is greater than 1.0, it increases the length, which makes the sequence slower.
If the `factor` is lower than 1.0, it decreases the length, which makes the sequence faster.
"""
function temporalstretch!(ns::NoteSequence, factor::Real)
    ns.isquantized && throw(error("Can only stretch unquantized NoteSequence"))

    factor == 1.0 && return ns

    # Stretch notes and total time of the NoteSequence
    for note in ns.notes
        note.start_time = round(note.start_time * factor)
        note.end_time = round(note.end_time * factor)
    end
    ns.total_time = round(ns.total_time * factor)

    # Stretch tempos
    for tempo in ns.tempos
        tempo.qpm /= factor
    end

    # Stretch event times for all other events
    for event in flatten((ns.timesignatures, ns.keysignatures))
        event.dT = round(event.dT * factor)
    end

    for event in flatten((ns.tempos, ns.pitchbends, ns.controlchanges))
        event.time = round(event.time * factor)
    end

    ns
end

"""
    transpose(ns::NoteSequence, amount::Int, minpitch::Int, maxpitch::Int)

Transpose the NoteSequence by `amount` half-steps
and return the transposed sequence along with the number of deleted notes.
The allowed range of values for a pitch is defined by `minpitch` and `maxpitch`.
If the transposed pitch for a note goes out of that range, the note is removed.
"""
function transpose(sequence::NoteSequence, amount::Int, minpitch::Int, maxpitch::Int)
    ns = deepcopy(sequence)
    new_notes = Vector{SeqNote}()
    num_deleted = 0
    end_time = 0

    for note in ns.notes
        note.pitch += amount
        if minpitch <= note.pitch <= maxpitch
            # Keep track of note ending times
            end_time = max(end_time, note.end_time)
            push!(new_notes, note)
        else
            # Don't include this pitch since it's out of range
            num_deleted += 1
        end
    end

    if num_deleted > 0
        ns.notes = new_notes
        # Update total time, since some notes were removed
        ns.total_time = end_time
    end

    # TODO: transpose key signatures

    ns, num_deleted
end

"""
    applysustainchanges!(sequence::NoteSequence, sustaincontrolnumber::Int=64)
    applysustainchanges(sequence::NoteSequence, sustaincontrolnumber::Int=64)

Apply the sustain pedal control changes and return a new NoteSequence.
`sustaincontrolnumber` is the MIDI controller number for the sustain pedal.

For events with the given `sustaincontrolnumber` and values 0-63,
it will be considered as sustain pedal off.
For events with the control number and values 64-127,
it will be considered as sustain pedal on events.
"""
function applysustainchanges!(ns::NoteSequence, sustaincontrolnumber::Int=64)
    ns.isquantized &&
        throw(error("Can only apply sustain changes to an unquantized NoteSequence"))

    # Priority for events (used when sorting)
    _SUS_ON, _SUS_OFF, _NOTE_ON, _NOTE_OFF = (0, 1, 2, 3)

    events = Any[]
    onsets = [(note.start_time, _NOTE_ON, note) for note in ns.notes]
    offsets = [(note.end_time, _NOTE_OFF, note) for note in ns.notes]
    append!(events, [onsets..., offsets...])

    for cc in ns.controlchanges
        if cc.controller == sustaincontrolnumber
            if cc.value >= 64
                push!(events, (cc.time, _SUS_ON, cc))
            elseif cc.value < 64
                push!(events, (cc.time, _SUS_OFF, cc))
            end
        end
    end

    # Sort the list of events by time and priority (to ensure the order of processing)
    sort!(events, by=event->getindex(event, [1, 2]))

    # Map instrument to its currently active notes
    active_notes = DefaultDict{Int, Vector{SeqNote}}(Vector{SeqNote})
    # Map instrument to its sustain status
    sustainactive = DefaultDict{Int, Bool}(()->false)
    end_time = 0
    for (time, event_type, event) in events
        if event_type == _SUS_ON
            sustainactive[event.instrument] = true
        elseif event_type == _SUS_OFF
            sustainactive[event.instrument] = false
            # End all the notes that were being extended for this instrument
            new_active_notes = Vector{SeqNote}()
            for note in active_notes[event.instrument]
                if note.end_time < time
                    # Update this note's end time since it was being sustained
                    note.end_time = time
                    ns.total_time = time
                else
                    # Keep this note since it is still active
                    push!(new_active_notes, note)
                end
            end
            active_notes[event.instrument] = new_active_notes
        elseif event_type == _NOTE_ON
            if sustainactive[event.instrument]
                # If sustain is on,
                # end all the previous notes with the same pitch
                new_active_notes = Vector{SeqNote}()
                for note in active_notes[event.instrument]
                    if note.pitch == event.pitch
                        note.end_time = time
                        if note.start_time == note.end_time
                            # Remove this note since it has no duration
                            # as another note of the same pitch has started at the same time
                            idx = findall(n->n == note)
                            deleteat!(ns.notes, idx)
                        end
                    else
                        push!(new_active_notes, note)
                    end
                end
                active_notes[event.instrument] = new_active_notes
            end
            push!(active_notes[event.instrument], event)
        elseif event_type == _NOTE_OFF
            # If sustain is active, ignore
            # Otherwise, remove this event from active_notes
            if !sustainactive[event.instrument]
                if event in active_notes[event.instrument]
                    idx = findall(e->e == event, active_notes[event.instrument])
                    deleteat!(active_notes[event.instrument], idx)
                end
            end
        end
        end_time = time
    end

    for instrument_notes in values(active_notes)
        for note in instrument_notes
            note.end_time = end_time
            ns.total_time = end_time
        end
    end

    ns
end

function applysustainchanges(sequence::NoteSequence, sustaincontrolnumber::Int=64)
    ns = deepcopy(sequence)
    applysustainchanges!(ns, sustaincontrolnumber)
end

# Extended by sub modules
function encode_event end
function decode_event end


text:
#### EGSHeader, EGSPhspIterator
struct EGSHeader{P<:EGSParticle}
    particlecount::Int32
    photoncount::Int32
    max_E_kin::Float32
    min_E_kin_electrons::Float32
    originalcount::Float32
end

function Base.show(io::IO, o::EGSHeader)
    kwshow(io, o, calle=typeof(o))
end

"""
    get_originalcount(phsp)

Get the number of electrons at the target used to generate the phase space.
"""
get_originalcount(o::EGSHeader)::Float64 = Float64(o.originalcount)

function read_ZLAST(io::IO)
    mode = prod([read(io, Char) for _ in 1:5])
    if mode == "MODE0"
        Nothing
    elseif mode == "MODE2"
        Float32
    else
        error("Unknown mode $mode")
    end
end

function consume_egs_header(io::IO)
    ZLAST = read_ZLAST(io)
    P = EGSParticle{ZLAST}
    nphsp      = read(io, Int32)
    nphotphsp  = read(io, Int32)
    ekmaxphsp  = read(io, Float32)
    ekminphsp  = read(io, Float32)
    nincphsp   = read(io, Float32)
    EGSHeader{P}(nphsp, nphotphsp, ekmaxphsp, ekminphsp, nincphsp)
end

ptype(::Type{EGSHeader{P}}) where {P} = P

struct EGSPhspIterator{P <: EGSParticle, I <:IO} <: AbstractPhspIterator
    io::I
    header::EGSHeader{P}
    buffer::Base.RefValue{P}
    length::Int64
end

Base.eltype(::Type{<:EGSPhspIterator{P}}) where {P} = P

function egs_iterator(io::IO)
    h = consume_egs_header(io)
    total_size  = bytelength(io)
    P = ptype(h)
    float_len = (total_size / sizeof(P)) - 1
    len = floor(Int64, float_len)
    if float_len != len
        @warn "File size indicates $(float_len) particles, which is not an integer. Assuming $(len) particles instead."
    end
    if len != h.particlecount
        @warn "Particle count according to the header is $(h.particlecount), while there are actually $(float_len) particles stored in the file."
    end
    buffer = Base.RefValue{P}()
    EGSPhspIterator(io, h, buffer, len)
end

function Base.iterate(iter::EGSPhspIterator, nread=0)
    # skip header
    if nread == 0
        pos = sizeof(ptype(iter.header))
        seek(iter.io, pos)
    end
    if nread == length(iter)
        return nothing
    else
        p = read!(iter.io, iter.buffer)[]
        nread += 1
        p, nread
    end
end

zlast_type(P::Type{EGSParticle{ZLAST}}) where {ZLAST} = ZLAST
zlast_type(p::EGSParticle{ZLAST}) where {ZLAST} = ZLAST
zlast_type(o) = zlast_type(ptype(o))

function write_header(io::IO, h::EGSHeader)
    ZLAST = zlast_type(h)
    mode = if ZLAST == Nothing
        "MODE0"
    else
        @assert ZLAST == Float32
        "MODE2"
    end

    ret = 0
    ret += write(io, mode)
    ret += write(io, h.particlecount)
    ret += write(io, h.photoncount)
    ret += write(io, h.max_E_kin)
    ret += write(io, h.min_E_kin_electrons)
    ret += write(io, h.originalcount)
    psize = sizeof(ptype(h))
    while (ret < psize)
        ret += write(io, '\0')
    end
    @assert ret == psize
    ret
end

mutable struct EGSWriter{P <: EGSParticle, I <: IO}
    io::I
    particlecount::Int32
    photoncount::Int32
    max_E_kin::Float32
    min_E_kin_electrons::Float32
    originalcount::Float32
    buffer::Base.RefValue{P}
    function EGSWriter{P}(io::I, particlecount, photoncount,
                       max_E_kin, min_E_kin_electrons,
                       originalcount, buffer) where {P,I}

        w = new{P,I}(io, particlecount, photoncount,
                       max_E_kin, min_E_kin_electrons,
                       originalcount, buffer)
        finalizer(close, w)
        w
    end
end

function Base.write(w::EGSWriter{P}, p::P) where {P <: EGSParticle}
    w.particlecount += 1
    if isphoton(p)
        w.photoncount += 1
    end
    w.max_E_kin = max(w.max_E_kin, p.E)
    if iselectron(p)
        w.min_E_kin_electrons = min(w.min_E_kin_electrons, p.E)
    end
    w.buffer[] = p
    write(w.io, w.buffer)
end

function create_header(w::EGSWriter{P}) where {P}
    h = EGSHeader{P}(
        w.particlecount::Int32,
        w.photoncount::Int32,
        w.max_E_kin::Float32,
        w.min_E_kin_electrons::Float32,
        w.originalcount::Float32,
   )
end

"""
    egs_writer(f, path, P::Type{<:EGSParticle})

Write particles in EGS format to `path`:
```jldoctest
julia> using PhaseSpaceIO

julia> p = EGSParticle(
           E=1.0, weight=5.2134085,
           x=2.0, y=-3.0,
           u=0.99884456, v=-0.03811472, w=0.029271236,
           new_history=false, zlast=0.23947382f0, latch=Latch(charge=0));

julia> path = tempname() * ".egsphsp1";

julia> egs_writer(path, typeof(p)) do w
           write(w, p)
       end
32

julia> phsp_iterator(collect, path)
1-element Array{EGSParticle{Float32},1}:
 EGSParticle(latch=Latch(multicross=false, charge=0, creation=0, visited=@BitRegions(), brems=fals
e, ), new_history=false, E=1.0f0, x=2.0f0, y=-3.0f0, u=0.99884456f0, v=-0.03811472f0, w=0.02927123
6f0, weight=5.2134085f0, zlast=0.23947382f0, )
```
"""
function egs_writer(f, path, P)
    w = egs_writer(path, P)
    ret = call_fenced(f,w)
    close(w)
    ret
end

function egs_writer(path::AbstractString, P)
    io = open(path, "w")
    egs_writer(io, P)
end

function egs_writer(io::IO, ::Type{P}) where {P <: EGSParticle}
    buffer = Base.RefValue{P}()
    w = EGSWriter{P}(io, Int32(0),Int32(0),Float32(-Inf),Float32(Inf),Float32(1.), buffer)
    h = create_header(w)
    write_header(io, h)
    w
end

function Base.flush(w::EGSWriter)
    if isopen(w.io)
        h = create_header(w)
        pos = Base.position(w.io)
        seekstart(w.io)
        write_header(w.io, h)
        seek(w.io, pos)
    end
    flush(w.io)
end

function Base.close(w::EGSWriter)
    flush(w)
    close(w.io)
end


text:
module FireBarrier

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Fire Barrier" => Ahorn.EntityPlacement(
        Maple.FireBarrier,
        "rectangle"
    ),
)

Ahorn.minimumSize(entity::Maple.FireBarrier) = 8, 8
Ahorn.resizable(entity::Maple.FireBarrier) = true, true

Ahorn.selection(entity::Maple.FireBarrier) = Ahorn.getEntityRectangle(entity)

edgeColor = (246, 98, 18, 255) ./ 255
centerColor = (209, 9, 1, 102) ./ 255

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.FireBarrier, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, centerColor, edgeColor)
end

end


text:
using Test
using HTTP

import Base.==

==(a::HTTP.Response,b::HTTP.Response) = (a.status  == b.status)    &&
                                        (a.version == b.version)   &&
                                        (a.headers == b.headers)   &&
                                        (a.body    == b.body)

@testset "HTTP.Handler" begin

f = HTTP.HandlerFunction((req) -> HTTP.Response(200))
@test HTTP.handle(f, HTTP.Request()) == HTTP.Response(200)

r = HTTP.Router()
@test length(methods(r.func)) == 1
@test HTTP.handle(r, HTTP.Request()) == HTTP.Response(404)

HTTP.register!(r, "/path/to/greatness", f)
@test length(methods(r.func)) == 2
req = HTTP.Request()
req.target = "/path/to/greatness"
@test HTTP.handle(r, req) == HTTP.Response(200)

p = "/next/path/to/greatness"
f2 = HTTP.HandlerFunction((req) -> HTTP.Response(201))
HTTP.register!(r, p, f2)
@test length(methods(r.func)) == 3
req = HTTP.Request()
req.target = "/next/path/to/greatness"
@test HTTP.handle(r, req) == HTTP.Response(201)

r = HTTP.Router()
HTTP.register!(r, "GET", "/sget", f)
HTTP.register!(r, "POST", "/spost", f)
HTTP.register!(r, "POST", "/tpost", f)
req = HTTP.Request("GET", "/sget")
@test HTTP.handle(r, req) == HTTP.Response(200)
req = HTTP.Request("POST", "/sget")
@test HTTP.handle(r, req) == HTTP.Response(404)
req = HTTP.Request("GET", "/spost")
@test HTTP.handle(r, req) == HTTP.Response(404)
req = HTTP.Request("POST", "/spost")
@test HTTP.handle(r, req) == HTTP.Response(200)
req = HTTP.Request("GET", "/tpost")
@test HTTP.handle(r, req) == HTTP.Response(404)
req = HTTP.Request("POST", "/tpost")
@test HTTP.handle(r, req) == HTTP.Response(200)

r = HTTP.Router()
HTTP.register!(r, "/test", f)
HTTP.register!(r, "/test/*", f2)
f3 = HTTP.HandlerFunction((req) -> HTTP.Response(202))
HTTP.register!(r, "/test/sarv/ghotra", f3)
f4 = HTTP.HandlerFunction((req) -> HTTP.Response(203))
HTTP.register!(r, "/test/*/ghotra/seven", f4)

req = HTTP.Request()
req.target = "/test"
@test HTTP.handle(r, req) == HTTP.Response(200)

req.target = "/test/sarv"
@test HTTP.handle(r, req) == HTTP.Response(201)

req.target = "/test/sarv/ghotra"
@test HTTP.handle(r, req) == HTTP.Response(202)

req.target = "/test/sar/ghotra/seven"
@test HTTP.handle(r, req) == HTTP.Response(203)

end


text:
cd(@__DIR__)
str = read("raw.jl", String);

exprlins = split(str, '\n');

imports = filter(x-> startswith(x, "using"), exprlins)
rest = filter(x-> !startswith(x, "using"), exprlins)
N = 6
open("makie_precompile.jl", "w") do io
    println(io, "using Pkg, Test, LinearAlgebra, Random, Statistics, Dates, BinaryProvider")
    for elem in imports
        println(io, elem)
    end
    for line in rest
        x = Meta.parse(line, raise = true) # is parseable?
        Meta.isexpr(x, :incomplete) && continue
        line = string("try;", line, "; catch; end")
        println(io, line)
    end
end

include("makie_precompile.jl")


text:
using Hexagons

struct HexagonalBinGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement
    tag::Symbol
end
HexagonalBinGeometry(default_statistic; tag=empty_tag) =
        HexagonalBinGeometry(default_statistic, tag)
HexagonalBinGeometry(; xbincount=200, ybincount=200, tag=empty_tag) =
        HexagonalBinGeometry(Gadfly.Stat.hexbin(xbincount=xbincount, ybincount=ybincount), tag)

"""
    Geom.hexbin[(; xbincount=200, ybincount=200)]

Bin the `x` and `y` aesthetics into tiled hexagons and color by count.
`xbincount` and `ybincount` specify the number of bins.  This behavior relies
on the default use of [`Stat.hexbin`](@ref).

Alternatively, draw hexagons of size `xsize` and `ysize` at positions `x` and
`y` by passing [`Stat.identity`](@ref) to `plot` and manually binding the `color`
aesthetic.
"""
const hexbin = HexagonalBinGeometry

default_statistic(geom::HexagonalBinGeometry) = geom.default_statistic
element_aesthetics(geom::HexagonalBinGeometry) = [:x, :y, :xsize, :ysize, :color]

function render(geom::HexagonalBinGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    default_aes = Gadfly.Aesthetics()
    default_aes.color = discretize_make_ia(RGB{Float32}[theme.default_color])
    default_aes.xsize = [1.0]
    default_aes.ysize = [1.0]
    aes = inherit(aes, default_aes)

    Gadfly.assert_aesthetics_defined("Geom.hexbin", aes, :x, :y)
    Gadfly.assert_aesthetics_equal_length("Geom.hexbin", aes, :x, :y)

    n = length(aes.x)
    visibility = Bool[!ismissing(c) for c in takestrict(cycle(aes.color), n)]
    xs = aes.x[visibility]
    ys = aes.y[visibility]
    xsizes = collect(eltype(aes.xsize), takestrict(cycle(aes.xsize), n))[visibility]
    ysizes = collect(eltype(aes.ysize), takestrict(cycle(aes.ysize), n))[visibility]
    cs = collect(eltype(aes.color), takestrict(cycle(aes.color), n))[visibility]
    n = length(xs)

    return compose!(
        context(),
        Compose.polygon([hexpoints(xs[i], ys[i], xsizes[i], ysizes[i]) for i in 1:n], geom.tag),
        fill(cs),
        svgclass("geometry"))
end


text:
@testset "Full codebook update" begin

  d, n, m, h, V, rho = 32, 10_000, 4, 256, false, 1e-4
  X, _, B = generate_random_dataset(Float32, Int16, d, n, m, h)

  C1 = Rayuela.update_codebooks_fast_bin(X, B, h, V, rho)
  C2 = Rayuela.update_codebooks_fast_bin2(X, B, h, V, rho)
  @test isapprox(C1, C2)

end

# Make sure the fast version of codebook update is still okay
@testset "Chain codebook update" begin
  d, n, m, h, V, rho = 32, 10_000, 4, 256, false, 1e-4
  X, _, B = generate_random_dataset(Float64, Int16, d, n, m, h)

  # These two methods are equivalent, but the second should be faster
  C1, _ = Rayuela.update_codebooks_chain(X, B, h, V)
  C2, _ = Rayuela.update_codebooks_chain_bin(X, B, h, V, rho)
  @test isapprox(C1, C2)
end

# Chain quantization
@testset "Chain encoding" begin
  d, n, m, h = 32, Int(1e3), 4, 256
  X, C, B = generate_random_dataset(Float32, Int16, d, n, m, h)

  B1, _ = Rayuela.quantize_chainq(X, C) # Julia

  use_cuda, use_cpp = true, false
  B2, _ = Rayuela.quantize_chainq(X, C, use_cuda, use_cpp) # Cuda

  use_cuda, use_cpp = false, true
  B3, _ = Rayuela.quantize_chainq(X, C, use_cuda, use_cpp) # C++
  @test all(B1 .== B2 .== B3) # C++ implememtation
end


text:
fn = "day20/input.txt"
tiles_text = split(read(fn, String), "\n\n")

struct Tile
    id::Int
    geo::Array{Bool, 2}
end

tiles = []
for tile_text in tiles_text
    tile_lines = split(tile_text, "\n")
    tile_id = parse(Int, match(r"Tile (\d+):", tile_lines[1]).captures[1])
    tile_geo = [[s == '#' for s in tile_line] for tile_line in tile_lines[2:end]]
    tile_geo = vcat(permutedims.(tile_geo)...)
    push!(tiles, Tile(tile_id, tile_geo))
end

function edges(t::Tile)
    return [t.geo[1, 1:end], t.geo[end, 1:end], t.geo[1:end, 1], t.geo[1:end, end]]
end

corners = []
for tile in tiles
    matching_edges = -4
    my_edges = edges(tile)
    for tile2 in tiles
        matching_edges += count(e -> e in my_edges || reverse(e) in my_edges, edges(tile2))
    end
    if matching_edges == 2
        push!(corners, tile.id)
    end
end

println(prod(corners))


text:
using LazyStack

# for stack_iter
LazyStack.no_wraps(a::KeyedArray) = LazyStack.no_wraps(NamedDims.unname(parent(a)))

function LazyStack.rewrap_like(A, a::KeyedArray)
    B = LazyStack.rewrap_like(A, parent(a))
    KeyedArray(B, (axiskeys(a)..., ntuple(d -> keys_or_axes(B, d + ndims(a)), ndims(A) - ndims(a))...))
end

function LazyStack.rewrap_like(A, a::NamedTuple)
    KeyedArray(A, (collect(keys(a)), ntuple(d -> keys_or_axes(A, d + 1), ndims(A) - 1)...))
end

# tuple of arrays
function LazyStack.stack(x::Tuple{Vararg{<:KeyedArray}})
    KeyedArray(LazyStack.stack(map(parent, x)), stack_keys(x))
end

stack_keys(xs::Tuple{Vararg{<:KeyedArray}}) =
    (keys_or_axes(first(xs))..., Base.OneTo(length(xs)))

# array of arrays: first strip off outer containers...
function LazyStack.stack(xs::KeyedArray{<:AbstractArray})
    KeyedArray(stack(parent(xs)), stack_keys(xs))
end
function LazyStack.stack(xs::KeyedArray{<:AbstractArray,N,<:NamedDimsArray{L}}) where {L,N}
    data = stack(parent(parent(xs)))
    KeyedArray(LazyStack.ensure_named(data, LazyStack.getnames(xs)), stack_keys(xs))
end
function LazyStack.stack(xs::NamedDimsArray{L,<:AbstractArray,N,<:KeyedArray}) where {L,N}
    data = stack(parent(parent(xs)))
    LazyStack.ensure_named(KeyedArray(data, stack_keys(xs)), LazyStack.getnames(xs))
end

# ... then deal with inner ones:
function LazyStack.stack(x::AT) where {AT <: AbstractArray{<:KeyedArray{T,IN},ON}} where {T,IN,ON}
    KeyedArray(LazyStack.Stacked{T, IN+ON, AT}(x), stack_keys(x))
end
function LazyStack.stack(x::AT) where {AT <: AbstractArray{<:KeyedArray{T,IN,IT},ON}} where {IT<:NamedDimsArray} where {T,IN,ON}
    data = LazyStack.Stacked{T, IN+ON, AT}(x)
    KeyedArray(NamedDimsArray(data, LazyStack.getnames(x)), stack_keys(x))
end
function LazyStack.stack(x::AT) where {AT <: AbstractArray{<:NamedDimsArray{L,T,IN,IT},ON}} where {IT<:KeyedArray} where {T,IN,ON,L}
    data = LazyStack.Stacked{T, IN+ON, AT}(x)
    NamedDimsArray(KeyedArray(data, stack_keys(x)), LazyStack.getnames(x))
end

stack_keys(xs::AbstractArray{<:AbstractArray}) =
    (keys_or_axes(first(xs))..., keys_or_axes(xs)...)

function LazyStack.getnames(xs::KeyedArray{<:AbstractArray,N,OT}) where {N,OT}
    (dimnames(eltype(xs))..., dimnames(OT)...)
end
function LazyStack.getnames(xs::AbstractArray{<:KeyedArray{T,N,IT}}) where {T,N,IT}
    out_names = hasnames(xs) ? dimnames(xs) : NamedDims.dimnames(xs)
    (NamedDims.dimnames(IT)..., out_names...)
end


text:
import Compose: compose, context, polygon

function sierpinski(n)
    if n == 0
        Compose.compose(context(), polygon([(1,1), (0,1), (1/2, 0)]))
    else
        t = sierpinski(n - 1)
        Compose.compose(context(),
                (context(1/4,   0, 1/2, 1/2), t),
                (context(  0, 1/2, 1/2, 1/2), t),
                (context(1/2, 1/2, 1/2, 1/2), t))
    end
end

function main(window)
    push!(window.assets, "widgets")

    iterᵗ=Signal(0)

    vbox(title(2, "Sierpinski's Triangle"),
        vskip(1em),
        hbox("Iterations: ", slider(0:6) >>> iterᵗ),
        vskip(1em),
        map(iterᵗ) do iter
            sierpinski(iter)
        end
    ) |> Escher.pad(2em)
end


text:
module Halite

export Game, set_logger, DockedStatus, start_game, thrust, dock, undock, update_map,
    all_players, all_planets, all_ships, all_docked_ships, get_player, get_me, get_planet,
    get_docked_ship, isdocked, isowned, isfull, can_dock, navigate, send_command_queue,
    closest_point_to, nearest_unoccupied_planet

@enum DockedStatus UNDOCKED=0 DOCKING=1 DOCKED=2 UNDOCKING=3

include("constants.jl")
include("types.jl")
include("parse.jl")
include("utils.jl")
include("commands.jl")
include("logging.jl")

end #module

max_stars_repo_name: Cerbes/Halite-II

max_issues_repo_name: Cerbes/Halite-II

max_forks_repo_name: Cerbes/Halite-II


text: println("This should ring a bell.\a")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Each base as 'one hot' nucleotide encoding
#
# A   = 0b0001
# C   = 0b0010
# G   = 0b0100
# U   = 0b1000
# Gap = 0b0000

# Pairs encoded as adjacent sets of 'one hot' values
# so AA would be 0b00010001, and AU  0b00011000 and so forth.

abstract type NucleotidePair end

primitive type RNAPair<:NucleotidePair 8 end

const AU_PAIR = reinterpret(RNAPair, 0b00011000)
const UA_PAIR = reinterpret(RNAPair, 0b10000001)
const CG_PAIR = reinterpret(RNAPair, 0b00100100)
const GC_PAIR = reinterpret(RNAPair, 0b01000010)
const GU_PAIR = reinterpret(RNAPair, 0b01001000)
const UG_PAIR = reinterpret(RNAPair, 0b10000100)

primitive type RNAMismatch<:NucleotidePair 8 end

const AA_MISMATCH = reinterpret(RNAMismatch, 0b00010001)
const AC_MISMATCH = reinterpret(RNAMismatch, 0b00010010)
const AG_MISMATCH = reinterpret(RNAMismatch, 0b00010100)
const CA_MISMATCH = reinterpret(RNAMismatch, 0b00100001)
const CC_MISMATCH = reinterpret(RNAMismatch, 0b00100010)
const CU_MISMATCH = reinterpret(RNAMismatch, 0b00101000)
const GA_MISMATCH = reinterpret(RNAMismatch, 0b01000001)
const GG_MISMATCH = reinterpret(RNAMismatch, 0b01000100)
const UC_MISMATCH = reinterpret(RNAMismatch, 0b10000010)
const UU_MISMATCH = reinterpret(RNAMismatch, 0b10001000)

primitive type RNABulge<:NucleotidePair 8 end

const AB_BULGE = reinterpret(RNABulge, 0b00010000)
const CB_BULGE = reinterpret(RNABulge, 0b00100000)
const GB_BULGE = reinterpret(RNABulge, 0b01000000)
const UB_BULGE = reinterpret(RNABulge, 0b10000000)

const BA_BULGE = reinterpret(RNABulge, 0b00000001)
const BC_BULGE = reinterpret(RNABulge, 0b00000010)
const BG_BULGE = reinterpret(RNABulge, 0b00000100)
const BU_BULGE = reinterpret(RNABulge, 0b00001000)

primitive type InvalidPair<:NucleotidePair 8 end

const XX_INVALID = reinterpret(InvalidPair, 0b11111111)

const XA_INVALID = reinterpret(InvalidPair, 0b11110001)
const XC_INVALID = reinterpret(InvalidPair, 0b11110010)
const XG_INVALID = reinterpret(InvalidPair, 0b11110100)
const XU_INVALID = reinterpret(InvalidPair, 0b11111000)

const AX_INVALID = reinterpret(InvalidPair, 0b00011111)
const CX_INVALID = reinterpret(InvalidPair, 0b00101111)
const GX_INVALID = reinterpret(InvalidPair, 0b01001111)
const UX_INVALID = reinterpret(InvalidPair, 0b10001111)

primitive type InvalidBulge<:NucleotidePair 8 end

const BX_INVALID = reinterpret(InvalidBulge, 0b00001111)
const XB_INVALID = reinterpret(InvalidBulge, 0b11110000)

const BULGES = [AB_BULGE, CB_BULGE,
                GB_BULGE, UB_BULGE,
                BA_BULGE, BC_BULGE,
                BG_BULGE, BU_BULGE]

const MISMATCHES = [AA_MISMATCH, AC_MISMATCH,
                    AG_MISMATCH, CA_MISMATCH,
                    CC_MISMATCH, CU_MISMATCH,
                    GA_MISMATCH, GG_MISMATCH,
                    UC_MISMATCH, UU_MISMATCH]

const PAIRS = [AU_PAIR, UA_PAIR,
               CG_PAIR, GC_PAIR,
               GU_PAIR, UG_PAIR]

const INVALIDS = [XA_INVALID, XC_INVALID,
                  XG_INVALID, XU_INVALID,
                  AX_INVALID, CX_INVALID,
                  GX_INVALID, UX_INVALID,
                  XX_INVALID]

Base.convert{NP <: NucleotidePair}(::Type{NP}, x::UInt8) = reinterpret(NP, x)
Base.convert(::Type{UInt8}, x::NucleotidePair)           = reinterpret(UInt8, x)

function Base.convert{NP <: NucleotidePair, X <: NucleicAcid, Y <: NucleicAcid}(::Type{NP}, first::X,
                                                                                last::Y)
   enc  = reinterpret(UInt8, first) << 4
   enc |= reinterpret(UInt8, last)
   reinterpret(NP, enc)
end

function Base.convert{NP <: NucleotidePair}(::Type{NP}, first::UInt8, last::UInt8 )
   enc  = first << 4
   enc |= last
   reinterpret(NP, enc)
end

# convert two adjacent 'one hot' nucleotide 4-bit codes
# into two adjacent 2-bit encoding 0b00 = A, 0b01 = C, 0b10 = G, 0b11 = U
# for both left and right nucleotides, so 0b00010001 becomes 0b0000
# and 0b10001000 becomes 0b1111, then add 1 to bring to 1:16 range
function index(x::NucleotidePair)
   idx   = trailing_zeros(convert(UInt8, x) >> 4) << 2
   right = trailing_zeros(convert(UInt8, x) & 0x0F)
   (idx | right) + 1
end

function index(x::RNABulge)
   trailing_zeros(reinterpret(UInt8, x)) + 1
end

isgap{N <: NucleicAcid}( nuc::N ) = reinterpret(UInt8, nuc) == 0 ? true : false

isbulge(x::RNAPair)     = false
isbulge(x::RNAMismatch) = false
isbulge(x::RNABulge)    = true

ispair(x::RNAPair)     = true
ispair(x::RNAMismatch) = false
ispair(x::RNABulge)    = false

ismismatch(x::RNAPair)     = false
ismismatch(x::RNAMismatch) = true
ismismatch(x::RNABulge)    = false

function type_cnt{NP <: NucleotidePair}( v::Vector{NP}, func::Function )
   cnt = 0
   @inbounds for i in 1:length(v)
      func(v[i]) && (cnt += 1)
   end
   cnt
end

npairs{NP <: NucleotidePair}( v::Vector{NP} )      = type_cnt( v, ispair )
nmismatches{NP <: NucleotidePair}( v::Vector{NP} ) = type_cnt( v, ismismatch )
nbulges{NP <: NucleotidePair}( v::Vector{NP} )     = type_cnt( v, isbulge )

isfiveprime(x::RNABulge) = reinterpret(UInt8, x) & 0x0F == 0 ? true : false

isbulgefive(x::RNABulge)    = isfiveprime(x)
isbulgefive(x::RNAPair)     = false
isbulgefive(x::RNAMismatch) = false

isbulgethree(x::RNABulge)    = !isfiveprime(x)
isbulgethree(x::RNAPair)     = false
isbulgethree(x::RNAMismatch) = false

function flip{NP <: NucleotidePair}(x::NP)
   left   = convert(UInt8, x) >> 4
   retval = (convert(UInt8, x) & 0x0F) << 4
   reinterpret( NP, retval | left )
end

function Base.split(x::NucleotidePair)
   left  = trailing_zeros(convert(UInt8, x) >> 4) + 1
   right = trailing_zeros(convert(UInt8, x) & 0x0F) + 1
   (left, right)
end

function Base.split(x::RNABulge)
   x_int = reinterpret(UInt8, x)
   isfiveprime(x) ? trailing_zeros(x_int >> 4) + 1 : trailing_zeros(x_int) + 1
end

function Base.split(x::InvalidPair)
   left  = convert(UInt8, x) >> 4
   lret  = left == 15 ? 15 : trailing_zeros(left) + 1
   right = convert(UInt8, x) & 0x0F
   rret  = right == 15 ? 15 : trailing_zeros(right) + 1
   (lret, rret)
end

# this function returns the first, or last if flag is true
# five prime or three prime indexes (1:4) in a path of NucleotidePairs
function five_three{NP <: NucleotidePair}( path::Vector{NP}, range::UnitRange;
                                           last::Bool=false )
   five  = 0
   three = 0
   for i in (last ? reverse(range) : range)
      if isa( path[i], RNABulge )
         if isfiveprime( path[i] )
            five = split( path[i] )
         else
            three = split( path[i] )
         end
      else
         x,y   = split( path[i] )
         five  = five > 0 ? five : x
         three = three > 0 ? three : y
      end
      if five > 0 && three > 0
         break
      end
   end
   five, three
end

function is_purine_pyrimidine(x::RNAPair)
   if x == AU_PAIR || x == GC_PAIR || x == GU_PAIR
      return true
   else
      return false
   end
end

const PAIR_DICT = Dict{Tuple{Int,Int},NucleotidePair}()

for i in [PAIRS; MISMATCHES; INVALIDS]
   PAIR_DICT[split(i)] = i
end
for i in BULGES
   if isfiveprime(i)
      PAIR_DICT[(0,split(i))] = i
   else
      PAIR_DICT[(split(i),0)] = i
   end
end


text:
export sync, sync!

function get_sync_t(s::Union{String,DateTime}, t::Array{Int64,1})
  isa(s, DateTime) && return floor(Int64, d2u(s)*sμ)
  if s == "first"
    return minimum(t)
  elseif s == "last"
    return maximum(t)
  else
    return floor(Int64, d2u(DateTime(s))*sμ)
  end
end

function get_sync_inds(t::AbstractArray{Int64,1}, Ω::Bool, t₀::Int64, t₁::Int64)
  if Ω == true
    return union(findall(t.<t₀), findall(t.>t₁))
  else
    return findall(t.<t₀)
  end
end

@doc """
sync!(S::GphysData)

Synchronize the start times of all data in S to begin at or after the last
start time in S.

sync!(S, [s=ST, t=ET, v=VV])

Synchronize all data in S to start at `ST` and terminate at `EN` with verbosity level VV.

For regularly-sampled channels, gaps between the specified and true times
are filled with the mean; this isn't possible with irregularly-sampled data.

#### Specifying start time
* s="last": (Default) sync to the last start time of any channel in `S`.
* s="first": sync to the first start time of any channel in `S`.
* A numeric value is treated as an epoch time (`?time` for details).
* A DateTime is treated as a DateTime. (see Dates.DateTime for details.)
* Any string other than "last" or "first" is parsed as a DateTime.

#### Specifying end time (t)
* t="none": (Default) end times are not synchronized.
* t="last": synchronize all channels to end at the last end time in `S`.
* t="first" synchronize to the first end time in `S`.
* numeric, datetime, and non-reserved strings are treated as for `-s`.

Related functions: time, Dates.DateTime, parsetimewin
""" sync!
function sync!(S::GphysData;
                s="last"::Union{String,DateTime},
                t="none"::Union{String,DateTime},
                v::Int64=KW.v,
                )

  # Dekete empty traces
  deleteat!(S, findall([isempty(S.x[i]) for i =1:S.n]))     # delete empty channels
  S.n == 0 && return nothing                                # pointless to continue

  do_end = t=="none" ? false : true

  # Do not edit order of operations -------------------------------------------
  start_times = zeros(Int64, S.n)
  if do_end
    end_times = zeros(Int64, S.n)
  end

  # non-timeseries data
  irr = falses(S.n)
  non_ts = findall(S.fs .== 0)
  irr[non_ts] .= true
  z = zero(Int64)
  for i = 1:S.n
    start_times[i] = S.t[i][1,2]
    if do_end
      end_times[i] = S.t[i][end,2] + (irr[i] ? z : start_times[i] + sum(S.t[i][2:end,2]) +
        round(Int64, (length(S.x[i])-1)/(μs*S.fs[i])))
        # ts data: start time in μs from epoch + sum of time gaps in μs + length of trace in μs
        # non-ts data: time of last sample
    end
  end
  # Do not edit order of operations -------------------------------------------

  # Start and end times
  t_start = get_sync_t(s, start_times)
  sstr = string(u2d(t_start*μs))

  if do_end
    t_end = get_sync_t(t, end_times)
    t_end > t_start || error("No time overlap with given start & end times!")
    tstr = string(u2d(t_end*μs))
    if v > 0
      @info(@sprintf("Synchronizing %.2f seconds of data\n", (t_end - t_start)*μs))
      if v > 1
        @info(string("t_start = ", u2d(t_start*μs)))
        @info(string("t_end = ", tstr))
    end
    elseif v > 0
      @info(string("Synchronizing to start at ", u2d(t_start*μs)))
    end
  else
    t_end = z
  end

  # Loop over non-timeseries data
  dflag = falses(S.n)
  for i in non_ts
    t = view(S.t[i], :, 2)
    Lt = length(t)

    j = get_sync_inds(t, do_end, t_start, t_end)
    Lj = length(j)
    if Lj ≥ Lt
      dflag[i] = true
      continue
    elseif do_end
      note!(S, i, @sprintf("sync!, s = %s, t = %s, removed %i samples (out of time range) from :x", sstr, tstr, Lj))
    else
      note!(S, i, @sprintf("sync! s = %s, t = none, removed %i samples (out of time range) from :x", sstr, Lj))
    end
    ti = collect(1:Lt)
    deleteat!(S.x[i], j)
    deleteat!(ti, j)
    S.t[i] = S.t[i][ti,:]
  end

  # Loop over timeseries data
  for i = 1:S.n
    if irr[i] == false
      sync_str = Array{String,1}(undef,0)
      desc_str = ""

      # truncate X to values within bounds
      t = t_expand(S.t[i], S.fs[i])
      j = get_sync_inds(t, do_end, t_start, t_end)
      deleteat!(S.x[i], j)
      deleteat!(t, j)

      #= length 0 traces _can_ happen with resampled :x where (lx*f_rat) < 0.5
      requries short series of high-frequency data resampled to too low fs =#
      if length(S.x[i]) == 0
        dflag[i] = true
        continue
      end

      if !isempty(j)
        push!(sync_str, string("deleted ", length(j), " samples from :x"))
      end

      S.t[i] = t_collapse(t, S.fs[i])

      # prepend points to time series data that begin late
      T = eltype(S.x[i])
      fμs = S.fs[i]*μs
      dtμ = round(Int64, 1/fμs)
      μ_x = T(mean(S.x[i]))

      if (start_times[i] - t_start) ≥ dtμ
        ni = div(start_times[i]-t_start, dtμ)
        prepend!(S.x[i], ones(T, ni).*μ_x)

        # corrected 2019-02-28
        S.t[i][1,2] = S.t[i][1,2] - ni*dtμ

        push!(sync_str, string("prepended ", ni, " samples to :x."))
      end

      # append points to time series data that end early
      if do_end
        end_times[i] = S.t[i][1,2] +
                        sum(S.t[i][2:end,2]) +
                        round(Int64, length(S.x[i])/fμs)
        if (t_end - end_times[i]) ≥ dtμ
          ni = div(t_end-end_times[i], dtμ)
          append!(S.x[i], ones(T,ni).*μ_x)

          push!(sync_str, string("appended ", ni, " samples to :x."))
        end

        # Correct for length aberration if necessary
        S.t[i][end,1] = length(S.x[i])

        desc_str = string("sync!, s = ", sstr, ", t = ", tstr)
      else
        desc_str = string("sync!, s = ", sstr, ", t = none")
      end

      if length(sync_str) > 0
        desc_str *= string(", ", join(sync_str, ";"))
      end
      note!(S, i, desc_str)
    end
  end
  del_flagged!(S, dflag, "length 0 after sync")
  return nothing
end

function sync!(C::SeisChannel;
                s::DateTime,
                t="none"::Union{String,DateTime},
                v::Int64=KW.v )
  S = SeisData(C)
  sync!(S, s=s, t=t, v=v)
  return nothing
end

@doc (@doc sync!)
function sync(S::GphysData;
                s="last"::Union{String,DateTime},
                t="none"::Union{String,DateTime},
                v::Int64=KW.v )

  T = deepcopy(S)
  sync!(T, s=s, t=t, v=v)
  return T
end

function sync(C::SeisChannel;
                s::DateTime,
                t="none"::Union{String,DateTime},
                v::Int64=KW.v )

  U = deepcopy(C)
  S = SeisData(U)
  sync!(S, s=s, t=t, v=v)
  return S[1]
end


text:
"""
    SeisFKFilter(d ; <keyword arguments>)

Decimate a multidimensional array input. Input and output have the same dimension

# Arguments
* `d`: 2D Input data.

# Keyword arguments
* `dt=0.002`: time sampling interval.
* `dx=10`: space sampling interval.
* `va=-2000,vb=-3000,vc=3000,vd=2000`: corner velocities to be filtered

# Output
* `out`: Filtered data.

"""

function SeisFKFilter(d;dt=0.002,dx=10,va=-2000,vb=-3000,vc=3000,vd=2000)

    nt = size(d,1)
    nx = size(d,2)
    pa = 1./va
    pb = 1./vb
    pc = 1./vc
    pd = 1./vd
    m = fft(d,1)
    m = fft(m,2)
    nf = size(m,1)
    dw = 1./nf/dt
    nw = isodd(nf) ? round(Int,floor(nf/2)) + 1 : round(Int,floor(nf/2))
    nk = size(m,2)
    dk = 1./nk/dx
    for iw=1:nw
	w = dw*(iw-1)
	for ik=1:nk
	    k = ik < round(Int,floor(nk/2)) + 1 ? dk*(ik-1) : -(dk*nk - dk*(ik-1))
	    p = k/w;
	    if p < pa
		m[iw,ik] *= 0.
	    elseif p >= pa && p < pb
		m[iw,ik] *= (p-pa)/(pb-pa)
	    elseif p >= pb && p <= pc
				m[iw,ik] *= 1.
	    elseif p > pc && p <= pd
		m[iw,ik] *= 1. - (p-pc)/(pd-pc)
	    elseif p > pd
				m[iw,ik] *= 0.
	    end
	end
    end
    # symmetries
    m = ifft(m,2)
    for iw=nw+1:nf
	m[iw,:] = conj(m[nf-iw+2,:])
    end
    d = real(ifft(m,1))

    return d[1:nt,1:nx];
end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
module UDUnits
import Base:+, -, *, / , ^, inv, log, log10, √, show, cconvert

# package code goes here
const depfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(depfile)
    include(depfile)
else
    error("libudunits2 not properly installed. Please run Pkg.build(\"UDUnits\")")
end

using Libdl

const UT_encoding_t = Cint
const UT_UTF8 = 2
const UT_NAMES = 4
const UT_DEFINITION = 8

const buffer_size = 100

# default encoding is UTF8
const UT_ENC = UT_UTF8

# Unit System

mutable struct System
    ptr :: Ptr{Nothing}
end

cconvert(::Type{Ptr{Nothing}}, sys::System) = sys.ptr
_free_system(system::System) = ccall((:ut_free_system,libudunits2),Ptr{Nothing},(Ptr{Nothing},),system)

"""
    system = System()

Load the default SI unit system. `system` behaves similar to a dictionary as it
supports indexing and the function `haskey` (to determine if a unit is a valid):

```julia
using UDUnits
sys = System()
m = sys["meter"]
haskey(sys,"μm") # returns true
```

"""
function System()
    sys = System(ccall((:ut_read_xml,libudunits2),Ptr{Nothing},(Ptr{Nothing},),C_NULL))
    finalizer(_free_system,sys)
    return sys
end

# Unit

mutable struct Unit
    ptr :: Ptr{Nothing}
end

cconvert(t::Type{Ptr{Nothing}}, unit::Unit) = unit.ptr

_free_unit(unit::Unit) = ccall((:ut_free,libudunits2),Ptr{Nothing},(Ptr{Nothing},),unit)
_ut_parse(system::System,unit::AbstractString) =
    ccall((:ut_parse ,libudunits2),
          Ptr{Nothing},(Ptr{Nothing},Ptr{UInt8},UT_encoding_t),system,unit,UT_ENC)

function Unit(system::System,unit::AbstractString)
    ptr = _ut_parse(system,unit)
    if ptr == C_NULL
        error("UDUnits cannot parse $(unit)")
    end
    unit = Unit(ptr)
    finalizer(_free_unit,unit)
    return unit
end

"""
    s = format(unit::Unit; names = false, definition = false)

Format the unit `unit` as a string. If names is true, then definition uses
the unit names (e.g. meter) instead of symbols (e.g. m). If definition is
true, then the unit should be expressed in terms of basic units (m²·kg·s⁻³
instead of W).

"""
function format(unit::Unit; names = false, definition = false)
    buffer = Vector{UInt8}(undef,buffer_size)
    opts = UT_ENC

    if names
        opts |= UT_NAMES
    end

    if definition
        opts |= UT_DEFINITION
    end

    len = ccall((:ut_format,libudunits2),Cint,(Ptr{Nothing},Ptr{UInt8},Csize_t,Cint),unit,buffer,length(buffer),opts)
    if len == -1
        return "unknown"
    else
        return unsafe_string(pointer(buffer))
    end
end

Base.string(unit::Unit) = format(unit)

function Base.show(io::IO,unit::Unit)
    print(io,"<Unit: ")
    printstyled(io, format(unit); color = :blue)
    print(io,">\n")
end

Base.haskey(system::System,unit::AbstractString) = _ut_parse(system,unit) != C_NULL
Base.getindex(system::System,unit::AbstractString) = Unit(system,unit)

"""
    s = symbol(unit::Unit)

Returns the symbol of `unit`.
"""
symbol(unit::Unit) =
    unsafe_string(ccall((:ut_get_symbol,libudunits2),Ptr{UInt8},(Ptr{Nothing},UT_encoding_t),unit,UT_ENC))

"""
    s = name(unit::Unit)

Returns the name of `unit`.
"""
name(unit::Unit) =
    unsafe_string(ccall((:ut_get_name,libudunits2),Ptr{UInt8},(Ptr{Nothing},UT_encoding_t),unit,UT_ENC))

+(unit::Unit,offset::Number) = Unit(ccall((:ut_offset,libudunits2),Ptr{Nothing},(Ptr{Nothing},Float64),unit,offset))
-(unit::Unit,offset::Number) = unit+(-offset)
*(factor::Number,unit::Unit) = Unit(ccall((:ut_scale,libudunits2),Ptr{Nothing},(Float64,Ptr{Nothing}),factor,unit))
*(unit1::Unit,unit2::Unit) = Unit(ccall((:ut_multiply,libudunits2),Ptr{Nothing},(Ptr{Nothing},Ptr{Nothing}),unit1,unit2))
inv(unit::Unit) = Unit(ccall((:ut_invert,libudunits2),Ptr{Nothing},(Ptr{Nothing},),unit))
/(unit1::Unit,unit2::Unit) = Unit(ccall((:ut_divide,libudunits2),Ptr{Nothing},(Ptr{Nothing},Ptr{Nothing}),unit1,unit2))
/(factor::Number,unit::Unit) = factor * inv(unit)

root(unit::Unit,power::Integer) = Unit(ccall((:ut_root,libudunits2),Ptr{Nothing},(Ptr{Nothing},Cint),unit,power))

^(unit::Unit,power::Integer) = Unit(ccall((:ut_raise,libudunits2),Ptr{Nothing},(Ptr{Nothing},Cint),unit,power))
^(unit::Unit,power::Rational{<:Integer}) = root(unit^numerator(power),denominator(power))

√(unit::Unit) = root(unit,2)

log(base::Number,unit::Unit) = Unit(ccall((:ut_log,libudunits2),Ptr{Nothing},(Float64,Ptr{Nothing}),base,unit))
log(unit::Unit) = log(ℯ,unit)
log10(unit::Unit) = log(10,unit)

# Converter

mutable struct Converter
    ptr :: Ptr{Nothing}
end

cconvert(t::Type{Ptr{Nothing}}, conv::Converter) = conv.ptr
_free_converter(converter::Converter) = ccall((:cv_free,libudunits2),Ptr{Nothing},(Ptr{Nothing},),converter)

"""
    converter = Converter(from_unit::Unit,to_unit::Unit)

Creates a converter function of numeric values in the `from_unit` to equivalent
values in the `to_unit`.

```julia
using UDUnits
sys = System()
conv = Converter(sys["m/s"],sys["km/h"])
speed_in_m_per_s = 100.
speed_in_km_per_h = conv(speed_in_m_per_s)
```
"""
function Converter(from_unit::Unit,to_unit::Unit)
    if !areconvertible(from_unit,to_unit)
        error("UDUnits cannot convert from $(from_unit) to $(to_unit)")
    end

    ptr = ccall((:ut_get_converter ,libudunits2),
                Ptr{Nothing},(Ptr{Nothing},Ptr{Nothing}),from_unit,to_unit)

    # ptr should not be null since the units are convertible
    @assert ptr != C_NULL

    converter = Converter(ptr)

    finalizer(_free_converter,converter)
    return converter
end

"""
    bool = areconvertible(from_unit::Unit,to_unit::Unit)

Return true if `from_unit` and `to_unit` are convertible, otherwise false.
"""
areconvertible(unit1::Unit,unit2::Unit) = ccall((:ut_are_convertible ,libudunits2),Cint,(Ptr{Nothing},Ptr{Nothing}),unit1,unit2) == 1

convert(conv::Converter,v::Float64) = ccall((:cv_convert_double,libudunits2),Float64,(Ptr{Nothing},Float64),conv,v)
(conv::Converter)(v::Number) = convert(conv,v)

"""
    s = expression(conv::Converter; variable = "x")
Return a string representation of the converter `conv` using the
variable `variable`.
"""
function expression(conv::Converter; variable = "x")
    buffer = Vector{UInt8}(undef,buffer_size)
    len = ccall((:cv_get_expression,libudunits2),Cint,(Ptr{Nothing},Ptr{UInt8},Csize_t,Ptr{UInt8}),conv,buffer,length(buffer),variable)

    @assert len != -1

    return unsafe_string(pointer(buffer))
end

function Base.show(io::IO,conv::Converter)
    print(io,"<Converter: ")
    printstyled(io, expression(conv); color = :blue)
    print(io,">\n")
end

export System, Unit, Converter, areconvertible, name, symbol, expression, format

end # module