text:
# This file hosts functions that mimic some GEOS functions and interfaces the use off
# GMT and GDAL data types.

"""
    setproj!(type, proj)

Set a referencing system to the `type` object. This object can be a `GMTgrid`, a `GMTimage`,
a `GMTdataset` or an `AbstractDataset`.

- `proj` Is either a Proj4 string or a WKT. Alternatively, it can also be another grid, image or dataset
         type, in which case its referencing system is copied into `type`
"""
function setproj!(tipo::AbstractArray, proj::String="")
	(!isa(tipo, GMTgrid) && !isa(tipo, GMTimage) && !isa(tipo, GMT.GMTdataset) && !isa(tipo, Vector{GMT.GMTdataset})) &&
		error("Wrong data type for this function. Must be a grid, image or dataset")
	(proj == "") && error("the projection string cannot obviously be empty")
	isproj4 = (startswith(proj, "+proj") !== nothing)
	obj = (isa(tipo, Vector{GMT.GMTdataset})) ? tipo[1] : tipo
	(isproj4) ? (obj.proj4 = proj) : (obj.wkt = proj)
	return nothing
end
function setproj!(tipo::AbstractArray, ref)
	(!isa(ref, GMTgrid) && !isa(ref, GMTimage) && !isa(ref, GMT.GMTdataset) && !isa(ref, Vector{GMT.GMTdataset})) &&
		error("Wrong REFERENCE data type for this function. Must be a grid, image or dataset")
	obj = (isa(ref, Vector{GMT.GMTdataset})) ? ref[1] : ref
	((prj = obj.proj4) == "") && (prj = obj.wkt)
	(prj == "") && error("The REFERENCE type is not referenced with either PROJ4 or WKT string")
	setproj!(tipo, prj)
end
function setproj!(dataset::AbstractDataset, projstring::AbstractString)
	result = GDALSetProjection(dataset.ptr, projstring)
	@cplerr result "Could not set projection"
	return dataset
end

# ---------------------------------------------------------------------------------------------------
"""
    buffer(geom, dist::Real, quadsegs::Integer=30; gdataset=false)

### Parameters
* `geom`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `dist`: the buffer distance to be applied. Should be expressed into the
    same unit as the coordinates of the geometry.
* `quadsegs`: the number of segments used to approximate a 90 degree (quadrant) of curvature.
* `gdataset`: Returns a GDAL IGeometry even when input is a GMTdataset or Matrix

Compute buffer of geometry.

Builds a new geometry containing the buffer region around the geometry on which
it is invoked. The buffer is a polygon containing the region within the buffer
distance of the original geometry.

Some buffer sections are properly described as curves, but are converted to
approximate polygons. The `quadsegs` parameter can be used to control how many
segments should be used to define a 90 degree curve - a quadrant of a circle.
A value of 30 is a reasonable default. Large values result in large numbers of
vertices in the resulting buffer geometry while small numbers reduce the
accuracy of the result.

### Returns
A GMT dataset when input is a Matrix or a GMT type (except if `gdaset=true`), or a GDAL IGeometry otherwise
"""
buffer(geom::AbstractGeometry, dist::Real, quadsegs::Integer=30) = IGeometry(OGR_G_Buffer(geom.ptr, dist, quadsegs))
function buffer(D::AbstractArray, dist::Real, quadsegs::Integer=30; gdataset=false)
	geom = helper_1geom(D)
	ig = IGeometry(OGR_G_Buffer(geom.ptr, dist, quadsegs))
	return (gdataset) ? ig : gd2gmt(ig)
end

# ---------------------------------------------------------------------------------------------------
"""
    centroid(geom; gdataset=false)

### Parameters
* `geom`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `gdataset`: Returns a GDAL IGeometry even when input is a GMTdataset or Matrix

Compute the geometry centroid.

The centroid is not necessarily within the geometry.

(This method relates to the SFCOM ISurface::get_Centroid() method however the current implementation
based on GEOS can operate on other geometry types such as multipoint, linestring, geometrycollection
such as multipolygons. OGC SF SQL 1.1 defines the operation for surfaces (polygons). SQL/MM-Part 3
defines the operation for surfaces and multisurfaces (multipolygons).)

### Returns
A GMT dataset when input is a Matrix or a GMT type (except if `gdaset=true`), or a GDAL IGeometry otherwise
"""
function centroid(geom::AbstractGeometry)
	point = createpoint()
	centroid!(geom, point)
	return point
end
centroid(D; gdataset=false) = helper_geoms_run_fun(centroid, D; gdataset=gdataset)

function centroid!(geom::AbstractGeometry, centroid::AbstractGeometry)
	result = OGR_G_Centroid(geom.ptr, centroid.ptr)
	@ogrerr result "Failed to compute the geometry centroid"
	return centroid
end

# ---------------------------------------------------------------------------------------------------
"""
    intersection(geom1, geom2; gdataset=false)

### Parameters
* `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type
* `gdataset`: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

Returns a new geometry representing the intersection of the geometries, or NULL
if there is no intersection or an error occurs.

Generates a new geometry which is the region of intersection of the two geometries operated on.
The `intersects` function can be used to test if two geometries intersect.

### Returns
A GMT dataset when input is a Matrix or a GMT type (except if `gdaset=true`), or a GDAL IGeometry otherwise
"""
intersection(g1::AbstractGeometry, g2::AbstractGeometry) = IGeometry(OGR_G_Intersection(g1.ptr, g2.ptr))
intersection(D1, D2; gdataset=false) = helper_geoms_run_fun(intersection, D1, D2; gdataset=gdataset)

# ---------------------------------------------------------------------------------------------------
"""
    polyunion(geom1, geom2; gdataset=false)

### Parameters
* `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type
* `gdataset`: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

Computes a new geometry representing the union of the geometries.

### Returns
A GMT dataset when input is a Matrix or a GMT type (except if `gdaset=true`), or a GDAL IGeometry otherwise
"""
polyunion(g1::AbstractGeometry, g2::AbstractGeometry) = IGeometry(OGR_G_Union(g1.ptr, g2.ptr))
polyunion(D1, D2; gdataset=false) = helper_geoms_run_fun(polyunion, D1, D2; gdataset=gdataset)

# ---------------------------------------------------------------------------------------------------
"""
    difference(geom1, geom2; gdataset=false)

### Parameters
* `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type
* `gdataset`: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

Generates a new geometry which is the region of this geometry with the region of the other geometry removed.

### Returns
A GMT dataset when input is a Matrix or a GMT type (except if `gdaset=true`), or a GDAL IGeometry otherwise
"""
difference(g1::AbstractGeometry, g2::AbstractGeometry) = IGeometry(OGR_G_Difference(g1.ptr, g2.ptr))
difference(D1, D2; gdataset=false) = helper_geoms_run_fun(difference, D1, D2; gdataset=gdataset)

# ---------------------------------------------------------------------------------------------------
"""
    symdifference(geom1, geom2; gdataset=false)

### Parameters
* `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type
* `gdataset`: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

Generates a new geometry representing the symmetric difference of the geometries
or NULL if the difference is empty or an error occurs.

### Returns
A GMT dataset when input is a Matrix or a GMT type (except if `gdaset=true`), or a GDAL IGeometry otherwise
"""
symdifference(g1::AbstractGeometry, g2::AbstractGeometry) = IGeometry(OGR_G_SymDifference(g1.ptr, g2.ptr))
symdifference(D1, D2; gdataset=false) = helper_geoms_run_fun(symdifference, D1, D2; gdataset=gdataset)

# ---------------------------------------------------------------------------------------------------
"""
    distance(geom1, geom2)

### Parameters
* `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type

Returns the distance between the geometries or -1 if an error occurs.
"""
distance(g1::AbstractGeometry, g2::AbstractGeometry) = OGR_G_Distance(g1.ptr, g2.ptr)
distance(D1, D2) = helper_geoms_run_fun(distance, D1, D2, false)

# ---------------------------------------------------------------------------------------------------
"""
    geomarea(geom)

### Parameters
* `geom`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

Returns the area of the geometry or 0.0 for unsupported geometry types.
"""
geomarea(geom::AbstractGeometry) = OGR_G_Area(geom.ptr)
geomarea(D) = helper_geoms_run_fun(geomarea, D, false)

# ---------------------------------------------------------------------------------------------------
"""
    geomlength(geom)

### Parameters
* `geom`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

Returns the length of the geometry, or 0.0 for unsupported geometry types.
"""
geomlength(geom::AbstractGeometry) = OGR_G_Length(geom.ptr)
geomlength(D) = helper_geoms_run_fun(geomlength, D, false)

# ---------------------------------------------------------------------------------------------------
"""
    envelope(geom)

### Parameters
* `geom`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

Computes and returns the bounding envelope for this geometry.
"""
function envelope(geom::AbstractGeometry)
    envelope = Ref{OGREnvelope}(OGREnvelope(0, 0, 0, 0))
    OGR_G_GetEnvelope(geom.ptr, envelope)
    return envelope[]
end
envelope(D) = helper_geoms_run_fun(envelope, D, false)

# ---------------------------------------------------------------------------------------------------
"""
    envelope3d(geom)

### Parameters
* `geom`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

Computes and returns the bounding envelope (3D) for this geometry
"""
function envelope3d(geom::AbstractGeometry)
    envelope = Ref{OGREnvelope3D}(OGREnvelope3D(0, 0, 0, 0, 0, 0))
    OGR_G_GetEnvelope3D(geom.ptr, envelope)
    return envelope[]
end
envelope3d(D) = helper_geoms_run_fun(envelope3d, D, false)

# ---------------------------------------------------------------------------------------------------
"""
    simplify(geom::AbstractGeometry, tol::Real)

Compute a simplified geometry.

### Parameters
* `geom`: the geometry.
* `tol`: the distance tolerance for the simplification.
"""
simplify(geom::AbstractGeometry, tol::Real) = IGeometry(OGR_G_Simplify(geom.ptr, tol))
function simplify(D, tol::Real; gdataset=false)
	ig = simplify(helper_1geom(D), tol)
	return (gdataset) ? ig : gd2gmt(ig)
end

# ---------------------------------------------------------------------------------------------------
"""
    delaunay(geom::AbstractGeometry, tol::Real, onlyedges::Bool)

### Parameters
* `geom`: the geometry.
* `tol`: optional snapping tolerance to use for improved robustness
* `onlyedges`: if `true`, will return a MULTILINESTRING, otherwise it
    will return a GEOMETRYCOLLECTION containing triangular POLYGONs.

Return a Delaunay triangulation of the vertices of the geometry.
"""
delaunaytriangulation(geom::AbstractGeometry, tol::Real, edges::Bool) = IGeometry(OGR_G_DelaunayTriangulation(geom.ptr, tol, edges))
function delaunaytriangulation(D, tol::Real, edges::Bool; gdataset=false)
	geom = helper_1geom(D)
	ig = delaunaytriangulation(geom, tol, edges)
	return (gdataset) ? ig : gd2gmt(ig)
end

# ---------------------------------------------------------------------------------------------------
"""
    intersects(geom1, geom2)

### Parameters
* `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type

Returns whether the geometries intersect

Determines whether two geometries intersect. If GEOS is enabled, then this is done in rigorous fashion
otherwise `true` is returned if the envelopes (bounding boxes) of the two geometries overlap.
"""
intersects(g1::AbstractGeometry, g2::AbstractGeometry) = Bool(OGR_G_Intersects(g1.ptr, g2.ptr))
function intersects(D1, D2)::Bool
	g1, g2 = helper_2geoms(D1, D2)
	intersects(g1, g2)
end

# ---------------------------------------------------------------------------------------------------
"""
    equals(geom1, geom2)

### Parameters
* `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type

Returns `true` if the geometries are equivalent.
"""
equals(g1::AbstractGeometry, g2::AbstractGeometry) = Bool(OGR_G_Equals(g1.ptr, g2.ptr))
equals(D1, D2) = helper_geoms_run_fun(equals, D1, D2, false)

# ---------------------------------------------------------------------------------------------------
"""
    disjoint(geom1, geom2)

### Parameters
* `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type

Returns `true` if the geometries are disjoint.
"""
disjoint(g1::AbstractGeometry, g2::AbstractGeometry) = Bool(OGR_G_Disjoint(g1.ptr, g2.ptr))
disjoint(D1, D2) = helper_geoms_run_fun(disjoint, D1, D2, false)

# ---------------------------------------------------------------------------------------------------
"""
    touches(geom1, geom2)

### Parameters
* `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type

Returns `true` if the geometries are touching.
"""
touches(g1::AbstractGeometry, g2::AbstractGeometry) = Bool(OGR_G_Touches(g1.ptr, g2.ptr))
touches(D1, D2) = helper_geoms_run_fun(touches, D1, D2, false)

# ---------------------------------------------------------------------------------------------------
"""
    crosses(geom1, geom2)

### Parameters
* `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type

Returns `true` if the geometries are crossing.
"""
crosses(g1::AbstractGeometry, g2::AbstractGeometry) = Bool(OGR_G_Crosses(g1.ptr, g2.ptr))
crosses(D1, D2) = helper_geoms_run_fun(crosses, D1, D2, false)

# ---------------------------------------------------------------------------------------------------
"""
    within(geom1, geom2)

### Parameters
* `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type

Returns `true` if g1 is contained within g2.
"""
within(g1::AbstractGeometry, g2::AbstractGeometry) = Bool(OGR_G_Within(g1.ptr, g2.ptr))
within(D1, D2) = helper_geoms_run_fun(within, D1, D2, false)

# ---------------------------------------------------------------------------------------------------
"""
    contains(geom1, geom2)

### Parameters
* `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type

Returns `true` if g1 contains g2.
"""
contains(g1::AbstractGeometry, g2::AbstractGeometry) = Bool(OGR_G_Contains(g1.ptr, g2.ptr))
contains(D1, D2) = helper_geoms_run_fun(contains, D1, D2, false)

# ---------------------------------------------------------------------------------------------------
"""
    overlaps(geom1, geom2)

### Parameters
* `geom1`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `geom2`: Second geometry. AbstractGeometry if `geom1::AbstractGeometry` or a GMTdataset/Matrix if `geom1` is GMT type

Returns `true` if the geometries overlap.
"""
overlaps(g1::AbstractGeometry, g2::AbstractGeometry) = Bool(OGR_G_Overlaps(g1.ptr, g2.ptr))
overlaps(D1, D2) = helper_geoms_run_fun(overlaps, D1, D2, false)

# ---------------------------------------------------------------------------------------------------
"""
    boundary(geom; gdataset=false)

### Parameters
* `geom`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `gdataset`: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

A new geometry object is created and returned containing the boundary of the
geometry on which the method is invoked.

### Returns
A GMT dataset when input is a Matrix or a GMT type (except if `gdaset=true`), or a GDAL IGeometry otherwise
"""
boundary(geom::AbstractGeometry) = IGeometry(OGR_G_Boundary(geom.ptr))
boundary(D; gdataset=false) = helper_geoms_run_fun(boundary, D; gdataset=gdataset)

# ---------------------------------------------------------------------------------------------------
"""
    convexhull(geom; gdataset=false)

### Parameters
* `geom`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `gdataset`: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

A new geometry object is created and returned containing the convex hull of the geometry on which the method is invoked.

### Returns
A GMT dataset when input is a Matrix or a GMT type (except if `gdaset=true`), or a GDAL IGeometry otherwise
"""
convexhull(geom::AbstractGeometry) = IGeometry(OGR_G_ConvexHull(geom.ptr))
convexhull(D; gdataset=false) = helper_geoms_run_fun(convexhull, D; gdataset=gdataset)

# ---------------------------------------------------------------------------------------------------
"""
    pointalongline(geom, distance::Real; gdataset=false)

### Parameters
* `geom`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `distance`: distance along the curve at which to sample position. This
    distance should be between zero and geomlength() for this curve.
* `gdataset`: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

Fetch point (or NULL) at given distance along curve.

### Returns
A GMT dataset when input is a Matrix or a GMT type (except if `gdaset=true`), or a GDAL IGeometry otherwise
"""
pointalongline(geom::AbstractGeometry, distance::Real) = IGeometry(OGR_G_Value(geom.ptr, distance))
function pointalongline(D, distance::Real; gdataset=false)
	geom = helper_1geom(D)
	ig = pointalongline(geom, distance)
	return (gdataset) ? ig : gd2gmt(ig)
end

# ---------------------------------------------------------------------------------------------------
"""
    polygonize(geom; gdataset=false)

### Parameters
* `geom`: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix
* `gdataset`: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

Polygonizes a set of sparse edges.

A new geometry object is created and returned containing a collection of
reassembled Polygons: NULL will be returned if the input collection doesn't
correspond to a MultiLinestring, or when reassembling Edges into Polygons is
impossible due to topological inconsistencies.

### Returns
A GMT dataset when input is a Matrix or a GMT type (except if `gdaset=true`), or a GDAL IGeometry otherwise
"""
polygonize(geom::AbstractGeometry) = IGeometry(OGR_G_Polygonize(geom.ptr))
polygonize(D; gdataset=false) = helper_geoms_run_fun(polygonize, D; gdataset=gdataset)

# ---------------------------------------------------------------------------------------------------
function helper_geoms_run_fun(f::Function, D1, D2, retds::Bool=true; gdataset=false)
	# Helper function that checks if geoms intersect and then run the F function
	# If GDATASET is true, return the GDAL dataset
	g1, g2 = helper_2geoms(D1, D2)
	ig = f(g1, g2)
	(gdataset) && (retds = false)
	return (retds) ? gd2gmt(ig) : ig
end
function helper_geoms_run_fun(f::Function, D, retds::Bool=true; gdataset=false)
	geom = helper_1geom(D)
	ig = f(geom)
	(gdataset) && (retds = false)
	return (retds) ? gd2gmt(ig) : ig
end

function helper_2geoms(D1, D2)
	# Helpr function that deals with arg checking and data type conversions, This is common to several funs here
	(!isa(D1, Matrix{<:Real}) && !isa(D2, Matrix{<:Real}) && !isa(D1, GMT.GMTdataset) && !isa(D1, Vector{<:GMT.GMTdataset}) && !isa(D2, GMT.GMTdataset) && !isa(D2, Vector{<:GMT.GMTdataset})) && error("Input mut be GMTdatset or Matrix{Real}")
	ds1 = gmt2gd(D1)
	ds2 = gmt2gd(D2)
	g1  = getgeom(unsafe_getfeature(getlayer(ds1, 0),0))
	g2  = getgeom(unsafe_getfeature(getlayer(ds2, 0),0))
	return g1, g2
end

function helper_1geom(D)
	(!isa(D, Matrix{<:Real}) && !isa(D, Matrix{<:Real}) && !isa(D, GMT.GMTdataset) && !isa(D, Vector{<:GMT.GMTdataset}) ) && error("Input mut be GMTdatset or Matrix{Real}")
	ds = gmt2gd(D)
	getgeom(unsafe_getfeature(getlayer(ds, 0),0))
end


text:
macro asserteq(x, y)
  xs = sprint(print, x)
  ys = sprint(print, y)

  quote
    xv = $(esc(x))
    yv = $(esc(y))
    if xv ≠ yv
      msg = sprint(print, "(", $xs, " = ", xv, ") ≠ (", $ys, " = ", yv, ")")
      throw(AssertionError(msg))
    end
  end
end

symbol_n(cstr::Union{Vector{UInt8}, SubArray{UInt8, 1}}, len) =
  ccall(:jl_symbol_n, Symbol, (Ptr{UInt8}, Int), cstr, len)

struct BSONElem
  tag::BSONType
  pos::Int
end

BSONElem(tag::BSONType, io::IO) = BSONElem(tag, position(io))

mutable struct ParseCtx
  refindx::Vector{BSONElem}
  refs::Vector{Any}
  curref::Int32
end

ParseCtx() = ParseCtx([], [], -1)

struct Document{IOT <: IO, DET}
  io::IOT
  ctx::ParseCtx
  index::Dict{Symbol, Int}
  elems::Vector{Pair{Symbol, BSONElem}}
end

function Document(io::IO, ::Type{DET}, ctx::ParseCtx) where DET
  len = read(io, Int32)
  index = Dict{Symbol, Int}()
  elems = Pair{Symbol, BSONElem}[]

  while (tag = read(io, BSONType)) ≠ eof
    name = parse_cstr_unsafe(io)
    nsym = symbol_n(name, length(name))

    push!(elems, nsym => BSONElem(tag, position(io)))
    index[nsym] = length(elems)

    skip_over(io, tag)
  end

  Document{typeof(io), DET}(io, ctx, index, elems)
end

function Document(io::IO, ::Type{DET}) where DET
  doc = Document(io, DET, ParseCtx())

  build_refs_indx!(io, doc, doc.ctx)
  doc
end
Document(io::IO) = Document(io, Any)

Base.getindex(doc::Document{IOT, DET}, key::Symbol) where {IOT, DET} = doc[DET, key]
function Base.getindex(doc::Document, ::Type{ET}, key::Symbol)::ET where ET
  elem = doc.elems[doc.index[key]].second
  seek(doc.io, elem.pos)
  parse_specific(doc.io, ET, elem.tag, doc.ctx)
end

Base.iterate(doc::Document) = iterate(doc, 1)
function Base.iterate(doc::Document{IOT, DET}, i::Int) where {IOT, DET}
  if i > length(doc.elems)
    return nothing
  end

  (name, elem) = doc.elems[i]
  if name == :_backrefs
    iterate(doc, i + 1)
  else
    seek(doc.io, elem.pos)
    (name => parse_specific(doc.io, DET, elem.tag, doc.ctx), i + 1)
  end
end

Base.length(doc::Document) = length(doc.elems)

struct ParseArrayIter{IOT <: IO, ET}
  io::IOT
  ctx::ParseCtx
  len::Int
end

ParseArrayIter(io::IOT, ::Type{ET}, len, ctx) where {IOT, ET} =
  ParseArrayIter{IOT, ET}(io, ctx, len)

ParseArrayIter(io::IOT, ::Type{ET}, ctx) where {IOT, ET} =
  ParseArrayIter{IOT, ET}(io, ctx, Int(read(io, Int32)) - 4)

Base.iterate(itr::ParseArrayIter) = iterate(itr, 0)

function Base.iterate(itr::ParseArrayIter{IOT, E},
                      len::Int)::Union{Nothing, Tuple{E, Int}} where {IOT, E}
  startpos = position(itr.io)
  tag = read(itr.io, BSONType)
  if tag == eof
    @asserteq (len + sizeof(BSONType)) itr.len
    return nothing
  end

  while read(itr.io, UInt8) != 0x00 end

  obj = parse_specific(itr.io, E, tag, itr.ctx)::E
  (obj, len + position(itr.io) - startpos)
end

parse_array_len(io::IO, ctx::ParseCtx) = Int(read(io, Int32)) - 4

function parse_array_tag(io::IO, ctx::ParseCtx)::BSONType
  tag = read(io, BSONType)
  tag == eof && return eof

  while read(io, UInt8) ≠ 0x00 end

  tag
end

Base.isempty(itr::ParseArrayIter) = itr.len == 1

function skip_over(io::IO, tag::BSONType)
  len = if tag == document || tag == array
    read(io, Int32) - 4
  elseif tag == string
    read(io, Int32)
  elseif tag == binary
    read(io, Int32) + 1
  elseif tag == null
    0
  else
    sizeof(jtype(tag))
  end

  seek(io, position(io) + len)
#  @info "Skipped" tag len position(io)
end

"Create an index into the _backrefs entry in the root document"
function build_refs_indx!(io::IO, doc::Document, ctx::ParseCtx)
  i = get(doc.index, :_backrefs, nothing)
  i == nothing && return
  elem = doc.elems[i].second

  if elem.tag != array
    error("_backrefs is not an array; tag = $(elem.tag)")
  end

  seek(io, elem.pos)
  len = read(io, Int32)
  #@info "Processing _backrefs" position(io) len

  while (tag = read(io, BSONType)) ≠ eof
    while read(io, UInt8) ≠ 0x00 end

    push!(ctx.refindx, BSONElem(tag, position(io)))
    push!(ctx.refs, nothing)
    skip_over(io, tag)
  end
end

function setref(obj, ctx::ParseCtx)
  if ctx.curref ≠ -1
    @asserteq ctx.refs[ctx.curref] nothing
    ctx.refs[ctx.curref] = obj
    ctx.curref = -1
  end
end

function parse_bin(io::IO, ctx::ParseCtx)::Vector{UInt8}
  len = read(io, Int32)
  subtype = read(io, 1)
  bin = read(io, len)
  setref(bin, ctx)
  bin
end

parse_bin_unsafe(io, ctx) = parse_bin(io, ctx)
function parse_bin_unsafe(io::IOBuffer, ctx::ParseCtx)::SubArray{UInt8, 1}
  len = read(io, Int32)
  subtype = read(io, 1)
  bin = view(io.data, (io.ptr):(io.ptr + len - 1))
  setref(bin, ctx)
  bin
end

function parse_tag(io::IO, tag::BSONType, ctx::ParseCtx)
  if tag == null
    @asserteq ctx.curref -1
  elseif tag == document
    parse_doc(io, ctx)
  elseif tag == array
    parse_any_array(io, ctx)
  elseif tag == string
    @asserteq ctx.curref -1
    len = read(io, Int32)-1
    s = String(read(io, len))
    eof = read(io, 1)
    s
  elseif tag == binary
    parse_bin(io, ctx)
  else
    @asserteq ctx.curref -1
    read(io, jtype(tag))
  end
end

function parse_symbol(io::IO, name::BSONElem, ctx::ParseCtx)
  @asserteq name.tag string
  seek(io, name.pos)
  len = read(io, Int32)-1
  cstr = parse_cstr_unsafe(io)
  @asserteq len length(cstr)
  symbol_n(cstr, len)
end

for (T, expected) in (Nothing => null,
                      Bool => boolean,
                      Int32 => int32,
                      Int64 => int64,
                      Float64 => double)
  @eval begin
    function parse_specific(io::IO, ::Type{$T}, tag::BSONType,
                            ctx::ParseCtx)::$T
      @asserteq tag $expected
      read(io, $T)
    end
  end
end

parse_specific(io::IO, ::Type{Any}, tag::BSONType, ctx::ParseCtx) =
  parse_tag(io, tag, ctx)

function parse_specific(io::IO, ::Type{String},
                        tag::BSONType, ctx::ParseCtx)::String
  @asserteq tag string
  len = read(io, Int32)-1
  s = String(read(io, len))
  eof = read(io, 1)
  #@info "Parse specific string" s
  s
end

function parse_specific(io::IO, ::Type{Symbol},
                        tag::BSONType, ctx::ParseCtx)::Symbol
  @asserteq tag document

  len = read(io, Int32)
  local name

  for _ in 1:3
    if (tag = read(io, BSONType)) == eof
      break
    end

    k = parse_cstr_unsafe(io)
    if k == b"tag"
      @asserteq tag string
    elseif k == b"name"
      name = BSONElem(tag, io)
    else
      error("Expected tag or name, but got '$(String(k))'")
    end

    skip_over(io, tag)
  end
  endpos = position(io)

  ret = parse_symbol(io, name, ctx)
  seek(io, endpos)
  ret
end

function parse_specific(io::IO, ::Type{Vector{UInt8}},
                        tag::BSONType, ctx::ParseCtx)::Vector{UInt8}
  len = read(io, Int32)
  subtype = read(io, 1)
  read(io, len)
end

function parse_specific(io::IO, ::Type{BSONArray},
                        tag::BSONType, ctx::ParseCtx)::BSONArray
  @asserteq tag array
  parse_any_array(io, ctx)
end

function parse_specific(io::IO, ::Type{BSONDict},
                        tag::BSONType, ctx::ParseCtx)::BSONDict
  @asserteq tag document
  parse_doc(io, Any, ctx)
end

function parse_specific(io::IO, ::Type{DataType},
                        tag::BSONType, ctx::ParseCtx)::DataType
  @asserteq tag document
  len = read(io, Int32)
  local name, params
  ref = nothing

  for _ in 1:4
    if (tag = read(io, BSONType)) == eof
      break
    end

    k = parse_cstr_unsafe(io)
    if k == b"tag"
      @asserteq tag string
    elseif k == b"name"
      name = BSONElem(tag, io)
    elseif k == b"params"
      params = BSONElem(tag, io)
    elseif k == b"ref"
      ref = BSONElem(tag, io)
    else
      error("Expected tag, name, ref or params, but got '$(String(k))'")
    end

    skip_over(io, tag)
  end
  endpos = position(io)

  T = if ref ≠ nothing
    parse_specific_ref(io, DataType, ref, ctx)::DataType
  else
    parse_type(io, name, params, ctx)
  end
  seek(io, endpos)
  T
end

function parse_specific(io::IO, ::Type{T}, tag::BSONType,
                        ctx::ParseCtx)::T where T
  @asserteq tag document
  #@info "Parse specfic" T

  if @generated
    if !isconcretetype(T)
      return :(parse_tag(io, tag, ctx))
    end

    quote
      startpos = position(io)
      len = read(io, Int32)
      local data::BSONElem
      ref = nothing

      for _ in 1:4
        if (tag = read(io, BSONType)) == eof
          break
        end

        k = parse_cstr_unsafe(io)
        if k == b"tag"
          @asserteq tag string
        elseif k == b"ref"
          ref = BSONElem(tag, io)
        elseif k == b"type"
          @asserteq tag document
        elseif k == b"data"
          @asserteq tag array
          data = BSONElem(tag, io)
        end

        skip_over(io, tag)
      end

      endpos = position(io)
      @asserteq (startpos + len) endpos

      if ref ≠ nothing
        ret = parse_specific_ref(io, $T, ref, ctx)::$T
        seek(io, endpos)
        return ret
      end

      seek(io, data.pos)
      ret = load_struct(io, $T, data.tag, ctx)
      seek(io, endpos)
      ret
    end
  else
    parse_tag(io::IO, tag, ctx)::T
  end
end

function parse_specific(io::IO, ::Type{Array{T, N}}, tag::BSONType,
                        ctx::ParseCtx)::Array{T, N} where {T, N}
  if @generated
    expr = quote
      startpos = position(io)
      len = read(io, Int32)
      local data::BSONElem, sizes

      if tag == array
        @asserteq $N 1
        return $T[ParseArrayIter(io, $T, len - 4, ctx)...]
      end
      @asserteq tag document

      for _ in 1:5
        if (tag = read(io, BSONType)) == eof
          break
        end

        k = parse_cstr_unsafe(io)
        if k == b"tag"
          @asserteq tag string
          skip_over(io, tag)
        elseif k == b"type"
          @asserteq tag document
          skip_over(io, tag)
        elseif k == b"size"
          @asserteq tag array
          sizes = (ParseArrayIter(io, Int64, ctx)...,)
        elseif k == b"data"
          @assert tag == array || tag == binary
          data = BSONElem(tag, io)
          skip_over(io, tag)
        end
      end

      endpos = position(io)
      @asserteq (startpos + len) endpos
      @asserteq length(sizes) $N
      seek(io, data.pos)
    end

    expr = if isbitstype(T)
      :($expr; ret = load_bits_array(io, $T, sizes, data, ctx))
    else
      :($expr; ret = load_array(io, $T, sizes, ctx))
    end

    :($expr;
      seek(io, endpos);
      ret)
  else
    parse_tag(io, tag, ctx)
  end
end

parse_specific(::IO, ::Type{Function}, ::BSONType, ::ParseCtx) =
  error("Functions are not supported, use load_compat")

function parse_type(io::IO,
                    name::BSONElem, params::BSONElem, ctx::ParseCtx)::DataType
  @asserteq name.tag array
  @asserteq params.tag array
  curref = ctx.curref
  ctx.curref = -1

  seek(io, name.pos)
  T = resolve(ParseArrayIter(io, String, ctx))
  #@info "Union all type" T
  seek(io, params.pos)
  ctx.curref = curref
  p = constructtype(T, ParseArrayIter(io, Any, ctx))
  setref(p, ctx)
  p
end

function parse_any_array(io::IO, ctx::ParseCtx)
  len = read(io, Int32)

  # If this array has a reference, then it might reference itself.
  if ctx.curref ≠ -1
    ps = BSONArray()
    setref(ps, ctx)

    while (tag = parse_array_tag(io, ctx)) ≠ eof
      push!(ps, parse_tag(io, tag, ctx))
    end

    ps
  else
    # If it is not referenced then we can recreate the array with a wider type
    # as we discover new types (see base/array.jl)
    tag = parse_array_tag(io, ctx)
    tag ≠ eof || return BSONArray()
    ps = [parse_tag(io, tag, ctx)]

    parse_any_array!(io::IO, ps, ctx)
  end
end

function parse_any_array!(io::IO, ps::Vector{T}, ctx::ParseCtx) where T
  while (tag = parse_array_tag(io, ctx)) ≠ eof
    e = parse_tag(io, tag, ctx)

    if e isa T || typeof(e) === T
      push!(ps, e::T)
    else
      new = sizehint!(empty(ps, Base.typejoin(T, typeof(e))), length(ps))

      append!(new, ps)
      push!(new, e)
      return parse_any_array!(io, new, ctx)
    end
  end

  ps
end

function load_bits_array(io::IO, ::Type{T}, sizes,
                         data::BSONElem, ctx::ParseCtx) where T
  arr = if sizeof(T) == 0
      fill(T(), sizes...)
    else
      @asserteq data.tag binary
      reshape(reinterpret_(T, parse_bin_unsafe(io, ctx)), sizes...)
    end

    setref(arr, ctx)
    arr
end

function load_array(io::IO, ::Type{T}, sizes, ctx::ParseCtx) where T
    arr = Array{T}(undef, sizes...)
    setref(arr, ctx)
    bsonarr = ParseArrayIter(io, T, ctx)

    if (itr = iterate(bsonarr)) ≠ nothing
      for (i, len) = enumerate(sizes), j = 1:len
        (elem, s) = itr
        arr[i, j] = elem
        itr = iterate(bsonarr, s)
      end
    end

    arr
end

function parse_array(io::IO, ttype::BSONElem, size::BSONElem,
                     data::BSONElem, ctx::ParseCtx)::AbstractArray
  # Save current ref incase T is a backref
  curref = ctx.curref
  ctx.curref = -1

  seek(io, ttype.pos)
  T = parse_specific(io, DataType, ttype.tag, ctx)::DataType
  #@info "New array type" T

  ctx.curref = curref

  seek(io, size.pos)
  sizes = (ParseArrayIter(io, Int64, ctx)...,)

  seek(io, data.pos)
  if isbitstype(T)
    load_bits_array(io, T, sizes, data, ctx)
  else
    load_array(io, T, sizes, ctx)
  end
end

function parse_specific_ref(io::IO, ::Type{T}, ref::BSONElem,
                            ctx::ParseCtx)::T where T
  seek(io, ref.pos)
  id = if ref.tag == int64
    convert(Int32, read(io, Int64))
  elseif ref.tag == int32
    read(io, Int32)
  else
    error("Expecting Int type found: $ref_tag")
  end

  ret = ctx.refs[id]::Union{Nothing, T}
  if ret ≠ nothing
    ret
  else
    ctx.curref = id
    obj = ctx.refindx[id]
    seek(io, obj.pos)
    ctx.refs[id] = parse_specific(io, T, obj.tag, ctx)::T
  end
end

function parse_backref(io::IO, ref::BSONElem, ctx::ParseCtx)
  parse_specific_ref(io, Any, ref, ctx)
end

function load_dict!(io::IO, d::Dict{K, V},
                    ctx::ParseCtx) where {K, V}
  setref(d, ctx)

  parse_array_len(io, ctx)
  tag = parse_array_tag(io, ctx)
  ks = parse_specific(io, Vector{K}, tag, ctx)::Vector{K}

  tag = parse_array_tag(io, ctx)
  vs = parse_specific(io, Vector{V}, tag, ctx)::Vector{V}

  for (k, v) in zip(ks, vs)
      d[k] = v
  end

  d
end

function load_dict!(io::IO, d::Dict{K, Nothing},
                    ctx::ParseCtx) where K
  setref(d, ctx)

  parse_array_len(io, ctx)
  tag = parse_array_tag(io, ctx)

  for k in parse_specific(io, Vector{K}, tag, ctx)::Vector{K}
      d[k] = nothing
  end

  d
end

function load_struct(io::IO, ::Type{T}, dtag::BSONType, ctx::ParseCtx)::T where T
  #@info "Load struct" T

  if @generated
    n = fieldcount(T)

    if isprimitive(T)
      @assert isbitstype(T)
      quote
        @asserteq dtag binary
        bits = parse_bin_unsafe(io, ctx)
        ccall(:jl_new_bits, Any, (Any, Ptr{Cvoid}), $T, bits)
      end
    elseif T <: Dict
      :(load_dict!(io, $T(), ctx))
    elseif n < 1
      :($T())
    elseif !ismutable(T)
      block = :(x = Vector{Any}(undef, $n);
                parse_array_len(io, ctx))

      for i in 1:n
        FT = fieldtype(T, i)
        block = :($block;
                  tag = parse_array_tag(io, ctx);
                  x[$i] = parse_specific(io, $FT, tag, ctx)::$FT)
      end

      :($block;
        ccall(:jl_new_structv, Any, (Any,Ptr{Cvoid},UInt32), $T, x, $n))
    else
      block = :(x = ccall(:jl_new_struct_uninit, Any, (Any,), $T);
                setref(x, ctx);
                parse_array_len(io, ctx))

      for i in 1:n
        FT = fieldtype(T, i)
        block = :($block;
                  tag = parse_array_tag(io, ctx);
                  f = parse_specific(io, $FT, tag, ctx)::$FT;
                  ccall(:jl_set_nth_field, Nothing, (Any, Csize_t, Any), x, $i-1, f))
      end

      :($block; x)
    end
  else
    n = fieldcount(T)

    if isprimitive(T)
      @assert isbitstype(T)
      @asserteq dtag binary
      bits = parse_bin_unsafe(io, ctx)
      ccall(:jl_new_bits, Any, (Any, Ptr{Cvoid}), T, bits)
    elseif T <: Dict
      load_dict!(io, T(), ctx)
    elseif n < 1
      T()
    elseif !ismutable(T)
      x = Vector{Any}(undef, n)
      parse_array_len(io, ctx)

      for i in 1:n
        FT = fieldtype(T, i)
        tag = parse_array_tag(io, ctx)
        x[i] = parse_specific(io, FT, tag, ctx)::FT
      end

      ccall(:jl_new_structv, Any, (Any,Ptr{Cvoid},UInt32), T, x, n)
    else
      @asserteq dtag array

      x = ccall(:jl_new_struct_uninit, Any, (Any,), T)
      setref(x, ctx)

      for i in 1:n
        tag = parse_array_tag(io, ctx)
        FT = fieldtype(T, i)
        f = parse_specific(io, FT, tag, ctx)::FT
        ccall(:jl_set_nth_field, Nothing, (Any, Csize_t, Any), x, i-1, f)
      end

      x
    end
  end
end

function parse_struct(io::IO, ttype::BSONElem, data::BSONElem, ctx::ParseCtx)
  # Save current ref incase T is a backref
  curref = ctx.curref
  ctx.curref = -1

  seek(io, ttype.pos)
  T = parse_specific(io, DataType, ttype.tag, ctx)::DataType
  #@info "New struct type" T

  seek(io, data.pos)
  ctx.curref = curref

  load_struct(io, T, data.tag, ctx)
end

function parse_doc(io::IO, ::Type{V},
                   ctx::ParseCtx)::Dict{Symbol, V} where V
  len = read(io, Int32)
  dic = Dict{Symbol, V}()
  setref(dic, ctx)

  while (tag = read(io, BSONType)) ≠ eof
    cstr = parse_cstr_unsafe(io)
    k = symbol_n(cstr, length(cstr))
    dic[k] = parse_specific(io, V, tag, ctx)
  end

  dic
end

function parse_doc(io::IO, ctx::ParseCtx)
  startpos = position(io)
  len = read(io, Int32)

  seen::Int64 = 0
  see(it::Int64) = seen = seen | it
  saw(it::Int64)::Bool = seen & it != 0
  only_saw(it::Int64)::Bool = seen == it

  # First decide if this document is tagged with a Julia type. Saving the BSON tag types
  local tref, tdata, ttype, ttypename, ttag, tname, tparams, tpath,
        tsize, tvar, tbody
  local k::AbstractVector{UInt8}

  for _ in 1:6
    if (tag = read(io, BSONType)) == eof
      break
    end
    k = parse_cstr_unsafe(io)
    #@info "Read key" String(k)

    if k == b"tag"
      if tag == string && (dtag = parse_doc_tag(io)) isa Int64
        #@info "Read tag" dtag
        see(SEEN_TAG | dtag)
        continue
      else
        @goto FALLBACK
      end
    end

    if k == b"ref"
      see(SEEN_REF)
      tref = BSONElem(tag, io)
    elseif k == b"data"
      see(SEEN_DATA)
      tdata = BSONElem(tag, io)
    elseif k == b"type"
      see(SEEN_TYPE)
      ttype = BSONElem(tag, io)
    elseif k == b"typename"
      see(SEEN_TYPENAME)
      ttypename = BSONElem(tag, io)
    elseif k == b"name"
      see(SEEN_NAME)
      tname = BSONElem(tag, io)
    elseif k == b"params"
      see(SEEN_PARAMS)
      tparams = BSONElem(tag, io)
    elseif k == b"path"
      see(SEEN_PATH)
      tpath = BSONElem(tag, io)
    elseif k == b"size"
      see(SEEN_SIZE)
      tsize = BSONElem(tag, io)
    elseif k == b"var"
      see(SEEN_VAR)
      tvar = BSONElem(tag, io)
    elseif k == b"body"
      see(SEEN_BODY)
      tbody = BSONElem(tag, io)
    elseif k == b"_backrefs"
      nothing
    else
      @goto FALLBACK
    end

    skip_over(io, tag)
  end
  endpos = position(io)

  ret = if only_saw(SEEN_TAG | SEEN_REF | SEEN_TAG_BACKREF)
    #@info "Found backref" tref
    parse_backref(io, tref, ctx)
  elseif only_saw(SEEN_TAG | SEEN_TYPE | SEEN_DATA | SEEN_TAG_STRUCT)
    #@info "Found Struct" ttype tdata
    parse_struct(io, ttype, tdata, ctx)
  elseif only_saw(SEEN_TAG | SEEN_NAME | SEEN_PARAMS | SEEN_TAG_DATATYPE)
    #@info "Found Type" tname tparams
    parse_type(io, tname, tparams, ctx)
  elseif only_saw(SEEN_TAG | SEEN_NAME | SEEN_TAG_SYMBOL)
    #@info "Found Symbol" tname
    parse_symbol(io, tname, ctx)
  elseif only_saw(SEEN_TAG | SEEN_DATA | SEEN_TAG_TUPLE)
    #@info "Found Tuple" tdata
    seek(io, tdata.pos)
    (ParseArrayIter(io, Any, ctx)...,)
  elseif only_saw(SEEN_TAG | SEEN_DATA | SEEN_TAG_SVEC)
    #@info "Found svec" tdata
    seek(io, tdata.pos)
    Core.svec(ParseArrayIter(io, Any, ctx)...)
  elseif only_saw(SEEN_TAG | SEEN_TAG_UNION)
    Union{}
  elseif only_saw(SEEN_TAG | SEEN_TYPENAME | SEEN_PARAMS | SEEN_TAG_ANON)
    error("Functions are not supported, use load_compat")
  elseif only_saw(SEEN_TAG | SEEN_PATH | SEEN_TAG_REF)
    error("Refs not implemented")
  elseif only_saw(SEEN_TAG | SEEN_TYPE | SEEN_SIZE | SEEN_DATA | SEEN_TAG_ARRAY)
    parse_array(io, ttype, tsize, tdata, ctx)
  elseif only_saw(SEEN_TAG | SEEN_VAR | SEEN_BODY | SEEN_TAG_UNIONALL)
    error("Unionall not implemented")
  else
    @goto FALLBACK
  end

  seek(io, endpos)
  return ret

  @label FALLBACK
  #@info "Found plain dictionary"
  # This doc doesn't appear to be tagged with all the necessay julia type info
  seek(io, startpos)
  parse_doc(io, Any, ctx)
end

function parse(io::IO, ctx::ParseCtx)
  doc = Document(io, Any, ctx)
  build_refs_indx!(io, doc, ctx)

  # The root document is usually a plain BSONDict, but it can also be a Julia
  # struct or some other type if the user figures out how to cause that
  if haskey(doc.index, :tag)
    seek(io, 0)
    parse_doc(io, ctx)
  else
    BSONDict(doc)
  end
end

parse(path::String, ctx::ParseCtx; mmap=false, opts...) = open(path) do io
  if mmap
    parse(IOBuffer(Mmap.mmap(io; shared=false)), ctx; opts...)
  else
    parse(io, ctx; opts...)
  end
end

load(x; opts...) = parse(x, ParseCtx(); opts...)

function directtrip(ting::T) where {T}
  io = IOBuffer()
  bson(io, Dict(:stuff => ting))
  seek(io, 0)
  doc = Document(io, T)
  try
    doc[:stuff]::T
  catch e
    @error "Error during parsing" io doc.ctx
    rethrow(e)
  end
end

ext: jl


text:
# cuda21: Kernels for Array->Array reduction.

fp = open("cuda22.cu","w")
#using Knet: reduction_ops

function cuda22src(f, j, op, f1, v0; BLK=128, THR=128) # BLK not used, determined by ny
    sprint() do s
        for (T,F) in [("float","$(f)_32"),("double","$(f)_64")]
            print(s,
"""
__device__ void _$(F)_22_0(volatile $T *x, int i) {
//for optimizing warps, volatile must be used as register optimization will lead to wrong answers
  $T ai, xi;
  ai=x[i]; xi=x[i+32]; x[i]=$op;
  ai=x[i]; xi=x[i+16]; x[i]=$op;
  ai=x[i]; xi=x[i+ 8]; x[i]=$op;
  ai=x[i]; xi=x[i+ 4]; x[i]=$op;
  ai=x[i]; xi=x[i+ 2]; x[i]=$op;
  ai=x[i]; xi=x[i+ 1]; x[i]=$op;
}

__global__ void _$(F)_22(int nx, int xd1, $T *x, int s1, int s2, int ny, $T *y) {
  // x[i] goes into y[(i/sy)%ny]
  __shared__ $T buffer[$THR];
  int t = threadIdx.x;
  int b = blockIdx.x;
  int d = blockDim.x;
  $T ai, xi;

  for (int bi = b; bi < ny; bi += blockDim.x) {
    // sum the elements assigned to this thread
    ai = $v0;
    int lower = (bi%s1) + (bi/s1)*s2;
    int upper = lower + s1*xd1;
    for (int i=lower+t*s1; i <= upper; i+=d*s1) {
      xi = x[i]; xi=$f1; ai=$op;
    }
    buffer[t] = ai;
    __syncthreads();

    // help sum the entries in the block
    for(int stride=$THR/2; stride>32; stride>>=1) {
      if(t < stride) {
        ai=buffer[t]; xi=buffer[stride+t]; buffer[t]=$op;
      }
      __syncthreads();   // Q: can this be outside the for loop?
    }

    if(t<32) {
      _$(F)_22_0(buffer,t);  // This reuses warpSum from 20 scalar reduction.
    }
    __syncthreads();

    if(t==0) {  // the first thread in the block writes the block result to y
      y[bi]=buffer[0];
    }
  }
}

extern "C" {
$DLLEXPORT void $(F)_22(int nx, int xd1, $T *x, int s1, int s2, int ny, $T *y) {
  _$(F)_22<<<$BLK,$THR>>>(nx,xd1,x,s1,s2,ny,y);
}
$DLLEXPORT void $(F)_22_stream(int nx, int xd1, $T *x, int s1, int s2, int ny, $T *y, cudaStream_t STR) {
  _$(F)_22<<<$BLK,$THR,0,STR>>>(nx,xd1,x,s1,s2,ny,y);
}
}

""")
        end
    end
end

for a in reduction_ops
    if !isa(a,Tuple); a=(a,); end
    print(fp,cuda22src(a...))
end

close(fp)


text:
function _spawn_runme(gwt::GWTask)

    runme = _runme_file(gwt)
    chmod(rscfile, 0o755)
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
# Pentaquark decay chain: Lb -> Jψ p K
#   in the s-channel scattering kinematics
#   Lb Jψ -> p K
#
# it is check for an example described in arXiv:1805.02113
# https://arxiv.org/pdf/1805.02113.pdf
# at the page 29

using ThreeBodyDecay

mLb = 5.62; mJψ = 3.09; mp=0.938; mK = 0.49367
const tbs = ThreeBodySystem(mJψ,mp,mKmLb; two_jps = ([2,1,0,1],['-','+','-','+']))

# helicity couplings
# Lb Jψ
H1(two_λ0,two_λ1;two_LS=(4,3),two_s=3) = CG_doublearg(1,two_λ0,2,-two_λ1,two_LS[2],two_λ0-two_λ1) *
                                         CG_doublearg(two_LS[1],0,two_LS[2],two_λ0-two_λ1,two_s,two_λ0-two_λ1)
# pK
H2(two_λ2;two_l=4,two_s=3) = CG_doublearg(1,two_λ2,0,0,1,two_λ2) *
                             CG_doublearg(two_l,0,1,two_λ2,two_s,two_λ2)

# amplitude
function A(z,two_λs;two_s=3,two_LS=(4,3),two_l=4)
    return H1(two_λs[4],two_λs[1]; two_LS=two_LS, two_s=two_s) *
        H2(two_λs[2]; two_l=two_l, two_s=two_s) *
        # (two_λs[1]!=0 ? 1.0 : 2.0) * # weight of A(λ_J/ψ = 0) could be different with E/m factor
        wignerd_doublearg(two_s,two_λs[4]-two_λs[1],two_λs[2],z)
end

# amplitude squared
Asq_summed(z;two_s=3,two_LS=(4,3),two_l=4) =
    sum(abs2(A(z,two_λs;two_s=two_s,two_LS=two_LS,two_l=two_l))
        for two_λs in Iterators.product(-tbs.two_js[1]:2:tbs.two_js[1],
                                        -tbs.two_js[2]:2:tbs.two_js[2],
                                        -tbs.two_js[3]:2:tbs.two_js[3],
                                        -tbs.two_js[4]:2:tbs.two_js[4]))
I(σs) = Asq_summed(cosθ23(σs,tbs.msq))

## plotting
using Plots
# S=3/2, D-wave
plot(z->Asq_summed(z),-1,1, lab="", ylim=(0,1))
# P-wave
plot(z->Asq_summed(z;two_s=3,two_LS=(2,3),two_l=4),-1,1, lab="", ylim=(0,1))
# s = 1/2, LS=(0,1/2), l=0, angle dependence is constant again
plot(z->Asq_summed(z;two_s=1,two_LS=(0,1),two_l=0)+1e-5*z,-1,1, lab="", ylim=(0,3))
#

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
using Base.Iterators: Iterators, product
using ..Datasets: getsavefolder, Cube

"""
    savecube(cube,name::String)

Save a [`YAXArray`](@ref) to the folder `name` in the YAXArray working directory.
"""
function savecube(
    c,
    name::AbstractString;
    chunksize = Dict(),
    max_cache = 1e8,
    backend = :zarr,
    backendargs...,
)
    allax = caxes(c)
    if !isa(chunksize, AbstractDict) && !isa(chunksize, NamedTuple)
        @warn "Chunksize must be provided as a Dict mapping axis names to chunk size in the future"
        chunksize = OrderedDict(axname(i[1]) => i[2] for i in zip(allax, chunksize))
    end
    firstaxes = sort!([findAxis(String(k), allax) for k in keys(chunksize)])
    lastaxes = setdiff(1:length(allax), firstaxes)
    allax = allax[[firstaxes; lastaxes]]
    dl = cumprod(length.(caxes(c)))
    isplit = findfirst(i -> i > max_cache / sizeof(eltype(c)), dl)
    isplit isa Nothing && (isplit = length(dl) + 1)
    forcesingle = (isplit + 1) < length(firstaxes)
    axin = allax[1:isplit-1]
    axn = axname.(axin)
    indims = InDims(axn...)
    outdims = OutDims(
        axn...,
        backend = backend,
        chunksize = chunksize,
        path = name;
        backendargs...,
    )
    function cop(xout, xin)
        xout .= xin
    end
    if forcesingle
        nprocs() > 1 && println("Forcing single core processing because of bad chunk size")
        o = mapCube(
            cop,
            c,
            indims = indims,
            outdims = outdims,
            ispar = false,
            max_cache = max_cache,
        )
    else
        o = mapCube(cop, c, indims = indims, outdims = outdims, max_cache = max_cache)
    end
end

function loadcube(s)
    Cube(getsavefolder(s, true))
end

function rmcube(s)
    p = getsavefolder(s, true)
    if isfile(p) || isdir(p)
        rm(p, recursive = true)
    end
end

function writefun(
    xout,
    xin::AbstractArray{Union{Missing,T}},
    a,
    nd,
    cont_loop,
    filename,
) where {T}

    #TODO replace -9999.0 this with a default missing value
    x = map(ix -> ismissing(ix) ? convert(T, -9999.0) : ix, xin)

    count_vec = fill(-1, nd)
    start_vec = fill(1, nd)

    used_syms = Symbol[]
    for (k, v) in cont_loop
        count_vec[k] = 1
        start_vec[k] = a[Symbol(v)][1]
        push!(used_syms, Symbol(v))
    end

    splinds = Iterators.filter(i -> !in(i, used_syms), keys(a))
    vn = join(string.([a[s][2] for s in splinds]), "_")
    isempty(vn) && (vn = "layer")
    ncwrite(x, filename, vn, start = start_vec, count = count_vec)
end

global const projection = Dict(
    "grid_mapping_name" => "transverse_mercator",
    "longitude_of_central_meridian" => -9.0,
    "false_easting" => 500000.0,
    "false_northing" => 0.0,
    "latitude_of_projection_origin" => 0.0,
    "scale_factor_at_central_meridian" => 0.9996,
    "long_name" => "CRS definition",
    "longitude_of_prime_meridian" => 0.0,
    "semi_major_axis" => 6378137.0,
    "inverse_flattening" => 298.257223563,
    "spatial_ref" => "PROJCS[\"WGS 84 / UTM zone 29N\",GEOGCS[\"WGS 84\",
    DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,
    AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],
    PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],
    UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],
    AUTHORITY[\"EPSG\",\"4326\"]],
    PROJECTION[\"Transverse_Mercator\"],
    PARAMETER[\"latitude_of_origin\",0],
    PARAMETER[\"central_meridian\",-9],
    PARAMETER[\"scale_factor\",0.9996],
    PARAMETER[\"false_easting\",500000],
    PARAMETER[\"false_northing\",0],
    UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],
    AXIS[\"Easting\",EAST],
    AXIS[\"Northing\",NORTH],
    AUTHORITY[\"EPSG\",\"32629\"]]",
    "GeoTransform" => "722576.2320803495 20 0 4115483.464603715 0 -20 ",
)

global const epsg4326 = Dict(
    "grid_mapping_name" => "latitude_longitude",
    "long_name" => "CRS definition",
    "longitude_of_prime_meridian" => 0.0,
    "semi_major_axis" => 6378137.0,
    "inverse_flattening" => 298.257223563,
    "spatial_ref" =>
        "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]",
    "GeoTransform" =>
        "-98.68712696894481 0.0001796630568239077 0 20.69179551612753 0 -0.0001796630568239077 ",
)


text:
using LinearAlgebra
eye(n) = Matrix{Float64}(I, n, n)

function gaussnewton(x0, f, tol=1e-6, maxit = 100)
    alpha = 1.0e-4
    it = 1
    xc = x0
    fc, jac, gc = f(xc)
    n_f = 1
    n_g = 1
    n_h = 0
    while norm(gc) > tol && it < maxit
        dc = (jac'*jac)\gc
        lambda = 1.0
        xt = xc - lambda*dc
        ft, = f(xt)
        n_f = n_f
        iarm = 0
        it += 1
        # goal for sufficient decrease
        fgoal = fc - alpha*lambda*(gc'*dc)
        while ft > fgoal
            iarm += 1
            lambda = lambda/2
            fgoal = fc - alpha*lambda*(gc'*dc)
            xt = xc - lambda*dc
            ft, = f(xt)
            n_f += 1
            if iarm > 10
                println("Armijo error in Gauss-Newton")
            end
        end
        xc = xt
        fc, jac, gc = f(xc)
        n_f += 1
        n_g += 1
    end
    return xc
end

function levenbergmarquardt(x0, f, tol=1e-6, maxit=100)
    it = 1
    xc = x0
    fc, jac, gc, rc = f(xc)
    nu0 = 0.001
    m = length(gc)
    n = length(xc)
    nu = norm(gc)
    while norm(gc) > tol && it < maxit
        it += 1
        hc = [jac; sqrt(nu).*eye(n)]
        rhs = [rc; zeros(n)]
        dc = -hc\rhs
        xt = xc + dc
        xp, nup, id = trustregiontest(f, xc, xt, fc, jac, gc, nu, nu0, rc)
        if id > 30
            error("To many iteraion in trust region solver")
        end
        xc = xp
        nu = nup
        fc, jac, gc, rc = f(xc)
    end
    return xc
end

function trustregiontest(f, xc, xt, fc, jac, gc, nu, nu0, rc)
    mu0 = 0.0
    mulow = 0.25
    muhigh = 0.75
    mup = 2.0
    mdown = 0.5
    z = xc
    n = length(z)
    it = 0
    while z == xc && it < 30
        ft,  = f(xt)
        it += 1
        s = xt - xc
        diff = fc - ft
        pred = - (gc'*s)/2.0
        rat = diff/pred
        if rat < mu0
            nu = max(nu*mup, nu0)
            hc = [jac; sqrt(nu).*eye(n)]
            rhs = [rc; zeros(n)]
            dc = -hc\rhs
            xt = xc + dx
        elseif rat < mulow
            z = xt
            nu = max(nu*mup, nu0)
        else
            z = xt
            if rat > muhigh
                nu = mdown * nu
                if nu < nu0
                    nu = 0.0
                end
            end
        end
    end
    return z, nu, it
end


text:
#==============
1D Earth models
==============#
"""
Abstract type for radially-symmetric models of the Earth.  1D Earth models should
be subtypes of this.
"""
abstract type EarthModel1D <: EarthModel end

const model_variables_EarthModel1D = (:vp, :vs, :density, :vph, :vpv, :vsh, :vsv,
    :eta, :Qμ, :Qκ)
const model_names_EarthModel1D = ("isotropic average Vp", "isotropic average Vs",
    "density", "Vph", "Vpv", "Vsh", "Vsv", "η", "Qμ", "Qκ")
const model_units_EarthModel1D = (" (km/s)", " (km/s)", " (g/cm^3)", " (km/s)",
    " (km/s)", " (km/s)", " (km/s)", "", "", "")

"""
    findlayer(m::EarthModel1D, r) -> layer_index

Return the layer index in which radius `r` km is located in the EarthModel1D `m`.
"""
function findlayer(m::EarthModel1D, r::Real)
    r < 0. && throw(DomainError(r, "radius cannot be negative"))
    r > m.a && throw(DomainError(r, "radius is greater than Earth radius for model ($(m.a) km)"))
    for i in 1:length(m.r)
        r < m.r[i] && return i
    end
    length(m.r)
end

"""
    PREMPolyModel <: EarthModel1D

Type describing the Earth as a set of layers within which properties vary
according to a set of polynomials.

Physical parameters are represented by arrays of size `(order+1, n)`, where `n`
is the number of layers, and `order` is the order of polynomial which is
used to represent the parameter.  Hence a constant layer has size `(1,n)`,
and to compute the value of `x` in layer `i`, for an Earth radius of `a` km,
at a radius of `r` km, the expression is:

    val_x = x[i,1] + (r/a)*x[i,2] + (r/a)^2*x[i,3] ... (r/a)^order*x[i,order+1]
"""
struct PREMPolyModel <: EarthModel1D
    "Earth radius in km"
    a :: Float64
    "Number of distinct layers"
    n :: Int
    "Radii in km of top of each layer"
    r :: Vector{Float64}
    "Physical parameters"
    vp :: Array{Float64}
    vs :: Array{Float64}
    density :: Array{Float64}
    "Optional anisotropic parameters"
    aniso :: Bool
    vph :: Array{Float64}
    vpv :: Array{Float64}
    vsh :: Array{Float64}
    vsv :: Array{Float64}
    eta :: Array{Float64}
    "Optional attenuation parameters"
    attenuation :: Bool
    Qμ :: Array{Float64}
    Qκ :: Array{Float64}
end

# Evaluation routines--all documented here
# TODO: Use Horner's method à la Base.@evalpoly
for (sym, name, unit) in zip(model_variables_EarthModel1D, model_names_EarthModel1D, model_units_EarthModel1D)
    @eval begin
        """
            $(split(string($sym), ".")[end])(m::EarthModel1D, r; depth=false) -> $($name)

        Return the value of $($name)$($unit) for model `m` at radius `r` km.

        If `depth` is `true`, then `r` is given as a depth in km instead.
        """
        function ($sym)(m::PREMPolyModel, r::Real; depth::Bool=false)
            length(m.$sym) > 0 || throw(ArgumentError("$($name) not defined for model"))
            depth && (r = radius(m, r))
            ir = findlayer(m, r)
            x = r/m.a
            val = m.$(sym)[1,ir]
            for k in 2:size(m.$sym, 1)
                val = val + m.$(sym)[k,ir]*x^(k-1)
            end
            val
        end
    end
end
# Alternative names
const ρ = density
const Qmu = Qμ
const Qkappa = Qκ

"""
    evaluate(m::EarthModel1D, field::Symbol, r) -> vals

Evaluate the model `m` at radius `r` km for the different property/ies in `field`,
returning a scalar for scalar input, and an array for array input.
"""
function evaluate(m::PREMPolyModel, field::Symbol, r; depth::Bool=false)
    y = getfield(m, field)
    length(y) > 0 || throw(ArgumentError("'$field' not defined for model"))
    depth && (r = radius(m, r))
    ir = findlayer(m, r)
    x = r/m.a
    val = y[1,ir]
    for k in 2:size(y, 1)
        val = val + y[k,ir]*x^(k-1)
    end
    val
end

"""
    SteppedLayeredModel <: EarthModel1D

A `SteppedLayeredModel` contains `n` layers with maximum radius `r` km,
each with a constant velocity.
"""
struct SteppedLayeredModel <: EarthModel1D
    "Earth radius in km"
    a :: Float64
    "Number of layers"
    n :: Int
    "Radii in km of top of each layer"
    r :: Vector{Float64}
    "Physical parameters"
    vp :: Vector{Float64}
    vs :: Vector{Float64}
    density :: Vector{Float64}
    "Optional anisotropic parameters"
    aniso :: Bool
    vph :: Vector{Float64}
    vpv :: Vector{Float64}
    vsh :: Vector{Float64}
    vsv :: Vector{Float64}
    eta :: Vector{Float64}
    "Optional attenuation parameters"
    attenuation :: Bool
    Qμ :: Vector{Float64}
    Qκ :: Vector{Float64}
end

for sym in model_variables_EarthModel1D
    @eval function ($sym)(m::SteppedLayeredModel, r::Real; depth::Bool=false)
        length(m.$sym) > 0 ||
            throw(ArgumentError("'$split(string($sym), ".")[end]' not defined for model"))
        depth && (r = radius(m, r))
        m.$(sym)[findlayer(m, r)]
    end
end

function evaluate(m::SteppedLayeredModel, field::Symbol, r; depth::Bool=false)
    length(getfield(m, field)) > 0 ||
        throw(ArgumentError("$field' not defined for model"))
    depth && (r = radius(m, r))
    getfield(m, field)[findlayer(m, r)]
end

"""
    LinearLayeredModel <: EarthModel1D

A `LinearLayeredModel` contains `n` points at which velocities are defined, with
linear interpolation between them.  Hence there are `n - 1` layers.

Discontinuities are represented by two layers with the same radii.
"""
struct LinearLayeredModel <: EarthModel1D
    "Earth radius in km"
    a :: Float64
    "Number of layers"
    n :: Int
    "Radii in km of top of each layer"
    r :: Vector{Float64}
    "Physical parameters"
    vp :: Vector{Float64}
    vs :: Vector{Float64}
    density :: Vector{Float64}
    "Optional anisotropic parameters"
    aniso :: Bool
    vph :: Vector{Float64}
    vpv :: Vector{Float64}
    vsh :: Vector{Float64}
    vsv :: Vector{Float64}
    eta :: Vector{Float64}
    "Optional attenuation parameters"
    attenuation :: Bool
    Qμ :: Vector{Float64}
    Qκ :: Vector{Float64}
end

function findlayer(m::LinearLayeredModel, r::Real)
    r < 0. && error("Radius cannot be negative")
    r > m.a && error("Radius $r km is greater than Earth radius for model ($(m.a) km)")
    for i in 1:(length(m.r) - 1)
        r < m.r[i+1] && return i
    end
    return length(m.r) - 1
end

for sym in model_variables_EarthModel1D
    @eval function ($sym)(m::LinearLayeredModel, r::Real; depth::Bool=false)
        length(m.$sym) > 0 ||
            throw(ArgumentError("'$(split(string($sym), ".")[end])' not defined for model"))
        depth && (r = radius(m, r))
        ir = findlayer(m, r)
        r0 = m.r[ir]
        r1 = m.r[ir+1]
        (m.$sym[ir]*(r1 - r) + m.$sym[ir+1]*(r - r0))/(r1 - r0)
    end
end

function evaluate(m::LinearLayeredModel, field::Symbol, r; depth::Bool=false)
    y = getfield(m, field)
    length(y) > 0 ||
        throw(ArgumentError("'$field' not defined for model"))
    depth && (r = radius(m, r))
    ir = findlayer(m, r)
    r0 = m.r[ir]
    r1 = m.r[ir+1]
    (y[ir]*(r1 - r) + y[ir+1]*(r - r0))/(r1 - r0)
end

## Derived quantities
"""
    bulk_modulus(m::EarthModel1D, r) -> K

Return the bulk modulus `K` in GPa at radius `r` in the model `m`.
"""
bulk_modulus(m::EarthModel1D, r) = density(m, r)*vp(m, r)^2 - 4/3*shear_modulus(m, r)

const NewtonG = 6.67428e-11
"""
    gravity(m::EarthModel1D, r) -> g

Return the acceleration due to gravity, `g`, in m/s^2 at radius `r` km.
"""
gravity(m::EarthModel1D, r) = (r == 0) ? 0. : NewtonG*mass(m, r)/(r*1.e3)^2

"""
    mass(m::EarthModel1D, r) -> mass

Return the mass in kg between the centre of the model and the radius `r` km.
"""
function mass(m::PREMPolyModel, r)
	l = findlayer(m, r)
	r *= 1.e3 # SI
	M = 0.
	for i = 1:l
		if i == 1
			R0 = 0.
		else
			R0 = m.r[i-1]*1.e3
		end
		if i < l
			R = m.r[i]*1.e3
		else
			R = r
		end
        for k in 1:size(m.density, 1)
            rho = m.density[k,i]/m.a^(k-1)*1.e3^(2-k)/(k+2)
            M += rho*(R^(k+2) - R0^(k+2))
        end
	end
	4*π*M
end

function mass(m::SteppedLayeredModel, r)
    l = findlayer(m, r)
    r *= 1.e3
    M = 0.
    for i in 1:l
        R0 = i == 1 ? 0. : m.r[i-1]*1.e3
        R = i < l ? m.r[i]*1.e3 : r
        M += m.density[i]*(R^3 - R0^3)
    end
    4/3*π*M
end

function mass(m::LinearLayeredModel, r)
    l = findlayer(m, r)
    r *= 1.e3
    M = 0.
    for i in 1:l
        m.r[i] == m.r[i+1] && continue
        R0 = m.r[i]*1.e3
        R = i == l ? r : m.r[i+1]*1.e3
        dρ_dr = 1.e3*(m.density[i+1] - m.density[i])/(m.r[i+1]*1.e3 - R0)
        ρ0 = 1.e3*m.density[i] - dρ_dr*R0
        M += ρ0*(R^3 - R0^3)/3 + dρ_dr*(R^4 - R0^4)/4
    end
    4*π*M
end

"""
    pressure(m::EarthModel1D, r) -> p

Return the pressure `p` in Pa at radius `r` km.
"""
function pressure(m::EarthModel1D, r)
    f(r) = pressure_integration_func(m, r)
    1.e3*quadgk(f, r, m.a)[1]
end
pressure_integration_func(m::EarthModel1D, r) = 1.e3*density(m, r)*gravity(m, r)

"""
    poissons_ratio(m, r) -> ν

Return the Poisson's ratio for the model `m` given a radius `r` in km.
The calculation uses the isotropic average velocities at `r`.
"""
function poissons_ratio(m::EarthModel1D, r)
    G = shear_modulus(m, r)
    K = bulk_modulus(m, r)
    (3K - 2G)/(6K + 2G)
end

"""
    shear_modulus(m::EarthModel1D, r) -> μ

Return the shear modulus `μ` (often also called G) in GPa at radius `r` in the
model `m`.
"""
shear_modulus(m::EarthModel1D, r) = vs(m, r)^2*density(m, r)

"""
    surface_mass(m::EarthModel1D, r) -> mass

Return the mass in kg betwen radius `r` km and the surface.
"""
surface_mass(m::EarthModel1D, r) = mass(m, m.a) - mass(m, r)

"""
    youngs_modulus(m, r) -> E

Return the Young's modulus `E` in GPa for the model `m` given a radius `r` in km.
"""
youngs_modulus(m::EarthModel1D, r) = 2*shear_modulus(m, r)*(1 + poissons_ratio(m, r))

"""
    moment_of_inertia(m, r0=0, r1=surface_radius(m)) -> I

Return the moment of interia `I` in kg m² for the model `m` between radii `r0`
and `r1` in km.
"""
moment_of_inertia(m::EarthModel1D, r0=0, r1=surface_radius(m)) =
    8/3*π*quadgk(r->1e3*density(m, r/1e3)*r^4, r0*1e3, r1*1e3)[1]


text:
@testset "quat basic tests" begin

    @testset "type conversion" for T in (Int32, Int64, Float32, Float64, BigFloat, BigInt)
        @test real(quat(T)) == T
        @test complex(quat(T)) == complex(T)
        @test widen(Quaternion{T}) == Quaternion{widen(T)}
        @test quat(Complex{T}) == Quaternion{T}
        @test quat(Quaternion{T}) == Quaternion{T}
        @test float(Quaternion{T}(1,1,1,1)) == quat(float(T(1)), float(T(1)), float(T(1)), float(T(1)))
        @test big(Quaternion{T}) == Quaternion{big(T)}
    end

    @test promote_rule(QuatF16, Float32) == QuatF32
    @test promote_rule(QuatF64, ComplexF32) == QuatF64
    @test promote_rule(Quaternion{BigFloat}, QuatF64) == Quaternion{BigFloat}

    @test big(quat(1,1,1,1)) == quat(big(1), big(1), big(1), big(1))

    @test Quaternion{Float64}(true) == quat(1.0)

    @test real(quat(1.0)) == 1.0
    @test real(quat(1.0+2.0im)) == 1.0
    @test convert(Quaternion, true) == quat(true)
    @test convert(Quaternion{Float64}, true) == Quaternion{Float64}(1.0)
    @test convert(Quaternion, 1.0) == quat(1.0)
    @test convert(Quaternion, 1.0+2.0im) == quat(1.0+2.0im)
    @test convert(Quaternion, 1+im+jm+km) == quat(1,1,1,1)
    @test convert(Quaternion{Float32}, 1.0) == quat(1.0)
    @test convert(Quaternion{Float32}, 1.0+2.0im) == quat(1.0+2.0im)
    @test convert(Quaternion{Int32}, 1+im+jm+km) == quat(1,1,1,1)
    @test imag(quat(1.0+2.0im)) == 2.0
    @test jmag(1.0+2.0im+3.0jm+4.0km) == 3.0
    @test kmag(1.0+2.0im+3.0jm+4.0km) == 4.0
    @test quat(vec([0.0 2.0 3.0 4.0])) == 2im + 3jm + 4km
    @test isinteger(quat(1,0,0,0))
    @test !isinteger(quat(1.1, 0.0, 0.0, 0.0))
    @test !isinteger(jm)
    @test isfinite(quat(1.0, 10.0, 100.0, 1000.0))
    @test !isfinite(quat(1.0, Inf, Inf, Inf))
    @test isnan(quat(NaN, NaN, NaN, NaN))
    @test !isnan(1 + im + jm + km)
    @test isinf(quat(1.0, Inf, Inf, Inf))
    @test !isinf(quat(1.0, 10.0, 100.0, 1000.0))
    @test iszero(quat(0.0, 0.0, 0.0, 0.0))
    @test iszero(quat(0, 0, 0, 0))
    @test !iszero(im + jm + km)
    @test quat(jm+km) == jm+km
    @test complex(1+im+jm+km) == 1+im
    @test vec(1+2im+3jm+4km) == [1; 2; 3; 4]

    @test conj(1+2im+3jm+4km) == 1-2im-3jm-4km

    @test quat(false,true,false,false) == im
    @test quat(false,false,true,false) == jm
    @test quat(false,false,false,true) == km
    @test quat(1,0,0,0) == 1
    @test quat(0,1,0,0) == im
    @test quat(0,0,1,0) == jm
    @test quat(0,0,0,1) == km
    @test quat(1.0,0.0,0.0,0.0) == 1
    @test quat(0.0,1.0,0.0,0.0) == im
    @test quat(0.0,0.0,1.0,0.0) == jm
    @test quat(0.0,0.0,0.0,1.0) == km
    @test quat(-1,0,0,0) == -1
    @test quat(0,-1,0,0) == -im
    @test quat(0,0,-1,0) == -jm
    @test quat(0,0,0,-1) == -km
    @test quat(-1.0,0.0,0.0,0.0) == -1
    @test quat(0.0,-1.0,0.0,0.0) == -im
    @test quat(0.0,0.0,-1.0,0.0) == -jm
    @test quat(0.0,0.0,0.0,-1.0) == -km
    @test im*jm ==  km == -jm*im
    @test jm*im == -km == -im*jm
    @test jm*km ==  im == -km*jm
    @test km*jm == -im == -jm*km
    @test km*im ==  jm == -im*km
    @test im*km == -jm == -km*im
    @test im*im == jm*jm == km*km == -1

    @test true + jm == quat(1,0,1,0)
    @test km + true == quat(1,0,0,1)
    @test true - jm == quat(1,0,-1,0)
    @test km - true == quat(-1,0,0,1)
    @test true * jm == jm
    @test jm * true == jm
    @test km * false == zero(Quaternion{Bool})
    @test jm + im == quat(false, true, true, false)
    @test im + km == quat(false, true, false, true)

    @test complex(true, true) - quat(true, true, false, false) == 0
    @test quat(true, true, false, false) - complex(true, true) == 0

    @test true + quat(0.0,1.0,1.0,1.0) == quat(1.0, 1.0, 1.0, 1.0)
    @test quat(0.0,1.0,1.0,1.0) + true == quat(1.0, 1.0, 1.0, 1.0)
    @test true - quat(0.0,1.0,1.0,1.0) == conj(quat(1.0, 1.0, 1.0, 1.0))
    @test quat(0.0,1.0,1.0,1.0) - true == quat(-1.0, 1.0, 1.0, 1.0)
    @test 1.0 + jm == quat(1.0, 0.0, 1.0, 0.0)
    @test jm + 1.0 == quat(1.0, 0.0, 1.0, 0.0)
    @test 1.0 - jm == quat(1.0, 0.0, -1.0, 0.0)
    @test jm - 1.0 == quat(-1.0, 0.0, 1.0, 0.0)
    @test jm * 1.0 == quat(0.0, 0.0, 1.0, 0.0)
    @test 1.0 * jm == quat(0.0, 0.0, 1.0, 0.0)

    @test 1.0 - quat(0.0, -1.0, -1.0, -1.0) == quat(1.0, 1.0, 1.0, 1.0)

    q, c, r = rand(QuatF64), rand(ComplexF64), rand(Float64)
    @test q * true == q
    @test false * q == quat(0.0)
    @test r + q == quat(r + q.re, q.im, q.jm, q.km)
    @test q + r == quat(q.re + r, q.im, q.jm, q.km)
    @test c + q == quat(c.re + q.re, c.im + q.im, q.jm, q.km)
    @test q + c == quat(q.re + c.re, q.im + c.im, q.jm, q.km)
    @test c - q == quat(c.re - q.re, c.im - q.im, -q.jm, -q.km)
    @test q - c == quat(q.re - c.re, q.im - c.im, q.jm, q.km)

    @test complex(1, 1) * quat(1, 1, 1, 1) == quat(0, 2, 0, 2)
    @test quat(1, 1, 1, 1) * complex(1, 1) == quat(0, 2, 2, 0)

    @test isone(quat(1,0,0,0))
    @test isone(quat(1.0))
    @test !isone(jm)

    @test bswap(quat(1,1,1,1)) == quat(bswap(1), bswap(1), bswap(1), bswap(1))

    @test quat(complex(1,2), complex(3,4)) == quat(1,2,3,4)
    @test quat(quat(1,2,3,4)) == quat(1,2,3,4)
    @test quat(jm) == jm
    @test Quaternion(km) == km
end

# These mirror tests from Complex.jl
@testset "arithmetic" begin
    @testset for T in (Float16, Float32, Float64, BigFloat) # BigFloat was failing
        t = true
        f = false
        u = quat(T(+1.0), T(+1.0), T(+1.0), T(+1.0))

        @testset "add and subtract" begin

            @test isequal(T(+0.0) + im + jm + km, quat(T(+0.0),T(+1.0),T(+1.0),T(+1.0)))
            @test isequal(T(-0.0) + im + jm + km, quat(T(-0.0),T(+1.0),T(+1.0),T(+1.0)))
            @test isequal(T(+0.0) - im - jm - km, quat(T(+0.0),T(-1.0),T(-1.0),T(-1.0)))
            @test isequal(T(-0.0) - im - jm - km, quat(T(-0.0),T(-1.0),T(-1.0),T(-1.0)))
            @test isequal(T(+1.0) + im + jm + km, quat(T(+1.0),T(+1.0),T(+1.0),T(+1.0)))
            @test isequal(T(-1.0) + im + jm + km, quat(T(-1.0),T(+1.0),T(+1.0),T(+1.0)))
            @test isequal(T(+1.0) - im - jm - km, quat(T(+1.0),T(-1.0),T(-1.0),T(-1.0)))
            @test isequal(T(-1.0) - im - jm - km, quat(T(-1.0),T(-1.0),T(-1.0),T(-1.0)))
            @test isequal(jm + T(+0.0), quat(T(+0.0),T(+0.0),T(+1.0),T(+0.0)))
            @test isequal(jm + T(-0.0), quat(T(-0.0),T(+0.0),T(+1.0),T(+0.0)))
            @test isequal(jm - T(+0.0), quat(T(+0.0),T(+0.0),T(+1.0),T(+0.0)))
            @test isequal(jm - T(-0.0), quat(T(+0.0),T(+0.0),T(+1.0),T(+0.0)))
            @test isequal(km + T(+0.0), quat(T(+0.0),T(+0.0),T(+0.0),T(+1.0)))
            @test isequal(km + T(-0.0), quat(T(-0.0),T(+0.0),T(+0.0),T(+1.0)))
            @test isequal(km - T(+0.0), quat(T(+0.0),T(+0.0),T(+0.0),T(+1.0)))
            @test isequal(km - T(-0.0), quat(T(+0.0),T(+0.0),T(+0.0),T(+1.0)))
            @test isequal(im + jm + km + T(+1.0), quat(T(+1.0),T(+1.0),T(+1.0),T(+1.0)))
            @test isequal(im + jm + km + T(-1.0), quat(T(-1.0),T(+1.0),T(+1.0),T(+1.0)))
            @test isequal(im + jm + km - T(+1.0), quat(T(-1.0),T(+1.0),T(+1.0),T(+1.0)))
            @test isequal(im + jm + km - T(-1.0), quat(T(+1.0),T(+1.0),T(+1.0),T(+1.0)))
            @test T(f) + im + jm + km == quat(T(+0.0),T(+1.0),T(+1.0),T(+1.0))
            @test T(t) + im + jm + km == quat(T(+1.0),T(+1.0),T(+1.0),T(+1.0))
            @test T(f) - im - jm - km == quat(T(+0.0),T(-1.0),T(-1.0),T(-1.0))
            @test T(t) - im - jm - km == quat(T(+1.0),T(-1.0),T(-1.0),T(-1.0))
            @test im + jm + km + T(f) == quat(T(+0.0),T(+1.0),T(+1.0),T(+1.0))
            @test im + jm + km + T(t) == quat(T(+1.0),T(+1.0),T(+1.0),T(+1.0))
            @test im + jm + km - T(f) == quat(T(+0.0),T(+1.0),T(+1.0),T(+1.0))
            @test im + jm + km - T(t) == quat(T(-1.0),T(+1.0),T(+1.0),T(+1.0))
        end

        @testset "multiply" begin
            oq  = quat(T(+1.0),T(+0.0),T(+0.0),T(+0.0))
            iq  = quat(T(+0.0),T(+1.0),T(+0.0),T(+0.0))
            jq  = quat(T(+0.0),T(+0.0),T(+1.0),T(+0.0))
            kq  = quat(T(+0.0),T(+0.0),T(+0.0),T(+1.0))
            oqn = quat(T(-1.0),T(-0.0),T(-0.0),T(-0.0))
            iqn = quat(T(-0.0),T(-1.0),T(-0.0),T(-0.0))
            jqn = quat(T(-0.0),T(-0.0),T(-1.0),T(-0.0))
            jkn = quat(T(-0.0),T(-0.0),T(-0.0),T(-1.0))
            @test isequal(T(+1.0) * u,  u)
            @test isequal(u * T(+1.0),  u)
            @test isequal(T(-1.0) * u, -u)
            @test isequal(u * T(-1.0), -u)
            @test isequal(T(+0.0) * u, quat(T(+0.0),T(+0.0),T(+0.0),T(+0.0)))
            @test isequal(u * T(+0.0), quat(T(+0.0),T(+0.0),T(+0.0),T(+0.0)))
            @test isequal(T(-0.0) * u, quat(T(-0.0),T(-0.0),T(-0.0),T(-0.0)))
            @test isequal(u * T(-0.0), quat(T(-0.0),T(-0.0),T(-0.0),T(-0.0)))
            @test isequal(oq * oq, oq)
            @test isequal(iq * iq, quat(T(-1.0),T(+0.0),T(+0.0),T(+0.0)))
            @test isequal(jq * jq, quat(T(-1.0),T(+0.0),T(+0.0),T(+0.0)))
            @test isequal(kq * kq, quat(T(-1.0),T(+0.0),T(+0.0),T(+0.0)))
            @test iq*jq ==  kq == -jq*iq
            @test jq*iq == -kq == -iq*jq
            @test jq*kq ==  iq == -kq*jq
            @test kq*jq == -iq == -jq*kq
            @test kq*iq ==  jq == -iq*kq
            @test iq*kq == -jq == -kq*iq
        end

        @testset "divide" begin
            @test u/T(+1.0) == u
            @test u/quat(T(+1.0)) == u
            @test T(+1.0)/u == 0.25conj(u)
        end

        @testset "rounding" begin
            for i=1:10
                q = rand(Quaternion{T})
                @test round(q) == quat(round(q.re), round(q.im), round(q.jm), round(q.km))
                @test round(q, RoundNearest, RoundToZero) == quat(round(q.re, RoundNearest),
                                                                  round(q.im, RoundToZero),
                                                                  round(q.jm, RoundToZero),
                                                                  round(q.km, RoundToZero))
                @test float(q) == q
            end
        end
    end
end

@testset "printing" begin
    @test sprint(show,  1 + 2im + 3jm + 4km) == "1 + 2im + 3jm + 4km"
    @test sprint(show, -1 - 2im - 3jm - 4km) == "-1 - 2im - 3jm - 4km"
    @test sprint(show,  1.0+2.0im+3.0jm+4.0km, context=:compact=>true) == "1.0+2.0im+3.0jm+4.0km"
    @test sprint(show, -1.0-2.0im-3.0jm-4.0km, context=:compact=>true) == "-1.0-2.0im-3.0jm-4.0km"
    @test sprint(show, NaN + NaN*im + NaN*jm + NaN*km) == "NaN + NaN*im + NaN*jm + NaN*km"
    @test sprint(show, quat(im)) == "im"
    @test sprint(show, jm) == "jm"
    @test sprint(show, km) == "km"
    @test sprint(show, Quaternion(true,true,true,true)) == "Quaternion(true,true,true,true)"

    @test sprint(show, "text/html",  1 + 2im + 3jm + 4km) == "1 + 2i + 3j + 4k"
    @test sprint(show, "text/html", -1 - 2im - 3jm - 4km) == "-1 - 2i - 3j - 4k"
    @test sprint(show, "text/html",  1.0+2.0im+3.0jm+4.0km, context=:compact=>true) == "1.0+2.0i+3.0j+4.0k"
    @test sprint(show, "text/html", -1.0-2.0im-3.0jm-4.0km, context=:compact=>true) == "-1.0-2.0i-3.0j-4.0k"
    @test sprint(show, "text/html", NaN + NaN*im + NaN*jm + NaN*km) == "NaN + NaN*i + NaN*j + NaN*k"
    @test sprint(show, "text/html", quat(im)) == "i"
    @test sprint(show, "text/html", jm) == "j"
    @test sprint(show, "text/html", km) == "k"
    @test sprint(show, "text/html", Quaternion(true,true,true,true)) == "Quaternion(true,true,true,true)"

    @testset "io" begin
        mktemp() do path, file
            q = quat(1,2,3,4)
            @test write(file, quat(1,2,3,4)) > 0
            seekstart(file)
            @test read(file, typeof(q)) == q
        end
    end
end

@testset "Matrix tests" begin

    Q = [1+2im+3jm+4jm -im+km; -jm+km -5+6im-7jm-8km]

    @test qmatrix(cmatrix(1+2im+3jm+4km))[1,1] == 1+2im+3jm+4km
    @test qmatrix(cmatrix(Q)) == Q
    @test quat(Q) == Q
    @test cmatrix(1+im+jm+km) == [1+im 1+im; im-1 1-im]
    @test cmatrix([1 1im; 1jm 1km]) == [1 im 0 0; 0 0 1 im; 0 0 1 (-im); (-1) im 0 0]

    for i=1:10
        Q = randn(QuatF64, 2, 2)
        W = randn(QuatF64, 2, 2)
        @test qmatrix(cmatrix(Q)) == Q
        @test cmatrix(Q + W) == cmatrix(Q) + cmatrix(W)
        @test cmatrix(Q * W) ≈ cmatrix(Q) * cmatrix(W)
    end

end

@testset "math" begin

    @test abs(quat(1,1,1,1)) == 2.0
    @test abs(quat(-1.0,2.0,-3.0,4.0)) == sqrt(1+4+9+16)
    @test abs2(quat(1,1,1,1)) == 4
    @test abs2(quat(-1.0,2.0,-3.0,4.0)) == 1.0+4.0+9.0+16.0

    @test inv(quat(1.0,0,0,0)) == quat(1.0,0,0,0)
    @test inv(quat(1,0,0,0)) == quat(1.0,0,0,0)
    @test inv(quat(0,1,0,0)) == quat(-1.0im)
    @test inv(jm) == -1.0jm
    @test inv(4km) == -0.25km

    @test norm(quat(1,2,3,4), 1) == 10.0
    @test norm(quat(1,2,3,4), 2) == sqrt(1+4+9+16)
    @test norm(quat(1,2,3,4), Inf) == 4.0
    @test normalize(quat(1,1,1,1), 1) == quat(1/4, 1/4, 1/4, 1/4)
    @test normalize(quat(-1,-1,-1,-1), 1) == -quat(1/4, 1/4, 1/4, 1/4)
    @test normalize(quat(1,1,1,1), 2) == quat(1/2, 1/2, 1/2, 1/2)
    @test normalize(quat(-8,-4,-4,-2), 2) == -quat(0.8, 0.4, 0.4, 0.2)

    @test exp(quat(1.0)) == quat(exp(1.0))
    @test exp(log(quat(1,1,1,1))) ≈ quat(1.0, 1.0, 1.0, 1.0)
    @test sqrt(quat(1.0)) ≈ quat(1.0)
    @test sqrt(quat(-1.0)) ≈ quat(1.0)
    # @test quat(im)^2.0 ≈ quat(jm)^2.0 ≈ quat(km)^2.0 ≈ quat(-1)

    for u in [1, im, jm, km]
        # @test exp(log(quat(u))) ≈ quat(u)
        @test log(exp(quat(u))) ≈ quat(u)
    end

    for i=1:100
        q = (rand(QuatF64) - quat(1/2,1/2,1/2,1/2))
        @test exp(log(q)) ≈ q
        @test log(exp(q)) ≈ q
        @test exp(log(quat(q.re))) ≈ quat(abs(q.re))
        @test log(exp(quat(q.re))) ≈ quat(q.re)
        @test exp(quat(q.re, q.im, 0, 0)) ≈ exp(complex(q.re, q.im))
        @test log(quat(q.re, q.im, 0, 0)) ≈ log(complex(q.re, q.im))
        # @test sqrt(q^2) ≈ q
        @test sqrt(q)^2.0 ≈ q
    end
end


text:
module Types

# types
export Arg, Option, Flag, EntryCommand, NodeCommand, LeafCommand, AbstractCommand

# interfaces
export print_cmd,
    cmd_doc, cmd_name, cmd_sym, set_brief_length!, max_brief_length, check_first_sentence_length

const MAX_DOC_WIDTH = 28
const INDENT = 2
const HELP_FLAG = "-h, --help"
const VERSION_FLAG = "-V, --version"
const HELP_FLAG_DOC = "print this help message"
const VERSION_FLAG_DOC = "print version information"
const PRESERVED = ["h", "help"]
const MAX_BRIEF_LENGTH = Ref{Int}(80)

"""
    set_brief_length!(n::Int=120)

Set the maximum length of a brief description. Default is 120.

!!! note

    This is only effective when called before [`@main`](@ref) or [`codegen`](@ref).
"""
function set_brief_length!(n::Int = 120)
    MAX_BRIEF_LENGTH[] = n
end

"""
    max_brief_length()

Get the maximum length allowed for a brief description.
"""
max_brief_length() = MAX_BRIEF_LENGTH[]

"""
    AbstractCommand

abstract type for commands.
"""
abstract type AbstractCommand end

struct CommandDoc
    first::String
    rest::String

    function CommandDoc(name::String, lineinfo::LineNumberNode, doc::String)
        doc = strip(doc)
        brief = check_first_sentence_length(name, lineinfo, doc)
        new(brief, doc[length(brief)+1:end])
    end

    function CommandDoc()
        new("", "")
    end
end

function Base.iterate(doc::CommandDoc, st = 1)
    if st == 1
        return doc.first, 2
    elseif st == 2
        return doc.rest, 3
    else
        return nothing
    end
end

function Base.endswith(doc::CommandDoc, b::String)
    if isempty(doc.rest)
        return endswith(doc.first, b)
    else
        return endswith(doc.rest, b)
    end
end

function Base.show(io::IO, doc::CommandDoc)
    print(io, "  ", doc.first, doc.rest)
end

Base.@kwdef struct Arg
    name::String = "arg"
    line::LineNumberNode = LineNumberNode(0)
    doc::CommandDoc = CommandDoc(name, line, "positional argument. ")
    require::Bool = true
    vararg::Bool = false
    type = Any
    default::Union{String,Nothing} = nothing
end

function Arg(
    name::String,
    line::LineNumberNode,
    doc::String,
    require::Bool,
    vararg::Bool,
    type,
    default,
)
    Arg(name, line, CommandDoc(name, line, doc), require, vararg, type, default)
end

Base.@kwdef struct Option
    name::String
    line::LineNumberNode = LineNumberNode(0)
    doc::CommandDoc = CommandDoc()
    arg::Arg = Arg()
    short::Bool = false
end

function Option(name::String, line::LineNumberNode, doc::String, arg::Arg, short::Bool)
    Option(name, line, CommandDoc(name, line, doc), arg, short)
end

Base.@kwdef struct Flag
    name::String
    line::LineNumberNode = LineNumberNode(0)
    doc::CommandDoc = CommandDoc()
    short::Bool = true
end

function Flag(name::String, line::LineNumberNode, doc::String, short::Bool)
    Flag(name, line, CommandDoc(name, line, doc), short)
end

"""
    EntryCommand <: AbstractCommand

`EntryCommand` describes the entry of the CLI. It contains
an actual root command (either [`LeafCommand`](@ref) or [`NodeCommand`](@ref))
of the entire CLI and a version number. The version number is `v"0.0.0"` by default.
"""
Base.@kwdef struct EntryCommand <: AbstractCommand
    root::Any
    version::VersionNumber = v"0.0.0"
    line::LineNumberNode = LineNumberNode(0)
end

"""
    NodeCommand <: AbstractCommand

`NodeCommand` describes the command in the middle of a CLI, e.g
in the following `remote` is a `NodeCommand`, it will dispatch
the call to its sub-command `show`. See also [`LeafCommand`](@ref).

```sh
git remote show origin
```
"""
Base.@kwdef struct NodeCommand <: AbstractCommand
    name::String
    line::LineNumberNode = LineNumberNode(0)
    doc::CommandDoc = CommandDoc()
    subcmds::Vector{Any}
end

function NodeCommand(name::String, line::LineNumberNode, doc::String, subcmds::Vector)
    NodeCommand(name, line, CommandDoc(name, line, doc), subcmds)
end

function NodeCommand(name::String, subcmds::Vector; kwargs...)
    NodeCommand(; name = name, subcmds = subcmds, kwargs...)
end

"""
    LeafCommand <: AbstractCommand

`LeafCommand` describes the command at the end of a CLI, e.g
in the following `show` is `LeafCommand`, it is the command
that actually executes things. See also [`NodeCommand`](@ref).

```sh
git remote show origin
```
"""
struct LeafCommand <: AbstractCommand
    entry::Any # a callable
    name::String
    args::Vector{Arg}
    nrequire::Int
    options::Vector{Option}
    flags::Vector{Flag}
    doc::CommandDoc
    line::LineNumberNode

    function LeafCommand(entry, name, args, options, flags, doc, line)
        check_duplicate_short_options(options, flags)
        check_varargs_position(args)
        nrequire = check_required_args(args)
        new(entry, name, args, nrequire, options, flags, doc, line)
    end
end

function LeafCommand(entry, name, args, options, flags, doc::String, line)
    LeafCommand(entry, name, args, options, flags, CommandDoc(name, line, doc), line)
end

Arg(name; kwargs...) = Arg(; name = name, kwargs...)
Option(name, arg = Arg(); kwargs...) = Option(; name = name, arg = arg, kwargs...)
Flag(name; kwargs...) = Flag(; name = name, kwargs...)
EntryCommand(root; kwargs...) = EntryCommand(; root = root, kwargs...)

function LeafCommand(
    entry;
    name::String = string(nameof(entry)),
    args::Vector{Arg} = Arg[],
    options::Vector{Option} = Option[],
    flags::Vector{Flag} = Flag[],
    doc = CommandDoc(),
    line::LineNumberNode = LineNumberNode(0),
)

    LeafCommand(entry, name, args, options, flags, doc, line)
end

cmd_name(cmd::EntryCommand) = cmd_name(cmd.root)
cmd_name(cmd) = cmd.name

cmd_doc(cmd::EntryCommand) = cmd_doc(cmd.root)
cmd_doc(cmd) = cmd.doc

cmd_sym(cmd) = Symbol(replace(cmd_name(cmd), "-" => "_"))
cmd_lineinfo(cmd) = cmd.line

function default_str(x::Arg)
    return "default is \e[36m" * x.default * "\e[39m. "
end
# Printings

print_option_or_flag(io::IO, xs...) = printstyled(io, xs...; color = :light_cyan)
print_args(io::IO, xs...) = printstyled(io, xs...; color = :light_magenta)

"""
    print_cmd([io, ]cmd)

Print a command object. This is used to generate command help.
"""
function print_cmd end

function Base.show(io::IO, cmd::NodeCommand)
    printstyled(io, cmd_name(cmd); color = :light_blue, bold = true)
    printstyled(io, " <command>", color = :light_blue)
end

function Base.show(io::IO, entry::EntryCommand)
    printstyled(io, cmd_name(entry); color = :light_blue, bold = true)

    if entry.version > v"0.0.0"
        print(io, " v", entry.version)
    end
end

function Base.show(io::IO, cmd::LeafCommand)
    printstyled(io, cmd_name(cmd); color = :light_blue, bold = true)

    if !isempty(cmd.options)
        printstyled(io, " [options]"; color = :light_cyan)
    end

    if !isempty(cmd.flags)
        printstyled(io, " [flags]"; color = :light_cyan)
    end

    if !isempty(cmd.args)
        print(io, " ")
        join(io, cmd.args, " ")
    end
end

function Base.show(io::IO, opt::Option)
    if opt.short
        print_option_or_flag(io, "-", first(opt.name), ", ")
    end
    print_option_or_flag(io, "--", opt.name)
    print(io, " ", opt.arg)
end

function Base.show(io::IO, opt::Flag)
    if opt.short
        print_option_or_flag(io, "-", first(opt.name), ", ")
    end
    print_option_or_flag(io, "--", opt.name)
end

function Base.show(io::IO, arg::Arg)
    notation = get(io, :notation, true)
    show_notation = notation && !arg.require
    show_notation && print_args(io, "[")

    if ignore_type(arg.type)
        print_args(io, "<", arg.name, ">")
    else
        print_args(io, "<", arg.name, "::", arg.type, ">")
    end

    show_notation && print_args(io, "]")
end

ignore_type(type) = type in [Any, String] || type <: AbstractString

print_cmd(x) = print_cmd(stdout, x)

function print_cmd(io::IO, cmd::EntryCommand)
    print_title(io, cmd)

    print_section(io, "Usage")
    println(io, " "^INDENT, cmd.root, "\n")
    print_body(io, cmd.root, true)
end

function print_cmd(io::IO, cmd::NodeCommand)
    if get(io, :inline, false)
        partition(io, cmd, cmd_doc(cmd).first)
    else
        print_title(io, cmd)
        print_body(io, cmd)
    end
end

function print_cmd(io::IO, cmd::LeafCommand)
    if get(io, :inline, false)
        partition(io, cmd, cmd_doc(cmd).first)
    else
        print_title(io, cmd)
        print_body(io, cmd)
    end
end

function print_cmd(io::IO, x::Flag)
    partition(io, x, cmd_doc(x)...)
end

function print_cmd(io::IO, x::Option)
    if x.arg.default === nothing
        partition(io, x, cmd_doc(x)...)
    else
        doc = cmd_doc(x)
        partition(io, x, doc..., _hint(doc, default_str(x.arg)))
    end
end

function print_cmd(io::IO, x::Arg)
    contents = String[]

    if !x.require
        push!(contents, "optional. ")
    end

    doc = cmd_doc(x)
    push!(contents, doc...)

    if x.default !== nothing
        push!(contents, _hint(doc, default_str(x)))
    end

    return partition(io, x, contents)
end

function _hint(doc::CommandDoc, hint::String)
    if endswith(doc, ".")
        return " " * hint
    elseif endswith(doc, ". ")
        return hint
    else
        return ". " * hint
    end
end

function print_help(io::IO)
    doc_indent = first_line_doc_indent(io, HELP_FLAG)
    print(io, " "^get(io, :indent, 0))
    print_option_or_flag(io, HELP_FLAG)
    print(io, " "^doc_indent, HELP_FLAG_DOC, "\n\n")
end

function print_version(io::IO)
    doc_indent = first_line_doc_indent(io, VERSION_FLAG)
    print(io, " "^get(io, :indent, 0))
    print_option_or_flag(io, VERSION_FLAG)
    print(io, " "^doc_indent, VERSION_FLAG_DOC, "\n\n")
end

function partition(io::IO, cmd, xs...; kwargs...)
    return partition(io, cmd, join(xs); kwargs...)
end

function partition(io::IO, cmd, xs::Vector{String}; kwargs...)
    return partition(io, cmd, join(xs); kwargs...)
end

function partition(io::IO, cmd, doc::String; width = get(io, :terminal_width, 80))
    doc_indent = get(io, :doc_indent, -1)
    doc_width = width - doc_indent
    first_line_indent = first_line_doc_indent(io, cmd)
    indent = " "^get(io, :indent, 0)
    print(io, indent, cmd)
    isempty(doc) && return

    print_doc(io, cmd, doc, doc_width, first_line_indent, indent, doc_indent)
end

# only ommit the description if they can be displayed in full later
function print_doc(
    io::IO,
    cmd::AbstractCommand,
    doc::String,
    width::Int,
    first_line_indent::Int,
    indent::String,
    doc_indent::Int,
)
    brief = first_sentence(doc)
    lines = splitlines(brief, width)
    print(io, " "^first_line_indent, first(lines))

    for i in 2:min(3, length(lines))
        print(io, "\n", indent, " "^doc_indent, lines[i])
    end

    if length(lines) > 3
        print(io, "...")
    end
end

function print_doc(
    io::IO,
    cmd,
    doc::String,
    width::Int,
    first_line_indent::Int,
    indent::String,
    doc_indent::Int,
)
    lines = splitlines(doc, width)
    print(io, " "^first_line_indent, first(lines))
    for i in 2:length(lines)
        print(io, "\n", indent, " "^doc_indent, lines[i])
    end
end

function first_sentence(content)
    index = findfirst(". ", content)
    if index === nothing
        index = findfirst(".", content)
    end

    if index === nothing
        return content
    else
        return content[1:first(index)]
    end
end

"""
    splittext(s)

Split the text in string `s` into an array, but keep all the separators
attached to the preceding word.

!!! note

    this is copied from Luxor/text.jl
"""
function splittext(s::String)
    # split text into array, keeping all separators
    # hyphens stay with first word
    result = Array{String,1}()
    iobuffer = IOBuffer()
    for c in s
        if isspace(c)
            push!(result, String(take!(iobuffer)))
            iobuffer = IOBuffer()
        elseif c == '-' # hyphen splits words but needs keeping
            print(iobuffer, c)
            push!(result, String(take!(iobuffer)))
            iobuffer = IOBuffer()
        else
            print(iobuffer, c)
        end
    end
    push!(result, String(take!(iobuffer)))
    return result
end

function splitlines(s, width = 80)
    words = splittext(s)
    lines = String[]
    current_line = String[]
    space_left = width
    for word in words
        word == "" && continue
        word_width = length(word)

        if space_left < word_width
            # start a new line
            push!(lines, strip(join(current_line)))
            current_line = String[]
            space_left = width
        end

        if endswith(word, "-")
            push!(current_line, word)
            space_left -= word_width
        else
            push!(current_line, word * " ")
            space_left -= word_width + 1
        end
    end
    isempty(current_line) || push!(lines, strip(join(current_line)))
    return lines
end

function scan_indent(x::NodeCommand, min_indent = 4)
    indent = min_indent
    indent = max(indent, length(HELP_FLAG) + min_indent)
    indent = max(indent, length(VERSION_FLAG) + min_indent)

    for cmd in x.subcmds
        indent = max(indent, length(string(cmd)) + min_indent)
    end
    return indent
end

function scan_indent(x::LeafCommand, min_indent = 4)
    indent = min_indent
    indent = max(indent, length(HELP_FLAG) + min_indent)
    indent = max(indent, length(VERSION_FLAG) + min_indent)

    for arg in x.args
        arg_str = sprint(print, arg; context = :notation => false)
        indent = max(indent, length(arg_str) + min_indent)
    end

    for opt in x.options
        indent = max(indent, length(string(opt)) + min_indent)
    end

    for flag in x.flags
        indent = max(indent, length(string(flag)) + min_indent)
    end
    return indent
end

function first_line_doc_indent(io::IO, x)
    doc_indent = get(io, :doc_indent, -1)
    notation = get(io, :notation, true)
    x_str = sprint(print, x; context = :notation => notation)
    return doc_indent > 0 ? doc_indent - length(x_str) : 4
end

function wrap_io(io::IO, x, notation = true)
    return IOContext(
        io,
        :indent => INDENT,
        :doc_indent => scan_indent(x),
        :inline => true,
        :notation => notation,
    )
end

function print_title(io::IO, x)
    println(io)
    println(io, "  ", x)
    println(io)
    println(io, cmd_doc(x))
    println(io)
end

function print_body(io::IO, cmd::NodeCommand, isentry = false)
    io = wrap_io(io, cmd)
    print_list(io, "Commands", cmd.subcmds)

    print_section(io, "Flags")
    print_help(io)
    isentry && print_version(io)
    println(io)
end

function print_body(io::IO, cmd::LeafCommand, isentry = false)
    io = wrap_io(io, cmd, false)
    if !isempty(cmd.args)
        print_list(io, "Args", cmd.args)
    end

    if !isempty(cmd.options)
        print_list(io, "Options", cmd.options)
    end

    if !isempty(cmd.flags)
        print_list(io, "Flags", cmd.flags)
        print_help(io)
        isentry && print_version(io)
    else
        print_section(io, "Flags")
        print_help(io)
        isentry && print_version(io)
    end

    println(io)
end

function print_section(io::IO, sec)
    printstyled(io, sec; bold = true)
    print(io, "\n\n")
end

function print_list(io::IO, title, list)
    print_section(io, title)

    for x in sort(list; by = x -> first(cmd_name(x)))
        print_cmd(io, x)
        print(io, "\n\n")
    end
end

# validate
function check_duplicate_short_options(options, flags)
    flags_and_options = Iterators.flatten((options, flags))

    for cmd in flags_and_options
        if cmd_name(cmd) in PRESERVED
            error("$cmd is preserved")
        end

        n_duplicate = count(flags_and_options) do x
            cmd_name(x) == cmd_name(cmd)
        end

        if n_duplicate > 1
            error("$cmd is duplicated, found $n_duplicate")
        end

        if cmd.short
            first_letter = string(first(cmd_name(cmd)))
            if first_letter in PRESERVED
                error("$cmd cannot use short version since -$first_letter is preserved.")
            end

            n_duplicate = count(flags_and_options) do x
                x.short && string(first(cmd_name(x))) == first_letter
            end

            if n_duplicate > 1
                error("the short version of $cmd is duplicated, $n_duplicate found")
            end
        end
    end
    return
end

function check_required_args(args)
    count = 0
    prev_require = 0
    for (i, arg) in enumerate(args)
        if arg.require && !(arg.vararg)
            prev_require + 1 == i || error("optional positional arguments must occur at the end")
            count += 1
            prev_require = i
        end
    end

    return count
end

function check_varargs_position(args)
    first_vararg_position = 0
    first_vararg = nothing
    for (i, arg) in enumerate(args)
        if arg.vararg
            first_vararg_position = i
            break
        end
    end
    if first_vararg_position > 0 && first_vararg_position != length(args)
        throw(Meta.ParseError("syntax: invalid \"...\" on non-final argument around $(first_vararg.line)"))
    end
    return
end

function check_first_sentence_length(name, lineinfo, doc)
    brief = first_sentence(doc)
    if length(brief) > max_brief_length()
        error(
            "the first sentence of doc should not be larger than $(max_brief_length()). ",
            "please revise the doc string for [$name] at $(lineinfo.file):$(lineinfo.line)",
        )
    end
    return brief
end

end # module


text:
@testset "split_mesh.jl" begin

end