text:
export histogram

const histogram_plot_types = [BarPlot, Heatmap, Volume]

function convert_arguments(P::Type{<:AbstractPlot}, h::StatsBase.Histogram{<:Any, N}) where N
    ptype = plottype(P, histogram_plot_types[N])
    f(edges) = edges[1:end-1] .+ diff(edges)./2
    kwargs = N == 1 ? (; width = step(h.edges[1])) : NamedTuple()
    to_plotspec(ptype, convert_arguments(ptype, map(f, h.edges)..., Float64.(h.weights)); kwargs...)
end

to_weights(v) = StatsBase.weights(v)
to_weights(v::StatsBase.AbstractWeights) = v

bool2options(t::Bool, a, b) = t ? a : b
bool2options(t, a, b) = t

function _histogram(args...; edges = automatic, weights = automatic, normalize = false, kwargs...)
    ea = edges === automatic ? () : (to_tuple(edges),)
    wa = weights === automatic ? () : (to_weights(weights),)
    ha = length(args) == 1 ? args[1] : args
    attr = Dict(kwargs)
    isempty(ea) || pop!(attr, :nbins, nothing)
    h = fit(StatsBase.Histogram, to_tuple(ha), wa..., ea...; attr...)
    norm_option = bool2options(normalize, :pdf, :none)
    StatsBase.normalize(h, mode = norm_option)
end

const histogram = Analysis(_histogram)

function adjust_globally(hist::Analysis{typeof(_histogram)}, traces)
    (length(traces) == 1 || get(hist.kwargs, :edges, automatic) !== automatic) && return hist
    global_output = map(concatenate, (trace.output for trace in traces)...)
    h = hist(global_output...)
    Analysis(hist.f; merge(hist.kwargs, (edges = h.edges,))...)
end


text:
"""
A Julia package for solving large-scale electromagnetic scattering problems in
two dimensions; specifically, those containing a large number of penetrable
smooth particles. Provides the ability to optimize over the particle parameters
for various design problems.
"""
module ParticleScattering

#Core functionality
using SpecialFunctions, IterativeSolvers, LinearMaps, Optim, SparseArrays, LinearAlgebra
import Statistics: mean
import LineSearches
#For plotting with PyPlot - actual import is done at runtime
using PyPlot, PyCall
const patch = PyNULL()
function __init__()
    copy!(patch, pyimport("matplotlib.patches")) #circles, polygons
end

include("PS_types.jl")
include("shapes.jl")
include("scattering.jl")
include("multipole.jl")
include("visualization.jl")
include("fmm_matrices.jl")
include("fmm_mvp.jl")
include("fmm_main.jl")
include("optimize_phis.jl")
include("optimize_rs.jl")
include("utilities.jl")
include("minimum_N_P.jl")
include("incident.jl")

#methods, shapes.jl
export rounded_star, squircle, ellipse, square_grid, rect_grid, hex_grid,
    randpoints, luneburg_grid, verify_min_distance
#methods, multipole.jl
export solve_particle_scattering, scattered_field_multipole
#methods, minimum_N_P.jl
export minimumN, minimumP
#methods, visualization.jl
export plot_far_field, plot_near_field, calc_near_field, draw_shapes
#methods, fmm_main.jl
export solve_particle_scattering_FMM
#methods, scattering.jl
export get_potential, get_potentialPW, scatteredfield
#types, PS_types.jl
export ScatteringProblem, OptimBuffer, FMMoptions, R_multipole,
    ShapeParams, CircleParams, AbstractShapeParams, PlaneWave, LineSource,
    CurrentSource, Einc
#methods, utilities.jl
export find_border, uniqueind
#methods, optimize_phis.jl
export optimize_φ
#methods, optimize_rs.jl
export optimize_radius
#methods, incident.jl
export uinc, hxinc, hyinc
#consts, incident.jl
export eta0

#temp
export divideSpace, FMMtruncation, particleExpansion, FMMbuildMatrices
export FMMbuffer, FMMmatrices

include("optimize_rs_mf.jl")
export optimize_rs_mf

include("poynting.jl")
export calc_power

include("optimize_rs_pwr.jl")
export PowerBuffer, OptimProblemBuffer, optMatrixPwr, optimize_pwr_rs_f, optimize_pwr_rs_g!
include("optimize_phis_pwr.jl")
export optimize_pwr_φ_f, optimize_pwr_φ_g!
end

ext: jl


text: # LU factorization for unsymmetric systems


text:
using Test, Gantner, Dates

name = joinpath(@__DIR__, "..", "example", "data.dat")

@testset "gantnerinfo" begin
    # check gantnerinfo
    nc, nv, fs, chanlegendtext, starttime, finishtime = gantnerinfo(name)
    @test nc == 5
    @test nv == 200_000
    @test fs == 10_000.0
    @test chanlegendtext == ["Trig Condition", "DI_01a", "FI_01a", "2nd MUX AI_01", "1st MUX AI_03"]
    #@test isapprox(starttime, unix2datetime(ctime(name)), atol=1e-3) # isapprox not valid with these types
end

@testset "gantnerread" begin
    #check gantnerread for all channels
    (t, data, fs, chanTextLegend) = gantnerread(name)
    @test t == 0.0:0.0001:19.9999
    @test data[2,5] == 0.001806761370971799

    #check gantnerread for all channels with scaling array
    scalefactors = [1.0; 2.0; 3.0; 4.0; 5.0]
    (t, data, fs, chanTextLegend) = gantnerread(name, scale=scalefactors)
    @test t == 0.0:0.0001:19.9999
    @test data[2,5] == 5.0 * 0.001806761370971799

    #check gantnerread for all channels with scaling Float64
    scalefactors = 5.0
    (t, data, fs, chanTextLegend) = gantnerread(name, scale=scalefactors)
    @test t == 0.0:0.0001:19.9999
    @test data[2,5] == 5.0 * 0.001806761370971799

    # check gantnerread where the sampling time is read, rather than recreated
    (t, data, fs, chanTextLegend) = gantnerread(name; lazytime = false)
    @test t[1] == 44210.79703523496

    # check gantnerread, where only one channel is read and scale factor is used
    (t, data, fs, chanTextLegend) = gantnerread(name, 5, scale=2.0)
    @test typeof(data) == Vector{Float64}
    @test length(data) == 200_000
    @test data[5] == 2 * -0.004325194749981165

    # check gantnerread, where a number of channels are read and scale factor is used
    (t, data, fs, chanTextLegend) = gantnerread(name, 2:4, scale=2.0)
    @test typeof(data) == Matrix{Float64}
    @test size(data) == (200_000, 3)
    @test data[5,3] == 0.028291847556829453
end

@testset "gantnermask" begin
    a = [0.0,3,0,3,0,3,0,3,0,3,0,3,0,3,0,3,0,3,0,3,0,3,0,3,0,3,0,3,0,3,0,3,0,3,0,3]
    b = gantnermask(a,2)
    @test typeof(a) == typeof(b)
    @test all(a .>= b)
    @test all(a .== 3 .* b)
    @test all(gantnermask(a, 2) .== gantnermask(a, 1))
    c = gantnermask(a,1,1)
    @test all(b .== c)
    @test all(gantnermask(a, 2, 5) .== gantnermask(a, 1, 5))
    @test all(gantnermask(a, 2, 3) .== repeat(0:1, inner=3, outer=6))

    c = [0.0, 1, 2, 4, 8, 16, 32]
    @test all(gantnermask(c, 1) .== [0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    @test all(gantnermask(c, 2) .== [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0])
    @test all(gantnermask(c, 3) .== [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0])
    @test all(gantnermask(c, 4) .== [0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0])
    @test all(gantnermask(c, 5) .== [0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0])
    @test all(gantnermask(c, 6) .== [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0])
    @test all(gantnermask(c, 7) .== [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    @test all(gantnermask(c, 8) .== [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    @test all(gantnermask(c, 9) .== [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    @test all(gantnermask(c, 10) .== [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])

    @test all(gantnermask(c, 1, 1) .== [0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    @test all(gantnermask(c, 2, 1) .== [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0])
    @test all(gantnermask(c, 3, 1) .== [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0])
    @test all(gantnermask(c, 4, 1) .== [0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0])
    @test all(gantnermask(c, 5, 1) .== [0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0])
    @test all(gantnermask(c, 6, 1) .== [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0])
    @test all(gantnermask(c, 7, 1) .== [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    @test all(gantnermask(c, 8, 1) .== [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    @test all(gantnermask(c, 9, 1) .== [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    @test all(gantnermask(c, 10, 1) .== [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
end


text:
@testset "cosineknn" begin
    A = load_matrix_network("bfs_example")
    S = cosineknn(A,2)

    # A = speye(Int64,4)
    A = sparse(1.0I,4,4)
    A[4,1] = 1
    CKNN = cosineknn(A,2)
    @test ((CKNN[4,1] *2)^2 - 2) <= 1e-7

    # A = speye(Int64,4)
    A = sparse(1.0I,4,4)
    A[4,1] = 1
    CKNN = cosineknn(MatrixNetwork(A),2)
    @test ((CKNN[4,1] *2)^2 - 2) <= 1e-7
end


text:
#=
Model/data deviations for calibration

Intended workflow:
* Set up deviations while defining model. Load in the data.
* After solving the model: fill in model values.
* Compute deviations for calibration.
* Report deviations.
* Show model fit. (All the last steps can simply work of the deviation vector)

# Todo

* Make an iterator for deviations
=#

# Numeric values are stored as this type (in containers where type is unknown; e.g. Dict)
const DevType = Float64

## -----------  Types

"""
	AbstractDeviation

Abstract type for Deviation objects. A Deviation object implements:

    * `set_model_values`
    * `scalar_dev`

Contains:

    * `stdV`: std errors of data. Defaults to []. Used when showing deviations.
    * `wtV`: weights within a deviation (e.g. 1/std error)
    * `scalarWt`: used by `DevVector` to weight each `scalar_dev`
    * `showFct`: function that takes an `AbstractDeviation` as input and produces a model/data comparison
    * `showPath`: file name where `showFct` files its output. `stdout` is used if empty.
"""
abstract type AbstractDeviation{F1} end

function Base.isempty(d :: AbstractDeviation)
    return d.name == :empty
end

"""
    Deviation

Holds numeric arrays. The default for deviations.

Contains model values, data values, data std deviations (for display only), weights (for constructing scalar deviation).

Can specify that only `modelV[idxV...]` is to be compared with `dataV`.
"""
Base.@kwdef mutable struct Deviation{F1 <: AbstractFloat} <:
    AbstractDeviation{F1}

    name  :: Symbol     # eg 'fracEnterIq'
    modelV  :: Array{F1} = zeros(F1, 1)  # model values
    dataV  :: Array{F1} = zeros(F1, 1)   # data values
    stdV :: Array{F1} = zeros(F1, 1)
    # relative weights, sum to user choice
    wtV  :: Array{F1} = ones(F1, size(dataV))
    # Indices such that `modelV[idxV...]` matches `dataV`
    # Default is to use all
    idxV :: Vector{Any} = []
    # Scales the scalar deviation.
    scalarWt :: F1 = one(F1)
    # Using weighted normp. By default: sum of abs deviations.
    normP :: F1 = one(F1)
    shortStr  :: String = String(name)      # eg 'enter/iq'
    # eg 'fraction entering college by iq quartile'
    longStr  :: String = shortStr
    # For displaying the deviation. Compatible with `Formatting.sprintf1`
    # E.g. "%.2f"
    fmtStr  :: String = "%.2f"
    showFct = deviation_show_fct
    showPath :: String = ""
end

"""
    ScalarDeviation

Here the `wtV` field is intended to hold 1 / std error of the moment.
The scalar deviation is always the absolute deviation. For scalars, the exponent in the norm makes no difference.
"""
Base.@kwdef mutable struct ScalarDeviation{F1 <: AbstractFloat} <: AbstractDeviation{F1}
    name  :: Symbol     # eg 'fracEnterIq'
    modelV  :: F1 = zero(F1)  # model values
    dataV  :: F1 = zero(F1)   # data values
    stdV :: F1 = zero(F1)
    # Used when a std error of the data moment is known
    wtV :: F1 = one(F1)
    scalarWt :: F1 = one(F1)
    # Using weighted normp. By default: sum of abs deviations.
    normP :: F1 = one(F1)
    shortStr  :: String = String(name)      # eg 'enter/iq'
    longStr  :: String = shortStr
    fmtStr  :: String = "%.2f"
    showFct = scalar_show_fct
    showPath :: String = ""
end

"""
	RegressionDeviation

Holds model and data in the form of `RegressionTable` objects.
Some regressors may be excluded from the computation of the scalar deviation.
"""
Base.@kwdef mutable struct RegressionDeviation{F1} <: AbstractDeviation{F1}
    name  :: Symbol
    modelV  :: RegressionTable = RegressionTable()
    dataV  :: RegressionTable = RegressionTable()
    excludeV :: Vector{Symbol} = Vector{Symbol}()
    scalarWt :: F1 = one(F1)
    normP :: F1 = one(F1)
    shortStr  :: String = String(name)      # eg 'enter/iq'
    longStr  :: String = shortStr
    fmtStr  :: String = "%.2f"
    showFct = regression_show_fct
    showPath :: String = ""
end

"""
	$(SIGNATURES)

Bounds deviation. Returns zero scalar deviation until model values get out of bounds.
"""
Base.@kwdef mutable struct BoundsDeviation{F1 <: AbstractFloat} <: AbstractDeviation{F1}
    name  :: Symbol
    modelV  :: Array{F1} = zeros(F1, 1)
    # Bounds
    lbV  :: Array{F1} = zeros(F1, 1)
    ubV  :: Array{F1} = ones(F1, 1)
    # relative weights, sum to user choice
    wtV  :: Array{F1} = ones(F1, size(lbV))
    scalarWt :: F1 = one(F1)
    shortStr  :: String = String(name)
    # eg 'fraction entering college by iq quartile'
    longStr  :: String = shortStr
    # For displaying the deviation. Compatible with `Formatting.sprintf1`
    # E.g. "%.2f"
    fmtStr  :: String = "%.2f"
    showFct = bounds_show_fct
    showPath :: String = ""
end

"""
	$(SIGNATURES)

Penalty deviation. Calls a function on model values to return scalar deviation.
"""
Base.@kwdef mutable struct PenaltyDeviation{F1 <: AbstractFloat} <: AbstractDeviation{F1}
    name  :: Symbol
    modelV  :: Array{F1} = zeros(F1, 1)
    scalarDevFct :: Function
    scalarWt :: F1 = one(F1)
    shortStr  :: String = String(name)
    longStr  :: String = shortStr
    showFct = penalty_show_fct
    showPath :: String = ""
end

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


text:
#### ODE Steppers for y'(t) = f(t, y(t)) #####

##### Secondary Methods ######

function Derivative(x, f, d=1, h=10.0^-8)
    """
    Estimate the derivative of f at x using the center
    finite difference technique with step size h.

    x and f can be vector valued, in which case
    d must be provided. d should be a list of only
    zeros and ones.

    For example, if x and f(x) are vectors of length four
    and we want the derivative with respect to the second
    component in x, then we call Derivative(x, f, [0, 1, 0, 0]).
    """
    (f(x+h*d) - f(x-h*d)) / (2*h)
end

function NewtonsMethod(x, f, error_tol=10.0^-10)
    """
    Estimate the root of f starting from x.

    x and f(x) are vector valued. Newton's method becomes

    \vec{x}_{n+1} = \vec{x}_{n} - inv(Jacobian) * \vec{f(\vec{x}_n)}

    I create the Jacobian by estimating the derivatives
    using a center finite difference method.
    """
    h, n = 10.0^-5, length(x)
    J = zeros(n, n)
    for j=1:n
        d = zeros(n)
        d[j] = 1
        J[:, j] = Derivative(x, f, d)
    end
    count = 0
    while norm(f(x)) > error_tol
        count += 1
        x -= inv(J) * f(x)
        if count > 50
            println("Newton's Method took unusually long so it was terminated prematurely.")
            return x
        end
    end
    x
end

##### End Secondary Methods #####

###### Single Steppers ######

## Step functions that do not include an error estimate ##
# Call these with Stepper(f, method, t0, y0, tf, h)

function EulerStep(f, t, y, h)
    h * f(t, y)
end

function ImprovedEulerStep(f, t, y, h)
    f0 = f(t, y)
    (h / 2.0) * (f0 + f(t + h, y + h * f0))
end

function BackwardEulerStep(f, t, y, h)
    """ Implicit Backward Euler Step """
    #Use Euler step as starting point for NewtonsMethod.
    NewtonsMethod(f(t, y), f1 -> h*f(t+h, y+f1) - f1)
end

function GL2Step(f, t, y, h)
    """ Implicit Gauss-Legendre Order 2 Method """
    #Use Euler step as starting point for NewtonsMethod.
    NewtonsMethod(f(t, y), f1 -> h*f(t+h/2.0, y+f1/2.0) - f1)
end

function MidpointStep(f, t, y, h)
    h * f(t + h / 2, y + (h / 2) * f(t, y))
end

function CrankNicholsonStep(f, t, y, h)
    """ Second order implicit """
    f1 = f(t, y)
    NewtonsMethod(f1, f2 -> (h/2.0) * (f1 + f(t+h, y+f2)) - f2)
end

function RK4Step(f, t, y, h)
    f1 = f(t, y)
    f2 = f(t + h / 2, y + (h / 2) * f1)
    f3 = f(t + h / 2, y + (h / 2) * f2)
    f4 = f(t + h, y + h * f3)
    (h / 6) * (f1 + f4 + 2 * (f2 + f3))
end

function Stepper(f, method, t0, y0, tf, h)
    """
    The step function `method` returns just the step,
    not an error estimate.
    """
    ts, ys = Float64[], Vector{Float64}[]
    while t0 <= tf
        push!(ts, t0)
        push!(ys, y0)
        y0 += method(f, t0, y0, h)
        t0 += h
    end
    ts, ys
end

## End step functions that do not include an error estimate ##

## Step functions that include an error estimate ##
# Call these with AdaptiveStepper(f, method, t0, y0, tf)

function SimpleErrorStep(f, t, y, h)
    """ Compare implicit Crank Nicholson to Midpoint to get error """
    step = CrankNicholsonStep(f, t, y, h)
    step, norm(step - MidpointStep(f, t, y, h))
end

function RKTrapStep(f, t, y, h)
    """ Implicit Runge-Kutta Trapazoidal Method """
    f1 = h * f(t, y)
    f2 = NewtonsMethod(f1, f2 -> h*f(t+h, y+f1/2.0+f2/2.0) - f2)

    soln = f1/2.0 + f2/2.0
    soln, norm(soln - f1)
end

function DopriStep(f, t, y, h)
    """ Dormand-Prince Method """
    f1 = h * f(t, y)
    f2 = h * f(t + h/5.0, y + f1/5.0)
    f3 = h * f(t + 3*h/10.0, y + 3*f1/40.0 + 9*f2/40.0)
    f4 = h * f(t + 4*h/5.0, y + 44*f1/45.0 - 56*f2/15.0 + 32*f3/9.0)
    f5 = h * f(t + 8*h/9.0, y + 19372*f1/6561.0 - 25360*f2/2187.0 + 64448*f3/6561.0 - 212*f4/729.0)
    f6 = h * f(t + h, y + 9017*f1/3168.0 - 355*f2/33.0 + 46732*f3/5247.0 + 49*f4/176.0 - 5103*f5/18656.0)
    soln = 35*f1/384.0 + 500*f3/1113.0 + 125*f4/192.0 - 2187*f5/6784.0 + 11*f6/84.0
    f7 = h * f(t + h, y + soln)
    alt = 5179*f1/57600.0 + 7571*f3/16695.0 + 393*f4/640.0 - 92097*f5/339200.0 + 187*f6/2100.0 + f7/40.0

    soln, norm(soln - alt)
end

function CarpStep(f, t, y, h)
    """ Cash-Karp Method """
    f1 = h * f(t, y)
    f2 = h * f(t + h/5.0, y + f1/5.0)
    f3 = h * f(t + 3*h/10, y + 3*f1/40.0 + 9*f2/40.0)
    f4 = h * f(t + 3*h/5.0, y + 3*f1/10.0 - 9*f2/10.0 + 6*f3/5.0)
    f5 = h * f(t + h, y - 11*f1/54.0 + 5*f2/2.0 - 70*f3/27.0 + 35*f4/27.0)
    f6 = h * f(t + 7*h/8.0, y + 1631*f1/55296 + 175*f2/512 + 575*f3/13824 + 44275*f4/110592.0 + 253*f5/4096.0)
    soln = 37*f1/378.0 + 250*f3/621.0 + 125*f4/594.0 + 512*f6/1771.0
    alt = 2825*f1/27648.0 + 18575*f3/48384.0 + 13525*f4/55296.0 + 277*f5/14336.0 + f6/4.0

    soln, norm(soln - alt)
end

function RKFStep(f, t, y, h)
    """ Runge-Kutta-Fehlberg Method """
    f1 = h * f(t, y)
    f2 = h * f(t + h/4.0, y + f1/4.0)
    f3 = h * f(t + 3*h/8.0, y + 3*f1/32.0 + 9*f2/32.0)
    f4 = h * f(t + 12*h/13.0, y + 1932*f1/2197.0 - 7200*f2/2197.0 + 7296*f3/2197.0)
    f5 = h * f(t + h, y + 439*f1/216.0 - 8*f2 + 3680*f3/513.0 - 845*f4/4104.0)
    f6 = h * f(t + h/2.0, y - 8*f1/27.0 + 2*f2 - 3544*f3/2565.0 + 1859*f4/4104.0 - 11*f5/40.0)
    soln = 16*f1/135.0 + 6656*f3/12825.0 + 28561*f4/56430.0 - 9*f5/50.0 + 2*f6/55.0
    alt = 25*f1/216.0 + 1408*f3/2565.0 + 2197*f4/4104.0 - f5/5.0

    soln, norm(soln - alt)
end

function AdaptiveStepper(f, method, t0, y0, tf, error_tol=10.0^-7)
    """
    The step function `method` must return the tuple,
    (step, error estimate).
    """
    ts, ys = Float64[t0], Vector{Float64}[y0]
    h, totalerror = 0.001, 0
    while t0 < tf
        if h < 10.0^-8 && h < tf - t0
            println("Step size effectively zero at t = ", t0)
            h = 0.001
        else
            if h > tf - t0
                h = tf - t0
            end
            m, error_est = method(f, t0, y0, h)
            if error_est > error_tol
                h *= .75
            else
                t0 += h
                y0 += m
                push!(ts, t0)
                push!(ys, y0)
                if error_est < error_tol / 10.0
                    h *= 1.2
                end
                totalerror += error_est
            end
        end
    end
    println("Estimated upper bound on total error of ", method, ": ", totalerror)
    ts, ys
end

## End step functions that include an error estimate ##

##### End Single Steppers #####

##### Multi Steppers ######
#Call these with MultiStepper(f, method, t0, y0, tf, h)

function AdamsBash2(f, t, y, h)
    """ Adams–Bashforth Order 2 explicit method """
    y0n = Array[y, y + RK4Step(f, t, y, h)]
    function g(t, yn)
        #yn = Array[y0, y1]
        y0, y1 = yn
        y2 = y1 + h * (1.5 * f(t + h, y1) - 0.5 * f(t, y0))
        Array[y1, y2]
    end
    g, y0n
end

function AdamsBash3(f, t, y, h)
    """ Adams–Bashforth Order 3 explicit method """
    y1 = y + RK4Step(f, t, y, h)
    y2 = y1 + RK4Step(f, t+h, y1, h)
    y0n = Array[y, y1, y2]
    function g(t, yn)
        #yn = Array[y0, y1, y2]
        y0, y1, y2 = yn
        y3 = y2 + h * (23.0 / 12.0 * f(t+2*h, y2) - 4.0 / 3.0 * f(t+h, y1) + 5.0 / 12.0 * f(t, y0))
        Array[y1, y2, y3]
    end
    g, y0n
end

function AdamsBash4(f, t, y, h)
    """ Adams–Bashforth Order 4 explicit method """
    y1 = y + RK4Step(f, t, y, h)
    y2 = y1 + RK4Step(f, t+h, y1, h)
    y3 = y2 + RK4Step(f, t+2*h, y2, h)
    y0n = Array[y, y1, y2, y3]
    function g(t, yn)
        #yn = Array[y0, y1, y2, y3]
        y0, y1, y2, y3 = yn
        y4 = y3+h * (55.0 / 24.0 * f(t+3*h, y3) - 59.0 / 24.0 * f(t+2*h, y2) + 37.0 / 24.0 * f(t+h, y1) - 3.0 / 8.0 * f(t, y0))
        Array[y1, y2, y3, y4]
    end
    g, y0n
end

function AdamsMoult2(f, t, y, h)
    """ Adams-Moulton Order 2 implicit method """
    y1 = y + RK4Step(f, t, y, h)
    y0n = Array[y, y1]
    function g(t, yn)
        y0, y1 = yn
        f0, f1 = f(t, y0), f(t+h, y1)
        y2 = NewtonsMethod(y1, y2 -> y1+h*(5*f(t+2*h, y2)/12.0 + 2*f1/3.0 - f0/12.0) - y2)
        Array[y1, y2]
    end
    g, y0n
end

function AdamsMoult3(f, t, y, h)
    """ Adams-Moulton Order 3 implicit method """
    y1 = y + RK4Step(f, t, y, h)
    y2 = y1 + RK4Step(f, t+h, y1, h)
    y0n = Array[y, y1, y2]
    function g(t, yn)
        y0, y1, y2 = yn
        f0, f1, f2 = f(t, y0), f(t+h, y1), f(t+2*h, y2)
        y3 = NewtonsMethod(y2, y3 -> y2+h*(3*f(t+3*h, y3)/8.0 + 19*f2/24.0 - 5*f1/24.0 + f0/24.0)-y3)
        Array[y1, y2, y3]
    end
    g, y0n
end

function AdamsMoult4(f, t, y, h)
    """ Adams-Moulton Order 4 implicit method """
    y1 = y + RK4Step(f, t, y, h)
    y2 = y1 + RK4Step(f, t+h, y1, h)
    y3 = y2 + RK4Step(f, t+2*h, y2, h)
    y0n = Array[y, y1, y2, y3]
    function g(t, yn)
        y0, y1, y2, y3 = yn
        f0, f1, f2, f3 = f(t, y0), f(t+h, y1), f(t+2*h, y2), f(t+3*h, y3)
        y4 = NewtonsMethod(y3, y4 -> y3+h*(251*f(t+4*h, y4)/720.0+646*f3/720.0 - 264*f2/720.0 + 106*f1/720.0 - 19*f0/720.0)-y4)
        Array[y1, y2, y3, y4]
    end
    g, y0n
end

function MultiStepper(f, method, t0, y0, tf, h)
    """
    method returns g(t, yn), y0n with y0n is
    an array  of previous values of y.
    """
    g, y0n = method(f, t0, y0, h)
    #y0n = Array[y0, y1, ..., yn]
    ts, ys = Float64[], Vector{Float64}[]
    while t0 <= tf + h
        push!(ts, t0)
        push!(ys, y0n[1])
        y0n = g(t0, y0n)
        t0 += h
    end
    ts, ys
end

##### End Multi Steppers #####

#### Test ####
function func(t, y)
    y1, ydot = y
    [ydot, -t*y1]
end

t, tmax, h = 0.0, 20, 0.01
y0 = [-1.0, 0.0]

#xs, ys = MultiStepper(func, AdamsMoult4, t, copy(y0), tmax, h)
xs, ys = Stepper(func, GL2Step, t, copy(y0), tmax, h)
#xs, ys = AdaptiveStepper(func, RKTrapStep, t, copy(y0), tmax, h)
figure(0)
plot(xs, [ys[i][1] for i=1:length(ys)], "o-", label="Position")
plot(xs, [ys[i][2] for i=1:length(ys)], "x--", label="Velocity")
legend()

xs, ys = AdaptiveStepper(func, RKFStep, t, copy(y0), tmax)
figure(1)
plot(xs, [ys[i][1] for i=1:length(ys)], "o-", label="Position")
plot(xs, [ys[i][2] for i=1:length(ys)], "x--", label="Velocity")
legend()

max_stars_repo_name: jiosue/QuadratureCandJulia

max_issues_repo_name: jiosue/QuadratureCandJulia

max_forks_repo_name: jiosue/QuadratureCandJulia


text:
module availability

export Supply, Demand, calculate_Avail, applySupplyAdjusters

include("model/supply.jl")
include("model/demand.jl")
include("serialization.jl")
include("adjuster.jl")
include("exception/exceptions.jl")
# include("router.jl")

function calculate_Avail(supply::Supply, demand::Demand)::Int128
    @info ("calculate availability")
    supply.quantity - demand.quantity
end

end # module

max_stars_repo_name: kawadesumit/availability

max_issues_repo_name: kawadesumit/availability

max_forks_repo_name: kawadesumit/availability


text:
# options checking specific to the Euler Equations

"""
  This function checks the options in the options dictionary after default
  values have been supplied and throws exceptions if unsupported/incompatable
  options are specified

  **Inputs**

   * opts: options dictionary
   * comm: MPI communicator

  **Outputs**

    none
"""
function checkOptions(opts, comm=MPI.COMM_WORLD)

  # this is a problem for nested physics
#  if opts["physics"] != PhysicsName
#    error("physics not specified as $PhysicsName, are you lost?")
#  end

  if opts["use_edge_res"]
    error("use_edge_res does not work correctly, do not use")
  end

  if opts["use_staggered_grid"]
    if opts["volume_integral_type"] != 2
      error("cannot use staggered grids with non entropy stable volume integrals")
    end
#=
    if opts["face_integral_type"] != 2
      error("cannot use staggered grids with non entropy stable face integrals")
    end
=#
    if opts["operator_type2"] == "SBPNone"
      error("must specify a real SBP operator for staggered grid")
    end
  end

  run_type = opts["run_type"]
  #TODO: what to do for jac_type 4?
  if opts["use_DG"] &&
     opts["Flux_name"] in keys(FluxDict_diff) &&
     (run_type == 5 || run_type == 40 || run_type == 20 || run_type == 50) &&
     opts["jac_type"] >= 2 #&&
#     !opts["use_src_term"]
    get!(opts, "calc_jac_explicit", true)
  else
    get!(opts, "calc_jac_explicit", false)
    if opts["calc_jac_explicit"]
      error("cannot calculate jacobian explicitly with given options")
    end
  end

  if opts["calc_jac_explicit"]
    get!(opts, "preallocate_jacobian_coloring", false)
  else
    get!(opts, "preallocate_jacobian_coloring", true)
  end

  get!(opts, "homotopy_function", "FirstOrderDissipation")
  get!(opts, "freeze_viscosity", false)
  get!(opts, "shock_capturing_variables", "IR")
  get!(opts, "sensor_pp_use_filtered", false)
  get!(opts, "sensor_pp_s0_offset", 0)

  return nothing
end

import PDESolver: _getSparsityPattern

function _getSparsityPattern(mesh::AbstractMesh, sbp::AbstractOperator,
                             eqn::EulerData, opts)

  if !opts["addShockCapturing"]
    return INVISCID
  else
    return getSparsityPattern(eqn.shock_capturing)
  end

end


text:
mutable struct Graph
    num_levels::Int64
    acts
    μacts
    props
    μprops
    states
    leveled
    initprops
    goalprops
    function Graph()
        new(0, Dict(0=>[]), Dict(0=>[]),  Dict(0=>[]), Dict(0=>[]),
            Dict(), false, [], [])
    end
end

mutable struct NoOp
    name
    args
    pos_prec
    neg_prec
    pos_eff
    neg_eff
    function NoOp(prop)
        new(:NoOp, [], [prop], [], [prop], [])
    end
end

mutable struct Action
    name
    args
    pos_prec
    neg_prec
    pos_eff
    neg_eff
    function Action(name)
        new(name, [], [], [], [], [])
    end
end


text:
using JSON
using JLD2

module SemanticWords
import SimilaritySearch: L2Distance, CosineDistance, DenseCosine, Sequential, KnnResult
import NearNeighborGraph: optimize!, fit!, LocalSearchIndex, BeamSearch, FixedNeighborhood, LogNeighborhood, LogSatNeighborhood

include("codebook.jl")

#export compute_cluster, read_vocabulary

#function compute_cluster(X, numcenters, k, kind; maxiter=10, tol=0.001, recall=0.9)
#    if kind == "random"
#        centroids, codes, rlist = codebook(X, numcenters, k, maxiter=maxiter, tol=tol, recall=recall)
#    elseif kind == "fft"
#        centroids, codes, rlist = fftcodebook(X, numcenters, k)
#    else # kind == "dnet"
#        centroids, codes, rlist = dnetcodebook(X, numcenters, k)
#    end
#
#    centroids, codes
#end

end


text:
"Your optional docstring here"
function distance(a, b)
    if length(a) != length(b)
        throw(ArgumentError("Lengths of the sequencess differ"))
    end

    count = 0
    for i in 1:length(a)
        if a[i] != b[i]
            count += 1
        end
    end
    return count
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
using Reexport

@reexport module GridEstimators
    import ..isboundarycondition
    import ..invariantmeasure, InvariantDistribution
    import ..transferoperator
    import ..GridBasedTransferOperator

    include("binvisits.jl")
    include("SingleGrid.jl")

end


text:
module likelihoods

using ..gpflow
import ..gpflow: instantiate!
export
    Bernoulli,
    Beta,
    Exponential,
    Gamma,
    Gaussian,
    Likelihood,
    MonteCarloLikelihood,
    MultiClass,
    Ordinal,
    Poisson,
    RobustMax,
    SoftMax,
    StudentT,
    SwitchedLikelihood,
    instantiate!

mutable struct Bernoulli<:AbstractLikelihood
    # invlink # TODO: add capability to add custom invlink
    o::Union{PyObject, Nothing}
end

function Bernoulli()
    out = Bernoulli(nothing)
    instantiate!(out)
    return out
end

function instantiate!(o::Union{Bernoulli, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.likelihoods.Bernoulli()
end

"""
    This uses a reparameterisation of the Beta density. We have the mean of the
    Beta distribution given by the transformed process:

        m = sigma(f)

    and a scale parameter. The familiar alpha, beta parameters are given by

        m     = alpha / (alpha + beta)
        scale = alpha + beta

    so:
        alpha = scale * m
        beta  = scale * (1-m)
"""
mutable struct Beta{T}<:AbstractLikelihood
    # invlink # TODO: add capability to add custom invlink
    scale::T
    o::Union{PyObject, Nothing}
end

function Beta(;scale::Real=1.0)
    out = Beta(scale, nothing)
    instantiate!(out)
    return out
end

function instantiate!(o::Union{Beta, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.likelihoods.Beta(;scale=o.scale)
end

mutable struct Exponential<:AbstractLikelihood
    # invlink # TODO: add capability to add custom invlink
    o::Union{PyObject, Nothing}
end

function Exponential()
    out = Exponential(nothing)
    instantiate!(out)
    return out
end

function instantiate!(o::Union{Exponential, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.likelihoods.Exponential()
end

"""
    Use the transformed GP to give the *scale* (inverse rate) of the Gamma
"""
mutable struct Gamma<:AbstractLikelihood
    # invlink # TODO: add capability to add custom invlink
    o::Union{PyObject, Nothing}
end

function Gamma()
    out = Gamma(nothing)
    instantiate!(out)
    return out
end

function instantiate!(o::Union{Gamma, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.likelihoods.Gamma()
end

mutable struct Gaussian{T}<:AbstractLikelihood
    variance::T
    o::Union{PyObject, Nothing}
end

function Gaussian(;variance::T=1.0) where T <: Real
    out = Gaussian(variance, nothing)
    instantiate!(out)
    return out
end

function instantiate!(o::Union{Gaussian, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.likelihoods.Gaussian(;variance=o.variance)
end

# TODO GaussianMC not implemented in GPFlow
"""
    A likelihood for doing ordinal regression.

    The data are integer values from 0 to K, and the user must specify (K-1)
    'bin edges' which define the points at which the labels switch. Let the bin
    edges be [a_0, a_1, ... a_{K-1}], then the likelihood is

    p(Y=0|F) = phi((a_0 - F) / sigma)
    p(Y=1|F) = phi((a_1 - F) / sigma) - phi((a_0 - F) / sigma)
    p(Y=2|F) = phi((a_2 - F) / sigma) - phi((a_1 - F) / sigma)
    ...
    p(Y=K|F) = 1 - phi((a_{K-1} - F) / sigma)

    where phi is the cumulative density function of a Gaussian (the inverse probit
    function) and sigma is a parameter to be learned. A reference is:

    @article{chu2005gaussian,
      title={Gaussian processes for ordinal regression},
      author={Chu, Wei and Ghahramani, Zoubin},
      journal={Journal of Machine Learning Research},
      volume={6},
      number={Jul},
      pages={1019--1041},
      year={2005}
    }
"""
function GaussianMC()
    throw("GaussianMC not implemented")
end

mutable struct Likelihood<:AbstractLikelihood
    o::Union{PyObject, Nothing}
end

function Likelihood()
    out = Likelihood(nothing)
    instantiate!(out)
    return out
end

function instantiate!(o::Union{Likelihood, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.likelihoods.Likelihood()
end

mutable struct MonteCarloLikelihood<:AbstractLikelihood
    o::Union{PyObject, Nothing}
end

function MonteCarloLikelihood()
    out = MonteCarloLikelihood(nothing)
    instantiate!(out)
    return out
end

function instantiate!(o::Union{MonteCarloLikelihood, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.likelihoods.MonteCarloLikelihood()
end

mutable struct MultiClass<:AbstractLikelihood
    num_classes::Integer
    # invlink # TODO: add capability to add custom invlink
    o::Union{PyObject, Nothing}
end

function MultiClass(num_classes::Integer)
    out = MultiClass(num_classes, nothing)
    instantiate!(out)
    return out
end

function instantiate!(o::Union{MultiClass, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.likelihoods.MultiClass(o.num_classes)
end

"""
    A likelihood for doing ordinal regression.

    The data are integer values from 0 to K, and the user must specify (K-1)
    'bin edges' which define the points at which the labels switch. Let the bin
    edges be [a_0, a_1, ... a_{K-1}], then the likelihood is

    p(Y=0|F) = phi((a_0 - F) / sigma)
    p(Y=1|F) = phi((a_1 - F) / sigma) - phi((a_0 - F) / sigma)
    p(Y=2|F) = phi((a_2 - F) / sigma) - phi((a_1 - F) / sigma)
    ...
    p(Y=K|F) = 1 - phi((a_{K-1} - F) / sigma)

    where phi is the cumulative density function of a Gaussian (the inverse probit
    function) and sigma is a parameter to be learned. A reference is:

    @article{chu2005gaussian,
      title={Gaussian processes for ordinal regression},
      author={Chu, Wei and Ghahramani, Zoubin},
      journal={Journal of Machine Learning Research},
      volume={6},
      number={Jul},
      pages={1019--1041},
      year={2005}
    }
"""
mutable struct Ordinal{T}<:AbstractLikelihood
    bin_edges::T
    o::Union{PyObject, Nothing}
end

function Ordinal(bin_edges)
    out = Ordinal(bin_edges, nothing)
    instantiate!(out)
    return out
end

function instantiate!(o::Union{Ordinal, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.likelihoods.Ordinal(o.bin_edges)
end

"""
    Poisson likelihood for use with count data, where the rate is given by the (transformed) GP.

    let g(.) be the inverse-link function, then this likelihood represents

    p(y_i | f_i) = Poisson(y_i | g(f_i) * binsize)

    Note:binsize
    For use in a Log Gaussian Cox process (doubly stochastic model) where the
    rate function of an inhomogeneous Poisson process is given by a GP.  The
    intractable likelihood can be approximated by gridding the space (into bins
    of size 'binsize') and using this Poisson likelihood.
"""
mutable struct Poisson{T}<:AbstractLikelihood
    # invlink # TODO: add capability to add custom invlink
    binsize::T
    o::Union{PyObject, Nothing}
end

function Poisson(;binsize::Real=1.0)
    out = Poisson(binsize, nothing)
    instantiate!(out)
    return out
end

function instantiate!(o::Union{Poisson, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.likelihoods.Poisson(;binsize=o.binsize)
end

"""
    This class represent a multi-class inverse-link function. Given a vector
    f=[f_1, f_2, ... f_k], the result of the mapping is

    y = [y_1 ... y_k]

    with

    y_i = (1-eps)  i == argmax(f)
          eps/(k-1)  otherwise.
"""
mutable struct RobustMax{T}<:AbstractLikelihood where T <: Real
    num_classes::Integer
    epsilon::T
    o::Union{PyObject, Nothing}
end

function RobustMax(num_classes::Integer; epsilon::Real=1e-3)
    out = RobustMax(num_classes, epsilon, nothing)
    instantiate!(out)
    return out
end

function instantiate!(o::Union{RobustMax, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.likelihoods.RobustMax(o.num_classes; epsilon=o.epsilon)
end

"""
    The soft-max multi-class likelihood.
"""
mutable struct SoftMax<:AbstractLikelihood
    num_classes::Integer
    o::Union{PyObject, Nothing}
end

function SoftMax(num_classes::Integer)
    out = SoftMax(num_classes, nothing)
    instantiate!(out)
    return out
end

function instantiate!(o::Union{SoftMax, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.likelihoods.SoftMax(o.num_classes)
end

mutable struct StudentT{T1,T2}<:AbstractLikelihood where {T1 <: Real, T2 <: Real}
    scale::T1
    df::T2
    o::Union{PyObject, Nothing}
end

function StudentT(;scale::Real=1.0, df::Real=3.0)
    out = StudentT(scale, df, nothing)
    instantiate!(out)
    return out
end

function instantiate!(o::Union{StudentT, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.likelihoods.StudentT(;scale=o.scale, df=o.df)
end

mutable struct SwitchedLikelihood{T}<:AbstractLikelihood
    likelihood_list::T
    o::Union{PyObject, Nothing}
end

function SwitchedLikelihood(likelihood_list)
    out = SwitchedLikelihood(likelihood_list, nothing)
    instantiate!(out)
    return out
end

function instantiate!(o::Union{SwitchedLikelihood, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.likelihoods.SwitchedLikelihood(likelihood_list)
end

end # module


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=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])) :
                     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]))
    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) in (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) in (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."
    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."
    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,outer=dot(k,w),(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
        nneighbor>0 && mapreduce(isequal(Inf),|,result) && @warn "minimumlengths warning: larger(>$coordination) coordination or smaller(<$nneighbor) nneighbor may be needed."
    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
                    isapprox(len,dist,atol=atol) && (push!(result,(nb,i,(j-1)%size(cluster,2)+1,disps[:,j]));break)
                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
                isapprox(len,dist,atol=atol) && (push!(result,(nb,i,j));break)
            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
    isequal(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)
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
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))
Point(pid::PID,rcoord::AbstractVector{<:Real},icoord::AbstractVector{<:Real}=zero(SVector{length(rcoord),Float}))=Point{length(rcoord)}(pid,rcoord,icoord)
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 join(string.(p.rcoord),",") join(string.(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],",")
        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
    isequal(lattice1::AbstractLattice,lattice2::AbstractLattice) -> Bool

Overloaded equivalent operator.
"""
Base.:(==)(lattice1::AbstractLattice,lattice2::AbstractLattice) = ==(efficientoperations,lattice1,lattice2)
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
    isequal(bonds1::Bonds,bonds2::Bonds) -> Bool

Overloaded equivalent operator.
"""
Base.:(==)(bonds1::Bonds,bonds2::Bonds) = ==(efficientoperations,bonds1,bonds2)
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:
"""
   compute the gradient of objective function with respect to velocity model for one
frequency slice, r is the residue obtained by subtracting observed data from synthetic,
u is source side wavefield with PML boundary part, H is the LU factorization
of Helmhotz operator.
"""
function velocity_gradient(r::MonochromaticRecordings, u::Vector{Tv},
         H, params::FdParams) where {Tv<:Complex{Float64}}

    # insert the conplex conjugate of residue to snapshot
    b = Vector{Complex{Float64}}(undef, params.Nz*params.Nx)
    inject_rec2spt!(b, r)

    # square of radian frequency
    omega2 = r.omega * r.omega

    # some coefficients of mixed-grid
    w1 = 0.6287326
    w2 = 0.0928166
    w3 = (1.0 - w1 - 4.0*w2) / 4.0

    # create a table save the coefficients for mass acceleration average
    C  = OffsetArray{Float64}(undef, -1:1, -1:1)

    C[-1,-1] = w3; C[-1,0] = w2; C[-1,1] = w3;
    C[ 0,-1] = w2; C[ 0,0] = w1; C[ 0,1] = w2;
    C[ 1,-1] = w3; C[ 1,0] = w2; C[ 1,1] = w3;

    # compute adjoint wavefield
    u_adj = zeros(Complex{Float64}, params.Nz*params.Nx)
    ldiv!(u_adj, transpose(H), b)

    # apply the imaging condition to get the gradient
    g = zeros(Complex{Float64}, params.nz * params.nx)

    # both u and u_adj are vectors include PML layers
    for i2 = 1 : params.nx
        idx2 = (i2-1) * params.nz

        for i1 = 1 : params.nz
            idx1 = idx2 + i1

            # the constant part
            tmp   = params.spt2wfd[idx1]
            alpha = -2.0 * omega2 * u[tmp]
            alpha = alpha / (params.rho[i1,i2] * (params.vel[i1,i2])^3)

            for d2 = -1 : 1
                l2 = (i2+params.npml+d2-1) * params.Nz

                for d1 = -1 : 1
                    if 1 <= i1+d1 || !params.free_surface
                       l1 = l2 + i1 + d1 + params.ntop

                       g[idx1] = g[idx1] + C[d1,d2] * u_adj[l1]
                    end
                end
            end

            g[idx1] = g[idx1] * alpha
        end
    end

    return real(g)
end

# some other test the transpose of LU decomposition
# A = sprand(Complex{Float64}, 10, 10, 0.3) + spdiagm(0=>rand(Complex{Float64}, 10));
# B = Matrix(transpose(A))
# B = Sparse(B
# H = lu(B)
# P = lu(A)
#
# x0 = rand(Complex{Float64}, 10);
# b  = Vector{Complex{Float64}}(undef, 10);
# mul!(b, transpose(A), x0);
#
# x = Vector{Complex{Float64}}(undef, 10);
# ldiv!(x, H, b)
#
# x1 = Vector{Complex{Float64}}(undef, 10);
# ldiv!(x1, transpose(P), b)