text:
#########################
# Main types for meshes #
#########################
"""
    Node{dim, T}

A `Node` is a point in space.

# Fields
- `x::Vec{dim,T}`: stores the coordinates
"""
struct Node{dim,T}
    x::Vec{dim,T}
end
Node(x::NTuple{dim,T}) where {dim,T} = Node(Vec{dim,T}(x))
getcoordinates(n::Node) = n.x

abstract type AbstractCell{dim,N,M} end
"""
    Cell{dim,N,M} <: AbstractCell{dim,N,M}
A `Cell` is a sub-domain defined by a collection of `Node`s as it's vertices.
However, a `cell` is not defined by the nodes but rather by the global node ids.
The parameter `dim` refers here to the geometrical/ambient dimension, i.e. the dimension of the `nodes` in the grid and **not** the topological dimension of the cell.
A `Cell` has `N` nodes and `M` faces.

# Fields
- `nodes::Ntuple{N,Int}`: N-tuple that stores the node ids. The ordering defines a cell's and its subentities' orientations.
"""
struct Cell{dim,N,M} <: AbstractCell{dim,N,M}
    nodes::NTuple{N,Int}
end
nfaces(c::C) where {C<:AbstractCell} = nfaces(typeof(c))
nfaces(::Type{<:AbstractCell{dim,N,M}}) where {dim,N,M} = M
nnodes(c::C) where {C<:AbstractCell} = nnodes(typeof(c))
nnodes(::Type{<:AbstractCell{dim,N,M}}) where {dim,N,M} = N

# Typealias for commonly used cells
const Line  = Cell{1,2,2}
const Line2D = Cell{2,2,1}
const Line3D = Cell{3,2,0}
const QuadraticLine = Cell{1,3,2}

const Triangle = Cell{2,3,3}
const QuadraticTriangle = Cell{2,6,3}

const Quadrilateral = Cell{2,4,4}
const Quadrilateral3D = Cell{3,4,1}
const QuadraticQuadrilateral = Cell{2,9,4}

const Tetrahedron = Cell{3,4,4}
const QuadraticTetrahedron = Cell{3,10,4}

const Hexahedron = Cell{3,8,6}

"""
A `CellIndex` wraps an Int and corresponds to a cell with that number in the mesh
"""
struct CellIndex
    idx::Int
end

"""
A `FaceIndex` wraps an (Int, Int) and defines a local face by pointing to a (cell, face).
"""
struct FaceIndex <: BoundaryIndex
    idx::Tuple{Int,Int} # cell and side
end

"""
A `EdgeIndex` wraps an (Int, Int) and defines a local edge by pointing to a (cell, edge).
"""
struct EdgeIndex <: BoundaryIndex
    idx::Tuple{Int,Int} # cell and side
end

"""
A `VertexIndex` wraps an (Int, Int) and defines a local vertex by pointing to a (cell, vert).
"""
struct VertexIndex <: BoundaryIndex
    idx::Tuple{Int,Int} # cell and side
end

struct EntityNeighborhood{T<:Union{BoundaryIndex,CellIndex}}
    neighbor_info::Vector{T}
end

EntityNeighborhood(info::T) where T <: BoundaryIndex = EntityNeighborhood([info])
Base.zero(::Type{EntityNeighborhood{T}}) where T = EntityNeighborhood(T[])
Base.zero(::Type{EntityNeighborhood}) = EntityNeighborhood(BoundaryIndex[])
Base.length(n::EntityNeighborhood) = length(n.neighbor_info)
Base.getindex(n::EntityNeighborhood,i) = getindex(n.neighbor_info,i)
Base.firstindex(n::EntityNeighborhood) = 1
Base.lastindex(n::EntityNeighborhood) = length(n.neighbor_info)
Base.:(==)(n1::EntityNeighborhood, n2::EntityNeighborhood) = n1.neighbor_info == n2.neighbor_info
Base.iterate(n::EntityNeighborhood, state=1) = iterate(n.neighbor_info,state)

function Base.:+(n1::EntityNeighborhood, n2::EntityNeighborhood)
    neighbor_info = [n1.neighbor_info; n2.neighbor_info]
    return EntityNeighborhood(neighbor_info)
end

function Base.show(io::IO, ::MIME"text/plain", n::EntityNeighborhood)
    if length(n) == 0
        println(io, "No EntityNeighborhood")
    elseif length(n) == 1
        println(io, "$(n.neighbor_info[1])")
    else
        println(io, "$(n.neighbor_info...)")
    end
end

"""
    face_npoints(::AbstractCell{dim,N,M)
Specifies for each subtype of AbstractCell how many nodes form a face
"""
face_npoints(::Cell{2,N,M}) where {N,M} = 2
face_npoints(::Cell{3,4,1}) = 4 #not sure how to handle embedded cells e.g. Quadrilateral3D
"""
    edge_npoints(::AbstractCell{dim,N,M)
Specifies for each subtype of AbstractCell how many nodes form an edge
"""
edge_npoints(::Cell{3,4,1}) = 2 #not sure how to handle embedded cells e.g. Quadrilateral3D
face_npoints(::Cell{3,N,6}) where N = 4
face_npoints(::Cell{3,N,4}) where N = 3
edge_npoints(::Cell{3,N,M}) where {N,M} = 2

getdim(::Cell{dim}) where dim = dim

abstract type AbstractTopology end

"""
    ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell
`ExclusiveTopology` saves topological (connectivity) data of the grid. The constructor works with an `AbstractCell`
vector for all cells that dispatch `vertices`, `faces` and in 3D `edges` as well as the utility functions
`face_npoints` and `edge_npoints`.
The struct saves the highest dimensional neighborhood, i.e. if something is connected by a face and an
 edge only the face neighborhood is saved. The lower dimensional neighborhood is recomputed, if needed.

# Fields
- `vertex_to_cell::Dict{Int,Vector{Int}}`: global vertex id to all cells containing the vertex
- `cell_neighbor::Vector{EntityNeighborhood{CellIndex}}`: cellid to all connected cells
- `face_neighbor::SparseMatrixCSC{EntityNeighborhood,Int}`: `face_neighbor[cellid,local_face_id]` -> neighboring face
- `vertex_neighbor::SparseMatrixCSC{EntityNeighborhood,Int}`: `vertex_neighbor[cellid,local_vertex_id]` -> neighboring vertex
- `edge_neighbor::SparseMatrixCSC{EntityNeighborhood,Int}`: `edge_neighbor[cellid_local_vertex_id]` -> neighboring edge
- `vertex_vertex_neighbor::Dict{Int,EntityNeighborhood{VertexIndex}}`: global vertex id -> all connected vertices by edge or face
- `face_skeleton::Vector{FaceIndex}`: list of unique faces in the grid
"""
struct ExclusiveTopology <: AbstractTopology
    # maps a global vertex id to all cells containing the vertex
    vertex_to_cell::Dict{Int,Vector{Int}}
    # index of the vector = cell id ->  all other connected cells
    cell_neighbor::Vector{EntityNeighborhood{CellIndex}}
    # face_neighbor[cellid,local_face_id] -> exclusive connected entities (not restricted to one entity)
    face_neighbor::SparseMatrixCSC{EntityNeighborhood,Int}
    # vertex_neighbor[cellid,local_vertex_id] -> exclusive connected entities to the given vertex
    vertex_neighbor::SparseMatrixCSC{EntityNeighborhood,Int}
    # edge_neighbor[cellid,local_edge_id] -> exclusive connected entities of the given edge
    edge_neighbor::SparseMatrixCSC{EntityNeighborhood,Int}
    # maps global vertex id to all directly (by edge or face) connected vertices (no diagonal connection considered)
    vertex_vertex_neighbor::Dict{Int,EntityNeighborhood{VertexIndex}}
    # list of unique faces in the grid given as FaceIndex
    face_skeleton::Vector{FaceIndex}
end

function ExclusiveTopology(cells::Vector{C}) where C <: AbstractCell
    cell_vertices_table = vertices.(cells) #needs generic interface for <: AbstractCell
    vertex_cell_table = Dict{Int,Vector{Int}}()

    for (cellid, cell_nodes) in enumerate(cell_vertices_table)
       for node in cell_nodes
            if haskey(vertex_cell_table, node)
                push!(vertex_cell_table[node], cellid)
            else
                vertex_cell_table[node] = [cellid]
            end
        end
    end

    I_face = Int[]; J_face = Int[]; V_face = EntityNeighborhood[]
    I_edge = Int[]; J_edge = Int[]; V_edge = EntityNeighborhood[]
    I_vertex = Int[]; J_vertex = Int[]; V_vertex = EntityNeighborhood[]
    cell_neighbor_table = Vector{EntityNeighborhood{CellIndex}}(undef, length(cells))

    for (cellid, cell) in enumerate(cells)
        #cell neighborhood
        cell_neighbors = getindex.((vertex_cell_table,), cell_vertices_table[cellid]) # cell -> vertex -> cell
        cell_neighbors = unique(reduce(vcat,cell_neighbors)) # non unique list initially
        filter!(x->x!=cellid, cell_neighbors) # get rid of self neighborhood
        cell_neighbor_table[cellid] = EntityNeighborhood(CellIndex.(cell_neighbors))

        for neighbor in cell_neighbors
            neighbor_local_ids = findall(x->x in cell.nodes, cells[neighbor].nodes)
            cell_local_ids = findall(x->x in cells[neighbor].nodes, cell.nodes)
            # vertex neighbor
            if length(cell_local_ids) == 1
                _vertex_neighbor!(V_vertex, I_vertex, J_vertex, cellid, cell, neighbor_local_ids, neighbor, cells[neighbor])
            # face neighbor
            elseif length(cell_local_ids) == face_npoints(cell)
                _face_neighbor!(V_face, I_face, J_face, cellid, cell, neighbor_local_ids, neighbor, cells[neighbor])
            # edge neighbor
            elseif getdim(cell) > 2 && length(cell_local_ids) == edge_npoints(cell)
                _edge_neighbor!(V_edge, I_edge, J_edge, cellid, cell, neighbor_local_ids, neighbor, cells[neighbor])
            end
        end
    end

    face_neighbor = sparse(I_face,J_face,V_face)
    vertex_neighbor = sparse(I_vertex,J_vertex,V_vertex)
    edge_neighbor = sparse(I_edge,J_edge,V_edge)

    vertex_vertex_table = Dict{Int,EntityNeighborhood}()
    vertex_vertex_global = Dict{Int,Vector{Int}}()
    # Vertex Connectivity
    for global_vertexid in keys(vertex_cell_table)
        #Cellset that contains given vertex
        cellset = vertex_cell_table[global_vertexid]
        vertex_neighbors_local = VertexIndex[]
        vertex_neighbors_global = Int[]
        for cell in cellset
            neighbor_boundary = getdim(cells[cell]) == 2 ? [faces(cells[cell])...] : [edges(cells[cell])...] #get lowest dimension boundary
            neighbor_connected_faces = neighbor_boundary[findall(x->global_vertexid in x, neighbor_boundary)]
            neighbor_vertices_global = getindex.(neighbor_connected_faces, findfirst.(x->x!=global_vertexid,neighbor_connected_faces))
            neighbor_vertices_local= [VertexIndex(cell,local_vertex) for local_vertex in findall(x->x in neighbor_vertices_global, vertices(cells[cell]))]
            append!(vertex_neighbors_local, neighbor_vertices_local)
            append!(vertex_neighbors_global, neighbor_vertices_global)
        end
        vertex_vertex_table[global_vertexid] =  EntityNeighborhood(vertex_neighbors_local)
        vertex_vertex_global[global_vertexid] = vertex_neighbors_global
    end

    # Face Skeleton
    face_skeleton_global = Set{NTuple}()
    face_skeleton_local = Vector{FaceIndex}()
    fs_length = length(face_skeleton_global)
    for (cellid,cell) in enumerate(cells)
        for (local_face_id,face) in enumerate(faces(cell))
            push!(face_skeleton_global, sortface(face))
            fs_length_new = length(face_skeleton_global)
            if fs_length != fs_length_new
                push!(face_skeleton_local, FaceIndex(cellid,local_face_id))
                fs_length = fs_length_new
            end
        end
    end
    return ExclusiveTopology(vertex_cell_table,cell_neighbor_table,face_neighbor,vertex_neighbor,edge_neighbor,vertex_vertex_table,face_skeleton_local)
end

function _vertex_neighbor!(V_vertex, I_vertex, J_vertex, cellid, cell, neighbor, neighborid, neighbor_cell)
    vertex_neighbor = VertexIndex((neighborid, neighbor[1]))
    cell_vertex_id = findfirst(x->x==neighbor_cell.nodes[neighbor[1]], cell.nodes)
    push!(V_vertex,EntityNeighborhood(vertex_neighbor))
    push!(I_vertex,cellid)
    push!(J_vertex,cell_vertex_id)
end

function _edge_neighbor!(V_edge, I_edge, J_edge, cellid, cell, neighbor, neighborid, neighbor_cell)
    neighbor_edge = neighbor_cell.nodes[neighbor]
    if getdim(neighbor_cell) < 3
        neighbor_edge_id = findfirst(x->issubset(x,neighbor_edge), faces(neighbor_cell))
        edge_neighbor = FaceIndex((neighborid, neighbor_edge_id))
    else
        neighbor_edge_id = findfirst(x->issubset(x,neighbor_edge), edges(neighbor_cell))
        edge_neighbor = EdgeIndex((neighborid, neighbor_edge_id))
    end
    cell_edge_id = findfirst(x->issubset(x,neighbor_edge),edges(cell))
    push!(V_edge, EntityNeighborhood(edge_neighbor))
    push!(I_edge, cellid)
    push!(J_edge, cell_edge_id)
end

function _face_neighbor!(V_face, I_face, J_face, cellid, cell, neighbor, neighborid, neighbor_cell)
    neighbor_face = neighbor_cell.nodes[neighbor]
    if getdim(neighbor_cell) == getdim(cell)
        neighbor_face_id = findfirst(x->issubset(x,neighbor_face), faces(neighbor_cell))
        face_neighbor = FaceIndex((neighborid, neighbor_face_id))
    else
        neighbor_face_id = findfirst(x->issubset(x,neighbor_face), edges(neighbor_cell))
        face_neighbor = EdgeIndex((neighborid, neighbor_face_id))
    end
    cell_face_id = findfirst(x->issubset(x,neighbor_face),faces(cell))
    push!(V_face, EntityNeighborhood(face_neighbor))
    push!(I_face, cellid)
    push!(J_face, cell_face_id)
end

getcells(neighbor::EntityNeighborhood{T}) where T <: BoundaryIndex = first.(neighbor.neighbor_info)
getcells(neighbor::EntityNeighborhood{CellIndex}) = getproperty.(neighbor.neighbor_info, :idx)
getcells(neighbors::Vector{T}) where T <: EntityNeighborhood = reduce(vcat, getcells.(neighbors))

abstract type AbstractGrid{dim} end

ExclusiveTopology(grid::AbstractGrid) = ExclusiveTopology(getcells(grid))

"""
    Grid{dim, C<:AbstractCell, T<:Real} <: AbstractGrid}

A `Grid` is a collection of `Cells` and `Node`s which covers the computational domain, together with Sets of cells, nodes and faces.
There are multiple helper structures to apply boundary conditions or define subdomains. They are gathered in the `cellsets`, `nodesets`,
`facesets`, `edgesets` and `vertexsets`.

# Fields
- `cells::Vector{C}`: stores all cells of the grid
- `nodes::Vector{Node{dim,T}}`: stores the `dim` dimensional nodes of the grid
- `cellsets::Dict{String,Set{Int}}`: maps a `String` key to a `Set` of cell ids
- `nodesets::Dict{String,Set{Int}}`: maps a `String` key to a `Set` of global node ids
- `facesets::Dict{String,Set{FaceIndex}}`: maps a `String` to a `Set` of `Set{FaceIndex} (global_cell_id, local_face_id)`
- `edgesets::Dict{String,Set{EdgeIndex}}`: maps a `String` to a `Set` of `Set{EdgeIndex} (global_cell_id, local_edge_id`
- `vertexsets::Dict{String,Set{VertexIndex}}`: maps a `String` key to a `Set` of local vertex ids
- `boundary_matrix::SparseMatrixCSC{Bool,Int}`: optional, only needed by `onboundary` to check if a cell is on the boundary, see, e.g. Helmholtz example
"""
mutable struct Grid{dim,C<:AbstractCell,T<:Real} <: AbstractGrid{dim}
    cells::Vector{C}
    nodes::Vector{Node{dim,T}}
    # Sets
    cellsets::Dict{String,Set{Int}}
    nodesets::Dict{String,Set{Int}}
    facesets::Dict{String,Set{FaceIndex}}
    edgesets::Dict{String,Set{EdgeIndex}}
    vertexsets::Dict{String,Set{VertexIndex}}
    # Boundary matrix (faces per cell × cell)
    boundary_matrix::SparseMatrixCSC{Bool,Int}
end

function Grid(cells::Vector{C},
              nodes::Vector{Node{dim,T}};
              cellsets::Dict{String,Set{Int}}=Dict{String,Set{Int}}(),
              nodesets::Dict{String,Set{Int}}=Dict{String,Set{Int}}(),
              facesets::Dict{String,Set{FaceIndex}}=Dict{String,Set{FaceIndex}}(),
              edgesets::Dict{String,Set{EdgeIndex}}=Dict{String,Set{EdgeIndex}}(),
              vertexsets::Dict{String,Set{VertexIndex}}=Dict{String,Set{VertexIndex}}(),
              boundary_matrix::SparseMatrixCSC{Bool,Int}=spzeros(Bool, 0, 0)) where {dim,C,T}
    return Grid(cells, nodes, cellsets, nodesets, facesets, edgesets, vertexsets, boundary_matrix)
end

##########################
# Grid utility functions #
##########################
"""
    getneighborhood(top::ExclusiveTopology, grid::Grid{dim,C,T}, cellidx::CellIndex, include_self=false)
    getneighborhood(top::ExclusiveTopology, grid::Grid{dim,C,T}, faceidx::FaceIndex, include_self=false)
    getneighborhood(top::ExclusiveTopology, grid::Grid{dim,C,T}, vertexidx::VertexIndex, include_self=false)
    getneighborhood(top::ExclusiveTopology, grid::Grid{dim,C,T}, edgeidx::EdgeIndex, include_self=false)

Returns all directly connected entities of the same type, i.e. calling the function with a `VertexIndex` will return
a list of directly connected vertices (connected via face/edge). If `include_self` is true, the given `*Index` is included
in the returned list.

!!! warning
    This feature is highly experimental and very likely subjected to interface changes in the future.
"""
function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, cellidx::CellIndex, include_self=false)
    patch = getcells(top.cell_neighbor[cellidx.idx])
    if include_self
        return [patch; cellidx.idx]
    else
        return patch
    end
end

function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, faceidx::FaceIndex, include_self=false)
    if include_self
        return [top.face_neighbor[faceidx[1],faceidx[2]].neighbor_info; faceidx]
    else
        return top.face_neighbor[faceidx[1],faceidx[2]].neighbor_info
    end
end

function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, vertexidx::VertexIndex, include_self=false)
    cellid, local_vertexid = vertexidx[1], vertexidx[2]
    cell_vertices = vertices(getcells(grid,cellid))
    global_vertexid = cell_vertices[local_vertexid]
    if include_self
        return [top.vertex_vertex_neighbor[global_vertexid].neighbor_info; vertexidx]
    else
        return top.vertex_vertex_neighbor[global_vertexid].neighbor_info
    end
end

function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid{3}, edgeidx::EdgeIndex, include_self=false)
    if include_self
        return [top.edge_neighbor[edgeidx[1],edgeidx[2]].neighbor_info; edgeidx]
    else
        return top.edge_neighbor[edgeidx[1],edgeidx[2]].neighbor_info
    end
end

"""
    faceskeleton(grid) -> Vector{FaceIndex}
Returns an iterateable face skeleton. The skeleton consists of `FaceIndex` that can be used to `reinit`
`FaceValues`.
"""
faceskeleton(top::ExclusiveTopology, grid::AbstractGrid) =  top.face_skeleton

toglobal(grid::Grid,vertexidx::VertexIndex) = vertices(getcells(grid,vertexidx[1]))[vertexidx[2]]
toglobal(grid::Grid,vertexidx::Vector{VertexIndex}) = unique(toglobal.((grid,),vertexidx))

@inline getdim(::AbstractGrid{dim}) where {dim} = dim
"""
    getcells(grid::AbstractGrid)
    getcells(grid::AbstractGrid, v::Union{Int,Vector{Int}}
    getcells(grid::AbstractGrid, setname::String)

Returns either all `cells::Collection{C<:AbstractCell}` of a `<:AbstractGrid` or a subset based on an `Int`, `Vector{Int}` or `String`.
Whereas the last option tries to call a `cellset` of the `grid`. `Collection` can be any indexable type, for `Grid` it is `Vector{C<:AbstractCell}`.
"""
@inline getcells(grid::AbstractGrid) = grid.cells
@inline getcells(grid::AbstractGrid, v::Union{Int, Vector{Int}}) = grid.cells[v]
@inline getcells(grid::AbstractGrid, setname::String) = grid.cells[collect(getcellset(grid,setname))]
"Returns the number of cells in the `<:AbstractGrid`."
@inline getncells(grid::AbstractGrid) = length(grid.cells)
"Returns the celltype of the `<:AbstractGrid`."
@inline getcelltype(grid::AbstractGrid) = eltype(grid.cells)
@inline getcelltype(grid::AbstractGrid, i::Int) = typeof(grid.cells[i])

"""
    getnodes(grid::AbstractGrid)
    getnodes(grid::AbstractGrid, v::Union{Int,Vector{Int}}
    getnodes(grid::AbstractGrid, setname::String)

Returns either all `nodes::Collection{N}` of a `<:AbstractGrid` or a subset based on an `Int`, `Vector{Int}` or `String`.
The last option tries to call a `nodeset` of the `<:AbstractGrid`. `Collection{N}` refers to some indexable collection where each element corresponds
to a Node.
"""
@inline getnodes(grid::AbstractGrid) = grid.nodes
@inline getnodes(grid::AbstractGrid, v::Union{Int, Vector{Int}}) = grid.nodes[v]
@inline getnodes(grid::AbstractGrid, setname::String) = grid.nodes[collect(getnodeset(grid,setname))]
"Returns the number of nodes in the grid."
@inline getnnodes(grid::AbstractGrid) = length(grid.nodes)
"Returns the number of nodes of the `i`-th cell."
@inline nnodes_per_cell(grid::AbstractGrid, i::Int=1) = nnodes(grid.cells[i])

"""
    getcellset(grid::AbstractGrid, setname::String)

Returns all cells as cellid in a `Set` of a given `setname`.
"""
@inline getcellset(grid::AbstractGrid, setname::String) = grid.cellsets[setname]
"""
    getcellsets(grid::AbstractGrid)

Returns all cellsets of the `grid`.
"""
@inline getcellsets(grid::AbstractGrid) = grid.cellsets

"""
    getnodeset(grid::AbstractGrid, setname::String)

Returns all nodes as nodeid in a `Set` of a given `setname`.
"""
@inline getnodeset(grid::AbstractGrid, setname::String) = grid.nodesets[setname]
"""
    getnodesets(grid::AbstractGrid)

Returns all nodesets of the `grid`.
"""
@inline getnodesets(grid::AbstractGrid) = grid.nodesets

"""
    getfaceset(grid::AbstractGrid, setname::String)

Returns all faces as `FaceIndex` in a `Set` of a given `setname`.
"""
@inline getfaceset(grid::AbstractGrid, setname::String) = grid.facesets[setname]
"""
    getfacesets(grid::AbstractGrid)

Returns all facesets of the `grid`.
"""
@inline getfacesets(grid::AbstractGrid) = grid.facesets

"""
    getedgeset(grid::AbstractGrid, setname::String)

Returns all edges as `EdgeIndex` in a `Set` of a given `setname`.
"""
@inline getedgeset(grid::AbstractGrid, setname::String) = grid.edgesets[setname]
"""
    getedgesets(grid::AbstractGrid)

Returns all edge sets of the grid.
"""
@inline getedgesets(grid::AbstractGrid) = grid.edgesets

"""
    getedgeset(grid::AbstractGrid, setname::String)

Returns all vertices as `VertexIndex` in a `Set` of a given `setname`.
"""
@inline getvertexset(grid::AbstractGrid, setname::String) = grid.vertexsets[setname]
"""
    getvertexsets(grid::AbstractGrid)

Returns all vertex sets of the grid.
"""
@inline getvertexsets(grid::AbstractGrid) = grid.vertexsets

n_faces_per_cell(grid::Grid) = nfaces(eltype(grid.cells))

"""
    function compute_vertex_values(grid::AbstractGrid, f::Function)
    function compute_vertex_values(grid::AbstractGrid, v::Vector{Int}, f::Function)
    function compute_vertex_values(grid::AbstractGrid, set::String, f::Function)

Given a `grid` and some function `f`, `compute_vertex_values` computes all nodal values,
 i.e. values at the nodes,  of the function `f`.
The function implements two dispatches, where only a subset of the grid's node is used.

```julia
    compute_vertex_values(grid, x -> sin(x[1]) + cos([2]))
    compute_vertex_values(grid, [9, 6, 3], x -> sin(x[1]) + cos([2])) #compute function values at nodes with id 9,6,3
    compute_vertex_values(grid, "right", x -> sin(x[1]) + cos([2])) #compute function values at nodes belonging to nodeset right
```

"""
@inline function compute_vertex_values(nodes::Vector{Node{dim,T}}, f::Function) where{dim,T}
    map(n -> f(getcoordinates(n)), nodes)
end

@inline function compute_vertex_values(grid::AbstractGrid, f::Function)
    compute_vertex_values(getnodes(grid), f::Function)
end

@inline function compute_vertex_values(grid::AbstractGrid, v::Vector{Int}, f::Function)
    compute_vertex_values(getnodes(grid, v), f::Function)
end

@inline function compute_vertex_values(grid::AbstractGrid, set::String, f::Function)
    compute_vertex_values(getnodes(grid, set), f::Function)
end

# Transformations
"""
    transform!(grid::Abstractgrid, f::Function)

Transform all nodes of the `grid` based on some transformation function `f`.
"""
function transform!(g::AbstractGrid, f::Function)
    c = similar(g.nodes)
    for i in 1:length(c)
        c[i] = Node(f(g.nodes[i].x))
    end
    copyto!(g.nodes, c)
    g
end

# Sets

_check_setname(dict, name) = haskey(dict, name) && throw(ArgumentError("there already exists a set with the name: $name"))
_warn_emptyset(set, name) = length(set) == 0 && @warn("no entities added to the set with name: $name")

"""
    addcellset!(grid::AbstractGrid, name::String, cellid::Union{Set{Int}, Vector{Int}})
    addcellset!(grid::AbstractGrid, name::String, f::function; all::Bool=true)

Adds a cellset to the grid with key `name`.
Cellsets are typically used to define subdomains of the problem, e.g. two materials in the computational domain.
The `MixedDofHandler` can construct different fields which live not on the whole domain, but rather on a cellset.

```julia
addcellset!(grid, "left", Set((1,3))) #add cells with id 1 and 3 to cellset left
addcellset!(grid, "right", x -> norm(x[1]) < 2.0 ) #add cell to cellset right, if x[1] of each cell's node is smaller than 2.0
```
"""
function addcellset!(grid::AbstractGrid, name::String, cellid::Union{Set{Int},Vector{Int}})
    _check_setname(grid.cellsets,  name)
    cells = Set(cellid)
    _warn_emptyset(cells, name)
    grid.cellsets[name] = cells
    grid
end

function addcellset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true)
    _check_setname(grid.cellsets, name)
    cells = Set{Int}()
    for (i, cell) in enumerate(getcells(grid))
        pass = all
        for node_idx in cell.nodes
            node = grid.nodes[node_idx]
            v = f(node.x)
            all ? (!v && (pass = false; break)) : (v && (pass = true; break))
        end
        pass && push!(cells, i)
    end
    _warn_emptyset(cells, name)
    grid.cellsets[name] = cells
    grid
end

"""
    addfaceset!(grid::AbstractGrid, name::String, faceid::Union{Set{FaceIndex},Vector{FaceIndex}})
    addfaceset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true)

Adds a faceset to the grid with key `name`.
A faceset maps a `String` key to a `Set` of tuples corresponding to `(global_cell_id, local_face_id)`.
Facesets are used to initialize `Dirichlet` structs, that are needed to specify the boundary for the `ConstraintHandler`.

```julia
addfaceset!(gird, "right", Set(((2,2),(4,2))) #see grid manual example for reference
addfaceset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0) #see incompressible elasticity example for reference
```
"""
addfaceset!(grid::Grid, name::String, set::Union{Set{FaceIndex},Vector{FaceIndex}}) =
    _addset!(grid, name, set, grid.facesets)
addedgeset!(grid::Grid, name::String, set::Union{Set{EdgeIndex},Vector{EdgeIndex}}) =
    _addset!(grid, name, set, grid.edgesets)
addvertexset!(grid::Grid, name::String, set::Union{Set{VertexIndex},Vector{VertexIndex}}) =
    _addset!(grid, name, set, grid.vertexsets)
function _addset!(grid::AbstractGrid, name::String, _set, dict::Dict)
    _check_setname(dict, name)
    set = Set(_set)
    _warn_emptyset(set, name)
    dict[name] = set
    grid
end

addfaceset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true) =
    _addset!(grid, name, f, Ferrite.faces, grid.facesets, FaceIndex; all=all)
addedgeset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true) =
    _addset!(grid, name, f, Ferrite.edges, grid.edgesets, EdgeIndex; all=all)
addvertexset!(grid::AbstractGrid, name::String, f::Function; all::Bool=true) =
    _addset!(grid, name, f, Ferrite.vertices, grid.vertexsets, VertexIndex; all=all)
function _addset!(grid::AbstractGrid, name::String, f::Function, _ftype::Function, dict::Dict, _indextype::Type; all::Bool=true)
    _check_setname(dict, name)
    _set = Set{_indextype}()
    for (cell_idx, cell) in enumerate(getcells(grid))
        for (face_idx, face) in enumerate(_ftype(cell))
            pass = all
            for node_idx in face
                v = f(grid.nodes[node_idx].x)
                all ? (!v && (pass = false; break)) : (v && (pass = true; break))
            end
            pass && push!(_set, _indextype(cell_idx, face_idx))
        end
    end
    _warn_emptyset(_set, name)
    dict[name] = _set
    grid
end

"""
    addnodeset!(grid::AbstractGrid, name::String, nodeid::Union{Vector{Int},Set{Int}})
    addnodeset!(grid::AbstractGrid, name::String, f::Function)

Adds a `nodeset::Dict{String, Set{Int}}` to the `grid` with key `name`. Has the same interface as `addcellset`.
However, instead of mapping a cell id to the `String` key, a set of node ids is returned.
"""
function addnodeset!(grid::AbstractGrid, name::String, nodeid::Union{Vector{Int},Set{Int}})
    _check_setname(grid.nodesets, name)
    grid.nodesets[name] = Set(nodeid)
    _warn_emptyset(grid.nodesets[name], name)
    grid
end

function addnodeset!(grid::AbstractGrid, name::String, f::Function)
    _check_setname(grid.nodesets, name)
    nodes = Set{Int}()
    for (i, n) in enumerate(getnodes(grid))
        f(n.x) && push!(nodes, i)
    end
    grid.nodesets[name] = nodes
    _warn_emptyset(grid.nodesets[name], name)
    grid
end

"""
    getcoordinates!(x::Vector{Vec{dim,T}}, grid::AbstractGrid, cell::Int)
Fills the vector `x` with the coordinates of a cell, defined by its cell id.
"""
@inline function getcoordinates!(x::Vector{Vec{dim,T}}, grid::AbstractGrid, cell::Int) where {dim,T}
    #@assert length(x) == N
    @inbounds for i in 1:length(x)
        x[i] = grid.nodes[grid.cells[cell].nodes[i]].x
    end
end
@inline getcoordinates!(x::Vector{Vec{dim,T}}, grid::AbstractGrid, cell::CellIndex) where {dim, T} = getcoordinates!(x, grid, cell.idx)
@inline getcoordinates!(x::Vector{Vec{dim,T}}, grid::AbstractGrid, face::FaceIndex) where {dim, T} = getcoordinates!(x, grid, face.idx[1])

"""
    getcoordinates(grid::AbstractGrid, cell)
Return a vector with the coordinates of the vertices of cell number `cell`.
"""
@inline function getcoordinates(grid::AbstractGrid, cell::Int)
    # TODO pretty ugly, worth it?
    dim = typeof(grid.cells[cell]).parameters[1]
    T = typeof(grid).parameters[3]
    nodeidx = grid.cells[cell].nodes
    return [grid.nodes[i].x for i in nodeidx]::Vector{Vec{dim,T}}
end
@inline getcoordinates(grid::AbstractGrid, cell::CellIndex) = getcoordinates(grid, cell.idx)
@inline getcoordinates(grid::AbstractGrid, face::FaceIndex) = getcoordinates(grid, face.idx[1])

function Base.show(io::IO, ::MIME"text/plain", grid::Grid)
    print(io, "$(typeof(grid)) with $(getncells(grid)) ")
    typestrs = sort!(collect(Set(celltypes[typeof(x)] for x in grid.cells)))
    str = join(io, typestrs, '/')
    print(io, " cells and $(getnnodes(grid)) nodes")
end

const celltypes = Dict{DataType, String}(Cell{1,2,2}  => "Line",
                                         Cell{2,2,2}  => "2D-Line",
                                         Cell{3,2,0}  => "3D-Line",
                                         Cell{1,3,2}  => "QuadraticLine",
                                         Cell{2,3,3}  => "Triangle",
                                         Cell{2,6,3}  => "QuadraticTriangle",
                                         Cell{2,4,4}  => "Quadrilateral",
                                         Cell{3,4,1}  => "3D-Quadrilateral",
                                         Cell{2,9,4}  => "QuadraticQuadrilateral",
                                         Cell{3,4,4}  => "Tetrahedron",
                                         Cell{3,10,4} => "QuadraticTetrahedron",
                                         Cell{3,8,6}  => "Hexahedron",
                                         Cell{3,20,6} => "Cell{3,20,6}")

# Functions to uniquely identify vertices, edges and faces, used when distributing
# dofs over a mesh. For this we can ignore the nodes on edged, faces and inside cells,
# we only need to use the nodes that are vertices.
# 1D: vertices
faces(c::Union{Line,QuadraticLine}) = (c.nodes[1], c.nodes[2])
vertices(c::Union{Line,Line2D,Line3D,QuadraticLine}) = (c.nodes[1], c.nodes[2])
# 2D: vertices, faces
faces(c::Line2D) = ((c.nodes[1],c.nodes[2]),)
vertices(c::Union{Triangle,QuadraticTriangle}) = (c.nodes[1], c.nodes[2], c.nodes[3])
faces(c::Union{Triangle,QuadraticTriangle}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[1]))
vertices(c::Union{Quadrilateral,Quadrilateral3D,QuadraticQuadrilateral}) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4])
faces(c::Union{Quadrilateral,QuadraticQuadrilateral}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[4]), (c.nodes[4],c.nodes[1]))
# 3D: vertices, edges, faces
edges(c::Line3D) = ((c.nodes[1],c.nodes[2]),)
vertices(c::Union{Tetrahedron,QuadraticTetrahedron}) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4])
edges(c::Union{Tetrahedron,QuadraticTetrahedron}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[1]), (c.nodes[1],c.nodes[4]), (c.nodes[2],c.nodes[4]), (c.nodes[3],c.nodes[4]))
faces(c::Union{Tetrahedron,QuadraticTetrahedron}) = ((c.nodes[1],c.nodes[2],c.nodes[3]), (c.nodes[1],c.nodes[2],c.nodes[4]), (c.nodes[2],c.nodes[3],c.nodes[4]), (c.nodes[1],c.nodes[4],c.nodes[3]))
vertices(c::Union{Hexahedron,Cell{3,20,6}}) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4], c.nodes[5], c.nodes[6], c.nodes[7], c.nodes[8])
edges(c::Union{Hexahedron,Cell{3,20,6}}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[4]), (c.nodes[4],c.nodes[1]), (c.nodes[5],c.nodes[6]), (c.nodes[6],c.nodes[7]), (c.nodes[7],c.nodes[8]), (c.nodes[8],c.nodes[5]), (c.nodes[1],c.nodes[5]), (c.nodes[2],c.nodes[6]), (c.nodes[3],c.nodes[7]), (c.nodes[4],c.nodes[8]))
faces(c::Union{Hexahedron,Cell{3,20,6}}) = ((c.nodes[1],c.nodes[4],c.nodes[3],c.nodes[2]), (c.nodes[1],c.nodes[2],c.nodes[6],c.nodes[5]), (c.nodes[2],c.nodes[3],c.nodes[7],c.nodes[6]), (c.nodes[3],c.nodes[4],c.nodes[8],c.nodes[7]), (c.nodes[1],c.nodes[5],c.nodes[8],c.nodes[4]), (c.nodes[5],c.nodes[6],c.nodes[7],c.nodes[8]))
edges(c::Union{Quadrilateral3D}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[4]), (c.nodes[4],c.nodes[1]))
faces(c::Union{Quadrilateral3D}) = ((c.nodes[1],c.nodes[2],c.nodes[3],c.nodes[4]),)

# random stuff
default_interpolation(::Union{Type{Line},Type{Line2D},Type{Line3D}}) = Lagrange{1,RefCube,1}()
default_interpolation(::Type{QuadraticLine}) = Lagrange{1,RefCube,2}()
default_interpolation(::Type{Triangle}) = Lagrange{2,RefTetrahedron,1}()
default_interpolation(::Type{QuadraticTriangle}) = Lagrange{2,RefTetrahedron,2}()
default_interpolation(::Union{Type{Quadrilateral},Type{Quadrilateral3D}}) = Lagrange{2,RefCube,1}()
default_interpolation(::Type{QuadraticQuadrilateral}) = Lagrange{2,RefCube,2}()
default_interpolation(::Type{Tetrahedron}) = Lagrange{3,RefTetrahedron,1}()
default_interpolation(::Type{QuadraticTetrahedron}) = Lagrange{3,RefTetrahedron,2}()
default_interpolation(::Type{Hexahedron}) = Lagrange{3,RefCube,1}()
default_interpolation(::Type{Cell{3,20,6}}) = Serendipity{3,RefCube,2}()

boundaryfunction(::Type{FaceIndex}) = Ferrite.faces
boundaryfunction(::Type{EdgeIndex}) = Ferrite.edges
boundaryfunction(::Type{VertexIndex}) = Ferrite.vertices

for INDEX in (:VertexIndex, :EdgeIndex, :FaceIndex)
    @eval begin
        #Constructor
        ($INDEX)(a::Int, b::Int) = ($INDEX)((a,b))

        Base.getindex(I::($INDEX), i::Int) = I.idx[i]

        #To be able to do a,b = faceidx
        Base.iterate(I::($INDEX), state::Int=1) = (state==3) ?  nothing : (I[state], state+1)

        #For (cellid, faceidx) in faceset
        Base.in(v::Tuple{Int, Int}, s::Set{$INDEX}) = in($INDEX(v), s)
    end
end


text:
using Flux.Data.MNIST
using Images: channelview, colorview, RGB
using Metalhead

using EmacsREPL: viewrepl

function test()
    # MNIST
    imgs = MNIST.images();
    size(imgs[1])               # 28,28
    typeof(imgs[1][1])          # ColorTypes.Gray
    viewrepl(imgs[1])

    # Converting to float
    size(float.(imgs[1]))       # (28,28)
    typeof(float.(imgs[1])[1])     # Float64
    viewrepl(float.(imgs[1]))
    size(float.(imgs[1])[:,:,:]) # 28,28,1
    viewrepl(float.(imgs[1])[:,:,:])

    # CIFAR
    imgs = trainimgs(CIFAR10);
    size(imgs[1].img)           # 32,32
    typeof(imgs[1].img[1])      # RGB

    viewrepl(imgs[1].img)

    getarray(X) = Float32.(permutedims(channelview(X), (2, 3, 1)))
    getRGB(X) = colorview(RGB, permutedims(X, (3,1,2)))

    # converting to float and channel view
    floatimg = getarray(imgs[1].img);
    size(floatimg)              # 32,32,3
    typeof(floatimg[1])         # Float32
    viewrepl(floatimg)          # this is not working

    # converting back to RGB
    rgbimg = getRGB(floatimg);
    viewrepl(rgbimg)
end

ext: jl

max_stars_repo_name: lihebi/julia-repl

max_issues_repo_name: lihebi/julia-repl

max_forks_repo_name: lihebi/julia-repl


text:
# This file is part of RLinearAlgebra.jl

module ProceduralTestLSLogOracle

using Test, RLinearAlgebra, Random, LinearAlgebra

Random.seed!(1010)

@testset "LS Oracle Log -- Procedural" begin

    # Verify appropriate super type
    @test supertype(LSLogOracle) == LinSysSolverLog

    # Verify Required Fields
    @test :iterations in fieldnames(LSLogOracle)
    @test :converged in fieldnames(LSLogOracle)

    # Verify log_update initialization
    let
        A = rand(2,2)
        x = rand(2)
        b = A * x
        z = rand(2)

        sampler = LinSysVecRowOneRandCyclic()
        log = LSLogOracle(x)

        RLinearAlgebra.log_update!(log, sampler, z, (), 0, A, b)

        @test length(log.error_hist) == 1
        @test log.error_hist[1] == norm(z - x)
        @test length(log.resid_hist) == 1
        @test log.resid_hist[1] == norm(A * z - b)
        @test log.iterations == 0
        @test log.converged == false
    end

    # Verify collection rate
    let
        A = rand(2,2)
        x = rand(2)
        b = A * x
        z = rand(2)

        sampler = LinSysVecRowOneRandCyclic()
        log = LSLogOracle(x, 3)

        # Test initialization of log
        RLinearAlgebra.log_update!(log, sampler, z, (), 0, A, b)

        @test length(log.error_hist) == 1
        @test log.error_hist[1] == norm(z - x)
        @test length(log.resid_hist) == 1
        @test log.resid_hist[1] == norm(A * z - b)
        @test log.iterations == 0
        @test log.converged == false

        # Test collection rate of log
        for i = 1:10
            RLinearAlgebra.log_update!(log, sampler, x + i*(z-x), (), i, A, b)
        end

        @test length(log.error_hist) == 4 # Record at 0, 3, 6, 9
        @test norm(log.error_hist - [1.0, 3.0, 6.0, 9.0] * norm(z - x)) < 1e-15
        @test length(log.resid_hist) == 4 # Record at 0, 3, 6, 9
        @test norm(log.resid_hist - [1.0, 3.0, 6.0, 9.0] * norm(A * z - b)) < 1e-15
        @test log.iterations == 10
        @test log.converged == false
    end

end

end # End Module


text:
# This file was generated, do not modify it. # hide
b2 = select(boston, Not(:NOx))
first(b2, 2)


text:
julia> using LinearAlgebraicRepresentation
julia> using Plasm
julia> using DataStructures
julia> Lar = LinearAlgebraicRepresentation

julia> v0 = [0. 1. 2. 3. 4. 5. 6. 7. 8.]
1×9 Array{Float64,2}:
 0.0  1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0

julia> v1 = collect(0.:12)
13-element Array{Float64,1}:
  0.0
  1.0
  2.0
  3.0
  4.0
  5.0
  6.0
  7.0
  8.0
  9.0
 10.0
 11.0
 12.0

julia> c0 = Lar.larGridSkeleton([1,1,1])(0)
8-element Array{Array{Int64,1},1}:
 [1]
 [2]
 [3]
 [4]
 [5]
 [6]
 [7]
 [8]

julia> m0 = (v0, c0)
([0.0 1.0 … 7.0 8.0], Array{Int64,1}[[1], [2], [3], [4], [5], [6], [7], [8]])

julia> c1 = [[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8]]
12-element Array{Array{Int64,1},1}:
 [1, 2]
 [3, 4]
 [5, 6]
 [7, 8]
 [1, 3]
 [2, 4]
 [5, 7]
 [6, 8]
 [1, 5]
 [2, 6]
 [3, 7]
 [4, 8]

julia> v1 = hcat(collect(0.:12)...)
1×13 Array{Float64,2}:
 0.0  1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  9.0  10.0  11.0  12.0

julia> m0 = (v0,c0)
([0.0 1.0 … 7.0 8.0], Array{Int64,1}[[1], [2], [3], [4], [5], [6], [7], [8]])

julia> m1 = (v1,c1)
([0.0 1.0 … 11.0 12.0], Array{Int64,1}[[1, 2], [3, 4], [5, 6], [7, 8], [1, 3], [2, 4], [5, 7], [6, 8], [1, 5], [2, 6], [3, 7], [4, 8]])

julia> Plasm.view(Lar.larModelProduct( m1,m0 ))

julia> m11 = Lar.larModelProduct( m0, m1 )
([0.0 0.0 … 8.0 8.0; 0.0 1.0 … 11.0 12.0], Array{Int64,1}[[1, 2], [3, 4], [5, 6], [7, 8], [1, 3], [2, 4], [5, 7], [6, 8], [1, 5], [2, 6]  …  [96, 97], [98, 99], [92, 94], [93, 95], [96, 98], [97, 99], [92, 96], [93, 97], [94, 98], [95, 99]])

julia> m12 = Lar.larModelProduct( m1, m0 )
([0.0 0.0 … 12.0 12.0; 0.0 1.0 … 7.0 8.0], Array{Int64,1}[[1, 10], [2, 11], [3, 12], [4, 13], [5, 14], [6, 15], [7, 16], [8, 17], [19, 28], [20, 29]  …  [25, 61], [26, 62], [28, 64], [29, 65], [30, 66], [31, 67], [32, 68], [33, 69], [34, 70], [35, 71]])

julia> Lar.Struct([ m11,m12 ])
LinearAlgebraicRepresentation.Struct(Tuple{Array{Float64,2},Array{Array{Int64,1},1}}[([0.0 0.0 … 8.0 8.0; 0.0 1.0 … 11.0 12.0], [[1, 2], [3, 4], [5, 6], [7, 8], [1, 3], [2, 4], [5, 7], [6, 8], [1, 5], [2, 6]  …  [96, 97], [98, 99], [92, 94], [93, 95], [96, 98], [97, 99], [92, 96], [93, 97], [94, 98], [95, 99]]), ([0.0 0.0 … 12.0 12.0; 0.0 1.0 … 7.0 8.0], [[1, 10], [2, 11], [3, 12], [4, 13], [5, 14], [6, 15], [7, 16], [8, 17], [19, 28], [20, 29]  …  [25, 61], [26, 62], [28, 64], [29, 65], [30, 66], [31, 67], [32, 68], [33, 69], [34, 70], [35, 71]])], Array{Float64,2}[[0.0; 0.0], [12.0; 12.0]], "10562118164415904967", 2, "feature")

julia> m2 = Lar.struct2lar(Lar.Struct([ m11,m12 ]))
([0.0 0.0 … 7.0 7.0; 0.0 1.0 … 6.0 7.0], Array{Int64,1}[[1, 2], [3, 4], [5, 6], [7, 8], [1, 3], [2, 4], [5, 7], [6, 8], [1, 5], [2, 6]  …  [23, 55], [24, 56], [25, 57], [26, 58], [27, 59], [28, 60], [29, 61], [30, 62], [31, 63], [32, 64]])

julia> Plasm.view(m2)

julia> m21 = Lar.larModelProduct( m2,m1 )
([0.0 0.0 … 7.0 7.0; 0.0 0.0 … 7.0 7.0; 0.0 1.0 … 11.0 12.0], Array{Int64,1}[[1, 2, 14, 15], [3, 4, 16, 17], [5, 6, 18, 19], [7, 8, 20, 21], [1, 3, 14, 16], [2, 4, 15, 17], [5, 7, 18, 20], [6, 8, 19, 21], [1, 5, 14, 18], [2, 6, 15, 19]  …  [408, 409, 824, 825], [410, 411, 826, 827], [404, 406, 820, 822], [405, 407, 821, 823], [408, 410, 824, 826], [409, 411, 825, 827], [404, 408, 820, 824], [405, 409, 821, 825], [406, 410, 822, 826], [407, 411, 823, 827]])

julia> Plasm.view(m21)

julia> m20 = Lar.larModelProduct( m1,m1 )
([0.0 0.0 … 12.0 12.0; 0.0 1.0 … 11.0 12.0], Array{Int64,1}[[1, 2, 14, 15], [3, 4, 16, 17], [5, 6, 18, 19], [7, 8, 20, 21], [1, 3, 14, 16], [2, 4, 15, 17], [5, 7, 18, 20], [6, 8, 19, 21], [1, 5, 14, 18], [2, 6, 15, 19]  …  [44, 45, 96, 97], [46, 47, 98, 99], [40, 42, 92, 94], [41, 43, 93, 95], [44, 46, 96, 98], [45, 47, 97, 99], [40, 44, 92, 96], [41, 45, 93, 97], [42, 46, 94, 98], [43, 47, 95, 99]])

julia> Plasm.view(m20)

julia> m20 = Lar.larModelProduct( m20,m0 )
([0.0 0.0 … 12.0 12.0; 0.0 0.0 … 12.0 12.0; 0.0 1.0 … 7.0 8.0], Array{Int64,1}[[1, 10, 118, 127], [2, 11, 119, 128], [3, 12, 120, 129], [4, 13, 121, 130], [5, 14, 122, 131], [6, 15, 123, 132], [7, 16, 124, 133], [8, 17, 125, 134], [19, 28, 136, 145], [20, 29, 137, 146]  …  [376, 412, 844, 880], [377, 413, 845, 881], [379, 415, 847, 883], [380, 416, 848, 884], [381, 417, 849, 885], [382, 418, 850, 886], [383, 419, 851, 887], [384, 420, 852, 888], [385, 421, 853, 889], [386, 422, 854, 890]])

julia> m20 = Lar.larModelProduct( m1,m1 )
([0.0 0.0 … 12.0 12.0; 0.0 1.0 … 11.0 12.0], Array{Int64,1}[[1, 2, 14, 15], [3, 4, 16, 17], [5, 6, 18, 19], [7, 8, 20, 21], [1, 3, 14, 16], [2, 4, 15, 17], [5, 7, 18, 20], [6, 8, 19, 21], [1, 5, 14, 18], [2, 6, 15, 19]  …  [44, 45, 96, 97], [46, 47, 98, 99], [40, 42, 92, 94], [41, 43, 93, 95], [44, 46, 96, 98], [45, 47, 97, 99], [40, 44, 92, 96], [41, 45, 93, 97], [42, 46, 94, 98], [43, 47, 95, 99]])

julia> m22 = Lar.larModelProduct( m20,m0 )
([0.0 0.0 … 12.0 12.0; 0.0 0.0 … 12.0 12.0; 0.0 1.0 … 7.0 8.0], Array{Int64,1}[[1, 10, 118, 127], [2, 11, 119, 128], [3, 12, 120, 129], [4, 13, 121, 130], [5, 14, 122, 131], [6, 15, 123, 132], [7, 16, 124, 133], [8, 17, 125, 134], [19, 28, 136, 145], [20, 29, 137, 146]  …  [376, 412, 844, 880], [377, 413, 845, 881], [379, 415, 847, 883], [380, 416, 848, 884], [381, 417, 849, 885], [382, 418, 850, 886], [383, 419, 851, 887], [384, 420, 852, 888], [385, 421, 853, 889], [386, 422, 854, 890]])

julia> Plasm.view(m22)

julia> Lar.struct2lar(Lar.Struct([ m22,m21 ]))
([0.0 0.0 … 7.0 7.0; 0.0 1.0 … 6.0 7.0; 0.0 0.0 … 7.0 7.0], Array{Int64,1}[[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16], [17, 18, 19, 20], [21, 22, 23, 24], [25, 26, 27, 28], [29, 30, 31, 32], [33, 34, 35, 36], [37, 38, 39, 40]  …  [244, 248, 500, 504], [252, 256, 508, 512], [228, 236, 484, 492], [232, 240, 488, 496], [244, 252, 500, 508], [248, 256, 504, 512], [228, 244, 484, 500], [232, 248, 488, 504], [236, 252, 492, 508], [240, 256, 496, 512]])

julia> m32 = Lar.struct2lar(Lar.Struct([ m22,m21 ]))
([0.0 0.0 … 7.0 7.0; 0.0 1.0 … 6.0 7.0; 0.0 0.0 … 7.0 7.0], Array{Int64,1}[[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16], [17, 18, 19, 20], [21, 22, 23, 24], [25, 26, 27, 28], [29, 30, 31, 32], [33, 34, 35, 36], [37, 38, 39, 40]  …  [244, 248, 500, 504], [252, 256, 508, 512], [228, 236, 484, 492], [232, 240, 488, 496], [244, 252, 500, 508], [248, 256, 504, 512], [228, 244, 484, 500], [232, 248, 488, 504], [236, 252, 492, 508], [240, 256, 496, 512]])

julia> Plasm.view(m32)


text:
"PDDL planning problem."
abstract type Problem end

"Returns the name of a problem."
get_name(problem::Problem) = error("Not implemented.")

"Returns the name of a problem's domain."
get_domain_name(problem::Problem) = error("Not implemented.")

"Returns an iterator over objects in a `problem`."
get_objects(problem::Problem) = error("Not implemented.")

"Returns a map from problem objects to their types."
get_objtypes(problem::Problem) = error("Not implemented.")

"Returns a list of terms that determine the initial state."
get_init_terms(problem::Problem) = error("Not implemented.")

"Returns the goal specification of a problem."
get_goal(problem::Problem) = error("Not implemented.")

"Returns the metric specification of a problem."
get_metric(problem::Problem) = error("Not implemented.")

"Returns the constraint specification of a problem."
get_constraints(problem::Problem) = error("Not implemented.")


text:
# This file contains all linear algebra code
# All functions are defined and called from the C library for interoperability

# Vector2
const Vector2 = NamedTuple{(:x, :y), Tuple{Float64, Float64}}

function Vector2(x::Float64, y::Float64)
  return Vector2((x, y))
end

function Base.:+(a::Vector2, b::Vector2)
  return ccall((:add2, @fullLibraryPath), Vector2, (Vector2, Vector2), a, b)
end

function Base.:-(a::Vector2, b::Vector2)
  return ccall((:subtract2, @fullLibraryPath), Vector2, (Vector2, Vector2), a, b)
end

function dot(a::Vector2, b::Vector2)
  return ccall((:dot2, @fullLibraryPath), Cdouble, (Vector2, Vector2), a, b)
end

function Base.:*(v::Vector2, s::Float64)
  return ccall((:scalar2, @fullLibraryPath), Vector2, (Vector2, Cdouble), v, s)
end

function magnitude(v::Vector2)
  return ccall((:magnitude2, @fullLibraryPath), Cdouble, (Vector2,), v)
end

function normalize(v::Vector2)
  return ccall((:normalize2, @fullLibraryPath), Vector2, (Vector2,), v)
end

# Vector3
const Vector3 = NamedTuple{(:x, :y, :z), Tuple{Float64, Float64, Float64}}

function Vector3(x::Float64, y::Float64, z::Float64)
  return Vector3((x, y, z))
end

function Base.:+(a::Vector3, b::Vector3)
  return ccall((:add3, @fullLibraryPath), Vector3, (Vector3, Vector3), a, b)
end

function Base.:-(a::Vector3, b::Vector3)
  return ccall((:subtract3, @fullLibraryPath), Vector3, (Vector3, Vector3), a, b)
end

function dot(a::Vector3, b::Vector3)
  return ccall((:dot3, @fullLibraryPath), Cdouble, (Vector3, Vector3), a, b)
end

function Base.:*(a::Vector3, b::Vector3)
  return ccall((:cross, @fullLibraryPath), Vector3, (Vector3, Vector3), a, b)
end

function Base.:*(v::Vector3, s::Float64)
  return ccall((:scalar3, @fullLibraryPath), Vector3, (Vector3, Cdouble), v, s)
end

function magnitude(v::Vector3)
  return ccall((:magnitude3, @fullLibraryPath), Cdouble, (Vector3,), v)
end

function normalize(v::Vector3)
  return ccall((:normalize3, @fullLibraryPath), Vector3, (Vector3,), v)
end

# Vector4
const Vector4 = NamedTuple{(:x, :y, :z, :w), Tuple{Float64, Float64, Float64, Float64}}

function Vector4(x::Float64, y::Float64, z::Float64, w::Float64)
  return Vector4((x, y, z, w))
end

function Base.:+(a::Vector4, b::Vector4)
  return ccall((:add4, @fullLibraryPath), Vector4, (Vector4, Vector4), a, b)
end

function Base.:-(a::Vector4, b::Vector4)
  return ccall((:subtract4, @fullLibraryPath), Vector4, (Vector4, Vector4), a, b)
end

function dot(a::Vector4, b::Vector4)
  return ccall((:dot4, @fullLibraryPath), Cdouble, (Vector4, Vector4), a, b)
end

function Base.:*(v::Vector4, s::Float64)
  return ccall((:scalar4, @fullLibraryPath), Vector4, (Vector4, Float64), v, s)
end

function magnitude(v::Vector4)
  return ccall((:magnitude4, @fullLibraryPath), Cdouble, (Vector4,), v)
end

function normalize(v::Vector4)
  return ccall((:normalize4, @fullLibraryPath), Vector4, (Vector4,), v)
end

# Matrices
const Matrix4x4 = NamedTuple{(:a11, :a12, :a13, :a14, :a21, :a22, :a23, :a24, :a31, :a32, :a33, :a34, :a41, :a42, :a43, :a44),
  Tuple{Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64}}

function Base.:*(a::Matrix4x4, b::Matrix4x4)
  return ccall((:multiplyMatrix4x4, @fullLibraryPath), Matrix4x4, (Matrix4x4, Matrix4x4), a, b)
end

function identityMatrix()
  return ccall((:identityMatrix4x4, @fullLibraryPath), Matrix4x4, ())
end

function invertMatrix(m)
  return ccall((:inverseMatrix4x4, @fullLibraryPath), Matrix4x4, (Matrix4x4,), m)
end

function translationMatrix(translation::Vector3)
  return ccall((:translationMatrix, @fullLibraryPath), Matrix4x4, (Vector3,), translation)
end

function scalingMatrix(scale::Vector3)
  return ccall((:scalingMatrix, @fullLibraryPath), Matrix4x4, (Vector3,), scale)
end

function rotationXMatrix(angle::Float64)
  return ccall((:rotationXMatrix, @fullLibraryPath), Matrix4x4, (Cdouble,), angle)
end

function rotationYMatrix(angle::Float64)
  return ccall((:rotationYMatrix, @fullLibraryPath), Matrix4x4, (Cdouble,), angle)
end

function rotationZMatrix(angle::Float64)
  return ccall((:rotationZMatrix, @fullLibraryPath), Matrix4x4, (Cdouble,), angle)
end

function Base.:*(a::Matrix4x4, b::Vector4)
  return ccall((:matrix4x4timesVector4, @fullLibraryPath), Vector4, (Matrix4x4, Vector4), a, b)
end

function perspective(zNear::Float64, zFar::Float64, fovy::Float64, aspect::Float64)
  return ccall((:perspectiveProjectionMatrix, @fullLibraryPath), Matrix4x4, (Float64, Float64, Float64, Float64), zNear, zFar, fovy, aspect)
end

max_stars_repo_name: jonathanbuchanan/PhysicsEngine

max_issues_repo_name: jonathanbuchanan/PhysicsEngine

max_forks_repo_name: jonathanbuchanan/PhysicsEngine


text:
module Todo

# Julia package for using and working with the todo.txt format (defined [here](https://github.com/todotxt/todo.txt))

export TodoTask, isComplete, priority, completed, created, description, projects, contexts, tags, load_todos, save_todos
export @todo_str, todos, todo_hits, reset_hits, clear_todos

using Dates: Dates, Date, format, @dateformat_str

#todo"more accurately parse description into a vector of strings, projects, contexts, and tags"
struct TodoTask
	isComplete::Bool
	priority::Union{Char, Nothing}
	completed::Union{Date, Nothing}
	created::Union{Date, Nothing}
	description::String

	function TodoTask(isComplete::Bool, priority::Union{Char, Nothing}, completed::Union{Date, Nothing}, created::Union{Date, Nothing}, description::AbstractString)
		completed !== nothing && created !== nothing && completed < created && error("TodoTask cannot be completed before it was created")
		completed !== nothing && !isComplete && error("TodoTask cannot have a completion date yet be marked incomplete")
		priority === nothing || 'A' <= priority <= 'Z' || error("TodoTask priority must a capital letter [A-Z]")

		desc = String(description)
		match(r"^x ", desc) === nothing || error("TodoTask description may not begin with a completion marking, e.g. `x `")
		match(r"^\([A-Z]\)", desc) === nothing || error("TodoTask description may not begin with a priority, e.g. `(A) `")
		match(r"^[\d]{4}-[\d]{2}-[\d]{2} ", desc) === nothing || error("TodoTask description may not begin with a date, e.g. `2019-01-23 `")

		#todo"ensure desc doesnt contain control chars"
		#todo"ensure tags, contexts, projects contain no whitespace chars"

		return new(isComplete, priority, completed, created, desc)
	end
end

function TodoTask(line::AbstractString = ""; kwargs...)
	foreach(k -> k in fieldnames(TodoTask) || error("Keyword argument $(k) is not a valid field for TodoTask ($(fieldnames(TodoTask)...))"), keys(kwargs))

	m = match(r"^(?:(x)[ ]+)?(?:\(([A-Z])\)[ ]+)?(?:([\d]{4}-[\d]{2}-[\d]{2})[ ]+)?(?:([\d]{4}-[\d]{2}-[\d]{2})[ ]+)?(.*)$", escape_string(line))
	m === nothing && error("Invalid line used to construct a TodoTask `$(line)`")

	(isC, pri, com, cre, des) = m.captures
	isC = isC === nothing ? false : true
	pri = pri === nothing ? nothing : pri[1]
	(com, cre) = cre === nothing ? (nothing, com) : (com, cre)
	com = com === nothing ? nothing : Date(com, dateformat"yyyy-mm-dd")
	cre = cre === nothing ? nothing : Date(cre, dateformat"yyyy-mm-dd")
	des = unescape_string(des)

	return TodoTask(
		get(kwargs, :isComplete, isC),
		get(kwargs, :priority, pri),
		get(kwargs, :completed, com),
		get(kwargs, :created, cre),
		get(kwargs, :description, des),
	)
end

isComplete(task::TodoTask) = task.isComplete
priority(task::TodoTask) = task.priority
completed(task::TodoTask) = task.completed
created(task::TodoTask) = task.created
description(task::TodoTask) = task.description
projects(task::TodoTask) = map(m -> m.captures[1], collect(eachmatch(r"(?:^| )\+([^ ]+)", task.description)))
contexts(task::TodoTask) = map(m -> m.captures[1], collect(eachmatch(r"(?:^| )@([^ ]+)", task.description)))
tags(task::TodoTask) = map(m -> m.captures[1] => m.captures[2], collect(eachmatch(r"(?:^| )([^ :+@][^ :]*):([^ ]+)", task.description)))

Base.string(task::TodoTask) = join((
	task.isComplete ? "x " : "",
	task.priority !== nothing ? "($(task.priority)) " : "",
	task.completed !== nothing ? format(task.completed, dateformat"yyyy-mm-dd ") : "",
	task.created !== nothing ? format(task.created, dateformat"yyyy-mm-dd ") : "",
	escape_string(task.description)
))

function load_todos(filePath)
	return open(filePath) do file
		return map(line -> TodoTask(unescape_string(line)), filter(!isempty, readlines(file)))
	end
end

function save_todos(filePath, todos; format = :todo, kwargs...)
	format in (:todo, :markdown) || error("`$(format)` is not a valid format for saving.")

	open(filePath, "w+") do file
		if format === :todo
			foreach(task::TodoTask -> println(file, string(task)), todos)
		elseif format === :markdown
			println(file, get(kwargs, :heading, "## Todos"))
			println(file)
			foreach(todos) do todo::TodoTask
				print(file, "- [$(todo.isComplete ? 'x' : ' ')] ")
				println(file, replace(todo.description, '\n' => "  \n"))
			end
			println(file)
		end
	end
	return nothing
end

_todos = Dict{TodoTask, Int}()

todos() = collect(keys(_todos))
todo_hits() = filter(((todo, hits),) -> hits > 0, _todos)
reset_hits() = foreach(key -> _todos[key] = 0, keys(_todos))

macro todo_str(line::String) ; return _add_todo(line, __source__, __module__) ; end

function _add_todo(line::String, src::LineNumberNode, mod::Module)
	# do nothing when file field is not provided or when in the REPL
	src.file === nothing && return quote end
	m = match(r"^REPL\[\d+\]$", string(src.file))
	m === nothing || return quote end

	task = TodoTask(line)

	#todo"escape spaces in profile and context strings"

	project = string(mod)
	project = replace(project, r"^Main\." => "")
	task = project in projects(task) ? task : TodoTask(unescape_string(string(task)), description = "$(task.description) +$(project)")

	context = "$(relpath(string(src.file))):$(src.line)"
	task = context in contexts(task) ? task : TodoTask(unescape_string(string(task)), description = "$(task.description) @$(context)")

	line = unescape_string(string(task))
	task = TodoTask(line)
	haskey(Todo._todos, task) || push!(Todo._todos, task => 0)

	return Todo.TRACK_TODOS ? :(let ; Todo._todos[TodoTask($(line))] += 1 ; nothing ; end) : nothing
end

function __init__()
	global TRACK_TODOS = isdefined(Main, :TRACK_TODOS) ? (Main.TRACK_TODOS === true) : isinteractive()
end

end


text:
module synthtest
using JSON

include("../src/Ebic.jl")
include("metrics.jl")

using .Ebic: run_ebic
using .metrics: eval_metrics

const OUTPUT_DIR = "results/Ebic.jl"

function test_dataset(dataset_path)
    out_path = joinpath(OUTPUT_DIR, splitpath(dataset_path)[end])
    println(out_path)
    isdir(out_path) || mkpath(out_path)

    for (root, _, files) in walkdir(dataset_path)
        isempty(files) && continue

        println("""
        ####################################
        Starting test case: '$(basename(root))'
        ####################################""")

        input_paths = Vector()
        biclusters_paths = Vector()
        for file in files
            if occursin("hidden", file)
                push!(biclusters_paths, joinpath(root, file))
            else
                push!(input_paths, joinpath(root, file))
            end
        end

        if length(input_paths) != length(biclusters_paths)
            throw(ErrorException("Something is wrong in '$root'."))
        end

        test_case_results = Vector()
        for (input_path, ground_truth_path) in zip(input_paths, biclusters_paths)
            println("##################################")
            println("Testing: '$(basename(input_path))'")
            println("Ground truth: '$(basename(ground_truth_path))'")

            ground_truth::Vector = JSON.parsefile(ground_truth_path)

            result = run_ebic(
                input_path,
                verbose = true,
                max_iterations = 20_000,
                max_biclusters = length(ground_truth),
                overlap_threshold = 0.75,
                negative_trends = true,
                approx_trends_ratio = 0.85f0,
                best_bclrs_stats = false,
            )

            result["input_data"] = input_path
            result["ground_truth"] = ground_truth_path

            relevance, recovery, ce =
                eval_metrics(result["biclusters"], input_path, ground_truth)

            result["relevance"] = relevance
            result["recovery"] = recovery
            result["ce"] = ce

            println("Prelic relevance   : $(relevance)")
            println("Prelic recovery    : $(recovery)")
            println("Clustering error   : $(ce)")
            println("Last iter tabu hits: $(result["last_iter_tabu_hits"])")

            push!(test_case_results, result)
        end

        open(joinpath(out_path, "$(basename(root))_res.json"), "w") do f
            JSON.print(f, test_case_results)
        end
    end
end

test_unibic() = test_dataset("data/unibic/")
test_recbic_sup() = test_dataset("data/recbic_sup/")
test_recbic_main() = test_dataset("data/recbic_maintext/")

function main()
    test_unibic()
    test_recbic_main()
    test_recbic_sup()
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end
end


text:
export OpenStreetMapSpace, OSMSpace, TEST_MAP
export osm_random_road_position,
    osm_plan_route,
    osm_map_coordinates,
    osm_road_length,
    osm_is_stationary,
    osm_random_route!,
    osm_latlon,
    osm_intersection,
    osm_road

struct OpenStreetMapSpace <: DiscreteSpace
    m::OpenStreetMapX.MapData
    s::Vector{Vector{Int}}
end

"""
    OpenStreetMapSpace(path::AbstractString; kwargs...)
Create a space residing on the Open Street Map (OSM) file provided via `path`.

The abbreviation `OSMSpace` may be used interchangeably.

Much of the functionality of this space is provided by interfacing with
[OpenStreetMapX.jl](https://github.com/pszufe/OpenStreetMapX.jl), for example the two
keyword arguments `use_cache = false` and `trim_to_connected_graph = true` can be
passed into the `OpenStreetMapX.get_map_data` function.

For details on how to obtain an OSM file for your use case, consult the OpenStreetMapX.jl
README. We provide a variable `TEST_MAP` to use as a `path` for testing.

This space represents the underlying map as a *continuous* entity choosing accuracy over
performance. An example of its usage can be found in [Zombie Outbreak](@ref).

If your solution can tolerate routes to and from intersections only, a
faster implementation can be achieved by using the
[graph representation](https://pszufe.github.io/OpenStreetMapX.jl/stable/reference/#OpenStreetMapX.MapData)
of your map provided by OpenStreetMapX.jl. For tips on how to implement this, see our
integration example: [Social networks with LightGraphs.jl](@ref).

## The OSMAgent

The base properties for an agent residing on an `OSMSpace` are as follows:
```julia
mutable struct OSMAgent <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int,Float64}
    route::Vector{Int}
    destination::Tuple{Int,Int,Float64}
end
```

Current `pos`ition and `destination` tuples are represented as
`(start intersection index, finish intersection index, distance travelled in meters)`.
The `route` is an ordered list of intersections, providing a path to reach `destination`.

Further details can be found in [`OSMAgent`](@ref).

## Routing

There are two ways to generate a route, depending on the situation.
1. [`osm_plan_route`](@ref), which provides `:shortest` and `:fastest` paths (with the option of a `return_trip`) between intersections or positions.
2. [`osm_random_route!`](@ref), choses a new `destination` an plans a new path to it; overriding the current route (if any).
"""
function OpenStreetMapSpace(
    path::AbstractString;
    use_cache = false,
    trim_to_connected_graph = true,
)
    m = get_map_data(path; use_cache, trim_to_connected_graph)
    agent_positions = [Int[] for i in 1:nv(m.g)]
    return OpenStreetMapSpace(m, agent_positions)
end

function Base.show(io::IO, s::OpenStreetMapSpace)
    print(
        io,
        "OpenStreetMapSpace with $(length(s.m.roadways)) roadways and $(length(s.m.intersections)) intersections",
    )
end

const OSMSpace = OpenStreetMapSpace

const TEST_MAP =
    joinpath(dirname(pathof(OpenStreetMapX)), "..", "test", "data", "reno_east3.osm")

#######################################################################################
# Custom functions for OSMSpace
#######################################################################################

# EXPORTED
"""
    osm_random_road_position(model::ABM{OpenStreetMapSpace})

Similar to `random_position`, but rather than providing only intersections, this method
returns a location somewhere on a road heading in a random direction.
"""
function osm_random_road_position(model::ABM{<:OpenStreetMapSpace})
    ll = generate_point_in_bounds(model.space.m)
    return osm_road(ll, model)
end

"""
    osm_random_route!(agent, model::ABM{<:OpenStreetMapSpace})

Plan a new random route for the agent, by selecting a random destination and
planning a route from the agent's current position. Overwrite any current route.
"""
function osm_random_route!(agent, model::ABM{<:OpenStreetMapSpace})
    agent.destination = osm_random_road_position(model)
    agent.route = osm_plan_route(agent.pos, agent.destination, model)
    return nothing
end

"""
    osm_plan_route(start, finish, model::ABM{<:OpenStreetMapSpace};
                   by = :shortest, return_trip = false, kwargs...)

Generate a list of intersections between `start` and `finish` points on the map.
`start` and `finish` can either be intersections (`Int`) or positions
(`Tuple{Int,Int,Float64}`).

When either point is a position, the associated intersection index will be removed from
the route to avoid double counting.

Route is planned via the shortest path by default (`by = :shortest`), but can also be
planned `by = :fastest`. Road speeds are needed for this method which can be passed in via
extra keyword arguments. Consult the OpenStreetMapX documentation for more details.

If `return_trip = true`, a route will be planned from start -> finish -> start.
"""
function osm_plan_route(
    start::Int,
    finish::Int,
    model::ABM{<:OpenStreetMapSpace};
    by = :shortest,
    return_trip = false,
    kwargs...,
)
    @assert by ∈ (:shortest, :fastest) "Can only plan route by :shortest or :fastest"
    planner = by == :shortest ? shortest_route : fastest_route
    route = if return_trip
        planner(
            model.space.m,
            model.space.m.n[start],
            model.space.m.n[finish],
            model.space.m.n[start];
            kwargs...,
        )[1]
    else
        planner(model.space.m, model.space.m.n[start], model.space.m.n[finish]; kwargs...)[1]
    end
    map(p -> getindex(model.space.m.v, p), route)
end

function osm_plan_route(
    start::Tuple{Int,Int,Float64},
    finish::Int,
    model::ABM{<:OpenStreetMapSpace};
    return_trip = false,
    kwargs...,
)
    path = if return_trip
        plan_return_route(start[2], finish, start[1], model; kwargs...)
    else
        osm_plan_route(start[2], finish, model; kwargs...)
    end
    ## Since we start on an edge, there are two possibilities here.
    ## 1. The route wants us to turn around, thus next id en-route will
    ## be pos[1]. That's fine.
    ## 2. The route wants us to move on, but start will be in the list,
    ## so we need to drop that.
    if !isempty(path) && path[1] == start[2]
        popfirst!(path)
    end
    return path
end

function osm_plan_route(
    start,
    finish::Tuple{Int,Int,Float64},
    model::ABM{<:OpenStreetMapSpace};
    kwargs...,
)
    path = osm_plan_route(start, finish[1], model; kwargs...)
    isempty(path) || pop!(path)
    return path
end

function plan_return_route(
    start::Int,
    middle::Int,
    finish::Int,
    model::ABM{<:OpenStreetMapSpace};
    by = :shortest,
    kwargs...,
)
    planner = by == :shortest ? shortest_route : fastest_route
    route = planner(
        model.space.m,
        model.space.m.n[start],
        model.space.m.n[middle],
        model.space.m.n[finish];
        kwargs...,
    )[1]
    map(p -> getindex(model.space.m.v, p), route)
end

"""
    osm_latlon(pos, model)
    osm_latlon(agent, model)

Return (latitude, longitude) of current road or intersection position.
"""
osm_latlon(pos::Int, model::ABM{<:OpenStreetMapSpace}) = latlon(model.space.m, pos)

function osm_latlon(pos::Tuple{Int,Int,Float64}, model::ABM{<:OpenStreetMapSpace})
    if pos[1] != pos[2]
        start = get_EastNorthUp_coordinate(pos[1], model)
        finish = get_EastNorthUp_coordinate(pos[2], model)
        travelled = pos[3] / osm_road_length(pos, model)
        enu_coord = ENU(
            getX(start) * (1 - travelled) + getX(finish) * travelled,
            getY(start) * (1 - travelled) + getY(finish) * travelled,
            getZ(start) * (1 - travelled) + getZ(finish) * travelled,
        )
        lla = LLA(enu_coord, model.space.m.bounds)
        return (lla.lat, lla.lon)
    else
        return latlon(model.space.m, pos[1])
    end
end

osm_latlon(agent::A, model::ABM{<:OpenStreetMapSpace,A}) where {A<:AbstractAgent} =
    osm_latlon(agent.pos, model)

"""
    osm_intersection(latlon::Tuple{Float64,Float64}, model::ABM{<:OpenStreetMapSpace})

Return the nearest intersection position to (latitude, longitude).
Quicker, but less precise than [`osm_road`](@ref).
"""
function osm_intersection(ll::Tuple{Float64,Float64}, model::ABM{<:OpenStreetMapSpace})
    idx = getindex(model.space.m.v, point_to_nodes(ll, model.space.m))
    return (idx, idx, 0.0)
end

"""
    osm_road(latlon::Tuple{Float64,Float64}, model::ABM{<:OpenStreetMapSpace})

Return a location on a road nearest to (latitude, longitude).
Slower, but more precise than [`osm_intersection`](@ref).
"""
function osm_road(ll::Tuple{Float64,Float64}, model::ABM{<:OpenStreetMapSpace})
    ll_enu = ENU(LLA(ll...), model.space.m.bounds)
    P = (ll_enu.east, ll_enu.north, ll_enu.up)

    # This is one index, close to the position.
    idx = getindex(model.space.m.v, point_to_nodes(ll, model.space.m))
    idx_enu = get_EastNorthUp_coordinate(idx, model)

    candidates = Tuple{Tuple{Int,Int,Float64},Float64}[]
    # This separation is only useful for one-way street situations.
    # In case of two way streets, either side may be returned with
    # little penalty.
    if abs(ll_enu.east - idx_enu.east) > abs(ll_enu.north - idx_enu.north)
        # idx is the first position
        nps = nearby_positions(idx, model; neighbor_type = :out)
        isempty(nps) && return (idx, idx, 0.0)
        np_enus = map(np -> get_EastNorthUp_coordinate(np, model), nps)
        A = (idx_enu.east, idx_enu.north, idx_enu.up)
        for (np_enu, np) in zip(np_enus, nps)
            B = (np_enu.east, np_enu.north, np_enu.up)
            closest = orthognonal_projection(A, B, P)
            candidate = (idx, np, distance(np_enu, ENU(closest...)))
            push!(candidates, (candidate, sum(abs.(osm_latlon(candidate, model) .- ll))))
        end
    else
        # idx is the second position
        nps = nearby_positions(idx, model; neighbor_type = :in)
        isempty(nps) && return (idx, idx, 0.0)
        np_enus = map(np -> get_EastNorthUp_coordinate(np, model), nps)
        B = (idx_enu.east, idx_enu.north, idx_enu.up)
        for (np_enu, np) in zip(np_enus, nps)
            A = (np_enu.east, np_enu.north, np_enu.up)
            closest = orthognonal_projection(A, B, P)
            candidate = (np, idx, distance(np_enu, ENU(closest...)))
            push!(candidates, (candidate, sum(abs.(osm_latlon(candidate, model) .- ll))))
        end
    end
    bestidx = findmin(last.(candidates))[2]
    return first(candidates[bestidx])
end

function orthognonal_projection(A, B, P)
    M = B .- A
    t0 = dot(M, P .- A) / dot(M, M)
    return A .+ t0 .* M
end

"""
    osm_map_coordinates(agent, model::ABM{OpenStreetMapSpace})

Return a set of coordinates for an agent on the underlying map. Useful for plotting.
"""
function osm_map_coordinates(agent, model)
    if agent.pos[1] != agent.pos[2]
        start = get_EastNorthUp_coordinate(agent.pos[1], model)
        finish = get_EastNorthUp_coordinate(agent.pos[2], model)
        travelled = agent.pos[3] / osm_road_length(agent.pos, model)
        (
            getX(start) * (1 - travelled) + getX(finish) * travelled,
            getY(start) * (1 - travelled) + getY(finish) * travelled,
        )
    else
        position = get_EastNorthUp_coordinate(agent.pos[1], model)
        (getX(position), getY(position))
    end
end

"""
    osm_road_length(start::Int, finish::Int, model)
    osm_road_length(pos::Tuple{Int,Int,Float64}, model)

Return the road length (in meters) between two intersections given by intersection ids.
"""
osm_road_length(pos::Tuple{Int,Int,Float64}, model::ABM{<:OpenStreetMapSpace}) =
    osm_road_length(pos[1], pos[2], model)
osm_road_length(p1::Int, p2::Int, model::ABM{<:OpenStreetMapSpace}) =
    model.space.m.w[p1, p2]

"""
    osm_is_stationary(agent)

Return `true` if agent has no route left to follow and is therefore standing still.
"""
function osm_is_stationary(agent::A) where {A<:AbstractAgent}
    return agent.pos == agent.destination && isempty(agent.route)
end

#HELPERS, NOT EXPORTED

"""
    get_EastNorthUp_coordinate(pos::Int, model)

Return an East-North-Up coordinate value for index `pos`.
"""
get_EastNorthUp_coordinate(pos::Int, model) = model.space.m.nodes[model.space.m.n[pos]]

#######################################################################################
# Agents.jl space API
#######################################################################################

function random_position(model::ABM{<:OpenStreetMapSpace})
    ll = generate_point_in_bounds(model.space.m)
    return osm_intersection(ll, model)
end

function add_agent_to_space!(
    agent::A,
    model::ABM{<:OpenStreetMapSpace,A},
) where {A<:AbstractAgent}
    push!(model.space.s[agent.pos[1]], agent.id)
    return agent
end

function remove_agent_from_space!(
    agent::A,
    model::ABM{<:OpenStreetMapSpace,A},
) where {A<:AbstractAgent}
    prev = model.space.s[agent.pos[1]]
    ai = findfirst(i -> i == agent.id, prev)
    deleteat!(prev, ai)
    return agent
end

function move_agent!(
    agent::A,
    pos::Tuple{Int,Int,Float64},
    model::ABM{<:OpenStreetMapSpace,A},
) where {A}
    remove_agent_from_space!(agent, model)
    agent.pos = pos
    add_agent_to_space!(agent, model)
end

"""
    move_agent!(agent, model::ABM{<:OpenStreetMapSpace}, distance::Real)

Move an agent by `distance` in meters along its planned route.
"""
function move_agent!(
    agent::A,
    model::ABM{<:OpenStreetMapSpace,A},
    distance::Real,
) where {A<:AbstractAgent}

    if osm_is_stationary(agent)
        return nothing
    end

    dist_to_intersection = osm_road_length(agent.pos, model) - agent.pos[3]

    if isempty(agent.route) && agent.pos[1:2] == agent.destination[1:2]
        # Last one or two moves before destination
        to_travel = agent.destination[3] - agent.pos[3]
        if distance >= to_travel
            pos = agent.destination
        else
            pos = (agent.destination[1:2]..., agent.pos[3] + distance)
        end
    elseif distance >= dist_to_intersection
        if !isempty(agent.route)
            pos = travel!(
                agent.pos[2],
                popfirst!(agent.route),
                distance - dist_to_intersection,
                agent,
                model,
            )
        else
            # Now moving to the final destination
            pos = park(distance - dist_to_intersection, agent, model)
        end
    else
        # move up current path
        pos = (agent.pos[1], agent.pos[2], agent.pos[3] + distance)
    end

    move_agent!(agent, pos, model)
end

function travel!(start, finish, distance, agent, model)
    # Assumes we have just reached the intersection of `start` and `finish`,
    # and have `distance` left to travel.
    edge_distance = osm_road_length(start, finish, model)
    if edge_distance <= distance
        if !isempty(agent.route)
            return travel!(
                finish,
                popfirst!(agent.route),
                distance - edge_distance,
                agent,
                model,
            )
        else
            #######################################
            # TODO: The code here can be simplified by using the existing `park`.
            # alright, so here we're in a situation where the agent is imagined to be at 'start' with 'distance left to travel.
            # the route is empty, but (start,finish) does not equal agent.destination[1:2]
            # what follows is the srouce code of the function "park", but the 'virtual' agent in it has position
            # pos=(start,finish,distance-edge_distance)
            # so I replaced that and nothing else.
            distance-=edge_distance
            if finish != agent.destination[1]
                # At the end of the route, we must travel
                last_distance = osm_road_length(finish, agent.destination[1], model)
                if distance >= last_distance + agent.destination[3]
                    # We reach the destination
                    return agent.destination
                elseif distance >= last_distance
                    # We reach the final road, but not the destination
                    return (agent.destination[1:2]..., distance - last_distance)
                else
                    # We travel the final leg
                    return (finish, agent.destination[1], distance)
                end
            else
                # Reached final road
                if distance >= agent.destination[3]
                    return agent.destination
                else
                    return (agent.destination[1:2]..., distance)
                end
            end
            #######################################
        end
    else
        return (start, finish, distance)
    end
end

function park(distance, agent, model)
    # We have no route left but have not quite yet arrived at our destination.
    # Assumes that when this is called, we have just completed the current leg
    # in `agent.pos`, and we have `distance` left to travel.
    if agent.pos[2] != agent.destination[1]
        # At the end of the route, we must travel
        last_distance = osm_road_length(agent.pos[2], agent.destination[1], model)
        if distance >= last_distance + agent.destination[3]
            # We reach the destination
            return agent.destination
        elseif distance >= last_distance
            # We reach the final road, but not the destination
            return (agent.destination[1:2]..., distance - last_distance)
        else
            # We travel the final leg
            return (agent.pos[2], agent.destination[1], distance)
        end
    else
        # Reached final road
        if distance >= agent.destination[3]
            return agent.destination
        else
            return (agent.destination[1:2]..., distance)
        end
    end
end

# Nearby positions must be intersections, since edges imply a direction.
# nearby agents/ids can be on an intersection or on a road X m away.
# We cannot simply use nearby_positions to collect nearby ids then.

# Default is searching both backwards and forwards.
# I assume it would be useful to turn off one or the other at some point.
function nearby_ids(
    pos::Tuple{Int,Int,Float64},
    model::ABM{<:OpenStreetMapSpace},
    distance::Real,
    args...;
    kwargs...,
)
    current_road = osm_road_length(pos, model)
    nearby = Int[]

    close = ids_on_road(pos[1], pos[2], model)
    pos_search = pos[3] + distance
    if distance > current_road - pos[3]
        # Local search distance
        pos_search = current_road

        # Check outgoing
        search_distance = distance - (current_road - pos[3])
        search_outward_ids!(nearby, search_distance, pos[1], pos[2], model)
    end
    # Check anyone close in the forward direction
    if !isempty(close)
        for (id, dist) in close
            if pos[3] <= dist <= pos_search
                push!(nearby, id)
            end
        end
    end
    rev_search = pos[3] - distance
    if rev_search < 0
        # Local search distance
        rev_search = 0.0

        # Check incoming
        search_distance = distance - pos[3]
        search_inward_ids!(nearby, search_distance, pos[1], pos[2], model)
    end
    # Check anyone close in the reverse direction
    if !isempty(close)
        for (id, dist) in close
            if rev_search <= dist <= pos[3]
                push!(nearby, id)
            end
        end
    end
    return nearby
end

function ids_on_road(
    pos1::Int,
    pos2::Int,
    model::ABM{<:OpenStreetMapSpace},
    reverse = false,
)
    distance = osm_road_length(pos1, pos2, model)
    dist_front(d) = reverse ? distance - d : d
    dist_back(d) = reverse ? d : distance - d
    # Agents listed in the current position, filtered to current road
    # (id, distance)
    res = [
        (model[i].id, dist_front(model[i].pos[3]))
        for i in model.space.s[pos1] if model[i].pos[2] == pos2
    ]
    # Opposite direction. We must invert the distances here to obtain a relative distance from `pos`
    #NOTE: I think a complication happens here when we're looking at divided roads. They don't end up having
    #the same ids, so this comparison isn't possible. Unsure if we can do something about that using
    #OpenStreetMapX or not.
    append!(
        res,
        [
            (model[i].id, dist_back(model[i].pos[3]))
            for i in model.space.s[pos2] if model[i].pos[2] == pos1
        ],
    )
    res
end

function search_outward_ids!(
    nearby::Vector{Int},
    distance::Real,
    pos1::Int,
    pos2::Int,
    model,
)
    # find all intersections the current end position connects to
    outgoing = filter(i -> i != pos1, nearby_positions(pos2, model; neighbor_type = :out))
    # Distances for each road
    outdist = [osm_road_length(pos2, o, model) for o in outgoing]
    # Identify roads that are shorter than the search distance
    go_deeper = findall(i -> i < distance, outdist)
    # Those in the `go_deeper` category are completely covered by the distance search.
    # Any agent found in the heading or opposite directions along this road are counted.
    for i in go_deeper
        finish = outgoing[i]
        for (id, _) in ids_on_road(pos2, finish, model)
            push!(nearby, id)
        end
        # We must recursively look up branches until the search distance is met
        search_outward_ids!(nearby, distance - outdist[i], pos2, finish, model)
    end
    # Those not in the `go_deeper` category are searched in the positive direction
    setdiff!(outgoing, outgoing[go_deeper])
    to_collect = map(i -> ids_on_road(pos2, i, model), outgoing)
    filter!(!isempty, to_collect)
    if !isempty(to_collect)
        for (id, dist) in vcat(to_collect...)
            if dist <= distance
                push!(nearby, id)
            end
        end
    end
end

function search_inward_ids!(
    nearby::Vector{Int},
    distance::Real,
    pos1::Int,
    pos2::Int,
    model,
)
    # find all intersections the current start position connects to
    incoming = filter(i -> i != pos2, nearby_positions(pos1, model; neighbor_type = :in))
    # Distances for each road
    indist = [osm_road_length(i, pos1, model) for i in incoming]
    # Identify roads that are shorter than the search distance
    go_deeper = findall(i -> i < distance, indist)
    # Those in the `go_deeper` category are completely covered by the distance search.
    # Any agent found in the heading or opposite directions along this road are counted.
    for i in go_deeper
        start = incoming[i]
        for (id, _) in ids_on_road(start, pos1, model)
            push!(nearby, id)
        end
        # We must recursively look up branches until the search distance is met
        search_inward_ids!(nearby, distance - indist[i], start, pos1, model)
    end
    # Those not in the `go_deeper` category are searched in the positive direction
    setdiff!(incoming, incoming[go_deeper])
    to_collect = map(i -> ids_on_road(i, pos1, model, true), incoming)
    filter!(!isempty, to_collect)
    if !isempty(to_collect)
        for (id, dist) in vcat(to_collect...)
            if dist <= distance
                push!(nearby, id)
            end
        end
    end
end

function nearby_ids(
    agent::A,
    model::ABM{<:OpenStreetMapSpace,A},
    args...;
    kwargs...,
) where {A<:AbstractAgent}
    all = nearby_ids(agent.pos, model, args...; kwargs...)
    filter!(i -> i ≠ agent.id, all)
end

nearby_positions(pos::Tuple{Int,Int,Float64}, model, args...; kwargs...) =
    nearby_positions(pos[1], model, args...; kwargs...)

function nearby_positions(
    position::Int,
    model::ABM{<:OpenStreetMapSpace};
    neighbor_type::Symbol = :default,
)
    @assert neighbor_type ∈ (:default, :all, :in, :out)
    neighborfn = if neighbor_type == :default
        LightGraphs.neighbors
    elseif neighbor_type == :in
        LightGraphs.inneighbors
    elseif neighbor_type == :out
        LightGraphs.outneighbors
    else
        LightGraphs.all_neighbors
    end
    neighborfn(model.space.m.g, position)
end


text:
#InspectDR: Functionality for user-overwritable styles (ex: used in stylesheets)
#-------------------------------------------------------------------------------

#==Main types
===============================================================================#
abstract type AbstractStyle end
struct StyleDefaultType end #Keep style default
const StyleDefault = StyleDefaultType()

#Defines a "style" object (to implement user-overwritable styles):
mutable struct StyleType{T<:AbstractStyle}
	values::T #Active values
	defaults::T #Defaults for this style
	overwrite::Set{Symbol} #Masks out which elements user wishes to overwrite
end
StyleType(defaults::T) where T<:AbstractStyle =
	StyleType(deepcopy(defaults), defaults, Set{Symbol}([]))

#==Main functions
===============================================================================#
function overwritestyleattrib(s::Set{Symbol}, key::Symbol, value::StyleDefaultType)
	if key in s
		pop!(s, key)
	end
	return
end
overwritestyleattrib(s::Set{Symbol}, key::Symbol, value) = push!(s, key)

overwrite!(s::AbstractStyle, key::Symbol, value::StyleDefaultType) = nothing
function overwrite!(s::AbstractStyle, key::Symbol, value)
	value = convert(typeof(getfield(s, key)), value) #Need proper type
	setfield!(s, key, value)
end

function overwrite!(s::StyleType, key::Symbol, value)
	overwrite!(s.values, key, value)
	overwritestyleattrib(s.overwrite, key, value)
end

function refresh!(s::StyleType)
	proplist = fieldnames(typeof(s.defaults))
	for key in proplist
		if !(key in s.overwrite)
			value = getfield(s.defaults, key)
			setfield!(s.values, key, value)
		end
	end
	return
end

function setstyle!(s::StyleType{T}, value::T; refresh::Bool=true) where T<:AbstractStyle
	s.defaults = value
	if refresh
		refresh!(s)
	end
	return
end

#==Mutators
===============================================================================#
Base.setindex!(s::AbstractStyle, value, key::Symbol) = overwrite!(s, key, value)
Base.setindex!(s::StyleType, value, key::Symbol) = overwrite!(s, key, value)
Base.getindex(s::AbstractStyle, key::Symbol) = getfield(s, key)
Base.getindex(s::StyleType, key::Symbol) = getfield(s.values, key)

#Last line


text:
using Test
using NQCDynamics
using ComponentArrays: ComponentVector
using StatsBase: mean
using Distributions: Normal

model = Harmonic()
sim = Simulation(Atoms(1), model; temperature=5)
rp_sim = RingPolymerSimulation(Atoms(1), model, 10; temperature=5)

vector = [ComponentVector(r=fill(1, size(sim)), v=randn(size(sim))) for _=1:100]
rp_vector = [ComponentVector(r=fill(1, size(rp_sim)), v=randn(size(rp_sim))) for _=1:100]

dist = DynamicalDistribution(Normal(), [randn(size(sim)) for _=1:100], size(sim))

@testset "@estimate" begin
    avg = mean(Estimators.potential_energy.(sim, vector))
    est = Estimators.@estimate potential_energy(sim, vector)
    @test avg ≈ est

    Estimators.@estimate potential_energy(sim, dist)
end

@testset "potential_energy" begin
    function f(u)
        r = DynamicsUtils.get_positions(u)
        sum(model.m .* model.ω.^2 .* r.^2) / 2
    end
    u = rand(vector)
    @test Estimators.potential_energy(sim, u) ≈ f(u)
    u = rand(rp_vector)
    @test Estimators.potential_energy(rp_sim, u) ≈ f(u) / nbeads(rp_sim)
end

@testset "kinetic_energy" begin
    function f(u)
        v = DynamicsUtils.get_velocities(u)
        sum(v.^2) / 2
    end
    u = rand(vector)
    @test Estimators.kinetic_energy(sim, u) ≈ f(u)
end


text:
module Correlation

include("wrapper/CAPI.jl")
import .CAPI: Problem, Parameter, Score

include("mine.jl")
include("utils.jl")

end


text:
using LookingGlass
using Test

module MF
    foo(x) = 2x
    bar(x) = foo(x) + 1
    module Inner
        f() = 1
    end

    g = 2
end

@test Set(LookingGlass.module_functions_names(MF)) == Set([:foo, :bar])
@test Set(LookingGlass.module_functions(MF)) == Set([MF.foo, MF.bar])
@test Set(LookingGlass.module_objects(MF)) ⊇ Set([MF.g, MF.Inner, typeof(MF.foo), typeof(MF.bar), MF.foo, MF.bar])

@testset "func_specializations()" begin
    @test length(LookingGlass.func_specializations(MF.foo)) == 0
    MF.foo(2)
    @test length(LookingGlass.func_specializations(MF.foo)) == 1
end

@testset "func_backedges()" begin
    @test length(first(LookingGlass.func_backedges(MF.foo))[2]) == 0
    MF.bar(2)
    @test length(first(LookingGlass.func_backedges(MF.foo))[2]) == 1
end

module Outer
    const g_outer = 1
end
module MV
    gv = 2
    const cv = 2
    vec = []
    module A
        const g_a = 1
    end
    module Inner
        i_x = 3
        const i_c = 2
        const i_vec = [2]

        import ...Outer
    end
end

@test LookingGlass.module_recursive_globals_names(MV) ==
    Dict(
        MV => sort([:gv, :cv, :vec]),
        MV.Inner => sort([:i_x, :i_c, :i_vec]),
        MV.A => sort([:g_a]),
        )

@test LookingGlass.module_recursive_globals_names(MV,
                     constness = :const, mutability = :mutable) ==
    Dict(
        MV => sort([]),
        MV.Inner => sort([:i_vec]),
        MV.A => sort([]),
        )

@test LookingGlass.module_recursive_globals_names(MV,
                     constness = :const, imported = true) ==
    Dict(
        MV => sort([:cv]),
        MV.Inner => sort([:i_c, :i_vec]),
        MV.A => sort([:g_a]),
        Outer => sort([:g_outer]),  # imported via MV.Inner
        )

@test LookingGlass.module_recursive_globals(MV) ==
    Dict(
        (MV, :gv) => MV.gv,
        (MV, :cv) => MV.cv,
        (MV, :vec) => MV.vec,
        (MV.Inner, :i_x) => MV.Inner.i_x,
        (MV.Inner, :i_c) => MV.Inner.i_c,
        (MV.Inner, :i_vec) => MV.Inner.i_vec,
        (MV.A, :g_a) => MV.A.g_a,
        )

@test LookingGlass.module_recursive_globals(MV, imported=true) ==
    Dict(
        (MV, :gv) => MV.gv,
        (MV, :cv) => MV.cv,
        (MV, :vec) => MV.vec,
        (MV.Inner, :i_x) => MV.Inner.i_x,
        (MV.Inner, :i_c) => MV.Inner.i_c,
        (MV.Inner, :i_vec) => MV.Inner.i_vec,
        (MV.A, :g_a) => MV.A.g_a,
        (Outer, :g_outer) => MV.Inner.Outer.g_outer,  # imported via MV.Inner
        )

@test LookingGlass.module_recursive_globals(MV, imported=true, constness=:const) ==
    Dict(
        (MV, :cv) => MV.cv,
        (MV.Inner, :i_c) => MV.Inner.i_c,
        (MV.Inner, :i_vec) => MV.Inner.i_vec,
        (MV.A, :g_a) => MV.A.g_a,
        (Outer, :g_outer) => MV.Inner.Outer.g_outer,  # imported via MV.Inner
        )


text:
using Peacock
using PyPlot, LinearAlgebra
p = TE
cutoff = 7
#epf(x, y) = x^2 + y^2 < r^2 ? 11.43 : 1
r0 = 0.3; d = 0.32
epf(x, y) = abs(x+d) <=r0/2  && abs(y+d) <= r0/2 || abs(x-d) <=r0/2  && abs(y+d) <= r0/2 || abs(x+d) <=r0/2  && abs(y-d) <= r0/2 || abs(x-d) <=r0/2  && abs(y-d) <= r0/2 ? 1 : 25
muf(x,y) = 1
a1 = [1, 0]; a2 = [0, 1]
d1 = 0.01*100/121; d2 = 0.01*100/121
geometry = Geometry(epf, muf, a1, a2, d1, d2)
solver = Solver(geometry, cutoff, cutoff)

G = BrillouinZoneCoordinate(  0,   0, "Γ")
X = BrillouinZoneCoordinate(1/2,   0, "X")
M = BrillouinZoneCoordinate(1/2, 1/2, "M")
ks = [G,X,M,G]
figure(); plot_band_diagram(solver, ks, p, bands=1:3, frequency_scale=1/2pi)
grid()

#compute transmission and reflection
kym = (-1:0.1:1)*0.5
omega0 = 0.005:0.005:0.3
RT = zeros(ComplexF64, length(kym), length(omega0))
TT = zeros(ComplexF64, length(kym), length(omega0))
jj = 1
for ky in kym
    ii = 1
    for omega in omega0
        kT, dT, VT = solver_EPWE(solver, ky, omega, p)
        RT[jj, ii], TT[jj, ii], _, _ = abcd(kT, dT, VT, omega, ky)
        ii = ii + 1
    end
    jj = jj +1
end

kymm = repeat(kym, 1, length(omega0))
omega00 = repeat(omega0, 1, length(kym))

pcolor(kymm, transpose(omega00),abs.(RT))


text:
#
# A simple toy atmospheric model that uses upwinded finite differences
# to simulate transport of a puff of non-reactive tracer gas.
#
# I wrote this to teach myself Juila.
# It's based on what I can remember of my PhD disseration from 12 years ago.
# This is absolutely _not_ how you go about solving this problem in the real world.
# I've no doubt there are much better packages out there for CFD.  Use those.
#
# Visualization is based on code from
# https://lazarusa.github.io/BeautifulMakie/surfWireLines/volumeTransparent/
#
# John Linford <john@redhpc.com>
#
# BSD 3-Clause License
#

using Dates
using GLMakie

"""
    advec_diff(...)

Advection and diffusion for one grid cell via upwind-biased finite differences.
Concentrations and wind values for the two neighboring cells are required.
"""
function advec_diff(
            Δx::Float64,
            c2l::Float64,
            c1l::Float64, w1l::Float64,
            c00::Float64, w00::Float64, d00::Float64,
            c1r::Float64, w1r::Float64,
            c2r::Float64)
    wind = (w1l + w00) / 2.0
    if wind > 0
        advecL = 1/6 * ( -c2l + 5.0*c1l + 2.0*c00 )
    else
        advecL = 1/6 * ( 2.0*c1l + 5.0*c00 - c1r )
    end
    advecL *= wind

    wind = (w1r + w00) / 2.0
    if wind > 0.0
        advecR = 1/6 * ( -c1l + 5.0*c00 + 2.0*c1r )
    else
        advecR = 1/6 * ( 2.0*c00 + 5.0*c1r - c2r )
    end
    advecR *= wind

    advec = (advecL - advecR) / Δx

    diffL = d00 * (c1l - c00)
    diffR = d00 * (c00 - c1r)
    diff = (diffL  -  diffR) / Δx^2

    return advec + diff
end

"""
    discretize1d(...)

Calculates advection/diffusion in 1D via `advec_diff` with boundaries wrapped around.

# Arguments
 - `Δx::Float64`: Space delta
 - `Δt::Float64`: Time delta
 - `d::Float64`: Diffusion coefficient
 - `c::Vector{Float64}`: Concentrations in ea
 - `w::Vector{Float64})`: Wind coefficients for each concentration
"""
function discretize1d(
            Δx::Float64,
            Δt::Float64,
            d::Float64,
            c::Vector{Float64},
            w::Vector{Float64})

    n = length(c)
    dcdx = Vector{Float64}(undef, n)

    # Boundaries
    dcdx[1] = advec_diff(Δx,
                         c[n-1],        # 2-left neighbor, wrapped around
                         c[n], w[n],    # 1-left neighbor, wrapped around
                         c[1], w[1], d, # Origin cell
                         c[2], w[2],    # 1-right neighbor
                         c[3])          # 2-right neighbor
    dcdx[2] = advec_diff(Δx,
                         c[n],          # 2-left neighbor, wrapped around
                         c[1], w[1],    # 1-left neighbor
                         c[2], w[2], d, # Origin cell
                         c[3], w[3],    # 1-right neighbor
                         c[4])          # 2-right neighbor

    # Interior
    for i in 3:n-2
        dcdx[i] = advec_diff(Δx,
                             c[i-2],            # 2-left neighbor
                             c[i-1], w[i-1],    # 1-left neighbor
                             c[i], w[i], d,     # Origin cell
                             c[i+1], w[i+1],    # 1-right neighbor
                             c[i+2])            # 2-right neighbor
    end

    # Boundaries
    dcdx[n-1] = advec_diff(Δx,
                           c[n-3],              # 2-left neighbor
                           c[n-2], w[n-2],      # 1-left neighbor
                           c[n-1], w[n-1], d,   # Origin cell
                           c[n], w[n],          # 1-right neighbor
                           c[1])                # 2-right neighbor, wrapped around
    dcdx[n] = advec_diff(Δx,
                         c[n-2],            # 2-left neighbor
                         c[n-1], w[n-1],    # 1-left neighbor
                         c[n], w[n], d,     # Origin cell
                         c[1], w[1],        # 1-right neighbor, wrapped around
                         c[2])              # 2-right neighbor, wrapped around

    c1 = c + Δt*dcdx
    # Mass loss here, if you care about that sort of thing
    return (c1 .* (c1 .> 0))
end

"""
    transport(...)

Discretization by first order upwind-biased finite differences on a possibly non-uniform grid.
3D discretization is achived via a 1D dimension split: x -> y -> z -> y -> x

# Arguments
 - `Δx`, `Δy`, `Δz`: Grid cell size in `x`, `y`, and `z` dimensions
 - `Δt`: timestep size
 - `dH`: Diffusion coefficient in x and y dimensions
 - `dV`: Diffusion coefficient in z dimension
 - `conc`: Concentration field
 - `wind`: Wind field
"""
function transport(
            Δx::Float64,
            Δy::Float64,
            Δz::Float64,
            Δt::Float64,
            dH::Float64,
            dV::Float64,
            conc::Array{Float64,4},
            wind::Array{Float64,4})

    nx, ny, nz, nspec = size(conc)

    for iy ∈ 1:ny, iz ∈ 1:nz, is ∈ 1:nspec
        conc[:,iy,iz,is] = discretize1d(Δx, Δt/2, dH, conc[:,iy,iz,is], wind[:,iy,iz,1])
    end
    for ix ∈ 1:nx, iz ∈ 1:nz, is ∈ 1:nspec
        conc[ix,:,iz,is] = discretize1d(Δy, Δt/2, dH, conc[ix,:,iz,is], wind[ix,:,iz,2])
    end
    for ix ∈ 1:nx, iy ∈ 1:ny, is ∈ 1:nspec
        conc[ix,iy,:,is] = discretize1d(Δz, Δt, dV, conc[ix,iy,:,is], wind[ix,iy,:,3])
    end
    for ix ∈ 1:nx, iz ∈ 1:nz, is ∈ 1:nspec
        conc[ix,:,iz,is] = discretize1d(Δy, Δt/2, dH, conc[ix,:,iz,is], wind[ix,:,iz,2])
    end
    for iy ∈ 1:ny, iz ∈ 1:nz, is ∈ 1:nspec
        conc[:,iy,iz,is] = discretize1d(Δx, Δt/2, dH, conc[:,iy,iz,is], wind[:,iy,iz,1])
    end
end

"""
    update_plot!(...)

Update Observables to advance the visualization.
"""
function update_plot!(points, colors, conc, spec, maxconc=8.61e9, minconc=1e-6)
    ptmp = Vector{Point3f0}()
    ctmp = Vector{Float64}()
    nx, ny, nz, _ = size(conc)
    for ix in 1:nx, iy in 1:ny, iz in 1:nz
        if conc[ix,iy,iz,spec] > minconc
            push!(ptmp, Point3f0(ix,iy,iz))
            push!(ctmp, conc[ix,iy,iz,spec]/maxconc)
        end
    end
    # Holy perf pain batman!  Don't update those observables until you absolutely have to!
    # One update outside the loop is at least 100x faster than updating in place.
    points[] = ptmp
    colors[] = ctmp
end

# Probably should be a main() function or something....
let
    # Time grid
    tstart = Dates.DateTime(2020, 01, 01)
    tend = tstart + Dates.Hour(24)
    Δt = Dates.Minute(15)

    # Space grid
    nx, ny, nz = 50, 50, 25

    # Grid cell dimensions (m)
    Δx, Δy, Δz = 1000.0, 1000.0, 1000.0

    # Horizontal and vertical diffusion coefficients (m²/s)
    dH, dV = 100.0, 50.0

    # Number of chemical species
    nspec = 1

    # Concentration field
    # Concentrations of the same species are contiguous
    conc = Array{Float64,4}(undef, nx, ny, nz, nspec)
    let
        conc .= 0.0
        # Dummy up a tracer point source in the center of the concentration field
        source_coords = size(conc)[1:3] .÷ 2
        conc[source_coords...] = 8.61e9
    end

    # Wind field
    wind = Array{Float64,4}(undef, nx, ny, nz, 3)
    wind_init = "simple"
    if wind_init == "simple"
        # 0.5m/sec cross wind in x-dimension
        wind[:,:,:,1] .= 0.5
        wind[:,:,:,2] .= 0.0
        wind[:,:,:,3] .= 0.0
    elseif wind_init == "cyclone"
        # Cyclonic wind
        for ix ∈ 1:nx, iy ∈ 1:ny, iz ∈ 1:nz
            wind[ix,iy,iz,:] = 0.5 .* [cos(nx/2 - ix), sin(ny/2 - iy), 0.25]
        end
    else
        error("Invalid wind_init: $wind_init")
    end

    # Temperature field
    # 297.15 °K == 24 °C == 75.2 °F
    temp = Array{Float64,3}(undef, nx, ny, nz)
    temp .= 297.15

    # Visualization
    fig = Figure(resolution=(1200,1200))
    ax = Axis3(fig[1,1], perspectiveness=0.5, azimuth=π/3, elevation=π/6, aspect=:data)
    points = Node(Vector{Point3f0}())
    colors = Node(Vector{Float64}())

    # Initial plot
    limits!(ax, 1, nx, 1, ny, 1, nz)
    update_plot!(points, colors, conc, 1)
    meshscatter!(ax, points, color=colors, colormap=RGBAf0.(to_colormap(:turbo, 10), LinRange(0,1,10)))

    # Time iteration with in-situ visualization
    record(fig, "plume.mp4", tstart:Δt:tend-Δt) do tnow
        println(tnow)
        transport(Δx, Δy, Δz, Float64(Dates.value(Dates.Second(Δt))), dH, dV, conc, wind)
        update_plot!(points, colors, conc, 1)
    end

end #let

max_stars_repo_name: jlinford/plume

max_issues_repo_name: jlinford/plume

max_forks_repo_name: jlinford/plume