text:
using HTTP
using DataFrames
using Dates

"""
    weather(; apikey, info, loc, units=:metric)

Returns weather information for the place identifed with loc keyword.

# Arguments
apikey::String  : API key to access HERE
loc::NamedTuple : a named tuple with lat and long, zipcode, or name of the place
info::String    : takes weather observation, weather forecast, or astronomy forecast as string

For more information visit: https://developer.here.com/documentation/destination-weather/dev_guide/topics/user-guide.html

Example request: https://developer.here.com/documentation/destination-weather/dev_guide/topics/examples.html
"""
function weather(;apikey, info, loc, units=:metric)
    # Request
    baseurl = "https://weather.cc.api.here.com"
    path = "weather/1.0"
    resource = "report.xml"

    if info == "weather observation" product = "observation"
    elseif info == "weather forecast" product = "forecast_7days_simple"
    elseif info == "weather alert" product = "alerts"
    elseif info == "astronomy forecast" product = "forecast_astronomy"
    else product = ""
    end

    if :latitude ∈ keys(loc) && :longitude ∈ keys(loc) place = "latitude=$(loc[:lat])&longitude=$(loc[:long])"
    elseif :zipcode ∈ keys(loc) place = "zipcode=$(loc[:zipcode])"
    elseif :name ∈ keys(loc) place = "name=$(loc[:name])"
    else place = ""
    end

    url = "$baseurl/$path/$resource?apiKey=$apikey&product=$product&$place&metric=$(units==:metric)"
    r = HTTP.request(:GET, url)
    f = parsexml(String(r.body))
    rootnode = root(f)

    # Collecting data
    if info == "weather observation"
        # TODO: Add funciton arguments and date of observation as metadata
        # Create a blank dataframe
        df = DataFrame(sid = Int64[], latitude = Float64[], longitude = Float64[], distance = Float64[])
        for node in iterate(rootnode)
            if nodename(node) == "observation"
                for childnode in eachelement(node) df[!, nodename(childnode)] = [] end  # TODO: Add types
                break
            end
        end
        # Fill the dataframe
        for node in iterate(rootnode)
            if nodename(node) == "location"
                push!(df[!, :sid], length(df[!, :sid]) + 1)
                push!(df[!, :latitude], parse(Float64, node["latitude"]))
                push!(df[!, :longitude], parse(Float64, node["longitude"]))
                push!(df[!, :distance], parse(Float64, node["distance"]))
            end
            if nodename(node) == "observation"
                for childnode in eachelement(node) push!(df[!, nodename(childnode)], nodecontent(childnode)) end
            end
        end
    elseif info == "weather forecast"
        # TODO: Add function argument and station information as metadata
        # Create a blank dataframe
        df = DataFrame(id = Int64[], datetime = String[], day = String[])
        for node in iterate(rootnode)
            if nodename(node) == "forecast"
                for childnode in eachelement(node) df[!, nodename(childnode)] = [] end  # TODO: Add types
                break
            end
        end
        # Fill the dataframe
        for node in iterate(rootnode)
            if nodename(node) == "forecast"
                push!(df[!, :id], length(df[!, :id]) + 1)
                push!(df[!, :datetime], node["utcTime"])
                push!(df[!, :day], node["weekday"])
                for childnode in eachelement(node) push!(df[!, nodename(childnode)], nodecontent(childnode)) end
            end
        end
    elseif info == "weather alert"
        # TODO: Add approporiate metadata
        # Create a dataframe
        seg = Dict("M" => "morning", "A" => "afternoon", "E" => "evening", "N" => "night")
        df = DataFrame(id = Int64[], type = Int64[], description = String[], daytime=String[])
        for node in iterate(rootnode)
            if nodename(parentnode(node)) == nodename(node) == "alerts"
                push!(df[!, :id], length(df[!, :id]) + 1)
                dt = ""
                for childnode in eachelement(node)
                    if nodename(childnode) == "type" push!(df[!, :type], parse(Int64,nodecontent(childnode))) end
                    if nodename(childnode) == "description" push!(df[!, :description], nodecontent(childnode)) end
                    if nodename(childnode) == "timeSegment" dt *= "$(dayname(parse(Int64, childnode["day_of_week"]))) - $(seg[childnode["segment"]]); " end
                end
                push!(df[!, :daytime], dt)
            end
        end
    elseif info == "astronomy forecast"
        # TODO: Add function argument information as metadata
        # Create a blank dataframe
        df = DataFrame(id = Int64[], datetime = String[])
        for node in iterate(rootnode)
            if nodename(parentnode(node)) == nodename(node) == "astronomy"
                for childnode in eachelement(node) df[!, nodename(childnode)] = [] end  # TODO: Add types
                break
            end
        end
        # Fill the dataframe
        for node in iterate(rootnode)
            if nodename(parentnode(node)) == nodename(node) == "astronomy"
                push!(df[!, :id], length(df[!, :id]) + 1)
                push!(df[!, :datetime], node["utcTime"])
                for childnode in eachelement(node) push!(df[!, nodename(childnode)], nodecontent(childnode)) end
            end
        end
    else df = DataFrame()
    end

    return df
end


text:
# This method is called when a SemidiscretizationHyperbolic is constructed.
# It constructs the basic `cache` used throughout the simulation to compute
# the RHS etc.
function create_cache(mesh::UnstructuredQuadMesh, equations,
                      dg::DG, RealT, uEltype)

  polydeg_ = polydeg(dg.basis)
  nvars = nvariables(equations)

  elements = init_elements(RealT, uEltype, mesh, dg.basis.nodes, nvars, polydeg_)

  interfaces = init_interfaces(uEltype, mesh, nvars, polydeg_)

  boundaries = init_boundaries(RealT, uEltype, mesh, elements, nvars, polydeg_)

  cache = (; elements, interfaces, boundaries)

  # perform a check on the sufficient metric identities condition for free-stream preservation
  # and halt computation if it fails
  if !isapprox(max_discrete_metric_identities(dg, cache), 0, atol=1e-12)
    error("metric terms fail free-stream preservation check with maximum error $(max_discrete_metric_identities(dg, cache))")
  end

  return cache
end

function rhs!(du, u, t,
              mesh::UnstructuredQuadMesh, equations,
              initial_condition, boundary_conditions, source_terms,
              dg::DG, cache)
  # Reset du
  @timed timer() "reset ∂u/∂t" du .= zero(eltype(du))

  # Calculate volume integral
  @timed timer() "volume integral" calc_volume_integral!(
    du, u, mesh,
    have_nonconservative_terms(equations), equations,
    dg.volume_integral, dg, cache)

  # Prolong solution to interfaces
  @timed timer() "prolong2interfaces" prolong2interfaces!(
    cache, u, mesh, equations, dg)

  # Calculate interface fluxes
  @timed timer() "interface flux" calc_interface_flux!(
    cache.elements.surface_flux_values, mesh,
    have_nonconservative_terms(equations), equations,
    dg, cache)

  # Prolong solution to boundaries
  @timed timer() "prolong2boundaries" prolong2boundaries!(
    cache, u, mesh, equations, dg)

  # Calculate boundary fluxes
  @timed timer() "boundary flux" calc_boundary_flux!(
    cache, t, boundary_conditions, mesh, equations, dg)

  # Calculate surface integrals
  @timed timer() "surface integral" calc_surface_integral!(
    du, mesh, equations, dg, cache)

  # Apply Jacobian from mapping to reference element
  #  Note! this routine is reused from dg_curved/dg_2d.jl
  @timed timer() "Jacobian" apply_jacobian!(
    du, mesh, equations, dg, cache)

  # Calculate source terms
  @timed timer() "source terms" calc_sources!(
    du, u, t, source_terms, equations, dg, cache)

  return nothing
end

# flux differencing volume integral on curvilinear quadrilateral elements. Averaging of the
# mapping terms, stored in `contravariant_vectors`, is peeled apart from the evaluation of
# the physical fluxes in each Cartesian direction
function calc_volume_integral!(du, u,
                               mesh::Union{CurvedMesh{2}, UnstructuredQuadMesh},
                               nonconservative_terms, equations,
                               volume_integral::VolumeIntegralFluxDifferencing,
                               dg::DGSEM, cache)
  @threaded for element in eachelement(dg, cache)
    split_form_kernel!(du, u, nonconservative_terms, volume_integral.volume_flux, element,
                       mesh, equations, dg, cache)
  end
end

@inline function split_form_kernel!(du::AbstractArray{<:Any,4}, u,
                                    nonconservative_terms::Val{false}, volume_flux, element,
                                    mesh::Union{CurvedMesh{2}, UnstructuredQuadMesh}, equations,
                                    dg::DGSEM, cache, alpha=true)
  @unpack derivative_split = dg.basis
  @unpack contravariant_vectors = cache.elements

  # Calculate volume integral in one element
  for j in eachnode(dg), i in eachnode(dg)
    u_node = get_node_vars(u, equations, dg, i, j, element)

    # compute the fluxes in the x and y directions
    flux1 = flux(u_node, 1, equations)
    flux2 = flux(u_node, 2, equations)

    # first direction: use consistency of the volume flux to make this evaluation cheaper
    # pull the contravariant vector
    Ja11_node, Ja12_node = get_contravariant_vector(1, contravariant_vectors, i, j, element)
    # compute the two contravariant fluxes
    fluxtilde1 = Ja11_node * flux1 + Ja12_node * flux2
    integral_contribution = alpha * derivative_split[i, i] * fluxtilde1
    add_to_node_vars!(du, integral_contribution, equations, dg, i, j, element)

    # second direction: use consistency of the volume flux to make this evaluation cheaper
    # pull the contravariant vector
    Ja21_node, Ja22_node = get_contravariant_vector(2, contravariant_vectors, i, j, element)
    fluxtilde2 = Ja21_node * flux1 + Ja22_node * flux2
    integral_contribution = alpha * derivative_split[j, j] * fluxtilde2
    add_to_node_vars!(du, integral_contribution, equations, dg, i, j, element)

    # use symmetry of the volume flux for the remaining terms in the first direction
    for ii in (i+1):nnodes(dg)
      u_node_ii = get_node_vars(u, equations, dg, ii, j, element)
      flux1 = volume_flux(u_node, u_node_ii, 1, equations)
      flux2 = volume_flux(u_node, u_node_ii, 2, equations)
      # pull the contravariant vectors and compute the average
      Ja11_node_ii, Ja12_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j, element)
      Ja11_avg = 0.5 * (Ja11_node + Ja11_node_ii)
      Ja12_avg = 0.5 * (Ja12_node + Ja12_node_ii)
      # compute the contravariant sharp flux
      fluxtilde1 = Ja11_avg * flux1 + Ja12_avg * flux2
      integral_contribution = alpha * derivative_split[i, ii] * fluxtilde1
      add_to_node_vars!(du, integral_contribution, equations, dg, i,  j, element)
      integral_contribution = alpha * derivative_split[ii, i] * fluxtilde1
      add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, element)
    end

    # use symmetry of the volume flux for the remaining terms in the second direction
    for jj in (j+1):nnodes(dg)
      u_node_jj = get_node_vars(u, equations, dg, i, jj, element)
      flux1 = volume_flux(u_node, u_node_jj, 1, equations)
      flux2 = volume_flux(u_node, u_node_jj, 2, equations)
      # pull the contravariant vectors and compute the average
      Ja21_node_jj, Ja22_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj, element)
      Ja21_avg = 0.5 * (Ja21_node + Ja21_node_jj)
      Ja22_avg = 0.5 * (Ja22_node + Ja22_node_jj)
      # compute the contravariant sharp flux
      fluxtilde2 = Ja21_avg * flux1 + Ja22_avg * flux2
      integral_contribution = alpha * derivative_split[j, jj] * fluxtilde2
      add_to_node_vars!(du, integral_contribution, equations, dg, i, j,  element)
      integral_contribution = alpha * derivative_split[jj, j] * fluxtilde2
      add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, element)
    end
  end
end

# prolong the solution into the convenience array in the interior interface container
# Note! this routine is for quadrilateral elements with "right-handed" orientation
function prolong2interfaces!(cache, u,
                             mesh::UnstructuredQuadMesh,
                             equations, dg::DG)
  @unpack interfaces = cache

  @threaded for interface in eachinterface(dg, cache)
    primary_element   = interfaces.element_ids[1, interface]
    secondary_element = interfaces.element_ids[2, interface]

    primary_side   = interfaces.element_side_ids[1, interface]
    secondary_side = interfaces.element_side_ids[2, interface]

    if primary_side == 1
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[1, v, i, interface] = u[v, i, 1, primary_element]
      end
    elseif primary_side == 2
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[1, v, i, interface] = u[v, nnodes(dg), i, primary_element]
      end
    elseif primary_side == 3
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[1, v, i, interface] = u[v, i, nnodes(dg), primary_element]
      end
    else # primary_side == 4
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[1, v, i, interface] = u[v, 1, i, primary_element]
      end
    end

    if secondary_side == 1
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[2, v, i, interface] = u[v, i, 1, secondary_element]
      end
    elseif secondary_side == 2
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[2, v, i, interface] = u[v, nnodes(dg), i, secondary_element]
      end
    elseif secondary_side == 3
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[2, v, i, interface] = u[v, i, nnodes(dg), secondary_element]
      end
    else # secondary_side == 4
      for i in eachnode(dg), v in eachvariable(equations)
        interfaces.u[2, v, i, interface] = u[v, 1, i, secondary_element]
      end
    end
  end

  return nothing
end

# compute the numerical flux interface coupling between two elements on an unstructured
# quadrilateral mesh
function calc_interface_flux!(surface_flux_values,
                              mesh::UnstructuredQuadMesh,
                              nonconservative_terms::Val{false}, equations, dg::DG, cache)
  @unpack surface_flux = dg
  @unpack u, start_index, index_increment, element_ids, element_side_ids = cache.interfaces
  @unpack normal_directions = cache.elements

  @threaded for interface in eachinterface(dg, cache)
    # Get neighboring elements
    primary_element   = element_ids[1, interface]
    secondary_element = element_ids[2, interface]

    # Get the local side id on which to compute the flux
    primary_side   = element_side_ids[1, interface]
    secondary_side = element_side_ids[2, interface]

    # initial index for the coordinate system on the secondary element
    secondary_index = start_index[interface]

    # loop through the primary element coordinate system and compute the interface coupling
    for primary_index in eachnode(dg)
      # pull the primary and secondary states from the boundary u values
      u_ll = get_one_sided_surface_node_vars(u, equations, dg, 1, primary_index, interface)
      u_rr = get_one_sided_surface_node_vars(u, equations, dg, 2, secondary_index, interface)

      # pull the outward pointing (normal) directional vector
      #   Note! this assumes a conforming approximation, more must be done in terms of the normals
      #         for hanging nodes and other non-conforming approximation spaces
      outward_direction = get_surface_normal(normal_directions, primary_index, primary_side,
                                             primary_element)

      # Call pointwise numerical flux with rotation. Direction is normalized inside this function
      flux = surface_flux(u_ll, u_rr, outward_direction, equations)

      # Copy flux back to primary/secondary element storage
      # Note the sign change for the normal flux in the secondary element!
      for v in eachvariable(equations)
        surface_flux_values[v, primary_index  , primary_side  , primary_element  ] =  flux[v]
        surface_flux_values[v, secondary_index, secondary_side, secondary_element] = -flux[v]
      end

      # increment the index of the coordinate system in the secondary element
      secondary_index += index_increment[interface]
    end
  end

  return nothing
end

# move the approximate solution onto physical boundaries within a "right-handed" element
function prolong2boundaries!(cache, u,
                             mesh::UnstructuredQuadMesh,
                             equations, dg::DG)
  @unpack boundaries = cache

  @threaded for boundary in eachboundary(dg, cache)
    element = boundaries.element_id[boundary]
    side    = boundaries.element_side_id[boundary]

    if side == 1
      for l in eachnode(dg), v in eachvariable(equations)
        boundaries.u[v, l, boundary] = u[v, l, 1, element]
      end
    elseif side == 2
      for l in eachnode(dg), v in eachvariable(equations)
        boundaries.u[v, l, boundary] = u[v, nnodes(dg), l, element]
      end
    elseif side == 3
      for l in eachnode(dg), v in eachvariable(equations)
        boundaries.u[v, l, boundary] = u[v, l, nnodes(dg), element]
      end
    else # side == 4
      for l in eachnode(dg), v in eachvariable(equations)
        boundaries.u[v, l, boundary] = u[v, 1, l, element]
      end
    end
  end

  return nothing
end

# TODO: Taal dimension agnostic
function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic,
                             mesh::UnstructuredQuadMesh, equations, dg::DG)
  @assert isempty(eachboundary(dg, cache))
end

function calc_boundary_flux!(cache, t, boundary_conditions,
                             mesh::UnstructuredQuadMesh, equations, dg::DG)
  @unpack boundary_condition_types, boundary_indices = boundary_conditions

  calc_boundary_flux_by_type!(cache, t, boundary_condition_types, boundary_indices,
                              mesh, equations, dg)
  return nothing
end

# Iterate over tuples of boundary condition types and associated indices
# in a type-stable way using "lispy tuple programming".
function calc_boundary_flux_by_type!(cache, t, BCs::NTuple{N,Any},
                                     BC_indices::NTuple{N,Vector{Int}},
                                     mesh::UnstructuredQuadMesh, equations, dg::DG) where {N}
  # Extract the boundary condition type and index vector
  boundary_condition = first(BCs)
  boundary_condition_indices = first(BC_indices)
  # Extract the remaining types and indices to be processed later
  remaining_boundary_conditions = Base.tail(BCs)
  remaining_boundary_condition_indices = Base.tail(BC_indices)

  # process the first boundary condition type
  calc_boundary_flux!(cache, t, boundary_condition, boundary_condition_indices, mesh, equations, dg)

  # recursively call this method with the unprocessed boundary types
  calc_boundary_flux_by_type!(cache, t, remaining_boundary_conditions,
                              remaining_boundary_condition_indices, mesh, equations, dg)

  return nothing
end

# terminate the type-stable iteration over tuples
function calc_boundary_flux_by_type!(cache, t, BCs::Tuple{}, BC_indices::Tuple{},
                                     mesh::UnstructuredQuadMesh, equations, dg::DG)
  nothing
end

function calc_boundary_flux!(cache, t, boundary_condition, boundary_indexing,
                             mesh::UnstructuredQuadMesh, equations, dg::DG)
  @unpack surface_flux_values = cache.elements
  @unpack element_id, element_side_id = cache.boundaries

  @threaded for local_index in eachindex(boundary_indexing)
    # use the local index to get the global boundary index from the pre-sorted list
    boundary = boundary_indexing[local_index]

    # get the element and side IDs on the boundary element
    element = element_id[boundary]
    side    = element_side_id[boundary]

    # calc boundary flux on the current boundary interface
    for node in eachnode(dg)
      calc_boundary_flux!(surface_flux_values, t, boundary_condition, mesh, equations, dg, cache,
                          node, side, element, boundary)
    end
  end
end

# inlined version of the boundary flux calculation along a physical interface where the
# boundary flux values are set according to a particular `boundary_condition` function
@inline function calc_boundary_flux!(surface_flux_values, t, boundary_condition,
                                     mesh::UnstructuredQuadMesh, equations, dg::DG, cache,
                                     node_index, side_index, element_index, boundary_index)
  @unpack normal_directions = cache.elements
  @unpack u, node_coordinates = cache.boundaries
  @unpack surface_flux = dg

  # pull the inner solution state from the boundary u values on the boundary element
  u_inner = get_node_vars(u, equations, dg, node_index, boundary_index)

  # pull the outward pointing (normal) directional vector
  outward_direction = get_surface_normal(normal_directions, node_index, side_index, element_index)

  # get the external solution values from the prescribed external state
  x = get_node_coords(node_coordinates, equations, dg, node_index, boundary_index)

  # Call pointwise numerical flux function in the rotated direction on the boundary
  #    Note! the direction is normalized inside this function
  flux = boundary_condition(u_inner, outward_direction, x, t, surface_flux, equations)

  for v in eachvariable(equations)
    surface_flux_values[v, node_index, side_index, element_index] = flux[v]
  end
end

# Note! The local side numbering for the unstructured quadrilateral element implementation differs
#       from the structured TreeMesh or CurvedMesh local side numbering:
#
#      TreeMesh/CurvedMesh sides   versus   UnstructuredMesh sides
#                  4                                  3
#          -----------------                  -----------------
#          |               |                  |               |
#          | ^ eta         |                  | ^ eta         |
#        1 | |             | 2              4 | |             | 2
#          | |             |                  | |             |
#          | ---> xi       |                  | ---> xi       |
#          -----------------                  -----------------
#                  3                                  1
# Therefore, we require a different surface integral routine here despite their similar structure.
function calc_surface_integral!(du, mesh::UnstructuredQuadMesh,
                                equations, dg::DGSEM, cache)
  @unpack boundary_interpolation = dg.basis
  @unpack surface_flux_values = cache.elements

  @threaded for element in eachelement(dg, cache)
    for l in eachnode(dg), v in eachvariable(equations)
      # surface contribution along local sides 2 and 4 (fixed x and y varies)
      du[v, 1,          l, element] += ( surface_flux_values[v, l, 4, element]
                                          * boundary_interpolation[1, 1] )
      du[v, nnodes(dg), l, element] += ( surface_flux_values[v, l, 2, element]
                                          * boundary_interpolation[nnodes(dg), 2] )
      # surface contribution along local sides 1 and 3 (fixed y and x varies)
      du[v, l, 1,          element] += ( surface_flux_values[v, l, 1, element]
                                          * boundary_interpolation[1, 1] )
      du[v, l, nnodes(dg), element] += ( surface_flux_values[v, l, 3, element]
                                          * boundary_interpolation[nnodes(dg), 2] )
    end
  end

  return nothing
end

# This routine computes the maximum value of the discrete metric identities necessary to ensure
# that the approxmiation will be free-stream preserving (i.e. a constant solution remains constant)
# on a curvilinear mesh.
#   Note! Independent of the equation system and is only a check on the discrete mapping terms.
#         Can be used for a metric identities check on CurvedMesh{2} or UnstructuredQuadMesh
function max_discrete_metric_identities(dg::DGSEM, cache)
  @unpack derivative_matrix = dg.basis
  @unpack contravariant_vectors = cache.elements

  ndims_ = size(contravariant_vectors, 1)

  metric_id_dx = zeros(eltype(contravariant_vectors), nnodes(dg), nnodes(dg))
  metric_id_dy = zeros(eltype(contravariant_vectors), nnodes(dg), nnodes(dg))

  max_metric_ids = zero(dg.basis.nodes[1])

  for i in 1:ndims_, element in eachelement(dg, cache)
    # compute D*Ja_1^i + Ja_2^i*D^T
    @views mul!(metric_id_dx, derivative_matrix, contravariant_vectors[i, 1, :, :, element])
    @views mul!(metric_id_dy, contravariant_vectors[i, 2, :, :, element], derivative_matrix')
    local_max_metric_ids = maximum( abs.(metric_id_dx + metric_id_dy) )

    max_metric_ids = max( max_metric_ids, local_max_metric_ids )
  end

  return max_metric_ids
end

ext: jl


text:
function F_(u::T, v::T, w::T) where T<:Real
	R1 = sqrt(v^2 + w^2)
	R2 = sqrt(u^2 + w^2)
	R3 = sqrt(u^2 + v^2)
	R  = sqrt(u^2 + v^2 + w^2)
	return sign(u*v) *
		(
			abs(u) * log( (abs(v) + R3) * R2 / (abs(u) * (abs(v) + R))) +
			abs(v) * log( (abs(u) + R3) * R1 / (abs(v) * (abs(u) + R))) +
			w * atan(abs(u * v) / w / R)
 		)
end

"""
   gbox(x, y, z, xextent, yextent, zextent)

Compute the vertical component of the gravitational attraction
rectangular box bounded by the axis-parallel planes defined
by `xextent`, `yextent`, and `zextent`.
Observation point is given by the cartesian coordinates `x, y, z`.

Output is given in units of m/s^2.
"""
function gbox(x::T, y::T, z::T, x1::Array{T}, y1::Array{T}, z1::Array{T}) where T<:Real
	Vz = zero(eltype(T))
	for i in 1:2
        for j in 1:2
            for k in 1:2
                f = (-1.0)^(i+j+k)
                Vz += f * F_(x - x1[i], y - y1[j], z - z1[k])
            end
        end
    end
	return Vz * 6.6742e-11
end

max_stars_repo_name: ruboerner/GravMag

max_issues_repo_name: ruboerner/GravMag

max_forks_repo_name: ruboerner/GravMag


text:
module TimeDag

using AbstractTrees
using AssociativeWindowAggregation
using Bijections
using DataStructures
using Dates
using LightGraphs
using LinearAlgebra
using PrettyTables
using Random
using RecipesBase
using StaticArrays
using Statistics
using Tables
using TeaFiles
using TimeZones

# Types.
export Block
# Singletons.
export INTERSECT, LEFT, UNION
# Node creation.
export wrap, wrapb
# Source nodes.
export block_node, constant, empty_node, iterdates, pulse, tea_file
# Other nodes.
export active_count, align, coalign, ema, first_knot, history, lag, right, left
# Evaluation & other utilities.
export evaluate, value_type

include("maybe.jl")

include("core.jl")
include("identity_map.jl")
include("constant.jl")  # Constant nodes are special, so we need to know about them first.

include("output_type.jl")

include("block.jl")
include("graph.jl")
include("evaluation.jl")
include("plotting.jl")

include("alignment.jl")

include("ops/common.jl")  # This defines macros useful for definition of other nodes.

include("ops/align.jl")
include("ops/array.jl")
include("ops/conditional.jl")
include("ops/history.jl")
include("ops/lagging.jl")
include("ops/random.jl")
include("ops/simple.jl")
include("ops/sources.jl")
include("ops/window.jl")

end


text:
export WreathProduct, WreathProductElem

###############################################################################
#
#   WreathProduct / WreathProductElem
#
###############################################################################

@doc doc"""
    WreathProduct(N, P) <: Group
> Implements Wreath product of a group `N` by permutation group $P = S_n$,
> usually written as $N \wr P$.
> The multiplication inside wreath product is defined as
>  >  `(n, σ) * (m, τ) = (n*σ(m), στ)`
> where `σ(m)` denotes the action (from the right) of the permutation group on
> `n-tuples` of elements from `N`

# Arguments:
* `N::Group` : the single factor of the group $N$
* `P::Generic.PermGroup` : full `PermutationGroup`
"""
struct WreathProduct{N, T<:Group, PG<:Generic.PermGroup} <: Group
   N::DirectPowerGroup{N, T}
   P::PG

   function WreathProduct(Gr::T, P::PG) where {T, PG<:Generic.PermGroup}
      N = DirectPowerGroup(Gr, Int(P.n))
      return new{Int(P.n), T, PG}(N, P)
   end
end

struct WreathProductElem{N, T<:GroupElem, P<:Generic.perm} <: GroupElem
   n::DirectPowerGroupElem{N, T}
   p::P

   function WreathProductElem(n::DirectPowerGroupElem{N,T}, p::P,
      check::Bool=true) where {N, T, P<:Generic.perm}
      if check
         N == length(p.d) || throw(DomainError(
            "Can't form WreathProductElem: lengths differ"))
      end
      return new{N, T, P}(n, p)
   end
end

###############################################################################
#
#   Type and parent object methods
#
###############################################################################

elem_type(::Type{WreathProduct{N, T, PG}}) where {N, T, PG} = WreathProductElem{N, elem_type(T), elem_type(PG)}

parent_type(::Type{WreathProductElem{N, T, P}}) where {N, T, P} =
   WreathProduct{N, parent_type(T), parent_type(P)}

parent(g::WreathProductElem) = WreathProduct(parent(g.n[1]), parent(g.p))

###############################################################################
#
#   Parent object call overloads
#
###############################################################################

function (G::WreathProduct{N})(g::WreathProductElem{N}) where {N}
   n = G.N(g.n)
   p = G.P(g.p)
   return WreathProductElem(n, p)
end

@doc doc"""
    (G::WreathProduct)(n::DirectPowerGroupElem, p::Generic.perm)
> Creates an element of wreath product `G` by coercing `n` and `p` to `G.N` and
> `G.P`, respectively.
"""
(G::WreathProduct)(n::DirectPowerGroupElem, p::Generic.perm) = WreathProductElem(n,p)

(G::WreathProduct)() = WreathProductElem(G.N(), G.P(), false)

@doc doc"""
    (G::WreathProduct)(p::Generic.perm)
> Returns the image of permutation `p` in `G` via embedding `p -> (id,p)`.
"""
(G::WreathProduct)(p::Generic.perm) = G(G.N(), p)

@doc doc"""
    (G::WreathProduct)(n::DirectPowerGroupElem)
> Returns the image of `n` in `G` via embedding `n -> (n,())`. This is the
> embedding that makes the sequence `1 -> N -> G -> P -> 1` exact.
"""
(G::WreathProduct)(n::DirectPowerGroupElem) = G(n, G.P())

(G::WreathProduct)(n,p) = G(G.N(n), G.P(p))

###############################################################################
#
#   Basic manipulation
#
###############################################################################

function hash(G::WreathProduct, h::UInt)
   return hash(G.N, hash(G.P, hash(WreathProduct, h)))
end

function hash(g::WreathProductElem, h::UInt)
   return hash(g.n, hash(g.p, hash(WreathProductElem, h)))
end

###############################################################################
#
#   String I/O
#
###############################################################################

function show(io::IO, G::WreathProduct)
   print(io, "Wreath Product of $(G.N.group) by $(G.P)")
end

function show(io::IO, g::WreathProductElem)
   print(io, "($(g.n)≀$(g.p))")
end

###############################################################################
#
#   Comparison
#
###############################################################################

function (==)(G::WreathProduct, H::WreathProduct)
   G.N == H.N || return false
   G.P == H.P || return false
   return true
end

function (==)(g::WreathProductElem, h::WreathProductElem)
   g.n == h.n || return false
   g.p == h.p || return false
   return true
end

###############################################################################
#
#   Group operations
#
###############################################################################

(p::perm)(n::DirectPowerGroupElem) = DirectPowerGroupElem(n.elts[p.d])

@doc doc"""
    *(g::WreathProductElem, h::WreathProductElem)
> Return the wreath product group operation of elements, i.e.
>
> `g*h = (g.n*g.p(h.n), g.p*h.p)`,
>
> where `g.p(h.n)` denotes the action of `g.p::Generic.perm` on
> `h.n::DirectPowerGroupElem` via standard permutation of coordinates.
"""
function *(g::WreathProductElem, h::WreathProductElem)
   return WreathProductElem(g.n*g.p(h.n), g.p*h.p, false)
end

^(g::WreathProductElem, n::Integer) = Base.power_by_squaring(g, n)

@doc doc"""
    inv(g::WreathProductElem)
> Returns the inverse of element of a wreath product, according to the formula
>   `g^-1 = (g.n, g.p)^-1 = (g.p^-1(g.n^-1), g.p^-1)`.
"""
function inv(g::WreathProductElem)
   pinv = inv(g.p)
   return WreathProductElem(pinv(inv(g.n)), pinv, false)
end

###############################################################################
#
#   Misc
#
###############################################################################

matrix_repr(g::WreathProductElem) = Any[matrix_repr(g.p) g.n]

function iterate(G::WreathProduct)
   n, state_N = iterate(G.N)
   p, state_P = iterate(G.P)
   return G(n,p), (state_N, p, state_P)
end

function iterate(G::WreathProduct, state)
   state_N, p, state_P = state
   res = iterate(G.N, state_N)

   if res == nothing
      resP = iterate(G.P, state_P)
      if resP == nothing
         return nothing
      else
         n, state_N = iterate(G.N)
         p, state_P = resP
      end
   else
      n, state_N = res
   end

   return G(n,p), (state_N, p, state_P)
end

eltype(::Type{WreathProduct{N,G,PG}}) where {N,G,PG} = WreathProductElem{N, elem_type(G), elem_type(PG)}

order(G::WreathProduct) = order(G.P)*order(G.N)
length(G::WreathProduct) = order(G)


text:
# An ImmutableDict constructor based on a list of Pair arguments has been introduced in Julia 1.6.0 and backported to Julia 1.5
if VERSION < v"1.5"
    function Base.ImmutableDict(KV::Pair{K,V}, KVs::Pair{K,V}...) where {K,V}
        d = Base.ImmutableDict(KV)
        for p in KVs
            d = Base.ImmutableDict(d, p)
        end
        return d
    end
end

"""
    @convert(<T1>::<T2>,
        <conversions>
    )

Generate `convert` functions both ways between ArchGDAL Enum of typeids (e.g. `ArchGDAL.OGRFieldType`)
and other types or typeids.

ArchGDAL uses Enum types, listing typeids of various data container used in GDAL/OGR object model.
Some of these types are used to implement concrete types in julia through parametric composite types
based on those Enum of typeids (e.g. `Geometry` and `IGeometry` types with `OGRwkbGeometryType`)

Other types or typeids can be:
- GDAL CEnum.Cenum typeids (e.g. `GDAL.OGRFieldType`),
- Base primitive DataType types (e.g. `Bool`),
- other parametric composite types (e.g. `ImageCore.Normed`)

# Arguments
- `(<T1>::<T2>)::Expr`: source and target supertypes, where `T1<:Enum`  and `T2<:CEnum.Cenum || T2::Type{DataType} || T2::UnionAll}``
- `(<stype1>::<stype2>)::Expr`: source and target subtypes or type ids with `stype1::T1` and
    - `stype2::T2 where T2<:CEnum.Cenum` or
    - `stype2::T2 where T2::Type{DataType}` or
    - `stype2<:T2`where T2<:UnionAll
- ...

**Note:** In the case where the mapping is not bijective, the last declared typeid of subtype is used.
Example:
```
@convert(
    OGRFieldType::DataType,
    OFTInteger::Bool,
    OFTInteger::Int16,
    OFTInteger::Int32,
)
```
will generate a `convert` functions giving:
- `Int32` type for `OFTInteger` and not `Ìnt16`
- `OFTInteger` OGRFieldType typeid for both `Int16` and `Int32`

# Usage
### General case:
```
@convert(GDALRWFlag::GDAL.GDALRWFlag,
    GF_Read::GDAL.GF_Read,
    GF_Write::GDAL.GF_Write,
)
```
does the equivalent of
```
const GDALRWFlag_to_GDALRWFlag_map = ImmutableDict(
    GF_Read => GDAL.GF_Read,
    GF_Write => GDAL.GF_Write
)
Base.convert(::Type{GDAL.GDALRWFlag}, ft::GDALRWFlag) =
    GDALRWFlag_to_GDALRWFlag_map[ft]

const GDALRWFlag_to_GDALRWFlag_map = ImmutableDict(
    GDAL.GF_Read => GF_Read,
    GDAL.GF_Write => GF_Write
)
Base.convert(::Type{GDALRWFlag}, ft::GDAL.GDALRWFlag) =
    GDALRWFlag_to_GDALRWFlag_map[ft]
```
### Case where 1st type `<: Enum` and 2nd type `== DataType` or `ìsa UnionAll`:
```
@convert(OGRFieldType::DataType,
    OFTInteger::Bool,
    OFTInteger::Int16,
)
```
does the equivalent of
```
const OGRFieldType_to_DataType_map = ImmutableDict(
    OFTInteger => Bool,
    OFTInteger => Int16,
)
Base.convert(::Type{DataType}, ft::OGRFieldType) =
    OGRFieldType_to_DataType_map[ft]

Base.convert(::Type{OGRFieldType}, ft::Type{Bool}) = OFTInteger
Base.convert(::Type{OGRFieldType}, ft::Type{Int16}) = OFTInteger
```

"""
macro convert(args...)
    @assert length(args) > 1
    @assert args[1].head == :(::)
    (T1, T2) = args[1].args

    # Types and type ids / subtypes checks
    stypes1, stypes2 = zip((eval.(a.args) for a in args[2:end])...)
    @assert(eval(T1) <: Enum && all(isa.(stypes1, eval(T1))))
    @assert(
        ((eval(T2) <: CEnum.Cenum) && all(isa.(stypes2, eval(T2)))) ||
        ((eval(T2) isa Type{DataType}) && all(isa.(stypes2, eval(T2)))) ||
        ((eval(T2) isa UnionAll) && all((<:).(stypes2, eval(T2))))
    )

    # Types other representations
    (T1_string, T2_string) = replace.(string.((T1, T2)), "." => "_")
    (type1, type2) = esc.((T1, T2))

    # Subtypes forward and backward mapping
    fwd_map = [Expr(:call, esc(:Pair), esc.(a.args)...) for a in args[2:end]]
    rev_map =
        [Expr(:call, esc(:Pair), esc.(reverse(a.args))...) for a in args[2:end]]

    #! Convert functions generation
    result_expr = Expr(:block)

    #* Forward conversions from ArchGDAL typeids
    T1_to_T2_map_name = esc(Symbol(T1_string * "_to_" * T2_string * "_map"))
    push!(
        result_expr.args,
        :(const $T1_to_T2_map_name = Base.ImmutableDict([$(fwd_map...)]...)),
    )
    push!(
        result_expr.args,
        :(function Base.convert(::Type{$type2}, ft::$type1)
            return $T1_to_T2_map_name[ft]
        end),
    )

    #* Reverse conversions to ArchGDAL typeids
    # Optimization for conversion from types
    if !(eval(T2) <: CEnum.Cenum)
        for stypes in [Tuple(esc.(reverse(a.args))) for a in args[2:end]]
            push!(
                result_expr.args,
                :(
                    function Base.convert(
                        ::Type{$type1},
                        ft::Type{$(stypes[1])},
                    )
                        return $(stypes[2])
                    end
                ),
            )
        end
        # Conversion from typeids
    else
        T2_to_T1_map_name = esc(Symbol(T2_string * "_to_" * T1_string * "_map"))
        push!(
            result_expr.args,
            :(
                const $T2_to_T1_map_name =
                    Base.ImmutableDict([$(rev_map...)]...)
            ),
        )
        push!(
            result_expr.args,
            :(function Base.convert(::Type{$type1}, ft::$type2)
                return $T2_to_T1_map_name[ft]
            end),
        )
    end

    return result_expr
end

macro gdal(args...)
    @assert length(args) > 0
    @assert args[1].head == :(::)
    fhead = (args[1].args[1], GDAL.libgdal)
    returntype = args[1].args[2]
    argtypes = Expr(:tuple, [esc(a.args[2]) for a in args[2:end]]...)
    args = [esc(a.args[1]) for a in args[2:end]]
    return quote
        ccall($fhead, $returntype, $argtypes, $(args...))
    end
end

macro ogrerr(code, message)
    return quote
        if $(esc(code)) != GDAL.OGRERR_NONE
            error($message)
        end
    end
end

macro cplerr(code, message)
    return quote
        if $(esc(code)) != GDAL.CE_None
            error($message)
        end
    end
end

macro cplwarn(code, message)
    return quote
        if $(esc(code)) != GDAL.CE_None
            @warn $message
        end
    end
end

macro cplprogress(progressfunc)
    @cfunction($(esc(progressfunc)), Cint, (Cdouble, Cstring, Ptr{Cvoid}))
end

# """
# Load a `NULL`-terminated list of strings

# That is it expects a "StringList", in the sense of the CPL functions, as a
# NULL terminated array of strings.
# """
# function unsafe_loadstringlist(pstringlist::Ptr{Cstring})
#     stringlist = Vector{String}()
#     (pstringlist == C_NULL) && return stringlist
#     i = 1
#     item = unsafe_load(pstringlist, i)
#     while item != C_NULL
#         push!(stringlist, unsafe_string(item))
#         i += 1
#         item = unsafe_load(pstringlist, i)
#     end
#     stringlist
# end

"""
    metadatadomainlist(obj)

Fetch list of (non-empty) metadata domains.
"""
metadatadomainlist(obj)::Vector{String} =
    GDAL.gdalgetmetadatadomainlist(obj.ptr)

"""
    metadata(obj; domain::AbstractString = "")

Fetch metadata. Note that relatively few formats return any metadata.
"""
metadata(obj; domain::AbstractString = "")::Vector{String} =
    GDAL.gdalgetmetadata(obj.ptr, domain)

"""
    metadataitem(obj, name::AbstractString, domain::AbstractString)

Fetch single metadata item.

### Parameters
* `name` the name of the metadata item to fetch.
* `domain` (optional) the domain to fetch for.

### Returns
The metadata item on success, or an empty string on failure.
"""
function metadataitem(
    obj,
    name::AbstractString;
    domain::AbstractString = "",
)::String
    item = GDAL.gdalgetmetadataitem(obj.ptr, name, domain)
    # Use `=== nothing` instead of isnothing() for performance.
    # See https://github.com/JuliaLang/julia/pull/36444 for context.
    return item === nothing ? "" : item
end

"""
    setconfigoption(option::AbstractString, value)

Set a configuration option for GDAL/OGR use.

Those options are defined as a (key, value) couple. The value corresponding to a
key can be got later with the `getconfigoption()` method.

### Parameters
* `option`  the key of the option
* `value`   the value of the option, or NULL to clear a setting.

This mechanism is similar to environment variables, but options set with
`setconfigoption()` overrides, for `getconfigoption()` point of view, values
defined in the environment.

If `setconfigoption()` is called several times with the same key, the value
provided during the last call will be used.
"""
function setconfigoption(option::AbstractString, value)::Nothing
    GDAL.cplsetconfigoption(option, value)
    return nothing
end

"""
    clearconfigoption(option::AbstractString)

This function can be used to clear a setting.

Note: it will not unset an existing environment variable; it will
just unset a value previously set by `setconfigoption()`.
"""
function clearconfigoption(option::AbstractString)::Nothing
    setconfigoption(option, C_NULL)
    return nothing
end

"""
    getconfigoption(option::AbstractString, default = C_NULL)

Get the value of a configuration option.

The value is the value of a (key, value) option set with `setconfigoption()`.
If the given option was not defined with `setconfigoption()`, it tries to find
it in environment variables.

### Parameters
* `option`  the key of the option to retrieve
* `default` a default value if the key does not match existing defined options

### Returns
the value associated to the key, or the default value if not found.
"""
function getconfigoption(option::AbstractString, default = C_NULL)::String
    result =
        @gdal(CPLGetConfigOption::Cstring, option::Cstring, default::Cstring)
    return (result == C_NULL) ? "" : unsafe_string(result)
end

"""
    setthreadconfigoption(option::AbstractString, value)

Set a configuration option for GDAL/OGR use.

Those options are defined as a (key, value) couple. The value corresponding to a
key can be got later with the `getconfigoption()` method.

### Parameters
* `option`  the key of the option
* `value`   the value of the option

This function sets the configuration option that only applies in the current
thread, as opposed to `setconfigoption()` which sets an option that applies on
all threads.
"""
function setthreadconfigoption(option::AbstractString, value)::Nothing
    GDAL.cplsetthreadlocalconfigoption(option, value)
    return nothing
end

"""
    clearthreadconfigoption(option::AbstractString)

This function can be used to clear a setting.

Note: it will not unset an existing environment variable; it will
just unset a value previously set by `setthreadconfigoption()`.
"""
function clearthreadconfigoption(option::AbstractString)::Nothing
    setthreadconfigoption(option, C_NULL)
    return nothing
end

"""
    getthreadconfigoption(option::AbstractString, default = C_NULL)

Same as `getconfigoption()` but with settings from `setthreadconfigoption()`.
"""
function getthreadconfigoption(option::AbstractString, default = C_NULL)::String
    result = @gdal(
        CPLGetThreadLocalConfigOption::Cstring,
        option::Cstring,
        default::Cstring
    )
    return (result == C_NULL) ? "" : unsafe_string(result)
end


text:
perp(x) = [-x[2]; x[1]]

function axis_aligned_bounding_box(rep::VRepresentation)
    n = fulldim(rep)
    ps = points(rep)
    l = [minimum(p[i] for p in ps) for i = 1 : n]
    u = [maximum(p[i] for p in ps) for i = 1 : n]
    l, u
end

function axis_aligned_bounding_box(rep::HRepresentation)
    axis_aligned_bounding_box(vrep(polyhedron(rep, CDDLibrary())))
end

function polyhedron_constraints(model::Model, rep::HRepresentation, vars::AbstractVector{Variable})
    for ineq in ineqs(rep)
        @constraint(model, ineq.a ⋅ vars <= ineq.β)
    end
    for eq in eqs(rep)
        @constraint(model, eq.a ⋅ vars == eq.β)
    end
end

function polyhedron_constraints(model::Model, rep::VRepresentation, vars::AbstractVector{Variable})
    hr = SimpleHRepresentation(hrep(polyhedron(rep, CDDLibrary())))
    polyhedron_constraints(model, hr, vars)
end

JuMP.getvalue(A::AxisArray) = AxisArray(JuMP.getvalue(A.data), A.axes)

# TODO: use in test:
# function mccormick_envelope_constraints(model::Model, w::Variable, u::Variable, v::Variable, umin, umax)
#     # convex approximation of w == u * v for u ∈ [umin, umax], v ∈ [0, 1]
#     # 4.21
#     @constraints(model, begin
#         w >= umin * v
#         w >= umax * v + u - umax
#         w <= umin * v + u - umin
#         w <= umax * v
#         umin <= u
#         u <= umax
#     end)
# end

# function mccormick_envelope_constraints(model::Model, w::Variable, u::Variable, v::Variable)
#     # convex approximation of w == u * v for u ∈ [0, 1], v ∈ [0, 1]
#     # 4.20
#     for u ∈ [0, 1], v ∈ [0, 1]
#     @constraints(model, begin
#         w <= v
#         w <= u
#         w >= u + v - 1
#         0 <= u
#         u <= 1
#         0 <= v
#         v <= 1
#         0 <= w
#     end)
# end


text:
module Escher

using Reexport

using WebIO
@reexport using Interact
@reexport using Mux

export webio_serve

end


text:
#!/usr/bin/env julia

using Coverage

run(`julia --code-coverage=all --inline=no test/runtests.jl`)
coverage = process_folder();

mkpath("coverage")
LCOV.writefile("coverage/cov.info", coverage);

run(`rm -rf coverage/src`)
run(`cp -r src coverage`)
run(`genhtml -o coverage coverage/cov.info`)

clean_folder("src")
clean_folder("test")


text:
module ImageInterpLast

using CachedArrays, ImageCore, ImageMetadata
using AxisArrays: AxisArrays, AxisArray

import CachedArrays: update_cache!, cache, AbstractCachedArray
import Base: size, getindex, show

export interp_last, ImgItpLast

include("util.jl")
include("interp_last.jl")

end


text:
function updateSP(M::Model)
    for s in Endogenous(M.variables)
		sid = find(State(M.variables),s.name)
		xid = find(State(M.variables,false),s.name)
		nG 	= length(M.G)
		for j = 1:size(M.ProbWeights,2)
			for i = 1:nG
				@inbounds M.SP[i+(j-1)*nG,sid] = M.X[i,xid]
			end
		end
    end
end

function updateXP(M::Model)
    futureloc = Int[v.val.targXid for v in Future(M.variables)]
    M.XP[:,:] = M.G(M.X[:,futureloc],M.SP)
    for i in 1:length(futureloc)
        clamp!(M.XP[:,i],Future(M.variables)[futureloc[i]].bounds[1],Future(M.variables)[futureloc[i]].bounds[2])
    end
end

function updatePolicy(M::Model,ϕ=0.8)
    M.F(M)
    ns = size(M.S,2)
    pid = (M.temp.ploc[1]:M.temp.ploc[end])-ns
    n = length(pid)
    for i =1:length(M.G)
        M.J(M,i)
        x = M.temp.J\M.temp.error[i,:]
        @fastmath @simd for j = 1:n
            @inbounds M.X[i,j] -= ϕ*x[j]
        end
    end
end

function forcebounds(M)
    ns = size(M.S,2)
    pid = (M.temp.ploc[1]:M.temp.ploc[end])-ns
    n = length(pid)
    for j = 1:n
        lb,ub = M.variables[ns+j].bounds[1],M.variables[ns+j].bounds[2]
        for i = 1:length(M.G)
            if M.X[i,j] < lb
                M.X[i,j] = lb
            elseif M.X[i,j] > ub
                M.X[i,j] = ub
            end
        end
    end
end

function solve(M::Model,n::Int=1000,ϕ::Float64=0.8;disp::Int=10000,crit=1e-6)
    MaxError = maximum(abs(M.temp.error),1)
    SError = sum(abs(M.temp.error).>crit*5,1)/length(M.G)
    M.temp.updateD(M)
    updateSP(M)
    updateXP(M)
    updatePolicy(M)
    forcebounds(M)

    for iter = 2:n
        mod(iter,20)==0 && (MaxError = maximum(abs(M.temp.error),1);SError = sum(abs(M.temp.error).>crit*5,1)/length(M.G))
        updateSP(M)
        updateXP(M)
        updatePolicy(M,ϕ)
        forcebounds(M)
        mod(iter,disp)==0 && println(round(log10(MaxError),2),"  ",round(SError,2))
        iter>min(n,5)  && (maximum(MaxError)<crit || all(SError.<0.05)) && (print(iter);break)
    end
    M.temp.updateD(M)
end


text:
using Test
using JDF
using RDatasets
using DataFrames: DataFrame
using CategoricalArrays

@testset "JDF.jl categorical arrays" begin
    df = DataFrame([collect(1:100) for i = 1:3])
    df[!, :x1] = categorical(df[!, :x1])
    df[!, :x2] = categorical(string.(df[!, :x2]))
    JDF.save(df, "a3cate.jdf")
    df2 = DataFrame(JDF.load("a3cate.jdf", cols = [:x2, :x1]); copycols=false)
    @test size(df2, 2) == 2
    @test size(df2, 1) == 100
    @time df2[!, :x1] isa CategoricalVector{Int}
    @time df2[!, :x2] isa CategoricalVector{String}

    rm("a3cate.jdf", force = true, recursive = true)
end

@testset "Guard against Github #27" begin
    iris = dataset("datasets", "iris")
    JDF.save(iris, "iris.jdf")
    JDF.load("iris.jdf")

    rm("iris.jdf", force = true, recursive = true)
end


text:
import .ArchGDAL: RasterDataset, AbstractRasterBand,
  getgeotransform, width, height, getname, getcolorinterp,
  getband, nraster, getdataset, ArchGDAL
using .ArchGDAL.DiskArrays: GridChunks, DiskArrays, eachchunk
const AG = ArchGDAL

struct GDALBand{T} <: AG.DiskArrays.AbstractDiskArray{T,2}
    filename::String
    band::Int
    size::Tuple{Int,Int}
    attrs::Dict{String,Any}
    cs::GridChunks{2}
end
function GDALBand(b, filename, i)
    s = size(b)
    atts = getbandattributes(b)
    GDALBand{AG.pixeltype(b)}(filename, i, s, atts, eachchunk(b))
end
Base.size(b::GDALBand) = b.size
DiskArrays.eachchunk(b::GDALBand) = b.cs
DiskArrays.haschunks(::GDALBand) = DiskArrays.Chunked()
function DiskArrays.readblock!(b::GDALBand, aout, r::AbstractUnitRange...)
    AG.read(b.filename) do ds
        AG.getband(ds, b.band) do bh
            DiskArrays.readblock!(bh, aout, r...)
        end
    end
end
function DiskArrays.writeblock!(b::GDALBand, ain, r::AbstractUnitRange...)
    AG.read(b.filename, flags=AG.OF_Update) do ds
        AG.getband(ds, b.band) do bh
            DiskArrays.writeblock!(bh, ain, r...)
        end
    end
end

struct GDALDataset
    filename::String
    bandsize::Tuple{Int,Int}
    projection::String
    trans::Vector{Float64}
    bands::OrderedDict{String}
end

function GDALDataset(filename)
    AG.read(filename) do r
        nb = AG.nraster(r)
        allbands = map(1:nb) do iband
            b = AG.getband(r, iband)
            gb = GDALBand(b, filename, iband)
            name = AG.GDAL.gdalgetdescription(b.ptr)
            if isempty(name)
                name = AG.getname(AG.getcolorinterp(b))
            end
            name => gb
        end
        proj = AG.getproj(r)
        trans = AG.getgeotransform(r)
        s = AG._common_size(r)
        allnames = first.(allbands)
        if !allunique(allnames)
            allbands = ["Band$i"=>last(v) for (i,v) in enumerate(allbands)]
        end
        GDALDataset(filename, s[1:end-1], proj, trans, OrderedDict(allbands))
    end
end
Base.haskey(ds::GDALDataset, k) = in(k, ("X", "Y")) || haskey(ds.bands, k)
#Implement Dataset interface
function get_var_handle(ds::GDALDataset, name)
    if name == "X"
        range(ds.trans[1], length = ds.bandsize[1], step = ds.trans[2])
    elseif name == "Y"
        range(ds.trans[4], length = ds.bandsize[2], step = ds.trans[6])
    else
        ds.bands[name]
    end
end

get_varnames(ds::GDALDataset) = collect(keys(ds.bands))

function get_var_dims(ds::GDALDataset, d)
    if d === "X"
        return ("X",)
    elseif d==="Y"
        return ("Y",)
    else
        return ("X", "Y")
    end
end

get_global_attrs(ds::GDALDataset) = Dict("projection"=>ds.projection)

function get_var_attrs(ds::GDALDataset, name)
    if name in ("Y", "X")
        Dict{String,Any}()
    else
        merge(ds.bands[name].attrs,get_global_attrs(ds))
    end
end

const colornames = AG.getname.(AG.GDALColorInterp.(0:16))

islat(s) = startswith(uppercase(s), "LAT")
islon(s) = startswith(uppercase(s), "LON")
isx(s) = uppercase(s) == "X"
isy(s) = uppercase(s) == "Y"

function totransform(x, y)
    xstep = diff(x)
    ystep = diff(y)
    if !all(isapprox(first(xstep)), xstep) || !all(isapprox(first(ystep)), ystep)
        throw(ArgumentError("Grid must have regular spacing"))
    end
    Float64[first(x), first(xstep), 0.0, first(y), 0.0, first(ystep)]
end
totransform(x::AbstractRange, y::AbstractRange) =
    Float64[first(x), step(x), 0.0, first(y), 0.0, step(y)]
getproj(userproj::String, attrs) = AG.importPROJ4(userproj)
getproj(userproj::AG.AbstractSpatialRef, attrs) = userproj
function getproj(::Nothing, attrs)
    if haskey(attrs, "projection_PROJ4")
        return AG.importPROJ4(attrs["projection_PROJ4"])
    elseif haskey(attrs, "projection_WKT")
        return AG.importWKT(attrs["projection_WKT"])
    else
        error(
            "Could not determine output projection from attributes, please specify userproj",
        )
    end
end

function create_dataset(
    ::Type{<:GDALDataset},
    outpath,
    gatts,
    dimnames,
    dimvals,
    dimattrs,
    vartypes,
    varnames,
    vardims,
    varattrs,
    varchunks;
    userproj = nothing,
    kwargs...,
)
    @assert length(dimnames) == 2
    proj, trans = if islon(dimnames[1]) && islat(dimnames[2])
        #Lets set srs to EPSG:4326
        proj = AG.importEPSG(4326)
        trans = totransform(dimvals[1], dimvals[2])
        proj, trans
    elseif isx(dimnames[1]) && isy(dimnames[2])
        #Try to find out srs
        proj = getproj(userproj, merge(gatts,varattrs...))
        trans = totransform(dimvals[1], dimvals[2])
        proj, trans
    else
        error("Did not find x, y or lon, lat dimensions in dataset")
    end
    cs = first(varchunks)
    @assert all(isequal(varchunks[1]), varchunks)

    driver = AG.getdriver(AG.extensiondriver(outpath))

    nbands = length(varnames)
    dtype = promote_type(vartypes...)
    s = (length.(dimvals)...,)
    bands = AG.create(
        outpath;
        driver = driver,
        width = length(dimvals[1]),
        height = length(dimvals[2]),
        nbands = nbands,
        dtype = dtype,
        options = ["BLOCKXSIZE=$(cs[1])", "BLOCKYSIZE=$(cs[2])"]
    ) do ds
        AG.setgeotransform!(ds, trans)
        bands = map(1:length(varnames)) do i
            b = AG.getband(ds, i)
            icol = findfirst(isequal(varnames[i]), colornames)
            if isnothing(icol)
                AG.setcolorinterp!(b, AG.GDAL.GDALColorInterp(0))
            else
                AG.setcolorinterp!(b, AG.GDAL.GDALColorInterp(icol - 1))
            end
            AG.GDAL.gdalsetdescription(b.ptr, varnames[i])
            atts = varattrs[i]
            haskey(atts, "missing_value") && AG.setnodatavalue!(b, atts["missing_value"])
            if haskey(atts, "labels")
                labeldict = atts[labels]
                maxlabel = maximum(keys(labeldict))
                kt = keytype(labeldict)
                labelvec = [haskey(labeldict, et(i)) ? labeldict[et(i)] : "" for i = 0:maxlabel]
                AG.setcategorynames!(b, labelvec)
            end
            haskey(atts, "units") && AG.setunittype!(b, atts["units"])
            haskey(atts, "scale_factor") && AG.setscale!(b, atts["scale_factor"])
            haskey(atts, "add_offset") && AG.setoffset!(b, atts["add_offset"])
            GDALBand{dtype}(outpath, i, s, atts, AG.DiskArrays.GridChunks(s,cs))
        end
    end
    return GDALDataset(outpath, s, AG.toPROJ4(proj), trans, OrderedDict(vn=>b for (vn,b) in zip(varnames, bands)))
end

allow_parallel_write(::Type{<:GDALDataset}) = false
allow_parallel_write(::GDALDataset) = false

allow_missings(::Type{<:GDALDataset}) = false
allow_missings(::GDALDataset) = false

backendlist[:gdal] = GDALDataset
push!(backendregex,r".tif$"=>GDALDataset)
push!(backendregex,r".gtif$"=>GDALDataset)
push!(backendregex,r".tiff$"=>GDALDataset)
push!(backendregex,r".gtiff$"=>GDALDataset)


text:
#=
Get the frequency of occurrence of each contiguous block of unique elements in an array.

Note: This returns 3 (not 2) values for [a,b,a].
=#
function get_distrib_sortedarr(arr::Array)
    vals = []
    freqs = Int[]

    val_curr = arr[1]
    num_curr = 1
    for i in 2:length(arr)
        if arr[i] != val_curr
            push!(vals,val_curr)
            push!(freqs,num_curr)
            val_curr = arr[i]
            num_curr = 1
        else
            num_curr += 1
        end
    end
    push!(vals,val_curr)
    push!(freqs,num_curr)

    return vals, freqs
end

#=
Get the cumulative frequencies.

Note: length(cum) == length(freqs)+1
=#
function get_cum(freqs::Array)
    cum = zeros(Int,length(freqs)+1)
    for i in 2:length(freqs)+1
        cum[i] = cum[i-1] + freqs[i-1]
    end
    return cum
end

function get_cum_sortedarr(arr::Array)
    vals, freqs = get_distrib_sortedarr(arr)
    return vals, get_cum(freqs)
end


text:
#md# 中文需要分词，须先配置python环境和安装结巴分词
#md# ### 安装PyCall
#md# > `ENV["PYTHON"] = "" #使用内嵌的python`
#md# > `using Pkg`
#md# > `Pkg.build("PyCall")`
#md# ### 安装结巴分词
#md# > `using Conda`
#md# > `Conda.pip_interop(true)`
#md# > `Conda.pip("install","jieba")`
#md#
#md# 安装完成后运行以下示例
using WordCloud
using PyCall

@pyimport jieba

TheInternationale = "起来，饥寒交迫的奴隶！\n起来，全世界受苦的人！\n满腔的热血已经沸腾，\n要为真理而斗争！\n旧世界打个落花流水，\n奴隶们起来，起来！\n不要说我们一无所有，\n我们要做天下的主人！\n\n这是最后的斗争，\n团结起来到明天，\n英特纳雄耐尔\n就一定要实现！\n这是最后的斗争，\n团结起来到明天，\n英特纳雄耐尔\n就一定要实现！\n\n从来就没有什么救世主，\n也不靠神仙皇帝！\n要创造人类的幸福，\n全靠我们自己！\n我们要夺回劳动果实，\n让思想冲破牢笼！\n快把那炉火烧得通红，\n趁热打铁才能成功！\n\n是谁创造了人类世界？\n是我们劳动群众！\n一切归劳动者所有，\n哪能容得寄生虫？！\n最可恨那些毒蛇猛兽，\n吃尽了我们的血肉！\n一旦把它们消灭干净，\n鲜红的太阳照遍全球！\n"

jieba.add_word("英特纳雄耐尔")

wc = wordcloud(
    processtext(jieba.lcut(TheInternationale)),
    colors="#DE2910",
#     mask = WordCloud.randommask(400, color="#FFDE00"),
    mask=pkgdir(WordCloud) * "/res/heart_mask.png",
    maskcolor="#FFDE00",
    density=0.65) |> generate!

println("结果保存在 中文.png")
paint(wc, "中文.png")
wc
#eval# runexample(:中文)
#md# ![](中文.png)

max_stars_repo_name: guo-yong-zhi/WordCloud

max_issues_repo_name: guo-yong-zhi/WordCloud

max_forks_repo_name: guo-yong-zhi/WordCloud


text:
module PeriodicGraphs

using LinearAlgebra
using StaticArrays
using Graphs

export PeriodicVertex, PeriodicEdge, PeriodicGraph,
       PeriodicVertex1D, PeriodicEdge1D, PeriodicGraph1D,
       PeriodicVertex2D, PeriodicEdge2D, PeriodicGraph2D,
       PeriodicVertex3D, PeriodicEdge3D, PeriodicGraph3D,
       coordination_sequence, cellgraph, periodiccellgraph,
       offset_representatives!, swap_axes!, find_edges,
       dimensionality, change_dimension

import Base: (==), isless, convert, show, showerror, eltype, iterate, zero,
             length, in, ndims, print, cmp
import Base.Order: Forward, Lt

"""
    PeriodicVertex{N}

Vertex type for an `N`-periodic graph.

A vertex is uniquely determined by the identifier of its representative in the
a fixed initial cell, and the offset of the cell containing the vertex compared
to the the initial cell.
Vertex identifiers start at 1.
"""
struct PeriodicVertex{N}
    v::Int
    ofs::SVector{N,Int}
end
PeriodicVertex{N}(n::Integer) where {N} = PeriodicVertex{N}(n, zero(SVector{N,Int}))
PeriodicVertex(v::Integer, ofs::Union{SVector{N},NTuple{N}}) where {N} = PeriodicVertex{N}(v, ofs)

function show(io::IO, x::PeriodicVertex{N}) where N
    if get(io, :typeinfo, Any) != PeriodicVertex{N}
        print(io, PeriodicVertex{N})
    end
    print(io, '(', x.v, ", (", join(x.ofs, ','))
    N == 1 && print(io, ',')
    print(io, ')', ')')
end
function convert(::Type{PeriodicVertex{N}}, (dst, offset)::Tuple{Integer,Any}) where N
    PeriodicVertex{N}(dst, offset)
end
function convert(::Type{PeriodicVertex}, (dst, offset)::Tuple{Integer,Union{SVector{N},NTuple{N}}}) where N
    PeriodicVertex{N}(dst, offset)
end
function cmp(x::PeriodicVertex{N}, y::PeriodicVertex{N}) where N
    c = cmp(x.v, y.v)
    iszero(c) || return c
    return cmp(x.ofs, y.ofs)
end
isless(x::PeriodicVertex{N}, y::PeriodicVertex{N}) where {N} = cmp(x, y) < 0
==(x::PeriodicVertex{N}, y::PeriodicVertex{M}) where {N,M} = N == M && iszero(cmp(x, y))

const PeriodicVertex1D = PeriodicVertex{1}
const PeriodicVertex2D = PeriodicVertex{2}
const PeriodicVertex3D = PeriodicVertex{3}

ZtoN(x::Signed) = -(x<0) + 2*abs(x)
function hash_position((x1,x2,x3)::SVector{3,<:Integer})
    x1 = ZtoN(x1); x2 = ZtoN(x2); x3 = ZtoN(x3)
    _b = x1 >= x2
    b1 = _b & (x1 >= x3)
    b2 = (!_b) & (x2 >= x3)
    b3 = (!b1) & (!b2)
    return b1*(x2*(x1 + 1) + x3 + x1^3) +
           b2*((x2 + 1)*(x1 + x2 + 1) + x3 + x2^3) +
           b3*((x3 + 1)*(2x3 + 1) + x3*(x1 + x3^2) + x2)
end

function hash_position((x1,x2)::SVector{2,<:Integer})
    x1 = ZtoN(x1); x2 = ZtoN(x2)
    b = x1 >= x2
    return b*(x2 + x1^2) + (!b)*(x1 + x2*(x2 + 1) + 1)
end

function hash_position((x,)::SVector{1,<:Integer})
    return ZtoN(x)
end

function hash_position(::SVector{0,<:Integer})
    return 0
end

"""
    hash_position(x::PeriodicVertex{N}, n::Integer) where {N}

Given a `PeriodicVertex{N}` and the number n of vertex identifiers in a graph, compute a
unique hash for the given vertex.

This hash function is a bijection between the set of all the vertices of the graph and
the set of positive integers. Its value is an integer between `1+n\\*(2d-1)^N`
(or `1` if `d == 0`) and `n*(2d+1)^N`, where `d = max.(abs.(x.ofs))`.
This means that when one unit cell A is further than another B (for the Manhattan
distance), all vertices in A will have a larger hash than all vertices in B.
"""
function hash_position(x::PeriodicVertex, n::Integer)
    return x.v + n*hash_position(x.ofs)
end

"""
    LoopException <: Exception

Error type for constructing an invalid `PeriodicEdge{N}` of the form (u, u, zeros(Int,N)).
Loops are not expected in the algorithms implemented in PeriodicGraphs.jl. If you
still want to construct them, use the `unsafe_edge{N}` constructor instead of
`PeriodicEdge{N}``.
"""
struct LoopException <: Exception
    src::Int
end
showerror(io::IO, e::LoopException) = print(io, "LoopException: a loop from vertex $(e.src) to itself in the same unit cell is a forbidden edges. Maybe the offset is wrong?")
@noinline __throw_loopexception(src) = throw(LoopException(src))

"""
    unsafe_edge{N}

Internal constructor for `PeriodicEdge{N}` that bypasses the loop check.
This function is not part of the official API, use with caution.
"""
struct unsafe_edge{N} end

"""
    PeriodicEdge{N} <: Graphs.SimpleGraphs.AbstractSimpleEdge{Int}

Edge type for an `N`-periodic graph.

An edge is uniquely determined by the vertex identifiers of its source and
destination, and the cell offset between the source vertex and the destination vertex.
"""
struct PeriodicEdge{N} <: Graphs.SimpleGraphs.AbstractSimpleEdge{Int}
    src::Int
    dst::PeriodicVertex{N}
    function (::Type{unsafe_edge{N}})(src, dst) where N
        return new{N}(src, dst)
    end
end
unsafe_edge{N}(src, dst, ofs) where {N} = unsafe_edge{N}(src, PeriodicVertex{N}(dst, ofs))
function PeriodicEdge{N}(src, dst::PeriodicVertex{N}) where N
    src == dst.v && all(iszero.(dst.ofs)) && __throw_loopexception(src)
    return unsafe_edge{N}(src, dst)
end
PeriodicEdge(src, dst::PeriodicVertex{N}) where {N} = PeriodicEdge{N}(src, dst)
function PeriodicEdge{N}(src, dst, offset) where N
    return PeriodicEdge{N}(src, PeriodicVertex{N}(dst, offset))
end
function PeriodicEdge(src, dst, offset::Union{SVector{N,T},NTuple{N,T}}) where {N,T}
    PeriodicEdge{N}(src, dst, offset)
end
function PeriodicEdge{N}((src, dst, offset)::Tuple{Any,Any,Union{SVector{N,T},NTuple{N,T}}}) where {N,T<:Integer}
    PeriodicEdge{N}(src, dst, offset)
end
function PeriodicEdge((src, dst, offset)::Tuple{Any,Any,Union{SVector{N,T},NTuple{N,T}}}) where {N,T<:Integer}
    PeriodicEdge{N}(src, dst, offset)
end

function convert(::Type{PeriodicEdge{N}}, (src, dst, offset)::Tuple{Any,Any,Any}) where {N}
    PeriodicEdge{N}(src, dst, offset)
end
function convert(::Type{PeriodicEdge}, (src, dst, offset)::Tuple{Any,Any,Union{SVector{N,T},NTuple{N,T}}}) where {N,T<:Integer}
    PeriodicEdge{N}(src, dst, offset)
end
function convert(::Union{Type{PeriodicEdge},Type{PeriodicEdge{N}}}, (src, v)::Tuple{Any,PeriodicVertex{N}}) where N
    PeriodicEdge{N}(src, v)
end

function show(io::IO, x::PeriodicEdge{N}) where N
    if get(io, :typeinfo, Any) != PeriodicEdge{N}
        print(io, PeriodicEdge{N})
    end
    print(io, '(', x.src, ", ", x.dst.v, ", (", join(x.dst.ofs, ','))
    N == 1 && print(io, ',')
    print(io, ')', ')')
end

const PeriodicEdge1D = PeriodicEdge{1}
const PeriodicEdge2D = PeriodicEdge{2}
const PeriodicEdge3D = PeriodicEdge{3}

function Graphs.reverse(e::PeriodicEdge{N}) where N
    return unsafe_edge{N}(e.dst.v, e.src, .-e.dst.ofs)
end
Graphs.src(e::PeriodicEdge) = e.src
Graphs.dst(e::PeriodicEdge) = e.dst.v

"""
    ofs(v::PeriodicVertex)
    ofs(e::PeriodicEdge)

The offset of a vertex or an edge.
"""
ofs(v::PeriodicVertex) = v.ofs
ofs(e::PeriodicEdge) = ofs(e.dst)
isindirectedge(e::PeriodicEdge) = e.src > e.dst.v || (e.src == e.dst.v && e.dst.ofs < zero(e.dst.ofs))
function cmp(x::PeriodicEdge{N}, y::PeriodicEdge{N}) where N
    c = cmp(x.src, y.src)
    iszero(c) || return c
    return cmp(x.dst, y.dst)
end
isless(x::PeriodicEdge{N}, y::PeriodicEdge{N}) where {N} = cmp(x, y) < 0
==(x::PeriodicEdge{N}, y::PeriodicEdge{M}) where {N,M} = N == M && iszero(cmp(x, y))

"""
    PeriodicGraph{N} <: AbstractGraph{Int}

Type representing an undirected `N`-periodic graph.
"""
struct PeriodicGraph{N} <: AbstractGraph{Int}
    ne::Base.RefValue{Int}
    nlist::Vector{Vector{PeriodicVertex{N}}} # For each vertex, the sorted list of its neighbors
    directedgestart::Vector{Int}
    width::Base.RefValue{Rational{Int}}
end

const PeriodicGraph1D = PeriodicGraph{1}
const PeriodicGraph2D = PeriodicGraph{2}
const PeriodicGraph3D = PeriodicGraph{3}

function PeriodicGraph{N}(nv::Integer, t, s) where N
    return PeriodicGraph{N}(Ref(nv), t, s, Ref(-1//1))
end

"""
    PeriodicGraph{N}(nv::Integer=0)

Construct a `PeriodicGraph{N}` with `nv` vertices and 0 edge.
"""
function PeriodicGraph{N}(n::Integer = 0) where N
    @assert n >= 0
    return PeriodicGraph{N}(0, [PeriodicVertex{N}[] for _ in 1:n], [1 for _ in 1:n])
end
"""
    PeriodicGraph([nv::Integer, ]edge_list::AbstractVector{PeriodicEdge{N}})
    PeriodicGraph{N}([nv::Integer, ]edge_list::AbstractVector{PeriodicEdge{N}})

Construct a `PeriodicGraph{N}` from a vector of edges.
If `nv` is unspecified, the number of vertices is the highest that is used
in an edge in `edge_list`.

### Implementation Notes
This constructor works the fastest when `edge_list` is sorted by the lexical ordering
and does not contain any duplicates.

## Examples
```julia
julia> el = PeriodicEdge2D[(1, 1, (1,0)), (1, 3, (0,1)), (3, 1, (0,-1))];

julia> g = PeriodicGraph(el)
PeriodicGraph2D(3, PeriodicEdge2D[(1, 1, (1, 0)), (1, 3, (0, 1))])

julia> ne(g)
2
```
"""
function PeriodicGraph{N}(nv::Integer, t::AbstractVector{PeriodicEdge{N}}) where N
    sort!(t); unique!(t)
    ne = length(t)
    g = PeriodicGraph{N}(nv)
    for e in t
        add_edge!(g, e)
    end
    return g
end
PeriodicGraph(nv::Integer, t::AbstractVector{PeriodicEdge{N}}) where {N} = PeriodicGraph{N}(nv, t)
function PeriodicGraph{N}(t::AbstractVector{PeriodicEdge{N}}) where N
    @static if VERSION < v"1.6-"
        isempty(t) && return PeriodicGraph(0, t)
        return PeriodicGraph(maximum(max(e.src, e.dst.v) for e in t), t)
    else
        return PeriodicGraph(maximum(max(e.src, e.dst.v) for e in t; init=0), t)
    end
end
PeriodicGraph(t::AbstractVector{PeriodicEdge{N}}) where {N} = PeriodicGraph{N}(t)

struct KeyString{T,S<:AbstractString}
    x::S
    start::Base.RefValue{Int}
end
function KeyString{T}(x) where T
    KeyString{T,typeof(x)}(x, Ref(firstindex(x)))
end
iterate(::KeyString, ::Nothing) = nothing
function iterate(k::KeyString{T}, (next_char, idx)) where T
    start = idx
    state = nothing
    while isspace(next_char)
        state = iterate(k.x, idx)
        state isa Nothing && return nothing
        start = idx
        next_char, idx = state
    end
    stop = start
    tmp = idx
    while !isspace(next_char)
        state = iterate(k.x, idx)
        state isa Nothing && break
        stop = tmp
        tmp = idx
        next_char, idx = state
    end
    ret = tryparse(T, SubString(k.x, start:stop))
    if ret isa T
        return (ret, state)
    end
    throw(ArgumentError("Input string does not represent a graph"))
end
function iterate(k::KeyString)
    iterate(k, (' ', k.start[]))
end
function Base.popfirst!(k::KeyString)
    next = iterate(k)
    next isa Nothing && throw(ArgumentError("Input string does not represent a graph"))
    char, state = next
    k.start[] = state isa Nothing ? (ncodeunits(k.x) + 1) : last(state)
    return char
end
Base.isempty(k::KeyString) = isempty(SubString(k.x, k.start[]))
Base.IteratorSize(::Type{<:KeyString}) = Base.SizeUnknown()

"""
    PeriodicGraph(key::AbstractString)
    PeriodicGraph{N}(key::AbstractString)

Construct a `PeriodicGraph{N}` from a `key`, which is a string of whitespace-separated
values of the form
`"N src1 dst1 ofs1_1 ofs1_2 ... ofs1_N src2 dst2 ofs2_1 ofs2_2 ... ofs2_N  ...  srcm dstm ofsm_1 ofsm_2 ... ofsm_N"`
where `N` is the number of repeating dimensions of the graph, `m` is the number of edges
and for all `i` between `1` and `m`, the number of edges, the `i`-th edge is described as
- `srci`, the vertex identifier of the source vertex,
- `dsti`, the vertex identifier of the destination vertex and
- `(ofsi_1, ofsi_2, ..., ofsi_N)` the offset of the edge.

This compact representation of a graph can be obtained simply by `print`ing the graph
or with `string`.

## Examples
```julia
julia> PeriodicGraph("2  1 2 0 0  2 1 1 0  1 1 0 -1")
PeriodicGraph2D(2, PeriodicEdge2D[(1, 1, (0,1)), (1, 2, (-1,0)), (1, 2, (0,0))])

julia> PeriodicGraph3D("3  1 1 0 0 1  1 1 0 1 0  1 1 1 0 0")
PeriodicGraph3D(1, PeriodicEdge3D[(1, 1, (0,0,1)), (1, 1, (0,1,0)), (1, 1, (1,0,0))])

julia> string(ans)
"3 1 1 0 0 1 1 1 0 1 0 1 1 1 0 0"

julia> string(PeriodicGraph(ans)) == ans
true
```
"""
function PeriodicGraph{N}(s::AbstractString) where N
    key = KeyString{Int}(s)
    M = popfirst!(key)
    M != N && throw(DimensionMismatch("Cannot construct a $N-dimensional graph from a $M-dimensional key"))
    edgs = PeriodicEdge{N}[]
    while !isempty(key)
        src = popfirst!(key)
        dst = popfirst!(key)
        ofs = SVector{N,Int}([popfirst!(key) for _ in 1:N])
        push!(edgs, PeriodicEdge{N}(src, dst, ofs))
    end
    return PeriodicGraph{N}(edgs)
end
function PeriodicGraph(s::AbstractString)
    N = first(KeyString{Int}(s))
    return PeriodicGraph{N}(s)
end
function PeriodicGraph{N}(g::PeriodicGraph{N}) where N
    PeriodicGraph{N}(g.ne[], [copy(x) for x in g.nlist], copy(g.directedgestart))
end
PeriodicGraph(g::PeriodicGraph{N}) where {N} = PeriodicGraph{N}(g)

function show(io::IO, g::PeriodicGraph{N}) where N
    if get(io, :typeinfo, Any) != PeriodicGraph{N}
        print(io, PeriodicGraph{N})
    end
    print(io, '(', nv(g), ',', ' ', collect(edges(g)), ')')
end
function print(io::IO, g::PeriodicGraph{N}) where N
    print(io, N)
    for e in edges(g)
        print(io, ' ', e.src, ' ', e.dst.v, ' ', join(e.dst.ofs, ' '))
    end
end

"""
    ==(g1::PeriodicGraph, g2::PeriodicGraph)

Test whether two `PeriodicGraph`s have the same representations. Beware that
two `PeriodicGraph`s may be isomorphic but still have different representations.
"""
function ==(g1::PeriodicGraph{N}, g2::PeriodicGraph{M}) where {N,M}
    return N == M && nv(g1) == nv(g2) && edges(g1) == edges(g2)
end

ndims(::PeriodicVertex{N}) where {N} = N
ndims(::PeriodicEdge{N}) where {N} = N
ndims(::PeriodicGraph{N}) where {N} = N

"""
    ndims(::PeriodicVertex)
    ndims(::PeriodicEdge)
    ndims(::PeriodicGraph)

Return the number of repeating dimensions.

## Examples
```julia
julia> ndims(PeriodicVertex{5}(2))
5

julia> ndims(PeriodicEdge3D(1, 3, (0,0,1)))
3

julia> ndims(PeriodicGraph2D(7))
2
```
"""
ndims

Graphs.ne(g::PeriodicGraph) = g.ne[]
Graphs.nv(g::PeriodicGraph) = length(g.nlist)
Graphs.vertices(g::PeriodicGraph) = Base.OneTo(nv(g))
Graphs.edges(g::PeriodicGraph{N}) where {N} = PeriodicEdgeIter{N}(g)
eltype(g::PeriodicGraph{N}) where {N} = PeriodicVertex{N}
Graphs.edgetype(::PeriodicGraph{N}) where {N} = PeriodicEdge{N}
function Graphs.has_edge(g::PeriodicGraph, s, d)
    ((s < 1) | (s > nv(g))) && return false
    #=@inbounds=# begin
        start = g.directedgestart[s]
        lo, hi = s > d ? (1, start-1) : (start, lastindex(g.nlist[s]))
        i = searchsortedfirst(g.nlist[s], d, lo, hi, Lt((x,y)->isless(x.v, y)))
        return i <= length(g.nlist[s]) && g.nlist[s][i].v == d
    end
end

"""
    find_edges(g::PeriodicGraph, s::Int, d::Int)

Return the set of PeriodicVertex `v` of graph `g` such that there is an edge
between a source vertex of identifier `s` and `v`, and the identifier of `v` is `d`.
"""
function find_edges(g::PeriodicGraph{N}, s::Int, d::Int) where N
    ((s < 1) | (s > nv(g))) && return false
    #=@inbounds=# begin
        start = g.directedgestart[s]
        lo, hi = s > d ? (1, start-1) : (start, lastindex(g.nlist[s]))
        rng = searchsorted(g.nlist[s], d, lo, hi, Lt((x,y)->isless(x isa Integer ? x : x.v, y isa Integer ? y : y.v)))
        if s == d
            rng = (2*first(rng) - last(rng) - 1):last(rng)
        end
        return g.nlist[s][rng]
    end
end
function Graphs.has_edge(g::PeriodicGraph, e::PeriodicEdge)
    s, d = e.src, e.dst
    ((s < 1) | (s > nv(g))) && return false
    #=@inbounds=# begin
        start = g.directedgestart[s]
        lo, hi = isindirectedge(e) ? (1, start-1) : (start, lastindex(g.nlist[s]))
        i = searchsortedfirst(g.nlist[s], d, lo, hi, Forward)
        return i <= length(g.nlist[s]) && g.nlist[s][i] == d
    end
end
Graphs.has_edge(g::PeriodicGraph, i, x::PeriodicVertex) = has_edge(g, PeriodicEdge(i, x))
Graphs.outneighbors(g::PeriodicGraph, v::Integer) = g.nlist[v]
Graphs.inneighbors(g::PeriodicGraph, v::Integer) = outneighbors(g, v)
zero(::Type{PeriodicGraph{N}}) where N = PeriodicGraph{N}(0)
Graphs.is_directed(::Type{<:PeriodicGraph}) = false
@static if isdefined(Graphs, :has_contiguous_vertices)
    @inline Graphs.has_contiguous_vertices(::Type{<:PeriodicGraph}) = true
end
Graphs.has_vertex(g::PeriodicGraph, v::Integer) = 1 <= v <= nv(g)
function Graphs.SimpleGraphs.add_vertices!(g::PeriodicGraph{N}, n::Integer) where N
    append!(g.nlist, [PeriodicVertex{N}[] for _ in 1:n])
    append!(g.directedgestart, [1 for _ in 1:n])
    g.width[] = -1
    return n
end
function Graphs.SimpleGraphs.add_vertex!(g::PeriodicGraph{N}) where N
    push!(g.nlist, PeriodicVertex{N}[])
    push!(g.directedgestart, 1)
    g.width[] = -1
    true
end

function _add_edge!(g::PeriodicGraph, e::PeriodicEdge, ::Val{check}) where check
    #=@inbounds=# begin
        s, dst = e.src, e.dst
        neigh = g.nlist[s]
        start = g.directedgestart[s]
        indirectedge = isindirectedge(e)
        lo, hi = indirectedge ? (1, start-1) : (start, lastindex(neigh))
        i = searchsortedfirst(neigh, dst, lo, hi, Forward)
        if check
            i <= length(neigh) && neigh[i] == dst && return false
        end
        g.directedgestart[s] += indirectedge
        insert!(neigh, i, dst)
        return true
    end
end
function Graphs.add_edge!(g::PeriodicGraph, e::PeriodicEdge)
    (src(e) < 1 || src(e) > nv(g) || dst(e) < 1 || dst(e) > nv(g)) && return false
    success = _add_edge!(g, e, Val(true)) && _add_edge!(g, reverse(e), Val(false))
    if success
        g.ne[] += 1
        g.width[] = -1
    end
    return success
end
Graphs.add_edge!(g::PeriodicGraph, i, x::PeriodicVertex) = add_edge!(g, PeriodicEdge(i, x))

function _rem_edge!(g::PeriodicGraph, e::PeriodicEdge, ::Val{check}) where check
    #=@inbounds=# begin
        s, dst = e.src, e.dst
        neigh = g.nlist[s]
        start = g.directedgestart[s]
        indirectedge = isindirectedge(e)
        lo, hi = indirectedge ? (1, start-1) : (start, lastindex(neigh))
        i = searchsortedfirst(neigh, dst, lo, hi, Forward)
        if check
            i <= length(neigh) && neigh[i] == dst || return false
        end
        g.directedgestart[s] -= indirectedge
        deleteat!(neigh, i)
        return true
    end
end
function Graphs.rem_edge!(g::PeriodicGraph, e::PeriodicEdge)
    (src(e) < 1 || src(e) > nv(g) || dst(e) < 1 || dst(e) > nv(g)) && return false
    success = _rem_edge!(g, e, Val(true)) && _rem_edge!(g, reverse(e), Val(false))
    if success
        g.ne[] -= 1
        g.width[] = -1
    end
    return success
end
Graphs.rem_edge!(g::PeriodicGraph, i, x::PeriodicVertex) = rem_edge!(g, PeriodicEdge(i, x))

function Graphs.SimpleGraphs.rem_vertices!(g::PeriodicGraph{N}, t::AbstractVector{<:Integer}, keep_order::Bool=false) where N
    isempty(t) && return collect(1:nv(g))
    sort!(t)
    (first(t) < 1 || last(t) > nv(g)) && throw(ArgumentError("Vertices to be removed must be in the range 1:nv(g)."))

    bt = falses(nv(g))
    bt[t] .= true

    vmap = Int[]
    rev_vmap = zeros(Int, nv(g))

    if keep_order
        append!(vmap, collect(1:nv(g)))
        deleteat!(vmap, bt)
        rev_vmap[vmap] .= 1:length(vmap)
        deleteat!(g.nlist, bt)
        deleteat!(g.directedgestart, bt)
    else
        i_next_vertex_to_del = 1
        next_vertex_to_del = t[i_next_vertex_to_del]
        t_end = length(t)
        g_end = length(g.nlist)
        i = 1
        while i <= g_end
            if i == next_vertex_to_del
                i_next_vertex_to_del += 1
                if i_next_vertex_to_del > t_end
                    next_vertex_to_del = nv(g) + 1
                else
                    next_vertex_to_del = t[i_next_vertex_to_del]
                end
                while t_end >= i_next_vertex_to_del && g_end == t[t_end]
                    t_end -= 1
                    g_end -= 1
                end
                if i < g_end
                    g.nlist[i], g.nlist[g_end] = g.nlist[g_end], g.nlist[i]
                    push!(vmap, g_end)
                    rev_vmap[g_end] = length(vmap)
                end
                g_end -= 1
            else
                push!(vmap, i)
                rev_vmap[i] = length(vmap)
            end
            i += 1
        end
        resize!(g.nlist, g_end)
        resize!(g.directedgestart, g_end)
    end

    counter_edges = 0

    for i in vertices(g)
        neighbors = g.nlist[i]
        remove_edges = falses(length(neighbors))
        startoffset = 1
        for (k, x) in enumerate(neighbors)
            if bt[x.v]
                remove_edges[k] = true
            else
                neigh = PeriodicVertex{N}(rev_vmap[x.v], x.ofs)
                neighbors[k] = neigh
                startoffset += isindirectedge(PeriodicEdge{N}(i, neigh))
            end
        end
        deleteat!(neighbors, remove_edges)
        sort!(neighbors)
        g.directedgestart[i] = startoffset
        counter_edges += startoffset - 1
    end

    g.ne[] = counter_edges
    g.width[] = -1

    return vmap
end

function Graphs.SimpleGraphs.rem_vertex!(g::PeriodicGraph, v::Integer)
    n = nv(g)
    return length(rem_vertices!(g, [v])) == n - 1
end

"""
    PeriodicEdgeIter{N} <: AbstractEdgeIter

Edge iterator type for undirected `N`-periodic graphs.

The iterator only yields edges in the form `(u, v, ofs)` with either `u < v` or
`u == v && ofs > zero(ofs)`.
This is possible because `PeriodicGraph`s are undirected, hence to each edge
`(u, v, ofs)` in the graph corresponds its reverse edge `(v, u, .-ofs)`. The iterator
thus yields each edge of the graph exactly once.
"""
struct PeriodicEdgeIter{N} <: AbstractEdgeIter
    g::PeriodicGraph{N}
end

eltype(::Type{PeriodicEdgeIter{N}}) where {N} = PeriodicEdge{N}
length(iter::PeriodicEdgeIter) = iter.g.ne[]

function iterate(iter::PeriodicEdgeIter{N}, (vertex, neigh)) where N
    nlists = iter.g.nlist
    n = length(nlists)
    #=@inbounds=# while vertex <= n
        if iszero(neigh)
            neigh = iter.g.directedgestart[vertex]
        end
        neighbors = nlists[vertex]
        if neigh > length(neighbors)
            vertex += 1
            neigh = 0
            continue
        end
        return (unsafe_edge{N}(vertex, neighbors[neigh]), (vertex, neigh+1))
    end
    return nothing
end

iterate(iter::PeriodicEdgeIter) = iterate(iter, (1, 0))

function in(edge, iter::PeriodicEdgeIter{N}) where N
    has_edge(iter.g, edge)
end

function cmp(it1::PeriodicEdgeIter{N}, it2::PeriodicEdgeIter{N}) where N
    n = length(it1)
    m = length(it2)
    n == m || return cmp(n,m)
    n == 0 && return 0
    cmpofs = 0
    e1::PeriodicEdge{N}, st1::Tuple{Int,Int} = iterate(it1)
    e2::PeriodicEdge{N}, st2::Tuple{Int,Int} = iterate(it2)
    for _ in 1:n-1
        c = cmp((e1.src, e1.dst.v), (e2.src, e2.dst.v))
        if iszero(c)
            if iszero(cmpofs)
                cmpofs = cmp(e1.dst.ofs, e2.dst.ofs)
            end
            e1, st1 = iterate(it1, st1)
            e2, st2 = iterate(it2, st2)
            continue
        end
        return c
    end
    c = cmp((e1.src, e1.dst.v), (e2.src, e2.dst.v))
    return iszero(c) ? (iszero(cmpofs) ? cmp(e1.dst.ofs, e2.dst.ofs) : cmpofs) : c
end

function isless(it1::PeriodicEdgeIter{N}, it2::PeriodicEdgeIter{N}) where N
    return cmp(it1, it2) < 0
end
function ==(it1::PeriodicEdgeIter{N}, it2::PeriodicEdgeIter{N}) where N
    return iszero(cmp(it1, it2))
end

"""
    vertex_permutation(g::PeriodicGraph, vlist)

Return the `PeriodicGraph` corresponding to `g` with its vertices identifiers
permuted according to `vlist`. `isperm(vlist)` must hold and will not be checked.

See also `induced_subgraph` for the more general case where `vlist` is not a
permutation.
"""
function vertex_permutation(g::PeriodicGraph{N}, vlist) where N
    n = length(vlist)
    newvid = Vector{Int}(undef, n)
    for i in 1:n
        newvid[vlist[i]] = i
    end
    edges = Vector{Vector{PeriodicVertex{N}}}(undef, n)
    startoffsets = [1 for _ in 1:n]
    #=@inbounds=# for i in 1:n
        startoffset = 1
        neighs = copy(g.nlist[vlist[i]])
        edges[i] = neighs
        for j in 1:length(neighs)
            dst = neighs[j]
            neigh = PeriodicVertex{N}(newvid[dst.v], dst.ofs)
            neighs[j] = neigh
            startoffsets[i] += isindirectedge(PeriodicEdge{N}(i, neigh))
        end
        sort!(neighs)
    end
    return PeriodicGraph{N}(g.ne[], edges, startoffsets)
end

function Graphs.induced_subgraph(g::PeriodicGraph{N}, vlist::AbstractVector{U}) where {N, U<:Integer}
    allunique(vlist) || __throw_unique_vlist()
    n = length(vlist)
    n == nv(g) && return (vertex_permutation(g, vlist), vlist)
    newvid = zeros(Int, nv(g))
    for i in 1:n
        newvid[vlist[i]] = i
    end

    ne = 0
    edges = Vector{Vector{PeriodicVertex{N}}}(undef, n)
    startoffsets = [1 for _ in 1:n]
    for i in 1:n
        edges[i] = PeriodicVertex{N}[]
        startne = ne
        for dst in g.nlist[vlist[i]]
            v = newvid[dst.v]
            iszero(v) && continue
            neigh = PeriodicVertex{N}(v, dst.ofs)
            push!(edges[i], neigh)
            ne += isindirectedge(PeriodicEdge{N}(i, neigh))
        end
        startoffsets[i] = 1 + ne - startne
        sort!(edges[i])
    end
    return (PeriodicGraph{N}(ne, edges, startoffsets), vlist)
end
@noinline __throw_unique_vlist() = throw(ArgumentError("Vertices in subgraph list must be unique"))

"""
    offset_representatives!(g::PeriodicGraph, offsets)

In-place modifies graph `g` so that the `i`-th vertex of the new initial cell corresponds
to the `i`-th vertex in cell `offsets[i]` compared to the previous initial cell.

Note that the resulting graph is isomorphic to the initial one, only the representation
has changed.
"""
function offset_representatives!(g::PeriodicGraph{N}, offsets) where N
    n = nv(g)
    length(offsets) == n || __throw_invalid_offsets()
    for i in 1:n
        neighs = g.nlist[i]
        ofsi = offsets[i]
        startoffset = 1
        for j in 1:length(neighs)
            x = neighs[j]
            neigh = PeriodicVertex{N}(x.v, x.ofs .+ ofsi .- offsets[x.v])
            neighs[j] = neigh
            startoffset += isindirectedge(PeriodicEdge{N}(i, neigh))
        end
        sort!(neighs)
        g.directedgestart[i] = startoffset
        g.width[] = -1
    end
    g
end
@noinline __throw_invalid_offsets() = throw(ArgumentError("The size of offsets does not match the number of vertices"))

"""
    swap_axes!(g::PeriodicGraph, t)

In-place modifies graph `g` so that the new initial cell corresponds to the previous
one with its axes swapped according to the permutation `t`.

Note that the resulting graph is isomorphic to the initial one, only the representation
has changed.
"""
function swap_axes!(g::PeriodicGraph{N}, t) where N
    length(t) == N || __throw_invalid_axesswap()
    tt::Vector{Int} = t isa Vector ? t : collect(t)
    for i in vertices(g)
        neighs = g.nlist[i]
        for j in 1:length(neighs)
            x = neighs[j]
            neigh = PeriodicVertex{N}(x.v, x.ofs[tt])
            neighs[j] = neigh
        end
        sort!(neighs)
        # Note: g.width[] is unchanged
    end
    g
end
@noinline __throw_invalid_axesswap() = throw(DimensionMismatch("The number of axes must match the dimension of the graph"))

function Graphs.connected_components(g::PeriodicGraph)
    nvg = nv(g)
    label = zeros(Int, nvg)

    for u in vertices(g)
        label[u] != 0 && continue
        label[u] = u
        Q = Int[]
        push!(Q, u)
        #=@inbounds=# while !isempty(Q)
            src = popfirst!(Q)
            for dst in outneighbors(g, src)
                vertex = dst.v
                if label[vertex] == 0
                    push!(Q, vertex)
                    label[vertex] = u
                end
            end
        end
    end
    c, d = Graphs.components(label)
    return c
end

"""
    graph_width!(g::PeriodicGraph)

Set the `width` internal field of the graph so that the for all `n` ∈ N\\*,
the `n`-th neighbor of any vertex `v` of the initial cell is in a cell
`(i,j,k)` such that `max(abs.((i,k,k))) ≤ 1 + fld((n - 1), width)`.

This function is meant for internal use and will be used whenever the `width` field
is required but unset. If you decide to modify the other internal fields of `g`,
it is probably a good idea to do `g.width[] = -1` so that this function gets
automatically called when needed, unless you are sure the width will not be
affected by your change.

It is never necessary to use this function or to touch the `width` field of the graph
if you only modify the graph using the official API (i.e. if you never directly touch
the fields).
"""
function graph_width!(g::PeriodicGraph{N}) where N
    distances = floyd_warshall_shortest_paths(cellgraph(g)).dists
    extremalpoints = NTuple{N,NTuple{2,Vector{Tuple{Int,Int}}}}([([],[]) for _ in 1:N])
    # a, x ∈ extremalpoints[i][j] where i ∈ ⟦1,N⟧ and j ∈ ⟦1,2⟧ means that
    # vertex x has a neighbor whose offset is a*(-1)^(j-1) along dimension i
    maxa = 1
    for e in edges(g)
        iszero(ofs(e)) && continue
        offset = ofs(e)
        for i in 1:N
            iszero(offset[i]) && continue
            j = signbit(offset[i]) + 1
            a = abs(offset[i])
            if a > maxa
                maxa = a
            end
            push!(extremalpoints[i][j], (a, src(e)))
            push!(extremalpoints[i][3-j], (a, dst(e)))
        end
    end

    width::Rational{Int} = Rational(nv(g)+1)
    for i in 1:N
        if any(isempty.(extremalpoints[i]))
            # TODO check this part
            if width > 1
                width = 1//1
            end
            continue
        end
        for (a1, x1) in extremalpoints[i][1], (a2, x2) in extremalpoints[i][2]
            dist = distances[x1, x2]
            if dist == typemax(Int)
                dist = 1
            end
            d = (dist + 1) // (a1 + a2)
            if d < width
                width = d
            end
        end
    end
    g.width[] = width == nv(g)+1 ? Rational(maxa) : width
end

function Graphs._neighborhood(g::Union{PeriodicGraph{0},PeriodicGraph{1},PeriodicGraph{2},PeriodicGraph{3}}, v::Integer, d::Real, distmx::AbstractMatrix{U}, neighborfn::Function) where U <: Real
    @assert typeof(neighborfn) === typeof(outneighbors)
    N = ndims(g)
    Q = Tuple{PeriodicVertex{N}, U}[]
    d < zero(U) && return Q
    start_vertex = PeriodicVertex{N}(v)
    push!(Q, (start_vertex, zero(U),) )
    n = nv(g)
    width = g.width[]
    if width == -1
        width = graph_width!(g)
    end
    seen_size = n*(2*(1 + fld(d-1, width)) + 1)^N
    seen = falses(seen_size)
    seen[hash_position(start_vertex, n)] = true
    #=@inbounds=# for (src, currdist) in Q
        currdist == d && continue # should be in Q but all its neighbours are too far
        for dst in outneighbors(g, src.v)
            dst = PeriodicVertex{N}(dst.v, dst.ofs .+ src.ofs)
            position = hash_position(dst, n)
            if !seen[position]
                seen[position] = true
                distance = currdist + distmx[src.v, dst.v]
                if distance <= d
                    push!(Q, (dst, distance))
                end
            end
        end
    end
    return Q
end

function Graphs._neighborhood(g::PeriodicGraph{N}, v::Integer, d::Real, distmx::AbstractMatrix{U}, neighborfn::Function) where {N,U <: Real}
    @assert typeof(neighborfn) === typeof(outneighbors)
    Q = Tuple{PeriodicVertex, U}[]
    d < zero(U) && return Q
    start_vertex = PeriodicVertex{N}(v)
    push!(Q, (start_vertex, zero(U),) )
    n = nv(g)
    seen = Set{PeriodicVertex{N}}()
    push!(seen, start_vertex)
    #=@inbounds=# for (src, currdist) in Q
        currdist == d && continue # should be in Q but all its neighbours are too far
        @simd for dst in outneighbors(g, src.v)
            dst = PeriodicVertex(dst.v, dst.ofs .+ src.ofs)
            if dst ∉ seen
                push!(seen, dst)
                distance = currdist + distmx[src.v, dst.v]
                if distance <= d
                    push!(Q, (dst, distance))
                end
            end
        end
    end
    return Q
end

"""
    coordination_sequence(g::PeriodicGraph, v::Integer, dmax)

Compute the list of numbers of `n`-th neighbors of vertex `v` in graph `g`, for
`1 ≤ n ≤ dmax`.
"""
function coordination_sequence(g::PeriodicGraph, v::Integer, dmax)
    Q = Graphs._neighborhood(g, v, dmax, weights(g), outneighbors)
    popfirst!(Q)
    ret = zeros(Int, dmax)
    for (_, d) in Q
        ret[d] += 1
    end
    return ret
end

"""
    cellgraph(g::PeriodicGraph)

Extract a simple graph from `g` by only keeping the edges that are strictly
within the initial cell.

See also `periodiccellgraph` to keep these edges.
"""
function cellgraph(g::PeriodicGraph)
    edgs = [Edge{Int}(src(x), dst(x)) for x in edges(g) if iszero(ofs(x))]
    ret = SimpleGraph(edgs)
    add_vertices!(ret, nv(g) - nv(ret))
    return ret
end

"""
    periodiccellgraph(g::PeriodicGraph)

Extract a simple graph from `g` by removing all indications of offset in the edges.
This means that edges that used to cross the boundaries of the initial cell now
bind the source vertex to the representative of the destination vertex that is in
the initial cell.

Note that these modified edges may turn into loops.

See also `cellgraph` to remove these edges.
"""
function periodiccellgraph(g::PeriodicGraph)
    ret = SimpleGraph([Edge{Int}(i, j.v) for i in vertices(g) for j in outneighbors(g, i)])
    add_vertices!(ret, nv(g) - nv(ret))
    return ret
end

include("dimensionality.jl")
include("precompile.jl")

end