text:
## ComponentArrays
# Similar
Base.similar(x::ComponentArray) = ComponentArray(similar(getdata(x)), getaxes(x)...)
Base.similar(x::ComponentArray, ::Type{T}) where T = ComponentArray(similar(getdata(x), T), getaxes(x)...)
# Base.similar(x::ComponentArray, ::Type{T}, ax::Tuple{Vararg{Int64,N}}) where {T,N} = similar(x, T, ax...)
# function Base.similar(x::ComponentArray, ::Type{T}, ax::Union{Integer, Base.OneTo}...) where T
#     A = similar(getdata(x), T, ax...)
#     if size(getdata(x)) == size(A)
#         return ComponentArray(A, getaxes(x))
#     else
#         return A
#     end
# end
function Base.similar(x::ComponentArray{T1,N,A,Ax}, ::Type{T}, dims::NTuple{N,Int}) where {T,T1,N,A,Ax}
    arr = similar(getdata(x), T, dims)
    return ComponentArray(arr, getaxes(x))
end
function Base.similar(x::ComponentArray{T1,N1,A,Ax}, ::Type{T}, dims::NTuple{N2,Int}) where {T,T1,N1,N2,A,Ax}
    return similar(getdata(x), T, dims)
end

## TODO: write length method for AbstractAxis so we can do this?
    # function Base.similar(::Type{CA}) where CA<:ComponentArray{T,N,A,Axes} where {T,N,A,Axes}
#     axs = getaxes(CA)
#     return ComponentArray(similar(A, length.(axs)...), axs...)
# end

Base.zero(x::ComponentArray) = zero.(x)

## FIXME: waiting on similar(::Type{<:ComponentArray})
# Base.zeros(CA::Type{<:ComponentArray}) = (similar(CA) .= 0)

# Base.ones(CA::Type{<:ComponentArray}) = (similar(CA) .= 1)

# Copy
Base.copy(x::ComponentArray) = ComponentArray(copy(getdata(x)), getaxes(x))

Base.copyto!(dest::AbstractArray, src::ComponentArray) = copyto!(dest, getdata(src))
function Base.copyto!(dest::ComponentArray, src::AbstractArray)
    copyto!(getdata(dest), src)
    return dest
end
function Base.copyto!(dest::ComponentArray, src::ComponentArray)
    copyto!(getdata(dest), getdata(src))
    return dest
end

Base.deepcopy(x::ComponentArray) = ComponentArray(deepcopy(getdata(x)), getaxes(x))

function Base.convert(::Type{ComponentArray{T,N,AA,Ax}}, A::AbstractArray) where {T,N,AA,Ax}
    return ComponentArray{Ax}(A)
end
function Base.convert(::Type{ComponentArray{T,N,A,Ax1}}, x::ComponentArray{T,N,A,Ax2}) where {T,N,A,Ax1,Ax2}
    return x
end
function Base.convert(::Type{ComponentArray{T1,N,A1,Ax1}}, x::ComponentArray{T2,N,A2,Ax2}) where {T1,T2,N,A1,A2,Ax1,Ax2}
    return T1.(x)
end
Base.convert(::Type{<:Array}, x::ComponentArray) = convert(Array, getdata(x))

# Conversion to from ComponentArray to NamedTuple (note, does not preserve numeric types of
# original NamedTuple)
function _namedtuple(x::ComponentVector)
    data = []
    idxmap = indexmap(getaxes(x)[1])
    for key in valkeys(x)
        idx = idxmap[getval(key)]
        # if idx isa AliasAxis
        #     val = idx.f
        # else
            val = getproperty(x, key) |> _namedtuple
        # end
        push!(data, getval(key) => val)
    end
    return (; data...)
end
_namedtuple(v::AbstractVector) = _namedtuple.(v)
_namedtuple(x) = x

Base.convert(::Type{NamedTuple}, x::ComponentVector) = _namedtuple(x)
Base.NamedTuple(x::ComponentVector) = _namedtuple(x)

## AbstractAxis conversion and promotion
Base.convert(::Type{Ax}, ax::AbstractAxis) where {Ax<:AbstractAxis} = ax #Ax()
# Base.convert(::Type{Axis}, ax::AbstractAxis) = ax #Ax()
Base.convert(::Type{Ax}, ax::AbstractAxis) where {Ax<:VarAxes} = convert.(typeof.(getaxes(Ax)), (ax,))

Base.promote_rule(AA::Type{<:Axis}, ::Type{NullAxis}) = AA
Base.promote_rule(AA::Type{<:Axis}, ::Type{FlatAxis}) = AA
Base.promote_rule(FA::Type{FlatAxis}, ::Type{NullAxis}) = FA
Base.promote_rule(FA::Type{FlatAxis}, ::Type{FlatAxis}) = FA
Base.promote_rule(NA::Type{NullAxis}, ::Type{NullAxis}) = NA
function Base.promote_rule(A1::Type{<:VarAxes}, A2::Type{<:VarAxes})
    # promote_type.(typeof.(getaxes(A1)), typeof.(getaxes(A2)))
    ax = typeof(first.(promote.(getaxes(A1), getaxes(A2))))
end

broadcast_promote_type(::Type{FlatAxis}, ax2::Type{<:Axis}) = ax2
broadcast_promote_type(ax1::Type{<:Axis}, ::Type{FlatAxis}) = ax1
broadcast_promote_type(ax1::Type{FlatAxis}, ::Type{NullAxis}) = ax1
broadcast_promote_type(::Type{NullAxis}, ax2::Type{FlatAxis}) = ax2
broadcast_promote_type(::Type{<:Axis}, ::Type{<:Axis}) = FlatAxis
broadcast_promote_type(ax1::Type{Ax}, ::Type{Ax}) where {Ax<:Axis} = ax1
function broadcast_promote_type(A1::Type{<:VarAxes}, A2::Type{<:VarAxes})
    ax = Tuple{broadcast_promote_typeof.(getaxes(A1), getaxes(A2))...}
end

broadcast_promote_typeof(ax1, ax2) = broadcast_promote_type(typeof(ax1), typeof(ax2))

broadcast_promote(ax1, ax2) = broadcast_promote_typeof(ax1, ax2)()


text:
using LinearAlgebra
import Simplices.simplexintersection
import Simplices.intersectingvertices

################################
# Exports
################################
export
AbstractSimplex,

# Size
npoints,
nvertices,
dimension,

# Properties of the simplex
orientation,
volume,
centroid,
radius,
issingular,

# Generating new points/simplices from existing ones
generate_interior_points,
generate_interior_simplex,
generate_exterior_point,
generate_exterior_points,
generate_exterior_simplex,
generate_intersecting_simplex,

# Intersections between simplices
intersect,
intersect_vertices,

# Plotting
connectvertices,
splitaxes

################################################################################
# Abstract simplex types
################################################################################
abstract type AbstractSimplex{D, T} end
abstract type AbstractImmutableSimplex{D, T} <: AbstractSimplex{D, T} end
abstract type AbstractMutableSimplex{D, T} <: AbstractSimplex{D, T} end

#######################
# Indexing
#######################
# Return the i-th point as column vector
Base.getindex(s::AbstractSimplex, i) = s.vertices[i]
Base.getindex(s::AbstractSimplex, i, j) = s.vertices[i][j]
Base.getindex(s::AbstractSimplex, i::Colon, j) = hcat(s.vertices[j]...,)
Base.getindex(s::AbstractSimplex, i::Colon, j::Colon) = hcat(s.vertices...,)

Base.length(s::AbstractSimplex) = length(s.vertices)
Base.size(s::AbstractSimplex) = length(s[1]), length(s)
Base.size(s::AbstractSimplex, i) = size(s)[i]
Base.IteratorSize(s::AbstractSimplex) = Base.HasLength()

Base.firstindex(s::AbstractSimplex) = 1
Base.lastindex(s::AbstractSimplex) = length(s)
Base.eachindex(s::AbstractSimplex) = Base.OneTo(length(s))
Base.iterate(s::AbstractSimplex, state = 1) = iterate(s.vertices, state)

###########################################################
# Vector{vertex} representation => Array representation
###########################################################
Base.Array(s::AbstractSimplex) = hcat(s.vertices...,)

#######################
# Sizes
#######################
dimension(s::AbstractSimplex) = length(s[1])
npoints(s::AbstractSimplex) = length(s)
nvertices(s::AbstractSimplex) = length(s)

export dimension, npoints, nvertices

#########################
# Orientation and volumes
#########################
# Orientation (convention: append rows of ones at top)
function orientation(s::AbstractSimplex)
    det([ones(1, dimension(s) + 1); s[:, :]])
end

""" Get the unscaled volume of the simplex. Divide by factorial(dim) to get the true volume."""
volume(s::AbstractSimplex) = abs(orientation(s))

#######################
# Centroid
#######################
function centroid(s::AbstractSimplex)
    D = dimension(s)
    # Results in a dim-by-1 matrix, but we just want a vector, so drop the last dimension
    centroid = dropdims(s[:, :] * (ones(D + 1, 1) / (D + 1)), dims = 2)
end

#######################
# Radius
#######################
function radius(s::AbstractSimplex)
    D = dimension(s)
    centroidmatrix = repeat(centroid(s), 1, D + 1)
    radius = sqrt(maximum(ones(1, D) * ((s[:, :] .- centroidmatrix) .^ 2)))
end

"""
    issingular(simplex::AbstractArray{T, 2}) where {T<:Number}

Determines if a simplex is singular by checking if any of its vertices are
identical.
"""
function issingular(simplex::AbstractArray{T, 2}) where {T<:Number}
    length(unique(simplex.vertices)) != length(simplex.vertices)
end

############################################
# Generating interior points/simplices
############################################

"""
    generate_interior_points(polytope, n::Int)

Generate `n` points in the interior of the polytope `p`, which
is represented by a dim-by-(# points) array.
"""
function generate_interior_points(polytope, n::Int)
    R = rand(size(polytope, 2), n)
    convex_coeffs = ((1 ./ sum(R, dims = 1)) .* R)
    polytope * convex_coeffs
end

"""
    generate_interior_point(simplex::T) where {T <: AbstractSimplex}

Generate a point that lies in the interior of `simplex`.
"""
function generate_interior_point(simplex::T) where {T <: AbstractSimplex}
    S = Array(simplex)
    dim = dimension(simplex)
    R = rand(dim + 1)
    convex_coeffs = ((1 ./ sum(R, dims = 1)) .* R)
    S * convex_coeffs
end

"""
    generate_interior_points(simplex::T, n::Int) where {T <: AbstractSimplex}

Generate `n` points that lies in the interior of `simplex`.
"""
function generate_interior_points(simplex::T, n::Int) where {T <: AbstractSimplex}
    S = Array(simplex)
    dim = dimension(simplex)

    # Random linear combination coefficients
    R = rand(dim + 1, n)

    # Normalise the coefficients so that they sum to one. We can then create the new point
    # as a convex linear combination of the vertices of the parent simplex.
    normalised_coeffs = (1 ./ sum(R, dims = 1)) .* R
    inside_pts = S * normalised_coeffs
    [inside_pts[:, i] for i = 1:n]
end

"""
    generate_interior_simplex(simplex::T) where {T <: AbstractSimplex}

Generate a new simplex contained in the interior of `simplex`.
"""
function generate_interior_simplex(simplex::T) where {T <: AbstractSimplex}
    # Create n_pts random points contained within the simplex
    n_pts = length(simplex)
    interior_pts = generate_interior_points(simplex, n_pts)

    # Return a simplex of the same type as the input.
    T([interior_pts[:, i] for i = 1:n_pts])
end

"""
    generate_exterior_point(simplex::T) where {T <: AbstractSimplex}

Generate a point in the exterior of `simplex`.
"""
function generate_exterior_point(simplex::T) where {T <: AbstractSimplex}
    parentsimplex = Array(simplex)
    dim = dimension(simplex)

    # Random linear combination coefficients
    R = rand(dim + 1)

    # Normalise the coefficients so that they sum to one. We can then create the new point
    # as a convex linear combination of the vertices of the parent simplex.
    normalised_coeffs = (1 ./ sum(R, dims = 1)) .* R
    normalised_coeffs[1] += 1.5

    parentsimplex * normalised_coeffs
end

"""
    generate_exterior_points(simplex::T, n::Int) where {T <: AbstractSimplex}

Generate `n` points in the exterior of `simplex`.
"""
function generate_exterior_points(simplex::T, n::Int) where {T <: AbstractSimplex}
    [generate_exterior_point(simplex) for i = 1:n]
end

"""
    generate_exterior_simplex(simplex::T) where {T <: AbstractSimplex}

Generate a new simplex lying in the exterior of `simplex`.
"""
function generate_exterior_simplex(simplex::T) where {T <: AbstractSimplex}
    T(generate_exterior_points(simplex, dimension(simplex) + 1))
end

#######################
# Subsampling simplices
#######################
include("simplex_subsampling.jl")

"""
    subsample(simplex::AbstractSimplex)

Subsample `simplex` with `n` points contained in its interior. If `sample_randomly = true`,
the points are randomly (uniformly) distributed within the simplex.
If `sample_randomly = false` (default), then a shape-preserving simplex splitting
routine is used to subdivide the simplex into subsimplices, and the centroids
of these simplices are taken as the sampling points. For the latter approach,
the number of sampling points is *at least `n`*, but might be somewhat higher
(the number of subsimplices generated by the splitting cannot be exactly `n`,
so we use a splitting  scheme that yields a minimum of `n` subsimplices).
"""
function subsample(simplex::AbstractSimplex; n::Int = 100, sample_randomly::Bool = false)
    dim = length(simplex[1])

    # Convex expansion coefficients to generate the centroids of the subsimplices
    # (if sample_randomly = false), or convex expansion coefficients to generate
    # random points.
    coeffs = subsample_coeffs(dim, n, sample_randomly)

    # Take the new points as as convex linear combinations of the vertices of
    # the simplex
    pts = simplex[:, :] * coeffs
    [SVector{dim, Float64}(pts[:, i]) for i = 1:size(pts, 2)]
end

export subsample

"""
    subsample(simplex::AbstractSimplex, coeffs::Array{Int, 2})

Subsample `simplex` using the precomputed convex coefficients in `coeffs`,
which contain instructions on how to combine the original vertices
to create the subsampled points. Will return as many points
as there are combinations/columns in `coeffs.`
"""
function subsample(simplex::AbstractSimplex, coeffs::Array{Float64, 2})
    dim = length(simplex[1])

    # Take the new points as as convex linear combinations of the vertices of
    # the simplex
    pts = simplex[:, :] * coeffs
    [SVector{dim, Float64}(pts[:, i]) for i = 1:size(pts, 2)]
end

"""
    subsample(simplex1::AbstractSimplex, simplex2::AbstractSimplex;
        n::Int = 100, sample_randomly::Bool = false)

Subsample `simplex` with `n` points contained in its interior.

If `sample_randomly = true`, the points are randomly (uniformly) distributed
within the simplex. If `sample_randomly = false` (default), then a
shape-preserving simplex splitting routine is used to subdivide the simplex
into subsimplices, and the centroids
of these simplices are taken as the sampling points. For the latter approach,
the number of sampling points is *at least `n`*, but might be somewhat higher
(the number of subsimplices generated by the splitting cannot be exactly `n`,
so we use a splitting  scheme that yields a minimum of `n` subsimplices).
"""
function subsample(simplex1::AbstractSimplex, simplex2::AbstractSimplex;
        n::Int = 100, sample_randomly::Bool = false)
    dim = length(simplex1[1])

    # Convex expansion coefficients to generate the centroids of the subsimplices
    # (if sample_randomly = false), or convex expansion coefficients to generate
    # random points.
    coeffs = subsample_coeffs(dim, n, sample_randomly)

    # Take the new points as as convex linear combinations of the vertices of
    # the simplex
    pts1 = simplex1[:, :] * coeffs
    pts2 = simplex2[:, :] * coeffs

    spts1 = [SVector{dim, Float64}(pts1[:, i]) for i = 1:size(pts1, 2)]
    spts2 = [SVector{dim, Float64}(pts2[:, i]) for i = 1:size(pts2, 2)]

    spts1, spts2
end

"""
    subsample(simplex::AbstractSimplex, coeffs::Array{Int, 2})

Subsample `simplex` using the precomputed convex coefficients in `coeffs`,
which contain instructions on how to combine the original vertices
to create the subsampled points. Will return as many points
as there are combinations/columns in `coeffs.`
"""
function subsample(simplex1::AbstractSimplex, simplex2::AbstractSimplex,
        coeffs::AbstractArray{Float64, 2})
    dim = length(simplex1[1])
    # Take the new points as as convex linear combinations of the vertices of
    # the simplex
    pts1 = simplex1[:, :] * coeffs
    pts2 = simplex2[:, :] * coeffs

    spts1 = [SVector{dim, Float64}(pts1[:, i]) for i = 1:size(pts1, 2)]
    spts2 = [SVector{dim, Float64}(pts2[:, i]) for i = 1:size(pts2, 2)]

    spts1, spts2
end

#######################
# Refine simplices
#######################
include("simplex_splitting.jl")

function sumvertices(vertices, k::Int, dim::Int)
    s = zeros(Float64, dim)
    for j = 1:k
        for i = 1:dim
            s[i] += vertices[j][i]
        end
    end
    return SVector{3, Float64}(s ./ k)
end

"""
    refine(simplex::Simplex; k::Int = 2)

Refine `simplex` with a shape-preserving sub-splitting of the simplex
into ``k^{dim}`` subsimplices.
"""
function refine(simplex::AbstractSimplex; k::Int = 2,
        splitting_rules = nothing)

    dim = length(simplex.vertices) - 1

    # Rules for forming the strictly new vertices of the subtriangulation
    if splitting_rules == nothing
        rules, subtriangulation = simplex_splitting_rules(k, dim)
    elseif splitting_rules isa Tuple{Array{Int64,2},Array{Int64,2}}
        # all is good
        rules, subtriangulation = splitting_rules
    else
        throw(ArgumentError("Either specify both rules and subtriangulation, or neither of them."))
    end

    # How many new vertices are created each split?
    nverts_persplit = size(rules, 1)

    # The original vertices and strictly new sub-vertices formed by the splitting.
    all_vertices = Vector{SVector{dim, Float64}}(undef, nverts_persplit)

    # Generate the strictly new vertices for each subsimplex
    for j = 1:nverts_persplit
        # Pick the corresponding original vertices with indices contained in rules[j, :]
        # and compute their combination, which will be the j-th new vertex.
        all_vertices[j] = sumvertices(simplex.vertices[rules[j, :]], k, dim)
    end

    # Generate the new simplices
    subsimplices = Vector{Simplex{3, Float64}}(undef, k^dim)
    for i = 1:k^dim
        subsimplices[i] = Simplex(all_vertices[subtriangulation[:, i]])
    end

    return subsimplices
end

"""
    refine(simplex1::Simplex, simplex2::Simplex; k::Int = 2,
        splitting_rules::Union{Nothing, Tuple{Array{Float64, 2}, Array{Int, 2}}} = nothing)

Refine `simplex1` and `simplex2` with a shape-preserving sub-splitting of
the simplices into ``k^{dim}`` subsimplices.

This method may be used, for example, to keep track of which simplices
correspond to each other when splitting a simplex and its image simplex
(under the forward linear map of the vertices of the original simplex).

`splitting_rules` may be generated beforehand to ensure that the same
coefficients are used if many pairs of simplices are refined using
the `simplex_splitting_rules` function.
"""
function refine(simplex1::AbstractSimplex, simplex2::AbstractSimplex;
        k::Int = 2,
        splitting_rules = nothing)

    dim = length(simplex1.vertices) - 1

    # Rules for forming the strictly new vertices of the subtriangulation
    if splitting_rules == nothing
        rules, subtriangulation = simplex_splitting_rules(k, dim)
    elseif splitting_rules isa Tuple{Array{Int64,2},Array{Int64,2}}
        # all is good
        rules, subtriangulation = splitting_rules
    else
        throw(ArgumentError("Either specify both rules and subtriangulation, or neither of them."))
    end

    # How many new vertices are created each split?
    nverts_persplit = size(rules, 1)

    # The original vertices and strictly new sub-vertices formed by the splitting.
    all_vertices1 = Vector{SVector{dim, Float64}}(undef, nverts_persplit)
    all_vertices2 = Vector{SVector{dim, Float64}}(undef, nverts_persplit)

    # Generate the strictly new vertices for each subsimplex
    for j = 1:nverts_persplit
        # Pick the corresponding original vertices with indices contained in rules[j, :]
        # and compute their combination, which will be the j-th new vertex.
        all_vertices1[j] = sumvertices(simplex1.vertices[rules[j, :]], k, dim)
        all_vertices2[j] = sumvertices(simplex2.vertices[rules[j, :]], k, dim)
    end

    # Generate the new simplices
    subsimplices1 = Vector{Simplex{3, Float64}}(undef, k^dim)
    subsimplices2 = Vector{Simplex{3, Float64}}(undef, k^dim)

    for i = 1:k^dim
        subsimplices1[i] = Simplex(all_vertices1[subtriangulation[:, i]])
        subsimplices2[i] = Simplex(all_vertices2[subtriangulation[:, i]])
    end

    return subsimplices1, subsimplices2
end

export refine

#######################
# Intersections
#######################
"""
    intersect(s1::AbstractSimplex, s2::AbstractSimplex)

Compute the volume intersection between two simplices.
"""
function Base.intersect(s1::AbstractSimplex, s2::AbstractSimplex)
    simplexintersection(s1[:, :], s2[:, :])
end

"""
    ∩(s1::AbstractSimplex, s2::AbstractSimplex)

Compute the volume intersection between two simplices.
"""
∩(s1::AbstractSimplex, s2::AbstractSimplex) = intersect(s1, s2)

"""
    intersect_vertices(s1::AbstractSimplex, s2::AbstractSimplex)

Return the vertices forming the convex hull of the volume intersection
between two simplices.
"""
function intersect_vertices(s1::AbstractSimplex, s2::AbstractSimplex)
    intersectingvertices(s1[:, :], s2[:, :])
end

"""
    generate_intersecting_simplex(simplex::T) where {T <: AbstractSimplex}

Generate a new simplex that intersects with `simplex` in a nontrivial way.
"""
function generate_intersecting_simplex(simplex::T) where {T <: AbstractSimplex}
    dim = dimension(simplex)

    n_inside = rand(1:dim)
    n_outside = dim + 1 - n_inside

    verts_inside = generate_interior_points(simplex, n_inside)
    verts_outside = generate_exterior_points(simplex, n_outside)

    T([verts_inside; verts_outside])
end

#######################
# Pretty printing
#######################
function summarise(s::AbstractSimplex)
    _type = typeof(s)
    n_vertices = npoints(s)
    D = dimension(s)

    matrixrepresentation = sprint(io -> show(IOContext(io, :limit=>true),
                MIME"text/plain"(), s[:, :]))
    matrixrepresentation = join(split(matrixrepresentation, '\n')[2:end], '\n')

    summary = join(["$D-dimensional $_type with $n_vertices vertices", "\n",
                    matrixrepresentation])
end

Base.show(io::IO, s::AbstractSimplex) = println(io, summarise(s))

################################
# Helper functions for plotting
################################
function connectvertices(s::T) where {T <: AbstractSimplex}
    # Convert to array
    s_arr = Array(s)
    hcat(s_arr, hcat([s_arr[:, i] for i in [1, 3, 1, 4, 2]]...,))
end

splitaxes(x) = ([x[k, :] for k = 1:size(x, 1)]...,)

ext: jl


text:
# All the IO related APIs.

export readdata, readlogdata, readtecdata, showhead

const tag = 4 # Fortran record tag

searchdir(path, key) = filter(x->occursin(key, x), readdir(path))

function Base.show(io::IO, data::Data)
   showhead(io, data)
   if data.list.bytes ≥ 1e9
      println(io, "filesize = $(data.list.bytes/1e9) GB")
   elseif data.list.bytes ≥ 1e6
      println(io, "filesize = $(data.list.bytes/1e6) MB")
   elseif data.list.bytes ≥ 1e3
      println(io, "filesize = $(data.list.bytes/1e3) KB")
   else
      println(io, "filesize = $(data.list.bytes) bytes")
   end
   println(io, "snapshots = $(data.list.npictinfiles)")
end

"""
	readdata(filenameIn; dir=".", npict=1, verbose=false)

Read data from BATSRUS output files. Stores the `npict` snapshot from an ascii or binary
data file into the arrays of coordinates `x` and data `w`.
Filename can be provided with wildcards.

# Examples
```
data = readdata("1d_raw*")
```
"""
function readdata(filenameIn::AbstractString; dir=".", npict=1, verbose=false)

   # Check the existence of files
   filename = searchdir(dir, Regex(replace(filenameIn, s"*" => s".*")))
   if isempty(filename)
      throw(ArgumentError("No matching filename was found for $(filenameIn)"))
   elseif length(filename) > 1
      throw(ArgumentError("Ambiguous filename $(filenameIn)!"))
   end

   filelist, fileID, pictsize = getfiletype(filename[1], dir)

   verbose &&
      @info "filename=$(filelist.name)\n"*"npict=$(filelist.npictinfiles)"

   if filelist.npictinfiles - npict < 0
      throw(ArgumentError("select snapshot $npict out of range $(filelist.npictinfiles)!"))
   end
   seekstart(fileID) # Rewind to start

   ## Read data from files
   # Skip npict-1 snapshots (because we only want the npict-th snapshot)
   skip(fileID, pictsize*(npict-1))

   filehead = getfilehead(fileID, filelist.type)

   # Read data
   fileType = lowercase(filelist.type)
   if fileType == "ascii"
      x, w = allocateBuffer(filehead, Float64) # why Float64?
      getascii!(x, w, fileID, filehead)
   else
      skip(fileID, tag) # skip record start tag.
      T = fileType == "real4" ? Float32 : Float64
      x, w = allocateBuffer(filehead, T)
      getbinary!(x, w, fileID, filehead, T)
   end

   close(fileID)

   #setunits(filehead,"")

	data = Data(filehead, x, w, filelist)

   verbose && @info "Finished reading $(filelist.name)"

   data
end

"Read information from log file."
function readlogdata(filename::AbstractString)

   f = open(filename, "r")
   nLine = countlines(f) - 2
   seekstart(f)
   headline  = readline(f)
   variables = split(readline(f))
   ndim      = 1
   it        = 0
   t         = 0.0
   gencoord  = false
   nx        = 1
   nw        = length(variables)

   data = zeros(nw,nLine)
   for i = 1:nLine
      line = split(readline(f))
      data[:,i] = parse.(Float64, line)
   end

   close(f)

   head = (ndim=ndim, headline=headline, it=it, time=t, gencoord=gencoord,
      nw=nw, nx=nx, variables=variables)

   head, data
end

"""
	readtecdata(filename; verbose=false)

Return header, data and connectivity from BATSRUS Tecplot outputs. Both 2D and
3D binary and ASCII formats are supported.
# Examples
```
filename = "3d_ascii.dat"
head, data, connectivity = readtecdata(filename)
```
"""
function readtecdata(filename::AbstractString; verbose=false)

   f = open(filename)

   nDim  = 3
   nNode = Int32(0)
   nCell = Int32(0)
   ET = ""

   # Read Tecplot header
   ln = readline(f) |> strip
   if startswith(ln, "TITLE")
      title = match(r"\"(.*?)\"", split(ln,'=', keepempty=false)[2])[1]
   else
      @warn "No title provided."
   end
   ln = readline(f) |> strip
   if startswith(ln, "VARIABLES")
      # Read until another keyword appears
      varline = split(ln,'=')[2]

      ln = readline(f) |> strip
      while !startswith(ln, "ZONE")
         varline *= strip(ln)
         ln = readline(f) |> strip
      end
      VARS = split(varline, '\"')
      deleteat!(VARS,findall(x -> x in (" ",", ") || isempty(x), VARS))
   else
      @warn "No variable names provided."
   end

   while !startswith(ln, "AUXDATA")
      if !startswith(ln, "ZONE") # ZONE allows multiple \n
         zoneline = split(ln, ", ", keepempty=false)
      else # if the ZONE line has nothing, this won't work!
         zoneline = split(ln[6:end], ", ", keepempty=false)
         replace(zoneline[1], '"'=>"") # Remove the quotes in T
      end
      for zline in zoneline
         name, value = split(zline,'=', keepempty=false)
         name = uppercase(name)
         if name == "T" # ZONE title
            T = value
         elseif name in ("NODES","N")
            nNode = parse(Int32, value)
         elseif name in ("ELEMENTS","E")
            nCell = parse(Int32, value)
         elseif name in ("ET","ZONETYPE")
            if uppercase(value) in ("BRICK","FEBRICK")
               nDim = 3
            elseif uppercase(value) in ("QUADRILATERAL", "FEQUADRILATERAL")
               nDim = 2
            end
            ET = uppercase(value)
         end
      end
      ln = readline(f) |> strip
   end

   auxdataname = String[]
   auxdata = Union{Int32, String}[]
   pt0 = position(f)

   while startswith(ln, "AUXDATA") || startswith(ln,"DT")
      name, value = split(ln,'"', keepempty=false)
      name = string(name[9:end-1])
      str = string(strip(value))
      if name in ("ITER","NPROC")
         str = parse(Int32, value)
      elseif name == "TIMESIM"
         sec = split(str,"=")
         str = string(strip(sec[2]))
      end
      push!(auxdataname, name)
      push!(auxdata, str)

      pt0 = position(f)
      ln = readline(f) |> strip
   end

   seek(f, pt0)

   data = Array{Float32,2}(undef, length(VARS), nNode)

   if nDim == 3
	   connectivity = Array{Int32,2}(undef,8,nCell)
   elseif nDim == 2
	   connectivity = Array{Int32,2}(undef,4,nCell)
   end

   IsBinary = false
   try
      parse.(Float32, split(readline(f)))
   catch
      IsBinary = true
      verbose && @info "reading binary file"
   end

   seek(f, pt0)

   if IsBinary
	   @inbounds for i = 1:nNode
		   read!(f, @view data[:,i])
	   end
	   @inbounds for i = 1:nCell
         read!(f, @view connectivity[:,i])
      end
   else
      @inbounds for i = 1:nNode
         x = readline(f)
         data[:,i] .= parse.(Float32, split(x))
      end
	   @inbounds for i = 1:nCell
		   x = readline(f)
		   connectivity[:,i] .= parse.(Int32, split(x))
	   end
   end

   close(f)

   head = (variables=VARS, nNode=nNode, nCell=nCell, nDim=nDim, ET=ET,
		title=title, auxdataname=auxdataname, auxdata=auxdata)

   head, data, connectivity
end

"Obtain file type."
function getfiletype(filename, dir)
   file = joinpath(dir, filename)
   fileID = open(file, "r")
   bytes = filesize(file)
   type  = ""

   # Check the appendix of file names
   # Gabor uses a trick: the first 4 bytes decides the file type
   if occursin(r"^.*\.(log)$", filename)
      type = "log"
      npictinfiles = 1
   elseif occursin(r"^.*\.(dat)$", filename)
      # Tecplot ascii format
      type = "dat"
      npictinfiles = 1
   else
      # Obtain filetype based on the length info in the first 4 bytes
      lenhead = read(fileID, Int32)

      if lenhead!=79 && lenhead!=500
         type = "ascii"
      else
         # The length of the 2nd line decides between real4 & real8
         # since it contains the time; which is real*8 | real*4
         skip(fileID, lenhead+tag)
         len = read(fileID, Int32)
         if len == 20
            type = "real4"
         elseif len == 24
            type = "binary"
         else
            throw(ArgumentError(
               "Error in getfiletype: incorrect formatted file: $filename"))
         end

         if lenhead == 500
            type = uppercase(type)
         end
      end
      # Obtain file size & number of snapshots
      seekstart(fileID)
      pictsize = getfilesize(fileID, type)
      npictinfiles = bytes ÷ pictsize
   end

   filelist = FileList(filename, type, dir, bytes, npictinfiles)

   filelist, fileID, pictsize
end

"""
	getfilehead(fileID, type) -> NameTuple

Obtain the header information from BATSRUS output file of `type` linked to `fileID`.
# Input arguments
- `fileID::IOStream`: file identifier.
- `type::String`: file type in ["ascii", "real4", "binary", "log"].
"""
function getfilehead(fileID::IOStream, type)

   ftype = string(lowercase(type))

   lenstr = ftype == type ? 79 : 500

   ## Read header
   pointer0 = position(fileID)

   if ftype == "ascii"
      headline = readline(fileID)
      line = readline(fileID)
      line = split(line)
      it = parse(Int,line[1])
      t = parse(Float64,line[2])
      ndim = parse(Int8,line[3])
      neqpar = parse(Int32,line[4])
      nw = parse(Int8,line[5])
      gencoord = ndim < 0
      ndim = abs(ndim)
      nx = parse.(Int64, split(readline(fileID)))
      if neqpar > 0
         eqpar = parse.(Float64, split(readline(fileID)))
      end
      varname = readline(fileID)
   elseif ftype ∈ ["real4", "binary"]
      skip(fileID, tag) # skip record start tag.
      headline = rstrip(String(read(fileID, lenstr)))
      skip(fileID, 2*tag) # skip record end/start tags.
      it = read(fileID, Int32)
      t = read(fileID, Float32)
      ndim = read(fileID, Int32)
      gencoord = (ndim < 0)
      ndim = abs(ndim)
      neqpar = read(fileID, Int32)
      nw = read(fileID, Int32)
      skip(fileID, 2*tag) # skip record end/start tags.
      nx = zeros(Int32, ndim)
      read!(fileID, nx)
      skip(fileID, 2*tag) # skip record end/start tags.
      if neqpar > 0
         eqpar = zeros(Float32, neqpar)
         read!(fileID, eqpar)
         skip(fileID, 2*tag) # skip record end/start tags.
      end
      varname = String(read(fileID, lenstr))
      skip(fileID, tag) # skip record end tag.
   end

   # Header length
   pointer1 = position(fileID)
   headlen = pointer1 - pointer0

   # Set variables array
   variables = split(varname) # returns a string array

	# Produce a wnames from the last file
   wnames = variables[ndim+1:ndim+nw]

   head = (ndim=ndim, headline=headline, it=it, time=t, gencoord=gencoord,
		neqpar=neqpar, nw=nw, nx=nx, eqpar=eqpar, variables=variables,
		wnames=wnames)
end

function skipline(s::IO)
   while !eof(s)
       c = read(s, Char)
       c == '\n' && break
   end
   nothing
end

"Return the size in bytes for one snapshot."
function getfilesize(fileID::IOStream, type::String)

   ftype = string(lowercase(type))

   if ftype == type lenstr = 79 else lenstr = 500 end

   # Read header
   pointer0 = position(fileID)

   if ftype == "ascii"
      skipline(fileID)
      line = readline(fileID)
      line = split(line)
      ndim = parse(Int32,line[3])
      neqpar = parse(Int32,line[4])
      nw = parse(Int8,line[5])
      gencoord = ndim < 0
      ndim = abs(ndim)
      nx = parse.(Int64, split(readline(fileID)))
      neqpar > 0 && skipline(fileID)
      skipline(fileID)
   elseif ftype ∈ ("real4", "binary")
      skip(fileID, tag)
      read(fileID, lenstr)
      skip(fileID, 2*tag)
      read(fileID, Int32)
      read(fileID, Float32)
      ndim = abs(read(fileID, Int32))
      tmp = read(fileID, Int32)
      nw = read(fileID, Int32)
      skip(fileID, 2*tag)
      nx = zeros(Int32, ndim)
      read!(fileID, nx)
      skip(fileID, 2*tag)
      if tmp > 0
         tmp2 = zeros(Float32, tmp)
         read!(fileID, tmp2)
         skip(fileID, 2*tag) # skip record end/start tags.
      end
      read(fileID, lenstr)
      skip(fileID, tag)
   end

   # Header length
   pointer1 = position(fileID)
   headlen = pointer1 - pointer0

   # Calculate the snapshot size = header + data + recordmarks
   nxs = prod(nx)
   if ftype == "log"
      pictsize = 1
   elseif ftype == "ascii"
      pictsize = headlen + (18*(ndim+nw)+1)*nxs
   elseif ftype == "real4"
      pictsize = headlen + 8*(1+nw) + 4*(ndim+nw)*nxs
   elseif ftype == "binary"
      pictsize = headlen + 8*(1+nw) + 8*(ndim+nw)*nxs
   end

   pictsize
end

"Create buffer for x and w."
function allocateBuffer(filehead::NamedTuple, T::DataType)

   if filehead.ndim == 1
      n1 = filehead.nx[1]
      x  = Array{T,2}(undef,n1,filehead.ndim)
      w  = Array{T,2}(undef,n1,filehead.nw)
   elseif filehead.ndim == 2
      n1, n2 = filehead.nx
      x  = Array{T,3}(undef,n1,n2,filehead.ndim)
      w  = Array{T,3}(undef,n1,n2,filehead.nw)
   elseif filehead.ndim == 3
      n1, n2, n3 = filehead.nx
      x  = Array{T,4}(undef,n1,n2,n3,filehead.ndim)
      w  = Array{T,4}(undef,n1,n2,n3,filehead.nw)
   end

   x, w
end

"Read ascii format data."
function getascii!(x, w, fileID::IOStream, filehead::NamedTuple)

   ndim, nx = filehead.ndim, filehead.nx

   # Read coordinates & values row by row
   if ndim == 1 # 1D
      for ix = 1:nx[1]
         temp = parse.(Float64, split(readline(fileID)))
         x[ix,:] .= temp[1]
         w[ix,:] .= temp[2:end]
      end
   elseif ndim == 2 # 2D
      for j = 1:nx[2], i = 1:nx[1]
         temp = parse.(Float64, split(readline(fileID)))
         x[i,j,:] .= temp[1:2]
         w[i,j,:] .= temp[3:end]
      end
   elseif ndim == 3 # 3D
      for k = 1:nx[3], j = 1:nx[2], i = 1:nx[1]
         temp = parse.(Float64, split(readline(fileID)))
         x[i,j,k,:] .= temp[1:3]
         w[i,j,k,:] .= temp[4:end]
      end
   end

   nothing
end

"Read binary format data."
function getbinary!(x, w, fileID::IOStream, filehead::NamedTuple, T::DataType)

   ndim = filehead.ndim
   nw   = filehead.nw

   # Read coordinates & values
   if ndim == 1 # 1D
      read!(fileID, x)
      skip(fileID, 2*tag) # skip record end/start tags.
      for iw = 1:nw
         read!(fileID, @view w[:,iw])
         skip(fileID, 2*tag) # skip record end/start tags.
      end
   elseif ndim == 2 # 2D
      read!(fileID, x)
      skip(fileID, 2*tag) # skip record end/start tags.
      for iw = 1:nw
         read!(fileID, @view w[:,:,iw])
         skip(fileID, 2*tag) # skip record end/start tags.
      end
   elseif ndim == 3 # 3D
      read!(fileID, x)
      skip(fileID, 2*tag) # skip record end/start tags.
      for iw = 1:nw
         read!(fileID, @view w[:,:,:,iw])
         skip(fileID, 2*tag) # skip record end/start tags.
      end
   end

   nothing
end

"""
	setunits(filehead, type; distance=1.0, mp=1.0, me=1.0)

Set the units for the output files.
If type is given as "SI", "CGS", "NORMALIZED", "PIC", "PLANETARY", "SOLAR", set
`typeunit = type`, otherwise try to guess from the fileheader. Based on `typeunit` set units
for distance [xSI], time [tSI], density [rhoSI], pressure [pSI], magnetic field [bSI] and
current density [jSI] in SI units. Distance unit [rplanet | rstar], ion and electron mass in
[amu] can be set with optional `distance`, `mp` and `me`.

Also calculate convenient constants ti0, cs0 ... for typical formulas.
This function is currently not used anywhere!
"""
function setunits(filehead, type; distance=1.0, mp=1.0, me=1.0)

   ndim      = filehead.ndim
   headline  = filehead.headline
   neqpar    = filehead.neqpar
   nw        = filehead.nw
   eqpar     = filehead.eqpar
   variables = filehead.variables

   mu0SI = 4π*1e-7      # H/m
   cSI   = 2.9978e8       # speed of light, [m/s]
   mpSI  = 1.6726e-27     # kg
   eSI   = 1.602e-19      # elementary charge, [C]
   AuSI  = 149597870700   # m
   RsSI  = 6.957e8        # m
   Mi    = 1.0            # Ion mass, [amu]
   Me    = 1.0/1836.15    # Electron mass, [amu]
   gamma = 5/3            # Adiabatic index for first fluid
   gammae= 5/3            # Adiabatic index for electrons
   kbSI  = 1.38064852e-23 # Boltzmann constant, [m2 kg s-2 K-1]
   e0SI  = 8.8542e-12     # [F/m]

   # This part is used to guess the units.
   # To be honest; I don`t understand the logic here. For example
   # nPa & m/s may appear in the same headline?
   if type !== ""
      typeunit = uppercase(type)
   elseif occursin("PIC", headline)
      typeunit = "PIC"
   elseif occursin(" AU ", headline)
      typeunit = "OUTERHELIO"
   elseif occursin(r"(kg/m3)|(m/s)", headline)
      typeunit = "SI"
   elseif occursin(r"(nPa)|( nT )", headline)
      typeunit = "PLANETARY"
   elseif occursin(r"(dyne)|( G)", headline)
      typeunit = "SOLAR"
   else
      typeunit = "NORMALIZED"
   end

   if typeunit == "SI"
      xSI   = 1.0             # m
      tSI   = 1.0             # s
      rhoSI = 1.0             # kg/m^3
      uSI   = 1.0             # m/s
      pSI   = 1.0             # Pa
      bSI   = 1.0             # T
      jSI   = 1.0             # A/m^2
   elseif typeunit == "CGS"
      xSI   = 0.01            # cm
      tSI   = 1.0             # s
      rhoSI = 1000.0          # g/cm^3
      uSI   = 0.01            # cm/s
      pSI   = 0.1             # dyne/cm^2
      bSI   = 1.0e-4          # G
      jSI   = 10*cSI          # Fr/s/cm^2
   elseif typeunit == "PIC"
      # Normalized PIC units
      xSI   = 1.0             # cm
      tSI   = 1.0             # s
      rhoSI = 1.0             # g/cm^3
      uSI   = 1.0             # cm/s
      pSI   = 1.0             # dyne/cm^2
      bSI   = 1.0             # G
      jSI   = 1.0             # Fr/s/cm^2
      c0    = 1.0             # speed of light always 1 for iPIC3D
   elseif typeunit == "NORMALIZED"
      xSI   = 1.0             # distance unit in SI
      tSI   = 1.0             # time unit in SI
      rhoSI = 1.0             # density unit in SI
      uSI   = 1.0             # velocity unit in SI
      pSI   = 1.0             # pressure unit in SI
      bSI   = √mu0SI     # magnetic unit in SI
      jSI   = 1/√(mu0SI)   # current unit in SI
      c0    = 1.0             # speed of light (for Boris correction)
   elseif typeunit == "PLANETARY"
      xSI   = 6378000         # Earth radius [default planet]
      tSI   = 1.0             # s
      rhoSI = mpSI*1e6        # mp/cm^3
      uSI   = 1e3             # km/s
      pSI   = 1e-9            # nPa
      bSI   = 1e-9            # nT
      jSI   = 1e-6            # muA/m^2
      c0    = cSI/uSI         # speed of light in velocity units
   elseif typeunit == "OUTERHELIO"
      xSI   = AuSI            # AU
      tSI   = 1.0             # s
      rhoSI = mpSI*1e6        # mp/cm^3
      uSI   = 1e3             # km/s
      pSI   = 1e-1            # dyne/cm^2
      bSI   = 1e-9            # nT
      jSI   = 1e-6            # muA/m^2
      c0    = cSI/uSI         # speed of light in velocity units
   elseif typeunit == "SOLAR"
      xSI   = RsSI            # radius of the Sun
      tSI   = 1.0             # s
      rhoSI = 1e3             # g/cm^3
      uSI   = 1e3             # km/s
      pSI   = 1e-1            # dyne/cm^2
      bSI   = 1e-4            # G
      jSI   = 1e-6            # muA/m^2
      c0    = cSI/uSI         # speed of light in velocity units
   else
      throw(ArgumentError("invalid typeunit=$(typeunit)"))
   end

   # Overwrite values if given by eqpar
   for ieqpar = 1:neqpar
      var = variables[ndim+nw+ieqpar]
      if var == "xSI"
         xSI   = eqpar[ieqpar]
      elseif var == "tSI"
         tSI   = eqpar[ieqpar]
      elseif var == "uSI"
         uSI   = eqpar[ieqpar]
      elseif var == "rhoSI"
         rhoSI = eqpar[ieqpar]
      elseif var == "mi"
         mi    = eqpar[ieqpar]
      elseif var == "m1"
         m1    = eqpar[ieqpar]
      elseif var == "me"
         me    = eqpar[ieqpar]
      elseif var == "qi"
         qi    = eqpar[ieqpar]
      elseif var == "q1"
         q1    = eqpar[ieqpar]
      elseif var == "qe"
         qe    = eqpar[ieqpar]
      elseif var == "g"
         gamma = eqpar[ieqpar]
      elseif var == "g1"
         gamma = eqpar[ieqpar]
      elseif var == "ge"
         ge    = eqpar[ieqpar]
      elseif var == "c"
         c     = eqpar[ieqpar]
      elseif var == "clight"
         clight= eqpar[ieqpar]
      elseif var == "r"
         r     = eqpar[ieqpar]
      elseif var == "rbody"
         rbody = eqpar[ieqpar]
      end
   end

   # Overwrite distance unit if given as an argument
   if !isempty(distance) xSI = distance end

   # Overwrite ion & electron masses if given as an argument
   if !isempty(mp)      Mi = mp end
   if !isempty(me) Me = me end

   # Calculate convenient conversion factors
   if typeunit == "NORMALIZED"
      ti0  = 1.0/Mi            # T      = p/rho*Mi           = ti0*p/rho
      cs0  = 1.0               # cs     = sqrt(gamma*p/rho)  = sqrt(gs*p/rho)
      mu0A = 1.0               # vA     = sqrt(b/rho)        = sqrt(bb/mu0A/rho)
      mu0  = 1.0               # beta   = p/(bb/2)           = p/(bb/(2*mu0))
      uH0  = Mi                # uH     = j/rho*Mi           = uH0*j/rho
      op0  = 1.0/Mi            # omegap = sqrt(rho)/Mi       = op0*sqrt(rho)
      oc0  = 1.0/Mi            # omegac = b/Mi               = oc0*b
      rg0  = √Mi               # rg = sqrt(p/rho)/b*sqrt(Mi) = rg0*sqrt(p/rho)/b
      di0  = c0*Mi             # di = c0/sqrt(rho)*Mi        = di0/sqrt(rho)
      ld0  = Mi                # ld = sqrt(p)/(rho*c0)*Mi    = ld0*sqrt(p)/rho
   elseif typeunit == "PIC"
      ti0  = 1.0/Mi            # T      = p/rho*Mi           = ti0*p/rho
      cs0  = 1.0               # cs     = sqrt(gamma*p/rho)  = sqrt(gs*p/rho)
      mu0A = 4*pi              # vA     = sqrt(b/(4*!pi*rho))= sqrt(bb/mu0A/rho)
      mu0  = 4*pi              # beta   = p/(bb/(8*!pi))     = p/(bb/(2*mu0))
      uH0  = Mi                # uH     = j/rho*Mi           = uH0*j/rho
      op0  = √(4π)/Mi          # omegap = sqrt(4*!pi*rho)/Mi = op0*sqrt(rho)
      oc0  = 1.0/Mi            # omegac = b/Mi               = oc0*b
      rg0  = √Mi               # rg = sqrt(p/rho)/b*sqrt(Mi) = rg0*sqrt(p/rho)/b
      di0  = 1.0/√(4π)         # di = 1/sqrt(4*!pi*rho)*Mi   = di0/sqrt(rho)
      ld0  = 1.0/√(4π)         # ld = sqrt(p/(4*!pi))/rho*Mi = ld0*sqrt(p)/rho
   else
      qom  = eSI/(Mi*mpSI); moq = 1/qom
      ti0  = mpSI/kbSI*pSI/rhoSI*Mi       # T[K]=p/(nk) = ti0*p/rho
      cs0  = pSI/rhoSI/uSI^2              # cs          = sqrt(gs*p/rho)
      mu0A = uSI^2*mu0SI*rhoSI*bSI^(-2)   # vA          = sqrt(bb/(mu0A*rho))
      mu0  = mu0SI*pSI*bSI^(-2)           # beta        = p/(bb/(2*mu0))
      uH0  = moq*jSI/rhoSI/uSI            # uH=j/(ne)   = uH0*j/rho
      op0  = qom*√(rhoSI/e0SI)*tSI        # omegap      = op0*sqrt(rho)
      oc0  = qom*bSI*tSI                  # omegac      = oc0*b
      rg0  = moq*√(pSI/rhoSI)/bSI/xSI/√(Mi)    # rg     = rg0*sqrt(p/rho)/b
      di0  = cSI/(op0/tSI)/xSI                 # di=c/omegap = di0/sqrt(rho)
      ld0  = moq*√(pSI)/rhoSI/xSI              # ld          = ld0*sqrt(p)/rho
   end
   return true
end

"""
	showhead(file, filehead)

Displaying file header information.
"""
function showhead(file::FileList, head, io=stdout)
   println(io, "filename : $(file.name)")
   println(io, "filetype : $(file.type)")
   println(io, "headline : $(head.headline)")
   println(io, "iteration: $(head.it)")
   println(io, "time     : $(head.time)")
   println(io, "gencoords: $(head.gencoord)")
   println(io, "ndim     : $(head.ndim)")
   println(io, "neqpar   : $(head.neqpar)")
   println(io, "nw       : $(head.nw)")
   println(io, "nx       : $(head.nx)")

   if head.neqpar > 0
      println(io, "parameters  = $(head.eqpar)")
      println(io, "coord names = $(head.variables[1:head.ndim])")
      println(io, "var   names = $(head.variables[head.ndim+1:head.ndim+head.nw])")
      println(io, "param names = $(head.variables[head.ndim+head.nw+1:end])")
   end
end

"""
	showhead(data)
   showhead(io, data)

Display file information of `data`.
"""
showhead(data::Data) = showhead(data.list, data.head)
showhead(io, data::Data) = showhead(data.list, data.head, io)


text:
planck(ν, T) =
    u"W*m^-2*cm*sr^-1"(c_1 * ν^3 / (exp(c_2 * ν / T) - 1) * 1u"sr^-1")

radiance_to_bt(ν, L) = u"K"(c_2 * ν / log(c_1 * ν^3 / L + 1))
radiance_to_bt(ν) = L -> radiance_to_bt(ν, L)

get_molecular_weight(gas) = LBLRTM_GAS_MOLECULAR_WEIGHTS[LBLRTM_GAS_INDEX[gas]]


text:
## We'll learn about control flow!

## Simple stuff
x = Inf # we can use Infinity as a number !
        # this makes a Float64
typeof(x)
##
if x < 1
    println("x is smaller than 1")
end
##
if x < 1
    println("x is smaller than 1")
else
    println("x is larger than or equal to 1")
end

##
x = -5
if x < 1
    println("x is smaller than 1")
elseif x < 0
    println("x is smaller than 0")
else
    println("x is larger than or equal to 1")
end

##
if x < 0
    println("x is smaller than 0")
elseif x < 1
    println("x is smaller than 1")
else
    println("x is larger than or equal to 1")
end

##
x = [1.0, 2.0]
if all(x >= 0) && all( x <= 5)
    println("all elements between 0 and 5")
end

##
x = [1.0, 2.0]
if all(x .>= 0) && all(x .<= 5)
    println("all elements between 0 and 5")
end
##
x = 0
y = 5
if x+y > 10 || x < 0
    y = 0
else
    y = 10
end

## That block returned something,  you can "save"
# the result of a comparison
z = if x+y > 10 || x < 0
    y = 0
else
    y = 10
end

## That's it for "if, elseif, else"

## On to for floops

for i=1:10 # for i = 1 to 10, starting a 1 , ending at 10
    @show i
    println("end loop iteration")
end

##
v = -5:5
for i=v # for i = 1 to 10, starting a 1 , ending at 10
    @show i
    println("end loop iteration")
end

##
for i in v
    @show i
    println("end loop iteration")
end
##
v = [1.0,Inf,-Inf,π]
for i in v
    @show i
    println("end loop iteration")
end
##
v = [1 2 3 4]
for i in v
    @show i
    println("end loop iteration")
end
## Julia iterates over matrices by columns
@show [v' v'.+1]
for i=[v' v'.+1]
    @show i
end
##
x = 1:5
x[3] = -3
##
x = collect(1:5)
x[3] = -3

##
x = copy(1:5)
x[3] = -3
## We want to turn this into a real array
x = collect(1:5)
x[3] = -3

##
@show A = [x' x'.+1]
for i in A
    @show i
end

##
##
@show A = [x x.+1]
for i in A
    @show i
end

##While loops
i = 0
while i < 10
    global i
    if i == 5
        break # this works in for loops too!
    end
    if i==4
        i = i+5
    end
    @show i
    i = i+1
end
@show i

##
function example_while()
    i = 0
    while i < 10
        if i == 5
            break # this works in for loops too!
        end
        if i==4
            i = i+5
        end
        @show i
        i = i+1
    end
    @show i
end
example_while()

## The place I usually run into that issue.
# Let's say I want to compute the number of random variables
# that are bigger than pi.
numsuccess = 0
numtrials = 10000
for i=1:numtrials
    numsuccess += randn() > pi
end

# This case hits that scoping issue :(
# In IJulia/Jupyter, they include a patch
# And working out something else is an active area of
# development for Julia, so this may change in the
# future.

##
numsuccess = 0
numtrials = 10000
for i=1:numtrials
    global numsuccess
    numsuccess += randn() > pi
end


text:
module Switch

export @switch

macro switch(value, body)
    cases = Any[]
    bodies = Any[]

    endlabel = gensym("switchend")

    hasdefault = false
    next_is_default = false
    default_label = nothing
    labels = Symbol[]

    for arg in body.args
        if isa(arg, Expr) && arg.head == :macrocall && arg.args[1] == symbol("@case")
            push!(cases, arg.args[2])
            push!(bodies, Any[])
            push!(labels, gensym("switchlabel"))
            next_is_default = false
        elseif isa(arg, Expr) && arg.head == :macrocall && arg.args[1] == symbol("@default")
            if hasdefault
                error("Multiple default cases in switch statement.")
            end
            hasdefault = true
            push!(cases, nothing)
            push!(bodies, Any[])
            default_label = gensym("switchlabel")
            push!(labels, default_label)
            next_is_default = true
        elseif !isempty(bodies)
            if isa(arg, Expr) && arg.head == :break
                body = Expr(:symbolicgoto, endlabel)
            else
                body = esc(arg)
            end
            push!(bodies[end], body)
        end
    end

    dispatch = Any[]
    for (case, label) in zip(cases, labels)
        if label == default_label
            continue
        end
        ex = quote
            if switch_value == $(case)
                $(Expr(:symbolicgoto, label))
            end
        end
        push!(dispatch, ex)
    end

    if hasdefault
        push!(dispatch, Expr(:symbolicgoto, default_label))
    end

    labeledbody = Any[]
    for (body, label) in zip(bodies, labels)
        push!(labeledbody, Expr(:symboliclabel, label))
        push!(labeledbody, body...)
    end

    quote
        begin
            local switch_value = $(esc(value))
            $(dispatch...)
            $(Expr(:symbolicgoto, endlabel))
            $(labeledbody...)
            $(Expr(:symboliclabel, endlabel))
        end
    end
end

end


text:
using ChainRulesCore: NO_FIELDS

# ChainRulesCore.frule(Δ, ::typeof(getproperty), x::ComponentArray, s::Symbol) = frule((_, Δ), getproperty, x, Val(s))
# function ChainRulesCore.frule(Δ, ::typeof(getproperty), x::ComponentArray, ::Val{s}) where s
#     zero_x = zero(x)
#     setproperty!(zero_x, s, Δ)
#     return getproperty(x, s), zero_x
# end
ChainRulesCore.rrule(::typeof(getproperty), x::ComponentArray, ::Val{s}) where s = ChainRulesCore.rrule(getproperty, x, s)
function ChainRulesCore.rrule(::typeof(getproperty), x::ComponentArray, s::Symbol)
    function getproperty_adjoint(Δ)
        zero_x = zero(x)
        setproperty!(zero_x, s, Δ)
        return (ChainRulesCore.NO_FIELDS, zero_x)
    end

    return getproperty(x, s), getproperty_adjoint
end

ChainRulesCore.rrule(::typeof(getdata), x::ComponentArray) = getdata(x), Δ->ComponentArray(Δ, getaxes(x))

ChainRulesCore.rrule(::Type{ComponentArray}, data, axes) = ComponentArray(data, axes), Δ->(NO_FIELDS, getdata(Δ), getaxes(Δ))

ChainRulesCore.rrule(::Type{ComponentArray}, data, axes) = ComponentArray(data, axes), Δ->(getdata(Δ), getaxes(Δ))


text:
module Walks

import ..WalksBase: get_position, get_x, get_y, get_z, get_time, get_nsteps, step!
import ..LatticeVars.init!
using ..WalksBase
using ..Actors
using ProgressMeter

export walk!, WalkPlan, trial!

function walk!(walk::AbstractWalkGeneral, actor::AbstractActor)
    init!(walk)
    Actors.init!(actor)
    while step!(walk) && Actors.act!(actor, walk)
    end
    return walk
end

struct WalkPlan{T, V}
    walk::T
    actor::V
end

for f in (:get_position, :get_x, :get_y, :get_z, :get_time, :get_nsteps)
    @eval ($f)(lwp::WalkPlan) = ($f)(lwp.walk)
end

function Base.show(io::IO, lp::WalkPlan)
    println(io, "WalkPlan")
    show(io, lp.walk)
    show(io, lp.actor)
end

walk!(lwp::WalkPlan) = walk!(lwp.walk, lwp.actor)

function sample!(walk_plan::WalkPlan, trial_actor)
    walk!(walk_plan)
    Actors.act!(trial_actor, walk_plan.walk)
    return nothing
end

function trial_body!(walk_plan::WalkPlan, sample_loop::SampleLoopActor, iter)
    for _ in iter
        sample!(walk_plan, sample_loop.actor)
    end
    return nothing
end

function trial_body!(walk_plan::WalkPlan, sample_loop::SampleLoopActor, progress_iter::ProgressIter)
    @showprogress 1 "Doing trials..." for i in progress_iter.iter
        sample!(walk_plan, sample_loop.actor)
    end
    return nothing
end

function trial!(walk_plan::WalkPlan, sample_loop::SampleLoopActor)
    trial_actor = sample_loop.actor
    Actors.init!(trial_actor)
    trial_body!(walk_plan, sample_loop, sample_loop.iter)
    Actors.finalize!(trial_actor)
    return nothing
end

end # module Walks


text:
instructions = readlines("day12/input.txt")

mutable struct Waypoint
    x::Int
    y::Int
end

mutable struct Ship
    x::Int
    y::Int
end

waypoint = Waypoint(10, 1)
ship = Ship(0,0)

function step_instruction(instruction)
    global waypoint
    command, amount = instruction[1], parse(Int, instruction[2:end])
    if command == 'N'
        waypoint.y += amount
    elseif command == 'S'
        waypoint.y -= amount
    elseif command == 'E'
        waypoint.x += amount
    elseif command == 'W'
        waypoint.x -= amount
    elseif command == 'L'
        for i in 1:(amount ÷ 90)
            waypoint.x, waypoint.y = (-waypoint.y, waypoint.x)
        end
    elseif command == 'R'
        for i in 1:(amount ÷ 90)
            waypoint.x, waypoint.y = (waypoint.y, -waypoint.x)
        end
    elseif command == 'F'
        for i in 1:amount
            ship.x += waypoint.x
            ship.y += waypoint.y
        end
    end
end

for ins in instructions
    step_instruction(ins)
end

println(abs(ship.x)+abs(ship.y))

# row row row your boat


text:
@testset "kronecker_generator" begin
  ij =  kronecker_generator(2,2)
  @test size(ij) == (2,8)
  ij1 = kronecker_generator(2,2;replicate=true,seed=[10,4,78,12,5,72])
  ij2 = kronecker_generator(2,2;replicate=true,seed=[10,4,78,12,5,72])
  @test ij1 == ij2
  @test get_min_type(4) == UInt8
  @test get_min_type(63) == UInt64
end


text:
# julia experiments\to_jld2.jl -i "..\..\..\lisa_data\nrn10.pklz" -o nrn10.jld2
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--input_path", "-i"
            help = "Input path"
			default = nothing
        "--output_path", "-o"
            help = "output path"
			default = nothing
    end

    return parse_args(s)
end
args = parse_commandline()

# using Test
using Logging
using Pio3d
using JLD2
ipath = args["input_path"]
opath = args["output_path"]

if opath == nothing
	opath = "$ipath.jld2"
end

@info "Reading $ipath..."
datap_in = Pio3d.read3d(ipath)

# @info "datap" datap_in

datap = Dict()
datap["data3d"] = datap_in["data3d"]
datap["voxelsize_mm"] = datap_in["voxelsize_mm"]

@info "Writing $opath..."
@JLD2.save opath datap


text:
function validedges(
    triangles::Array{Array{Point2D,1},1},
    del::DataFrame,
    scaler::Scaler,
    fd::Function,
    geps::Float64,
)::Array{GeometricalPredicates.Line2D{GeometricalPredicates.Point2D},1}
    inedges = Array{Array{Float64,2},1}()
    outedges = Array{Array{Float64,2},1}()
    for triangle in triangles
        center = unscaledcenterpoint(triangle, scaler)
        edges = vectorizededges(triangle)
        if fd([getx(center), gety(center)]) > -geps
            push!(inedges, edges...)
        else
            push!(outedges, edges...)
        end
    end
    return buildedges(inedges, outedges, del)
end

function unscaledcenterpoint(triangle::Array{Point2D,1}, scaler::Scaler)::Point2D
    return unscaledpoint2d(
        centroid(Primitive(triangle[1], triangle[2], triangle[3])),
        scaler,
    )
end

function vectorizededges(triangle::Array{Point2D,1})::Array{Array{Float64,2},1}
    a = triangle[1]
    b = triangle[2]
    c = triangle[3]
    ab = [getx(a) gety(a); getx(b) gety(b)]
    ba = [getx(b) gety(b); getx(a) gety(a)]
    bc = [getx(b) gety(b); getx(c) gety(c)]
    cb = [getx(c) gety(c); getx(b) gety(b)]
    ac = [getx(a) gety(a); getx(c) gety(c)]
    ca = [getx(c) gety(c); getx(a) gety(a)]
    return [ab, ba, bc, cb, ac, ca]
end

function buildedges(
    inedges::Array{Array{Float64,2},1},
    outedges::Array{Array{Float64,2},1},
    del::DataFrame,
)::Array{GeometricalPredicates.Line2D{GeometricalPredicates.Point2D},1}
    inedges = [edge for edge in inedges if !(edge in outedges)]
    return [
        Line(Point(r[1], r[2]), Point(r[3], r[4]))
        for r in eachrow(del) if !([r[1] r[2]; r[3] r[4]] in inedges)
    ]
end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
using GR

const x = [ 0, 1, 0, 1, 0.707 ]
const y = [ 0, 0, 1, 1, 0.707]
const z = [ 0, 0, 0, 0, 0.707 ]

const connections = [
    [ 4, 1, 2, 4, 3 ],  # 1 quadrilateral
    [ 3, 1, 2, 5 ],     # 4 triangles
    [ 3, 2, 4, 5 ],
    [ 3, 4, 3, 5 ],
    [ 3, 3, 1, 5 ]
]
#                               AABBGGRR
const colors = signed.(UInt32[0xdf753e80, 0xdf0068ff, 0xdfc0a35a, 0xdf2000c1, 0xdf62a2ce])

setviewport(0, 1, 0, 1)
setwindow(-1, 1, -1, 1)
setwindow3d(0, 1, 0, 1, 0, 1)
setspace3d(30, 80, 0, 0)

polygonmesh3d(x, y, z, vcat(connections...), colors)

updatews()


text:
_bls_by_shape(bls,s::Square) = fix_bls_by_shape!(bls,Rectangle(s))
function _bls_by_shape(bls,shape::Rectangle)
	stop = shape.x0-shape.a/2
	start = stop + bls[1].depth
	bls[1] = bls[1](startx=start, stopx=stop)

	stop = shape.x0+shape.a/2
	start = stop - bls[2].depth
	bls[2] = bls[2](startx=start, stopx=stop)

	stop = shape.y0-shape.b/2
	start = stop + bls[3].depth
	bls[3] = bls[3](starty=start, stopy=stop)

	stop = shape.y0+shape.b/2
	start = stop - bls[4].depth
	bls[4] = bls[4](starty=start, stopy=stop)
end
function _bls_by_shape(bls,c::Circle)
	stop=c.R
	start=stop-bls[1].depth
	bls[1] = bls[1](startx=start,stopx=stop)
end

#
# """
# 	plot(::Boundary)
#
# PML's plotted in blue patches (b/c absorbing, so cold), cPML's in red (b/c emitting, so hot)
#
# DirichletBC are solid black lines, NeumannBC dashed, FloquetBC dashdot, MatchedBC/noBC nothing
# """
# @recipe function f(bnd::Boundary)
# 	aspect_ratio --> 1
# 	legend --> false
#
# 	@series begin
# 		fillcolor --> BOUNDARY_COLOR
# 		color --> BOUNDARY_COLOR
# 		markershape := :none
# 		bnd.shape
# 	end
# 	@series (bnd,1)
# 	@series (bnd,1,1)
# end
# # plot boundary layers
# @recipe function f(bnd::Boundary,_1::Int)
# 	shape = bnd.shape
# 	cosθ = shape.cosθ
# 	sinθ = shape.sinθ
# 	typeof(shape)<:Square ? shape = Rectangle(shape) : nothing
# 	for i ∈ eachindex(bnd.bls)
# 		bls = bnd.bls[i]
# 		depth = bnd.bls[i].depth
# 		if typeof(shape)<:Rectangle
# 			if i==1
# 				x = -(shape.a-depth)/2
# 				y = 0
# 				bl_shape = Rectangle(depth,shape.b,shape.x0+cosθ*x-sinθ*y,shape.y0+sinθ*x+cosθ*y,shape.θ)
# 			elseif i==4
# 				x = 0
# 				y = (shape.b-depth)/2
# 				bl_shape = Rectangle(shape.a,depth,shape.x0+cosθ*x-sinθ*y,shape.y0+sinθ*x+cosθ*y,shape.θ)
# 			elseif i==2
# 				x = +(shape.a-depth)/2
# 				y = 0
# 				bl_shape = Rectangle(depth,shape.b,shape.x0+cosθ*x-sinθ*y,shape.y0+sinθ*x+cosθ*y,shape.θ)
# 			else
# 				x = 0
# 				y = -(shape.b-depth)/2
# 				bl_shape = Rectangle(shape.a,depth,shape.x0+cosθ*x-sinθ*y,shape.y0+sinθ*x+cosθ*y,shape.θ)
# 			end
# 		elseif typeof(shape)<:Circle
# 			bls = bnd.bls[1]
# 			depth = bls.depth
# 			bl_shape = Annulus(shape.R-depth,shape.R,shape.x0,shape.y0)
# 		end
# 		if typeof(bls)<:AbstractComplexBL
# 			@series begin
# 				alpha --> 0
# 				typeof(bls)<:PML ? fillcolor --> BOUNDARY_PML_COLOR : fillcolor --> BOUNDARY_CPML_COLOR
# 				typeof(bls)<:PML ? color --> BOUNDARY_PML_COLOR : color --> BOUNDARY_CPML_COLOR
# 				bl_shape
# 			end
# 		end
# 	end
# end
# # plot boundary conditions
# @recipe function f(bnd::Boundary,_1::Int,_2::Int)
# 	shape = bnd.shape
# 	for i ∈ eachindex(bnd.bcs)
# 		bcs = bnd.bcs[i]
# 		if typeof(shape)<:Rectangle
# 			if i==1
# 				x = [-1,-1]*shape.a/2
# 				y = [-1,+1]*shape.b/2
# 			elseif i==2
# 				x = [-1,+1]*shape.a/2
# 				y = [+1,+1]*shape.b/2
# 			elseif i==3
# 				x = [+1,+1]*shape.a/2
# 				y = [-1,+1]*shape.b/2
# 			else
# 				x = [-1,+1]*shape.a/2
# 				y = [-1,-1]*shape.b/2
# 			end
# 			if typeof(bcs)<:AbstractLocalBC
# 				@series begin
# 					if typeof(bcs)<:DirichletBC
# 						ls --> BOUNDARY_DIRICHLET_LINETYPE
# 					elseif typeof(bcs)<:NeumannBC
# 						ls --> BOUNDARY_NEUMANN_LINETYPE
# 					else
# 						ls --> BOUNDARY_FLOQUET_MATCHED_LINETYPE
# 					end
# 					seriestype --> :line
# 					markershape := :none
# 					color --> BOUNDARY_BC_COLOR
# 					shape.x0 .+ shape.cosθ*x-shape.sinθ*y, shape.y0 .+ shape.sinθ*x+shape.cosθ*y
# 				end
# 			end
# 		else
# 			if typeof(bcs)<:Union{AbstractLocalBC,FloquetBC}
# 				@series begin
# 					if typeof(bcs)<:DirichletBC
# 						ls --> BOUNDARY_DIRICHLET_LINETYPE
# 					elseif typeof(bcs)<:NeumannBC
# 						ls --> BOUNDARY_NEUMANN_LINETYPE
# 					else
# 						ls --> BOUNDARY_FLOQUET_MATCHED_LINETYPE
# 					end
# 					color --> BOUNDARY_BC_COLOR
# 					markershape := :none
# 					fillalpha --> 0
# 					alpha --> 1
# 					shape
# 				end
# 			end
# 		end
# 	end
# end
#


text:
"""
    Edge

An abstract edge for working with graph

* `getindex` : get the first and second nodes
* `string` : pretty print
* `reverse` : reverse the 2 ends

"""
abstract type Edge end

function Base.getindex(edge::Edge, i::Int64) #where {E<:Edge}
    return edge.vertices[i]
end

function Base.reverse(e::E) where {E<:Edge}
    return E(Base.reverse(e.vertices))
end

Base.string(e::Edge) = "$(e[1]) -> $(e[2])"
Base.show(io::IO, x::Edge) = print(io, string(x))

"""
    SimpleEdge <: Edge

A simple implementation of `Edge` which only has 2 points
"""
struct SimpleEdge <: Edge
    # u::Int64
    # v::Int64
    vertices::Tuple{Int64,Int64}
    function SimpleEdge(vertices)
        new(vertices)
    end
    function SimpleEdge(u,v)
        new((u,v))
    end
end

"""
    WeightedEdge <: Edge

A weighted version of `Edge` that has 2 points and a weight

* `weight`
* `reverse`
"""
struct WeightedEdge <: Edge
    vertices::Tuple{Int64,Int64}
    weight::Float64
    function WeightedEdge(vertices::Tuple{Int64,Int64}, weight)
        new(vertices, weight)
    end
    function WeightedEdge(u::Int64, v::Int64, weight)
        new((u,v), weight)
    end
end

Base.:(<)(x::WeightedEdge, y::WeightedEdge) = error("Unsupported operation")
Base.:(==)(x::WeightedEdge, y::WeightedEdge) = x.vertices == y.vertices #error("Unsupported operation")
# for sum(), Java.totalWeight()
Base.:(+)(x::WeightedEdge, y::WeightedEdge) = weight(x) + weight(y)
Base.:(+)(x::Float64, y::WeightedEdge) = x + weight(y)

function weight(e::WeightedEdge)
    e.weight
end

function Base.reverse(e::WeightedEdge)
    WeightedEdge(e[2], e[1], weight(e))
end


text:
# # Causal interventions
using ..Var, ..Basis
export |ᵈ, intervene, mergetags, Intervention

abstract type AbstractIntervention end

"The imperative that `x` should be replaced to value `v`"
struct Intervention{X, V} <: AbstractIntervention
  x::X
  v::V
end

Intervention(x::Pair{X, <:Number}) where X = Intervention(x.first, ω -> x.second)
Intervention(x::Pair) = Intervention(x.first, x.second)

"Multiple variables intervened"
struct MultiIntervention{XS} <: AbstractIntervention
  is::XS
end

"Intervened Variable: `x` had intervention `i` been the case"
struct Intervened{X, I}
  x::X
  i::I
end

"Merge Interventions"
mergeinterventions(i1::Intervention, i2::Intervention) = MultiIntervention((i1, i2))
mergeinterventions(i1::Intervention, i2::MultiIntervention) = MultiIntervention((i1, i2.is...))
mergeinterventions(i1::MultiIntervention, i2::Intervention) = MultiIntervention((i1.is..., i2))

"Merge Intervention Tags"
function mergetags(nt1::NamedTuple{K1, V1}, nt2::NamedTuple{K2, V2}) where {K1, K2, V1, V2}
  if K1 ∩ K2 == [:intervene]
    merge(merge(nt1, nt2), (intervene = mergeinterventions(nt2[:intervene], nt1[:intervene]),))
  else
    @assert false "Unimplemented"
  end
end

"intervened"
intervene(x, intervention::AbstractIntervention) = Intervened(x, intervention)
intervene(x, intervention::Pair) = Intervened(x, Intervention(intervention))
intervene(x, interventions::Tuple) =
  Intervened(x, MultiIntervention(map(Intervention, interventions)))

@inline x |ᵈ i = intervene(x, i)


text:
@testset "Minima and Kissing numbers" begin
  D = lattice_database()
  _minima = fmpq[ minimum(lattice(D, i)) for i in 1:300]
  @test _minima == fmpq[2, 1, 1, 1, 2, 2, 1, 2, 2, 1, 2, 3, 2, 3, 1, 2, 2, 3,
                        6, 1, 2, 1, 2, 3, 4, 2, 2, 1, 2, 1, 2, 2, 2, 1, 2, 2,
                        2, 1, 2, 3, 4, 2, 2, 4, 2, 4, 2, 2, 2, 4, 4, 2, 2, 4,
                        2, 3, 2, 4, 4, 6, 2, 2, 3, 1, 2, 5, 3, 4, 2, 4, 2, 2,
                        4, 4, 4, 7, 4, 2, 6, 4, 4, 4, 2, 4, 2, 4, 2, 3, 3, 2,
                        7, 2, 4, 2, 3, 2, 2, 4, 12, 2, 4, 6, 6, 8, 6, 8, 6, 6,
                        8, 6, 3, 6, 6, 4, 6, 6, 6, 4, 4, 4, 4, 4, 4, 4, 4, 2,
                        2, 4, 6, 4, 4, 10, 6, 4, 6, 4, 8, 2, 8, 2, 4, 2, 2, 2,
                        2, 4, 8, 4, 8, 3, 3, 4, 6, 6, 2, 9, 2, 4, 4, 16, 4, 4,
                        4, 3, 3, 2, 2, 10, 4, 6, 3, 2, 4, 4, 16, 4, 3, 4, 2,
                        11, 60, 2, 4, 4, 12, 4, 4, 4, 8, 4, 4, 8, 4, 2, 12, 2,
                        8, 2, 4, 6, 2, 2, 4, 2, 4, 2, 2, 2, 2, 2, 2, 15, 12, 4,
                        4, 4, 4, 4, 3, 4, 2, 13, 3, 5, 2, 4, 4, 9, 4, 4, 4, 3,
                        2, 14, 4, 4, 4, 4, 2, 4, 2, 4, 4, 4, 4, 4, 4, 7, 6, 6,
                        2, 15, 2, 3, 2, 4, 4, 4, 4, 4, 4, 4, 3, 3, 8, 6, 8, 12,
                        12, 8, 10, 6, 10, 4, 2, 16, 6, 8, 4, 4, 4, 4, 8, 6, 8,
                        4, 6, 2, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 2 ]

  _kissing = Int[Hecke.kissing_number(lattice(D, i)) for i in 1:300]
  @test _kissing == [2, 2, 2, 2, 6, 6, 2, 4, 4, 4, 12, 8, 12, 8, 6, 4, 12, 8,
                     8, 6, 6, 2, 4, 4, 6, 2, 2, 2, 2, 2, 8, 2, 2, 2, 2, 4, 4,
                     4, 4, 4, 4, 6, 6, 10, 20, 10, 24, 24, 24, 8, 6, 2, 2, 4,
                     24, 4, 6, 12, 4, 12, 6, 4, 8, 8, 30, 12, 20, 30, 40, 10,
                     12, 40, 12, 6, 30, 24, 6, 42, 14, 42, 44, 42, 60, 12, 72,
                     54, 72, 20, 8, 56, 16, 84, 14, 126, 56, 126, 126, 92, 2,
                     126, 64, 60, 56, 56, 58, 56, 60, 56, 56, 56, 56, 56, 56,
                     60, 56, 56, 56, 60, 60, 60, 58, 72, 56, 56, 56, 56, 84,
                     72, 64, 68, 64, 56, 120, 16, 48, 72, 30, 72, 18, 112, 16,
                     240, 240, 240, 240, 132, 12, 150, 42, 16, 16, 120, 48,
                     48, 90, 20, 144, 18, 180, 6, 198, 272, 270, 20, 20, 30,
                     110, 22, 220, 110, 40, 180, 20, 276, 54, 336, 80, 260,
                     132, 24, 24, 220, 22, 432, 82, 438, 432, 216, 270, 420,
                     336, 336, 360, 156, 26, 12, 126, 18, 126, 60, 24, 42,
                     756, 264, 24, 264, 60, 72, 36, 36, 36, 12, 2, 756, 756,
                     648, 632, 624, 64, 84, 182, 28, 52, 52, 312, 26, 918, 2,
                     906, 890, 888, 112, 210, 30, 378, 756, 210, 1260, 364,
                     28, 252, 1242, 1248, 1422, 1228, 1212, 1210, 156, 182,
                     546, 240, 32, 240, 240, 420, 30, 1746, 1644, 2340, 1872,
                     1630, 1596, 70, 240, 600, 960, 180, 1680, 480, 1200,
                     1440, 336, 336, 720, 272, 34, 1088, 360, 240, 4320, 4320,
                     440, 420, 2400, 360, 480, 360, 480, 32, 720, 2772, 2436,
                     2448, 4320, 2982, 2402, 2364, 512, 256, 306 ]
end


text:
using PQEcolaPoint
using Dates

# build a binary tree with three levels -
SF = 82.00
model = build(CRRLatticeModel; Sₒ = SF, number_of_levels=80, σ=0.4510, T=(37.0 / 365), μ=0.0045)

# what is the contract?
short_put_contract = PutContractModel()
short_put_contract.ticker = "XYZ"
short_put_contract.expiration_date = Date(2022, 05, 06)
short_put_contract.strike_price = (75.0)
short_put_contract.premium = 0.0
short_put_contract.number_of_contracts = 1
short_put_contract.direction = 1

# what is the premimum?
parray = premium(short_put_contract, model)


text:
myParse(x) = parse.(Int, x)
crabPositions = myParse.(split.(readline("input.txt"), ','))

#! Part 1

using Statistics

function findLeastFuelPositionPt1(positions)
    return round(median(positions))
end

totalFuelCostPt1(positions, destination) = sum([abs(x-destination) for x in positions])

print(totalFuelCostPt1(crabPositions, findLeastFuelPositionPt1(crabPositions)))

#! Part 2

function findLeastFuelCostPt2(positions)
    fuel = Inf
    approxPos = Int(round(mean(positions)))
    𝛜 = 3
    for pos in approxPos-𝛜:approxPos+𝛜
        dists = [abs(x - pos) for x in positions]
        f = sum([0.5*(d^2 + d) for d in dists])

        if fuel > f
            fuel = f
        end
    end

    return fuel
end

print(findLeastFuelCostPt2(crabPositions))

# Intuition for the least fuel Position functions: https://www.johnmyleswhite.com/notebook/2013/03/22/modes-medians-and-means-an-unifying-perspective/

max_stars_repo_path: Day 7/main.jl

max_issues_repo_path: Day 7/main.jl

max_forks_repo_path: Day 7/main.jl


text:
# --- PRIVATE METHODS --------------------------------------------------------------------------------------- #
function _compute_profit_loss_at_expiration(asset::VLAbstractAsset, assetPriceValueArray::Array{Float64,1})::Array{Float64,1}

    # initialize -
    pl_value_array = Array{Float64,1}()

    # compute the intrinsic value -
    for (index, asset_price_value) in enumerate(assetPriceValueArray)
        result = intrinsic_value(asset, asset_price_value)
        if (isa(result.value, Exception) == true)
            return result
        end
        pl_value = result.value.pl_value

        # grab -
        push!(pl_value_array, pl_value)
    end

    # call -
    return pl_value_array
end

function _compute_profit_loss_at_expiration(asset::VLAbstractAsset, assetPriceStart::Float64, assetPriceStop::Float64;
    number_of_price_steps::Int64=1000)::Array{Float64,1}

    # setup price range -
    asset_price_range = collect(range(assetPriceStart, assetPriceStop, length=number_of_price_steps))

    # return -
    return _compute_profit_loss_at_expiration(asset, asset_price_range)
end
# ----------------------------------------------------------------------------------------------------------- #

function expiration(assetSet::Set{VLAbstractAsset}, assetPriceArray::Array{Float64,1})::VLResult

    try

        # get -
        number_of_price_steps = length(assetPriceArray)

        # initialize -
        number_of_assets = length(assetSet)
        assetProfitLossArray = zeros(number_of_price_steps, number_of_assets)
        tradeProfitLossArray = zeros(number_of_price_steps, 2) # firstCol: price, secondCol: PL

        # process the set of assets -
        for (assetIndex, assetObject) in enumerate(assetSet)

            # process each asset type -
            profit_loss_array = _compute_profit_loss_at_expiration(assetObject, assetPriceArray)

            # add the PL for this asset to the overall array -
            for index = 1:number_of_price_steps
                assetProfitLossArray[index,assetIndex] = profit_loss_array[index]
            end
        end

        # ok, so the assetProfitLossArray holds the PL for each asset type in its cols -
        # to get the PL for the entire trade at a particular stock price, then sum across the cols -
        for price_index = 1:number_of_price_steps

            # grab col -
            tmp_row = assetProfitLossArray[price_index,:]

            # add -
            tradeProfitLossArray[price_index,1] = assetPriceArray[price_index]
            tradeProfitLossArray[price_index,2] = sum(tmp_row)
        end

        # return -
        return VLResult{Array{Float64,2}}(tradeProfitLossArray)
    catch error
        return VLResult(error)
    end
end


text:
################################################################################################
#
# demo 4: Live video capture (webcam)
#
################################################################################################

# Julia variables
device_index = CAP_ANY
winname = "Welcome"

# C++ OpenCV code
cxx"""
void videocapture(int device_index, const char *winname) {

   // convert const char* => string
   std::string wname = winname;

   // open Video device
   cv::VideoCapture capture(device_index);

   // check that device is opened
   if (!capture.isOpened())
    {
        std::cout<< "Can not open device!";
    }

 while(true) {
   // Create a Mat frame and show it in a window
    cv::Mat frame;
    cv::namedWindow(wname, 0);
    bool Success = capture.read(frame);

    if (!Success)
    {
        std::cout << "Failure to acquire stream!";
        break;
    }

    cv::imshow(wname, frame);

    if (cv::waitKey(30) == 27)
    {
        std::cout << "Press esc to end acquisition.";
        cv::destroyWindow(wname);
        break;
    }

  }
}
"""

@cxx videocapture(device_index, pointer(winname))

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


text:
# Generic methods for matrix correctio

using DataFrames
using PeriodicTable
using Roots

"""
`MatrixCorrection` is an abstract type for computing ϕ(ρz)-type matrix correction algorithms.  A sub-class
`MCA <: MatrixCorrection` should implement

    # Integral of the area under the ϕ(ρz)-curve
    F(mc::MCA)
    # Integral of the transmitted area under the ϕ(ρz)-curve
    Fχ(mc::MCA, xray::CharXRay, θtoa::Real)
    # Area under 0 to τ under the transmitted ϕ(ρz)-curve
    Fχp(mc::MCA, xray::CharXRay, θtoa::Real, τ::Real)
    # The ϕ(ρz)-curve
    ϕ(mc::MCA, ρz)
    # A factory method from MCA
    matrixcorrection(::Type{MCA}, mat::Material, ashell::AtomicSubShell, e0)::MCA

The algorithm should precompute as much as possible based on the `Material`, `AtomicSubShell` and beam
energy.  The class `MCA` should define member variables `subshell::AtomicSubShell`, `material::Material` and
`E0::AbstractFloat` to store the input values. See `XPP` for an example.

From these methods, other methods like `Z(...)`, `A(...)` are implemented.
"""
abstract type MatrixCorrection end

"""
    F(mc::MatrixCorrection)

Integral of the ϕ(ρz)-curve from ρz = 0 to ∞.
"""
F(mc::MatrixCorrection) = error("$mc does not implement F(mc::$mc)")
"""
    Fχ(mc::MatrixCorrection, xray::CharXRay, θtoa::Real)

Integral of the area under the absorption corrected ϕ(ρz)-curve from ρz = 0 to ∞.
"""
Fχ(mc::MatrixCorrection, xray::CharXRay, θtoa::Real)  = error("$mc does not implement Fχ(mc::$mc, xray::CharXRay, θtoa::Real)")

"""
    Fχp(mc::NeXLMatrixCorrection, xray::CharXRay, θtoa::Real, τ::Real)

The partial integral of the absorption corrected ϕ(ρz) curve from ρz = 0 to τ #
"""
Fχp(mc::MatrixCorrection, xray::CharXRay, θtoa::Real, τ::Real)  = error("$mc does not implement Fχp(mc::$mc, xray::CharXRay, θtoa::Real, τ::Real)")

"""
    ϕ(mc::MatrixCorrection, ρz)

The ϕ(ρz)-curve for visualization and other purposes.
"""
ϕ(mc::MatrixCorrection, ρz)  = error("$mc does not implement ϕ(mc::$mc, ρz)")

"""
    NeXLCore.atomicsubshell(mc::MatrixCorrection)

The sub-shell for which this `MatrixCorrection` has been calculated.
"""
NeXLCore.atomicsubshell(mc::MatrixCorrection) = mc.subshell

"""
    NeXLCore.material(mc::MatrixCorrection) = mc.material

The material for which this `MatrixCorrection` has been calculated.
"""
NeXLCore.material(mc::MatrixCorrection) = mc.material

"""
    beamEnergy(mc::MatrixCorrection) = mc.E0

The beam energy (eV) for which this `MatrixCorrection` has been calculated.
"""
beamEnergy(mc::MatrixCorrection) = mc.E0

"""
    χ(mat::Material, xray::CharXRay, θtoa)

Angle adjusted mass absorption coefficient.
"""
χ(mat::Material, xray::CharXRay, θtoa::AbstractFloat) = mac(mat, xray) * csc(θtoa)
χ(mat::Material, ea::Float64, θtoa::AbstractFloat) = mac(mat, ea) * csc(θtoa)

"""
    ϕabs(mc::MatrixCorection, ρz, xray::CharXRay, θtoa::AbstractFloat)

Computes the absorbed ϕ(ρz) curve according to the XPP algorithm.
"""
ϕabs(mc::MatrixCorrection, ρz, xray::CharXRay, θtoa::AbstractFloat) = ϕ(mc, ρz) * exp(-χ(material(mc), xray, θtoa) * ρz)

"""
    ZA(
      unk::MatrixCorrection,
      std::MatrixCorrection,
      xray::CharXRay,
      θunk::AbstractFloat,
      θstd::AbstractFloat
    )

The atomic number (Z) and absorption (A) correction factors.
"""
function ZA(unk::MatrixCorrection, std::MatrixCorrection, xray::CharXRay, θunk::AbstractFloat, θstd::AbstractFloat)
    @assert isequal(unk.subshell, inner(xray)) "Unknown and X-ray don't match in XPP"
    @assert isequal(std.subshell, inner(xray)) "Standard and X-ray don't match in XPP"
    return Fχ(unk, xray, θunk) / Fχ(std, xray, θstd)
end

"""
    Z(unk::MatrixCorrection, std::MatrixCorrection)

The atomic number correction factor.
"""
function Z(unk::MatrixCorrection, std::MatrixCorrection)
    @assert isequal(unk.subshell, std.subshell)
    "Unknown and standard matrix corrections don't apply to the same sub-shell."
    return F(unk) / F(std)
end

"""
    A(unk::MatrixCorrection, std::MatrixCorrection, xray::CharXRay, θunk::AbstractFloat, θstd::AbstractFloat)
    A(unk::MatrixCorrection, xray::CharXRay, θ::AbstractFloat)

The absorption correction factors.
"""
function A(unk::MatrixCorrection, std::MatrixCorrection, xray::CharXRay, θunk::AbstractFloat, θstd::AbstractFloat)
    @assert isequal(unk.subshell, inner(xray)) "Unknown and X-ray don't match in XPP"
    @assert isequal(std.subshell, inner(xray)) "Standard and X-ray don't match in XPP"
    return ZA(unk, std, xray, θunk, θstd) / Z(unk, std)
end

function A(unk::MatrixCorrection, xray::CharXRay, θunk::AbstractFloat)
    return Fχ(unk, xray, θunk)/F(unk)
end

correctcontinuum(mc::MatrixCorrection, θtoa::Real) = Fχ(mc, θtoa) / F(mc)

"""
    kcoating(ty::Type{<:MatrixCorrection}, subtrate::Material, coating::Material, cxr::CharXRay, e0::Real, toa::Real, τ::Real)

Estimate the k-ratio for a coating of mass-thickness τ (g/cm²) on the specified substrate. The standard for the coating
is assumed to be of the same material as the coating.
"""
function kcoating(
    ty::Type{<:MatrixCorrection},
    subtrate::Material,
    coating::Material,
    cxr::CharXRay,
    e0::Real,
    toa::Real,
    τ::Real,
)
    @assert element(cxr) in keys(coating) "$cxr must be produced by one of the elements in $(name(coating))"
    @assert e0 > energy(inner(cxr)) "The beam energy must exceed the edge energy for $cxr."
    submc, coatmc = matrixcorrection(ty, subtrate, inner(cxr), e0), matrixcorrection(ty, coating, inner(cxr), e0)
    return Fχp(submc, cxr, toa, τ) / Fχ(coatmc, cxr, toa)
end

""""
    massthickness(ty::Type{<:MatrixCorrection}, subtrate::Material, coating::Material, cxr::CharXRay, e0::Real, toa::Real, k::Real)

Estimate the mass-thickness of a ultra-thin layer of a `coating` material on a `substrate` from a measured k-ratio `k`
of a characteristic X-ray `cxr`.  Works for k-ratios of the order of 1 %.  The standard for the coating is assumed
to be of the same material as the coating.
"""
function NeXLCore.massthickness(
    ty::Type{<:MatrixCorrection},
    substrate::Material,
    coating::Material,
    cxr::CharXRay,
    e0::Real,
    toa::Real,
    k::Real,
)
    submc, coatingmc = matrixcorrection(ty, substrate, inner(cxr), e0), matrixcorrection(ty, coating, inner(cxr), e0)
    f(τ) = k - Fχp(submc, cxr, toa, τ) / Fχ(coatingmc, cxr, toa)
    return find_zero(f, 0.01 * range(ty, substrate, e0, false), Roots.Order1())
end

function coatingasfilm(
    ty::Type{<:MatrixCorrection},
    substrate::Material,
    coating::Material,
    cxr::CharXRay,
    e0::Real,
    toa::Real,
    k::Real,)
    @assert haskey(coating.properties, :Density)
    return Film(coating, massthickness(ty, substrate, coating, cxr, e0, toa, k) / coating[:Density])
end

"""
    correctkratios(krs::AbstractVector{KRatio}, coating::Material, θtoa::Real, ρz::Real)::Vector{KRatio}

This function is mainly for pedagogical purposes.  It takes a `KRatio[]`, a coating `Material` on the unknown,
and a mass-thickness (g/cm²) and creates a new `KRatio[]` that accounts for the intensity missing due to absorption
by the coating.  Favor the coating correction built into `ZAFCorrection` or `MultiZAF`.
"""
function correctkratios(krs::AbstractVector{KRatio}, coating::Material, ρz::Real)::Vector{KRatio}
    res = KRatio[]
    for kr in krs
        if !(kr.element in keys(coating))
            kratio = kr.kratio / exp(-χ(coating, brightest(kr.xrays), kr.unkProps[:TakeOffAngle]) * ρz)
            push!(res, KRatio(kr.xrays, kr.unkProps, kr.stdProps, kr.standard, kratio))
        else
            push!(res, kr)
        end
    end
    return res
end


text:
export Dirac

struct Dirac{X} <: PrimitiveMeasure
    x::X
end

sampletype(μ::Dirac{X}) where {X} = X

function (μ::Dirac{X})(s) where {X}
    μ.x ∈ s && return 1
    return 0
end

logdensity(μ::Dirac, x) = (x == μ.x) ? 0.0 : -Inf

Base.rand(::Random.AbstractRNG, T::Type, μ::Dirac) = μ.x

export dirac

dirac(d::AbstractMeasure) = Dirac(rand(d))

function logdensity(μ::Dirac{M}, ν::Dirac{M}, x) where {M}
    if μ.x == ν.x
        x == μ.x && return 0.0
        return -Inf
    elseif μ.x == x
        return Inf
    else
        return -Inf
    end
end

testvalue(d::Dirac) = d.x


text:
using NumbersFromText
using Base.Test

dir = joinpath(dirname(@__FILE__),"test_files/")

@testset "SpaceTokenizer" begin

end

array_int_float_filename = joinpath(dir, "arrays_int_float.txt")
array_int_float_contents = Any[
Vector{Int}([1 , 2 , 3 , 4 , 5 , 6 , -1, -2]),
Vector{Float64}([3.0, 4.0, 5.0, 6.0, -1.0, 1.0e18, -1.0e18, 1.23412341234123512341234e-000008])
]

@testset "readarray" begin
  data = [1/3*ones(500); pi*ones(500); nextfloat(0.0)*ones(500)]
  buf = IOBuffer()
  for i in data
    print(buf, i)
    print(buf, "  ")
  end
  seek(buf, 0) # reset buffer
  a = readarray!(buf, zeros(0))
  @test a == data

  data = convert(Array{UInt8}, "1.0"*("0"^2048)*"      2.0")
  buf = IOBuffer(data)
  @test_throws ArgumentError readarray(buf)
  seek(buf,0)
  @test_throws ArgumentError readarray(buf; maxbuf=2048)
  seek(buf,0)
  a = readarray(buf; maxbuf=2052)
  @test a == [1.0,2.0]

  data = [1/3*ones(500); pi*ones(500); nextfloat(0.0)*ones(500)]
  buf = IOBuffer()
  for i in data
    print(buf, i)
    print(buf, " ")
  end
  seek(buf, 0) # reset buffer
  a = readarray!(buf, zeros(0))
  @test a == data

  data = [1/3*ones(500); pi*ones(500); nextfloat(0.0)*ones(500)]
  buf = IOBuffer()
  for i in data
    print(buf, i)
    print(buf, "  ")
  end
  seek(buf, 0) # reset buffer
  a = readarray!(buf, zeros(0))
  @test a == data

  data = ones(Int, 500)
  buf = IOBuffer()
  foreach(x -> println(buf, x), data)
  seek(buf, 0)
  a = readarray!(buf, zeros(Int, 0))
  @test a == data
  seek(buf, 0)
  a = readarray(Int, buf)

  a = readarray(array_int_float_filename)
  @test a == vec(hcat(array_int_float_contents...)')

end

@testset "readarrays" begin
  data1 = 1/3*ones(500)
  data2 = ones(Int, 500)
  buf = IOBuffer()
  for i in eachindex(data1)
    print(buf, data1[i] )
    print(buf, "  ")
    print(buf, data2[i] )
    print(buf, "  ")
  end
  seek(buf, 0) # reset buffer
  a1, a2 = readarrays!(buf, zeros(0), zeros(Int, 0))
  @test a1 == data1
  @test a2 == data2

  seek(buf, 0) # reset buffer
  a1, a2 = readarrays(buf, Float64, Int)
  @test a1 == data1
  @test a2 == data2

  data = convert(Array{UInt8}, "1.0"*("0"^2048)*"      2.0")
  buf = IOBuffer(data)
  @test_throws ArgumentError readarrays(buf, Float64)
  seek(buf,0)
  @test_throws ArgumentError readarrays(buf, Float64; maxbuf=2048)
  seek(buf,0)
  a, = readarrays(buf, Float64; maxbuf=2052)
  @test a == [1.0,2.0]

  a1,a2 = readarrays(array_int_float_filename, Int, Float64)
  @test Any[a1,a2] == array_int_float_contents

end

@testset "readmatrix" begin
  data = randn(3,2)
  buf = IOBuffer()
  writedlm(buf, data, "\t ")
  seek(buf, 0) # reset
  a = readmatrix(buf)
  @test a == data

  buf = IOBuffer("1.0")
  a = readmatrix(buf)
  @test a == ones(1,1)

  data = convert(Array{UInt8}, "1.0"*("0"^2048))
  buf = IOBuffer(data)
  @test_throws ArgumentError a = readmatrix(buf)

  data = convert(Array{UInt8}, "1.0 1.0\n1.0"*("0"^2048)*"      2.0")
  buf = IOBuffer(data)
  a = readmatrix(buf; maxbuf = 2052)
  @test a == [1.0 1.0; 1.0 2.0]

  data = convert(Array{UInt8}, "1.0 1.0\n1.0"*("0"^2048)*"      2.0")
  buf = IOBuffer(data)
  a = readmatrix(buf; maxbuf = 2052, transpose=false)
  @test a == [1.0 1.0; 1.0 2.0]'

  a = readmatrix(joinpath(array_int_float_filename))
  @test a == hcat(array_int_float_contents...)
end

@testset "partition_buffer" begin
  data = collect(1:500)
  buf = IOBuffer()
  foreach(x -> println(buf, x), data)
  seek(buf, 0)
  dataarray = take!(buf)

  maxparts = 49
  bufs = allocate_buffers(maxparts)

  for nparts=1:maxparts
    n = partition_buffer(dataarray, length(dataarray), nparts, bufs, UInt8('\n'))
    newdata = Vector{UInt8}()
    for i=1:n
      push!(newdata, read(bufs[i])...)
    end
    @test dataarray == newdata
  end

end

@testset "parallel_readarrays" begin
  #=
  data = rand(Int, 32*1024*1024)
  buf = IOBuffer()
  foreach(x -> println(buf, x), data)
  seek(buf, 0)
  a, = readarrays!(Val{true}, buf, zeros(Int, 0))
  @test a == data
  =#
end