text:
# Container data structure (structure-of-arrays style) for DG elements
# TODO: Taal refactor, remove NVARS, POLYDEG?
mutable struct ElementContainer3D{RealT<:Real, NVARS, POLYDEG} <: AbstractContainer
  inverse_jacobian::Vector{RealT}      # [elements]
  node_coordinates::Array{RealT, 5}    # [orientation, i, j, k, elements]
  surface_flux_values::Array{RealT, 5} # [variables, i, j, direction, elements]
  cell_ids::Vector{Int}                # [elements]
  # internal `resize!`able storage
  _node_coordinates::Vector{RealT}
  _surface_flux_values::Vector{RealT}
end

# Only one-dimensional `Array`s are `resize!`able in Julia.
# Hence, we use `Vector`s as internal storage and `resize!`
# them whenever needed. Then, we reuse the same memory by
# `unsafe_wrap`ping multi-dimensional `Array`s around the
# internal storage.
function Base.resize!(elements::ElementContainer3D{RealT, NVARS, POLYDEG},
                      capacity) where {RealT, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  @unpack _node_coordinates, _surface_flux_values,
          inverse_jacobian, cell_ids = elements

  resize!(inverse_jacobian, capacity)

  resize!(_node_coordinates, 3 * n_nodes * n_nodes * n_nodes * capacity)
  elements.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates),
                                          (3, n_nodes, n_nodes, n_nodes, capacity))

  resize!(_surface_flux_values, NVARS * n_nodes * n_nodes * 2 * 3 * capacity)
  elements.surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values),
                                             (NVARS, n_nodes, n_nodes, 2 * 3, capacity))

  resize!(cell_ids, capacity)

  return nothing
end

function ElementContainer3D{RealT, NVARS, POLYDEG}(capacity::Integer) where {RealT<:Real, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  nan = convert(RealT, NaN)

  # Initialize fields with defaults
  inverse_jacobian = fill(nan, capacity)

  _node_coordinates = fill(nan, 3 * n_nodes * n_nodes * n_nodes * capacity)
  node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates),
                                 (3, n_nodes, n_nodes, n_nodes, capacity))

  _surface_flux_values = fill(nan, NVARS * n_nodes * n_nodes * 2 * 3 * capacity)
  surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values),
                                    (NVARS, n_nodes, n_nodes, 2 * 3, capacity))

  cell_ids = fill(typemin(Int), capacity)

  return ElementContainer3D{RealT, NVARS, POLYDEG}(
    inverse_jacobian, node_coordinates, surface_flux_values, cell_ids,
    _node_coordinates, _surface_flux_values)
end

# Return number of elements
nelements(elements::ElementContainer3D) = length(elements.cell_ids)
# TODO: Taal performance, 1:nelements(elements) vs. Base.OneTo(nelements(elements))
@inline eachelement(elements::ElementContainer3D) = Base.OneTo(nelements(elements))

# Create element container and initialize element data
function init_elements(cell_ids, mesh::TreeMesh3D,
                       equations::AbstractEquations{3, NVARS},
                       basis::LobattoLegendreBasis{T, NNODES}, ::Type{RealT}) where {RealT<:Real, NVARS, T, NNODES}
  # Initialize container
  n_elements = length(cell_ids)
  elements = ElementContainer3D{RealT, NVARS, NNODES-1}(n_elements)

  init_elements!(elements, cell_ids, mesh, basis.nodes)
  return elements
end

function init_elements!(elements, cell_ids, mesh::TreeMesh3D, nodes)
  n_nodes = length(nodes)

  # Store cell ids
  elements.cell_ids .= cell_ids

  # Calculate inverse Jacobian and node coordinates
  for element in eachelement(elements)
    # Get cell id
    cell_id = cell_ids[element]

    # Get cell length
    dx = length_at_cell(mesh.tree, cell_id)

    # Calculate inverse Jacobian as 1/(h/2)
    elements.inverse_jacobian[element] = 2/dx

    # Calculate node coordinates
    for k in 1:n_nodes, j in 1:n_nodes, i in 1:n_nodes
      elements.node_coordinates[1, i, j, k, element] = (
          mesh.tree.coordinates[1, cell_id] + dx/2 * nodes[i])
      elements.node_coordinates[2, i, j, k, element] = (
          mesh.tree.coordinates[2, cell_id] + dx/2 * nodes[j])
      elements.node_coordinates[3, i, j, k, element] = (
          mesh.tree.coordinates[3, cell_id] + dx/2 * nodes[k])
    end
  end

  return elements
end

# Container data structure (structure-of-arrays style) for DG interfaces
# TODO: Taal refactor, remove NVARS, POLYDEG?
mutable struct InterfaceContainer3D{RealT<:Real, NVARS, POLYDEG} <: AbstractContainer
  u::Array{RealT, 5}        # [leftright, variables, i, j, interfaces]
  neighbor_ids::Matrix{Int} # [leftright, interfaces]
  orientations::Vector{Int} # [interfaces]
  # internal `resize!`able storage
  _u::Vector{RealT}
  _neighbor_ids::Vector{Int}
end

# See explanation of Base.resize! for the element container
function Base.resize!(interfaces::InterfaceContainer3D{RealT, NVARS, POLYDEG},
                      capacity) where {RealT, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  @unpack _u, _neighbor_ids, orientations = interfaces

  resize!(_u, 2 * NVARS * n_nodes * n_nodes * capacity)
  interfaces.u = unsafe_wrap(Array, pointer(_u),
                             (2, NVARS, n_nodes, n_nodes, capacity))

  resize!(_neighbor_ids, 2 * capacity)
  interfaces.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids ),
                                        (2, capacity))

  resize!(orientations, capacity)

  return nothing
end

function InterfaceContainer3D{RealT, NVARS, POLYDEG}(capacity::Integer) where {RealT<:Real, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  nan = convert(RealT, NaN)

  # Initialize fields with defaults
  _u = fill(nan, 2 * NVARS * n_nodes * n_nodes * capacity)
  u = unsafe_wrap(Array, pointer(_u),
                  (2, NVARS, n_nodes, n_nodes, capacity))

  _neighbor_ids = fill(typemin(Int), 2 * capacity)
  neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids ),
                             (2, capacity))

  orientations = fill(typemin(Int), capacity)

  return InterfaceContainer3D{RealT, NVARS, POLYDEG}(
    u, neighbor_ids, orientations,
    _u, _neighbor_ids)
end

# Return number of interfaces
ninterfaces(interfaces::InterfaceContainer3D) = length(interfaces.orientations)

# Create interface container and initialize interface data in `elements`.
function init_interfaces(cell_ids, mesh::TreeMesh3D,
                         elements::ElementContainer3D{RealT, NVARS, POLYDEG}) where {RealT<:Real, NVARS, POLYDEG}
  # Initialize container
  n_interfaces = count_required_interfaces(mesh, cell_ids)
  interfaces = InterfaceContainer3D{RealT, NVARS, POLYDEG}(n_interfaces)

  # Connect elements with interfaces
  init_interfaces!(interfaces, elements, mesh)
  return interfaces
end

# Count the number of interfaces that need to be created
function count_required_interfaces(mesh::TreeMesh3D, cell_ids)
  count = 0

  # Iterate over all cells
  for cell_id in cell_ids
    for direction in eachdirection(mesh.tree)
      # Only count interfaces in positive direction to avoid double counting
      if direction % 2 == 1
        continue
      end

      # If no neighbor exists, current cell is small or at boundary and thus we need a mortar
      if !has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # Skip if neighbor has children
      neighbor_id = mesh.tree.neighbor_ids[direction, cell_id]
      if has_children(mesh.tree, neighbor_id)
        continue
      end

      count += 1
    end
  end

  return count
end

# Initialize connectivity between elements and interfaces
function init_interfaces!(interfaces, elements, mesh::TreeMesh3D)
  # Construct cell -> element mapping for easier algorithm implementation
  tree = mesh.tree
  c2e = zeros(Int, length(tree))
  for element in eachelement(elements)
    c2e[elements.cell_ids[element]] = element
  end

  # Reset interface count
  count = 0

  # Iterate over all elements to find neighbors and to connect via interfaces
  for element in eachelement(elements)
    # Get cell id
    cell_id = elements.cell_ids[element]

    # Loop over directions
    for direction in eachdirection(mesh.tree)
      # Only create interfaces in positive direction
      if direction % 2 == 1
        continue
      end

      # If no neighbor exists, current cell is small and thus we need a mortar
      if !has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # Skip if neighbor has children
      neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id]
      if has_children(mesh.tree, neighbor_cell_id)
        continue
      end

      # Create interface between elements (1 -> "left" of interface, 2 -> "right" of interface)
      count += 1
      interfaces.neighbor_ids[2, count] = c2e[neighbor_cell_id]
      interfaces.neighbor_ids[1, count] = element

      # Set orientation (x -> 1, y -> 2, z -> 3)
      if direction in (1, 2)
        interfaces.orientations[count] = 1
      elseif direction in (3, 4)
        interfaces.orientations[count] = 2
      else
        interfaces.orientations[count] = 3
      end
    end
  end

  @assert count == ninterfaces(interfaces) ("Actual interface count ($count) does not match " *
                                            "expectations $(ninterfaces(interfaces))")
end

# Container data structure (structure-of-arrays style) for DG boundaries
# TODO: Taal refactor, remove NVARS, POLYDEG?
mutable struct BoundaryContainer3D{RealT<:Real, NVARS, POLYDEG} <: AbstractContainer
  u::Array{RealT, 5}                # [leftright, variables, i, j, boundaries]
  neighbor_ids::Vector{Int}         # [boundaries]
  orientations::Vector{Int}         # [boundaries]
  neighbor_sides::Vector{Int}       # [boundaries]
  node_coordinates::Array{RealT, 4} # [orientation, i, j, elements]
  n_boundaries_per_direction::SVector{6, Int} # [direction]
  # internal `resize!`able storage
  _u::Vector{RealT}
  _node_coordinates::Vector{RealT}
end

# See explanation of Base.resize! for the element container
function Base.resize!(boundaries::BoundaryContainer3D{RealT, NVARS, POLYDEG},
                      capacity) where {RealT, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  @unpack _u, _node_coordinates,
          neighbor_ids, orientations, neighbor_sides = boundaries

  resize!(_u, 2 * NVARS * n_nodes * n_nodes * capacity)
  boundaries.u = unsafe_wrap(Array, pointer(_u),
                             (2, NVARS, n_nodes, n_nodes, capacity))

  resize!(_node_coordinates, 3 * n_nodes * n_nodes * capacity)
  boundaries.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates ),
                                            (3, n_nodes, n_nodes, capacity))

  resize!(neighbor_ids, capacity)

  resize!(orientations, capacity)

  resize!(neighbor_sides, capacity)

  return nothing
end

function BoundaryContainer3D{RealT, NVARS, POLYDEG}(capacity::Integer) where {RealT<:Real, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  nan = convert(RealT, NaN)

  # Initialize fields with defaults
  _u = fill(nan, 2 * NVARS * n_nodes * n_nodes * capacity)
  u = unsafe_wrap(Array, pointer(_u),
                  (2, NVARS, n_nodes, n_nodes, capacity))

  neighbor_ids = fill(typemin(Int), capacity)

  orientations = fill(typemin(Int), capacity)

  neighbor_sides = fill(typemin(Int), capacity)

  _node_coordinates = fill(nan, 3 * n_nodes * n_nodes * capacity)
  node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates),
                                 (3, n_nodes, n_nodes, capacity))

  n_boundaries_per_direction = SVector(0, 0, 0, 0, 0, 0)

  return BoundaryContainer3D{RealT, NVARS, POLYDEG}(
    u, neighbor_ids, orientations, neighbor_sides,
    node_coordinates, n_boundaries_per_direction,
    _u, _node_coordinates)
end

# Return number of boundaries
nboundaries(boundaries::BoundaryContainer3D) = length(boundaries.orientations)

# Create boundaries container and initialize boundary data in `elements`.
function init_boundaries(cell_ids, mesh::TreeMesh3D,
                         elements::ElementContainer3D{RealT, NVARS, POLYDEG}) where {RealT<:Real, NVARS, POLYDEG}
  # Initialize container
  n_boundaries = count_required_boundaries(mesh, cell_ids)
  boundaries = BoundaryContainer3D{RealT, NVARS, POLYDEG}(n_boundaries)

  # Connect elements with boundaries
  init_boundaries!(boundaries, elements, mesh)
  return boundaries
end

# Count the number of boundaries that need to be created
function count_required_boundaries(mesh::TreeMesh3D, cell_ids)
  count = 0

  # Iterate over all cells
  for cell_id in cell_ids
    for direction in eachdirection(mesh.tree)
      # If neighbor exists, current cell is not at a boundary
      if has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # If coarse neighbor exists, current cell is not at a boundary
      if has_coarse_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # No neighbor exists in this direction -> must be a boundary
      count += 1
    end
  end

  return count
end

# Initialize connectivity between elements and boundaries
function init_boundaries!(boundaries, elements, mesh::TreeMesh3D)
  # Reset boundaries count
  count = 0

  # Initialize boundary counts
  counts_per_direction = MVector(0, 0, 0, 0, 0, 0)

  # OBS! Iterate over directions first, then over elements, and count boundaries in each direction
  # Rationale: This way the boundaries are internally sorted by the directions -x, +x, -y etc.,
  #            obviating the need to store the boundary condition to be applied explicitly.
  # Loop over directions
  for direction in eachdirection(mesh.tree)
    # Iterate over all elements to find missing neighbors and to connect to boundaries
    for element in eachelement(elements)
      # Get cell id
      cell_id = elements.cell_ids[element]

      # If neighbor exists, current cell is not at a boundary
      if has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # If coarse neighbor exists, current cell is not at a boundary
      if has_coarse_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # Create boundary
      count += 1
      counts_per_direction[direction] += 1

      # Set neighbor element id
      boundaries.neighbor_ids[count] = element

      # Set neighbor side, which denotes the direction (1 -> negative, 2 -> positive) of the element
      if direction in (2, 4, 6)
        boundaries.neighbor_sides[count] = 1
      else
        boundaries.neighbor_sides[count] = 2
      end

      # Set orientation (x -> 1, y -> 2)
      if direction in (1, 2)
        boundaries.orientations[count] = 1
      elseif direction in (3, 4)
        boundaries.orientations[count] = 2
      else
        boundaries.orientations[count] = 3
      end

      # Store node coordinates
      enc = elements.node_coordinates
      if direction == 1 # -x direction
        boundaries.node_coordinates[:, :, :, count] .= enc[:, 1,   :,    :,   element]
      elseif direction == 2 # +x direction
        boundaries.node_coordinates[:, :, :, count] .= enc[:, end, :,    :,   element]
      elseif direction == 3 # -y direction
        boundaries.node_coordinates[:, :, :, count] .= enc[:, :,   1,    :,   element]
      elseif direction == 4 # +y direction
        boundaries.node_coordinates[:, :, :, count] .= enc[:, :,   end,  :,   element]
      elseif direction == 5 # -z direction
        boundaries.node_coordinates[:, :, :, count] .= enc[:, :,   :,    1,   element]
      elseif direction == 6 # +z direction
        boundaries.node_coordinates[:, :, :, count] .= enc[:, :,   :,    end, element]
      else
        error("should not happen")
      end
    end
  end

  @assert count == nboundaries(boundaries) ("Actual boundaries count ($count) does not match " *
                                            "expectations $(nboundaries(boundaries))")
  @assert sum(counts_per_direction) == count

  boundaries.n_boundaries_per_direction = SVector(counts_per_direction)

  return SVector(counts_per_direction)
end

# Container data structure (structure-of-arrays style) for DG L2 mortars
# Positions/directions for orientations = 1, large_sides = 2:
# mortar is orthogonal to x-axis, large side is in positive coordinate direction wrt mortar
#   /----------------------------\  /----------------------------\
#   |             |              |  |                            |
#   | upper, left | upper, right |  |                            |
#   |      3      |      4       |  |                            |
#   |             |              |  |           large            |
#   |-------------|--------------|  |             5              |
# z |             |              |  |                            |
#   | lower, left | lower, right |  |                            |
# ^ |      1      |      2       |  |                            |
# | |             |              |  |                            |
# | \----------------------------/  \----------------------------/
# |
# ⋅----> y
# Left and right are always wrt to a coordinate direction:
# * left is always the negative direction
# * right is always the positive direction
#
# Left and right are used *both* for the numbering of the mortar faces *and* for the position of the
# elements with respect to the axis orthogonal to the mortar.
mutable struct L2MortarContainer3D{RealT<:Real, NVARS, POLYDEG} <: AbstractContainer
  u_upper_left ::Array{RealT, 5} # [leftright, variables, i, j, mortars]
  u_upper_right::Array{RealT, 5} # [leftright, variables, i, j, mortars]
  u_lower_left ::Array{RealT, 5} # [leftright, variables, i, j, mortars]
  u_lower_right::Array{RealT, 5} # [leftright, variables, i, j, mortars]
  neighbor_ids::Matrix{Int}      # [position, mortars]
  # Large sides: left -> 1, right -> 2
  large_sides::Vector{Int}   # [mortars]
  orientations::Vector{Int}  # [mortars]
  # internal `resize!`able storage
  _u_upper_left::Vector{RealT}
  _u_upper_right::Vector{RealT}
  _u_lower_left::Vector{RealT}
  _u_lower_right::Vector{RealT}
  _neighbor_ids::Vector{Int}
end

# See explanation of Base.resize! for the element container
function Base.resize!(mortars::L2MortarContainer3D{RealT, NVARS, POLYDEG},
                      capacity) where {RealT, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  @unpack _u_upper_left, _u_upper_right, _u_lower_left, _u_lower_right,
          _neighbor_ids, large_sides, orientations = mortars

  resize!(_u_upper_left, 2 * NVARS * n_nodes * n_nodes * capacity)
  mortars.u_upper_left = unsafe_wrap(Array, pointer(_u_upper_left),
                                     (2, NVARS, n_nodes, n_nodes, capacity))

  resize!(_u_upper_right, 2 * NVARS * n_nodes * n_nodes * capacity)
  mortars.u_upper_right = unsafe_wrap(Array, pointer(_u_upper_right),
                                      (2, NVARS, n_nodes, n_nodes, capacity))

  resize!(_u_lower_left, 2 * NVARS * n_nodes * n_nodes * capacity)
  mortars.u_lower_left = unsafe_wrap(Array, pointer(_u_lower_left),
                                     (2, NVARS, n_nodes, n_nodes, capacity))

  resize!(_u_lower_right, 2 * NVARS * n_nodes * n_nodes * capacity)
  mortars.u_lower_right = unsafe_wrap(Array, pointer(_u_lower_right),
                                      (2, NVARS, n_nodes, n_nodes, capacity))

  resize!(_neighbor_ids, 5 * capacity)
  mortars.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids ),
                                        (5, capacity))

  resize!(large_sides, capacity)

  resize!(orientations, capacity)

  return nothing
end

function L2MortarContainer3D{RealT, NVARS, POLYDEG}(capacity::Integer) where {RealT<:Real, NVARS, POLYDEG}
  n_nodes = POLYDEG + 1
  nan = convert(RealT, NaN)

  # Initialize fields with defaults
  _u_upper_left = fill(nan, 2 * NVARS * n_nodes * n_nodes * capacity)
  u_upper_left = unsafe_wrap(Array, pointer(_u_upper_left),
                             (2, NVARS, n_nodes, n_nodes, capacity))

  _u_upper_right = fill(nan, 2 * NVARS * n_nodes * n_nodes * capacity)
  u_upper_right = unsafe_wrap(Array, pointer(_u_upper_right),
                              (2, NVARS, n_nodes, n_nodes, capacity))

  _u_lower_left = fill(nan, 2 * NVARS * n_nodes * n_nodes * capacity)
  u_lower_left = unsafe_wrap(Array, pointer(_u_lower_left),
                             (2, NVARS, n_nodes, n_nodes, capacity))

  _u_lower_right = fill(nan, 2 * NVARS * n_nodes * n_nodes * capacity)
  u_lower_right = unsafe_wrap(Array, pointer(_u_lower_right),
                              (2, NVARS, n_nodes, n_nodes, capacity))

  _neighbor_ids = fill(typemin(Int), 5 * capacity)
  neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids),
                             (5, capacity))

  large_sides = fill(typemin(Int), capacity)

  orientations = fill(typemin(Int), capacity)

  return L2MortarContainer3D{RealT, NVARS, POLYDEG}(
    u_upper_left, u_upper_right,
    u_lower_left, u_lower_right,
    neighbor_ids, large_sides, orientations,
    _u_upper_left, _u_upper_right,
    _u_lower_left, _u_lower_right,
    _neighbor_ids)
end

# Return number of L2 mortars
nmortars(l2mortars::L2MortarContainer3D) = length(l2mortars.orientations)

# Allow printing container contents
function Base.show(io::IO, ::MIME"text/plain", c::L2MortarContainer3D)
  println(io, '*'^20)
  for idx in CartesianIndices(c.u_upper_left)
    println(io, "c.u_upper_left[$idx] = $(c.u_upper_left[idx])")
  end
  for idx in CartesianIndices(c.u_upper_right)
    println(io, "c.u_upper_right[$idx] = $(c.u_upper_right[idx])")
  end
  for idx in CartesianIndices(c.u_lower_left)
    println(io, "c.u_lower_left[$idx] = $(c.u_lower_left[idx])")
  end
  for idx in CartesianIndices(c.u_lower_right)
    println(io, "c.u_lower_right[$idx] = $(c.u_lower_right[idx])")
  end
  println(io, "transpose(c.neighbor_ids) = $(transpose(c.neighbor_ids))")
  println(io, "c.large_sides = $(c.large_sides)")
  println(io, "c.orientations = $(c.orientations)")
  print(io,   '*'^20)
end

# Create mortar container and initialize mortar data in `elements`.
function init_mortars(cell_ids, mesh::TreeMesh3D,
                      elements::ElementContainer3D{RealT, NVARS, POLYDEG},
                      mortar::LobattoLegendreMortarL2) where {RealT<:Real, NVARS, POLYDEG}
  # Initialize containers
  n_mortars = count_required_mortars(mesh, cell_ids)
  mortars = L2MortarContainer3D{RealT, NVARS, POLYDEG}(n_mortars)

  # Connect elements with mortars
  init_mortars!(mortars, elements, mesh)
  return mortars
end

# Count the number of mortars that need to be created
function count_required_mortars(mesh::TreeMesh3D, cell_ids)
  count = 0

  # Iterate over all cells and count mortars from perspective of coarse cells
  for cell_id in cell_ids
    for direction in eachdirection(mesh.tree)
      # If no neighbor exists, cell is small with large neighbor or at boundary -> do nothing
      if !has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # If neighbor has no children, this is a conforming interface -> do nothing
      neighbor_id = mesh.tree.neighbor_ids[direction, cell_id]
      if !has_children(mesh.tree, neighbor_id)
        continue
      end

      count +=1
    end
  end

  return count
end

# Initialize connectivity between elements and mortars
function init_mortars!(mortars, elements, mesh::TreeMesh3D)
  # Construct cell -> element mapping for easier algorithm implementation
  tree = mesh.tree
  c2e = zeros(Int, length(tree))
  for element in eachelement(elements)
    c2e[elements.cell_ids[element]] = element
  end

  # Reset interface count
  count = 0

  # Iterate over all elements to find neighbors and to connect via interfaces
  for element in eachelement(elements)
    # Get cell id
    cell_id = elements.cell_ids[element]

    for direction in eachdirection(mesh.tree)
      # If no neighbor exists, cell is small with large neighbor -> do nothing
      if !has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # If neighbor has no children, this is a conforming interface -> do nothing
      neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id]
      if !has_children(mesh.tree, neighbor_cell_id)
        continue
      end

      # Create mortar between elements (3 possible orientations):
      #
      # mortar in x-direction:
      # 1 -> small element in lower, left position  (-y, -z)
      # 2 -> small element in lower, right position (+y, -z)
      # 3 -> small element in upper, left position  (-y, +z)
      # 4 -> small element in upper, right position (+y, +z)
      #
      # mortar in y-direction:
      # 1 -> small element in lower, left position  (-x, -z)
      # 2 -> small element in lower, right position (+x, -z)
      # 3 -> small element in upper, left position  (-x, +z)
      # 4 -> small element in upper, right position (+x, +z)
      #
      # mortar in z-direction:
      # 1 -> small element in lower, left position  (-x, -y)
      # 2 -> small element in lower, right position (+x, -y)
      # 3 -> small element in upper, left position  (-x, +y)
      # 4 -> small element in upper, right position (+x, +y)
      #
      # Always the case:
      # 5 -> large element
      #
      count += 1
      mortars.neighbor_ids[5, count] = element

      # Directions are from the perspective of the large element
      # ("Where are the small elements? Ah, in the ... direction!")
      if direction == 1 # -x
        mortars.neighbor_ids[1, count] = c2e[mesh.tree.child_ids[2, neighbor_cell_id]]
        mortars.neighbor_ids[2, count] = c2e[mesh.tree.child_ids[4, neighbor_cell_id]]
        mortars.neighbor_ids[3, count] = c2e[mesh.tree.child_ids[6, neighbor_cell_id]]
        mortars.neighbor_ids[4, count] = c2e[mesh.tree.child_ids[8, neighbor_cell_id]]
      elseif direction == 2 # +x
        mortars.neighbor_ids[1, count] = c2e[mesh.tree.child_ids[1, neighbor_cell_id]]
        mortars.neighbor_ids[2, count] = c2e[mesh.tree.child_ids[3, neighbor_cell_id]]
        mortars.neighbor_ids[3, count] = c2e[mesh.tree.child_ids[5, neighbor_cell_id]]
        mortars.neighbor_ids[4, count] = c2e[mesh.tree.child_ids[7, neighbor_cell_id]]
      elseif direction == 3 # -y
        mortars.neighbor_ids[1, count] = c2e[mesh.tree.child_ids[3, neighbor_cell_id]]
        mortars.neighbor_ids[2, count] = c2e[mesh.tree.child_ids[4, neighbor_cell_id]]
        mortars.neighbor_ids[3, count] = c2e[mesh.tree.child_ids[7, neighbor_cell_id]]
        mortars.neighbor_ids[4, count] = c2e[mesh.tree.child_ids[8, neighbor_cell_id]]
      elseif direction == 4 # +y
        mortars.neighbor_ids[1, count] = c2e[mesh.tree.child_ids[1, neighbor_cell_id]]
        mortars.neighbor_ids[2, count] = c2e[mesh.tree.child_ids[2, neighbor_cell_id]]
        mortars.neighbor_ids[3, count] = c2e[mesh.tree.child_ids[5, neighbor_cell_id]]
        mortars.neighbor_ids[4, count] = c2e[mesh.tree.child_ids[6, neighbor_cell_id]]
      elseif direction == 5 # -z
        mortars.neighbor_ids[1, count] = c2e[mesh.tree.child_ids[5, neighbor_cell_id]]
        mortars.neighbor_ids[2, count] = c2e[mesh.tree.child_ids[6, neighbor_cell_id]]
        mortars.neighbor_ids[3, count] = c2e[mesh.tree.child_ids[7, neighbor_cell_id]]
        mortars.neighbor_ids[4, count] = c2e[mesh.tree.child_ids[8, neighbor_cell_id]]
      elseif direction == 6 # +z
        mortars.neighbor_ids[1, count] = c2e[mesh.tree.child_ids[1, neighbor_cell_id]]
        mortars.neighbor_ids[2, count] = c2e[mesh.tree.child_ids[2, neighbor_cell_id]]
        mortars.neighbor_ids[3, count] = c2e[mesh.tree.child_ids[3, neighbor_cell_id]]
        mortars.neighbor_ids[4, count] = c2e[mesh.tree.child_ids[4, neighbor_cell_id]]
      else
        error("should not happen")
      end

      # Set large side, which denotes the direction (1 -> negative, 2 -> positive) of the large side
      if direction in (2, 4, 6)
        mortars.large_sides[count] = 1
      else
        mortars.large_sides[count] = 2
      end

      # Set orientation (x -> 1, y -> 2, z -> 3)
      if direction in (1, 2)
        mortars.orientations[count] = 1
      elseif direction in (3, 4)
        mortars.orientations[count] = 2
      else
        mortars.orientations[count] = 3
      end
    end
  end

  @assert count == nmortars(mortars) ("Actual mortar count ($count) does not match " *
                                      "expectations $(nmortars(mortars))")
end


text:
using HalfSine, Test

@test 1==1

ext: jl


text:
################################################################################
#
#  Type
#
################################################################################

embedding_type(::Type{AnticNumberField}) = NumFieldEmbNfAbs

embedding_type(::AnticNumberField) = embedding_type(AnticNumberField)

################################################################################
#
#  Field access
#
################################################################################

number_field(f::NumFieldEmbNfAbs) = f.K

isreal(f::NumFieldEmbNfAbs) = f.isreal

_absolute_index(f::NumFieldEmbNfAbs) = f.i

################################################################################
#
#  Construction
#
################################################################################

function complex_embeddings(K::AnticNumberField; conjugates::Bool = true)
  res = get_attribute!(K, :complex_embeddings) do
    return __complex_embeddings(K)
  end::Vector{embedding_type(K)}
  if conjugates
    return res
  else
    r, s = signature(K)
    return res[1:(r + s)]
  end
end

function __complex_embeddings(K::AnticNumberField)
  c = conjugate_data_arb(K)
  res = Vector{embedding_type(K)}(undef, degree(K))
  for i in 1:degree(K)
    res[i] = NumFieldEmbNfAbs(K, c, i)
  end
  return res
end

################################################################################
#
#  Conjugate
#
################################################################################

function conj(f::NumFieldEmbNfAbs)
  return complex_embeddings(number_field(f))[f.conjugate]
end

################################################################################
#
#  Equality
#
################################################################################

function Base.:(==)(f::NumFieldEmbNfAbs, g::NumFieldEmbNfAbs)
  return number_field(f) === number_field(g) &&
      _absolute_index(f) == _absolute_index(g)
end

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

function Base.show(io::IO, ::MIME"text/plain", f::NumFieldEmbNfAbs)
  print(io, "Embedding of\n", number_field(f), "\ncorresponding to ≈ ")
  _print_acb_neatly(io, f.r)
end

function Base.show(io::IO, f::NumFieldEmbNfAbs)
  print(io, "Embedding corresponding to ≈ ")
  _print_acb_neatly(io, f.r)
end

################################################################################
#
#  Evaluation
#
################################################################################

evaluate(x::nf_elem, f::NumFieldEmbNfAbs, p::Int) = f(x, p)

function (f::NumFieldEmbNfAbs)(x::nf_elem, abs_tol::Int = 32)
  K = parent(x)
  d = degree(K)
  r1, r2 = signature(K)
  target_tol = abs_tol
  abs_tol = Int(floor(abs_tol * 1.1))

  i = f.i

  while true
    prec_too_low = false
    c = conjugate_data_arb_roots(K, abs_tol)

    abs_tol < 0 && error("Precision overflow in evaluation of embedding")

    CC = AcbField(abs_tol, cached = false)
    RR = ArbField(abs_tol, cached = false)

    xpoly = arb_poly(parent(K.pol)(x), abs_tol)

    if i <= r1
      o = RR()
      ccall((:arb_poly_evaluate, libarb), Nothing,
            (Ref{arb}, Ref{arb_poly}, Ref{arb}, Int),
             o, xpoly, c.real_roots[i], abs_tol)

      if !isfinite(o) || !radiuslttwopower(o, -target_tol)
        abs_tol = 2*abs_tol
        prec_too_low = true
      else
        cc = CC(o)
        break
      end
    else
      tacb = CC()
      j = i <= r1 + r2 ? i - r1 : i - r1 - r2
      ccall((:arb_poly_evaluate_acb, libarb), Nothing,
            (Ref{acb}, Ref{arb_poly}, Ref{acb}, Int),
             tacb, xpoly, c.complex_roots[j], abs_tol)

      if !isfinite(tacb) || !radiuslttwopower(tacb, -target_tol)
        abs_tol = 2*abs_tol
        prec_too_low = true
      else
        if i <= r1 + r2
          cc = tacb
        else
          cc = conj(tacb)
        end
        break
      end
    end

    if prec_too_low
      continue
    end
  end

  cc = expand!(cc, -target_tol)
  return cc
end

#################################################################################
#
#  Restrict from bigger field
#
################################################################################

function restrict(e::NumFieldEmb, f::NumFieldMor{AnticNumberField, <: Any, <: Any})
  @req number_field(e) === codomain(f) "Number fields do not match"
  L = domain(f)
  emb = complex_embeddings(L)
  a = gen(L)
  # I need to find the embedding of the base_field of L
  b = e(f(a))
  cn = Bool[overlaps(b, e(a)) for e in emb]
  @assert count(cn) == 1
  i = findfirst(cn)
  return emb[i]
end


text:
using RandomLineChart
using Test

@testset "RandomLineChart.jl" begin
    # Write your tests here.
end


text:
"""
    LocationScale(μ,σ,ρ)

A location-scale transformed distribution with location parameter `μ`,
scale parameter `σ`, and given distribution `ρ`.

```math
f(x) = \\frac{1}{σ} ρ \\! \\left( \\frac{x-μ}{σ} \\right)
```

```julia
LocationScale(μ,σ,ρ) # location-scale transformed distribution
params(d)            # Get the parameters, i.e. (μ, σ, and the base distribution)
location(d)          # Get the location parameter
scale(d)             # Get the scale parameter
```

External links
[Location-Scale family on Wikipedia](https://en.wikipedia.org/wiki/Location%E2%80%93scale_family)
"""
struct LocationScale{T<:Real, D<:ContinuousUnivariateDistribution} <: ContinuousUnivariateDistribution
    μ::T
    σ::T
    ρ::D
    LocationScale{T, D}(μ::T,σ::T,ρ::D) where {T<:Real, D<:ContinuousUnivariateDistribution} = new{T, D}(μ, σ, ρ)
end

function LocationScale(μ::T,σ::T,ρ::D) where {T<:Real, D<:ContinuousUnivariateDistribution}
    @check_args(LocationScale, σ > zero(σ))
    return LocationScale{T,D}(μ,σ,ρ)
end

function LocationScale(μ::T,σ::T,ρ::D, ::NoArgCheck) where {T<:Real, D<:ContinuousUnivariateDistribution}
    return LocationScale{T,D}(μ,σ,ρ)
end

function LocationScale(μ::Integer, σ::Integer, ρ::ContinuousUnivariateDistribution)
    return LocationScale(float(μ), float(σ), ρ)
end

LocationScale(μ::Real, σ::Real, ρ::D) where {D<:ContinuousUnivariateDistribution} = LocationScale(promote(μ,σ)...,ρ)

minimum(d::LocationScale) = d.μ + d.σ * minimum(d.ρ)
maximum(d::LocationScale) = d.μ + d.σ * maximum(d.ρ)

#### Conversions

convert(::Type{LocationScale{T}}, μ::Real, σ::Real, ρ::D) where {T<:Real, D<:ContinuousUnivariateDistribution} = LocationScale(T(μ),T(σ),ρ)
convert(::Type{LocationScale{T}}, d::LocationScale{S}) where {T<:Real, S<:Real} = LocationScale(T(d.μ),T(d.σ),d.ρ, NoArgCheck())

#### Parameters

location(d::LocationScale) = d.μ
scale(d::LocationScale) = d.σ
params(d::LocationScale) = (d.μ,d.σ,d.ρ)
@inline partype(d::LocationScale{T}) where {T<:Real} = T

#### Statistics

mean(d::LocationScale) = d.μ + d.σ * mean(d.ρ)
median(d::LocationScale) = d.μ + d.σ * median(d.ρ)
mode(d::LocationScale) = d.μ + d.σ * mode(d.ρ)
modes(d::LocationScale) = d.μ .+ d.σ .* modes(d.ρ)

var(d::LocationScale) = d.σ^2 * var(d.ρ)
std(d::LocationScale) = d.σ * std(d.ρ)
skewness(d::LocationScale) = skewness(d.ρ)
kurtosis(d::LocationScale) = kurtosis(d.ρ)

isplatykurtic(d::LocationScale) = isplatykurtic(d.ρ)
isleptokurtic(d::LocationScale) = isleptokurtic(d.ρ)
ismesokurtic(d::LocationScale) = ismesokurtic(d.ρ)

entropy(d::LocationScale) = entropy(d.ρ) + log(d.σ)
mgf(d::LocationScale,t::Real) = exp(d.μ*t) * mgf(d.ρ,d.σ*t)

#### Evaluation & Sampling

pdf(d::LocationScale,x::Real) = pdf(d.ρ,(x-d.μ)/d.σ) / d.σ
logpdf(d::LocationScale,x::Real) = logpdf(d.ρ,(x-d.μ)/d.σ) - log(d.σ)
cdf(d::LocationScale,x::Real) = cdf(d.ρ,(x-d.μ)/d.σ)
logcdf(d::LocationScale,x::Real) = logcdf(d.ρ,(x-d.μ)/d.σ)
quantile(d::LocationScale,q::Real) = d.μ + d.σ * quantile(d.ρ,q)

rand(rng::AbstractRNG, d::LocationScale) = d.μ + d.σ * rand(rng, d.ρ)
cf(d::LocationScale, t::Real) = cf(d.ρ,t*d.σ) * exp(1im*t*d.μ)
gradlogpdf(d::LocationScale, x::Real) = gradlogpdf(d.ρ,(x-d.μ)/d.σ) / d.σ


text:
# The Fibonacci tiling

using SubstitutionTilings
using SubstitutionTilings.NumFields
using StructEquality
using Luxor
using Plots

@NumFields.simple_number_field_concrete Qτ [1,1] τ
Base.promote_rule(::Type{Qτ}, ::Type{<:Integer}) = Qτ

@def_structequal struct FibElem <: DGroupElem
    a :: Qτ
end

@enum FibTile begin
    A=1
    B=2
end

function Base.:*(x :: FibElem, y :: FibElem)
    return FibElem(x.a + y.a)
end

function SubstitutionTilings.dilate(λ :: Qτ, x :: FibElem)
    return FibElem(λ*x.a)
end

function Base.inv(x :: FibElem)
    return FibElem(-x.a)
end

function SubstitutionTilings.id(::Type{FibElem})
    return FibElem(0)
end

function embed_float(g :: FibElem)
    fτ = (1+sqrt(5))/2
    return embed_field(float,fτ, g.a)
end

function SubstitutionTilings.embed_aff(g :: FibElem)
    return [1, 0, 0, 1, embed_float(g), 0]
end

function SubstitutionTilings.draw(ptile::FibTile, action)
    fτ = (1+sqrt(5))/2
    if ptile == A
        return Luxor.poly([
            Point(-fτ*0.9,-10),
            Point(-fτ*0.9,10),
            Point(fτ*0.9,10),
            Point(fτ*0.9,-10)
        ], close = true, action)
    else
        return Luxor.poly([
            Point(-0.9,-10),
            Point(-0.9,10),
            Point(0.9,10),
            Point(0.9,-10)
        ], close = true, action)
    end
end

fib = SubSystem(Dict(A => [FibElem(Qτ(-1)//2) => A, FibElem(τ//2) => B], B => [FibElem(0) => A]),τ)

fib_tiling = substitute(fib, Dict([FibElem(0) => A]), 3)

function color(tile :: Pair{FibElem,FibTile})
    return (tile[2] == A) ? 1 : 2
end

width = 800
height = 20
sc = 20
@png begin
    colors = ["#DD93FC", "#E7977A"]
    patch = [FibElem(-τ//2) => A, FibElem(τ//2) => A]
    tiling = substitute(fib, patch, 4)
    setline(1)

    for tile in tiling
        origin()
        draw(tile, sc, colors[color(tile)], :fill)
    end

    origin()
    scale(sc)
    sethue("black")
    Luxor.poly([
            Point(-0.05,-10),
            Point(-0.05,10),
            Point(0.05,10),
            Point(0.05,-10)
        ], close = true, :fill)
end width height

function SubstitutionTilings.is_interior(tiling :: Dict, t :: FibElem)
    return haskey(tiling, t) && (haskey(tiling, FibElem(t.a-τ)) || haskey(tiling, FibElem(t.a-(1+τ)//2))) && (haskey(tiling, FibElem(t.a+τ)) || haskey(tiling, FibElem(t.a+(1+τ)//2)))
end

function SubstitutionTilings.collar_in(tiling :: Dict, t :: FibElem)
    if !is_interior(tiling, t)
        throw(UnrecognizedCollar)
    end
    collar = Dict{FibElem, FibTile}([])
    shifts = ([0, -τ, -(1+τ)//2, τ, (1+τ)//2])
    for shift in shifts
        if haskey(tiling, FibElem(t.a+shift))
            collar[FibElem(t.a+shift)] = tiling[FibElem(t.a+shift)]
        end
    end
    return collar
end

initial_collar = collar_in(fib_tiling, FibElem(0))
(collars, fib_c) = total_collaring(fib, initial_collar)

collars

width = 800
height = 40
sc = 20
@png begin
    colors = ["#DD93FC", "#E7977A"]
    first_tile = [FibElem(0) => A]
    tiling = collars[4]
    setline(1)

    for tile in tiling
        origin()
        draw(tile, sc, colors[color(tile)], :fill)
    end
end width height

frequency(fib, initial_collar, collars[4], 4)
nu = autocorrelation(fib, initial_collar, 9)
xs = Vector{Float64}()
ys = Vector{Float64}()
for delta in nu
    if abs(embed_float(delta[1])) < 10
        push!(xs, embed_float(delta[1]))
        push!(ys, delta[2])
    end
end
plot(xs, ys, seriestype = :scatter, ylims = (0,1.1))


text:
using Distributions
N=Normal(-0.940,7.38);
function f(x)
  return pdf(N,x);
end
(P1,err)=quadgk(f,20,100);
(P2,err)=quadgk(f,-100,-20);
P=P1+P2;
println(P);

(P3,err)=quadgk(f,-10,10);

println(P3);


text:
function seedSimulation!( mpSim::MPsim, seed::Integer, sysEnt::Bool )

    # Seed the simulation.
    mpSim.seedRNG = seed < 0 ? MersenneTwister() : MersenneTwister( seed )

    for name in sort( collect( keys( mpSim.recruitmentByName ) ) )
        for recruitment in mpSim.recruitmentByName[name]
            recSeed = sysEnt ? nothing : rand( mpSim.seedRNG, UInt64 )
            recruitment.recRNG = MersenneTwister( recSeed )
            recDist = recruitmentDists[recruitment.recruitmentDistType][3]
            recDistNodes = recruitment.recruitmentDistNodes
            recruitment.recruitmentDist = recDist( recruitment, recDistNodes,
                sort( collect( keys( recDistNodes ) ) ) )

            ageSeed = sysEnt ? nothing : rand( mpSim.seedRNG, UInt64 )
            recruitment.ageRNG = MersenneTwister( ageSeed )
            ageDist = recruitmentDists[recruitment.ageDistType][4]
            ageDistNodes = recruitment.ageDistNodes
            recruitment.ageDist = ageDist( recruitment, ageDistNodes,
                sort( collect( keys( ageDistNodes ) ) ) )
        end  # for recruitment in mpSim.recruitmentByName[name]
    end  # for name in keys( mpSim.recruitmentByName )

    # Seed the attributes.
    for name in sort( collect( keys( mpSim.attributeList ) ) )
        attSeed = sysEnt ? nothing : rand( mpSim.seedRNG, UInt64 )
        attribute = mpSim.attributeList[name]
        attribute.initRNG = MersenneTwister( attSeed )
    end  # for name in keys( mpSim.attributeList )

    # Seed the attrition times.
    for name in sort( collect( keys( mpSim.attritionSchemes ) ) )
        attSeed = sysEnt ? nothing : rand( mpSim.seedRNG, UInt64 )
        attrition = mpSim.attritionSchemes[name]
        attrition.timeRNG = MersenneTwister( attSeed )
    end  # for name in keys( mpSim.attritionSchemes )

    # Seed the transitions.
    for name in sort( collect( keys( mpSim.transitionsByName ) ) )
        for transition in mpSim.transitionsByName[name]
            transSeed = sysEnt ? nothing : rand( mpSim.seedRNG, UInt64 )
            transition.probRNG = MersenneTwister( transSeed )
        end  # for transition in mpSim.transitionsByName[name]
    end  # for name in keys( mpSim.transitionsByName )

end  # seedSimulation!( mpSim, seed, sysEnt )

setSimulationDBcommits!( mpSim::MPsim, nCommits::Integer ) =
    ( mpSim.nCommits = nCommits )

function wipeDatabase( mpSim::MPsim )

    # Wipe the tables if needed.
    DBInterface.execute( mpSim.simDB, string( "DROP TABLE IF EXISTS `",
        mpSim.histDBname, "`" ) )
    DBInterface.execute( mpSim.simDB, string( "DROP TABLE IF EXISTS `",
        mpSim.transDBname, "`" ) )
    DBInterface.execute( mpSim.simDB, string( "DROP TABLE IF EXISTS `",
        mpSim.persDBname, "`" ) )

    # Reset the population number of the simulation.
    mpSim.orgSize = 0
    mpSim.dbSize = 0
    mpSim.isVirgin = true

    # Create the tables.
    sqliteCmd = string( "CREATE TABLE `", mpSim.persDBname, "`(",
        "\n    `", mpSim.idKey, "` TEXT NOT NULL PRIMARY KEY,",
        "\n    timeEntered REAL,",
        "\n    timeExited REAL,",
        "\n    ageAtRecruitment REAL,",
        "\n    expectedAttritionTime REAL,",
        "\n    currentNode TEXT,",
        "\n    inNodeSince REAL,",
        join( string.( "\n    `", collect( keys( mpSim.attributeList ) )
            , "` TEXT," ) ),
        "\n    status TEXT )" )
    DBInterface.execute( mpSim.simDB, sqliteCmd )

    sqliteCmd = string( "CREATE TABLE `", mpSim.transDBname, "`(",
        "\n    `", mpSim.idKey, "` TEXT,",
        "\n    timeIndex REAL,",
        "\n    transition TEXT,",
        "\n    sourceNode TEXT,",
        "\n    targetNode TEXT,",
        "\n    FOREIGN KEY (`", mpSim.idKey, "`) REFERENCES `",
        mpSim.persDBname, "`(`", mpSim.idKey, "`) )" )
    DBInterface.execute( mpSim.simDB, sqliteCmd )

    sqliteCmd = string( "CREATE TABLE `", mpSim.histDBname, "`(",
        "\n    `", mpSim.idKey, "` TEXT,",
        "\n    timeIndex REAL,",
        "\n    attribute TEXT,",
        "\n    value TEXT,",
        "\n    FOREIGN KEY (`", mpSim.idKey, "`) REFERENCES `",
        mpSim.persDBname, "`(`", mpSim.idKey, "`) )" )
    DBInterface.execute( mpSim.simDB, sqliteCmd )

end  # wipeDatabase( mpSim )

function wipeDatabase( mpSim::MPsim, snapshotIDs::Vector{String} )

    # Reset the population number of the simulation.
    mpSim.orgSize = length(snapshotIDs)
    mpSim.dbSize = length(snapshotIDs)
    mpSim.isVirgin = false

    # Reset the tables.
    sqliteCmd = string( "DELETE FROM `", mpSim.transDBname,
        "` WHERE transition IS NOT 'Init'" )
    DBInterface.execute( mpSim.simDB, sqliteCmd )

    sqliteCmd = string( "DELETE FROM `", mpSim.histDBname,
        "` WHERE timeIndex > 0 OR `", mpSim.idKey, "` NOT IN ('",
        join( snapshotIDs, "', '" ), "')" )
    DBInterface.execute( mpSim.simDB, sqliteCmd )

    sqliteCmd = string( "DELETE FROM `", mpSim.persDBname,
        "` WHERE `", mpSim.idKey, "` NOT IN ('", join( snapshotIDs, "', '" ),
        "')" )
    DBInterface.execute( mpSim.simDB, sqliteCmd )

    # Set the "inNodeSince" for every node.
    queryCmd = string( "SELECT `", mpSim.idKey,
        "`, timeIndex, targetNode FROM `", mpSim.transDBname, "`" )
    inNodeInfo = DataFrame(DBInterface.execute( mpSim.simDB, queryCmd ))

    for ii in 1:size( inNodeInfo, 1 )
        node = mpSim.baseNodeList[inNodeInfo[ii, "targetNode"]]
        node.inNodeSince[inNodeInfo[ii, mpSim.idKey]] =
            inNodeInfo[ii, "timeIndex"]
    end  # for ii in 1:size( inNodeInfo, 1 )

end  # wipeDatabase( mpSim, snapshotIDs )

function resetSimulation( mpSim::MPsim, keepSnap::Bool=true )

    mpSim.sim = Simulation()

    snapshotIDs = Vector{String}()
    dbtables = SQLite.tables(mpSim.simDB)

    if keepSnap && haskey( dbtables, :name ) &&
        ( [mpSim.persDBname, mpSim.histDBname, mpSim.transDBname] ⊆
            dbtables.name )
        queryCmd = string( "SELECT `", mpSim.idKey, "` FROM `",
            mpSim.transDBname, "` WHERE transition IS 'Init'" )
        snapshotIDs = DataFrame(DBInterface.execute( mpSim.simDB,
            queryCmd ))[!, 1]
    end  # if haskey( dbtables, :name ) && ...

    # Clear the "inNodeSince" from every node.
    for name in keys( mpSim.baseNodeList )
        node = mpSim.baseNodeList[name]
        empty!( node.inNodeSince )
    end  # for name in keys( mpSim.baseNodeList )

    isempty(snapshotIDs) ? wipeDatabase(mpSim) :
        wipeDatabase( mpSim, snapshotIDs )

    # Set database to new style.
    mpSim.isOldDB = false
    mpSim.sNode = "sourceNode"
    mpSim.tNode = "targetNode"
    mpSim.valName = "value"

    # Reset timers.
    mpSim.simExecTime = Millisecond( 0 )
    mpSim.attritionExecTime = Millisecond( 0 )

end  # resetSimulation( mpSim )


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module ProfileViewSVG

function __init__()
    eval(Expr(:import, :ProfileView))
end

function view(data = Profile.fetch(); C = false, lidict = nothing, colorgc = true, fontsize = 12, combine = true, pruned = true)
    img, lidict, imgtags = ProfileView.prepare(data, C=C, lidict=lidict, colorgc=colorgc, combine=combine, pruned=pruned)
    ProfileView.ProfileData(img, lidict, imgtags, fontsize)
end

end


text:
include("../src/autodiff_module.jl")
include("../../__lib__/math/common/numder/src/numder_module.jl")

module dtest
using ..numder
using ..autodiff
using PyPlot
PyPlot.pygui(true)

t = autodiff.Tape()
x = autodiff.create_variable(t, pi)
y = autodiff.create_variable(t, 20.0)

# z = cos(x * y + sin(x))

z1 = x + y
z2 = x * y
z3 = exp(x)
z4 = sin(y)
z5 = abs(y)
z6 = z1^z2
z7 = 3 * z1
z8 = max(z1, z2)
z = z1 + z2 + z3 + z4 + z5 + z6 + z7

g = autodiff.grad(z)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# Adds the unprimed model states to the dictionary of ranges to properly
# fill the Jacobian matrix
# Note: This is not augmentation in the usual sense, which incorporates lags of
# model state variables after the transition equation has been solved for
# but rather, augmenting the model states with lags prior to solution since the
# solution method requires lags
function augment_model_states(endo::AbstractDict{Symbol, UnitRange}, n_model_states::Int64)
    endo_aug = deepcopy(endo)
    for (state::Symbol, inds::UnitRange) in endo
        unprimed_state = unprime(state)
        unprimed_inds::UnitRange  = reindex_unprimed_model_states(inds, n_model_states)
        endo_aug[unprimed_state] = unprimed_inds
    end

    # Ensure all ranges are consecutive
    endo_ranges = endo_aug.vals
  #=  for i in 1:length(endo_ranges)
        if i == 1
            @assert endo_ranges[1].start == 1
        else
            @assert endo_ranges[i].start == endo_ranges[i-1].stop + 1
        end
    end =#

    return endo_aug
end

# For the RANK equivalent of the HANK models
function augment_model_states(endo::AbstractDict{Symbol, Int64}, n_model_states::Int64)
    endo_aug = deepcopy(endo)
    for (state::Symbol, ind::Int) in endo
        unprimed_state           = unprime(state)
        endo_aug[unprimed_state] = ind + n_model_states
    end

    return endo_aug
end

# The reason for this function is that the canonical form for the Klein solution method is
# as follows:
# E_t A([x_{t+1}, y_{t+1}]) = B[x_t, y_t]
# Where we only need to track [x_{t+1}, y_{t+1}] as states when
# we transform the system from canonical form to state space form
# Hence because we only keep track of the t+1 indexed variables, which we denote
# with a ′, we need a way of calculating the indices in the Jacobian corresponding
# to the t indexed variables.
function reindex_unprimed_model_states(inds::UnitRange, n_model_states::Int64)
    return inds .+ n_model_states
end

# Returns a model state variable symbol without the prime
function unprime(state::Symbol)
    return Symbol(replace(string(state), "′" => ""))
end

######################################################################################
# For normalizing the states
function normalize_model_state_indices!(m::AbstractModel)
    normalized_model_states = m.normalized_model_states

    endo                     = m.endogenous_states
    model_state_keys         = collect(keys(endo))
    normalized_model_states = intersect(model_state_keys, normalized_model_states) #findall((in)(normalized_model_states), model_state_keys)

    state_inds                  = get_setting(m, :state_indices)
    jump_inds                   = get_setting(m, :jump_indices)
    state_normalization_factor  = get_setting(m, :backward_looking_states_normalization_factor)
    jump_normalization_factor   = get_setting(m, :jumps_normalization_factor)

    m.endogenous_states = normalize(m, endo, model_state_keys, normalized_model_states,
                                    state_inds, jump_inds,
                                    state_normalization_factor, jump_normalization_factor)
end

# Shift a UnitRange type down by an increment
# If this UnitRange is the first_range in the group being shifted, then do not
# subtract the increment from the start
# e.g.
# If you have 1:80 as your first range, then shift(1:80, -1; first_range = true)
# should return 1:79
# However, if you have 81:160 as range, that is not your first range, then the
# function should return 80:159
function shift(inds::UnitRange, increment::Int64; first_range::Bool = false)
    if first_range
        return UnitRange(inds.start, inds.stop + increment)
    else
        return UnitRange(inds.start + increment, inds.stop + increment)
    end
end

# normalize the distributional states located in indices specified by normalized_state_inds
function normalize(m::AbstractModel,
                   endo::AbstractDict{Symbol, UnitRange},
                   model_state_keys::Vector{Symbol},
                   normalized_model_states::Vector{Symbol},
                   state_inds::AbstractArray{Int64},
                   jump_inds::AbstractArray{Int64},
                   state_normalization_factor::Int64,
                   jump_normalization_factor::Int64)
    n_model_state_vars = length(endo)

    # For each model state that needs to be normalized...
    for state in normalized_model_states
        normalization_factor =
        if state in get_setting(m, :states)
            state_normalization_factor
        elseif state in get_setting(m, :jumps)
            jump_normalization_factor
        end
        # Subtract 1 from both the beginning and end of the UnitRange...
        for state in model_state_keys
            inds = endo[state]
            # except for first UnitRange in a given normalization
            first_range = (state==get_setting(m, :states)[1])
            endo[state] = shift(inds, -normalization_factor, first_range = first_range)
        end
        #=for j in i:n_model_state_vars
            inds = endo[model_state_keys[j]]
            # Except for the first UnitRange in a given normalization
            first_range = (j == i)
            endo[model_state_keys[j]] = shift(inds, -normalization_factor, first_range = first_range)
        end=#
    end

    # Ensure all ranges are consecutive
    endo_ranges = endo.vals
   #= for i in 1:n_model_state_vars
        if i == 1
            @assert endo_ranges[1].start == 1
        else
            @assert endo_ranges[i].start == endo_ranges[i-1].stop + 1
        end
    end =#

    return endo
end

"""
```
stack_indices(key_dict::OrderedDict, keys::Vector{Symbol})
```
Stacks the indices of a OrderedDict{Symbol, AbstractRange} into a Vector of Ints defined by the Ranges in the dictionary.
"""
function stack_indices(key_dict::AbstractDict, keys::Vector{Symbol})
    indices = Vector{Int64}(undef, 0)
    for i in getindex.(Ref(key_dict), keys)
        indices = [indices; i]
    end
    return indices
end


text:
;; ---------------------------------------------------------------------
;;	csv_create.jl
;;
;;					Jun/06/2013
;;
;; --------------------------------------------------------------------
(require 'tables)
(load-file "/var/www/data_base/common/rep_common/text_manipulate.jl")
;; --------------------------------------------------------------------

;; --------------------------------------------------------------------
(defun data_prepare_proc ()
	(setq dict_aa (make-table string-hash string=))
;
	(dict_append_proc dict_aa "t1271" "千葉" "78361" "1930-8-17")
	(dict_append_proc dict_aa "t1272" "勝浦" "41657" "1930-10-21")
	(dict_append_proc dict_aa "t1273" "市原" "94523" "1930-5-16")
	(dict_append_proc dict_aa "t1274" "流山" "51482" "1930-9-23")
	(dict_append_proc dict_aa "t1275" "八千代" "47851" "1930-8-9")
	(dict_append_proc dict_aa "t1276" "我孫子" "82753" "1930-5-19")
	(dict_append_proc dict_aa "t1277" "鴨川" "51748" "1930-9-12")
	(dict_append_proc dict_aa "t1278" "銚子" "23689" "1930-1-18")
	(dict_append_proc dict_aa "t1279" "市川" "91543" "1930-5-15")
)

;; --------------------------------------------------------------------
(format t "*** 開始 ***")
(setq file_out (cadr command-line-args))
(format t file_out)

(setq dict_aa (data_prepare_proc))
(table-walk dict_display_proc_single dict_aa)

(csv_write_proc file_out dict_aa)

(format t "*** 終了 ***")

;; ---------------------------------------------------------------------

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
"""
Write results into a csv file. This function takes as input the results from the
`learn_paths` and `build_paths` functions, including the information on gold paths
that is optionally returned as second output result.
"""
function write2csv end

"""
Reformat results into a dataframe. This function takes as input the results from the
`learn_paths` and `build_paths` functions, including the information on gold paths
that is optionally returned as second output result.
"""
function write2df end

"""
Write comprehension evaluation into a CSV file, include target and predicted
ids and indentifiers and their correlations.
"""
function write_comprehension_eval end

"""
    write2csv(res, data, cue_obj_train, cue_obj_val, filename)

Write results into csv file for the results from `learn_paths` and `build_paths`.

# Obligatory Arguments
- `res::Array{Array{Result_Path_Info_Struct,1},1}`: the results from `learn_paths` or `build_paths`
- `data::DataFrame`: the dataset
- `cue_obj_train::Cue_Matrix_Struct`: the cue object for training dataset
- `cue_obj_val::Cue_Matrix_Struct`: the cue object for validation dataset
- `filename::String`: the filename

# Optional Arguments
- `grams::Int64=3`: the number n in n-gram cues
- `tokenized::Bool=false`: if true, the dataset target is tokenized
- `sep_token::Union{Nothing, String, Char}=nothing`: separator
- `start_end_token::Union{String, Char}="#"`: start and end token in boundary cues
- `output_sep_token::Union{String, Char}=""`: output separator
- `path_sep_token::Union{String, Char}=":"`: path separator
- `target_col::Union{String, Symbol}=:Words`: the column name for target strings
- `root_dir::String="."`: dir path for project root dir
- `output_dir::String="."`: output dir inside root dir

# Examples
```julia
# writing results for training data
JudiLing.write2csv(
    res_train,
    latin_train,
     cue_obj_train,
    cue_obj_train,
    "res_latin_train.csv",
    grams=3,
    tokenized=false,
    sep_token=nothing,
    start_end_token="#",
    output_sep_token="",
    path_sep_token=":",
    target_col=:Word,
    root_dir=".",
    output_dir="test_out")

# writing results for validation data
JudiLing.write2csv(
    res_val,
    latin_val,
    cue_obj_train,
    cue_obj_val,
    "res_latin_val.csv",
    grams=3,
    tokenized=false,
    sep_token=nothing,
    start_end_token="#",
    output_sep_token="",
    path_sep_token=":",
    target_col=:Word,
    root_dir=".",
    output_dir="test_out")
```
"""
function write2csv(
    res,
    data,
    cue_obj_train,
    cue_obj_val,
    filename;
    grams = 3,
    tokenized = false,
    sep_token = nothing,
    start_end_token = "#",
    output_sep_token = "",
    path_sep_token = ":",
    target_col = :Words,
    root_dir = ".",
    output_dir = ".",
)

    output_path = joinpath(root_dir, output_dir)
    # create path if not exist
    mkpath(output_path)

    if isnothing(output_sep_token)
        output_sep_token = ""
    end

    df = JudiLing.write2df(
        res,
        data,
        cue_obj_train,
        cue_obj_val,
        grams = grams,
        tokenized = tokenized,
        sep_token = sep_token,
        start_end_token = start_end_token,
        output_sep_token = output_sep_token,
        path_sep_token = path_sep_token,
        target_col = target_col,
    )

    CSV.write(joinpath(output_path, filename), df, quotestrings = true)
    nothing
end

"""
    write2csv(gpi::Vector{Gold_Path_Info_Struct}, filename)

Write results into csv file for the gold paths' information optionally returned by
`learn_paths` and `build_paths`.

# Obligatory Arguments
- `gpi::Vector{Gold_Path_Info_Struct}`: the gold paths' information
- `filename::String`: the filename

# Optional Arguments
- `root_dir::String="."`: dir path for project root dir
- `output_dir::String="."`: output dir inside root dir

# Examples
```julia
# write gold standard paths to csv for training data
JudiLing.write2csv(
    gpi_train,
    "gpi_latin_train.csv",
    root_dir=".",
    output_dir="test_out"
    )

# write gold standard paths to csv for validation data
JudiLing.write2csv(
    gpi_val,
    "gpi_latin_val.csv",
    root_dir=".",
    output_dir="test_out"
    )
```
"""
function write2csv(gpi::Vector{Gold_Path_Info_Struct}, filename; root_dir = ".", output_dir = ".")
    output_path = joinpath(root_dir, output_dir)
    # create path if not exist
    mkpath(output_path)
    # open a file
    io = open(joinpath(output_path, filename), "w")

    # write header
    write(
        io,
        "\"utterance\",\"weakest_support\",\"weakest_support_timestep\",\"support\",\"gold_path\",\"timestep_support\"\n",
    )

    for (i, g) in enumerate(gpi)
        if isempty(g.ngrams_ind_support)
            ws = missing
            wst = missing
        else
            ws, wst = findmin(g.ngrams_ind_support)
        end
        write(
            io,
            "\"$i\",\"$(ws)\",\"$(wst)\",\"$(g.support)\",\"$(g.ngrams_ind)\",\"$(g.ngrams_ind_support)\"\n",
        )
    end
    # close file
    close(io)
end

"""
    write2csv(ts::Threshold_Stat_Struct, filename)

Write results into csv file for threshold and tolerance proportion for each timestep.

# Obligatory Arguments
- `gpi::Vector{Gold_Path_Info_Struct}`: the gold paths' information
- `filename::String`: the filename

# Optional Arguments
- `root_dir::String="."`: dir path for project root dir
- `output_dir::String="."`: output dir inside root dir

# Examples
```julia
JudiLing.write2csv(ts, "ts.csv", root_dir = @__DIR__, output_dir="out")
```
"""
function write2csv(ts::Threshold_Stat_Struct, filename; root_dir = ".", output_dir = ".")
    output_path = joinpath(root_dir, output_dir)

    # create path if not exist
    mkpath(output_path)

    df = JudiLing.write2df(ts)
    CSV.write(joinpath(output_path, filename), df, quotestrings = true)
    nothing
end

"""
    write2df(res, data, cue_obj_train, cue_obj_val)

Reformat results into a dataframe for the results form `learn_paths` and `build_paths`
functions.

# Obligatory Arguments
- `data::DataFrame`: the dataset

# Optional Arguments
- `grams::Int64=3`: the number n in n-gram cues
- `tokenized::Bool=false`: if true, the dataset target is tokenized
- `sep_token::Union{Nothing, String, Char}=nothing`: separator
- `start_end_token::Union{String, Char}="#"`: start and end token in boundary cues
- `output_sep_token::Union{String, Char}=""`: output separator
- `path_sep_token::Union{String, Char}=":"`: path separator
- `target_col::Union{String, Symbol}=:Words`: the column name for target strings

# Examples
```julia
# writing results for training data
JudiLing.write2df(
    res_train,
    latin_train,
    cue_obj_train,
    cue_obj_train,
    grams=3,
    tokenized=false,
    sep_token=nothing,
    start_end_token="#",
    output_sep_token="",
    path_sep_token=":",
    target_col=:Word)

# writing results for validation data
JudiLing.write2df(
    res_val,
    latin_val,
    cue_obj_train,
    cue_obj_val,
    grams=3,
    tokenized=false,
    sep_token=nothing,
    start_end_token="#",
    output_sep_token="",
    path_sep_token=":",
    target_col=:Word)
```
"""
function write2df(
    res,
    data,
    cue_obj_train,
    cue_obj_val;
    grams = 3,
    tokenized = false,
    sep_token = nothing,
    start_end_token = "#",
    output_sep_token = "",
    path_sep_token = ":",
    target_col = :Words,
)

    utterance_vec = Union{Int64,Missing}[]
    identifier_vec = Union{String,Missing}[]
    path_vec = Union{String,Missing}[]
    pred_vec = Union{String,Missing}[]
    num_tolerance_vec = Union{Int64,Missing}[]
    support_vec = Union{Float64,Missing}[]
    isbest_vec = Union{Bool,Missing}[]
    iscorrect_vec = Union{Bool,Missing}[]
    isnovel_vec = Union{Bool,Missing}[]

    i2f = cue_obj_train.i2f
    for (i, r) in enumerate(res)
        is_best = true
        if length(r) == 0
            utterance = i
            identifier = data[i, target_col]
            path = missing
            pred = missing
            num_tolerance = missing
            support = missing
            is_correct = missing
            is_novel = missing
            is_best = missing
            push!(utterance_vec, utterance)
            push!(identifier_vec, identifier)
            push!(path_vec, path)
            push!(pred_vec, pred)
            push!(num_tolerance_vec, num_tolerance)
            push!(support_vec, support)
            push!(isbest_vec, is_best)
            push!(iscorrect_vec, is_correct)
            push!(isnovel_vec, is_novel)
        end

        for p in r
            utterance = i
            identifier = data[i, target_col]
            path = translate_path(p.ngrams_ind, i2f, sep_token = path_sep_token)
            pred = translate(
                p.ngrams_ind,
                i2f,
                grams,
                tokenized,
                sep_token,
                start_end_token,
                output_sep_token,
            )
            num_tolerance = p.num_tolerance
            support = p.support
            is_correct = iscorrect(cue_obj_val.gold_ind[i], p.ngrams_ind)
            is_novel = isnovel(cue_obj_train.gold_ind, p.ngrams_ind)
            push!(utterance_vec, utterance)
            push!(identifier_vec, identifier)
            push!(path_vec, path)
            push!(pred_vec, pred)
            push!(num_tolerance_vec, num_tolerance)
            push!(support_vec, support)
            push!(isbest_vec, is_best)
            push!(iscorrect_vec, is_correct)
            push!(isnovel_vec, is_novel)
            is_best = false
        end
    end

    df = DataFrame()
    df.utterance = utterance_vec
    df.identifier = identifier_vec
    df.path = path_vec
    df.pred = pred_vec
    df.num_tolerance = num_tolerance_vec
    df.support = support_vec
    df.isbest = isbest_vec
    df.iscorrect = iscorrect_vec
    df.isnovel = isnovel_vec

    df
end

"""
    write2df(gpi::Vector{Gold_Path_Info_Struct})

Write results into a dataframe for the gold paths' information optionally returned by
`learn_paths` and `build_paths`.

# Obligatory Arguments
- `gpi::Vector{Gold_Path_Info_Struct}`: the gold paths' information

# Examples
```julia
# write gold standard paths to df for training data
JudiLing.write2csv(gpi_train)

# write gold standard paths to df for validation data
JudiLing.write2csv(gpi_val)
```
"""
function write2df(gpi::Vector{Gold_Path_Info_Struct})

    utterance_vec = Union{Int64,Missing}[]
    weakest_support_vec = Union{Float64,Missing}[]
    weakest_support_timestep_vec = Union{Int64,Missing}[]
    support_vec = Union{Float64,Missing}[]
    gold_path_vec = Union{Vector,Missing}[]
    timestep_support_vec = Union{Vector,Missing}[]

    for (i, g) in enumerate(gpi)
        ws, wst = findmin(g.ngrams_ind_support)
        push!(utterance_vec, i)
        push!(weakest_support_vec, ws)
        push!(weakest_support_timestep_vec, wst)
        push!(support_vec, g.support)
        push!(gold_path_vec, g.ngrams_ind)
        push!(timestep_support_vec, g.ngrams_ind_support)
    end

    df = DataFrame()
    df.utterance = utterance_vec
    df.weakest_support = weakest_support_vec
    df.weakest_support_timestep = weakest_support_timestep_vec
    df.support = support_vec
    df.gold_path = gold_path_vec
    df.timestep_support = timestep_support_vec

    df
end

"""
    write2df(ts::Threshold_Stat_Struct)

Write results into a dataframe for threshold and tolerance proportion for each timestep.

# Obligatory Arguments
- `ts::Threshold_Stat_Struct`: the threshold and tolerance proportion

# Examples
```julia
JudiLing.write2df(ts)
```
"""
function write2df(ts::Threshold_Stat_Struct)
    n = ts.timestep

    df = DataFrame(vcat([Int64, Float64, Float64], [Float64 for k in 1:2n]),
        vcat([:utterance, :threshold, :tolerance], [Symbol("thr_prop_$k") for k in 1:n], [Symbol("tlr_prop_$k") for k in 1:n]),
        0)

    for i in 1:ts.n_rows
        push!(df, vcat([i, ts.threshold, ts.tolerance], ts.threshold_prop[i,:], ts.tolerance_prop[i,:]))
    end

    df
end

"""
    save_L_matrix(L, filename)

Save lexome matrix into csv file.

# Obligatory Arguments
- `L::L_Matrix_Struct`: the lexome matrix struct
- `filename::String`: the filename/filepath

# Examples
```julia
JudiLing.save_L_matrix(L, joinpath(@__DIR__, "L.csv"))
```
"""
function save_L_matrix(L, filename)

    L_df = DataFrame(L.L)
    insertcols!(L_df, 1, :col_names => L.i2f)
    CSV.write(filename, L_df, quotestrings = true, header = false)
    nothing
end

"""
    load_L_matrix(filename)

Load lexome matrix from csv file.

# Obligatory Arguments
- `filename::String`: the filename/filepath

# Optional Arguments
- `header::Bool=false`: header in csv

# Examples
```julia
L_load = JudiLing.load_L_matrix(joinpath(@__DIR__, "L.csv"))
```
"""
function load_L_matrix(filename; header = false)

    L_df = DataFrame(CSV.File(filename, header = header))
    i2f = L_df[:, 1]
    f2i = Dict(v => i for (i, v) in enumerate(i2f))
    ncol = size(L_df, 2) - 1
    L = Array(select(L_df, Not(1)))

    L_Matrix_Struct(L, f2i, i2f, ncol)
end

"""
    save_S_matrix(S, filename, data, target_col)

Save S matrix into a csv file.

# Obligatory Arguments
- `S::Matrix`: the S matrix
- `filename::String`: the filename/filepath
- `data::DataFrame`: the data
- `target_col::Symbol`: the name of target column

# Optional Arguments
- `sep::Bool=" "`: separator in CSV file

# Examples
```julia
JudiLing.save_S_matrix(S, joinpath(@__DIR__, "S.csv"), latin, :Word)
```
"""
function save_S_matrix(S, filename, data, target_col; sep = " ")

    S_df = DataFrame(S)
    insertcols!(S_df, 1, :col_names => data[:,target_col])
    CSV.write(filename, S_df, quotestrings = true, header = false, delim=sep)
    nothing
end

"""
    load_S_matrix(filename)

Load S matrix from a csv file.

# Obligatory Arguments
- `filename::String`: the filename/filepath

# Optional Arguments
- `header::Bool=false`: header in csv
- `sep::Bool=" "`: separator in CSV file

# Examples
```julia
JudiLing.load_S_matrix(joinpath(@__DIR__, "S.csv"))
```
"""
function load_S_matrix(filename; header = false, sep = " ")

    S_df = DataFrame(CSV.File(filename, header = header, delim = sep))
    words = S_df[:, 1]
    S = Array(select(S_df, Not(1)))

    S, words
end

"""
    write_comprehension_eval(SChat, SC, data, target_col, filename)

Write comprehension evaluation into a CSV file, include target and predicted
ids and indentifiers and their correlations.

# Obligatory Arguments
- `SChat::Matrix`: the Shat/Chat matrix
- `SC::Matrix`: the S/C matrix
- `data::DataFrame`: the data
- `target_col::Symbol`: the name of target column
- `filename::String`: the filename/filepath

# Optional Arguments
- `k`: top k candidates
- `root_dir::String="."`: dir path for project root dir
- `output_dir::String="."`: output dir inside root dir

# Examples
```julia
JudiLing.write_comprehension_eval(Chat, cue_obj.C, latin, :Word, "output.csv",
    k=10, root_dir=@__DIR__, output_dir="out")
```
"""
function write_comprehension_eval(
    SChat,
    SC,
    data,
    target_col,
    filename;
    k = 10,
    root_dir = ".",
    output_dir = "."
    )

    output_path = joinpath(root_dir, output_dir)
    mkpath(output_path)
    total = size(SChat, 1)

    io = open(joinpath(output_path, filename), "w")

    write(
        io,
        "\"target_utterance\",\"target_identifier\",\"predicted_utterance\",\"predicted_identifier\",\"support\"\n",
    )

    rSC = cor(
        convert(Matrix{Float64}, SChat),
        convert(Matrix{Float64}, SC),
        dims = 2,
    )

    for i = 1:total
        p = sortperm(rSC[i, :], rev = true)
        p = p[1:k]
        for j in p
            write(
                io,
                "\"$i\",\"$(data[i,target_col])\",\"$j\",\"$(data[j,target_col])\",\"$(rSC[i,j])\"\n",
            )
        end
    end

    close(io)
end

"""
    write_comprehension_eval(SChat, SC, SC_rest, data, data_rest, target_col, filename)

Write comprehension evaluation into a CSV file for both training and validation
datasets, include target and predicted ids and indentifiers and their
correlations.

# Obligatory Arguments
- `SChat::Matrix`: the Shat/Chat matrix
- `SC::Matrix`: the S/C matrix
- `SC_rest::Matrix`: the rest S/C matrix
- `data::DataFrame`: the data
- `data_rest::DataFrame`: the rest data
- `target_col::Symbol`: the name of target column
- `filename::String`: the filename/filepath

# Optional Arguments
- `k`: top k candidates
- `root_dir::String="."`: dir path for project root dir
- `output_dir::String="."`: output dir inside root dir

# Examples
```julia
JudiLing.write_comprehension_eval(Shat_val, S_val, S_train, latin_val, latin_train,
    :Word, "all_output.csv", k=10, root_dir=@__DIR__, output_dir="out")
```
"""
function write_comprehension_eval(
    SChat,
    SC,
    SC_rest,
    data,
    data_rest,
    target_col,
    filename;
    k = 10,
    root_dir = ".",
    output_dir = "."
    )

    data_combined = copy(data)
    append!(data_combined, data_rest)

    write_comprehension_eval(
        SChat,
        vcat(SC, SC_rest),
        data_combined,
        target_col,
        filename,
        k = k,
        root_dir = root_dir,
        output_dir = output_dir
        )

end


text:
# Generate scenarios of the demands for use within the kNN-SAA framework of Bertsimas and Kallus

# find the indices of the k minimum elements of an array
function mink(arr::Array{Float64}, k_val::Int64)

	v = (sortperm(arr,rev=false))[1:k_val]
	return v
end

# compute distances between covariate samples for kNN
function computeSampCovDist(inputs::Array{Float64,2})

	numSamples::Int64 = size(inputs,1)
	dist_cov = zeros(Float64,numSamples,numSamples)
	for s = 1:numSamples
		for t = s+1:numSamples
			dist_cov[s,t] = norm(inputs[s,:]-inputs[t,:],2)
			dist_cov[t,s] = dist_cov[s,t]
		end
	end

	return dist_cov
end

function generatekNNScenarios(outputs::Array{Float64,2},inputs::Array{Float64,2},numFolds::Int64,minExponent::Float64,maxExponent::Float64,covariate_obs::Array{Float64,1})

	dist_cov = computeSampCovDist(inputs)

	numDependentVar::Int64 = size(outputs,2)
	numSamples::Int64 = size(outputs,1)

	min_error_est::Float64 = Inf
	numRepeatCV::Int64 = 1
	k_sweep = floor.(Int64,numSamples^minExponent):ceil.(Int64,numSamples^maxExponent)
	k_optimal::Int64 = k_sweep[1]
	error_k = zeros(Float64,length(k_sweep))

	# replicates of k-fold CV, if needed
	for rep = 1:numRepeatCV

		# first, randomly permute the given data points for k-fold CV
		rand_perm = randperm(numSamples)

		output_perm = zeros(Float64,numSamples,numDependentVar)
		for d = 1:numDependentVar
			output_perm[:,d] = outputs[rand_perm,d]
		end

		# next, compute distances for kNN
		dist_perm_tmp = zeros(Float64,numSamples,numSamples)
		dist_perm = zeros(Float64,numSamples,numSamples)
		for s = 1:numSamples
			dist_perm_tmp[:,s] = dist_cov[:,rand_perm[s]]
		end
		for s = 1:numSamples
			dist_perm[s,:] = dist_perm_tmp[rand_perm[s],:]
		end

		# then, compute the folds of data for k-fold CV
		numBaseElements::Int64 = floor.(Int64,numSamples*1.0/numFolds)
		numRemElements::Int64 = numSamples - numFolds*numBaseElements
		index_folds = zeros(Int64,numFolds,2)
		index_folds[1,1] = 1
		index_folds[1,2] = numBaseElements
		if(numRemElements > 0)
			index_folds[1,2] += 1
			numRemElements -= 1
		end
		for f = 2:numFolds
			index_folds[f,1] = index_folds[f-1,2] + 1
			index_folds[f,2] = index_folds[f,1] + (numBaseElements-1)
			if(numRemElements > 0)
				index_folds[f,2] += 1
				numRemElements -= 1
			end
		end

		# compute the k nearest neighbors
		ind_dist_sort = zeros(Int64,numSamples,k_sweep[end])
		for f = 1:numFolds
			for s = index_folds[f,1]:index_folds[f,2]
				dist_s_f = dist_perm[s,:]
				dist_s_f[index_folds[f,1]:index_folds[f,2]] = Inf
				ind_dist_sort[s,:] = mink(dist_s_f,k_sweep[end])
			end
		end

		# try different values of the parameter k and pick one that yields the smallest CV prediction error
		for i = 1:length(k_sweep)
			k_val::Int64 = k_sweep[i]
			for s = 1:numSamples
				est_output_s = mean(output_perm[ind_dist_sort[s,1:k_val],:],1)'
				error_k[i] += norm(est_output_s-output_perm[s,:],2)^2
			end
		end

	end

	for i = 1:length(k_sweep)
		if(error_k[i] < min_error_est)
			min_error_est = error_k[i]
			k_optimal = k_sweep[i]
		end
	end

	# finally, use kNN to predict the scenarios for the demands given a new realization of the covariates
	dist_obs = zeros(Float64,numSamples)
	for s = 1:numSamples
		dist_obs[s] = norm(inputs[s,:]-covariate_obs,2)
	end
	min_ind = mink(dist_obs,k_optimal)

	numSAAScenarios::Int64 = size(min_ind,1)
	demand_scen_kNN = zeros(Float64,numSAAScenarios,numDependentVar)
	for s = 1:numSAAScenarios
		demand_scen_kNN[s,:] = outputs[min_ind[s],:]
	end

	return demand_scen_kNN
end

max_stars_repo_name: rohitkannan/DD-SAA

max_issues_repo_name: rohitkannan/DD-SAA

max_forks_repo_name: rohitkannan/DD-SAA


text:
export SeisEvent

mutable struct SeisEvent
  hdr::SeisHdr
  source::SeisSrc
  data::EventTraceData

  SeisEvent(hdr::SeisHdr, source::SeisSrc, data::EventTraceData) = new(hdr, source, data)
end

function SeisEvent(;
                    hdr ::SeisHdr     = SeisHdr(),
                    source::SeisSrc   = SeisSrc(),
                    data::T           = EventTraceData()
                    ) where {T<:GphysData}

  if T != EventTraceData
    return SeisEvent(hdr, source, convert(EventTraceData, data))
  else
    return SeisEvent(hdr, source, data)
  end
end

# =============================================================================
# Methods from Base
isequal(S::SeisEvent, T::SeisEvent) = min(isequal(S.hdr, T.hdr),
                                          isequal(S.source, T.source),
                                          isequal(S.data, T.data))
==(S::SeisEvent, T::SeisEvent) = isequal(S,T)
sizeof(Ev::SeisEvent) = 24 + sizeof(Ev.hdr) + sizeof(Ev.source) + sizeof(Ev.data)

# SeisEvent
write(io::IO, W::SeisEvent) = ( write(io, getfield(W, :hdr));
                                write(io, getfield(W, :source));
                                write(io, getfield(W, :data))
                                )

read(io::IO, ::Type{SeisEvent}) = SeisEvent(read(io, SeisHdr),
                                            read(io, SeisSrc),
                                            read(io, EventTraceData))

summary(V::SeisEvent) = string("Event ", V.hdr.id, ": SeisEvent with ",
  V.data.n, " channel", V.data.n == 1 ? "" : "s")

function show(io::IO, S::SeisEvent)
  println(io, summary(S))
  println(io, "\n(.hdr)")
  show(io, getfield(S, :hdr))
  println(io, "\n(.source)")
  show(io, getfield(S, :source))
  println(io, "\n(.data)")
  println(io, summary(getfield(S, :data)))
  return nothing
end
show(S::SeisEvent) = show(stdout, S)


text:
#---------------------------------------------------------
# # [Groundwater discharge](@id groundwater-discharge)
#---------------------------------------------------------

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/6_groundwater_discharge.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/6_groundwater_discharge.ipynb)

#md # !!! note
#md #     All the AIBECS tutorials and how-to guides are available as Jupyter notebooks.
#md #     You can view them with [nbviewer](https://nbviewer.jupyter.org/)or execute them online with [binder](https://mybinder.org/) by clicking on the badges above!
#md #     (Note that binder can be slow to launch and its memory caps can be a problem when running.)

#md # !!! note
#md #     This tutorial is very similar to the [river discharge notebook](@id river-discharge)

# In this tutorial we will simulate a fictitious radioactive tracer that is injected into the ocean by groundwater discharge.
# The global groundwater discharge dataset for 40,000 coastal watersheds from [*Luijendijk et al.* (2020)](https://www.nature.com/articles/s41467-020-15064-8) is available from within the AIBECS.
# Once "born", our ficitious tracer decays with a parameter timescale $\tau$ as it flows through ocean basins.
#
# The 3D tracer equation is:
#
# $$\left[\frac{\partial}{\partial t} + \nabla \cdot (\boldsymbol{u} + \mathbf{K}\nabla)\right] x = s_\mathsf{gw} - x / \tau$$
#
# where $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \nabla \right]$ represents the ocean circulation transport.
# ([Tracer transport operators are described in the documentation](@ref tracer-transport-operators).)
# The source of the tracer is $s_\mathsf{gw}$, and $x / \tau$ is the decay rate.

# In AIBECS, we must recast this equation in the generic form
#
# $$\left[\frac{\partial}{\partial t} + \mathbf{T}(\boldsymbol{p})\right] \boldsymbol{x} = \boldsymbol{G}(\boldsymbol{x},\boldsymbol{p}).$$

# We start by telling Julia that we want to use the AIBECS and the OCIM2 transport matrix for the ocean circulation.

using AIBECS
grd, T_OCIM2 = OCIM2.load()

# For the radioactive decay, we simply use

function decay(x, p)
    @unpack τ = p
    return x / τ
end

# To build the groundwater sources, we will load the geographic locations and discharge of groundwater (in m³ yr⁻¹) from the *Luijendijk et al.* (2020)) dataset.

gws = GroundWaters.load()

# This is an array of groundwater sources, for which the type `GroundWaterSource{T}` contains the lat–lon coordinates and discharge in m³ yr⁻¹.
# For example, the first element is

gws[1]

# We can check the locations with

using Plots
scatter([x.lon for x in gws], [x.lat for x in gws],
        zcolor=log10.(ustrip.([x.VFR for x in gws] / u"m^3/s")),
        colorbartitle="log₁₀(discharge / (1 m³ s⁻¹))",
        clim=(0,10), fmt=:png)

# We can regrid these into the OCIM2 grid and return the corresponding vector with

gws_grd = ustrip.(u"m^3/s", regrid(gws, grd))

# (We convert years to seconds for AIBECS)

# We can control the global magnitude of the discharge by chosing a concentration of our tracer in groundwater, $C_\mathsf{gw}$, in mol m⁻³.
# For that, we need the volumes of the grid boxes

v = volumevec(grd)

# because the source is given by the mass (discharge × concentration) divided by the volume of each box:

function s_gw(p)
    @unpack C_gw = p
    return @. gws_grd * C_gw / v
end

# We then write the generic $\boldsymbol{G}$ function, which is

G_gw(x,p) = s_gw(p) - decay(x,p)

# ##### Parameters

# We specify some initial values for the parameters and also include units.

import AIBECS: @units, units
import AIBECS: @initial_value, initial_value
@initial_value @units struct GroundWatersParameters{U} <: AbstractParameters{U}
    τ::U    | 20.0 | u"yr"
    C_gw::U |  1.0 | u"mol/m^3"
end

# Finally, thanks to the initial values we provided, we can instantiate the parameter vector succinctly as

p = GroundWatersParameters()

# We build the state function `F`,

F = AIBECSFunction(T_OCIM2, G_gw)

# the steady-state problem `prob`,

nb = sum(iswet(grd))
x = ones(nb) # initial guess
prob = SteadyStateProblem(F, x, p)

# and solve it

sol = solve(prob, CTKAlg()).u * u"mol/m^3"

# Let's now run some visualizations using the plot recipes.
# Taking a horizontal slice of the 3D field at 200m gives

cmap = :viridis
plothorizontalslice(sol, grd, zunit=u"μmol/m^3", depth=200, color=cmap, clim=(0,100))

# and at 500m:

plothorizontalslice(sol, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,25))

# Or we can increase the decay timescale (×10) and decrease the groundwater concentration (÷10) to get a different (more well-mixed) tracer distribution:

p = GroundWatersParameters(τ = 200.0u"yr", C_gw = 0.1u"mol/m^3")
prob = SteadyStateProblem(F, x, p)
sol_τ50 = solve(prob, CTKAlg()).u * u"mol/m^3"
plothorizontalslice(sol_τ50, grd, zunit=u"μmol/m^3", depth=500, color=cmap, clim=(0,25))


text:
if !isdefined( Main, :PostgresQoL )
    include("../src/PostgresQoL.jl")
    using Main.PostgresQoL
end

conn = dbconnection( "mydb", user="juliadev", password="juliadev" )

dbname(conn)
gettables(conn)

DropTable( Var.( conn |> gettables ) ) |> execute(conn)

CreateTable("table1") |>
    addtablecolumn!( "id", "integer", ["PRIMARY KEY"] ) |>
    addtablecolumn!( "var1", "real" ) |>
    addtablecolumn!( "var2", "text" ) |>
    execute(conn)

gettables(conn)
getcolumnnames( conn, "table1" )

CreateTable("weather") |>
    addtablecolumn!( "city", "varchar(80)" ) |>
    addtablecolumn!( "temp_lo", "int" ) |>
    addtablecolumn!( "temp_hi", "int" ) |>
    addtablecolumn!( "prcp", "real" ) |>
    addtablecolumn!( "date", "date" ) |>
    execute(conn)

CreateTable("cities") |>
    addtablecolumn!( "name", "varchar(80)" ) |>
    addtablecolumn!( "location", "point" ) |>
    execute(conn)

InsertInto("weather") |>
    addentries!( ["San Francisco" |> SQLString, 46, 50, .25, "1994-11-27" |> SQLString] ) |>
    execute(conn)

InsertInto("cities") |>
    addentries!( ["San Francisco" |> SQLString, "(-194.0, 53.0)" |> SQLString] ) |>
    execute(conn)

InsertInto("weather", "city", "temp_lo", "temp_hi", "prcp", "date") |>
    addentries!( ["San Francisco" |> SQLString, 43, 57, 0.0, "1994-11-29" |> SQLString] ) |>
    execute(conn)

InsertInto("weather", "date", "city", "temp_hi", "temp_lo") |>
    addentries!( ["'1994-11-29'", "'Hayward'", 54, 37] ) |>
    execute(conn)

Select("weather") |>
    addqueryfield!("*") |>
    execute(conn) |> display

Select("weather") |>
    addqueryfield!("city") |>
    addqueryfield!( "(temp_hi+temp_lo)/2", "temp_avg" ) |>
    addqueryfield!("date") |>
    execute(conn) |> display

Select("weather") |>
    addqueryfield!("*") |>
    setwhereclause!( SQLAnd(
        SQLeq( "city", "San Francisco" |> SQLString ),
        SQLgt( "prcp", 0.0 )
    ) ) |>
    execute(conn) |> display

Select("weather") |>
    addqueryfield!("*") |>
    addorderby!("city") |>
    execute(conn) |> display

Select("weather") |>
    addqueryfield!("*") |>
    addorderby!("city") |>
    addorderby!("temp_lo") |>
    execute(conn) |> display

Select( "weather", distinct=true ) |>
    addqueryfield!("city") |>
    execute(conn) |> display

Select( "weather", distinct=true ) |>
    addorderby!("city") |>
    addqueryfield!("city") |>
    execute(conn) |> display

jt = JoinTable( "weather", "cities" ) |>
    setjoincondition!( SQLeq( "city", "name" ) )
Select(jt) |>
    addqueryfield!("*") |>
    execute(conn) |> display

jt = JoinTable( "weather", "cities" ) |>
    setjoincondition!( SQLeq( "city", "name" ) )
Select(jt) |>
    addqueryfield!("city") |>
    addqueryfield!("temp_lo") |>
    addqueryfield!("temp_hi") |>
    addqueryfield!("prcp") |>
    addqueryfield!("date") |>
    addqueryfield!("location") |>
    execute(conn) |> display

jt = JoinTable( "weather", "cities" ) |>
    setjoincondition!( SQLeq( Var( "weather", "city" ), Var( "cities", "name" ) ) )
Select(jt) |>
    addqueryfield!(Var( "weather", "city" )) |>
    addqueryfield!(Var( "weather", "temp_lo" )) |>
    addqueryfield!(Var( "weather", "temp_hi" )) |>
    addqueryfield!(Var( "weather", "prcp" )) |>
    addqueryfield!(Var( "weather", "date" )) |>
    addqueryfield!(Var( "cities", "location" )) |>
    execute(conn) |> display

jt = JoinTable( "weather", "cities", :left ) |>
    setjoincondition!( SQLeq( Var( "weather", "city" ), Var( "cities", "name" ) ) )
Select(jt) |>
    addqueryfield!("*") |>
    execute(conn) |> display

jt = JoinTable( "weather", "weather", alias="w1", jalias="w2" ) |>
    setjoincondition!( SQLAnd(
        SQLlt( Var( "w1", "temp_lo" ), Var("w2", "temp_lo" ) ),
        SQLgt( Var( "w1", "temp_hi" ), Var("w2", "temp_hi" ) )
    ) )
Select(jt) |>
    addqueryfield!( Var( "w1", "city" ) ) |>
    addqueryfield!( Var( "w1", "temp_lo" ), "low" ) |>
    addqueryfield!( Var( "w1", "temp_hi" ), "high" ) |>
    addqueryfield!( Var( "w2", "city" ) ) |>
    addqueryfield!( Var( "w2", "temp_lo" ), "low" ) |>
    addqueryfield!( Var( "w2", "temp_hi" ), "high" ) |>
    execute(conn) |> display

Select("weather") |>
    addqueryfield!("max(temp_lo)") |>
    execute(conn) |> display

Select("weather") |>
    addqueryfield!("city") |>
    setwhereclause!( SQLeq( "temp_lo",
        string( "(", Select("weather") |> addqueryfield!("max(temp_lo)"), ")" )
    ) ) |>
    execute(conn) |> display

Select("weather") |>
    addqueryfield!("city") |>
    addqueryfield!("max(temp_lo)") |>
    setgroupby!("city") |>
    execute(conn) |> display

Select("weather") |>
    addqueryfield!("city") |>
    addqueryfield!("max(temp_lo)") |>
    setgroupby!("city") |>
    sethavingclause!( SQLlt( "max(temp_lo)", 40 ) ) |>
    execute(conn) |> display

Select("weather") |>
    addqueryfield!("city") |>
    addqueryfield!("max(temp_lo)") |>
    setwhereclause!( SQLlike( "city", "S%" ) ) |>
    setgroupby!("city") |>
    sethavingclause!( SQLlt( "max(temp_lo)", 40 ) ) |>
    execute(conn) |> display

Update("weather") |>
    addcolumnupdate!( "temp_hi", "temp_hi - 2" ) |>
    addcolumnupdate!( "temp_lo", "temp_lo - 2" ) |>
    setwhereclause!( SQLgt( "date", "'1994-11-28'" ) ) |>
    execute(conn)

Select("weather") |>
    addqueryfield!("*") |>
    execute(conn) |> display

DeleteFrom("weather") |>
    setwhereclause!( SQLeq( "city", "Hayward" |> SQLString ) ) |>
    execute(conn)

Select("weather") |>
    addqueryfield!("*") |>
    execute(conn) |> display

# close(conn)


text:
__precompile__(true)
module Pkg3

const DEPOTS = [joinpath(homedir(), ".julia")]
depots() = DEPOTS

const USE_LIBGIT2_FOR_ALL_DOWNLOADS = false
const NUM_CONCURRENT_DOWNLOADS      = 8

iswindows() = @static VERSION < v"0.7-" ? Sys.is_windows() : Sys.iswindows()
isapple()   = @static VERSION < v"0.7-" ? Sys.is_apple()   : Sys.isapple()
islinux()   = @static VERSION < v"0.7-" ? Sys.is_linux()   : Sys.islinux()

# Backport of Equalto
if !isdefined(Base, :EqualTo)
    struct EqualTo{T} <: Function
        x::T
        EqualTo(x::T) where {T} = new{T}(x)
    end
    (f::EqualTo)(y) = isequal(f.x, y)
    const equalto = EqualTo
end

# load snapshotted dependencies
include("../ext/BinaryProvider/src/BinaryProvider.jl")
include("../ext/TOML/src/TOML.jl")
include("../ext/TerminalMenus/src/TerminalMenus.jl")

include("Pkg2/Pkg2.jl")
include("Types.jl")
include("Display.jl")
include("Operations.jl")
include("REPLMode.jl")
include("API.jl")

import .API: add, rm, up, test
const update = up

@enum LoadErrorChoice LOAD_ERROR_QUERY LOAD_ERROR_INSTALL LOAD_ERROR_ERROR

Base.@kwdef mutable struct GlobalSettings
    load_error_choice::LoadErrorChoice = LOAD_ERROR_QUERY # query, install, or error, when not finding package on import
end

GLOBAL_SETTINGS = GlobalSettings()

function __init__()
    push!(empty!(LOAD_PATH), dirname(dirname(@__DIR__)))
    isdefined(Base, :active_repl) && REPLMode.repl_init(Base.active_repl)
end

function Base.julia_cmd(julia::AbstractString)
    cmd = invoke(Base.julia_cmd, Tuple{Any}, julia)
    push!(cmd.exec, "-L$(abspath(@__DIR__, "require.jl"))")
    return cmd
end

if VERSION < v"0.7.0-DEV.2303"
    Base.find_in_path(name::String, wd::Void)   = _find_package(name)
    Base.find_in_path(name::String, wd::String) = _find_package(name)
else
    Base.find_package(name::String) = _find_package(name)
end

function _find_package(name::String)
    isabspath(name) && return name
    base = name
    if endswith(name, ".jl")
        base = name[1:end-3]
    else
        name = string(base, ".jl")
    end
    info = Pkg3.Operations.package_env_info(base, verb = "use")
    info == nothing && @goto find_global
    haskey(info, "uuid") || @goto find_global
    haskey(info, "hash-sha1") || @goto find_global
    uuid = Base.Random.UUID(info["uuid"])
    hash = Pkg3.Types.SHA1(info["hash-sha1"])
    path = Pkg3.Operations.find_installed(uuid, hash)
    ispath(path) && return joinpath(path, "src", name)

    # If we still haven't found the file, look if the package exists in the registry
    # and query the user (if we are interactive) to install it.
    @label find_global
    if isinteractive()
        env = Types.EnvCache()
        pkgspec = [Types.PackageSpec(base)]

        r = Operations.registry_resolve!(env, pkgspec)
        Types.has_uuid(r[1]) || return nothing

        GLOBAL_SETTINGS.load_error_choice == LOAD_ERROR_INSTALL && @goto install
        GLOBAL_SETTINGS.load_error_choice == LOAD_ERROR_ERROR   && return nothing

        choice = TerminalMenus.request("Could not find package \e[1m$(base)\e[22m, do you want to install it?",
                       TerminalMenus.RadioMenu(["yes", "yes (remember)", "no", "no (remember)"]))

        if choice == 3 || choice == 4
            choice == 4 && (GLOBAL_SETTINGS.load_error_choice = LOAD_ERROR_ERROR)
            return nothing
        end

        choice == 2 && (GLOBAL_SETTINGS.load_error_choice = LOAD_ERROR_INSTALL)
        @label install
        Pkg3.Operations.ensure_resolved(env, pkgspec, true)
        Pkg3.Operations.add(env, pkgspec)
        return _find_package(name)
    end
    return nothing
end

end # module


text:
export iauTporv
"""
In the tangent plane projection, given the rectangular coordinates
of a star and its direction cosines, determine the direction
cosines of the tangent point.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   xi,eta   double    rectangular coordinates of star image (Note 2)
   v        double[3] star's direction cosines (Note 3)

Returned:
   v01      double[3] tangent point's direction cosines, Solution 1
   v02      double[3] tangent point's direction cosines, Solution 2

Returned (function value):
               int     number of solutions:
                     0 = no solutions returned (Note 4)
                     1 = only the first solution is useful (Note 5)
                     2 = both solutions are useful (Note 5)

Notes:

   1. The tangent plane projection is also called the "gnomonic
   projection" and the "central projection".

   2. The eta axis points due north in the adopted coordinate system.
   If the direction cosines represent observed (RA,Dec), the tangent
   plane coordinates (xi,eta) are conventionally called the
   "standard coordinates".  If the direction cosines are with
   respect to a right-handed triad, (xi,eta) are also right-handed.
   The units of (xi,eta) are, effectively, radians at the tangent
   point.

   3. The vector v must be of unit length or the result will be wrong.

   4. Cases where there is no solution can arise only near the poles.
   For example, it is clearly impossible for a star at the pole
   itself to have a non-zero xi value, and hence it is meaningless
   to ask where the tangent point would have to be.

   5. Also near the poles, cases can arise where there are two useful
   solutions.  The return value indicates whether the second of the
   two solutions returned is useful;  1 indicates only one useful
   solution, the usual case.

   6. The basis of the algorithm is to solve the spherical triangle
   PSC, where P is the north celestial pole, S is the star and C is
   the tangent point.  Calling the celestial spherical coordinates
   of the star and tangent point (a,b) and (a0,b0) respectively, and
   writing rho^2 = (xi^2+eta^2) and r^2 = (1+rho^2), and
   transforming the vector v into (a,b) in the normal way, side c is
   then (pi/2-b), side p is sqrt(xi^2+eta^2) and side s (to be
   found) is (pi/2-b0), while angle C is given by sin(C) = xi/rho
   and cos(C) = eta/rho;  angle P (to be found) is (a-a0).  After
   solving the spherical triangle, the result (a0,b0) can be
   expressed in vector form as v0.

   7. This function is a member of the following set:

      spherical      vector         solve for

      iauTpxes      iauTpxev         xi,eta
      iauTpsts      iauTpstv          star
      iauTpors    > iauTporv <       origin

References:

   Calabretta M.R. & Greisen, E.W., 2002, "Representations of
   celestial coordinates in FITS", Astron.Astrophys. 395, 1077

   Green, R.M., "Spherical Astronomy", Cambridge University Press,
   1987, Chapter 13.

This revision:   2018 January 2

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauTporv(double xi, double eta, double v[3],
#              double v01[3], double v02[3])

function iauTporv(xi::Real, eta::Real, v::AbstractVector{<:Real})
   v01 = zeros(Float64, 3)
   v02 = zeros(Float64, 3)

   status = ccall((:iauTporv, libsofa_c), Cint,
       (Cdouble, Cdouble, Ref{Cdouble},
       Ptr{Cdouble}, Ptr{Cdouble}),
       convert(Float64, xi),
       convert(Float64, eta),
       convert(Array{Float64, 1}, v),
       v01, v02)

   return status, SVector{3}(v01), SVector{3}(v02)
end