text:
""" The discrete exterior calculus (DEC) for simplicial sets.

This module provides the dual complex associated with a delta set (the primal
complex), which is a discrete incarnation of Hodge duality, as well as the many
operators of the DEC that depend on it, such as the Hodge star, codifferential,
wedge product, interior product, and Lie derivative. The main reference for this
module is Hirani's 2003 PhD thesis.
"""
module DiscreteExteriorCalculus
export DualSimplex, DualV, DualE, DualTri, DualChain, DualForm,
  PrimalVectorField, DualVectorField,
  AbstractDeltaDualComplex1D, DeltaDualComplex1D,
  OrientedDeltaDualComplex1D, EmbeddedDeltaDualComplex1D,
  AbstractDeltaDualComplex2D, DeltaDualComplex2D,
  OrientedDeltaDualComplex2D, EmbeddedDeltaDualComplex2D,
  SimplexCenter, Barycenter, Circumcenter, Incenter, geometric_center,
  subsimplices, primal_vertex, elementary_duals, dual_boundary, dual_derivative,
  ⋆, hodge_star, inv_hodge_star, δ, codifferential, ∇², laplace_beltrami, Δ, laplace_de_rham,
  ♭, flat, ♯, sharp, ∧, wedge_product, interior_product, interior_product_flat,
  ℒ, lie_derivative, lie_derivative_flat,
  vertex_center, edge_center, triangle_center, dual_triangle_vertices,
  dual_point, dual_volume, subdivide_duals!, DiagonalHodge, GeometricHodge

import Base: ndims
using LinearAlgebra: Diagonal, dot, norm, cross
using SparseArrays
using StaticArrays: @SVector, SVector

using Catlab, Catlab.CategoricalAlgebra.CSets
using Catlab.CategoricalAlgebra.FinSets: deleteat
using Catlab.Graphs: HasGraph
import Catlab.Theories: Δ
using ..ArrayUtils, ..SimplicialSets
using ..SimplicialSets: DeltaCategory1D, DeltaCategory2D, CayleyMengerDet,
  operator_nz, ∂_nz, d_nz, cayley_menger, negate
import ..SimplicialSets: ∂, d, volume

abstract type DiscreteFlat end
struct DPPFlat <: DiscreteFlat end

abstract type DiscreteSharp end
struct PPSharp <: DiscreteSharp end

abstract type DiscreteHodge end
struct GeometricHodge <: DiscreteHodge end
struct DiagonalHodge  <: DiscreteHodge end

# 1D dual complex
#################

# Should be expressed using a coproduct of two copies of `DeltaCategory1D`.

@present SchemaDualComplex1D <: DeltaCategory1D begin
  # Dual vertices and edges.
  (DualV, DualE)::Ob
  (D_∂v0, D_∂v1)::Hom(DualE, DualV)

  # Centers of primal simplices are dual vertices.
  vertex_center::Hom(V, DualV)
  edge_center::Hom(E, DualV)

  # Every primal edge is subdivided into two dual edges.
  #
  # (∂v0_dual, ∂v1_dual)::Hom(E,DualE)
  #
  # ∂v0_dual ⋅ D_∂v1 == ∂v0 ⋅ vertex_center
  # ∂v1_dual ⋅ D_∂v1 == ∂v1 ⋅ vertex_center
  # ∂v0_dual ⋅ D_∂v0 == edge_center
  # ∂v1_dual ⋅ D_∂v0 == edge_center
  #
  # We could, and arguably should, track these through dedicated morphisms, as
  # in the commented code above. We don't because it scales badly in dimension:
  # an ``n``-simplex is subdivided into ``n!`` sub-simplices. So we would need 6
  # morphisms for triangles and 24 for tets, and an even larger number of
  # equations to say how everything fits together. Moreover, in practice, the
  # incidence data for the dual complex suffices to construct the dual cells of
  # primal simplices.
end

""" Abstract type for dual complex of a 1D delta set.
"""
@abstract_acset_type AbstractDeltaDualComplex1D <: HasDeltaSet1D

""" Dual complex of a one-dimensional delta set.

The data structure includes both the primal complex and the dual complex, as
well as the mapping between them.
"""
@acset_type DeltaDualComplex1D(SchemaDualComplex1D,
  index=[:∂v0,:∂v1,:D_∂v0,:D_∂v1]) <: AbstractDeltaDualComplex1D

""" Dual vertex corresponding to center of primal vertex.
"""
vertex_center(s::HasDeltaSet, args...) = s[args..., :vertex_center]

""" Dual vertex corresponding to center of primal edge.
"""
edge_center(s::HasDeltaSet1D, args...) = s[args..., :edge_center]

subsimplices(::Type{Val{1}}, s::HasDeltaSet1D, e::Int) =
  SVector{2}(incident(s, edge_center(s, e), :D_∂v0))

primal_vertex(::Type{Val{1}}, s::HasDeltaSet1D, e...) = s[e..., :D_∂v1]

elementary_duals(::Type{Val{0}}, s::AbstractDeltaDualComplex1D, v::Int) =
  incident(s, vertex_center(s,v), :D_∂v1)
elementary_duals(::Type{Val{1}}, s::AbstractDeltaDualComplex1D, e::Int) =
  SVector(edge_center(s,e))

""" Boundary dual vertices of a dual triangle.

This accessor assumes that the simplicial identities for the dual hold.
"""
function dual_triangle_vertices(s::HasDeltaSet1D, t...)
  SVector(s[s[t..., :D_∂e1], :D_∂v1],
          s[s[t..., :D_∂e0], :D_∂v1],
          s[s[t..., :D_∂e0], :D_∂v0])
end

# 1D oriented dual complex
#-------------------------

@present SchemaOrientedDualComplex1D <: SchemaDualComplex1D begin
  Orientation::AttrType
  edge_orientation::Attr(E, Orientation)
  D_edge_orientation::Attr(DualE, Orientation)
end

""" Oriented dual complex of an oriented 1D delta set.
"""
@acset_type OrientedDeltaDualComplex1D(SchemaOrientedDualComplex1D,
  index=[:∂v0,:∂v1,:D_∂v0,:D_∂v1]) <: AbstractDeltaDualComplex1D

dual_boundary_nz(::Type{Val{1}}, s::AbstractDeltaDualComplex1D, x::Int) =
  # Boundary vertices of dual 1-cell ↔
  # Dual vertices for cofaces of (edges incident to) primal vertex.
  d_nz(Val{0}, s, x)

dual_derivative_nz(::Type{Val{0}}, s::AbstractDeltaDualComplex1D, x::Int) =
  negatenz(∂_nz(Val{1}, s, x))

negatenz((I, V)) = (I, negate.(V))

""" Construct 1D dual complex from 1D delta set.
"""
function (::Type{S})(t::AbstractDeltaSet1D) where S <: AbstractDeltaDualComplex1D
  s = S()
  copy_parts!(s, t)
  make_dual_simplices_1d!(s)
  return s
end

""" Make dual vertice and edges for dual complex of dimension ≧ 1.

Although zero-dimensional duality is geometrically trivial (subdividing a vertex
gives back the same vertex), we treat the dual vertices as disjoint from the
primal vertices. Thus, a dual vertex is created for every primal vertex.

If the primal complex is oriented, an orientation is induced on the dual
complex. The dual edges are oriented relative to the primal edges they subdivide
(Hirani 2003, PhD thesis, Ch. 2, last sentence of Remark 2.5.1).
"""
function make_dual_simplices_1d!(s::HasDeltaSet1D)
  # Make dual vertices and edges.
  s[:vertex_center] = vcenters = add_parts!(s, :DualV, nv(s))
  s[:edge_center] = ecenters = add_parts!(s, :DualV, ne(s))
  D_edges = map((0,1)) do i
    add_parts!(s, :DualE, ne(s);
               D_∂v0 = ecenters, D_∂v1 = view(vcenters, ∂(1,i,s)))
  end

  # Orient elementary dual edges.
  if has_subpart(s, :edge_orientation)
    edge_orient = s[:edge_orientation]
    s[D_edges[1], :D_edge_orientation] = negate.(edge_orient)
    s[D_edges[2], :D_edge_orientation] = edge_orient
  end

  D_edges
end

# 1D embedded dual complex
#-------------------------

@present SchemaEmbeddedDualComplex1D <: SchemaOrientedDualComplex1D begin
  (Real, Point)::AttrType
  point::Attr(V, Point)
  length::Attr(E, Real)
  dual_point::Attr(DualV, Point)
  dual_length::Attr(DualE, Real)
end

""" Embedded dual complex of an embedded 1D delta set.

Although they are redundant information, the lengths of the primal and dual
edges are precomputed and stored.
"""
@acset_type EmbeddedDeltaDualComplex1D(SchemaEmbeddedDualComplex1D,
  index=[:∂v0,:∂v1,:D_∂v0,:D_∂v1]) <: AbstractDeltaDualComplex1D

""" Point associated with dual vertex of complex.
"""
dual_point(s::HasDeltaSet, args...) = s[args..., :dual_point]

struct PrecomputedVol end

volume(::Type{Val{n}}, s::EmbeddedDeltaDualComplex1D, x) where n =
  volume(Val{n}, s, x, PrecomputedVol())
dual_volume(::Type{Val{n}}, s::EmbeddedDeltaDualComplex1D, x) where n =
  dual_volume(Val{n}, s, x, PrecomputedVol())

volume(::Type{Val{1}}, s::HasDeltaSet1D, e, ::PrecomputedVol) = s[e, :length]
dual_volume(::Type{Val{1}}, s::HasDeltaSet1D, e, ::PrecomputedVol) =
  s[e, :dual_length]

dual_volume(::Type{Val{1}}, s::HasDeltaSet1D, e::Int, ::CayleyMengerDet) =
  volume(dual_point(s, SVector(s[e,:D_∂v0], s[e,:D_∂v1])))

hodge_diag(::Type{Val{0}}, s::AbstractDeltaDualComplex1D, v::Int) =
  sum(dual_volume(Val{1}, s, elementary_duals(Val{0},s,v)))
hodge_diag(::Type{Val{1}}, s::AbstractDeltaDualComplex1D, e::Int) =
  1 / volume(Val{1},s,e)

""" Compute geometric subdivision for embedded dual complex.

Supports different methods of subdivision through the choice of geometric
center, as defined by [`geometric_center`](@ref). In particular, barycentric
subdivision and circumcentric subdivision are supported.
"""
function subdivide_duals!(s::EmbeddedDeltaDualComplex1D, args...)
  subdivide_duals_1d!(s, args...)
  precompute_volumes_1d!(s)
end

function subdivide_duals_1d!(s::HasDeltaSet1D, alg)
  for v in vertices(s)
    s[vertex_center(s,v), :dual_point] = point(s, v)
  end
  for e in edges(s)
    s[edge_center(s,e), :dual_point] = geometric_center(
      point(s, edge_vertices(s, e)), alg)
  end
end

function precompute_volumes_1d!(s::HasDeltaSet1D)
  for e in edges(s)
    s[e, :length] = volume(1,s,e,CayleyMengerDet())
  end
  for e in parts(s, :DualE)
    s[e, :dual_length] = dual_volume(1,s,e,CayleyMengerDet())
  end
end

# 2D dual complex
#################

# Should be expressed using a coproduct of two copies of `DeltaCategory2D` or
# perhaps a pushout of `SchemaDualComplex2D` and `DeltaCategory1D`.

@present SchemaDualComplex2D <: DeltaCategory2D begin
  # Dual vertices, edges, and triangles.
  (DualV, DualE, DualTri)::Ob
  (D_∂v0, D_∂v1)::Hom(DualE, DualV)
  (D_∂e0, D_∂e1, D_∂e2)::Hom(DualTri, DualE)

  # Simplicial identities for dual simplices.
  D_∂e1 ⋅ D_∂v1 == D_∂e2 ⋅ D_∂v1
  D_∂e0 ⋅ D_∂v1 == D_∂e2 ⋅ D_∂v0
  D_∂e0 ⋅ D_∂v0 == D_∂e1 ⋅ D_∂v0

  # Centers of primal simplices are dual vertices.
  vertex_center::Hom(V, DualV)
  edge_center::Hom(E, DualV)
  tri_center::Hom(Tri, DualV)
end

""" Abstract type for dual complex of a 2D delta set.
"""
@abstract_acset_type AbstractDeltaDualComplex2D <: HasDeltaSet2D

""" Dual complex of a two-dimensional delta set.
"""
@acset_type DeltaDualComplex2D(SchemaDualComplex2D,
  index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:D_∂v0,:D_∂v1,:D_∂e0,:D_∂e1,:D_∂e2]) <: AbstractDeltaDualComplex2D

""" Dual vertex corresponding to center of primal triangle.
"""
triangle_center(s::HasDeltaSet2D, args...) = s[args..., :tri_center]

subsimplices(::Type{Val{2}}, s::HasDeltaSet2D, t::Int) =
  SVector{6}(incident(s, triangle_center(s,t), @SVector [:D_∂e1, :D_∂v0]))

primal_vertex(::Type{Val{2}}, s::HasDeltaSet2D, t...) =
  primal_vertex(Val{1}, s, s[t..., :D_∂e2])

elementary_duals(::Type{Val{0}}, s::AbstractDeltaDualComplex2D, v::Int) =
  incident(s, vertex_center(s,v), @SVector [:D_∂e1, :D_∂v1])
elementary_duals(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, e::Int) =
  incident(s, edge_center(s,e), :D_∂v1)
elementary_duals(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, t::Int) =
  SVector(triangle_center(s,t))

# 2D oriented dual complex
#-------------------------

@present SchemaOrientedDualComplex2D <: SchemaDualComplex2D begin
  Orientation::AttrType
  edge_orientation::Attr(E, Orientation)
  tri_orientation::Attr(Tri, Orientation)
  D_edge_orientation::Attr(DualE, Orientation)
  D_tri_orientation::Attr(DualTri, Orientation)
end

""" Oriented dual complex of an oriented 2D delta set.
"""
@acset_type OrientedDeltaDualComplex2D(SchemaOrientedDualComplex2D,
  index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:D_∂v0,:D_∂v1,:D_∂e0,:D_∂e1,:D_∂e2]) <: AbstractDeltaDualComplex2D

dual_boundary_nz(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, x::Int) =
  # Boundary vertices of dual 1-cell ↔
  # Dual vertices for cofaces of (triangles incident to) primal edge.
  negatenz(d_nz(Val{1}, s, x))
dual_boundary_nz(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, x::Int) =
  # Boundary edges of dual 2-cell ↔
  # Dual edges for cofaces of (edges incident to) primal vertex.
  d_nz(Val{0}, s, x)

dual_derivative_nz(::Type{Val{0}}, s::AbstractDeltaDualComplex2D, x::Int) =
  ∂_nz(Val{2}, s, x)
dual_derivative_nz(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, x::Int) =
  negatenz(∂_nz(Val{1}, s, x))

""" Construct 2D dual complex from 2D delta set.
"""
function (::Type{S})(t::AbstractDeltaSet2D) where S <: AbstractDeltaDualComplex2D
  s = S()
  copy_parts!(s, t)
  make_dual_simplices_2d!(s)
  return s
end

""" Make dual simplices for dual complex of dimension ≧ 2.

If the primal complex is oriented, an orientation is induced on the dual
complex. The elementary dual edges are oriented following (Hirani, 2003, Example
2.5.2) or (Desbrun et al, 2005, Table 1) and the dual triangles are oriented
relative to the primal triangles they subdivide.
"""
function make_dual_simplices_2d!(s::HasDeltaSet2D)
  # Make dual vertices and edges.
  D_edges01 = make_dual_simplices_1d!(s)
  s[:tri_center] = tri_centers = add_parts!(s, :DualV, ntriangles(s))
  D_edges12 = map((0,1,2)) do e
    add_parts!(s, :DualE, ntriangles(s);
               D_∂v0=tri_centers, D_∂v1=edge_center(s, ∂(2,e,s)))
  end
  D_edges02 = map(triangle_vertices(s)) do vs
    add_parts!(s, :DualE, ntriangles(s);
               D_∂v0=tri_centers, D_∂v1=vertex_center(s, vs))
  end

  # Make dual triangles.
  # Counterclockwise order in drawing with vertices 0, 1, 2 from left to right.
  D_triangle_schemas = ((0,1,1),(0,2,1),(1,2,0),(1,0,1),(2,0,0),(2,1,0))
  D_triangles = map(D_triangle_schemas) do (v,e,ev)
    add_parts!(s, :DualTri, ntriangles(s);
               D_∂e0=D_edges12[e+1], D_∂e1=D_edges02[v+1],
               D_∂e2=view(D_edges01[ev+1], ∂(2,e,s)))
  end

  if has_subpart(s, :tri_orientation)
    # Orient elementary dual triangles.
    tri_orient = s[:tri_orientation]
    rev_tri_orient = negate.(tri_orient)
    for (i, D_tris) in enumerate(D_triangles)
      s[D_tris, :D_tri_orientation] = isodd(i) ? rev_tri_orient : tri_orient
    end

    # Orient elementary dual edges.
    for e in (0,1,2)
      s[D_edges12[e+1], :D_edge_orientation] = relative_sign.(
        s[∂(2,e,s), :edge_orientation],
        isodd(e) ? rev_tri_orient : tri_orient)
    end
    # Remaining dual edges are oriented arbitrarily.
    s[lazy(vcat, D_edges02...), :D_edge_orientation] = one(eltype(tri_orient))
  end

  D_triangles
end

relative_sign(x, y) = sign(x*y)
relative_sign(x::Bool, y::Bool) = (x && y) || (!x && !y)

# 2D embedded dual complex
#-------------------------

@present SchemaEmbeddedDualComplex2D <: SchemaOrientedDualComplex2D begin
  (Real, Point)::AttrType
  point::Attr(V, Point)
  length::Attr(E, Real)
  area::Attr(Tri, Real)
  dual_point::Attr(DualV, Point)
  dual_length::Attr(DualE, Real)
  dual_area::Attr(DualTri, Real)
end

""" Embedded dual complex of an embedded 12 delta set.

Although they are redundant information, the lengths and areas of the
primal/dual edges and triangles are precomputed and stored.
"""
@acset_type EmbeddedDeltaDualComplex2D(SchemaEmbeddedDualComplex2D,
  index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2,:D_∂v0,:D_∂v1,:D_∂e0,:D_∂e1,:D_∂e2]) <: AbstractDeltaDualComplex2D

volume(::Type{Val{n}}, s::EmbeddedDeltaDualComplex2D, x) where n =
  volume(Val{n}, s, x, PrecomputedVol())
dual_volume(::Type{Val{n}}, s::EmbeddedDeltaDualComplex2D, x) where n =
  dual_volume(Val{n}, s, x, PrecomputedVol())

volume(::Type{Val{2}}, s::HasDeltaSet2D, t, ::PrecomputedVol) = s[t, :area]
dual_volume(::Type{Val{2}}, s::HasDeltaSet2D, t, ::PrecomputedVol) =
  s[t, :dual_area]

function dual_volume(::Type{Val{2}}, s::HasDeltaSet2D, t::Int, ::CayleyMengerDet)
  dual_vs = SVector(s[s[t, :D_∂e1], :D_∂v1],
                    s[s[t, :D_∂e2], :D_∂v0],
                    s[s[t, :D_∂e0], :D_∂v0])
  volume(dual_point(s, dual_vs))
end

hodge_diag(::Type{Val{0}}, s::AbstractDeltaDualComplex2D, v::Int) =
  sum(dual_volume(Val{2}, s, elementary_duals(Val{0},s,v)))
hodge_diag(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, e::Int) =
  sum(dual_volume(Val{1}, s, elementary_duals(Val{1},s,e))) / volume(Val{1},s,e)
hodge_diag(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, t::Int) =
  1 / volume(Val{2},s,t)

function ♭(s::AbstractDeltaDualComplex2D, X::AbstractVector, ::DPPFlat)
  # XXX: Creating this lookup table shouldn't be necessary. Of course, we could
  # index `tri_center` but that shouldn't be necessary either. Rather, we should
  # loop over incident triangles instead of the elementary duals, which just
  # happens to be inconvenient.
  tri_map = Dict{Int,Int}(triangle_center(s,t) => t for t in triangles(s))

  map(edges(s)) do e
    e_vec = (point(s, tgt(s,e)) - point(s, src(s,e))) * sign(1,s,e)
    dual_edges = elementary_duals(1,s,e)
    dual_lengths = dual_volume(1, s, dual_edges)
    mapreduce(+, dual_edges, dual_lengths) do dual_e, dual_length
      X_vec = X[tri_map[s[dual_e, :D_∂v0]]]
      dual_length * dot(X_vec, e_vec)
    end / sum(dual_lengths)
  end
end

function ♯(s::AbstractDeltaDualComplex2D, α::AbstractVector, ::PPSharp)
  α♯ = zeros(eltype(s[:dual_point]), nv(s))
  for t in triangles(s)
    area = volume(2,s,t)
    tri_center, tri_edges = triangle_center(s,t), triangle_edges(s,t)
    tri_point = dual_point(s, tri_center)
    for (i, (v₀, e₀)) in enumerate(zip(triangle_vertices(s,t), tri_edges))
      e_vec = point(s, tgt(s, e₀)) - point(s, src(s, e₀))
      e_vec /= norm(e_vec)
      e2_vec = point(s, v₀) - point(s, src(s, e₀))
      out_vec = e2_vec - dot(e2_vec, e_vec)*e_vec
      h = norm(out_vec)
      out_vec /= h^2 # length == 1/h
      for e in deleteat(tri_edges, i)
        v, sgn = src(s,e) == v₀ ? (tgt(s,e), -1) : (src(s,e), +1)
        dual_area = sum(dual_volume(2,s,d) for d in elementary_duals(0,s,v)
                        if s[s[d, :D_∂e0], :D_∂v0] == tri_center)
        α♯[v] += sgn * sign(1,s,e) * α[e] * (dual_area / area) * out_vec
      end
    end
  end
  α♯
end

function ∧(::Type{Tuple{1,1}}, s::HasDeltaSet2D, α, β, x::Int)
  # XXX: This calculation of the volume coefficients is awkward due to the
  # design decision described in `SchemaDualComplex1D`.
  dual_vs = vertex_center(s, triangle_vertices(s, x))
  dual_es = sort(SVector{6}(incident(s, triangle_center(s, x), :D_∂v0)),
                 by=e -> s[e,:D_∂v1] .== dual_vs, rev=true)[1:3]
  coeffs = map(dual_es) do e
    sum(dual_volume(2, s, SVector{2}(incident(s, e, :D_∂e1))))
  end / volume(2, s, x)

  # Wedge product of two primal 1-forms, as in (Hirani 2003, Example 7.1.2).
  # This formula is not the same as (Hirani 2003, Equation 7.1.2) but it is
  # equivalent.
  e0, e1, e2 = ∂(2,0,s,x), ∂(2,1,s,x), ∂(2,2,s,x)
  dot(coeffs, SVector(α[e2] * β[e1] - α[e1] * β[e2],
                      α[e2] * β[e0] - α[e0] * β[e2],
                      α[e1] * β[e0] - α[e0] * β[e1])) / 2
end

function subdivide_duals!(s::EmbeddedDeltaDualComplex2D, args...)
  subdivide_duals_2d!(s, args...)
  precompute_volumes_2d!(s)
end

function subdivide_duals_2d!(s::HasDeltaSet2D, alg)
  subdivide_duals_1d!(s, alg)
  for t in triangles(s)
    s[triangle_center(s,t), :dual_point] = geometric_center(
      point(s, triangle_vertices(s, t)), alg)
  end
end

function precompute_volumes_2d!(s::HasDeltaSet2D)
  precompute_volumes_1d!(s)
  for t in triangles(s)
    s[t, :area] = volume(2,s,t,CayleyMengerDet())
  end
  for t in parts(s, :DualTri)
    s[t, :dual_area] = dual_volume(2,s,t,CayleyMengerDet())
  end
end

# General operators
###################

""" Wrapper for dual simplex or simplices of dimension `D`.

See also: [`DualV`](@ref), [`DualE`](@ref), [`DualTri`](@ref).
"""
@parts_array_struct DualSimplex{D}

""" Vertex in simplicial set: alias for `Simplex{0}`.
"""
const DualV = DualSimplex{0}

""" Edge in simplicial set: alias for `Simplex{1}`.
"""
const DualE = DualSimplex{1}

""" Triangle in simplicial set: alias for `Simplex{2}`.
"""
const DualTri = DualSimplex{2}

""" Wrapper for chain of dual cells of dimension `n`.

In an ``N``-dimensional complex, the elementary dual simplices of each
``n``-simplex together comprise the dual ``(N-n)``-cell of the simplex. Using
this correspondence, a basis for primal ``n``-chains defines the basis for dual
``(N-n)``-chains.

!!! note

    In (Hirani 2003, Definition 3.4.1), the duality operator assigns a certain
    sign to each elementary dual simplex. For us, all of these signs should be
    regarded as positive because we have already incorporated them into the
    orientation of the dual simplices.
"""
@vector_struct DualChain{n}

""" Wrapper for form, aka cochain, on dual cells of dimension `n`.
"""
@vector_struct DualForm{n}

""" Wrapper for vector field on primal vertices.
"""
@vector_struct PrimalVectorField

""" Wrapper for vector field on dual vertices.
"""
@vector_struct DualVectorField

ndims(s::AbstractDeltaDualComplex1D) = 1
ndims(s::AbstractDeltaDualComplex2D) = 2

volume(s::HasDeltaSet, x::DualSimplex{n}, args...) where n =
  dual_volume(Val{n}, s, x.data, args...)
@inline dual_volume(n::Int, s::HasDeltaSet, args...) =
  dual_volume(Val{n}, s, args...)

""" List of dual simplices comprising the subdivision of a primal simplex.

A primal ``n``-simplex is always subdivided into ``n!`` dual ``n``-simplices,
not be confused with the [`elementary_duals`](@ref) which have complementary
dimension.

The returned list is ordered such that subsimplices with the same primal vertex
appear consecutively.
"""
subsimplices(s::HasDeltaSet, x::Simplex{n}) where n =
  DualSimplex{n}(subsimplices(Val{n}, s, x.data))
@inline subsimplices(n::Int, s::HasDeltaSet, args...) =
  subsimplices(Val{n}, s, args...)

""" Primal vertex associated with a dual simplex.
"""
primal_vertex(s::HasDeltaSet, x::DualSimplex{n}) where n =
  V(primal_vertex(Val{n}, s, x.data))
@inline primal_vertex(n::Int, s::HasDeltaSet, args...) =
  primal_vertex(Val{n}, s, args...)

""" List of elementary dual simplices corresponding to primal simplex.

In general, in an ``n``-dimensional complex, the elementary duals of primal
``k``-simplices are dual ``(n-k)``-simplices. Thus, in 1D dual complexes, the
elementary duals of...

- primal vertices are dual edges
- primal edges are (single) dual vertices

In 2D dual complexes, the elementary duals of...

- primal vertices are dual triangles
- primal edges are dual edges
- primal triangles are (single) dual triangles
"""
elementary_duals(s::HasDeltaSet, x::Simplex{n}) where n =
  DualSimplex{ndims(s)-n}(elementary_duals(Val{n}, s, x.data))
@inline elementary_duals(n::Int, s::HasDeltaSet, args...) =
  elementary_duals(Val{n}, s, args...)

""" Boundary of chain of dual cells.

Transpose of [`dual_derivative`](@ref).
"""
@inline dual_boundary(n::Int, s::HasDeltaSet, args...) =
  dual_boundary(Val{n}, s, args...)
∂(s::HasDeltaSet, x::DualChain{n}) where n =
  DualChain{n-1}(dual_boundary(Val{n}, s, x.data))

function dual_boundary(::Type{Val{n}}, s::HasDeltaSet, args...) where n
  operator_nz(Int, nsimplices(ndims(s)-n+1,s),
              nsimplices(ndims(s)-n,s), args...) do x
    dual_boundary_nz(Val{n}, s, x)
  end
end

""" Discrete exterior derivative of dual form.

Transpose of [`dual_boundary`](@ref). For more info, see (Desbrun, Kanso, Tong,
2008: Discrete differential forms for computational modeling, §4.5).
"""
@inline dual_derivative(n::Int, s::HasDeltaSet, args...) =
  dual_derivative(Val{n}, s, args...)
d(s::HasDeltaSet, x::DualForm{n}) where n =
  DualForm{n+1}(dual_derivative(Val{n}, s, x.data))

function dual_derivative(::Type{Val{n}}, s::HasDeltaSet, args...) where n
  operator_nz(Int, nsimplices(ndims(s)-n-1,s),
              nsimplices(ndims(s)-n,s), args...) do x
    dual_derivative_nz(Val{n}, s, x)
  end
end

""" Hodge star operator from primal ``n``-forms to dual ``N-n``-forms.

!!! note

    Some authors, such as (Hirani 2003) and (Desbrun 2005), use the symbol ``⋆``
    for the duality operator on chains and the symbol ``*`` for the Hodge star
    operator on cochains. We do not explicitly define the duality operator and
    we use the symbol ``⋆`` for the Hodge star.
"""
⋆(s::HasDeltaSet, x::SimplexForm{n}; kw...) where n =
  DualForm{ndims(s)-n}(⋆(Val{n}, s, x.data; kw...))
@inline ⋆(n::Int, s::HasDeltaSet, args...; kw...) = ⋆(Val{n}, s, args...; kw...)
@inline ⋆(::Type{Val{n}}, s::HasDeltaSet;
          hodge::DiscreteHodge=GeometricHodge()) where n = ⋆(Val{n}, s, hodge)
@inline ⋆(::Type{Val{n}}, s::HasDeltaSet, form::AbstractVector;
          hodge::DiscreteHodge=GeometricHodge()) where n = ⋆(Val{n}, s, form, hodge)

⋆(::Type{Val{n}}, s::HasDeltaSet, form::AbstractVector, ::DiagonalHodge) where n =
  applydiag(form) do x, a; a * hodge_diag(Val{n},s,x) end
⋆(::Type{Val{n}}, s::HasDeltaSet, ::DiagonalHodge) where n =
  Diagonal([ hodge_diag(Val{n},s,x) for x in simplices(n,s) ])

# Note that this cross product defines the positive direction for flux to
# always be in the positive z direction. This will likely not generalize to
# arbitrary meshes embedded in 3D space, and so will need to be revisited.
# Potentially this orientation can be provided by the simplicial triangle
# orientation?
crossdot(v1, v2) = begin
  v1v2 = cross(v1, v2)
  norm(v1v2) * (last(v1v2) == 0 ? 1.0 : sign(last(v1v2)))
end

""" Hodge star operator from primal 1-forms to dual 1-forms.

This specific hodge star implementation is based on the hodge star presented in
(Ayoub et al 2020), which generalizes the operator presented in (Hirani 2003).
This reproduces the diagonal hodge for a dual mesh generated under
circumcentric subdivision and provides off-diagonal correction factors for
meshes generated under other subdivision schemes (e.g. barycentric).
"""
function ⋆(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, ::GeometricHodge)

  vals = Dict{Tuple{Int64, Int64}, Float64}()
  I = Vector{Int64}()
  J = Vector{Int64}()
  V = Vector{Float64}()

  rel_orient = 0.0
  for t in triangles(s)
    e = reverse(triangle_edges(s, t))
    ev = point(s, tgt(s, e)) .- point(s, src(s,e))

    tc = dual_point(s, triangle_center(s, t))
    dv = map(enumerate(dual_point(s, edge_center(s, e)))) do (i,v)
      (tc - v) * (i == 2 ? -1 : 1)
    end

    diag_dot = map(1:3) do i
      dot(ev[i], dv[i]) / dot(ev[i], ev[i])
    end

    # This relative orientation needs to be redefined for each triangle in the
    # case that the mesh has multiple independent connected components
    rel_orient = 0.0
    for i in 1:3
      diag_cross = sign(Val{2}, s, t) * crossdot(ev[i], dv[i]) /
                      dot(ev[i], ev[i])
      if diag_cross != 0.0
        # Decide the orientation of the mesh relative to z-axis (see crossdot)
        # For optimization, this could be moved out of this loop
        if rel_orient == 0.0
          rel_orient = sign(diag_cross)
        end

        push!(I, e[i])
        push!(J, e[i])
        push!(V, diag_cross * rel_orient)
      end
    end

    for p ∈ ((1,2,3), (1,3,2), (2,1,3),
             (2,3,1), (3,1,2), (3,2,1))
      val = rel_orient * sign(Val{2}, s, t) * diag_dot[p[1]] *
              dot(ev[p[1]], ev[p[3]]) / crossdot(ev[p[2]], ev[p[3]])
      if val != 0.0
        push!(I, e[p[1]])
        push!(J, e[p[2]])
        push!(V, val)
      end
    end
  end
  sparse(I,J,V)
end

⋆(::Type{Val{0}}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) =
  ⋆(Val{0}, s, DiagonalHodge())
⋆(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) =
  ⋆(Val{2}, s, DiagonalHodge())

⋆(::Type{Val{0}}, s::AbstractDeltaDualComplex2D, form::AbstractVector, ::GeometricHodge) =
  ⋆(Val{0}, s, form, DiagonalHodge())
⋆(::Type{Val{1}}, s::AbstractDeltaDualComplex2D, form::AbstractVector, ::GeometricHodge) =
  ⋆(Val{1}, s, GeometricHodge()) * form
⋆(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, form::AbstractVector, ::GeometricHodge) =
  ⋆(Val{2}, s, form, DiagonalHodge())

⋆(::Type{Val{n}}, s::AbstractDeltaDualComplex1D, ::GeometricHodge) where n =
  ⋆(Val{n}, s, DiagonalHodge())
⋆(::Type{Val{n}}, s::AbstractDeltaDualComplex1D, form::AbstractVector, ::GeometricHodge) where n =
  ⋆(Val{n}, s, form, DiagonalHodge())

""" Alias for the Hodge star operator [`⋆`](@ref).
"""
const hodge_star = ⋆

""" Inverse Hodge star operator from dual ``N-n``-forms to primal ``n``-forms.

Confusingly, this is *not* the operator inverse of the Hodge star [`⋆`](@ref)
because it carries an extra global sign, in analogy to the smooth case
(Gillette, 2009, Notes on the DEC, Definition 2.27).
"""
@inline inv_hodge_star(n::Int, s::HasDeltaSet, args...; kw...) =
  inv_hodge_star(Val{n}, s, args...; kw...)
@inline inv_hodge_star(::Type{Val{n}}, s::HasDeltaSet;
                       hodge::DiscreteHodge=GeometricHodge()) where n =
  inv_hodge_star(Val{n}, s, hodge)
@inline inv_hodge_star(::Type{Val{n}}, s::HasDeltaSet, form::AbstractVector;
                       hodge::DiscreteHodge=GeometricHodge()) where n =
  inv_hodge_star(Val{n}, s, form, hodge)

function inv_hodge_star(::Type{Val{n}}, s::HasDeltaSet,
                        form::AbstractVector, ::DiagonalHodge) where n
  if iseven(n*(ndims(s)-n))
    applydiag(form) do x, a; a / hodge_diag(Val{n},s,x) end
  else
    applydiag(form) do x, a; -a / hodge_diag(Val{n},s,x) end
  end
end

function inv_hodge_star(::Type{Val{n}}, s::HasDeltaSet, ::DiagonalHodge) where n
  if iseven(n*(ndims(s)-n))
    Diagonal([ 1 / hodge_diag(Val{n},s,x) for x in simplices(n,s) ])
  else
    Diagonal([ -1 / hodge_diag(Val{n},s,x) for x in simplices(n,s) ])
  end
end

function inv_hodge_star(::Type{Val{1}}, s::AbstractDeltaDualComplex2D,
                        ::GeometricHodge)
  -1 * inv(Matrix(⋆(Val{1}, s, GeometricHodge())))
end
function inv_hodge_star(::Type{Val{1}}, s::AbstractDeltaDualComplex2D,
                        form::AbstractVector, ::GeometricHodge)
  -1 * (Matrix(⋆(Val{1}, s, GeometricHodge())) \ form)
end

inv_hodge_star(::Type{Val{0}}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) =
  inv_hodge_star(Val{0}, s, DiagonalHodge())
inv_hodge_star(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, ::GeometricHodge) =
  inv_hodge_star(Val{2}, s, DiagonalHodge())

inv_hodge_star(::Type{Val{0}}, s::AbstractDeltaDualComplex2D,
               form::AbstractVector, ::GeometricHodge) =
  inv_hodge_star(Val{0}, s, form, DiagonalHodge())
inv_hodge_star(::Type{Val{2}}, s::AbstractDeltaDualComplex2D,
               form::AbstractVector, ::GeometricHodge) =
  inv_hodge_star(Val{2}, s, form, DiagonalHodge())

inv_hodge_star(::Type{Val{n}}, s::AbstractDeltaDualComplex1D,
               ::GeometricHodge) where n =
  inv_hodge_star(Val{n}, s, DiagonalHodge())
inv_hodge_star(::Type{Val{n}}, s::AbstractDeltaDualComplex1D,
               form::AbstractVector, ::GeometricHodge) where n =
  inv_hodge_star(Val{n}, s, form, DiagonalHodge())

""" Codifferential operator from primal ``n`` forms to primal ``n-1``-forms.
"""
δ(s::HasDeltaSet, x::SimplexForm{n}; kw...) where n =
  SimplexForm{n-1}(δ(Val{n}, s, GeometricHodge(), x.data; kw...))
@inline δ(n::Int, s::HasDeltaSet, args...; kw...) =
  δ(Val{n}, s, args...; kw...)
@inline δ(::Type{Val{n}}, s::HasDeltaSet; hodge::DiscreteHodge=GeometricHodge(),
          matrix_type::Type=SparseMatrixCSC{Float64}) where n =
  δ(Val{n}, s, hodge, matrix_type)
@inline δ(::Type{Val{n}}, s::HasDeltaSet, form::AbstractVector;
          hodge::DiscreteHodge=GeometricHodge()) where n =
  δ(Val{n}, s, hodge, form)

function δ(::Type{Val{n}}, s::HasDeltaSet, ::DiagonalHodge, args...) where n
  # The sign of δ in Gillette's notes (see test file) is simply a product of
  # the signs for the inverse hodge and dual derivative involved.
  sgn = iseven((n-1)*(ndims(s)*(n-1) + 1)) ? +1 : -1
  operator_nz(Float64, nsimplices(n-1,s), nsimplices(n,s), args...) do x
    c = hodge_diag(Val{n}, s, x)
    I, V = dual_derivative_nz(Val{ndims(s)-n}, s, x)
    V = map(I, V) do i, a
      sgn * c * a / hodge_diag(Val{n-1}, s, i)
    end
    (I, V)
  end
end

function δ(::Type{Val{n}}, s::HasDeltaSet, ::GeometricHodge, matrix_type) where n
  inv_hodge_star(n-1, s) * dual_derivative(ndims(s)-n, s) * ⋆(n, s)
end

function δ(::Type{Val{n}}, s::HasDeltaSet, ::GeometricHodge, form::AbstractVector) where n
  Vector(inv_hodge_star(n - 1, s, dual_derivative(ndims(s)-n, s, ⋆(n, s, form))))
end

""" Alias for the codifferential operator [`δ`](@ref).
"""
const codifferential = δ

""" Laplace-Beltrami operator on discrete forms.

This linear operator on primal ``n``-forms defined by ``∇² α := -δ d α``, where
[`δ`](@ref) is the codifferential and [`d`](@ref) is the exterior derivative.

!!! note

    For following texts such as Abraham-Marsden-Ratiu, we take the sign
    convention that makes the Laplace-Beltrami operator consistent with the
    Euclidean Laplace operator (the divergence of the gradient). Other authors,
    such as (Hirani 2003), take the opposite convention, which has the advantage
    of being consistent with the Laplace-de Rham operator [`Δ`](@ref).
"""
∇²(s::HasDeltaSet, x::SimplexForm{n}; kw...) where n =
  SimplexForm{n}(∇²(Val{n}, s, x.data; kw...))
@inline ∇²(n::Int, s::HasDeltaSet, args...; kw...) = ∇²(Val{n}, s, args...; kw...)

∇²(::Type{Val{n}}, s::HasDeltaSet, form::AbstractVector; kw...) where n =
  -δ(n+1, s, d(Val{n}, s, form); kw...)
∇²(::Type{Val{n}}, s::HasDeltaSet; matrix_type::Type=SparseMatrixCSC{Float64}, kw...) where n =
  -δ(n+1, s; matrix_type=matrix_type, kw...) * d(Val{n}, s, matrix_type)

""" Alias for the Laplace-Beltrami operator [`∇²`](@ref).
"""
const laplace_beltrami = ∇²

""" Laplace-de Rham operator on discrete forms.

This linear operator on primal ``n``-forms is defined by ``Δ := δ d + d δ``.
Restricted to 0-forms, it reduces to the negative of the Laplace-Beltrami
operator [`∇²`](@ref): ``Δ f = -∇² f``.
"""
Δ(s::HasDeltaSet, x::SimplexForm{n}; kw...) where n =
  SimplexForm{n}(Δ(Val{n}, s, x.data; kw...))
@inline Δ(n::Int, s::HasDeltaSet, args...; kw...) = Δ(Val{n}, s, args...; kw...)

Δ(::Type{Val{0}}, s::HasDeltaSet, form::AbstractVector; kw...) =
  δ(1, s, d(Val{0}, s, form); kw...)
Δ(::Type{Val{0}}, s::HasDeltaSet; matrix_type::Type=SparseMatrixCSC{Float64}, kw...) =
  δ(1,s; matrix_type=matrix_type, kw...) * d(Val{0},s,matrix_type)

Δ(::Type{Val{n}}, s::HasDeltaSet, form::AbstractVector; kw...) where n =
  δ(n+1, s, d(Val{n}, s, form); kw...) + d(Val{n-1}, s, δ(n, s, form; kw...))
Δ(::Type{Val{n}}, s::HasDeltaSet; matrix_type::Type=SparseMatrixCSC{Float64}, kw...) where n =
  δ(n+1,s; matrix_type=matrix_type, kw...) * d(Val{n},s,matrix_type) +
		d(Val{n-1},s,matrix_type) * δ(n,s; matrix_type=matrix_type, kw...)

Δ(::Type{Val{1}}, s::AbstractDeltaDualComplex1D, form::AbstractVector; kw...) =
  d(Val{0}, s, δ(1, s, form; kw...))
Δ(::Type{Val{1}}, s::AbstractDeltaDualComplex1D; matrix_type::Type=SparseMatrixCSC{Float64}, kw...) =
  d(Val{0},s,matrix_type) * δ(1,s; matrix_type=matrix_type, kw...)

Δ(::Type{Val{2}}, s::AbstractDeltaDualComplex2D, form::AbstractVector; kw...) =
  d(Val{1}, s, δ(2, s, form; kw...))
Δ(::Type{Val{2}}, s::AbstractDeltaDualComplex2D; matrix_type::Type=SparseMatrixCSC{Float64}, kw...) =
  d(Val{1},s,matrix_type) * δ(2,s; matrix_type=matrix_type, kw...)
""" Alias for the Laplace-de Rham operator [`Δ`](@ref).
"""
const laplace_de_rham = Δ

""" Flat operator converting vector fields to 1-forms.

A generic function for discrete flat operators. Currently only the DPP-flat from
(Hirani 2003, Definition 5.5.2) and (Desbrun et al 2005, Definition 7.3) is
implemented.

See also: the sharp operator [`♯`](@ref).
"""
♭(s::HasDeltaSet, X::DualVectorField) = EForm(♭(s, X.data, DPPFlat()))

""" Alias for the flat operator [`♭`](@ref).
"""
const flat = ♭

""" Sharp operator for converting 1-forms to vector fields.

A generic function for discrete sharp operators. Currently only the
primal-primal flat from (Hirani 2003, Definition 5.8.1 and Remark 2.7.2) is
implemented.

!!! note

    A PP-flat is also defined in (Desbrun et al 2005, Definition 7.4) but
    differs in two ways: Desbrun et al's notation suggests a *unit* normal
    vector, whereas the gradient of Hirani's primal-primal interpolation
    function is not necessarily a unit vector. More importantly, Hirani's vector
    is a normal to a different face than Desbrun et al's, with further confusion
    created by the fact that Hirani's Figure 5.7 agrees with Desbrun et al's
    description rather than his own. That being said, to the best of our
    knowledge, our implementation is the correct one and agrees with Hirani's
    description, if not his figure.

See also: the flat operator [`♭`](@ref).
"""
♯(s::HasDeltaSet, α::EForm) = PrimalVectorField(♯(s, α.data, PPSharp()))

""" Alias for the sharp operator [`♯`](@ref).
"""
const sharp = ♯

""" Wedge product of discrete forms.

The wedge product of a ``k``-form and an ``l``-form is a ``(k+l)``-form.

The DEC and related systems have several flavors of wedge product. This one is
the discrete primal-primal wedge product introduced in (Hirani, 2003, Chapter 7)
and (Desbrun et al 2005, Section 8). It depends on the geometric embedding and
requires the dual complex.
"""
∧(s::HasDeltaSet, α::SimplexForm{k}, β::SimplexForm{l}) where {k,l} =
  SimplexForm{k+l}(∧(Tuple{k,l}, s, α.data, β.data))
@inline ∧(k::Int, l::Int, s::HasDeltaSet, args...) = ∧(Tuple{k,l}, s, args...)

function ∧(::Type{Tuple{k,l}}, s::HasDeltaSet, α, β) where {k,l}
  map(simplices(k+l, s)) do x
    ∧(Tuple{k,l}, s, α, β, x)
  end
end

∧(::Type{Tuple{0,0}}, s::HasDeltaSet, f, g, x::Int) = f[x]*g[x]
∧(::Type{Tuple{k,0}}, s::HasDeltaSet, α, g, x::Int) where k =
  wedge_product_zero(Val{k}, s, g, α, x)
∧(::Type{Tuple{0,k}}, s::HasDeltaSet, f, β, x::Int) where k =
  wedge_product_zero(Val{k}, s, f, β, x)

""" Wedge product of a 0-form and a ``k``-form.
"""
function wedge_product_zero(::Type{Val{k}}, s::HasDeltaSet,
                            f, α, x::Int) where k
  subs = subsimplices(k, s, x)
  vs = primal_vertex(k, s, subs)
  coeffs = map(x′ -> dual_volume(k,s,x′), subs) / volume(k,s,x)
  dot(coeffs, f[vs]) * α[x] / factorial(k)
end

""" Alias for the wedge product operator [`∧`](@ref).
"""
const wedge_product = ∧

""" Interior product of a vector field (or 1-form) and a ``n``-form.

Specifically, this operation is the primal-dual interior product defined in
(Hirani 2003, Section 8.2) and (Desbrun et al 2005, Section 10). Thus it takes a
primal vector field (or primal 1-form) and a dual ``n``-forms and then returns a
dual ``(n-1)``-form.
"""
interior_product(s::HasDeltaSet, X♭::EForm, α::DualForm{n}; kw...) where n =
  DualForm{n-1}(interior_product_flat(Val{n}, s, X♭.data, α.data); kw...)

""" Interior product of a 1-form and a ``n``-form, yielding an ``(n-1)``-form.

Usually, the interior product is defined for vector fields; this function
assumes that the flat operator [`♭`](@ref) (not yet implemented for primal
vector fields) has already been applied to yield a 1-form.
"""
@inline interior_product_flat(n::Int, s::HasDeltaSet, args...; kw...) =
  interior_product_flat(Val{n}, s, args...; kw...)

function interior_product_flat(::Type{Val{n}}, s::HasDeltaSet,
                               X♭::AbstractVector, α::AbstractVector;
                               kw...) where n
  # TODO: Global sign `iseven(n*n′) ? +1 : -1`
  n′ = ndims(s) - n
  hodge_star(n′+1,s, wedge_product(n′,1,s, inv_hodge_star(n′,s, α; kw...), X♭); kw...)
end

""" Lie derivative of ``n``-form with respect to a vector field (or 1-form).

Specifically, this is the primal-dual Lie derivative defined in (Hirani 2003,
Section 8.4) and (Desbrun et al 2005, Section 10).
"""
ℒ(s::HasDeltaSet, X♭::EForm, α::DualForm{n}; kw...) where n =
  DualForm{n}(lie_derivative_flat(Val{n}, s, X♭, α.data; kw...))

""" Alias for Lie derivative operator [`ℒ`](@ref).
"""
const lie_derivative = ℒ

""" Lie derivative of ``n``-form with respect to a 1-form.

Assumes that the flat operator [`♭`](@ref) has already been applied to the
vector field.
"""
@inline lie_derivative_flat(n::Int, s::HasDeltaSet, args...; kw...) =
  lie_derivative_flat(Val{n}, s, args...; kw...)

function lie_derivative_flat(::Type{Val{0}}, s::HasDeltaSet,
                             X♭::AbstractVector, α::AbstractVector; kw...)
  interior_product_flat(1, s, X♭, dual_derivative(0, s, α); kw...)
end

function lie_derivative_flat(::Type{Val{1}}, s::HasDeltaSet,
                             X♭::AbstractVector, α::AbstractVector; kw...)
  interior_product_flat(2, s, X♭, dual_derivative(1, s, α); kw...) +
    dual_derivative(0, s, interior_product_flat(1, s, X♭, α; kw...))
end

function lie_derivative_flat(::Type{Val{2}}, s::HasDeltaSet,
                             X♭::AbstractVector, α::AbstractVector; kw...)
  dual_derivative(1, s, interior_product_flat(2, s, X♭, α; kw...))
end

# Euclidean geometry
####################

""" A notion of "geometric center" of a simplex.

See also: [`geometric_center`](@ref).
"""
abstract type SimplexCenter end

""" Calculate the center of simplex spanned by given points.

The first argument is a list of points and the second specifies the notion of
"center", via an instance of [`SimplexCenter`](@ref).
"""
function geometric_center end

""" Barycenter, aka centroid, of a simplex.
"""
struct Barycenter <: SimplexCenter end

function geometric_center(points, ::Barycenter)
  sum(points) / length(points)
end

""" Circumcenter, or center of circumscribed circle, of a simplex.

The circumcenter is calculated by inverting the Cayley-Menger matrix, as
explained by
[Westdendorp](https://westy31.home.xs4all.nl/Circumsphere/ncircumsphere.htm).
This method of calculation is also used in the package
[AlphaShapes.jl](https://github.com/harveydevereux/AlphaShapes.jl).
"""
struct Circumcenter <: SimplexCenter end

function geometric_center(points, ::Circumcenter)
  CM = cayley_menger(points...)
  barycentric_coords = inv(CM)[1,2:end]
  mapreduce(*, +, barycentric_coords, points)
end

""" Incenter, or center of inscribed circle, of a simplex.
"""
struct Incenter <: SimplexCenter end

function geometric_center(points, ::Incenter)
  length(points) > 2 || return geometric_center(points, Barycenter())
  face_volumes = map(i -> volume(deleteat(points, i)), eachindex(points))
  barycentric_coords = face_volumes / sum(face_volumes)
  mapreduce(*, +, barycentric_coords, points)
end

end


text:
@testset "Maze" begin

UP = CartesianIndex(-1, 0)
DOWN = CartesianIndex(1, 0)
LEFT = CartesianIndex(0, -1)
RIGHT = CartesianIndex(0, 1)
DIRECTIONS = [UP, DOWN, LEFT, RIGHT]

manhattan(a::CartesianIndex, b::CartesianIndex) = sum(abs.((b - a).I))
getmazeneighbours(maze, state) = filter(x -> (1 <= x[1] <= size(maze)[1]) && (1 <= x[2] <= size(maze)[2]) && (!maze[x]), [state + d for d in DIRECTIONS])

@testset "Maze goal reachable" begin
  maze = [0 0 1 0 0;
          0 1 0 0 0;
          0 1 0 0 1;
          0 0 0 1 1;
          1 0 1 0 0] .== 1
  start = CartesianIndex(1, 1)
  goal = CartesianIndex(1, 5)

  neighbours(state) = getmazeneighbours(maze, state)

  res = astar(neighbours, start, goal, heuristic=manhattan)
  @test res.status == :success
  @test res.path ==  CartesianIndex{2}[CartesianIndex(1, 1), CartesianIndex(2, 1), CartesianIndex(3, 1), CartesianIndex(4, 1), CartesianIndex(4, 2), CartesianIndex(4, 3), CartesianIndex(3, 3), CartesianIndex(2, 3), CartesianIndex(2, 4), CartesianIndex(1, 4), CartesianIndex(1, 5)]
  @test res.cost == 10
end

@testset "Maze Goal Unreachable" begin
  maze = [0 0 1 0 0;
          0 1 0 0 0;
          0 1 0 0 1;
          1 0 0 1 1;
          1 0 1 0 0] .== 1
  start = CartesianIndex(1, 1)
  goal = CartesianIndex(1, 5)

  neighbours(state) = getmazeneighbours(maze, state)

  res = astar(neighbours, start, goal, heuristic=manhattan)
  @test res.status == :nopath
  @test res.path == CartesianIndex{2}[CartesianIndex(1, 1), CartesianIndex(1, 2)]
  @test res.cost == 1
end

end

ext: jl


text:
@testset "Substitute" begin
    @testset "Default" begin
        test_all(ImputorTester(Substitute))
    end

    @testset "reals" begin
        # Defining our missing datasets
        a = allowmissing(1.0:1.0:20.0)
        a[[2, 3, 7]] .= missing
        fill_val = median(skipmissing(a))

        result = impute(a, Substitute())
        expected = copy(a)
        expected[[2, 3, 7]] .= fill_val

        @test result == expected
        @test result == Impute.substitute(a)
    end

    @testset "counts" begin
        a = allowmissing([1, 12, 4, 6, 2, 5, 9, 19, 24, 35, 44, 99])
        a[[2, 3, 7]] .= missing

        # We should default to taking the  median because otherwise `mode` will
        # just return `1`
        fill_val = round(Int, median(skipmissing(a)))

        result = impute(a, Substitute())
        expected = copy(a)
        expected[[2, 3, 7]] .= fill_val

        @test result == expected
        @test result == Impute.substitute(a)
    end

    @testset "ratings" begin
        # Slightly imbalanced ratings
        a = allowmissing(vcat(repeat(1:5, 5), [1, 1, 5]))
        a[[2, 3, 7]] .= missing

        # We likely want to the mode because we only have a few unique values.
        fill_val = mode(skipmissing(a))
        @test fill_val == 1

        result = impute(a, Substitute())
        expected = copy(a)
        expected[[2, 3, 7]] .= fill_val

        @test result == expected
        @test result == Impute.substitute(a)
    end

    @testset "bools" begin
        a = allowmissing(vcat(falses(14), trues(6)))
        a[[2, 3, 7]] .= missing

        # For the same reason as for ratings we should probably just use the mode.
        # Though most of the time they'll give the same answer once rounded.
        fill_val = mode(skipmissing(a))
        @test fill_val == false

        result = impute(a, Substitute())
        expected = copy(a)
        expected[[2, 3, 7]] .= fill_val

        @test result == expected
        @test result == Impute.substitute(a)
    end

    @testset "non-real" begin
        a = allowmissing(DateTime(2000, 1, 1):Day(1):DateTime(2000, 1, 20))
        a[[2, 3, 7]] .= missing

        # Median of `DateTime`s doesn't apply, so we fallback to `mode`
        fill_val = mode(skipmissing(a))

        # In this case that's just going to take the first observation it finds
        fill_val == DateTime(2000, 1, 1)

        result = impute(a, Substitute())
        expected = copy(a)
        expected[[2, 3, 7]] .= fill_val

        @test result == expected
        @test result == Impute.substitute(a)
    end

    @testset "custom statistic" begin
        # Defining our missing datasets
        a = allowmissing(1.0:1.0:20.0)
        a[[2, 3, 7]] .= missing

        # We'll do mean - 1 std for some reason :)
        μ, σ = mean_and_std(skipmissing(a))
        fill_val = μ - σ

        expected = copy(a)
        expected[[2, 3, 7]] .= fill_val
        result = Impute.substitute(
            a;
            statistic=data -> -(mean_and_std(data)...)
        )
        @test result == expected
    end

    @testset "all missing" begin
        a = [missing 2; missing 3]
        @test_logs (:warn, r"all values are missing") Impute.substitute(a; dims=2)
    end
end

@testset "WeightedSubstitute" begin
    # Don't worry about the default tests since the weights are too coupled to the data.
    @testset "reals" begin
        # Defining our missing datasets
        a = allowmissing(1.0:1.0:20.0)
        wv = eweights(20, 0.3)
        mask = trues(20)
        mask[[2, 3, 7]] .= false
        a[[2, 3, 7]] .= missing

        fill_val = median(disallowmissing(a[mask]), wv[mask])

        result = impute(a, WeightedSubstitute(; weights=wv))
        expected = copy(a)
        expected[[2, 3, 7]] .= fill_val

        @test result == expected
        @test result == Impute.wsubstitute(a; weights=wv)
    end

    @testset "counts" begin
        a = allowmissing([1, 12, 4, 6, 2, 5, 9, 19, 24, 35, 44, 99])
        wv = eweights(12, 0.3)
        mask = trues(12)
        mask[[2, 3, 7]] .= false
        a[[2, 3, 7]] .= missing

        # We should default to taking the  median because otherwise `mode` will
        # just return `1`
        fill_val = round(Int, median(disallowmissing(a[mask]), wv[mask]))

        result = impute(a, WeightedSubstitute(; weights=wv))
        expected = copy(a)
        expected[[2, 3, 7]] .= fill_val

        @test result == expected
        @test result == Impute.wsubstitute(a; weights=wv)
    end

    @testset "ratings" begin
        # Slightly imbalanced ratings
        a = allowmissing(vcat(repeat(1:5, 5), [1, 1, 5]))
        wv = eweights(28, 0.3)
        mask = trues(28)
        mask[[2, 3, 7]] .= false
        a[[2, 3, 7]] .= missing

        # We likely want to the mode because we only have a few unique values.
        fill_val = Impute._mode(a[mask], wv[mask])

        result = impute(a, WeightedSubstitute(; weights=wv))
        expected = copy(a)
        expected[[2, 3, 7]] .= fill_val

        @test result == expected
        @test result == Impute.wsubstitute(a; weights=wv)
    end

    @testset "bools" begin
        a = allowmissing(vcat(falses(14), trues(6)))
        wv = eweights(20, 0.3)
        mask = trues(20)
        mask[[2, 3, 7]] .= false
        a[[2, 3, 7]] .= missing

        # For the same reason as for ratings we should probably just use the mode.
        # Though most of the time they'll give the same answer once rounded.
        fill_val = Impute._mode(a[mask], wv[mask])

        result = impute(a, WeightedSubstitute(; weights=wv))
        expected = copy(a)
        expected[[2, 3, 7]] .= fill_val

        @test result == expected
        @test result == Impute.wsubstitute(a; weights=wv)
    end

    @testset "non-real" begin
        a = allowmissing(DateTime(2000, 1, 1):Day(1):DateTime(2000, 1, 20))
        wv = eweights(20, 0.3)
        mask = trues(20)
        mask[[2, 3, 7]] .= false
        a[[2, 3, 7]] .= missing

        # Median of `DateTime`s doesn't apply, so we fallback to `mode`
        fill_val = Impute._mode(a[mask], wv[mask])

        result = impute(a, WeightedSubstitute(; weights=wv))
        expected = copy(a)
        expected[[2, 3, 7]] .= fill_val

        @test result == expected
        @test result == Impute.wsubstitute(a; weights=wv)
    end

    @testset "custom statistic" begin
        # Defining our missing datasets
        a = allowmissing(1.0:1.0:20.0)
        wv = eweights(20, 0.3)
        mask = trues(20)
        mask[[2, 3, 7]] .= false
        a[[2, 3, 7]] .= missing

        # We'll do mean - 1 std for some reason :)
        μ, σ = mean_and_std(disallowmissing(a[mask]), wv[mask])
        fill_val = μ - σ

        expected = copy(a)
        expected[[2, 3, 7]] .= fill_val
        result = Impute.wsubstitute(
            a;
            statistic=(data, weights) -> -(mean_and_std(data, weights)...),
            weights=wv
        )
        @test result == expected
    end

    @testset "all missing" begin
        a = [missing 2; missing 3]
        @test_logs (:warn, r"all values are missing") Impute.substitute(a; dims=2)
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    treewalk(f, tree::GitTree, post::Bool=false)

Traverse the entries in `tree` and its subtrees in post or pre order. Preorder
means beginning at the root and then traversing the leftmost subtree (and
recursively on down through that subtree's leftmost subtrees) and moving right
through the subtrees. Postorder means beginning at the bottom of the leftmost
subtree, traversing upwards through it, then traversing the next right subtree
(again beginning at the bottom) and finally visiting the tree root last of all.

The function parameter `f` should have following signature:

    (String, GitTreeEntry) -> Cint

A negative value returned from `f` stops the tree walk. A positive value means
that the entry will be skipped if `post` is `false`.
"""
function treewalk(f, tree::GitTree, post::Bool = false)
    ensure_initialized()
    # NOTE: don't use @check/GitError directly, because the code can be arbitrary
    payload = Any[ tree, f ]
    cbf = @cfunction(function treewalk_callback(root_cstr::Cstring, entry_ptr::Ptr{Cvoid}, payload::Vector{Any})::Cint
            # decode arguments
            root = unsafe_string(root_cstr)
            tree = payload[1]::GitTree
            f = payload[2]
            entry = GitTreeEntry(tree, entry_ptr, false)
            return f(root, entry)
        end, Cint, (Cstring, Ptr{Cvoid}, Ref{Vector{Any}}))
    err = ccall((:git_tree_walk, :libgit2), Cint,
                (Ptr{Cvoid}, Cint, Ptr{Cvoid}, Any),
                tree.ptr, post, cbf, payload)
    if err < 0
        err_class, _ = Error.last_error()
        if err_class != Error.Callback
            # now we now the code is valid
            throw(GitError(err))
        end
    end
    nothing
end

repository(tree::GitTree) = tree.owner
repository(te::GitTreeEntry) = repository(te.owner)

"""
    filename(te::GitTreeEntry)

Return the filename of the object on disk to which `te` refers.
"""
function filename(te::GitTreeEntry)
    ensure_initialized()
    str = ccall((:git_tree_entry_name, :libgit2), Cstring, (Ptr{Cvoid},), te.ptr)
    str != C_NULL && return unsafe_string(str)
    return nothing
end

"""
    filemode(te::GitTreeEntry) -> Cint

Return the UNIX filemode of the object on disk to which `te` refers as an integer.
"""
function filemode(te::GitTreeEntry)
    ensure_initialized()
    return ccall((:git_tree_entry_filemode, :libgit2), Cint, (Ptr{Cvoid},), te.ptr)
end

"""
    entrytype(te::GitTreeEntry)

Return the type of the object to which `te` refers. The result will be
one of the types which [`objtype`](@ref) returns, e.g. a `GitTree` or `GitBlob`.
"""
function entrytype(te::GitTreeEntry)
    ensure_initialized()
    otype = ccall((:git_tree_entry_type, :libgit2), Cint, (Ptr{Cvoid},), te.ptr)
    return objtype(Consts.OBJECT(otype))
end

"""
    entryid(te::GitTreeEntry)

Return the [`GitHash`](@ref) of the object to which `te` refers.
"""
function entryid(te::GitTreeEntry)
    ensure_initialized()
    GC.@preserve te begin
        oid_ptr = ccall((:git_tree_entry_id, :libgit2), Ptr{UInt8}, (Ptr{Cvoid},), te.ptr)
        oid = GitHash(oid_ptr)
    end
    return oid
end

function count(tree::GitTree)
    ensure_initialized()
    return ccall((:git_tree_entrycount, :libgit2), Csize_t, (Ptr{Cvoid},), tree.ptr)
end

function Base.getindex(tree::GitTree, i::Integer)
    if i < 1 || i > count(tree)
        throw(BoundsError(tree, i))
    end
    ensure_initialized()
    te_ptr = ccall((:git_tree_entry_byindex, :libgit2),
                   Ptr{Cvoid},
                   (Ptr{Cvoid}, Csize_t), tree.ptr, i-1)
    return GitTreeEntry(tree, te_ptr, false)
end

"""
    (::Type{T})(te::GitTreeEntry) where T<:GitObject

Get the git object to which `te` refers and return it as its actual type (the type
[`entrytype`](@ref) would show), for instance a `GitBlob` or `GitTag`.

# Examples
```julia
tree = LibGit2.GitTree(repo, "HEAD^{tree}")
tree_entry = tree[1]
blob = LibGit2.GitBlob(tree_entry)
```
"""
function GitObject(e::GitTreeEntry) end
function (::Type{T})(te::GitTreeEntry) where T<:GitObject
    ensure_initialized()
    repo = repository(te)
    obj_ptr_ptr = Ref{Ptr{Cvoid}}(C_NULL)
    @check ccall((:git_tree_entry_to_object, :libgit2), Cint,
                  (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Ref{Nothing}),
                   obj_ptr_ptr, repo.ptr, te.ptr)
    return T(repo, obj_ptr_ptr[])
end

function Base.show(io::IO, te::GitTreeEntry)
    println(io, "GitTreeEntry:")
    println(io, "Entry name: ", filename(te))
    println(io, "Entry type: ", entrytype(te))
    println(io, "Entry OID: ", entryid(te))
end

function Base.show(io::IO, tree::GitTree)
    println(io, "GitTree:")
    println(io, "Owner: ", repository(tree))
    println(io, "Number of entries: ", count(tree))
end

"""
    getindex(tree::GitTree, target::AbstractString) -> GitObject

Look up `target` path in the `tree`, returning a [`GitObject`](@ref) (a [`GitBlob`](@ref) in
the case of a file, or another [`GitTree`](@ref) if looking up a directory).

# Examples
```julia
tree = LibGit2.GitTree(repo, "HEAD^{tree}")
readme = tree["README.md"]
subtree = tree["test"]
runtests = subtree["runtests.jl"]
```
"""
function Base.getindex(tree::GitTree, target::AbstractString)
    if basename(target) == ""
        # get rid of any trailing separator
        target = dirname(target)
    end
    if target == "" || target == "/"
        return tree
    end

    local oid = nothing
    function _getindex_callback(root::String, entry::GitTreeEntry)::Cint
        path = joinpath(root, filename(entry))
        if path == target
            # we found the target, save the oid and stop the walk
            oid = entryid(entry)
            # workaround for issue: https://github.com/libgit2/libgit2/issues/4693
            ensure_initialized()
            ccall((:giterr_set_str, :libgit2), Cvoid,
                  (Cint, Cstring), Cint(Error.Callback),
                  "git_tree_walk callback returned -1")
            return -1
        elseif entrytype(entry) == GitTree && !startswith(target, path)
            # this subtree isn't relevant, so skip it
            return 1
        end
        return 0
    end
    treewalk(_getindex_callback, tree)
    oid === nothing && throw(KeyError(target))
    return GitObject(repository(tree), oid)
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
#######################################

"""
    function groebner(
                polynomials::Vector{Poly};
                reduced::Bool=true,
                ordering::Symbol=:input,
                randomized::Bool=true,
                rng::Rng=Random.MersenneTwister(42),
                loglevel::Logging.LogLevel=Logging.Warn
    ) where {Poly, Rng<:Random.AbstractRNG}

Computes a Groebner basis of the ideal generated by `polynomials`.

If `reduced` is set, returns the reduced basis, which is **unique**.

Uses the ordering on `polynomials` by default.
If `ordering` is explicitly specialized, it takes precedence.
*(On the fly ordering change not implemented yet :D)*

The algorithm is randomized by default, but
this can be changed with the `randomized` param.
*(Derandomized version not implemented yet :D)*

"""
function groebner(
            polynomials::Vector{Poly};
            reduced::Bool=true,
            ordering::Symbol=:input,
            randomized::Bool=true,
            rng::Rng=Random.MersenneTwister(42),
            loglevel::Logging.LogLevel=Logging.Warn
            ) where {Poly, Rng<:Random.AbstractRNG}

    #= set the logger =#
    prev_logger = Logging.global_logger(ConsoleLogger(stderr, loglevel))

    #= extract ring information, exponents and coefficients
       from input polynomials =#
    # Copies input, so that polys would not be changed itself.
    ring, exps, coeffs = convert_to_internal(polynomials, ordering)

    #= compute the groebner basis =#
    if ring.ch != 0
        # if finite field
        # Always returns UInt coefficients #
        bexps, bcoeffs = groebner_ff(ring, exps, coeffs, reduced, rng)
    else
        # if rational coefficients
        # Always returns rational coefficients #
        bexps, bcoeffs = groebner_qq(ring, exps, coeffs, reduced, randomized, rng)
    end

    #=
    Assuming ordering of `bexps` here matches `ring.ord`
    =#

    #= revert logger =#
    Logging.global_logger(prev_logger)

    # ring contains ordering of computation, it is the requested ordering
    #= convert result back to representation of input =#
    convert_to_output(ring, polynomials, bexps, bcoeffs)
end

#######################################
# Finite field groebner

function groebner_ff(
            ring::PolyRing,
            exps::Vector{Vector{Vector{UInt16}}},
            coeffs::Vector{Vector{UInt64}},
            reduced::Bool,
            rng::Rng) where {Rng<:Random.AbstractRNG}
    # specialize on ordering (not yet)
    # groebner_ff(ring, exps, coeffs, reduced, rng, Val(ring.ord))
    f4(ring, exps, coeffs, rng, reduced)
end

#######################################
# Rational field groebner

function modular_f4_step(
            ring::PolyRing,
            exps::Vector{Vector{Vector{UInt16}}},
            coeffs::Vector{Vector{UInt64}},
            rng::Rng,
            reduced::Bool
            ) where {Rng<:Random.AbstractRNG}

    f4(ring, exps, coeffs, rng, reduced)
end

function groebner_qq(
            ring::PolyRing,
            exps::Vector{Vector{Vector{UInt16}}},
            coeffs::Vector{Vector{Rational{BigInt}}},
            reduced::Bool,
            randomized::Bool,
            rng::Rng,
            ) where {Rng<:Random.AbstractRNG}

    # scale coefficients to integer ring inplace
    coeffs_zz = scale_denominators!(coeffs)

    gbcoeffs_accum = Vector{Vector{BigInt}}(undef, 0)
    gbexps = Vector{Vector{Vector{UInt16}}}(undef, 0)
    gbcoeffs_qq = Vector{Vector{Rational{BigInt}}}(undef, 0)

    prime::Int64 = 1
    modulo = BigInt(1)
    moduli = Int[]

    i = 1
    while true
        # lucky reduction prime
        prime  = nextluckyprime(coeffs_zz, prime)
        push!(moduli, prime)
        @info "$i: selected lucky prime $prime"

        # compute the image of coeffs_zz in GF(prime),
        # by coercing each coefficient into the finite field
        ring_ff, coeffs_ff = reduce_modulo(coeffs_zz, ring, prime)

        # groebner basis over finite field ideal
        # TODO: gbexps can be only computed once
        @info "Computing Groebner basis"
        gbexps, gbcoeffs_ff = modular_f4_step(
                                    ring_ff, exps, coeffs_ff,
                                    rng, reduced)

        # TODO: add majority rule based choice

        # reconstruct basis coeffs into integers
        # from the previously accumulated basis and the new one,
        @info "CRT modulo ($modulo, $(ring_ff.ch))"
        gbcoeffs_zz, modulo = reconstruct_crt!(
                            gbcoeffs_accum, modulo,
                            gbcoeffs_ff, ring_ff.ch)

        gbcoeffs_accum = gbcoeffs_zz

        # try to reconstruct basis coeffs from integers
        # into rationals
        @info "Reconstructing modulo $modulo"
        gbcoeffs_qq = reconstruct_modulo(gbcoeffs_zz, modulo)

        # run correctness checks to assure the reconstrction is correct
        if correctness_check(ring, exps, coeffs_zz, gbexps,
                                gbcoeffs_qq, moduli, randomized, rng)
            break
        end

        # not correct, goto next prime
        i += 1
        if i > 10000
            @error "Something probably went wrong in groebner.."
            return
        end
    end

    gbexps, gbcoeffs_qq
end

#######################################
# Generic isgroebner

"""
    function isgroebner(
                polynomials::Vector{Poly};
                ordering::Symbol=:input,
                randomized::Bool=true,
                rng::Rng=Random.MersenneTwister(42),
                loglevel::Logging.LogLevel=Logging.Warn
    ) where {Poly, Rng<:Random.AbstractRNG}

Checks if `polynomials` forms a Groebner basis.

Uses the ordering on `polynomials` by default.
If `ordering` is explicitly specialized, it takes precedence.
*(On the fly ordering change not implemented yet :D)*

The algorithm is randomized by default, but
this can be changed with the `randomized` param.
*(Derandomized version not implemented yet :D)*

"""
function isgroebner(
            polynomials::Vector{Poly};
            ordering=:input,
            randomized::Bool=true,
            rng::Rng=Random.MersenneTwister(42),
            loglevel::LogLevel=Logging.Warn
    ) where {Poly, Rng<:Random.AbstractRNG}

    #= set the logger =#
    prev_logger = Logging.global_logger(ConsoleLogger(stderr, loglevel))

    #= extract ring information, exponents and coefficients
       from input polynomials =#
    # Copies input, so that polys would not be changed itself.
    ring, exps, coeffs = convert_to_internal(polynomials, ordering)

    #= compute the groebner basis =#
    if ring.ch != 0
        # if finite field
        # Always returns UInt coefficients #
        flag = isgroebner_ff(ring, exps, coeffs, rng)
    else
        # if rational coefficients
        # Always returns rational coefficients #
        flag = isgroebner_qq(ring, exps, coeffs, randomized, rng)
    end

    #=
    Assuming ordering of `bexps` here matches `ring.ord`
    =#

    #= revert logger =#
    Logging.global_logger(prev_logger)

    flag
end

# UWU!
function isgroebner_ff(
            ring::PolyRing,
            exps::Vector{Vector{Vector{UInt16}}},
            coeffs::Vector{Vector{UInt64}},
            rng)

    isgroebner_f4(ring, exps, coeffs, rng)
end

# TODO
function isgroebner_qq(
            ring::PolyRing,
            exps::Vector{Vector{Vector{UInt16}}},
            coeffs::Vector{Vector{Rational{BigInt}}},
            randomized::Bool,
            rng)

    if randomized
        coeffs_zz = scale_denominators!(coeffs)
        goodprime = nextgoodprime(coeffs_zz, Int[], 2^30 + 3)
        ring_ff, coeffs_ff = reduce_modulo(coeffs_zz, ring, goodprime)
        isgroebner_f4(ring_ff, exps, coeffs_ff, rng)
    else
        error("Sorry, not randomized version is not implemented yet.")
    end
end

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


text:
#
println("Af) Test of the frequency-independent Breit interaction for an internally generated chlorine-like multiplet.")
#
wa = Atomic.Computation("xx",  Nuclear.Model(54.);
                        configs=[Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6"), Configuration("[Ne] 3s^2 3p^4 3d")],
                        asfSettings=AsfSettings(true, false, "meanDFS", "hydrogenic", "xxx", [1],    40, 1.0e-6, JAC.Subshell[],
                                                true, false, NoneQed(), "yyy", LSjjSettings(true),
                                                false, [1,2,3,4], false, JAC.LevelSymmetry[] )  )

wb = perform(wa)


text:
"""
    beamform(s::AbstractVector{Seis.Trace}, t1, t2, sx1, sx2, sy1, sy2, ds; maxima=1, method=:linear, n=1) -> ::BeamformGrid

Compute normalised beam power at each horizontal slowness point
`sx1:ds:sx2` × `sy1:ds:sy2` s/°, between times `t1` and `t2` s.  The slowness
grid is spaced by `ds` s/°.

The `maxima` maximum beam power points are set in the returned `BeamformGrid`.

    beamform(s, t1, t2, smax, ds; maxima=1, method=:linear, n=1) -> ::BeamformGrid

Use a slowness grid of `-smax:ds:smax` in both dimensions.

# Keyword arguments

- `maxima=1`: The top `maxima` points with the highest beam power are saved to the
  `BeamformGrid` returned.
- `method=:linear`: Set the stacking method to use.  See [`stack`](@ref) for available
  stacking methods.
- `n`: Power of the nth-root or phase-weighted stack.
- `wavefront=:plane`: Geometry of problem to consider.  If `:plane`, consider
  a plane wave moving across the array.  If `circle`, assumes a circular wavefront
  with backazimuth and absolute slowness determined by the slowness grid points.
"""
function beamform(s::TraceArray{T}, t1, t2, sx1, sx2, sy1, sy2, ds;
        maxima=1, method=:linear, n=nothing, wavefront=:plane) where T
    sx = s_per_km.(sx1:ds:sx2)
    sy = s_per_km.(sy1:ds:sy2)
    nx, ny = length.((sx, sy))
    beam_power = Array{T}(undef, nx, ny)
    lon, lat = s.sta.lon, s.sta.lat
    evt = first(s).evt
    x, y, z, mean_lon, mean_lat = array_geometry(lon, lat, zeros(length(s)))
    Δ̄ = wavefront === :circle ? delta(mean_lon, mean_lat, evt.lon, evt.lat) : nothing
    for (iy, isy) in enumerate(sy), (ix, isx) in enumerate(sx)
        align = shift_times(x, y, isx, isy, wavefront;
            lon=lon, lat=lat, mean_lon=mean_lon, mean_lat=mean_lat, mean_dist=Δ̄)
        S = stack(s, (t1,t2), align; method=method, n=n)
        beam_power[ix,iy] = sum(x->x^2, Seis.trace(S))
    # end
    end
    sx, sy = s_per_degree.(sx), s_per_degree.(sy)
    sx_max, sy_max, p, β = find_maxima(sx, sy, beam_power, maxima)
    BeamformGrid{T,eltype(s)}(mean_lon, mean_lat, x, y,
        collect(sx1:ds:sx2), collect(sy1:ds:sy2),
        beam_power./maximum(beam_power), maxima,
        sx_max, sy_max, p, β, s, method, n, t1, t2)
end
beamform(s::TraceArray, t1, t2, smax, ds; kwargs...) =
    beamform(s, t1, t2, -smax, smax, -smax, smax, ds; kwargs...)

"""
    vespagram(traces::AbstractVector{<:Seis.Trace}, t1, t2, s1, s2, ds; kwargs...) -> vesp::VespaGrid

Create a vespagram for a set of `traces`, between times `t1` and `t2` seconds,
and for the range of slownesses from `s1` to `s2` s/° in steps of `ds`
s/°.

# Keyword arguments
- `maxima=1`: The top `maxima` points with the highest beam power are saved to the
  `VespaGrid` returned.
- `method=:linear`: Set the stacking method to use.  See [`stack`](@ref) for available
  stacking methods.
- `n`: Power of the nth-root or phase-weighted stack.
- `wavefront=:plane`: Geometry of problem to consider.  If `:plane`, consider
  a plane wave moving across the array.  If `circle`, assumes a circular wavefront
  with backazimuth and absolute slowness determined by the slowness grid points.
- `envelope=false`: Whether or not to stack the envelopes of traces
"""
function vespagram(s::TraceArray{T}, t1, t2, s1, s2, ds;
        maxima=1, method=:linear, n=nothing, envelope=false, wavefront=:plane) where T
    slow = s_per_km.(s1:ds:s2)
    t = t1:first(s).delta:t2
    nslow = length(slow)
    nt = length(t)
    vespa = Array{T}(undef, nt, nslow)
    lon, lat = s.sta.lon, s.sta.lat
    evt = first(s).evt
    x, y, z, mean_lon, mean_lat = array_geometry(lon, lat, zeros(length(s)))
    Δ̄ = wavefront === :circle ? delta(mean_lon, mean_lat, evt.lon, evt.lat) : nothing
    ᾱ = Geodesics.azimuth(mean_lon, mean_lat, evt.lon, evt.lat) #+ 180
    for (i, islow) in enumerate(slow)
        sx, sy = islow.*sincos(deg2rad(ᾱ))
        align = shift_times(x, y, sx, sy, wavefront; lon=lon, lat=lat,
            mean_lon=mean_lon, mean_lat=mean_lat, mean_dist=Δ̄)
        S = stack(s, (t1, t2), align; method=method, n=n)
        envelope && (Seis.envelope!(S))
        vespa[:,i] .= Seis.trace(S)
    end
    slow = s_per_degree(slow)
    t_max, slow_max, _, _ = find_maxima(t, slow, vespa, maxima)
    VespaGrid{T,eltype(s)}(mean_lon, mean_lat, x, y, t, slow, vespa,
        maxima, t_max, slow_max, s, method, n, wavefront, envelope)
end

"""
    shift_times(x, y, sx, sy, wavefront=:plane; kwargs...) -> shifts::Vector

Compute a set of shift times in s, `shifts`, which represent the relative
time at which a wavefront travelling at x-slowness `sx`  and y-slowness `sy` s/km
occurs at each point in space `x` and `y` km.  By default, a plane wave is assumed;
set `wavefront=:circle` to use a circular wavefront.

# Keyword arguments
## Required if `wavefront == :circle`:
- `lon, lat`: The longitude and latitude (°) of each station.
- `mean_lon, mean_lat`: Reference longitude and latitude (°) of the array
- `mean_dist`: Distance from event to reference point
## Optional if `wavefront == :circle`:
- `earth_radius = $(R_EARTH_KM_DEFAULT)`: Planetary radius in km
"""
function shift_times(x, y, sx, sy, wavefront=:plane;
                     lon=nothing, lat=nothing, mean_lon=nothing, mean_lat=nothing,
                     mean_dist=nothing, earth_radius=R_EARTH_KM_DEFAULT)
    if wavefront === :plane
        -(x.*sx .+ y.*sy)
    elseif wavefront === :circle
        any(x -> x===nothing, (lon, lat, mean_lon, mean_lat, mean_dist)) &&
            throw(ArgumentError("lon, lat, mean_lon, mean_lat and mean_dist " *
                                "must be given when `wavefront` is `:circle`"))
        s = sqrt(sx^2 + sy^2) # s/km
        β = rad2deg(atan(sx, sy))
        # Location of new 'event' at this slowness and backazimuth
        new_evt_lon, new_evt_lat = geodesic_endpoint(mean_lon, mean_lat, β, mean_dist)
        # Distances to new 'event' from each station, in °
        Δs = delta.(new_evt_lon, new_evt_lat, lon, lat)
        # Distance from array reference to new 'event', in °
        Δ = delta(new_evt_lon, new_evt_lat, mean_lon, mean_lat)
        # Station delays are relative times compared to array centre
        (Δs .- Δ).*s_per_degree(s, earth_radius)
    else
        throw(ArgumentError("unknown `wavefront` type"))
    end
end

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


text:
using DataTypesBasic
using Suppressor

@test map(Identity(3)) do x
  x + 4
end == Identity(3 + 4)

@test repr(Identity(4)) == "Identity(4)"

@test eltype(Identity("hello")) == String
@test eltype(Identity{Int}) == Int
@test eltype(Identity) == Any

@test collect(Identity(1)) == [1]
@test length(Identity([1,2,3,4])) == 1

@test get(Identity(42)) == 42
@test Identity(42)[] == 42

@test isidentity(Identity(4))
@test !isidentity(Const(3))

@test Iterators.flatten(Identity(Identity(5))) == Identity(5)
# special support for Const:
@test Iterators.flatten(Identity(Const(5))) == Const(5)

@test @capture_out(foreach(print, Identity(5))) == "5"
@test @capture_out(for i in Identity(5)
  print(5)
end) == "5"

@test convert(Identity{Float32}, Identity(1)) == Identity(Float32(1))
@test convert(Identity{Number}, Identity(1)) == Identity(1)

@test promote_type(Identity{Int}, Identity{Number}) == Identity{Number}
@test promote_type(Identity{Int}, Identity{String}) == Identity
@test promote_type(Identity{Int}, Identity) == Identity
@test promote_type(Identity, Identity) == Identity


text:
mutable struct MCSimulation{RSG <: AbstractRandomSequenceGenerator, T <: MCTrait, P <: AbstractPathPricer, S <: StochasticProcess1D}
  antitheticVariate::Bool
  controlVariate::Bool
  rsg::RSG
  mcTrait::T
  mcModel::MonteCarloModel{P, RSG, S}

  # MCSimulation(antitheticVariate::Bool,
  #             controlVariate::Bool,
  #             rsg::RSG,
  #             mcTrait::T,
  #             mcModel::MonteCarloModel{P, RSG, S}) where {RSG, T, P, S} =
  #             new{RSG, T, P, S}(antitheticVariate, controlVariate, rsg, mcTrait, mcModel)
end

function MCSimulation(pe::MCEngine{S, RSG},
                      controlVariate::Bool,
                      inst::Instrument,
                      mcTrait::T) where {S <: StochasticProcess1D, RSG <: AbstractRandomSequenceGenerator, T <: MCTrait}
  # create model
  model = MonteCarloModel(path_generator(pe, inst), path_pricer(pe, inst), gen_RiskStatistics(), pe.antitheticVariate)
  return MCSimulation{RSG, T, typeof(model.pathPricer), S}(pe.antitheticVariate, controlVariate, pe.rsg, mcTrait, model)
end

max_error(err::Float64) = err

function get_value!(mcsim::MCSimulation, tolerance::Float64, maxSamples::Int = typemax(Int), minSamples::Int = 1023)
  sampleNumber = size(mcsim.mcModel.sampleAccumulator.samplesMatrix)[1]
  if sampleNumber < minSamples
    idxNum = sampleNumber + 1
    add_samples!(mcsim.mcModel, minSamples - sampleNumber, idxNum) # add an index so we can add to this
    sampleNumber = size(mcsim.mcModel.sampleAccumulator.samplesMatrix)[1]
  end

  err = error_estimate(mcsim.mcModel.sampleAccumulator)
  while max_error(err) > tolerance
    sampleNumber < maxSamples || error("max number of samples reached while error still above tolerance")
    # conservative estimate of how many samples are needed
    order = max_error(err * err) / tolerance / tolerance
    nextBatch = round(Int, floor(max(sampleNumber * order * 0.8 - sampleNumber, minSamples)))
    nextBatch = min(nextBatch, maxSamples - sampleNumber)
    idxNum = sampleNumber + 1
    sampleNumber += nextBatch
    add_samples!(mcsim.mcModel, nextBatch, idxNum)
    err = error_estimate(mcsim.mcModel.sampleAccumulator)
  end
  return stats_mean(mcsim.mcModel.sampleAccumulator)
end

function get_value_with_samples!(mcsim::MCSimulation, samples::Int)
  sampleNumber = size(mcsim.mcModel.sampleAccumulator.samplesMatrix)[1]

  samples >= sampleNumber || error("number of already simulated samples greater than requested samples")

  idxNum = sampleNumber + 1
  add_samples!(mcsim.mcModel, samples - sampleNumber, idxNum)

  return stats_mean(mcsim.mcModel.sampleAccumulator)
end

function _calculate!(mcsim::MCSimulation, requiredTolerance::Float64, requiredSamples::Int, maxSamples::Int)
  # TODO check if control variate
  # mcsim.mcModel = MonteCarloModel(path_generator(pe, inst), path_pricer(pe, inst), gen_RiskStatistics(), mcsim.antitheticVariate)

  if requiredTolerance != -1.0
    if maxSamples != -1
      get_value!(mcsim, requiredTolerance, maxSamples)
    else
      get_value!(mcsim, requiredTolerance)
    end
  else
    get_value_with_samples!(mcsim, requiredSamples)
  end

  return mcsim
end


text:
using BenchmarkTools

function expensive()
    X = randn(10_000, 1_000)
    return inv(X'X)
end

function go_serial()
    A = expensive()
    B = expensive()

    return A + B
end

function go_threads()
    A = Threads.@spawn expensive()
    B = Threads.@spawn expensive()

    # Wait for both to complete
    return fetch(A) + fetch(B)
end

display(@benchmark go_serial())
display(@benchmark go_threads())


text:
function Arc(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return Arc((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function check_point(obj::Arc, arg0::jdouble)
    return jcall(obj, "checkPoint", Region_Location, (jdouble,), arg0)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_barycenter(obj::Arc)
    return jcall(obj, "getBarycenter", jdouble, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_inf(obj::Arc)
    return jcall(obj, "getInf", jdouble, ())
end

function get_offset(obj::Arc, arg0::S1Point)
    return jcall(obj, "getOffset", jdouble, (S1Point,), arg0)
end

function get_offset(obj::Arc, arg0::jdouble)
    return jcall(obj, "getOffset", jdouble, (jdouble,), arg0)
end

function get_size(obj::Arc)
    return jcall(obj, "getSize", jdouble, ())
end

function get_sup(obj::Arc)
    return jcall(obj, "getSup", jdouble, ())
end

function get_tolerance(obj::Arc)
    return jcall(obj, "getTolerance", jdouble, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
using Gadfly

n=16
pts = repeat(1:n,inner=(1,n))
plot(x=pts, y=pts', color=string.(1:n*n), Theme(key_position=:none))


text:
punctuation = ['!','"','#','$','%','&','\'','(',')','*','+',',',
'-','.','/',':',';','<','=','>','?','@','[','\\',']','^','_','`',
'{','|','}','~']

function clean(string)
  string = Base.lowercase(string)
  sarray = convert(Array{typeof(string)}, Base.split(string))
  x = 1

  while x <= length(sarray)
    sarray[x] = Base.strip(sarray[x], punctuation)
    x += 1
  end

  return sarray
end

function parse_xml(doc_path)
  xdoc = parse_file(doc_path)
  xroot = root(xdoc)
  ces = get_elements_by_tagname(xroot, "EEBO")
  children = collect(child_elements(ces[1]))
  if(name(children[2]) == "GROUP")
    ces = get_elements_by_tagname(ces[1], "GROUP")
  else
    ces = get_elements_by_tagname(ces[1], "TEXT")
  end
  body = content(ces[1])
  text = string(body)

  free(xdoc)

  return text
end

function get_metadata(doc_path)
  metadata = Array(String,4)
  xdoc = parse_file(doc_path)
  xroot = root(xdoc)
  ces = get_elements_by_tagname(xroot, "HEADER")
  filedesc = get_elements_by_tagname(ces[1], "FILEDESC")
  source = get_elements_by_tagname(filedesc[1], "SOURCEDESC")
  bib = get_elements_by_tagname(source[1], "BIBLFULL")
  titlestmt = get_elements_by_tagname(bib[1], "TITLESTMT")
  author = get_elements_by_tagname(titlestmt[1], "AUTHOR")

  if length(author) == 0
    metadata[1] = "NA"
  else
    author = content(author[1])
    author = string(author)
    metadata[1] = author
  end

  title = get_elements_by_tagname(titlestmt[1], "TITLE")
  if length(title) == 0
    metadata[2] = "NA"
  else
    title = content(title[1])
    title = string(title)
    metadata[2]  = title
  end

  date = get_elements_by_tagname(bib[1], "PUBLICATIONSTMT")
  date = get_elements_by_tagname(date[1], "DATE")
  if length(date) == 0
    metadata[3] = "NA"
  else
    date = content(date[1])
    date = string(date)
    metadata[3] = date
  end

  profile = get_elements_by_tagname(ces[1], "PROFILEDESC")
  langusage = get_elements_by_tagname(profile[1], "LANGUSAGE")
  lang = get_elements_by_tagname(langusage[1], "LANGUAGE")
  if length(lang) == 0
    metadata[4] = "NA"
  else
    lang = content(lang[1])
    lang = string(lang)
    metadata[4] = lang
  end

  free(xdoc)

  return metadata
end

function load_featurevector(path)
  if isfile(path)
    fv = FeatureVector()
    parsed_file = parse_xml(path)
    clean_text = clean(parsed_file)
    for word in clean_text
      if word in keys(fv)
        fv[word] += 1
      else
        setindex!(fv,1,word)
      end
    end
    fv[""] = 0
    fv.mdata = get_metadata(path)
    return fv
  else
    Base.warn("$path is not to a valid file. Please check the path and try again.")
  end
end

function load_cluster(path)
  if isdir(path)
    curdir = pwd()
    cd(path)
    files = readdir(pwd())
    cl = Cluster()
    for file in files
      if endswith(file,"xml")
        fv = load_featurevector(file)
        cl[file] = fv
      else
        Base.warn("$file is not a valid file. Please enter a path to a directory of .xml files.")
      end
    end
    cd(curdir)
    return cl
  else
    Base.warn("$path is not to a valid directory. Please check the path and try again.")
  end
end

function load_dataset(path)
  if isdir(path)
    curdir = pwd()
    cd(path)
    dirs = readdir(pwd())
    ds = DataSet()
    for dir in dirs
      cl = Cluster()
      if isdir(dir)
        cl = load_cluster(dir)
        ds[dir] = cl
      else
        Base.warn("$dir is not a valid Directory. This will not be added.")
      end
    end
    cd(curdir)
    return ds
  else
    Base.warn("$path is not to a valid directory. Please check the path and try again.")
  end
end


text:
export d2u, j2md, md2j, parsetimewin, timestamp, u2d, timespec

# =====================================================================
# Time stamping
function tstr(t::DateTime)
  Y, M, D, h, m, s, μ = year(t), month(t), day(t), hour(t), minute(t), second(t), millisecond(t)
  Y = lpad(Y, 4, "0")
  M = lpad(M, 2, "0")
  D = lpad(D, 2, "0")
  h = lpad(h, 2, "0")
  m = lpad(m, 2, "0")
  s = lpad(s, 2, "0")
  return string(Y, "-", M, "-", D, "T", h, ":", m, ":", s)
end

"Alias to Dates.unix2datetime"
u2d(k::Real) = Dates.unix2datetime(k)
"Alias to Dates.datetime2unix"
d2u(k::DateTime) = Dates.datetime2unix(k)

@doc """
    timestamp()

Return current time formatted YYYY-mm-ddTHH:MM:SS.
""" timestamp
timestamp() = tstr(Dates.unix2datetime(time()))
timestamp(t::DateTime) = tstr(t)
timestamp(t::Real) = tstr(u2d(t))
timestamp(t::String) = tstr(Dates.DateTime(t))
tnote(s::String) = string(timestamp(), " ¦ ", s)

# =====================================================================
# Time and date converters
"""
    m,d = j2md(y,j)

Convert Julian day j of year y to month m, day d
"""
function j2md(y::T, j::T) where T<:Integer
  if T != Int32
    y = Int32(y)
    j = Int32(j)
  end
  z = zero(Int32)
  o = one(Int32)
  m = z
  d = o
  if j > Int32(31)
    leapyear = ((j > 59) && ((y % Int32(400) == z) || (y % Int32(4) == z && y % Int32(100) != z)))
    while j > z
      d = j
      m += o
      j -= days_per_month[m]
      if leapyear && m == 2
        j -= o
      end
    end
  else
    m = o
    d = j
  end
  return m,d
end

"""
    j = md2j(y, m, d)

Convert month `m`, day `d` of year `y` to Julian day (day of year)
"""
function md2j(y::T, m::T, d::T) where T<:Integer
  j = zero(Int32)
  if T != Int32
    y = Int32(y)
    m = Int32(m)
    d = Int32(d)
  end
  z = zero(Int32)
  i = one(Int8)
  while i < m
    j = j+getindex(days_per_month, i)
    i = i+1
  end
  j = j+d
  if m > 2 && ((y % Int32(400) == z) ||
               (y % Int32(4)   == z &&
                y % Int32(100) != z))
    j = j+1
  end
  return T(j)
end
md2j(y::AbstractString, m::AbstractString, d::AbstractString) = md2j(parse(Int, y), parse(Int, m), parse(Int, d))

@doc """
    t_arr!(B::Array{Int32,1}, t::Int64)

Convert `t` to [year, day of year, hour, minute, second, frac_second], overwriting the first 6 values in `B` with the result.
""" t_arr!
function t_arr!(tbuf::Array{Int32, 1}, t::Int64)
  dt = u2d(t*μs)
  tbuf[1] = Int32(year(dt))
  tbuf[2] = md2j(tbuf[1], Int32(month(dt)), Int32(day(dt)))
  tbuf[3] = Int32(hour(dt))
  tbuf[4] = Int32(minute(dt))
  tbuf[5] = Int32(second(dt))
  tbuf[6] = Int32(millisecond(dt))
  return nothing
end

function unpack_u8(v::UInt8)
  a = signed(div(v,0x10))*Int8(10)
  b = signed(rem(v,0x10))
  return Int64(a+b)
end

function datehex2μs!(a::Array{Int64,1}, datehex::Array{UInt8,1})
  a[1] = 100*unpack_u8(getindex(datehex, 1)) + unpack_u8(getindex(datehex, 2))
  a[2] = md2j(getindex(a,1), unpack_u8(getindex(datehex, 3)), unpack_u8(getindex(datehex, 4))) - 1
  setindex!(a, y2μs(getindex(a,1)), 1)
  a[3] = unpack_u8(getindex(datehex, 5))
  a[4] = unpack_u8(getindex(datehex, 6))
  a[5] = unpack_u8(getindex(datehex, 7))
  a[6] = unpack_u8(getindex(datehex, 8))
  return a[1] + a[2]*86400000000 + a[3]*3600000000 + a[4]*60000000 + a[5]*1000000 + a[6]*10000
end

function y2μs(y::T) where T<:Integer
  y = Int64(y)-1
  return 86400000000 * (y*365 + div(y,4) - div(y,100) + div(y,400)) - 62135596800000000
end

# ts = round(Int64, d2u(DateTime(iv[1], m, d, iv[3], iv[4], iv[5], iv[6]))*sμ
mktime(y::T, j::T, h::T, m::T, s::T, μ::T) where T<:Integer = (y2μs(y) +
          Int64(j-one(T))*86400000000 +
          Int64(h)*3600000000 +
          Int64(m)*60000000 +
          Int64(s)*1000000 +
          Int64(μ))
mktime(t::Array{T,1}) where T<:Integer =(y2μs(t[1]) +
          Int64(t[2]-one(T))*86400000000 +
          Int64(t[3])*3600000000 +
          Int64(t[4])*60000000 +
          Int64(t[5])*1000000 +
          Int64(t[6]))

# =====================================================================
# TimeSpec parsing
"""
    (str0, str1) = parsetimewin(ts1::TimeSpec, ts2::TimeSpec)

Convert times `s` and `t` to strings and sorts s.t. d0 < d1.

See also: `TimeSpec`
"""
function parsetimewin(s::DateTime, t::DateTime)
  if s < t
    return (string(s), string(t))
  else
    return (string(t), string(s))
  end
end
parsetimewin(s::DateTime, t::String) = parsetimewin(s, DateTime(t))
parsetimewin(s::DateTime, t::Real) = parsetimewin(s, u2d(d2u(s)+t))
parsetimewin(s::Real, t::DateTime) = parsetimewin(t, u2d(d2u(t)+s))
parsetimewin(s::String, t::Union{Real,DateTime}) = parsetimewin(DateTime(s), t)
parsetimewin(s::Union{Real,DateTime}, t::String) = parsetimewin(s, DateTime(t))
parsetimewin(s::String, t::String) = parsetimewin(DateTime(s), DateTime(t))
parsetimewin(s::Real, t::Real) = parsetimewin(u2d(60*floor(Int, time()/60) + s), u2d(60*floor(Int, time()/60) + t))

# convert a formatted time string to integer μs from the Unix epoch
function tstr2int(s::String)
  str = split(s, ".", limit=2)
  if length(str) < 2
    μ = 0
  else
    μ = parse(Int64, rpad(str[2], 6, '0'))
  end
  return DateTime(str[1]).instant.periods.value*1000 - dtconst + μ
end

# convert a time in integer μs (measured from the Unix epoch) to a string
function int2tstr(t::Int64)
  dt = unix2datetime(div(t, 1000000))
  v = 1000*getfield(getfield(getfield(dt, :instant), :periods), :value)
  r = string(t - v + dtconst)
  s = string(dt) * "." * lpad(r, 6, '0')
  return s
end

# =====================================================================
# Functions for SeisIO time matrices (:t field)

# Check whether a time matrix has gaps
function is_gapless(t::Array{Int64, 2})
  (length(t) == 4) || return false

  # Definition of a well-formed time matrix with no gaps
  return ((t[1,1] == 1) && (t[2,1] > 1) && (t[2,2] == 0))
end

function mk_t(nx::Integer, ts::Int64)
  t = Array{Int64, 2}(undef, 2, 2)
  setindex!(t, one(Int64), 1)
  setindex!(t, nx, 2)
  setindex!(t, ts, 3)
  setindex!(t, zero(Int64), 4)
  return t
end

"""
    tx = t_expand(t::Array{Int64,2}, fs::Float64)

Expand SeisIO time matrix `t` for data sampled at `fs` Hz. If `x` is a data
vector whose sample times are represented by `t`, then `tx` is a vector of
sample times where `tx[i]` is the sample time of `x[i]`.
"""
function t_expand(t::Array{Int64,2}, fs::Float64)
  fs == 0.0 && return t[:,2]
  t[end,1] == 1 && return [t[1,2]]
  dt = round(Int64, 1.0/(fs*μs))
  tt = dt.*ones(Int64, t[end,1])
  tt[1] -= dt
  for i = 1:size(t,1)
    tt[t[i,1]] += t[i,2]
  end
  cumsum!(tt, tt)
  return tt
end

"""
    t = t_collapse(tx::Array{Int64, 1}, fs::Float64)

Collapse vector of sample times `tx` sampled at `fs` Hz to compact SeisIO time
matrix representation `t`.
"""
function t_collapse(tt::Array{Int64,1}, fs::Float64)
  if fs == 0.0
    t = hcat(collect(1:1:length(tt)), tt)
  else
    dt = round(Int64, 1.0/(fs*μs))
    ts = Array{Int64,1}([dt; diff(tt)::Array{Int64,1}])
    L = length(tt)
    i = findall(ts .!= dt)
    t = Array{Int64,2}([[1 tt[1]];[i ts[i].-dt]])
    if isempty(i) || i[end] != L
      t = vcat(t, hcat(L,0))
    end
  end
  return t
end

function t_bounds(T::Array{Int64,2}, Δ::Int64)
  tmin = typemax(Int64)
  tmax = typemin(Int64)
  isempty(T) && return(tmin, tmax)
  t0 = 0
  t1 = 0
  n = size(T,1)-1
  if T[n+1,2] != 0
    T = vcat(T, [T[n+1,1] 0])
    n += 1
  end
  w0 = -(Δ)
  for i = 1:n
    t0 = T[i,2] + w0 + Δ
    t1 = t0 + Δ*(T[i+1,1]-T[i,1]-1)
    w0 = t1
    if t0 < tmin
      tmin = t0
    end
    if t1 > tmax
      tmax = t1
    end
  end
  t1 += Δ
  if t1 > tmax
    tmax = t1
  end
  return (tmin, tmax)
end

function t_win(T::Array{Int64,2}, Δ::Int64)
  isempty(T) && return(T)
  n = size(T,1)-1
  if T[n+1,2] != 0
    T = vcat(T, [T[n+1,1] 0])
    n += 1
  end
  w0 = -(Δ)
  W = Array{Int64,2}(undef,n,2)
  for i = 1:n
    W[i,1] = T[i,2] + w0 + Δ
    W[i,2] = W[i,1] + Δ*(T[i+1,1]-T[i,1]-1)
    w0 = W[i,2]
  end
  W[n,2] += Δ
  return W
end
t_win(T::Array{Int64,2}, fs::Float64) = t_win(T, round(Int64, sμ/fs))

function w_time(W::Array{Int64,2}, Δ::Int64)
  n = size(W,1)+1
  T = Array{Int64,2}(undef,n,2)
  T[1,1] = Int64(1)
  T[1,2] = W[1,1]
  for i = 2:n-1
    T[i,1] = T[i-1,1] + div(W[i-1,2]-W[i-1,1], Δ) + 1
    T[i,2] = W[i,1] - W[i-1,2] - Δ
  end
  T[n,1] = T[n-1,1] + div(W[n-1,2]-W[n-1,1], Δ)
  T[n,2] = 0
  if T[n,1] == T[n-1,1]
    T = T[1:n-1,:]
  end
  return T
end
w_time(W::Array{Int64,2}, fs::Float64) = w_time(W, round(Int64, sμ/fs))

# Sort based on start of each time window
function sort_segs!(W::Array{Int64, 2})
  (size(W, 1) < 2) && return
  j = sortperm(W[:,1])
  d = diff(j)
  if (maximum(d) > 1) || (minimum(d) < 1)
    W .= W[j,:]
  end
  return nothing
end

function sort_segs(t::Array{Int64, 2}, Δ::Int64)
  is_gapless(t) && return
  W = t_win(t, Δ)
  sort_segs!(W)
  return w_time(W, Δ)
end

"""
    te = endtime(t::Array{Int64,2}, Δ::Int64)

Compute the time of the last sample in `t`, a SeisIO time matrix sampled at
interval `Δ` [μs] or frequency `fs` [Hz]. Output is integer μs measured from
the Unix epoch.
"""
function endtime(t::Array{Int64,2}, Δ::Int64)
  if isempty(t)
    t_end = 0
  elseif is_gapless(t)
    t_end = t[1,2] + (t[2,1]-1)*Δ
  elseif minimum(t) < 0
    t_end = t_bounds(t, Δ)[2]
  else
    L = size(t,1)
    t_end = (t[L,1]-1)*Δ
    if L > 2
      t_end += getindex(sum(t, dims=1),2)
    else
      t_end += t[1,2] + t[2,2]
    end
  end
  return t_end
end
function endtime(t::Array{Int64,2}, fs::Float64)
  if fs == 0.0
    return isempty(t) ? 0 : t[size(t,1), 2]
  else
    return endtime(t, round(Int64, 1.0/(fs*μs)))
  end
end

"""
    ts = starttime(t::Array{Int64,2}, Δ::Int64)

Compute the time of the first sample in SeisIO time matrix `t`, sampled at
interval `Δ` [μs] or frequency `fs` [Hz]. Output is integer μs measured from
the Unix epoch.
"""
function starttime(t::Array{Int64,2}, Δ::Int64)
  return (if isempty(t)
            0
          elseif minimum(t) < 0
            t_bounds(t, Δ)[1]
          else
            t[1,2]
          end)
end
function starttime(t::Array{Int64,2}, fs::Float64)
  if fs == 0.0
    return isempty(t) ? 0 : minimum(t, dims=1)[2]
  else
    return starttime(t, round(Int64, sμ/fs))
  end
end

#=
Change to t_extend 2020-03-06

Bad behavior in old function for nx == 0 && nt > 0:
if nx == 0
  check for gap
    if gap
      increment T[end,2] by δt

Extending a channel with nx=0 would change gap at end but not add samples.
This should never happen.
=#
"""
    t_extend(T::Array{Int64,2}, t_new::Int64, n_new::Int64, Δ::Int64)

Extend SeisIO time matrix *T* sampled at interval *Δ* μs or frequency *fs* Hz. For matrix *Tᵢ*:
* *t_new* is the start time of the next segment in data vector *Xᵢ*
* *n_new* is the expected number of samples in the next segment of *Xᵢ*

`check_for_gap!` acts as a more specfic case of `t_extend` that operates on a
GphysData structure where `n_new` is known and no time gaps are possible
in the new segment.

This function has a mini-API in the time API (https://github.com/jpjones76/SeisIO.jl/blob/master/docs/DevGuides/time.md).

See also: `check_for_gap!`
"""
function t_extend(T::Array{Int64,2}, ts::Integer, nx::Integer, Δ::Int64)
  nt = div(length(T), 2)
  n0 = 0

  # channel has some data already
  if nt > 0
    (nx == 0) && (return nothing)
    n = T[nt, 1]
    t0 = endtime(T, Δ)
    δt = ts - t0 - Δ

    if abs(δt) > div(Δ, 2)
      if nx == 1
        if T[nt, 2] == 0
          # Case 3: normal, no gap from x[n] to x[n+1], nx=1
          T1 = copy(T)
          T1[nt, 1] = n+nx
          T1[nt, 2] = δt
          return T1

        else
          # Case 4: abnormal (gap before x[n]), gap from x[n] to x[n+1], nx=1
          return vcat(T, [n+1 δt])
        end

      elseif T[nt, 2] == 0
        # Case 5: normal, gap from x[n] to x[n+1], nx>1
        T1 = vcat(T, [n+nx 0])
        T1[nt, 1] = n+1
        T1[nt, 2] = δt
        return T1
      else
        # Case 6: abnormal (gap before x[n]), gap from x[n] to x[n+1], nx>1
        return vcat(T, [n+1 δt; n+nx 0])
      end

    elseif T[nt, 2] > 0
      # Case 2: abnormal (gap before x[n]), no gap from x[n] to x[n+1]
      return vcat(T, [n+nx 0])

    else
      # Case 1: normal (no gap before x[n]), no gap from x[n] to x[n+1]
      setindex!(T, n+nx, nt)
      return nothing
    end

  # extend t to end at ts (counterintuitive syntax)
  elseif nx == 0
    return mk_t(nx, ts)[1:1,:]
  # I really don't like how this behaves...is it even used anymore?

  # behavior for a new channel
  else
    return mk_t(nx, ts)
  end
end

# Change 2020-03-07: Account for fs = 0.0
function t_extend(T::Array{Int64,2}, ts::Integer, nx::Integer, fs::Float64)
  T1 = (if fs == 0.0
    nt = div(length(T), 2)
    T2 = zeros(Int64, nx, 2)
    T2[:,1].=nt.+(1:nx)
    T2[1,2] = ts
    vcat(T, T2)
  else
    t_extend(T, ts, nx, round(Int64, sμ/fs))
  end)
  return T1
end

function tx_float(t::Array{Int64,2}, fs::Float64)
  fs == 0.0 && return map(Float64, t[:,2])
  t[end,1] == 1 && return Float64[t[1,2]]
  Nt = t[end,1]
  dt = 1.0/fs
  tt = dt.*ones(Float64, Nt)
  tt[1] -= dt
  for i = 2:size(t,1)
    tt[t[i,1]] += t[i,2]
  end
  cumsum!(tt, tt)
  return tt
end

function x_inds(t::Array{Int64,2})
  nt = size(t, 1)-1
  inds = zeros(Int64, nt, 2)
  for i in 1:nt
    inds[i,1] = t[i,1]
    inds[i,2] = t[i+1,1] - (i == nt ? 0 : 1)
  end
  if t[nt+1,2] != 0
    inds[nt,2] -= 1
    inds = vcat(inds, [t[nt+1,1] t[nt+1,1]])
  end
  return inds
end


text:
using Geotherm
using Test

@testset "Geotherm.jl" begin
    include("Geometry.jl")
    include("Integrate.jl")
    include("Geotherm.jl")
end