text:
module Spatials

using LinearAlgebra: norm, dot, cross, inv
using StaticArrays: SVector, SMatrix
using Printf: @printf
using NearestNeighbors: KDTree, knn, inrange
using Base.Iterators: flatten, product
using ...Prerequisites: atol, rtol, Float
using ...Prerequisites.TypeTraits: efficientoperations
using ...Prerequisites.SimpleTrees: SimpleTree, simpletreedepth, isleaf
using ...Mathematics.Combinatorics: Combinations
using ...Mathematics.AlgebraOverFields: SimpleID

import ...Interfaces: decompose, rank, dimension, kind, expand, reset!

export distance, azimuthd, azimuth, polard, polar, volume
export isparallel, isonline, isintratriangle, issubordinate
export reciprocals, translate, rotate, tile, minimumlengths
export intralinks, interlinks
export PID, AbstractBond, pidtype, Point, Bond, rcoord, icoord, isintracell
export AbstractLattice, nneighbor, LatticeIndex, bonds!, bonds, LatticeBonds, Bonds, latticetype, bondtypes
export latticebondsstructure, allbonds, zerothbonds, insidebonds, acrossbonds, intrabonds, interbonds
export Lattice, SuperLattice, Cylinder

"""
    distance(p1::AbstractVector{<:Real}, p2::AbstractVector{<:Real}) -> Float

Get the distance between two points.

!!! note
    Compared to `norm(p1-p2)`, this function avoids the memory allocation for `p1-p2`, thus is more efficient.
"""
function distance(p1::AbstractVector{<:Real}, p2::AbstractVector{<:Real})
    @assert length(p1) == length(p2) "distance error: dismatched length of input vectors."
    result = 0
    for i = 1:length(p1)
        result = result + (p1[i]-p2[i])^2
    end
    sqrt(result)
end

"""
    azimuthd(v::AbstractVector{<:Real}) -> Float

Get the azimuth angle in degrees of a vector.
"""
function azimuthd(v::AbstractVector{<:Real})
    @assert length(v) ∈ (1, 2, 3) "azimuthd error: wrong dimensioned input vector."
    result = acosd(v[1]/((length(v) == 3) ? sqrt(v[1]^2+v[2]^2) : norm(v)))
    (length(v) > 1) && (v[2] < 0) && (result = 360 - result)
    return result
end

"""
    azimuth(v::AbstractVector{<:Real}) -> Float

Get the azimuth angle in radians of a vector.
"""
function azimuth(v::AbstractVector{<:Real})
    @assert length(v) ∈ (1, 2, 3) "azimuth error: wrong dimensioned input vector."
    result = acos(v[1]/((length(v) == 3) ? sqrt(v[1]^2+v[2]^2) : norm(v)))
    (length(v) > 1) && (v[2] < 0) && (result = 2pi - result)
    return result
end

"""
    polard(v::AbstractVector{<:Real}) -> Float

Get the polar angle in degrees of a vector.
"""
function polard(v::AbstractVector{<:Real})
    @assert length(v) == 3 "polard error: wrong dimensioned input vector."
    return acosd(v[3]/norm(v))
end

"""
    polar(v::AbstractVector{<:Real}) -> Float

Get the polar angle in radians of a vector.
"""
function polar(v::AbstractVector{<:Real})
    @assert length(v) == 3 "polard error: wrong dimensioned input vector."
    return acos(v[3]/norm(v))
end

"""
    volume(v1::AbstractVector{<:Real}, v2::AbstractVector{<:Real}, v3::AbstractVector{<:Real}) -> Real

Get the volume spanned by three vectors.
"""
function volume(v1::AbstractVector{<:Real}, v2::AbstractVector{<:Real}, v3::AbstractVector{<:Real})
    if (length(v1) ∈ (1, 2)) || (length(v2) ∈ (1, 2)) || (length(v3) ∈ (1, 2))
        result = zero(v1|>eltype)
    elseif (length(v1) == 3) && (length(v2) == 3) && (length(v3) == 3)
        result = v1[1]*(v2[2]*v3[3]-v2[3]*v3[2]) + v1[2]*(v2[3]*v3[1]-v2[1]*v3[3]) + v1[3]*(v2[1]*v3[2]-v2[2]*v3[1])
    else
        error("volume error: wrong dimensioned input vectors.")
    end
    return result
end

"""
    isparallel(v1::AbstractVector{<:Real}, v2::AbstractVector{<:Real}; atol::Real=atol, rtol::Real=rtol) -> Int

Judge whether two vectors are parallel to each other with the given tolerance, `0` for not parallel, `1` for parallel and `-1` for antiparallel.
"""
function isparallel(v1::AbstractVector{<:Real}, v2::AbstractVector{<:Real}; atol::Real=atol, rtol::Real=rtol)
    norm1, norm2 = norm(v1), norm(v2)
    if isapprox(norm1, 0.0, atol=atol, rtol=rtol) || isapprox(norm2, 0.0, atol=atol, rtol=rtol)
        result = 1
    elseif length(v1) == length(v2)
        temp = dot(v1, v2) / norm1 / norm2
        result = isapprox(temp, 1, atol=atol, rtol=rtol) ? 1 : isapprox(temp, -1, atol=atol, rtol=rtol) ? -1 : 0
    else
        error("isparallel error: shape dismatch of the input vectors.")
    end
    return result
end

"""
    isonline(p::AbstractVector{<:Real}, p1::AbstractVector{<:Real}, p2::AbstractVector{<:Real}; ends::Tuple{Bool, Bool}=(true, true), atol::Real=atol, rtol::Real=rtol) -> Bool

Judge whether a point is on a line segment whose end points are `p1` and `p2` with the given tolerance. `ends` defines whether the line segment should contain its ends.
"""
function isonline(p::AbstractVector{<:Real}, p1::AbstractVector{<:Real}, p2::AbstractVector{<:Real}; ends::Tuple{Bool, Bool}=(true, true), atol::Real=atol, rtol::Real=rtol)
    @assert length(p) == length(p1) == length(p2) "isonline error: shape dismatch of input point and line segment."
    d1, d2, d = distance(p, p1), distance(p, p2), distance(p1, p2)
    isapprox(d1, 0.0, atol=atol, rtol=rtol) && return ends[1]
    isapprox(d2, 0.0, atol=atol, rtol=rtol) && return ends[2]
    return isapprox(d1+d2, d, atol=atol, rtol=rtol)
end

"""
    decompose(v0::AbstractVector{<:Real}, v1::AbstractVector{<:Real}) -> Tuple{Float}
    decompose(v0::AbstractVector{<:Real}, v1::AbstractVector{<:Real}, v2::AbstractVector{<:Real}) -> Tuple{Float, Float}
    decompose(v0::AbstractVector{<:Real}, v1::AbstractVector{<:Real}, v2::AbstractVector{<:Real}, v3::AbstractVector{<:Real}) -> Tuple{Float, Float, Float}

Decompose a vector with respect to input basis vectors.
"""
function decompose(v0::AbstractVector{<:Real}, v1::AbstractVector{<:Real})
    @assert length(v0) == length(v1) "decompose error: dismatched length of input vectors."
    n0, n1 = norm(v0), norm(v1)
    sign = dot(v0, v1) / n0 / n1
    @assert isapprox(abs(sign), 1.0, atol=atol, rtol=rtol) "decompose error: insufficient basis vectors."
    return (sign*n0/n1,)
end
function decompose(v0::AbstractVector{<:Real}, v1::AbstractVector{<:Real}, v2::AbstractVector{<:Real})
    @assert length(v0) == length(v1) == length(v2) "decompose error: dismatched length of input vectors."
    @assert (length(v0) == 2) || (length(v0) == 3) "decompose error: unsupported dimension($(length(v0))) of input vectors."
    if length(v0) == 2
        det = v1[1]*v2[2] - v1[2]*v2[1]
        x1 = (v0[1]*v2[2]-v0[2]*v2[1]) / det
        x2 = (v0[2]*v1[1]-v0[1]*v1[2]) / det
    else
        v3 = SVector{3, Float}(v1[2]*v2[3]-v1[3]*v2[2], v1[3]*v2[1]-v1[1]*v2[3], v1[1]*v2[2]-v1[2]*v2[1])
        x1, x2, x3 = decompose(v0, v1, v2, v3)
        @assert isapprox(x3, 0.0, atol=atol, rtol=rtol) "decompose error: insufficient basis vectors."
    end
    return x1, x2
end
function decompose(v0::AbstractVector{<:Real}, v1::AbstractVector{<:Real}, v2::AbstractVector{<:Real}, v3::AbstractVector{<:Real})
    @assert length(v0) == length(v1) == length(v2) == length(v3) "decompose error: dismatched length of input vectors."
    @assert length(v0) == 3 "decompose error: unsupported dimension($(length(v0))) of input vectors."
    V = volume(v1, v2, v3)
    r1 = (v2[2]*v3[3]/V-v2[3]*v3[2]/V, v2[3]*v3[1]/V-v2[1]*v3[3]/V, v2[1]*v3[2]/V-v2[2]*v3[1]/V)
    r2 = (v3[2]*v1[3]/V-v3[3]*v1[2]/V, v3[3]*v1[1]/V-v3[1]*v1[3]/V, v3[1]*v1[2]/V-v3[2]*v1[1]/V)
    r3 = (v1[2]*v2[3]/V-v1[3]*v2[2]/V, v1[3]*v2[1]/V-v1[1]*v2[3]/V, v1[1]*v2[2]/V-v1[2]*v2[1]/V)
    return dot(r1, v0), dot(r2, v0), dot(r3, v0)
end

"""
    isintratriangle(p::AbstractVector{<:Real},
                    p1::AbstractVector{<:Real},
                    p2::AbstractVector{<:Real},
                    p3::AbstractVector{<:Real};
                    vertexes::NTuple{3, Bool}=(true, true, true),
                    edges::NTuple{3, Bool}=(true, true, true),
                    atol::Real=atol,
                    rtol::Real=rtol
                    ) -> Bool

Judge whether a point belongs to the interior of a triangle whose vertexes are `p1`, 'p2' and `p3` with the give tolerance. `vertexes` and `edges` define whether the interior should contain the vertexes or edges, respectively.
!!! note
    1. The vertexes are in the order (p1, p2, p3) and the edges are in the order (p1p2, p2p3, p3p1).
    2. The edges do not contain the vertexes.
"""
function isintratriangle(
        p::AbstractVector{<:Real},
        p1::AbstractVector{<:Real},
        p2::AbstractVector{<:Real},
        p3::AbstractVector{<:Real};
        vertexes::NTuple{3, Bool}=(true, true, true),
        edges::NTuple{3, Bool}=(true, true, true),
        atol::Real=atol,
        rtol::Real=rtol
        )
    @assert length(p) == length(p1) == length(p2) == length(p3) "isintratriangle error: shape dismatch of input point and triangle."
    @assert (length(p) == 2) || (length(p) == 3) "isintratriangle error: unsupported dimension($(length(p))) of input points."
    x = if length(p) == 2
        decompose(SVector(p[1]-p1[1], p[2]-p1[2]), SVector(p2[1]-p1[1], p2[2]-p1[2]), SVector(p3[1]-p1[1], p3[2]-p1[2]))
    else
        decompose(SVector(p[1]-p1[1], p[2]-p1[2], p[3]-p1[3]), SVector(p2[1]-p1[1], p2[2]-p1[2], p2[3]-p1[3]), SVector(p3[1]-p1[1], p3[2]-p1[2], p3[3]-p1[3]))
    end
    x1_approx_0, x2_approx_0 = isapprox(x[1], 0.0, atol=atol, rtol=rtol), isapprox(x[2], 0.0, atol=atol, rtol=rtol)
    x1_approx_1, x2_approx_1 = isapprox(x[1], 1.0, atol=atol, rtol=rtol), isapprox(x[2], 1.0, atol=atol, rtol=rtol)
    x12_approx_1 = isapprox(x[1]+x[2], 1.0, atol=atol, rtol=rtol)
    x1_approx_0 && x2_approx_0 && return vertexes[1]
    x1_approx_1 && x2_approx_0 && return vertexes[2]
    x1_approx_0 && x2_approx_1 && return vertexes[3]
    x2_approx_0 && (0 < x[1] < 1) && return edges[1]
    x12_approx_1 && (0 < x[1] < 1) && (0 < x[2] < 1) && return edges[2]
    x1_approx_0 && (0 < x[2] < 1) && return edges[3]
    return (0 < x[1] < 1) && (0 < x[2] < 1) && (x[1]+x[2] < 1)
end

"""
    issubordinate(rcoord::AbstractVector{<:Real}, vectors::AbstractVector{<:AbstractVector{<:Real}}; atol::Real=atol, rtol::Real=rtol) -> Bool

Judge whether a coordinate belongs to a lattice defined by `vectors` with the given tolerance.
"""
function issubordinate(rcoord::AbstractVector{<:Real}, vectors::AbstractVector{<:AbstractVector{<:Real}}; atol::Real=atol, rtol::Real=rtol)
    @assert length(rcoord) ∈ (1, 2, 3) "issubordinate error: only 1, 2 and 3 dimensional coordinates are supported."
    @assert length(vectors) ∈ (1, 2, 3) "issubordinate error: the number of input basis vectors must be 1, 2 or 3."
    fapprox = xi->isapprox(round(xi), xi, atol=atol, rtol=rtol)
    if length(vectors) == 1
        result = mapreduce(fapprox, &, decompose(rcoord, vectors[1]))
    elseif length(vectors) == 2
        result = mapreduce(fapprox, &, decompose(rcoord, vectors[1], vectors[2]))
    else
        result = mapreduce(fapprox, &, decompose(rcoord, vectors[1], vectors[2], vectors[3]))
    end
    return result
end

"""
    reciprocals(vectors::AbstractVector{AbstractVector{<:Real}}) -> Vector{Vector{Float}}

Get the reciprocals dual to the input vectors.
"""
function reciprocals(vectors::AbstractVector{<:AbstractVector{<:Real}})
    @assert length(vectors) < 4 "reciprocals error: the number of input vectors should not be greater than 3."
    @assert mapreduce(v->(length(v) ∈ (1, 2, 3)), &, vectors, init=true) "reciprocals error: all input vectors must be 1, 2 or 3 dimensional."
    result = Vector{Float}[]
    if length(vectors) == 1
        push!(result, 2pi/mapreduce(vi->vi^2, +, vectors[1])*vectors[1])
    elseif length(vectors) == 2
        v1, v2 = vectors[1], vectors[2]
        @assert length(v1) == length(v2) "reciprocals error: dismatched length of input vectors."
        if length(v1) == 2
            det = 2pi / (v1[1]*v2[2]-v1[2]*v2[1])
            push!(result, [det*v2[2], -det*v2[1]])
            push!(result, [-det*v1[2], det*v1[1]])
        else
            v3 = SVector{3, Float}(v1[2]*v2[3]-v1[3]*v2[2], v1[3]*v2[1]-v1[1]*v2[3], v1[1]*v2[2]-v1[2]*v2[1])
            V = 2pi / volume(v1, v2, v3)
            push!(result, [v2[2]*v3[3]*V-v2[3]*v3[2]*V, v2[3]*v3[1]*V-v2[1]*v3[3]*V, v2[1]*v3[2]*V-v2[2]*v3[1]*V])
            push!(result, [v3[2]*v1[3]*V-v3[3]*v1[2]*V, v3[3]*v1[1]*V-v3[1]*v1[3]*V, v3[1]*v1[2]*V-v3[2]*v1[1]*V])
        end
    elseif length(vectors) == 3
        v1, v2, v3 = vectors
        V = 2pi / volume(v1, v2, v3)
        push!(result, [v2[2]*v3[3]*V-v2[3]*v3[2]*V, v2[3]*v3[1]*V-v2[1]*v3[3]*V, v2[1]*v3[2]*V-v2[2]*v3[1]*V])
        push!(result, [v3[2]*v1[3]*V-v3[3]*v1[2]*V, v3[3]*v1[1]*V-v3[1]*v1[3]*V, v3[1]*v1[2]*V-v3[2]*v1[1]*V])
        push!(result, [v1[2]*v2[3]*V-v1[3]*v2[2]*V, v1[3]*v2[1]*V-v1[1]*v2[3]*V, v1[1]*v2[2]*V-v1[2]*v2[1]*V])
    end
    return result
end

"""
    translate(cluster::AbstractMatrix{<:Real}, vector::AbstractVector{<:Real}) -> Matrix{vector|>eltype}

Get the translated cluster of the original one by a vector.
"""
translate(cluster::AbstractMatrix{<:Real}, vector::AbstractVector{<:Real}) = cluster .+ reshape(vector, (vector|>length, 1))

"""
    rotate(cluster::AbstractMatrix{<:Real}, angle::Real; axis::Tuple{Union{AbstractVector{<:Real}, Nothing}, Tuple{<:Real, <:Real}}=(nothing, (0, 0))) -> Matrix{Float}

Get a rotated cluster of the original one by a certain angle around an axis.

The axis is determined by a point it gets through (`nothing` can be used to denote the origin), and its polar as well as azimuth angles in radians. The default axis is the z axis.
!!! note
    1. The result is given by the [Rodrigues' rotation formula](https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula).
    2. Only 2 and 3 dimensional vectors can be rotated.
    3. When the input vectors are 2 dimensional, both the polar and azimuth of the axis must be 0.
"""
function rotate(cluster::AbstractMatrix{<:Real}, angle::Real; axis::Tuple{Union{AbstractVector{<:Real}, Nothing}, Tuple{<:Real, <:Real}}=(nothing, (0, 0)))
    @assert size(cluster, 1) ∈ (2, 3) "rotate error: only 2 and 3 dimensional vectors can be rotated."
    center, theta, phi = ((axis[1] === nothing) ? zeros(size(cluster, 1)) : axis[1]), axis[2][1], axis[2][2]
    @assert length(center) == size(cluster, 1) "rotate error: dismatched shape of the input cluster and the point on axis."
    if length(center) == 2
        @assert isapprox(theta, 0, atol=atol) && isapprox(phi, 0, atol=atol) "rotate error: both the polar and azimuth of the axis for 2d vectors must be 0."
    end
    cosθ, sinθ = cos(angle), sin(angle)
    k, w = SVector{3, Float}(sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)), zeros(3)
    result = zeros(size(cluster))
    for i = 1:size(cluster, 2)
        for j = 1:size(cluster, 1)
            w[j] = cluster[j, i] - center[j]
        end
        inner = dot(k, w)
        outer = (k[2]*w[3]-k[3]*w[2], k[3]*w[1]-k[1]*w[3], k[1]*w[2]-k[2]*w[1])
        for j = 1:size(cluster, 1)
            result[j, i] = w[j]*cosθ + outer[j]*sinθ + k[j]*inner*(1-cosθ) + center[j]
        end
    end
    return result
end

"""
    tile(cluster::AbstractMatrix{<:Real}, vectors::AbstractVector{<:AbstractVector{<:Real}}, translations::NTuple{M, NTuple{N, <:Real}}=()) where {N, M} -> Matrix{Float}

Tile a supercluster by translations of the input cluster.

Basically, the final supercluster is composed of several parts, each of which is a translation of the original cluster, with the translation vectors specified by `vectors` and each set of the translation indices contained in `translations`. When translation vectors are empty, a copy of the original cluster will be returned.
"""
function tile(cluster::AbstractMatrix{<:Real}, vectors::AbstractVector{<:AbstractVector{<:Real}}, translations::NTuple{M, NTuple{N, <:Real}}=()) where {N, M}
    (N == 0) && return copy(cluster)
    @assert length(vectors) == N "tile error: dismatched shape of input vectors and translations."
    supercluster = zeros(size(cluster, 1), size(cluster, 2)*length(translations))
    disp = zeros(size(cluster, 1))
    for (i, translation) in enumerate(translations)
        for i = 1:length(disp)
            disp[i] = 0.0
            for j = 1:length(vectors)
                disp[i] += vectors[j][i] * translation[j]
            end
        end
        for j = 1:size(cluster, 2)
            col = (i-1)*size(cluster, 2) + j
            for row = 1:size(cluster, 1)
                supercluster[row, col] = cluster[row, j] + disp[row]
            end
        end
    end
    return supercluster
end

"""
    minimumlengths(cluster::AbstractMatrix{<:Real}, vectors::AbstractVector{<:AbstractVector{<:Real}}, nneighbor::Int=1; coordination::Int=8) -> Vector{Float}

Use kdtree to search the lowest several minimum bond lengths within a lattice translated by a cluster.

When the translation vectors are not empty, the lattice will be considered periodic in the corresponding directions. Otherwise the lattice will be open in all directions. To search for the bonds accorss the periodic boundaries, the cluster will be pretranslated to become a supercluster, which has open boundaries but is large enough to contain all the nearest neighbors within the required order. The `coordination` parameter sets the average number of each order of nearest neighbors. If it is to small, larger bond lengths may not be searched, and the result will contain `Inf`. This is a sign that you may need a larger `coordination`. Another situation that `Inf` appears in the result occurs when the minimum lengths are searched in open lattices. Indeed, the cluster may be too small so that the required order just goes beyond it. In this case the warning message can be safely ignored.
"""
function minimumlengths(cluster::AbstractMatrix{<:Real}, vectors::AbstractVector{<:AbstractVector{<:Real}}, nneighbor::Int=1; coordination::Int=8)
    @assert nneighbor >= 0 "minimumlengths error: input nneighbor must be non negative."
    result = [Inf for i = 1:nneighbor]
    if size(cluster, 2) > 0
        translations = reshape(product((-nneighbor:nneighbor for i = 1:length(vectors))...)|>collect, :)
        for translation in translations
            if (length(translation) > 0) && mapreduce(x->(x != 0), |, translation)
                remove = ntuple(i->-translation[i], length(translation))
                filter!(key->(key != remove), translations)
            end
        end
        supercluster = tile(cluster, vectors, Tuple(translations))
        for len in flatten(knn(KDTree(supercluster), cluster, (nneighbor > 0) ? min(nneighbor*coordination, size(supercluster, 2)) : 1, true)[2])
            for (i, minlen) in enumerate(result)
                if isapprox(len, minlen, atol=atol)
                    break
                elseif 0.0 < len < minlen
                    (nneighbor > 0) && (result[i+1:nneighbor] = result[i:nneighbor-1])
                    result[i] = len
                    break
                end
            end
        end
        if (nneighbor > 0) && mapreduce(isequal(Inf), |, result)
            @warn "minimumlengths warning: larger(>$coordination) coordination or smaller(<$nneighbor) nneighbor may be needed."
        end
    end
    return result
end

"""
    intralinks( cluster::AbstractMatrix{<:Real},
                vectors::AbstractVector{<:AbstractVector{<:Real}},
                neighbors::Dict{Int, Float},
                maxtranslations::NTuple{N, Int}=ntuple(i->length(neighbors), length(vectors))
                ) where N -> Vector{Tuple{Int, Int, Int, SubArray{<:Real, 1}}}

Use kdtree to get the intracluster nearest neighbors.

As is similar to [`minimumlengths`](@ref), when `vectors` is nonempty, the cluster assumes periodic boundaries. `neighbors` provides the map between the bond length and the order of nearest neighbors. Note only those with the lengths present in `neighbors` will be included in the result. `maxtranslations` determines the maximum number of translations along those directions specified by `vectors` when the tiled supercluster is construted (See [`minimumlengths`](@ref) for the explanation of the method for periodic lattices).
"""
function intralinks(cluster::AbstractMatrix{<:Real},
                    vectors::AbstractVector{<:AbstractVector{<:Real}},
                    neighbors::Dict{Int, Float},
                    maxtranslations::NTuple{N, Int}=ntuple(i->length(neighbors), length(vectors))
                    ) where N
    @assert length(vectors) == N "intralinks error: dismatched shape of input vectors and maxtranslations."
    result = Tuple{Int, Int, Int, Vector{Float}}[]
    (length(neighbors) == 0) && return result
    translations = reshape(product((-nnb:nnb for nnb in maxtranslations)...)|>collect, :)
    for translation in translations
        if (length(translation) > 0) && mapreduce(x->(x != 0), |, translation)
            remove = ntuple(i->-translation[i], length(translation))
            filter!(key->(key != remove), translations)
        end
    end
    sort!(translations, by=norm)
    translations = Tuple(translations)
    supercluster = tile(cluster, vectors, translations)
    disps = tile(zero(cluster), vectors, translations)
    for (i, indices) in enumerate(inrange(KDTree(supercluster), cluster, max(values(neighbors)...)+atol, true))
        for j in indices
            if i < j
                dist = 0.0
                for k = 1:size(cluster, 1)
                    dist = dist + (supercluster[k, j]-cluster[k, i])^2
                end
                dist = sqrt(dist)
                for (nb, len) in neighbors
                    if isapprox(len, dist, atol=atol)
                        push!(result, (nb, i, (j-1)%size(cluster, 2)+1, disps[:, j]))
                        break
                    end
                end
            end
        end
    end
    return result
end

"""
    interlinks(cluster1::AbstractMatrix{<:Real}, cluster2::AbstractMatrix{<:Real}, neighbors::Dict{Int, Float}) -> Vector{Tuple{Int, Int, Int}}

Use kdtree to get the intercluster nearest neighbors.
"""
function interlinks(cluster1::AbstractMatrix{<:Real}, cluster2::AbstractMatrix{<:Real}, neighbors::Dict{Int, Float})
    @assert size(cluster1, 1) == size(cluster2, 1) "interlinks error: dismatched space dimension of input clusters."
    result = Tuple{Int, Int, Int}[]
    (length(neighbors) == 0) && return result
    for (i, indices) in enumerate(inrange(KDTree(cluster2), cluster1, max(values(neighbors)...)+atol, true))
        for j in indices
            dist = 0.0
            for k = 1:size(cluster1, 1)
                dist = dist + (cluster2[k, j]-cluster1[k, i])^2
            end
            dist = sqrt(dist)
            for (nb, len) in neighbors
                if isapprox(len, dist, atol=atol)
                    push!(result, (nb, i, j))
                    break
                end
            end
        end
    end
    return result
end

"""
    PID(scope, site::Int)
    PID(site::Int)
    PID(; scope="tz", site::Int=1)

The id of a point.
"""
struct PID{S} <: SimpleID
    scope::S
    site::Int
end
PID(site::Int) = PID('T', site)
PID(; scope='T', site::Int=1) = PID(scope, site)
Base.fieldnames(pid::PID) = (:scope, :site)

"""
    AbstractBond{N, P<:PID, R}

Abstract bond.
"""
abstract type AbstractBond{N, P<:PID, R} end

"""
    ==(b1::AbstractBond{N, P, R}, b2::AbstractBond{N, P, R}) where {N, P<:PID, R} -> Bool

Overloaded equivalent operator.
"""
Base.:(==)(b1::AbstractBond{N, P, R}, b2::AbstractBond{N, P, R}) where {N, P<:PID, R} = ==(efficientoperations, b1, b2)

"""
    isequal(b1::AbstractBond{N, P, R}, b2::AbstractBond{N, P, R}) where {N, P<:PID, R} -> Bool

Overloaded equivalent function.
"""
Base.isequal(b1::AbstractBond{N, P, R}, b2::AbstractBond{N, P, R}) where {N, P<:PID, R} = isequal(efficientoperations, b1, b2)

"""
    length(bond::AbstractBond) -> Int
    length(::Type{<:AbstractBond{N, <:PID, R} where N}) where R -> Int

Get the number of points of a bond.
"""
Base.length(bond::AbstractBond) = bond |> typeof |> length
Base.length(::Type{<:AbstractBond{N, <:PID, R} where N}) where {R} = R

"""
    eltype(bond::AbstractBond)
    eltype(::Type{<:AbstractBond{N, P}}) where {N, P<:PID}

Get the eltype of a bond.
"""
Base.eltype(bond::AbstractBond) = bond |> typeof |> eltype
Base.eltype(::Type{<:AbstractBond{N, P}}) where {N, P<:PID} = Point{N, P}

"""
    dimension(bond::AbstractBond) -> Int
    dimension(::Type{<:AbstractBond{N}}) where N -> Int

Get the space dimension of a concrete bond.
"""
dimension(bond::AbstractBond) = bond |> typeof |> dimension
dimension(::Type{<:AbstractBond{N}}) where {N} = N

"""
    pidtype(bond::AbstractBond)
    pidtype(::Type{<:AbstractBond{N, P} where N}) where {P<:PID}

Get the pid type of a concrete bond.
"""
pidtype(bond::AbstractBond) = bond |> typeof |> pidtype
pidtype(::Type{<:AbstractBond{N, P} where N}) where {P<:PID} = P

"""
    rank(bond::AbstractBond) -> Int
    rank(::Type{<:AbstractBond{N, <:PID, R} where N}) where R -> Int

Get the rank of a bond.
"""
rank(bond::AbstractBond) = bond |> typeof |> rank
rank(::Type{<:AbstractBond{N, <:PID, R} where N}) where R = R

"""
    Point(pid::PID, rcoord::SVector{N, <:Real}, icoord::SVector{N, <:Real}) where N
    Point(pid::PID, rcoord::NTuple{N, <:Real}, icoord::NTuple{N, <:Real}=ntuple(i->0.0, N)) where N
    Point(pid::PID, rcoord::AbstractVector{<:Real}, icoord::AbstractVector{<:Real}=zero(SVector{length(rcoord), Float}))
    Point{N}(pid::PID, rcoord::AbstractVector{<:Real}, icoord::AbstractVector{<:Real}=zero(SVector{N, Float})) where N

Labeled point.
"""
struct Point{N, P<:PID} <: AbstractBond{N, P, 1}
    pid::P
    rcoord::SVector{N, Float}
    icoord::SVector{N, Float}
    Point(pid::PID, rcoord::SVector{N, <:Real}, icoord::SVector{N, <:Real}) where {N} = new{N, pid|>typeof}(pid, rcoord, icoord)
end
function Point(pid::PID, rcoord::NTuple{N, <:Real}, icoord::NTuple{N, <:Real}=ntuple(i->0.0, N)) where N
    Point(pid, convert(SVector{N, Float}, rcoord), convert(SVector{N, Float}, icoord))
end
function Point(pid::PID, rcoord::AbstractVector{<:Real}, icoord::AbstractVector{<:Real}=zero(SVector{length(rcoord), Float}))
    Point{length(rcoord)}(pid, rcoord, icoord)
end
function Point{N}(pid::PID, rcoord::AbstractVector{<:Real}, icoord::AbstractVector{<:Real}=zero(SVector{N, Float})) where N
    @assert length(rcoord) == length(icoord) == N "Point error: dismatched length of input rcoord and icoord."
    Point(pid, convert(SVector{N, Float}, rcoord), convert(SVector{N, Float}, icoord))
end

"""
    show(io::IO, p::Point)

Show a labeled point.
"""
Base.show(io::IO, p::Point) = @printf io "Point(%s, %s, %s)" p.pid p.rcoord p.icoord

"""
    iterate(p::Point, state=1)

Iterate over the point.
"""
Base.iterate(p::Point, state=1) = (state == 1) ? (p, state+1) : nothing

"""
    kind(::Point) -> 0
    kind(::Type{<:Point}) -> 0

Get the bond kind of a point, which is defined to be 0.
"""
kind(::Point) = 0
kind(::Type{<:Point}) = 0

"""
    Bond(neighbor::Int, spoint::Point, epoint::Point)

A bond in a lattice.
"""
struct Bond{N, P<:PID} <: AbstractBond{N, P, 2}
    neighbor::Int
    spoint::Point{N, P}
    epoint::Point{N, P}
end

"""
    show(io::IO, bond::Bond)

Show a bond.
"""
Base.show(io::IO, bond::Bond) = @printf io "Bond(%s, %s, %s)" bond.neighbor bond.spoint bond.epoint

"""
    reverse(bond::Bond) -> Bond

Get the reversed bond.
"""
Base.reverse(bond::Bond) = Bond(bond.neighbor, bond.epoint, bond.spoint)

"""
    iterate(bond::Bond, state=1)

Iterate over the points in a bond.
"""
Base.iterate(bond::Bond, state=1) = (state == 1) ? (bond.epoint, state+1) : (state == 2) ? (bond.spoint, state+1) : nothing

"""
    kind(bond::Bond) -> Int

Get the bond kind of a bond.
"""
kind(bond::Bond) = bond.neighbor

"""
    rcoord(bond::Bond) -> SVector

Get the rcoord of the bond.
"""
rcoord(bond::Bond) = bond.epoint.rcoord - bond.spoint.rcoord

"""
    icoord(bond::Bond) -> SVector

Get the icoord of the bond.
"""
icoord(bond::Bond) = bond.epoint.icoord - bond.spoint.icoord

"""
    isintracell(bond::Bond) -> Bool

Judge whether a bond is intra the unit cell of a lattice.
"""
isintracell(bond::Bond) = isapprox(bond|>icoord|>norm, 0.0, atol=atol)

"""
    AbstractLattice{P<:PID, N}

Abstract type for all lattices.

It should have the following attributes
- `name::String`: the name of the lattice
- `pids::Vector{P}`: the pids of the lattice
- `rcoords::Matrix{Float}`: the rcoords of the lattice
- `icoords::Matrix{Float}`: the icoords of the lattice
- `vectors::Vector{SVector{N, Float}}`: the translation vectors of the lattice
- `reciprocals::Vector{SVector{N, Float}}`: the reciprocals of the lattice
- `neighbors::Dict{Int, Float}`: the order-distance map of the nearest neighbors of the lattice
"""
abstract type AbstractLattice{N, P<:PID} end

"""
    show(io::IO, lattice::AbstractLattice)

Show a lattice.
"""
function Base.show(io::IO, lattice::AbstractLattice)
    @printf io "%s(%s)\n" lattice|>typeof|>nameof lattice.name
    if length(lattice) > 0
        @printf io "  with %s %s:\n" length(lattice) (length(lattice) == 1) ? "point" : "points"
        for i = 1:length(lattice)
            @printf io "    %s\n" lattice[LatticeIndex{'P'}(i)]
        end
    end
    if length(lattice.vectors) > 0
        @printf io "  with %s translation %s:\n" length(lattice.vectors) (length(lattice.vectors) == 1) ? "vector" : "vectors"
        for i = 1:length(lattice.vectors)
            #@printf io "    [%s]\n" join(lattice.vectors[i], ", ")
            @printf io "    %s\n" lattice.vectors[i]
        end
    end
    if length(lattice.neighbors) > 0
        @printf io "  with %s %s of nearest neighbors:\n" length(lattice.neighbors) (length(lattice.neighbors) == 1) ? "order" : "orders"
        for (order, distance) in lattice.neighbors
            @printf io "    %s => %s\n" order distance
        end
    end
end

"""
    length(lattice::AbstractLattice) -> Int

Get the number of points contained in a lattice.
"""
Base.length(lattice::AbstractLattice) = length(lattice.pids)

"""
    ==(lattice1::AbstractLattice, lattice2::AbstractLattice) -> Bool

Overloaded equivalent operator.
"""
Base.:(==)(lattice1::AbstractLattice, lattice2::AbstractLattice) = ==(efficientoperations, lattice1, lattice2)

"""
    isequal(lattice1::AbstractLattice, lattice2::AbstractLattice) -> Bool

Overloaded equivalent function.
"""
Base.isequal(lattice1::AbstractLattice, lattice2::AbstractLattice) = isequal(efficientoperations, lattice1, lattice2)

"""
    keytype(lattice::AbstractLattice)
    keytype(::Type{<:AbstractLattice{N, P} where N}) where {P<:PID}

Get the pid type of the lattice.
"""
Base.keytype(lattice::AbstractLattice) = lattice |> typeof |> keytype
Base.keytype(::Type{<:AbstractLattice{N, P} where N}) where {P<:PID} = P

"""
    valtype(lattice::AbstractLattice)
    valtype(::Type{<:AbstractLattice{N, P}}) where {N, P<:PID}

Get the point type of the lattice.
"""
Base.valtype(lattice::AbstractLattice) = lattice |> typeof |> valtype
Base.valtype(::Type{<:AbstractLattice{N, P}}) where {N, P<:PID} = Point{N, P}

"""
    dimension(lattice::AbstractLattice) -> Int
    dimension(::Type{<:AbstractLattice{N}}) where N -> Int

Get the space dimension of the lattice.
"""
dimension(lattice::AbstractLattice) = lattice |> typeof |> dimension
dimension(::Type{<:AbstractLattice{N}}) where {N} = N

"""
    nneighbor(lattice::AbstractLattice) -> Int

Get the highest order of nearest neighbors.
"""
nneighbor(lattice::AbstractLattice) = max(keys(lattice.neighbors)...)

"""
    LatticeIndex{Kind}(index::Union{PID, Int}) where Kind

Lattice index.

`Kind` must be one of the followings:
* 'R': for getting the rcoord of a lattice
* 'I': for getting the icoord of a lattice
* 'P': for getting the point of a lattice
"""
struct LatticeIndex{Kind, I<:Union{PID, Int}}
    index::I
    function LatticeIndex{Kind}(index::Union{PID, Int}) where Kind
        @assert Kind in ('R', 'I', 'P') "LatticeIndex error: wrong input Kind($Kind)."
        new{Kind, typeof(index)}(index)
    end
end

"""
    getindex(lattice::AbstractLattice, i::LatticeIndex{'R'}) -> SVector
    getindex(lattice::AbstractLattice, i::LatticeIndex{'I'}) -> SVector
    getindex(lattice::AbstractLattice, i::LatticeIndex{'P'}) -> Point

Get a rcoord, an icoord or a point of a lattice according to the type of the input index.
"""
Base.getindex(lattice::AbstractLattice, i::LatticeIndex{'R', Int}) = latticestaticcoords(lattice.rcoords, i.index, lattice|>dimension|>Val)
Base.getindex(lattice::AbstractLattice, i::LatticeIndex{'I', Int}) = latticestaticcoords(lattice.icoords, i.index, lattice|>dimension|>Val)
Base.getindex(lattice::AbstractLattice, i::LatticeIndex{'P', Int}) = Point(lattice.pids[i.index], lattice[LatticeIndex{'R'}(i.index)], lattice[LatticeIndex{'I'}(i.index)])
Base.getindex(lattice::AbstractLattice, i::LatticeIndex{'R', <:PID}) = lattice[LatticeIndex{'R'}(findfirst(isequal(i.index), lattice.pids))]
Base.getindex(lattice::AbstractLattice, i::LatticeIndex{'I', <:PID}) = lattice[LatticeIndex{'I'}(findfirst(isequal(i.index), lattice.pids))]
Base.getindex(lattice::AbstractLattice, i::LatticeIndex{'P', <:PID}) = lattice[LatticeIndex{'P'}(findfirst(isequal(i.index), lattice.pids))]
@generated function latticestaticcoords(coords::Matrix{Float}, i::Int, ::Val{N}) where N
    exprs = [:(coords[$j, i]) for j = 1:N]
    return :(SVector{N, Float}($(exprs...)))
end

abstract type LatticeBonds{R} end
Base.eltype(::Type{<:LatticeBonds{0}}) = AbstractBond
Base.eltype(::Type{<:LatticeBonds{1}}) = Point
Base.eltype(::Type{<:LatticeBonds{2}}) = Bond
Base.eltype(::Type{L}, ::B) where {L<:AbstractLattice, B<:LatticeBonds} = eltype(B){L|>dimension, L|>keytype}
Base.eltype(::Type{L}, ::Val{B}) where {L<:AbstractLattice, B} = eltype(typeof(B)){L|>dimension, L|>keytype}
struct AllBonds <: LatticeBonds{0} end
struct ZerothBonds <: LatticeBonds{1} end
struct InsideBonds <: LatticeBonds{2} end
struct AcrossBonds <: LatticeBonds{2} end
"""
    allbonds

Indicate that all bonds are inquired.
"""
const allbonds = AllBonds()
"""
    zerothbonds

Indicate that zeroth bonds, i.e. the points are inquired.
"""
const zerothbonds = ZerothBonds()
"""
    insidebonds

Indicate that bonds inside the unitcell are inquired, which do not contain those across the periodic boundaries.
"""
const insidebonds = InsideBonds()
"""
    acrossbonds

Indicate that bonds across the unitcell are inquired, which are in fact those across the periodic boundaries.
"""
const acrossbonds = AcrossBonds()

"""
    latticebondsstructure(::Type{<:AbstractLattice}) -> SimpleTree{LatticeBonds, Nothing}

The tree structure of the lattice bonds.
"""
@generated function latticebondsstructure(::Type{<:AbstractLattice})
    structure = SimpleTree{LatticeBonds, Nothing}()
    push!(structure, allbonds, nothing)
    push!(structure, allbonds, zerothbonds, nothing)
    push!(structure, allbonds, insidebonds, nothing)
    push!(structure, allbonds, acrossbonds, nothing)
    return structure
end

"""
    expand(::Type{L}, ::Val{VS}) where {L<:AbstractLattice, VS} -> Tuple

Expand the lattice bond types to the leaf level.
"""
@generated function expand(::Type{L}, ::Val{VS}) where {L<:AbstractLattice, VS}
    leaves = LatticeBonds[]
    structure = latticebondsstructure(L)
    for V in VS
        for latticebonds in keys(structure, simpletreedepth, V)
            isleaf(structure, latticebonds) && push!(leaves, latticebonds)
        end
    end
    return Tuple(leaves)
end

"""
    bonds(lattice::AbstractLattice, inquiry=allbonds) -> Vector{eltype(lattice|>typeof, inquiry)}
    bonds(lattice::AbstractLattice, inquiries...) -> Vector{mapreduce(inquiry->eltype(lattice|>typeof, inquiry), typejoin, inquiries)}

Generate the required bonds of a lattice.
"""
bonds(lattice::AbstractLattice) = bonds(lattice, allbonds)
bonds(lattice::AbstractLattice, inquiry) = bonds!(eltype(lattice|>typeof, inquiry)[], lattice, inquiry)
bonds(lattice::AbstractLattice, inquiries...) = bonds!(mapreduce(inquiry->eltype(lattice|>typeof, inquiry), typejoin, inquiries)[], lattice, inquiries...)

"""
    bonds!(bonds::Vector, lattice::AbstractLattice, inquiries::LatticeBonds...) -> Vector
    bonds!(bonds::Vector, lattice::AbstractLattice, ::Val{zerothbonds}) -> Vector
    bonds!(bonds::Vector, lattice::AbstractLattice, ::Val{insidebonds}) -> Vector
    bonds!(bonds::Vector, lattice::AbstractLattice, ::Val{acrossbonds}) -> Vector

Generate the required bonds of a lattice and append them to the input bonds.
"""
bonds!(bonds::Vector, lattice::AbstractLattice) = bonds!(bonds, lattice, allbonds)
function bonds!(bonds::Vector, lattice::AbstractLattice, inquiries::LatticeBonds...)
    leaves = expand(typeof(lattice), Val(inquiries))
    for i = 1:length(leaves)
        bonds!(bonds, lattice, leaves[i]|>Val)
    end
    return bonds
end
function bonds!(bonds::Vector, lattice::AbstractLattice, ::Val{zerothbonds})
    for i = 1:length(lattice)
        push!(bonds, lattice[LatticeIndex{'P'}(i)])
    end
    return bonds
end
function bonds!(bonds::Vector, lattice::AbstractLattice, ::Val{insidebonds})
    for (neighbor, sindex, eindex) in interlinks(lattice.rcoords, lattice.rcoords, lattice.neighbors)
        if sindex < eindex
            spoint = lattice[LatticeIndex{'P'}(sindex)]
            epoint = lattice[LatticeIndex{'P'}(eindex)]
            push!(bonds, Bond(neighbor, spoint, epoint))
        end
    end
    return bonds
end
function bonds!(bonds::Vector, lattice::AbstractLattice, ::Val{acrossbonds})
    nnb = lattice |> nneighbor
    translations = reshape(product((-nnb:nnb for i = 1:length(lattice.vectors))...)|>collect, :)
    for translation in translations
        remove = ntuple(i->-translation[i], length(translation))
        filter!(key->(key != remove), translations)
    end
    if length(translations) > 0
        dim = lattice |> dimension |> Val
        superrcoords = tile(lattice.rcoords, lattice.vectors, Tuple(translations))
        supericoords = tile(lattice.icoords, lattice.vectors, Tuple(translations))
        for (neighbor, sindex, eindex) in interlinks(lattice.rcoords, superrcoords, lattice.neighbors)
            spoint = Point(lattice.pids[sindex], latticestaticcoords(lattice.rcoords, sindex, dim), latticestaticcoords(lattice.icoords, sindex, dim))
            epoint = Point(lattice.pids[(eindex-1)%size(lattice.rcoords, 2)+1], latticestaticcoords(superrcoords, eindex, dim), latticestaticcoords(supericoords, eindex, dim))
            push!(bonds, Bond(neighbor, spoint, epoint))
        end
    end
    return bonds
end

"""
    Lattice{N}( name::String,
                pids::Vector{<:PID},
                rcoords::AbstractMatrix{<:Real},
                icoords::AbstractMatrix{<:Real},
                vectors::AbstractVector{<:AbstractVector{<:Real}},
                neighbors::Union{Dict{Int, <:Real}, Int}=1;
                coordination::Int=8
                ) where N
    Lattice(    name::String,
                points::AbstractVector{<:Point};
                vectors::AbstractVector{<:AbstractVector{<:Real}}=SVector{0, SVector{points|>eltype|>dimension, Float}}(),
                neighbors::Union{Dict{Int, <:Real}, Int}=1,
                coordination::Int=8
                )
    Lattice(    name::String,
                sublattices::AbstractVector{<:AbstractLattice};
                vectors::AbstractVector{<:AbstractVector{<:Real}}=SVector{0, SVector{sublattices|>eltype|>dimension, Float}}(),
                neighbors::Union{Dict{Int, <:Real}, Int}=1,
                coordination::Int=8
                )

Simplest lattice.

A simplest lattice can be construted from its contents, i.e. pids, rcoords and icoords, or from a couple of points, or from a couple of sublattices.
"""
struct Lattice{N, P<:PID} <: AbstractLattice{N, P}
    name::String
    pids::Vector{P}
    rcoords::Matrix{Float}
    icoords::Matrix{Float}
    vectors::Vector{SVector{N, Float}}
    reciprocals::Vector{SVector{N, Float}}
    neighbors::Dict{Int, Float}
    function Lattice{N}(name::String,
                        pids::Vector{<:PID},
                        rcoords::AbstractMatrix{<:Real},
                        icoords::AbstractMatrix{<:Real},
                        vectors::AbstractVector{<:AbstractVector{<:Real}},
                        neighbors::Union{Dict{Int, <:Real}, Int}=1;
                        coordination::Int=8
                        ) where N
        @assert (N == size(rcoords, 1) == size(icoords, 1)) && (length(pids) == size(rcoords, 2) == size(icoords, 2)) "Lattice error: dismatched shape of input pids, rcoords and icoords."
        isa(neighbors, Int) && (neighbors = Dict(i=>minlen for (i, minlen) in enumerate(minimumlengths(rcoords, vectors, neighbors, coordination=coordination))))
        vectors = convert(Vector{SVector{N, Float}}, vectors)
        recipls = convert(Vector{SVector{N, Float}}, reciprocals(vectors))
        new{N, pids|>eltype}(name, pids, rcoords, icoords, vectors, recipls, neighbors)
    end
end
function Lattice(   name::String,
                    points::AbstractVector{<:Point};
                    vectors::AbstractVector{<:AbstractVector{<:Real}}=SVector{0, SVector{points|>eltype|>dimension, Float}}(),
                    neighbors::Union{Dict{Int, <:Real}, Int}=1,
                    coordination::Int=8
                    )
    pids = Vector{points|>eltype|>pidtype}(undef, points|>length)
    rcoords = zeros(Float, points|>eltype|>dimension, points|>length)
    icoords = zeros(Float, points|>eltype|>dimension, points|>length)
    for i = 1:length(points)
        pids[i] = points[i].pid
        rcoords[:, i] = points[i].rcoord
        icoords[:, i] = points[i].icoord
    end
    Lattice{points|>eltype|>dimension}(name, pids, rcoords, icoords, vectors, neighbors, coordination=coordination)
end
function Lattice(   name::String,
                    sublattices::AbstractVector{<:AbstractLattice};
                    vectors::AbstractVector{<:AbstractVector{<:Real}}=SVector{0, SVector{sublattices|>eltype|>dimension, Float}}(),
                    neighbors::Union{Dict{Int, <:Real}, Int}=1,
                    coordination::Int=8
                    )
    if length(sublattices) > 1
        @assert mapreduce(sl->(dimension(sl) == dimension(sublattices[1])), &, sublattices) "Lattice error: all sublattices should have the same dimension."
        @assert mapreduce(sl->(keytype(sl) === keytype(sublattices[1])), &, sublattices) "Lattice error: all sublattices should have the same keytype."
    end
    len = mapreduce(length, +, sublattices)
    pids = Vector{sublattices|>eltype|>keytype}(undef, len)
    rcoords = zeros(Float, sublattices|>eltype|>dimension, len)
    icoords = zeros(Float, sublattices|>eltype|>dimension, len)
    pos = 1
    for sublattice in sublattices
        pids[pos:pos+length(sublattice)-1] = sublattice.pids
        rcoords[:, pos:pos+length(sublattice)-1] = sublattice.rcoords
        icoords[:, pos:pos+length(sublattice)-1] = sublattice.icoords
        pos += length(sublattice)
    end
    Lattice{sublattices|>eltype|>dimension}(name, pids, rcoords, icoords, vectors, neighbors, coordination=coordination)
end

"""
    SuperLattice(   name::String,
                    sublattices::AbstractVector{<:AbstractLattice};
                    vectors::AbstractVector{<:AbstractVector{<:Real}}=SVector{0, SVector{sublattices|>eltype|>dimension, Float}}(),
                    neighbors::Dict{Int, <:Real}=Dict{Int, Float}()
                    )

SuperLattice that is composed of serveral sublattices.
"""
struct SuperLattice{L<:AbstractLattice, N, P<:PID} <: AbstractLattice{N, P}
    sublattices::Vector{L}
    name::String
    pids::Vector{P}
    rcoords::Matrix{Float}
    icoords::Matrix{Float}
    vectors::Vector{SVector{N, Float}}
    reciprocals::Vector{SVector{N, Float}}
    neighbors::Dict{Int, Float}
    function SuperLattice(  name::String,
                            sublattices::AbstractVector{<:AbstractLattice};
                            vectors::AbstractVector{<:AbstractVector{<:Real}}=SVector{0, SVector{sublattices|>eltype|>dimension, Float}}(),
                            neighbors::Dict{Int, <:Real}=Dict{Int, Float}()
                            )
        if length(sublattices) > 1
            @assert mapreduce(sl->(dimension(sl) == dimension(sublattices[1])), &, sublattices) "SuperLattice error: all sublattices should have the same dimension."
            @assert mapreduce(sl->(keytype(sl) === keytype(sublattices[1])), &, sublattices) "SuperLattice error: all sublattices should have the same keytype."
        end
        @assert all(length(sublattice.vectors) == 0 for sublattice in sublattices) "SuperLattice error: all sublattices should assume open boundaries."
        len = mapreduce(length, +, sublattices)
        pids = Vector{sublattices|>eltype|>keytype}(undef, len)
        rcoords = zeros(Float, sublattices|>eltype|>dimension, len)
        icoords = zeros(Float, sublattices|>eltype|>dimension, len)
        pos = 1
        for sublattice in sublattices
            pids[pos:pos+length(sublattice)-1] = sublattice.pids
            rcoords[:, pos:pos+length(sublattice)-1] = sublattice.rcoords
            icoords[:, pos:pos+length(sublattice)-1] = sublattice.icoords
            pos += length(sublattice)
        end
        vectors = convert(Vector{SVector{sublattices|>eltype|>dimension, Float}}, vectors)
        recipls = convert(Vector{SVector{sublattices|>eltype|>dimension, Float}}, reciprocals(vectors))
        new{sublattices|>eltype, sublattices|>eltype|>dimension, sublattices|>eltype|>keytype}(sublattices, name, pids, rcoords, icoords, vectors, recipls, neighbors)
    end
end

"""
    latticetype(sl::SuperLattice)
    latticetype(::Type{<:SuperLattice{L}}) where {L<:AbstractLattice}

Get the sublattice type of a superlattice.
"""
latticetype(sl::SuperLattice) = sl |> typeof |> latticetype
latticetype(::Type{<:SuperLattice{L}}) where {L<:AbstractLattice} = L

struct IntraBonds <: LatticeBonds{2} end
struct InterBonds <: LatticeBonds{2} end
"""
    intrabonds

Indicate that bonds intra the sublattices are inquired.
!!! note
    These bonds do not contain those accorss the periodic boundaries.
"""
const intrabonds = IntraBonds()
"""
    interbonds

Indicate that bonds inter the sublattices are inquired.
!!! note
    These bonds do not contain those accorss the periodic boundaries.
"""
const interbonds = InterBonds()

"""
    latticebondsstructure(::Type{<:SuperLattice}) -> SimpleTree{LatticeBonds, Nothing}

The tree structure of the lattice bonds.
"""
@generated function latticebondsstructure(::Type{<:SuperLattice})
    structure = SimpleTree{LatticeBonds, Nothing}()
    push!(structure, allbonds, nothing)
    push!(structure, allbonds, zerothbonds, nothing)
    push!(structure, allbonds, insidebonds, nothing)
    push!(structure, allbonds, acrossbonds, nothing)
    push!(structure, insidebonds, intrabonds, nothing)
    push!(structure, insidebonds, interbonds, nothing)
    return structure
end

"""
    bonds!(bonds::Vector, lattice::SuperLattice, ::Val{intrabonds}) -> Vector
    bonds!(bonds::Vector, lattice::SuperLattice, ::Val{interbonds}) -> Vector

Generate the required bonds of a superlattice and append them to the input bonds.
"""
function bonds!(bonds::Vector, lattice::SuperLattice, ::Val{intrabonds})
    for sublattice in lattice.sublattices
        bonds!(bonds, sublattice, insidebonds)
    end
    return bonds
end
function bonds!(bonds::Vector, lattice::SuperLattice, ::Val{interbonds})
    dim = lattice |> dimension |> Val
    for (i, j) in Combinations{2}(1:length(lattice.sublattices))
        sub1, sub2 = lattice.sublattices[i], lattice.sublattices[j]
        for (neighbor, sindex, eindex) in interlinks(sub1.rcoords, sub2.rcoords, lattice.neighbors)
            spoint = Point(sub1.pids[sindex], latticestaticcoords(sub1.rcoords, sindex, dim), latticestaticcoords(sub1.icoords, sindex, dim))
            epoint = Point(sub2.pids[eindex], latticestaticcoords(sub2.rcoords, eindex, dim), latticestaticcoords(sub2.icoords, eindex, dim))
            push!(bonds, Bond(neighbor, spoint, epoint))
        end
    end
    return bonds
end

"""
    Cylinder{P}(name::String,
                block::AbstractMatrix{<:Real},
                translation::SVector{N, <:Real};
                vector::Union{AbstractVector{<:Real}, Nothing}=nothing,
                neighbors::Union{Dict{Int, <:Real}, Int}=1
                ) where {P<:PID, N}

Cylinder of 1d and quasi 2d lattices.
"""
mutable struct Cylinder{P<:PID, N} <: AbstractLattice{N, P}
    block::Matrix{Float}
    translation::SVector{N, Float}
    name::String
    pids::Vector{P}
    rcoords::Matrix{Float}
    icoords::Matrix{Float}
    vectors::Vector{SVector{N, Float}}
    reciprocals::Vector{SVector{N, Float}}
    neighbors::Dict{Int, Float}
    function Cylinder{P}(   name::String,
                            block::AbstractMatrix{<:Real},
                            translation::SVector{N, <:Real};
                            vector::Union{AbstractVector{<:Real}, Nothing}=nothing,
                            neighbors::Union{Dict{Int, <:Real}, Int}=1
                            ) where {P<:PID, N}
        pids = Vector{P}[]
        rcoords = zeros(Float, N, 0)
        icoords = zeros(Float, N, 0)
        vectors = (vector === nothing) ? SVector{N, Float}[] : [convert(SVector{N, Float}, vector)]
        recipls = convert(Vector{SVector{N, Float}}, reciprocals(vectors))
        isa(neighbors, Int) && (neighbors = Dict{Int, Float}(i=>Inf for i = 1:neighbors))
        new{P, N}(block, translation, name, pids, rcoords, icoords, vectors, recipls, neighbors)
    end
end

"""
    insert!(cylinder::Cylinder, ps::S...; cut::Int=length(cylinder)÷2+1, scopes::Union{<:AbstractVector{S}, Nothing}=nothing, coordination::Int=9) where S -> Cylinder

Insert a couple of blocks into a cylinder.

The position of the cut of the cylinder is specified by the keyword argument `cut`, which is the center of the cylinder by default. All pids corresponding to a same newly inserted block share the same scope, which is specified by the parameter `ps`. Optionally, the scopes of the old pids in the cylinder can be replaced if the parameter `scopes` is assigned other than `nothing`. Note the length of `ps` is equal to the number of newly inserted blocks, while that of `scopes` should be equal to the old length of the cylinder.
"""
function Base.insert!(cylinder::Cylinder, ps::S...; cut::Int=length(cylinder)÷2+1, scopes::Union{<:AbstractVector{S}, Nothing}=nothing, coordination::Int=8) where S
    @assert S <: fieldtype(cylinder|>keytype, :scope) "insert! error: dismatched type of input scopes and old scopes."
    @assert 1 <= cut <= length(cylinder)+1 "insert! error: wrong cut($cut), which should be in [1, $(length(cylinder)+1)]."
    @assert length(cylinder)%size(cylinder.block, 2) == 0 "insert! error: wrong cylinder length."
    (scopes === nothing) || @assert length(scopes) == length(cylinder) "insert! error: dismatched length of input scopes and cylinder."
    blcklen = size(cylinder.block, 2)
    blcknum = length(cylinder)÷blcklen + length(ps)
    rsltlen = blcklen * blcknum
    pids = keytype(cylinder)[]
    rcoords = zeros(Float, cylinder|>dimension, rsltlen)
    icoords = zeros(Float, cylinder|>dimension, rsltlen)
    for i = 1:cut-1
        push!(pids, (scopes === nothing) ? cylinder.pids[i] : replace(cylinder.pids[i], scope=scopes[i]))
    end
    for (i, p) in enumerate(ps)
        for j = 1:blcklen
            push!(pids, PID(p, (cut == length(cylinder)+1) ? j : (cylinder.pids[cut].site+j-2)%blcklen+1))
        end
    end
    for i = cut:length(cylinder)
        push!(pids, (scopes === nothing) ? cylinder.pids[i] : replace(cylinder.pids[i], scope=scopes[i]))
    end
    for i = 1:blcknum
        trannum = i - (blcknum+1)/2
        for j = 1:blcklen
            for k = 1:dimension(cylinder)
                rcoords[k, (i-1)*blcklen+j] = cylinder.block[k, j] + trannum*cylinder.translation[k]
            end
        end
    end
    cylinder.pids = pids
    cylinder.rcoords = rcoords
    cylinder.icoords = icoords
    if (length(cylinder.neighbors) > 0) && mapreduce(isequal(Inf), |, values(cylinder.neighbors))
        minlens = minimumlengths(rcoords, cylinder.vectors, cylinder|>nneighbor, coordination=coordination)
        cylinder.neighbors = Dict{Int, Float}(i=>len for (i, len) in enumerate(minlens))
    end
    cylinder
end

"""
    (cylinder::Cylinder)(scopes::Any...; coordination::Int=8) -> Lattice

Construct a lattice from a cylinder with the assigned scopes.
"""
function (cylinder::Cylinder)(scopes::Any...; coordination::Int=8)
    @assert fieldtype(cylinder|>keytype, :scope) === scopes|>eltype "cylinder call error: wrong scope type."
    name = cylinder.name * string(length(scopes))
    pids = keytype(cylinder)[PID(scope, i) for scope in scopes for i = 1:size(cylinder.block, 2)]
    rcoords = tile(cylinder.block, [cylinder.translation], Tuple((i,) for i = -(length(scopes)-1)/2:(length(scopes)-1)/2))
    icoords = zero(rcoords)
    vectors = cylinder.vectors
    neighbors = cylinder.neighbors
    if (length(neighbors) > 0) && mapreduce(isequal(Inf), |, values(neighbors))
        minlens = minimumlengths(rcoords, vectors, cylinder|>nneighbor, coordination=coordination)
        neighbors = Dict{Int, Float}(i=>len for (i, len) in enumerate(minlens))
    end
    Lattice{cylinder|>dimension}(name, pids, rcoords, icoords, vectors, neighbors)
end

"""
    Bonds{T, L}(bonds::Tuple{Vararg{Vector{<:AbstractBond}}}) where {T, L<:AbstractLattice}
    Bonds(lattice::AbstractLattice, types::LatticeBonds...)

A set of lattice bonds.

`Bonds` itself is an `AbstractVector` of `AbstractBond`. The need for such a struct is to ensure the type stability during the iteration over a set of different concrete bonds. Although the default `iterate` function does not achieve this goal, users can get it with the generated function trick. Besides, it provides a high level of management of different categories of bonds based on the `LatticeBonds` system.
"""
struct Bonds{T, L<:AbstractLattice, BS<:Tuple{Vararg{Vector{<:AbstractBond}}}, B<:AbstractBond} <: AbstractVector{B}
    bonds::BS
    function Bonds{T, L}(bonds::Tuple{Vararg{Vector{<:AbstractBond}}}) where {T, L<:AbstractLattice}
        @assert isa(T, Tuple{Vararg{LatticeBonds}}) && (length(T) == length(bonds)) "Bonds error: dismatched input types and bonds."
        new{T, L, typeof(bonds), mapreduce(eltype, typejoin, bonds)}(bonds)
    end
end
Bonds(lattice::AbstractLattice) = Bonds(lattice, allbonds)
Bonds(lattice::AbstractLattice, types::LatticeBonds...) = Bonds(lattice, Val(types))
@generated function Bonds(lattice::AbstractLattice, ::Val{VS}) where VS
    types, bonds = [], []
    for V in expand(lattice, Val(VS))
        push!(types, V)
        push!(bonds, :(bonds(lattice, Val($V))))
    end
    types = Tuple(types)
    bonds = Expr(:tuple, bonds...)
    return :(Bonds{$types, typeof(lattice)}($bonds))
end

"""
    ==(bonds1::Bonds, bonds2::Bonds) -> Bool

Overloaded equivalent operator.
"""
Base.:(==)(bonds1::Bonds, bonds2::Bonds) = ==(efficientoperations, bonds1, bonds2)

"""
    isequal(bonds1::Bonds, bonds2::Bonds) -> Bool

Overloaded equivalent function.
"""
Base.isequal(bonds1::Bonds, bonds2::Bonds) = isequal(efficientoperations, bonds1, bonds2)

"""
    size(bs::Bonds) -> Tuple{Int}

Get the size of the set of lattice bonds.
"""
Base.size(bs::Bonds) = (mapreduce(length, +, bs.bonds),)

"""
    summary(io::IO, bs::Bonds)

Print the brief description of a set of lattice bonds to an io.
"""
Base.summary(io::IO, bs::Bonds) = @printf io "%s-element Bonds" length(bs)

"""
    bondtypes(bs::Bonds) -> Tuple{Vararg{LatticeBonds}}
    bondtypes(::Type{<:Bonds{T}}) where T -> Tuple{Vararg{LatticeBonds}}

Get the bondtypes of a set of lattice bonds.
"""
bondtypes(bs::Bonds) = bs |> typeof |> bondtypes
bondtypes(::Type{<:Bonds{T}}) where {T} = T

"""
    latticetype(bs::Bonds)
    latticetype(::Type{<:Bonds{T, L} where T}) where {L<:AbstractLattice}
"""
latticetype(bs::Bonds) = bs |> typeof |> latticetype
latticetype(::Type{<:Bonds{T, L} where T}) where {L<:AbstractLattice} = L

"""
    rank(bs::Bonds) -> Int
    rank(::Type{<:Bonds{T}}) where T -> Int

Get the rank of a set of lattice bonds.
"""
rank(bs::Bonds) = bs |> typeof |> rank
rank(::Type{<:Bonds{T}}) where {T} = length(T)

"""
    iterate(bs::Bonds, state=(1, 0))

Iterate over the lattice bonds in the set.
"""
function Base.iterate(bs::Bonds, state=(1, 0))
    s1, s2 = state[1], state[2]
    while s1 <= length(bs.bonds)
        s1, s2 = (s2+1 > length(bs.bonds[s1])) ? (s1+1, 1) : (s1, s2+1)
        ((s1 > length(bs.bonds)) || (s2 <= length(bs.bonds[s1]))) && break
    end
    (s1 > length(bs.bonds)) && return nothing
    return bs.bonds[s1][s2], (s1, s2)
end

"""
    getindex(bs::Bonds, i::Int) -> eltype(bs)

Get the ith bond in the set.
"""
function Base.getindex(bs::Bonds, i::Int)
    r, k = 1, i
    while (r <= rank(bs)) && (k > length(bs.bonds[r]))
        k = k - length(bs.bonds[r])
        r = r + 1
    end
    (r > rank(bs)) && error("getindex error: attempt to access $(length(bs))-element Bonds at index [$k].")
    return bs.bonds[r][k]
end

"""
    filter(lbs::LatticeBonds, bs::Bonds, choice::Union{Val{:include}, Val{:exclude}}=Val(:include)) -> Bonds
    filter(lbs::Tuple{Vararg{LatticeBonds}}, bs::Bonds, choice::Union{Val{:include}, Val{:exclude}}=Val(:include)) -> Bonds

Get a subset of a set of lattice bonds.

When `choice=Val(:include)`, the lattice bonds indicated by `lbs` will be selected;
When `choice=Val(:exclude)`, the lattice bonds not indicated by `lbs` will be selected.
"""
Base.filter(lbs::LatticeBonds, bs::Bonds, choice::Union{Val{:include}, Val{:exclude}}=Val(:include)) = bondfilter(Val((lbs,)), bs, choice)
Base.filter(lbs::Tuple{Vararg{LatticeBonds}}, bs::Bonds, choice::Union{Val{:include}, Val{:exclude}} = Val(:include)) = bondfilter(Val(lbs), bs, choice)
@generated function bondfilter(::Val{VS}, bs::Bonds, ::Val{:include}) where VS
    types, bonds = [], []
    for V in expand(latticetype(bs), Val(VS))
        index = findfirst(isequal(V), bondtypes(bs))
        if isa(index, Int)
            push!(types, V)
            push!(bonds, :(bs.bonds[$index]))
        end
    end
    types = Tuple(types)
    bonds = Expr(:tuple, bonds...)
    return :(Bonds{$types, latticetype(bs)}($bonds))
end
@generated function bondfilter(::Val{VS}, bs::Bonds, ::Val{:exclude}) where VS
    types, bonds = [], []
    excludes = expand(latticetype(bs), Val(VS))
    for (i, V) in enumerate(bondtypes(bs))
        if V ∉ excludes
            push!(types, V)
            push!(bonds, :(bs.bonds[$i]))
        end
    end
    types = Tuple(types)
    bonds = Expr(:tuple, bonds...)
    return :(Bonds{$types, latticetype(bs)}($bonds))
end

"""
    filter(select::Function, bs::Bonds) -> Bonds

Get a filtered set of bonds by a select function.
"""
@generated function Base.filter(select::Function, bs::Bonds)
    bonds = Expr(:tuple, [:(filter(select, bs.bonds[$i])) for i = 1:rank(bs)]...)
    return :(Bonds{bondtypes(bs), latticetype(bs)}($bonds))
end

"""
    empty!(bs::Bonds) -> Bonds

Empty a set of lattice bonds.
"""
@generated function Base.empty!(bs::Bonds)
    exprs = [:(empty!(bs.bonds[$i])) for i = 1:rank(bs)]
    push!(exprs, :(return bs))
    return Expr(:block, exprs...)
end

"""
    empty(bs::Bonds) -> Bonds

Get an empty copy of a set of lattice bonds.
"""
@generated function Base.empty(bs::Bonds)
    bonds = Expr(:tuple, [:(empty(bs.bonds[$i])) for i = 1:rank(bs)]...)
    return :(Bonds{bondtypes(bs), latticetype(bs)}($bonds))
end

"""
    reset!(bs::Bonds, lattice::AbstractLattice) -> Bonds

Reset a set of lattice bonds by a new lattice.
"""
@generated function reset!(bs::Bonds, lattice::AbstractLattice)
    exprs = []
    push!(exprs, quote
        @assert latticetype(bs) >: typeof(lattice) "reset! error: dismatched bonds and lattice."
        empty!(bs)
    end)
    for i = 1:rank(bs)
        push!(exprs, :(bonds!(bs.bonds[$i], lattice, bondtypes(bs)[$i])))
    end
    push!(exprs, :(return bs))
    return Expr(:block, exprs...)
end

end #module


text:
module TestAdder

using Mimi
using Test

############################################
# adder component without a different name #
############################################

model1 = Model()
set_dimension!(model1, :time, 1:10)
add_comp!(model1, Mimi.adder)

x = collect(1:10)
y = collect(2:2:20)

update_param!(model1, :adder, :input, x)
update_param!(model1, :adder, :add, y)

run(model1)

for i in 1:10
    @test model1[:adder, :output][i] == 3i
end

##############################################
# test adder component with a different name #
##############################################

model2 = Model()
set_dimension!(model2, :time, 1:10)
add_comp!(model2, Mimi.adder, :compA)
update_param!(model2, :compA, :input, x)
update_param!(model2, :compA, :add, y)
run(model2)

for i in 1:10
    @test model2[:compA, :output][i] == 3i
end

end #module

ext: jl


text:
"""Evaluates a set of SCM configurations for a single parameter vector."""

@everywhere using Pkg
@everywhere Pkg.activate("../../..")
@everywhere using ArgParse
@everywhere using CalibrateEDMF
@everywhere using CalibrateEDMF.Pipeline
@everywhere using CalibrateEDMF.TurbulenceConvectionUtils

@everywhere src_dir = dirname(pathof(CalibrateEDMF))
@everywhere using CalibrateEDMF.HelperFuncs
@everywhere include(joinpath(src_dir, "parallel.jl"))
using JLD2

# Read iteration number of ensemble to be recovered
s = ArgParseSettings()
@add_arg_table s begin
    "--version"
    help = "Calibration process number"
    arg_type = Int
    "--job_dir"
    help = "Job output directory"
    arg_type = String
    default = "output"
    "--mode"
    help = "Forward model evaluation mode: `train` or `validation`"
    arg_type = String
    default = "train"
end
parsed_args = parse_args(ARGS, s)
outdir_path = parsed_args["job_dir"]
include(joinpath(outdir_path, "config.jl"))

config = get_config()
version = parsed_args["version"]
mode = parsed_args["mode"]
versioned_model_eval_parallel(version, outdir_path, mode, config)


text:
"""
Authors: Yaroslav Salii, 2020+
         Kara Ignatenko, 2021

Submodule responsible for preprocessing airport and flight data.

airports.jl
2021-08-13 v.0.1: First modular version
2021-08-17 v.0.2: Removed dependency on io.jl, created new func getProcessedAPs
                  to make the control flow explicit rather than rely on default arguments
2021-08-19 v.0.3: Made min_annual_boardings an optional kwarg for getProcessedAPs
2021-09-28 v.0.3.1: fix mkFlightMx2 using a not-quite-string which caused issues in visualization
"""
module Airports
using DataFrames
using FromFile

export FlightMx,
       grpBTS,
       mkFlightMx2,
       dropSmallAPs,
       getProcessedAPs

"""The minimum number of annual boardings for airports to be counted."""
const defaultMinAnnualBoardings = 2500

"""
Data type that describes passenger flows between all airports
in a set.
"""
struct FlightMx
    """An adjacency matrix, M[i1][i2] is the flow from i1 to i2"""
    M::Matrix{Float64}
    """Maps IATA codes to indices in M"""
    ix::Dict{String, Int}
    """Maps indices in M to IATA codes"""
    xi::Dict{Int, String}
    """Lists every IATA code in the set"""
    apCodes::Vector{String}
end

#===BTS FLIGHT DATA PROCESSING===#
"""
Sum the passengers on the flights with same (ORG,DST) pairs,
sort lexicographically in (ORG,DST) order
and output the air travel graph as a *list of edges*.
"""
function grpBTS(idf::DataFrame)
    gd = groupby(idf,[:ORG,:DST])
    out  = combine(gd,:PSG => sum, renamecols = false)
    sort!(out,[:ORG,:DST]) #sort by departure/arrival AP names
end

"""
Given `apCodes [:IATA_Code]`,
outputs a `DataFrame [:IATA_Code,:IN,:OUT,:TOUR]`,
with :IN=Σ_incoming PSG, :OUT=Σ_outgoing,:TOUR=Σ_(:ORG=:DST)
in :IN and :OUT sums, `missing` is non-absorbing and the diagonal is omitted.
(opt) :TTL=:IN+:OUT.
"""
function mkAggFlows2(apCodes::Array{String}, M::Matrix{Float64})
    out=DataFrame(IATA_Code=apCodes
#col-wise total sans the reflexive, `missing` if the arrivals are only reflexive
    ,IN=[ (sum(M[:,j]) == M[j,j]) ? missing : (sum(M[:,j]) - M[j,j]) for j ∈ 1:size(M)[2] ]
#row-wise total sans the reflexive, `missing` if the departures are only reflexive
    ,OUT=[ (sum(M[i,:]) == M[i,i]) ? missing : (sum(M[i,:]) - M[i,i]) for i ∈ 1:size(M)[1] ]
#no. reflexive travelers, or `missing`
    ,TOUR=[M[i,i] for i ∈ 1:size(M)[1]] )#just the reflexive travelers
    return out
end

#===CONSTRUCTING THE FLIGHT MATRIX===#
"""
Transform a *list of arcs* into *adjacency matrix*.
input: DF must have [:ORG,:DST,:PSG] cols
"""
function mkFlightMx2(fs::DataFrame; daily=false::Bool,babble=false::Bool, init_to=0.0::Number)
    #make a sorted list of ALL the APs in `fs`
    allAPs = sort( (fs.ORG |> unique) ∪ (fs.DST |> unique))
    #delegate to the method with EXPLICIT ground set (`iretAPs`); with retaintours = true
    #note the map(String, allAPs) - this is because DataFrames uses an internal
    # string-like type that causes problems down the line
    return mkFlightMx2(fs, map(String, allAPs); daily=daily,babble=babble,init_to=init_to,retaintours=true)
end

"""
Transform a *list of arcs* into *adjacency matrix*,
taking an *explicit* list of vertices `iretAPs`.
input: DF must have [:ORG,:DST,:PSG] cols
"""
function mkFlightMx2(fs::DataFrame, iretAPs::Array{String};
    daily=true::Bool # divide by 365, the input was annualized
    ,retaintours = false::Bool #flights with ORG==DST not retained by default
    ,babble=true::Bool, init_to=0.0::Number)

    retAPs = sort(iretAPs) #ensure AP codes are sorted, for indexing porposes
    pret = a -> a.ORG ∈ retAPs && a.DST ∈ retAPs
    if !retaintours
        pret = a -> a.ORG ≠ a.DST &&   a.ORG ∈ retAPs && a.DST ∈ retAPs
    end
    #pretain(a)  ? (retaintours) : ()
    #retain only flights (tuples :ORG,:DST,PSG) for APs ∈ retAPs;
    retFlows = filter(pret, eachrow(fs)) |> DataFrame

    if babble #report if there were *any* isolated APs
        isolatedAPs = filter( a -> a ∉ retFlows.ORG && a ∉ retFlows.DST, retAPs)
        println("Found ", isolatedAPs |> length," isolated APs: ",isolatedAPs)
    end

    dim = length(retAPs)
    ix_ = zip(retAPs, 1:dim) |> Dict #get index by name
    xi_ = zip(1:dim,retAPs) |> Dict #get name by index
    M_ = fill(init_to,(dim,dim)) #default is 0.0; screw `missing`
    #fill the matrix with the *known* values
    for row ∈ eachrow(retFlows)
        M_[ix_[row.ORG],ix_[row.DST]] = row.PSG
    end

    outM = daily ? map(x -> x/365,M_) : M_ #annual-to-daily, if requested
    return FlightMx(outM, ix_, xi_, retAPs)
end

#===LINKING BTS DATA WITH OPENFLIGHTS===#
#=====BTS=<-->=OPENFLIGHTS=======#
"""
Pick just the given APs (by IATA_Code) from all in dfAPinfo and
output columns as [:IATA_Code,:LAT,:LNG,:IN,:OUT,:TOUR,:Name,:City,:Country]
"""
function pickAPs(myAPs::DataFrame, dfAPinfo::DataFrame)
    out= innerjoin(myAPs, dfAPinfo,on=:IATA_Code)
    select!(out,[:IATA_Code,:LAT,:LNG,:IN,:OUT,:TOUR,:Name,:City,:Country])
end

"""
Pick just the “clean”  APs (no `missing`) from all in dfAPinfo.
CAVEAT: doesn't retain any APs that have :IN or:OUT missing or 0;
recall that 0 was identified with `missing` in mkAggFlows()=
Also remove missings from the :IN and :OUT.
"""
function pickCleanAPs(myAPs::DataFrame, dfAPinfo::DataFrame)
    pickAPs(myAPs,dfAPinfo) |> scrubAPs
end

"""Throw out the APs with missing IN or OUT enplanements"""
scrubAPs(APs::DataFrame) = filter(row -> !ismissing(row.OUT) && !ismissing(row.IN)
,eachrow(APs)) |> DataFrame

#TBD: enplanement is just OUT, so maybe not sum it with in.
"""Only retain the APs with at least `p` annual enplanements (sum :IN and :OUT)"""
function dropSmallAPs(idf::DataFrame, p=defaultMinAnnualBoardings::Number)
    if (["IN","OUT"] ⊈ names(idf)) error("wrong DF") end
    filter(row -> row.IN+row.OUT ≥ p, eachrow(idf)) |> DataFrame
end

"""
Given flight data `pBTS` and raw airport data `rawAPs`, calculate AP in- and outflows
and get rid of the APs with missing data or less than `min_annual_boardings` annual passengers.
Input:
`pBTS`: [:ORG, :DST]
`rawAPs`: [:IATA_Code,:LAT,:LNG,:Name,:City,:Country]
Output: [:IATA_Code,:LAT,:LNG,:IN,:OUT,:TOUR,:Name,:City,:Country]
"""
function getProcessedAPs(pBTS::DataFrame, rawAPs::DataFrame; min_annual_boardings=defaultMinAnnualBoardings)
    fmx = mkFlightMx2(pBTS)
    aggFlows = mkAggFlows2(fmx.apCodes, fmx.M)
    cleanAPs = pickCleanAPs(aggFlows, rawAPs)
    dropSmallAPs(cleanAPs, min_annual_boardings)
end

end


text:
include("../RLAlgorithms.jl")
using .RLAlgorithms
using Plots
using Flux
using Flux: Optimise.update!, Zygote.pullback
using Statistics: mean
using StatsBase: sample, Weights

function act(policy, s)
    logits = policy(s)
    a = sample(Weights(exp.(logits)))
    loga = gather(logits, a) - log.(sum(exp.(logits), dims=1))
    return a, loga
end

approxkl(logp, oldlogp) = mean(oldlogp .- logp)

function ploss(logp, oldlogp, adv, ϵ)
    rat = exp.(logp .- oldlogp)
    return -mean(min.(rat .* adv, min.(max.(rat, one(ϵ)-ϵ), one(ϵ)+ϵ) .* adv))
end

function run!(policy, popt, value, vopt, env;
        epochs=100, steps=300, maxt=100, ptrainiters=50, vtrainiters=80,
        ϵ=2f-1, γ=99f-2, λ=97f-2, targetkl=15f-3)
    rewards = zeros(Float32, epochs)
    memory = PGReplayMemory{Float32, Int64, Float32, Float32}(
        length(env.observationspace), 1, steps, steps, γ=γ, λ=λ
    )
    policyps = params(policy)
    valueps = params(value)

    for i in 1:epochs
        episodes!(memory, env, steps, maxt) do s
            a, loga = act(policy, s)
            v = value(s)
            return a, loga[1], v[1]
        end

        let (s, a, oldlogp, r, adv) = sample(memory)
            for _ in 1:ptrainiters
                loga, back = pullback(policyps) do
                    logits = policy(s)
                    return gather(logits, a) .- vec(log.(sum(exp.(logits), dims=1)))
                end
                approxkl(loga, oldlogp) > targetkl && break
                gs = gradient(logp -> ploss(logp, oldlogp, adv, ϵ), loga)[1]
                update!(popt, policyps, back(gs))
            end

            Flux.train!(valueps, Iterators.repeated((s, r), vtrainiters), vopt) do s, r
                return Flux.mse(vec(value(s)), r)
            end
        end
        rewards[i] = test(s -> Flux.onecold(policy(s)), env)
    end
    return rewards
end

env = CartPole{Float32}()

policy = Chain(
    Dense(length(env.observationspace), 32, tanh),
    Dense(32, length(env.actionspace))
)
value = Chain(
    Dense(length(env.observationspace), 32, tanh),
    Dense(32, 1)
)

r = run!(policy, ADAM(3e-4), value, ADAM(1e-3), env)
plt = plot(r, labels="Reward")
display(plt)
#savefig(plt, "ppodisc")


text:
# Model
include_model("hopper")

mb = 3.0 # body mass
ml = 0.3  # leg mass
Jb = 0.75 # body inertia
Jl = 0.075 # leg inertia

model = Hopper{Discrete, FixedTime}(n, m, d,
			   mb, ml, Jb, Jl,
			   0.25, g,
			   qL, qU,
			   uL, uU,
			   nq,
		       nu,
		       nc,
		       nf,
		       nb,
		   	   ns,
		       idx_u,
		       idx_λ,
		       idx_b,
		       idx_ψ,
		       idx_η,
		       idx_s)

model_ft = free_time_model(model)

# stair traj
@load joinpath(@__DIR__, "hopper_stair.jld2") qm um γm bm ψm ηm μm hm

# Stair
function ϕ_func(model::Hopper, q)
	k = kinematics(model, q)
	if k[1] < 0.125
    	return @SVector [k[2] - 3 * 0.25]
	else
		return @SVector [k[2]]
	end
end

# Horizon
T = 80

# Time step
# tf = 0.75
h = hm #tf / (T - 1)

# Bounds
_uu = Inf * ones(model_ft.m)
_uu[model_ft.idx_u] .= Inf#10.0#model_ft.uU
_uu[end] = h#3.0 * h
_ul = zeros(model_ft.m)
_ul[model_ft.idx_u] .= -Inf#10.0 #model_ft.uL
_ul[end] = h#0.2 * h
ul, uu = control_bounds(model_ft, T, _ul, _uu)

# Initial and final states
z_h = 3 * 0.25

q1 = [0.0, 0.5 + z_h, 0.0, 0.5]
q11 = [0.0, 0.5 + z_h + 0.125, 0.0, 0.25]
qm1 = [0.125, 0.5 + z_h + 0.25, -0.5 * π, 0.25]
qm2 = [0.25, 0.5 + z_h + 0.125, -1.5 * π, 0.25]
qm3 = [0.375, 0.5 + z_h + 0.0625, -2.0 * π, 0.5]
qT = [0.5, 0.5, -2.0 * π, 0.5]

ql1 = linear_interpolation(q1, q11, 14)
ql2 = linear_interpolation(q11, qm1, 15)
ql3 = linear_interpolation(qm1, qm2, 15)
ql4 = linear_interpolation(qm2, qm3, 15)
ql5 = linear_interpolation(qm3, qT, 14)
ql6 = linear_interpolation(qT, qT, 14)

q_ref = [ql1...,
		 ql2[2:end]...,
		 ql3[2:end]...,
		 ql4[2:end]...,
		 ql5[2:end]...,
		 ql6[2:end]...]

θr = range(0.0, stop = -2.0 * π, length = (58 - 14))
for (i, t) = enumerate(15:58)
	q_ref[t][3] = θr[i]
end
# model_ft.qU[2] = 2.5
xl, xu = state_bounds(model_ft, T,
		[model_ft.qL; model_ft.qL],
		[model_ft.qU; model_ft.qU],
        x1 = [q1; q1],
		xT = [qT; qT])

# Objective
include_objective(["velocity", "nonlinear_stage", "control_velocity"])
qp = 0.0 * [0.01; 0.01; 1.0; 1.0]
obj_tracking = quadratic_time_tracking_objective(
    [Diagonal([qp; qp]) for t = 1:T],
    [Diagonal([1.0e-1, 1.0e-2,
		1.0e-5 * ones(model_ft.nc)..., 1.0e-5 * ones(model_ft.nb)...,
		zeros(model_ft.m - model_ft.nu - model_ft.nc - model_ft.nb - 1)..., 0.0])
		for t = 1:T-1],
    [[qT; qT] for t = 1:T],
    [zeros(model_ft.m) for t = 1:T],
    1.0)

obj_contact_penalty = PenaltyObjective(1.0e5, model_ft.m - 1)

obj_velocity = velocity_objective(
    [(t > 20 && t < 60) ? Diagonal(1.0e-2 * [1.0; 1.0; 100.0; 100.0]) : Diagonal(1.0e-2 * [1.0; 1.0; 1.0; 1.0]) for t = 1:T-1],
    model_ft.nq,
    h = h,
    idx_angle = collect([3]))

obj_ctrl_vel = control_velocity_objective(Diagonal([1.0e-1 * ones(model_ft.nu);
	1.0e-3 * ones(model_ft.nc + model_ft.nb);
	zeros(model_ft.m - model_ft.nu - model_ft.nc - model_ft.nb)]))

function l_stage(x, u, t)
	J = 0.0

	_q1 = view(x, 1:4)
	p1 = kinematics(model, _q1)

	_q2 = view(x, 4 .+ (1:4))
	p2 = kinematics(model, _q2)

	v = (p2 - p1) ./ h

	if t < 40
		J += 1000.0 * v[1]^2.0
	end

	if t > 60
		J += 1000.0 * v[1]^2.0
	end

	if true#t > 5 #|| (t > 20 && t < T)
		J += (_q1 - q_ref[t])' * Diagonal([100.0; 100.0; 1000.0; 1000.0]) * (_q1 - q_ref[t])
	end

	return J
end

l_stage(x) = l_stage(x, nothing, T)

obj_stage = nonlinear_stage_objective(l_stage, l_stage)
obj = MultiObjective([obj_tracking,
	obj_contact_penalty,
	obj_velocity,
	obj_stage,
	obj_ctrl_vel])

# Constraints
include_constraints(["free_time", "contact", "stage"])
con_free_time = free_time_constraints(T)
con_contact = contact_constraints(model_ft, T)

p1_ref = kinematics(model, q1)
pT_ref = kinematics(model, qT)

function pinned1!(c, x, u, t)
    q = view(x, 1:4)
    c[1:2] = p1_ref - kinematics(model, q)
	nothing
end

function pinnedT!(c, x, u, t)
    q = view(x, 4 .+ (1:4))
	c[1:2] = pT_ref - kinematics(model, q)
	nothing
end

function no_foot_slip!(c, x, u, t)
	q = view(x, 1:4)

	c[1] = kinematics(model, q)[1]
end

n_stage = 2
t_idx1 = vcat([t for t = 1:10]...)
t_idxT = vcat([(T - 10 + 1):T]...)

con_pinned1 = stage_constraints(pinned1!, 2, (1:0), t_idx1)
con_pinnedT = stage_constraints(pinnedT!, 2, (1:0), t_idxT)
con_no_slip = stage_constraints(no_foot_slip!, 1, (1:1), collect(1:40))

con = multiple_constraints([con_free_time, con_contact, con_pinned1, con_pinnedT, con_no_slip])#, con_loop])

# Problem
prob = trajectory_optimization_problem(model_ft,
               obj,
               T,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
x0 = configuration_to_state(q_ref) # linear interpolation on state
u0 = [[1.0e-2 * rand(model_ft.m-1); h] for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
@time z̄, info = solve(prob, copy(z0),
	nlp = :ipopt,
	tol = 1.0e-3, c_tol = 1.0e-3, mapl = 5,
	time_limit = 60)
@show check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
tf, t, h̄ = get_time(ū)

q = state_to_configuration(x̄)
u = [u[model.idx_u] for u in ū]
γ = [u[model.idx_λ] for u in ū]
b = [u[model.idx_b] for u in ū]
ψ = [u[model.idx_ψ] for u in ū]
η = [u[model.idx_η] for u in ū]
h̄ = mean(h̄)
# @save joinpath(pwd(), "examples/trajectories/hopper_vertical_gait.jld2") z̄ x̄ ū h̄ q u γ b

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
open(vis)

visualize!(vis, model_ft,
	q,
	Δt = h̄[1],
	scenario = :vertical)

setobject!(vis["box"], GeometryBasics.HyperRectangle(Vec(0.0, 0.0, 0.0),
	Vec(0.25, 0.5, 3 * 0.25)), MeshPhongMaterial(color = RGBA(0.5, 0.5, 0.5, 1.0)))
settransform!(vis["box"], Translation(-0.125, -0.25, 0))

using Plots
plot(hcat(q_ref...)[1:4, :]', color = :black , width = 2.0)
plot!(hcat(q...)[1:4, :]', color = :red, width = 1.0)

plot(hcat(u...)', linetype = :steppost)

# Save
hm = h̄
μm = model.μ
qm, um, γm, bm, ψm, ηm = q, u, γ, b, ψ, η

@save joinpath(@__DIR__, "hopper_tall_flip.jld2") qm um γm bm ψm ηm μm hm

# composite stairs + flip
@load joinpath(@__DIR__, "hopper_stair.jld2") qm um γm bm ψm ηm μm hm

function step_repeat(q, u, γ, b, ψ, η, T; steps = 2)
	qm = [deepcopy(q)...]
	um = [deepcopy(u)...]
	γm = [deepcopy(γ)...]
	bm = [deepcopy(b)...]
	ψm = [deepcopy(ψ)...]
	ηm = [deepcopy(η)...]

	stride = zero(qm[1])
	for i = 1:(steps-1)
		@show stride[1] += q[T+1][1] - q[2][1]
		@show stride[2] += 0.25
		for t = 1:T-1
			push!(qm, q[t+2] + stride)
			push!(um, u[t])
			push!(γm, γ[t])
			push!(bm, b[t])
			push!(ψm, ψ[t])
			push!(ηm, η[t])
		end
	end

	return qm, um, γm, bm, ψm, ηm
end

qm, um, γm, bm, ψm, ηm = step_repeat(qm, um, γm, bm, ψm, ηm, T, steps = 3)

@save joinpath(@__DIR__, "hopper_stairs_3.jld2") qm um γm bm ψm ηm μm hm
# @load joinpath(@__DIR__, "hopper_stairs_3.jld2") qm um γm bm ψm ηm μm hm

setobject!(vis["box1"], GeometryBasics.HyperRectangle(Vec(0.0, 0.0, 0.0),
	Vec(0.25, 0.5, 0.25)), MeshPhongMaterial(color = RGBA(0.5, 0.5, 0.5, 1.0)))
settransform!(vis["box1"], Translation(0.125, -0.25, 0))

setobject!(vis["box2"], GeometryBasics.HyperRectangle(Vec(0.0, 0.0, 0.0),
	Vec(0.25, 0.5, 2 * 0.25)), MeshPhongMaterial(color = RGBA(0.5, 0.5, 0.5, 1.0)))
settransform!(vis["box2"], Translation(0.125 + 0.25, -0.25, 0))

setobject!(vis["box3"], GeometryBasics.HyperRectangle(Vec(0.0, 0.0, 0.0),
	Vec(0.25, 0.5, 3 * 0.25)), MeshPhongMaterial(color = RGBA(0.5, 0.5, 0.5, 1.0)))
settransform!(vis["box3"], Translation(0.125 + 2 * 0.25, -0.25, 0))

tall_flip = load(joinpath(@__DIR__, "hopper_tall_flip.jld2"))

qm_f, um_f, γm_f, bm_f, ψm_f, ηm_f, μm_f, hm_f = tall_flip["qm"], tall_flip["um"], tall_flip["γm"], tall_flip["bm"], tall_flip["ψm"], tall_flip["ηm"], tall_flip["μm"], tall_flip["hm"]

str = zero(qm[1])
str[1] = qm[end][1]

for i = 1:10
	t = 1
	push!(qm, qm_f[t+2] + str)
	push!(um, um_f[t])
	push!(γm, γm_f[t])
	push!(bm, bm_f[t])
	push!(ψm, ψm_f[t])
	push!(ηm, ηm_f[t])
end

for t = 1:length(um_f)
	push!(qm, qm_f[t+2] + str)
	push!(um, um_f[t])
	push!(γm, γm_f[t])
	push!(bm, bm_f[t])
	push!(ψm, ψm_f[t])
	push!(ηm, ηm_f[t])
end

@save joinpath(@__DIR__, "hopper_stairs_3_flip.jld2") qm um γm bm ψm ηm μm hm

visualize!(vis, model_ft,
	qm,
	Δt = h̄[1],
	scenario = :stairs)

setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 20)

settransform!(vis["/Cameras/default"],
	compose(Translation(0.0, -90.0, -1.0),LinearMap(RotZ(-0.5 * π))))

max_issues_repo_name: thowell/DirectMotionPlanning


text:
#=
	This is the main function, it takes as input a starting array, and the number of steps.
	Randomizing the initial grid is also built-in through the function arguments.
=#
using Base
using Distributed
include("mapping_functions.jl")
@everywhere Base.MainInclude.include("mapping_functions.jl")
using Plots
#using ProgressMeter
#TODO: Transfer to Images package instead of built-in imfilter
import Plots.@animate
@everywhere using ImageFiltering

#clibrary(:cmocean)

#module GasLattice

function GasLattice_visual(steps=50,startarr=0,dim=(1000,1000))
	if startarr == 0
		randomize = true

	else
		randomize = false
		dim = size(startarr)
	end
    randomize ? Grid=rand([1 2 3 5 6],dim) : Grid=startarr           #2:up 3:left 5:down 7:right 1:none
    filtertop  =   centered([0 0 0;0 0 0;0 1 0])		#Filter of top cells
    filterbot  =   centered([0 1 0;0 0 0;0 0 0])		#Filter of bottom cells
    filterleft =   centered([0 0 0;0 0 1;0 0 0])		#Filter of left cells
    filterright=   centered([0 0 0;1 0 0;0 0 0])		#Filter of right cells
    rang=1:dim[1]				#Range for limiting the convolution product result

    heatmap(map(countt,Grid),c=:grayC)			#Initializer for the heatmap

    #p1=Progress(steps)				#Progress bar initializer
    anim=  @animate for k=1:steps
    #for k=1:steps
        Grid=padarray(Grid,Fill(0,(1,1),(1,1)))
    	  ##==The following are mapping each cell to the incoming particle.==#
        # @sync begin
        #     top_=@spawnat 1 imfilter(Grid,filtertop)[rang,rang]
        #     bot_=@spawnat 2 imfilter(Grid,filterbot)[rang,rang]
        #     rig_=@spawnat 3 imfilter(Grid,filterright)[rang,rang]
        #     lef_=@spawnat 4 imfilter(Grid,filterleft)[rang,rang]
        # end

        Grid_top=@spawnat 1 map(hasbot,imfilter(Grid,filtertop)[rang,rang])
        Grid_left=@spawnat 2 map(hasright,imfilter(Grid,filterleft)[rang,rang])
        Grid_bot=@spawnat 3 map(hastop,imfilter(Grid,filterbot)[rang,rang])
        Grid_right=@spawnat 4 map(hasleft,imfilter(Grid,filterright)[rang,rang])
        ###
        Grid_post_filter=fetch(Grid_top) .* fetch(Grid_left) .* fetch(Grid_right) .* fetch(Grid_bot)
        #Grid_post_filter=fetch(Grid_top) .* fetch(Grid_left) .* fetch(Grid_right) .* fetch(Grid_bot)	#The product of the different maps
        Grid_post_filter=Grid_post_filter[rang,rang]
	    Grid_pre_collision=mapzeroes(Grid[rang,rang],Grid_post_filter)	#This line maps the zeros of the initial grid
								#to their corresponding cells in the product
								#It is necessary to support inert objects in the box

	    Grid_pre_wall=map(collision,Grid_pre_collision)		#Mapping particle-particle collisions.
        ###
        Grid=map(bordercollisions,Grid_pre_wall)		#Mapping particle-border collisions.

        Grid_count=map(countt,Grid)				#Counting the number of particles in a cell
				#For the heatmap,.
        heatmap(Grid_count,c=:grayC)					#Next frame.
	      #next!(p1)
	#end
	      				#Next step in the progress bar.
    end
    gif(anim,"randomgaslattice.gif",fps=3)			#Create the gif.

end
#export GasLattice_visual
#end


text:
function evol(fitness, lb, ub, numParticles, maxiter, verbose)

    sr = [(lb[i], ub[i]) for i=1:length(lb)]

    fopt = 10000
    xopt = []
    for i=1:15
        println(i)
        result = BlackBoxOptim.bboptimize(fitness; SearchRange = sr, NumDimensions = length(lb),
                    Method = :adaptive_de_rand_1_bin_radiuslimited, MaxSteps = 18000,
                    TraceMode = :silent)
        f_i = BlackBoxOptim.best_fitness(result)
        if f_i < fopt
            fopt = f_i
            println(fopt)
            xopt = BlackBoxOptim.best_candidate(result)
        end
    end

    sr = [(xopt[i]-0.05*abs(xopt[i]), xopt[i]+0.05*abs(xopt[i])) for i=1:length(lb)]
    sr[end] = (lb[end], lb[end])

    result = BlackBoxOptim.bboptimize(fitness; SearchRange = sr, NumDimensions = length(lb),
            Method = :adaptive_de_rand_1_bin_radiuslimited, MaxSteps = 30000,
            TraceMode = :silent)
    fopt = BlackBoxOptim.best_fitness(result)
    xopt = BlackBoxOptim.best_candidate(result)

    return xopt, fopt

end


text:
# Autogenerated wrapper script for SFML_jll for x86_64-apple-darwin
export libsfml_audio, libsfml_graphics, libsfml_network, libsfml_system, libsfml_window

using Libglvnd_jll
using Ogg_jll
using FLAC_jll
using FreeType2_jll
using libvorbis_jll
using Xorg_libXrandr_jll
using Xorg_libX11_jll
using OpenAL_jll
using eudev_jll
JLLWrappers.@generate_wrapper_header("SFML")
JLLWrappers.@declare_library_product(libsfml_audio, "@rpath/libsfml-audio.2.5.dylib")
JLLWrappers.@declare_library_product(libsfml_graphics, "@rpath/libsfml-graphics.2.5.dylib")
JLLWrappers.@declare_library_product(libsfml_network, "@rpath/libsfml-network.2.5.dylib")
JLLWrappers.@declare_library_product(libsfml_system, "@rpath/libsfml-system.2.5.dylib")
JLLWrappers.@declare_library_product(libsfml_window, "@rpath/libsfml-window.2.5.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Libglvnd_jll, Ogg_jll, FLAC_jll, FreeType2_jll, libvorbis_jll, Xorg_libXrandr_jll, Xorg_libX11_jll, OpenAL_jll, eudev_jll)
    JLLWrappers.@init_library_product(
        libsfml_audio,
        "lib/libsfml-audio.2.5.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsfml_graphics,
        "lib/libsfml-graphics.2.5.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsfml_network,
        "lib/libsfml-network.2.5.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsfml_system,
        "lib/libsfml-system.2.5.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsfml_window,
        "lib/libsfml-window.2.5.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Test generator models using a multiperiod optimization

# - Model: `_FP.BFARadPowerModel` is used to be able to test generators in both active and
#   reactive power, while keeping the model linear. It requires a distribution network.
# - Problem: `flex_tnep`.

## Settings

file = normpath(@__DIR__,"..","test","data","case2","case2_d_gen.m") # Input case. Here 2-bus distribution network having 1 dispatchable generator, 1 non-dispatchable generator and 1 fixed load, all on bus 1 (bus 2 is empty).
number_of_hours = 5 # Number of time periods

## Plot function

# Uncomment this part and the commented lines further down to display some nice plots when manually editing a testset
#=
using StatsPlots
function plot_bus(mn_data, result)
    res_gen(i,key) = [result["solution"]["nw"]["$n"]["gen"]["$i"][key] for n in 1:number_of_hours]
    data_gen(i,key) = [mn_data["nw"]["$n"]["gen"]["$i"][key] for n in 1:number_of_hours]
    data_load(i,key) = [mn_data["nw"]["$n"]["load"]["$i"][key] for n in 1:number_of_hours]
    gen_matrix = hcat(res_gen(1,"pgcurt"), res_gen(2,"pgcurt"), res_gen(1,"pg"), res_gen(2,"pg")) # Rows: hours; columns: power categories
    plt = groupedbar(gen_matrix;
        title = "Bus 1",
        yguide = "Power [p.u.]",
        xguide = "Time [h]",
        framestyle = :zerolines,
        bar_position = :stack,
        bar_width = 1,
        linecolor = HSLA(0,0,1,0),
        legend_position = :topleft,
        label = ["gen1 pgcurt" "gen2 pgcurt" "gen1 pg" "gen2 pg"],
        seriescolor = [HSLA(0,0.5,0.5,0.5) HSLA(0,0.75,0.25,0.5) HSLA(210,0.75,0.5,0.5) HSLA(210,1,0.25,0.5)],
    )
    plot!(plt, data_load(1,"pd"); label="demand", seriestype=:stepmid, linecolor=:black, linewidth=2, linestyle=:dot)
    display(plt)
end
function plot_gen(mn_data, result, i)
    res_gen(i,key) = [result["solution"]["nw"]["$n"]["gen"]["$i"][key] for n in 1:number_of_hours]
    data_gen(i,key) = [mn_data["nw"]["$n"]["gen"]["$i"][key] for n in 1:number_of_hours]
    gen_matrix = hcat(res_gen.(i,["pgcurt" "pg"])...) # Rows: hours; columns: power categories
    plt = groupedbar(gen_matrix;
        title = "Generator $i",
        yguide = "Power [p.u.]",
        xguide = "Time [h]",
        framestyle = :zerolines,
        bar_position = :stack,
        bar_width = 1,
        linecolor = HSLA(0,0,1,0),
        legend_position = :topleft,
        label = ["pgcurt" "pg"],
        seriescolor = [HSLA(0,0.75,0.25,0.5) HSLA(210,0.75,0.5,0.5)],
    )
    plot!(plt, data_gen(i,"pmax"); label="pmax", seriestype=:stepmid, linecolor=:black, linewidth=2, linestyle=:dot)
    plot!(plt, data_gen(i,"pmin"); label="pmin", seriestype=:stepmid, linecolor=:black, linewidth=1, linestyle=:dash)
    display(plt)
end
=#

## Test results

@testset "Generator model" begin

    # The power required by a fixed load linearly increases from 0 to 20 MW. Generator 2 is
    # non-dispatchable and its reference power decreases from 20 MW to 0 MW, so it is
    # curtailed in the first half of the time horizon and used at full power in the second
    # half. Generator 1, which is dispatchable and can range from 0 to 15 MW, covers the
    # rest of the demand in subsequent periods, until reaches its maximum power; after that,
    # the load is curtailed.
    data = _FP.parse_file(file)
    _FP.add_dimension!(data, :hour, number_of_hours)
    _FP.add_dimension!(data, :year, 1; metadata = Dict{String,Any}("scale_factor"=>1))
    _FP.scale_data!(data; cost_scale_factor=1e-6)
    loadprofile = collect(reshape(range(0,2;length=number_of_hours),:,1)) # Create a load profile: ramp from 0 to 2 times the rated value of load
    genprofile = hcat(1.5.*ones(number_of_hours), reverse(loadprofile; dims=1)) # Generator 1: 2 times the rated value; generator 2: ramp from 2 times the rated value to 0
    time_series = _FP.make_time_series(data; loadprofile, genprofile) # Compute time series by multiplying the rated value by the profile
    mn_data = _FP.make_multinetwork(data, time_series)
    result = _FP.flex_tnep(mn_data, _FP.BFARadPowerModel, cbc)
    #plot_bus(mn_data, result)

    @testset "Dispatchable generator" begin
        #plot_gen(mn_data, result, 1)
        @test result["solution"]["nw"]["1"]["gen"]["1"]["pg"]     ≈  0.0 atol=1e-3
        @test result["solution"]["nw"]["2"]["gen"]["1"]["pg"]     ≈  0.0 atol=1e-3
        @test result["solution"]["nw"]["3"]["gen"]["1"]["pg"]     ≈  0.0 atol=1e-3 # Here demand is covered by generator 2, which is non-dispatchable
        @test result["solution"]["nw"]["4"]["gen"]["1"]["pg"]     ≈ 10.0 rtol=1e-3
        @test result["solution"]["nw"]["5"]["gen"]["1"]["pg"]     ≈ 15.0 rtol=1e-3 # Must not exceed `pmax` even if the load requires more power
        @test result["solution"]["nw"]["1"]["gen"]["1"]["pgcurt"] ≈  0.0 atol=1e-3 # Dispatchable generators are not curtailable: `pgcurt` is always zero
        @test result["solution"]["nw"]["2"]["gen"]["1"]["pgcurt"] ≈  0.0 atol=1e-3
        @test result["solution"]["nw"]["3"]["gen"]["1"]["pgcurt"] ≈  0.0 atol=1e-3
        @test result["solution"]["nw"]["4"]["gen"]["1"]["pgcurt"] ≈  0.0 atol=1e-3
        @test result["solution"]["nw"]["5"]["gen"]["1"]["pgcurt"] ≈  0.0 atol=1e-3
    end

    @testset "Non-dispatchable generator" begin
        #plot_gen(mn_data, result, 2)
        @test result["solution"]["nw"]["1"]["gen"]["2"]["pg"]     ≈  0.0 atol=1e-3 # Curtailment is the only way to decrease generated power; here is completely exploited
        @test result["solution"]["nw"]["2"]["gen"]["2"]["pg"]     ≈  5.0 rtol=1e-3
        @test result["solution"]["nw"]["3"]["gen"]["2"]["pg"]     ≈ 10.0 rtol=1e-3
        @test result["solution"]["nw"]["4"]["gen"]["2"]["pg"]     ≈  5.0 rtol=1e-3 # Must not exceed `pmax`; the rest of the demand is covered by generator 1
        @test result["solution"]["nw"]["5"]["gen"]["2"]["pg"]     ≈  0.0 atol=1e-3 # Must not exceed `pmax` even if the load requires more power
        @test result["solution"]["nw"]["1"]["gen"]["2"]["pgcurt"] ≈ 20.0 rtol=1e-3 # Curtailment is the only way to decrease generated power; here is completely exploited
        @test result["solution"]["nw"]["2"]["gen"]["2"]["pgcurt"] ≈ 10.0 rtol=1e-3
        @test result["solution"]["nw"]["3"]["gen"]["2"]["pgcurt"] ≈  0.0 atol=1e-3
        @test result["solution"]["nw"]["4"]["gen"]["2"]["pgcurt"] ≈  0.0 atol=1e-3
        @test result["solution"]["nw"]["5"]["gen"]["2"]["pgcurt"] ≈  0.0 atol=1e-3
    end

end;


text:
const OS = Dict(:Darwin => 'm', :Linux => 'l')
const QARCH = string(OS[Sys.KERNEL], Sys.WORD_SIZE)
try
  global const QHOME = ENV["QHOME"]
catch
  global const QHOME = joinpath(homedir(), "q")
end
const EMBED_DIR = joinpath(@__DIR__, "..", "embed")
if isdir(joinpath(QHOME, QARCH))
  run(`make -C $(EMBED_DIR) QHOME=$(QHOME) QARCH=$(QARCH) install`)
  run(`make -C $(EMBED_DIR) QHOME=$(QHOME) QARCH=$(QARCH) clean`)
  info("Installed server components for $(QARCH) in $(QHOME).")
else
  warn("Could not find a suitable kdb+ installation.",
       "\nServer components are not installed.")
end


text:
function heqbos!(F, x)
    n = length(x)
    c = 1.0
    mu = 0.5:1:n-0.5
    mu = mu / n
    h = 1.0 / n
    cval = sqrt(1.0 - c)
    A = zeros(n, n)
    for j = 1:n
        for i = 1:n
            A[i, j] = mu[j] / (mu[i] + mu[j])
        end
    end
    A = (c / 2) * h * A
    F .= (A * x)
    for ig = 1:n
        F[ig] = 1.0 / (cval + F[ig])
    end
    F .= x - F
end

function chandprintbos(x, mu, c)
    muc = 0:0.05:1
    n = length(mu)
    nx = length(x)
    LC = zeros(21, n)
    for j = 1:n
        for i = 1:21
            LC[i, j] = muc[i] / (muc[i] + mu[j])
        end
    end
    p = c * 0.5 / n
    LC = p * LC
    hout = LC * x
    hout = ones(21, 1) ./ (ones(21, 1) - hout)
    return [muc hout]
end

max_stars_repo_name: ctkelley/MPResults


text:
#!/usr/bin/env julia

import Celeste.ParallelRun: one_node_single_infer, infer_init, BoundingBox
import Celeste.SDSSIO: RunCamcolField, load_field_images
import Celeste.Infer: find_neighbors
import Celeste.DeterministicVIImagePSF: infer_source_fft

const datadir = joinpath(Pkg.dir("Celeste"), "test", "data")
wd = pwd()
cd(datadir)
run(`make`)
cd(wd)

"""
test infer with a single (run, camcol, field).
This is basically just to make sure it runs at all.
"""
function benchmark_infer()
    # very small patch of sky that turns out to have 4 sources.
    # We checked that this patch is in the given field.
    box = BoundingBox(164.39, 164.41, 39.11, 39.13)
    rcfs = [RunCamcolField(3900, 6, 269),]

    catalog, target_sources = infer_init(rcfs, datadir; box=box)
    images = load_field_images(rcfs, datadir)
    neighbor_map = find_neighbors(target_sources, catalog, images)
    ctni = (catalog, target_sources, neighbor_map, images)

    # Warm up---this compiles the code
    one_node_single_infer(ctni...; infer_source_callback=infer_source_fft)

    # clear allocations in case julia is running with --track-allocations=user
    Profile.clear_malloc_data()

    if isempty(ARGS)
        @time one_node_single_infer(ctni...; infer_source_callback=infer_source_fft)
    elseif ARGS[1] == "--profile"
        Profile.clear_malloc_data()
	Profile.init(delay=0.01)
        # about half the run time is psf fitting, the other half is elbo evaluation
        @profile one_node_single_infer(ctni...; infer_source_callback=infer_source_fft)
        Profile.print(format=:flat, sortedby=:count)
    end
end

benchmark_infer()


text:
"""
    Phases
"""
module Phases

"""
    abstract type Phase

Abstract supertype for all phases. See `subtypes(FluxTraining.Phase)`.
"""
abstract type Phase end

"""
    abstract type AbstractTrainingPhase <: Phase

An abstract type for phases where parameter updates
are being made. This exists so callbacks can dispatch on it and work
with custom training phases.

The default implementation is [`TrainingPhase`](#).
"""
abstract type AbstractTrainingPhase <: Phase end

"""
    TrainingPhase()

A regular training phase. It iterates over batches
in `learner.data.training` and updates the model parameters
using `learner.optim` after calculating the gradients.

Throws the following events: [`EpochBegin`](#), [`BatchBegin`](#),
[`LossBegin`](#), [`BackwardBegin`](#), [`BackwardEnd`](#), [`BatchEnd`](#),
[`EpochEnd`](#).
"""
struct TrainingPhase <: AbstractTrainingPhase end

"""
    TrainingPhase()

A regular validation phase. It iterates over batches
in `learner.data.validation` and performs a forward pass.

Throws the following events: [`EpochBegin`](#), [`BatchBegin`](#),
[`LossBegin`](#), [`BatchEnd`](#), [`EpochEnd`](#).
"""
struct ValidationPhase <: Phase end
struct TestPhase <: Phase end

export
    Phase,
    AbstractTrainingPhase,
    TrainingPhase,
    ValidationPhase,
    TestPhase
end # module

using .Phases