text:
struct BerryRGrid{T, LT, MT, MT1}
    hami::TbHami{T, LT, MT}
    A ::Vector{Vec3{MT}} #A_a(R) = <0|r_a|R> is the Fourier transform of the Berrry connection A_a(k) = i<u|del_a u> (a=x,y,z)the berry connection
    B ::Vector{Vec3{MT}} #B_a(R)=<0n|H(r-R)|Rm> is the Fourier transform of B_a(k) = i<u|H|del_a u> (a=x,y,z)
    C ::Vector{MT1} #CC_ab(R) = <0|r_a.H.(r-R)_b|R> is the Fourier transform of CC_ab(k) = <del_a u|H|del_b u> (a,b=x,y,z)}
end

#A_a(R) = <0|r_a|R> is the Fourier transform of the Berrry connection A_a(k) = i<u|del_a u> (a=x,y,z)the berry connection
#B_a(R)=<0n|H(r-R)|Rm> is the Fourier transform of B_a(k) = i<u|H|del_a u> (a=x,y,z)
function BerryRGrid(ab_initio_grid::AbInitioKGrid{T}, hami::TbHami, chk) where {T}
    irvec = map(x->x.R_cryst, hami)

    n_wann = n_wannier_functions(ab_initio_grid)

    berry_vec = () -> Vec3([zeros(Complex{T}, n_wann, n_wann) for i = 1:3]...)
    berry_mat = () -> Mat3([zeros(Complex{T}, n_wann, n_wann) for i = 1:9]...)

    n_kpoints = length(ab_initio_grid)

    A_q = [berry_vec() for k=1:n_kpoints]
    B_q = [berry_vec() for k=1:n_kpoints]
    C_q = [berry_mat() for k=1:n_kpoints]

    n_nearest = n_nearest_neighbors(ab_initio_grid)
    neighbor_weights = ab_initio_grid.neighbor_weights
    Threads.@threads for i in 1:n_kpoints
        kpoint = ab_initio_grid.kpoints[i]
        for n in 1:n_nearest
            neighbor_bond = kpoint.neighbors[n]
            weight        = neighbor_weights[n]
            vr            = ustrip.(neighbor_bond.vr)
            overlap = kpoint.overlaps[n]
            h    = kpoint.hamis[n]
            for v in 1:3
                t_fac = 1im * vr[v] * weight
                A_q[i][v] .+= t_fac .* overlap
                B_q[i][v] .+= t_fac .* h
                for n2 in 1:n_nearest
                    weight2 = neighbor_weights[n2]
                    neighbor_bond2 = kpoint.neighbors[n2]
                    vr2 = ustrip.(neighbor_bond2.vr)
                    uHu = kpoint.uHu[n2, n]
                    for v2 in 1:3
                        t_fac2 = t_fac * -1im * weight2 * vr2[v2]
                        C_q[i][v, v2] .+= t_fac2 .* uHu
                    end
                end
            end
        end
        for v in 1:3
            A_q[i][v] .= (A_q[i][v] + A_q[i][v]')/2
            for v2 in 1:v
                C_q[i][v, v2] .= C_q[i][v2, v]'
            end
        end
    end
    A_R = [berry_vec() for k=1:length(irvec)]
    B_R = [berry_vec() for k=1:length(irvec)]
    C_R = [berry_mat() for k=1:length(irvec)]
    fourier_q_to_R(ab_initio_grid.kpoints, irvec) do iR, ik, phase
        for v in 1:3
            A_R[iR][v] .+= phase .* A_q[ik][v]
            B_R[iR][v] .+= phase .* B_q[ik][v]
            for v2 in 1:3
                C_R[iR][v2, v] .+= phase .* C_q[ik][v2, v]
            end
        end
    end
    ws_shifts, ws_nshifts = chk.ws_shifts_cryst, chk.ws_nshifts
    points, degens = chk.ws_R_cryst, chk.ws_degens
    A_R_out = [berry_vec() for k=1:length(irvec)]
    B_R_out = [berry_vec() for k=1:length(irvec)]
    C_R_out = [berry_mat() for k=1:length(irvec)]
    for (ir1, (h, shifts, nshifts, degen)) in enumerate(zip(hami, ws_shifts, ws_nshifts, degens))
        for i in eachindex(h.block)
            ns = nshifts[i]
            frac = 1/(ns * degen)
             for is in 1:ns
                rcryst = h.R_cryst + shifts[i][is]
                ir2 = findfirst(x->x.R_cryst == rcryst, hami)
                for v in 1:3
                    A_R_out[ir2][v][i] += A_R[ir1][v][i]*frac
                    B_R_out[ir2][v][i] += B_R[ir1][v][i]*frac
                    for v2 in 1:3
                        C_R_out[ir2][v, v2][i] += C_R[ir1][v, v2][i]*frac
                    end
                end
            end
        end
    end

    for i in 1:length(A_R_out)
        for v=1:3
            A_R_out[i][v] ./= n_kpoints
            B_R_out[i][v] ./= n_kpoints
            for v2 in 1:3
                C_R_out[i][v, v2] ./= n_kpoints
            end
        end
    end
    return BerryRGrid(hami, A_R_out, B_R_out, C_R_out)
end

struct BerryKGrid{T, MT, MT1} <: AbstractKGrid{T}
    hamiltonian_kgrid::HamiltonianKGrid{T, MT}
    J_plus ::Vector{Vec3{MT}}
    J_minus::Vector{Vec3{MT}}
    J::Vector{Vec3{MT}}
    A::Vector{Vec3{MT}}
    Ω::Vector{Vec3{MT}} #pseudo form of berry connection A
    B::Vector{Vec3{MT}}
    C::Vector{MT1}
    f::Vector{MT}
    g::Vector{MT}
end
function BerryKGrid(berry_R_grid::BerryRGrid, kpoints::Vector{<:Vec3}, fermi::AbstractFloat)
    tb_hami = berry_R_grid.hami
    hamiltonian_kgrid = HamiltonianKGrid(tb_hami, kpoints)
    nk = length(kpoints)
    J_plus  = [Vec3(zeros_block(tb_hami), zeros_block(tb_hami), zeros_block(tb_hami)) for i=1:nk]
    J_minus = [Vec3(zeros_block(tb_hami), zeros_block(tb_hami), zeros_block(tb_hami)) for i=1:nk]
    J = [Vec3(zeros_block(tb_hami), zeros_block(tb_hami), zeros_block(tb_hami)) for i=1:nk]
    A_k = [Vec3(zeros_block(tb_hami), zeros_block(tb_hami), zeros_block(tb_hami)) for i=1:nk]
    Ω_k = [Vec3(zeros_block(tb_hami), zeros_block(tb_hami), zeros_block(tb_hami)) for i=1:nk]
    B_k = [Vec3(zeros_block(tb_hami), zeros_block(tb_hami), zeros_block(tb_hami)) for i=1:nk]
    C_k = [Mat3([zeros_block(tb_hami) for j = 1:9]...) for i=1:nk]
    f = [zeros_block(tb_hami) for i=1:nk]
    g = [zeros_block(tb_hami) for i=1:nk]
    Threads.@threads for i=1:nk
        Uk = hamiltonian_kgrid.eigvecs[i]
        Ek = hamiltonian_kgrid.eigvals[i]
        ∇Hk = Vec3(zeros_block(tb_hami), zeros_block(tb_hami), zeros_block(tb_hami))
        A = A_k[i]
        B = B_k[i]
        Ω = Ω_k[i]
        C = C_k[i]
        fourier_transform(tb_hami, kpoints[i]) do n, iR, R_cart, b, fac
            Rcart = ustrip.(R_cart)
            for v=1:3
                ∇Hk[v][n] += Rcart[v] * 1im * fac * b.block[n]
                A[v][n]   += fac * berry_R_grid.A[iR][v][n]

                B[v][n]   += fac * berry_R_grid.B[iR][v][n]
                for v2 = 1:3
                    C[v, v2][n] += fac * berry_R_grid.C[iR][v, v2][n]
                end
            end
            Ω[1][n] += 1im * fac * (Rcart[2] * berry_R_grid.A[iR][3][n] - Rcart[3] * berry_R_grid.A[iR][2][n])
            Ω[2][n] += 1im * fac * (Rcart[3] * berry_R_grid.A[iR][1][n] - Rcart[1] * berry_R_grid.A[iR][3][n])
            Ω[3][n] += 1im * fac * (Rcart[1] * berry_R_grid.A[iR][2][n] - Rcart[2] * berry_R_grid.A[iR][1][n])
        end
        J_plus_k  = J_plus[i]
        J_minus_k = J_minus[i]
        J_k = J[i]
        occupations_H_gauge = map(x -> x < fermi ? 1 : 0, Ek) #acting like it's only an insulator for now
        n_wann = length(Ek)

        f[i] .= Uk * diagm(0 => occupations_H_gauge) * Uk'

        g[i] .= map(x -> -x, f[i])
        for j = 1:n_wann
            g[i][j, j] += 1
        end
        for v = 1:3
            Hbar = Uk' * ∇Hk[v] * Uk
            for m = 1:n_wann
                Ek_m = Ek[m]
                for n = 1:n_wann
                    Ek_n = Ek[n]
                    if Ek_n > fermi && Ek_m < fermi
                        J_plus_k[v][n, m]  =  1im * Hbar[n, m]/(Ek[m] - Ek[n])
                        J_minus_k[v][m, n] =  1im * Hbar[m, n]/(Ek[n] - Ek[m])
                        #else is already taken care of by initializing with zeros
                    end
                    J_k[v][m, n] =  n==m ? 0.0 : 1im * Hbar[m, n]/(Ek[n] - Ek[m])
                end
            end
            J_plus_k[v]  .= Uk * J_plus_k[v] * Uk'
            J_minus_k[v] .= Uk * J_minus_k[v] * Uk'
            J_k[v] .= Uk * J_k[v] * Uk'
        end
    end
    return BerryKGrid(hamiltonian_kgrid, J_plus, J_minus, J, A_k, Ω_k, B_k, C_k, f, g)
end

n_wannier_functions(bgrid::BerryKGrid) = size(bgrid.f[1], 1)
n_kpoints(bgrid::BerryKGrid)           = length(bgrid.f)

core_kgrid(bgrid::BerryKGrid) = core_kgrid(bgrid.hamiltonian_kgrid)

for f in (:Hk, :eigvecs, :eigvals)
    @eval $f(bgrid::BerryKGrid) = $f(bgrid.hamiltonian_kgrid)
end

"""
    fourier_q_to_R(f::Function, q_vectors, R_vectors)

Performs a fourier transform from the ab-initio kpoints to the wigner seitz unit cells.
The function will be executed inside the fourier transform loop, being called like
`f(iR, ik, phase)`
"""
function fourier_q_to_R(f::Function, q_vectors, R_vectors)
    for iR in 1:length(R_vectors)
        for ik in 1:length(q_vectors)
            phase = exp(-2im * π * (k_cryst(q_vectors[ik]) ⋅ R_vectors[iR]))
            f(iR, ik, phase)
        end
    end
end

const pseudo_α = Vec3(2,3,1)
const pseudo_β = Vec3(3,1,2)

function orbital_angular_momentum(berry_K_grid::BerryKGrid{T}, fermi) where {T}
    nwann = n_wannier_functions(berry_K_grid)
    nk = n_kpoints(berry_K_grid)

    M_local     = [Vec3([zeros(T, nwann, nwann) for i =1:3]...) for ik=1:nk]
    M_itinerant = [Vec3([zeros(T, nwann, nwann) for i =1:3]...) for ik=1:nk]

    Threads.@threads for ik in 1:nk
        f = berry_K_grid.f[ik]
        g = berry_K_grid.g[ik]
        A = berry_K_grid.A[ik]
        B = berry_K_grid.B[ik]
        H = berry_K_grid.hamiltonian_kgrid.Hk[ik]
        C = berry_K_grid.C[ik]
        J = berry_K_grid.J[ik]
        Ω = berry_K_grid.Ω[ik]
        for iv in 1:3
            α = pseudo_α[iv]
            β = pseudo_β[iv]

            Fαβ = real.(f * Ω[iv]) .- 2 .* imag.(f * A[α] * g * J[β] .+ f * J[α] * g * A[β] .+ f * J[α] * g * J[β])

            Hαβ = real.(f * H * f * Ω[iv]) .+ 2 .* imag.(f * H * f * A[α] * f * A[β] .- f * H * f * (A[α] * g * J[β] + J[α] * g * A[β] + J[α] * g * J[β]))

            J0 = real.(f * -1im*(C[α, β] - C[α, β]')) .- 2 .* imag.(f * H * A[α] * f * A[β])
            J1 = - 2 .* imag.(f * J[α] * g * B[β] .- f * J[β] * g * B[α])
            J2 = - 2 .* imag.(f * J[α] * g * H * g * J[β])

            Gαβ = J0 + J1 + J2

            M_local[ik][iv] .= Gαβ .- fermi .* Fαβ
            M_itinerant[ik][iv] .= Hαβ .- fermi .* Fαβ
        end
    end
    return M_local, M_itinerant
end

function orbital_angular_momentum_w90(berry_K_grid::BerryKGrid{T}) where {T}
    nwann = n_wannier_functions(berry_K_grid)
    nk = n_kpoints(berry_K_grid)

    non_traced_Fαβ = [Vec3([zeros(T, nwann, nwann) for i =1:3]...) for ik=1:nk]
    non_traced_Hαβ = [Vec3([zeros(T, nwann, nwann) for i =1:3]...) for ik=1:nk]
    non_traced_Gαβ = [Vec3([zeros(T, nwann, nwann) for i =1:3]...) for ik=1:nk]

    Threads.@threads for ik in 1:nk
        f = berry_K_grid.f[ik]
        g = berry_K_grid.g[ik]
        A = berry_K_grid.A[ik]
        B = berry_K_grid.B[ik]
        H = berry_K_grid.hamiltonian_kgrid.Hk[ik]
        C = berry_K_grid.C[ik]
        J_plus = berry_K_grid.J_plus[ik]
        J_minus = berry_K_grid.J_minus[ik]
        Ω = berry_K_grid.Ω[ik]
        for iv in 1:3
            α = pseudo_α[iv]
            β = pseudo_β[iv]

            non_traced_Fαβ[ik][iv] .+= real.(f * Ω[iv]) .- 2 .* imag.(A[α] * J_plus[β] .+ J_minus[α] * A[β] .+ J_minus[α] * J_plus[β])
            t_1 = H * A[α]
            t_3 = H * Ω[iv]

            t_2 = -1im*(C[α, β] - C[α, β]')

            t_4 = t_1 * f

            t_5 = t_4 * A[β]

            s = 2 .* imag.(f * t_5)

            non_traced_Gαβ[ik][iv] .+= real.(f * t_2) .- s
            non_traced_Hαβ[ik][iv] .+= real.(f * t_3) .+ s

            t_4 = H * J_minus[α]

            non_traced_Gαβ[ik][iv] .-= 2 .* imag.(J_minus[α] * B[β] .- J_minus[β] * B[α])

            non_traced_Hαβ[ik][iv] .-= 2 .* imag.(t_1 * J_plus[β] .+ t_4 * A[β])
            t_4 = J_minus[α] * H
            t_5 = H * J_minus[α]

            non_traced_Gαβ[ik][iv] .-= 2 .* imag.(t_4 * J_plus[β])

            non_traced_Hαβ[ik][iv] .-= 2 .* imag.(t_5 * J_plus[β])

        end
    end
    return non_traced_Fαβ./nk, non_traced_Hαβ./nk, non_traced_Gαβ./nk
end


text:
module AdventOfCode2021
using Printf

const CURRENT_DAY = 20

include("dispatch.jl")
for day in 1:CURRENT_DAY
    include(@sprintf("day%02d.jl", day))
end

export Day, CURRENT_DAY
export data, testdata
export test1, test2, testresult1, testresult2
export solve1, solve2

end # module

ext: jl


text:
ismacro(ex, name::Symbol, nargs::Integer=-1) =
    isa(ex, Expr) && ex.head==:macrocall && ex.args[1]==name &&
    (nargs == -1 || length(ex.args) == nargs+1)
ismacro(ex, name::String, nargs::Integer=-1) = ismacro(ex, Symbol(name), nargs)

iscall(ex, name::Symbol, nargs::Integer=-1) =
    isa(ex, Expr) && ex.head==:call && ex.args[1]==name &&
    (nargs == -1 || length(ex.args) == nargs+1)


text:
using Documenter, Ferrite, Pkg

Pkg.precompile()

# Generate examples
include("generate.jl")

GENERATEDEXAMPLES = [joinpath("examples", f) for f in (
    "heat_equation.md",
    "l2_projection.md",
    "helmholtz.md",
    "incompressible_elasticity.md",
    "hyperelasticity.md",
    "threaded_assembly.md",
    "plasticity.md",
    "transient_heat_equation.md",
    "landau.md"
    )]

# Build documentation.
makedocs(
    format = Documenter.HTML(prettyurls = haskey(ENV, "GITHUB_ACTIONS")), # disable for local builds
    sitename = "Ferrite.jl",
    doctest = false,
    # strict = VERSION.minor == 6 && sizeof(Int) == 8, # only strict mode on 0.6 and Int64
    strict = false,
    pages = Any[
        "Home" => "index.md",
        "manual/fe_intro.md",
        "Manual" => [
            "manual/degrees_of_freedom.md",
            "manual/assembly.md",
            "manual/boundary_conditions.md",
            "manual/grid.md",
            "manual/export.md"
            ],
        "Examples" => GENERATEDEXAMPLES,
        "API Reference" => [
            "reference/quadrature.md",
            "reference/interpolations.md",
            "reference/fevalues.md",
            "reference/dofhandler.md",
            "reference/assembly.md",
            "reference/boundary_conditions.md",
            "reference/grid.md",
            "reference/export.md"
            ]
        ],
)

# make sure there are no *.vtu files left around from the build
cd(joinpath(@__DIR__, "build", "examples")) do
    foreach(file -> endswith(file, ".vtu") && rm(file), readdir())
end

# Deploy built documentation from Travis.
deploydocs(
    repo = "github.com/Ferrite-FEM/Ferrite.jl.git",
    push_preview=true,
)


text:
include("../required.jl")

transect

all_ids = vcat(unique(transect.newPredID), unique(transect.newPreyID))

##Create Adjacency matrix
#create n x n matrix of zeros
M = zeros(Int16, (length(unique(transect.newPredID))+length(unique(transect.newPreyID)))+1, (length(unique(transect.newPredID))+length(unique(transect.newPreyID)))+1)
#assign names for cols
for y in 2:length(all_ids)+1
    M[y,1] = all_ids[y-1]
end

#assign names for rows
for x in 2:length(all_ids)+1
    M[1,x] = all_ids[x-1]
end
M

# Fill adjacency matrix with ones where there are interactions
for i in 1:length(all_ids)
    M[transect[i, 1]+1, transect[i, 2]+1] = 1
end
M

# # Test to make sure the adjacency matrix is correctly filled
# transectTest = [[1,6],[2,7], [3,8], [4,9],[5, 10]]
# MTest = zeros(Int16, 6, 6)
# for y in 2:6
#     MTest[y,1] = transectTest[y-1][1]
# end
#
# for x in 2:6
#     MTest[1,x] = transectTest[x-1][2]
# end
#
# for i in 1:5
#     MTest[i + 1, transectTest[i][2] - 5 + 1] = 1
# end

max_stars_repo_name: kiristern/reseaux

max_issues_repo_name: kiristern/reseaux

max_forks_repo_name: kiristern/reseaux


text:
using Gnuplot
x = y = -5:1:5
z = [x^2 + y^2 for x in x, y in y]
@gsp x y z "w pm3d notit" "set view map" "set auto fix" "set size square"
@gsp :- xlab = "x" ylab = "y" "set cblabel 'z'" palette(:plasma)
@gsp :- title = "3D Heatmap"
# please also see Heatmap discrete and make sure you are using the right one.
save(term ="pngcairo size 600,400", output ="heatmap1.png")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
function maketinfo(cinfo,n)
tinfo = Tinfo()

@unpack en,P,M,Ω,basis = cinfo
x,wx,Tx,y,wy,Ty = makealltrans(M,n,Ω,basis)
W = wx.*wy'
@pack tinfo = Tx,Ty,W

return tinfo
end


text:
# run all tests

using NBInclude

include("d2f_zeros_test.jl")
include("overapprox_nd_unittest.jl")

@testset "test methods notebook" begin
    @nbinclude(joinpath(dirname(@__DIR__), "plots", "Methods Section.ipynb"))
end


text:
function ten_to_two(n)
    b = zeros(Int8,10)
    temp_b = zeros(Int8,10)

    ite=1
    while n > 0
        b[ite] = n % 2
        n = Int((n-b[ite])/2)
        ite = ite +1
    end

    #=
    for i in 1:6
        temp_b[i] = b[6+1-i]
    end
    =#

    return b
end

function two_to_ten(b)
    n = Int8(0)

    for i in 1:length(b)
        n = n + 2^(i-1)*b[i]
    end
    return n
end

# ----------------------
# -- Inverse matrix   --
# ----------------------
function inverse_matrix(eu,ev)
# ----------------------
# -- A = ( eu[1]  ev[1])  --
# --     ( eu[2]  ev[2])  --
# ----------------------
    invA = zeros(2,2)
    detA = eu[1]*ev[2] - eu[2]*ev[1]

    invA[1,1] = (ev[2])/detA
    invA[1,2] = (-ev[1])/detA
    invA[2,1] = (-eu[2])/detA
    invA[2,2] = (eu[1])/detA
    return invA
end

function cal_morton(ulx,uly,lrx,lry,n_div) # 0<=s,tでreturn
    # 左上と右下の座標を整数で
    morton_belong = 0
    morton_num = 0

    upl_morton = zeros(Int8,10)
    lowr_morton = zeros(Int8,10)
    temp_morton = zeros(Int8,10)

    # 左上の点のモートン空間番号
    x_two = ten_to_two(ulx)
    y_two = ten_to_two(uly)

    for k in 1:5
        upl_morton[k*2-1] = x_two[k]
        upl_morton[k*2] = y_two[k]
    end

    # 右下の点のモートン空間番号
    x_two = ten_to_two(lrx)
    y_two = ten_to_two(lry)

    for k in 1:5
        lowr_morton[k*2-1] = x_two[k]
        lowr_morton[k*2] = y_two[k]
    end

    # 排他的論理和 xor
    for i in 1:10
        if upl_morton[i]+lowr_morton[i] == 2
            temp_morton[i] = 0
        else
            temp_morton[i] = upl_morton[i]+lowr_morton[i]
        end
    end

    # 所属空間のチェック

    ite = Int64(0)
    for i in 1:n_div
        s_even = 2*(n_div-(i-1))
        s_odd = 2*(n_div-(i-1)) - 1

        if temp_morton[s_even] == 1 || temp_morton[s_odd] == 1
            morton_belong = ite                             # ルート空間
            break
        end
        #=
        if temp_morton[s_even] == 1 || temp_morton[s_odd] == 1
            morton_belong = 0                             # ルート空間
        elseif temp_morton[4] == 1 || temp_morton[3] == 1
            morton_belong = 1                             # 親空間
        elseif temp_morton[2] == 1 || temp_morton[1] == 1
            morton_belong = 2                             # 子空間
        else
            morton_belong = 3                             # 孫空間
        end
        =#
        ite += 1
    end

    # 所属空間番号
    temp_morton = zeros(Int8,10)
    if morton_belong == 0                             # ルート空間
        morton_num = 0
    elseif morton_belong == 1                         # 親空間
        s_even = 2*(n_div)
        s_odd = 2*(n_div) - 1

        # 4右シフト
        temp_morton[1] = lowr_morton[s_odd]
        temp_morton[2] = lowr_morton[s_even]

        morton_num = two_to_ten(temp_morton)

    elseif morton_belong == 2                         # 子空間

        s_even = 2*(n_div)
        s_odd = 2*(n_div) - 1

        # 2右シフト
        temp_morton[1] = lowr_morton[s_odd-2]
        temp_morton[2] = lowr_morton[s_even-2]
        temp_morton[3] = lowr_morton[s_odd]
        temp_morton[4] = lowr_morton[s_even]

        morton_num = two_to_ten(temp_morton)
    elseif morton_belong == 3                         # 孫空間
        s_even = 2*(n_div)
        s_odd = 2*(n_div) - 1

        # 0右シフト
        temp_morton[1] = lowr_morton[s_odd-4]
        temp_morton[2] = lowr_morton[s_even-4]
        temp_morton[3] = lowr_morton[s_odd-2]
        temp_morton[4] = lowr_morton[s_even-2]
        temp_morton[5] = lowr_morton[s_odd]
        temp_morton[6] = lowr_morton[s_even]

        morton_num = two_to_ten(temp_morton)

    elseif morton_belong == 4                         # 孫空間
        temp_morton = copy(lowr_morton)
        morton_num = two_to_ten(temp_morton)
    end

    return morton_belong,morton_num
end

function liner_morton(n_div)
    # 両方の探査を行うため，煩雑
    s = zeros(Int64,5)

    for j in 1:5
        for i in 1:(n_div+1)
            s[j] += floor(4.0^(i-j))
        end
        s[j] += (j-1)
    end

    search0 = zeros(Int64,4^0,s[1])
    search1 = zeros(Int64,4^1,s[2])
    search2 = zeros(Int64,4^2,s[3])
    search3 = zeros(Int64,4^3,s[4])
    search4 = zeros(Int64,4^4,s[5])

    # ルート空間：s0
    for i in 1:1
        for j in 1:s[1]
            search0[i,j] = j
        end
    end

    # 親空間：s1
    for i in 1:4
        search1[i,1] = 1　　　　　　　　　　　　　　　　# ルート空間
        search1[i,2] = i+1                           # 親空間
        for j in 1:4
            search1[i,j+2] = 6 + 4*(i-1) + (j-1)     # 子空間
        end
        for j in 1:16
            search1[i,j+6] = 22 + 16*(i-1) + (j-1)   # 孫空間
        end

        if n_div ==4
            for j in 1:64
                search1[i,j+22] = 86 + 64*(i-1) + (j-1)   # ひ孫空間
            end
        end
    end

    # 子空間：s2
    for i in 1:16
        search2[i,1] = 1                             # ルート空間
        search2[i,2] = div(i-1,4) + 2                # 親空間
        search2[i,3] = i+5                           # 子空間
        for j in 1:4
            search2[i,j+3] = 22 + 4*(i-1) + (j-1)    # 孫空間
        end
        if n_div ==4
            for j in 1:16
                search2[i,j+7] = 86 + 16*(i-1) + (j-1)    # 孫空間
            end
        end
    end

    # 孫空間
    for i in 1:64
        search3[i,1] = 1                             # ルート空間
        search3[i,2] = div(i-1,16) + 2               # 親空間
        search3[i,3] = div(i-1,4) + 6                # 子空間
        search3[i,4] = 22 + (i-1)                    # 孫空間
        if n_div >=4
            for j in 1:4
                search3[i,j+4] = 86 + 4*(i-1) + (j-1)    # 孫空間
            end
        end
    end

    # ひ孫空間
    if n_div >=4
        for i in 1:256
            search4[i,1] = 1                             # ルート空間
            search4[i,2] = div(i-1,64) + 2               # 親空間
            search4[i,3] = div(i-1,16) + 6                # 子空間
            search4[i,4] = div(i-1,4) + 22                # 孫空間
            search4[i,5] = 86 + (i-1)                    # 孫空間
        end
    end

    return s,search0,search1,search2,search3,search4
end

max_stars_repo_name: hide-dog/Qtree

max_issues_repo_name: hide-dog/Qtree

max_forks_repo_name: hide-dog/Qtree


text:
"""
    RasterDataSource

Abstract supertype for raster data collections.
"""
abstract type RasterDataSource end

"""
    RasterDataSet

Abstract supertype for datasets that belong to a [`RasterDataSource`](@ref).
"""
abstract type RasterDataSet end

"""
    BioClim <: RasterDataSet

BioClim datasets. Usually containing layers from `1:19`.
These can also be accessed with `:bioX`, e.g. `:bio5`.

They do not usually use `month` or `date` keywords, but may use
`date` in past/future scenarios.

Currently implemented for WorldClim and CHELSA as `WorldClim{BioClim}`,
`CHELSA{BioClim}` and `CHELSA{Future{BioClim, args..}}`.

See the [`getraster`](@ref) docs for implementation details.
"""
struct BioClim <: RasterDataSet end

# Bioclim has standardised layers for all data sources
layers(::Type{BioClim}) = values(bioclim_lookup)
layerkeys(T::Type{BioClim}) = keys(bioclim_lookup)
layerkeys(T::Type{BioClim}, layer) = bioclim_key(layer)
layerkeys(T::Type{BioClim}, layers::Tuple) = map(l -> bioclim_key(l), layers)

const bioclim_lookup = (
    bio1 = 1,
    bio2 = 2,
    bio3 = 3,
    bio4 = 4,
    bio5 = 5,
    bio6 = 6,
    bio7 = 7,
    bio8 = 8,
    bio9 = 9,
    bio10 = 10,
    bio11 = 11,
    bio12 = 12,
    bio13 = 13,
    bio14 = 14,
    bio15 = 15,
    bio16 = 16,
    bio17 = 17,
    bio18 = 18,
    bio19 = 19,
)

# We allow a range of bioclim keys, as they are listed with
# a lot of variants on CHELSA and WorldClim
bioclim_key(k::Symbol) = bioclim_key(string(k))
bioclim_key(k::AbstractString) = Symbol(replace(lowercase(k), "_" => ""))
bioclim_key(k::Integer) = keys(bioclim_lookup)[k]

bioclim_int(k::Integer) = k
bioclim_int(k::Symbol) = bioclim_lookup[bioclim_key(k)]

"""
    Climate <: RasterDataSet

Climate datasets. These are usually months of the year, not specific dates,
and use a `month` keyword in `getraster`. They also use `date` in past/future scenarios.

Currently implemented for WorldClim and CHELSA as `WorldClim{Climate}`,
`CHELSA{Climate}` and `CHELSA{Future{Climate, args..}}`.

See the [`getraster`](@ref) docs for implementation details.
"""
struct Climate <: RasterDataSet end

months(::Type{Climate}) = ntuple(identity, Val{12})

"""
    Weather <: RasterDataSet

Weather datasets. These are usually large time-series of specific dates,
and use a `date` keyword in `getraster`.

Currently implemented for WorldClim and CHELSA as `WorldClim{Weather}`,
and `CHELSA{Weather}`

See the [`getraster`](@ref) docs for implementation details.
"""
struct Weather <: RasterDataSet end

"""
    LandCover <: RasterDataSet

Land-cover datasets.

Currently implemented for EarthEnv as `EarchEnv{LandCover}`.

See the [`getraster`](@ref) docs for implementation details.
"""
struct LandCover{X} <: RasterDataSet end

"""
    HabitatHeterogeneity <: RasterDataSet

Habitat heterogeneity datasets.

Currently implemented for EarchEnv as `EarchEnv{HabitatHeterogeneity}`.

See the [`getraster`](@ref) docs for implementation details.
"""
struct HabitatHeterogeneity <: RasterDataSet end

"""
    ClimateModel

Abstract supertype for climate models use in [`Future`](@ref) datasets.
"""
abstract type ClimateModel end

struct ACCESS1 <: ClimateModel end
struct BNUESM <: ClimateModel end
struct CCSM4 <: ClimateModel end
struct CESM1BGC <: ClimateModel end
struct CESM1CAM5 <: ClimateModel end
struct CMCCCMS <: ClimateModel end
struct CMCCCM <: ClimateModel end
struct CNRMCM5 <: ClimateModel end
struct CSIROMk3 <: ClimateModel end
struct CanESM2 <: ClimateModel end
struct FGOALS <: ClimateModel end
struct FIOESM <: ClimateModel end
struct GFDLCM3 <: ClimateModel end
struct GFDLESM2G <: ClimateModel end
struct GFDLESM2M <: ClimateModel end
struct GISSE2HCC <: ClimateModel end
struct GISSE2H <: ClimateModel end
struct GISSE2RCC <: ClimateModel end
struct GISSE2R <: ClimateModel end
struct HadGEM2AO <: ClimateModel end
struct HadGEM2CC <: ClimateModel end
struct IPSLCM5ALR <: ClimateModel end
struct IPSLCM5AMR <: ClimateModel end
struct MIROCESMCHEM <: ClimateModel end
struct MIROCESM <: ClimateModel end
struct MIROC5 <: ClimateModel end
struct MPIESMLR <: ClimateModel end
struct MPIESMMR <: ClimateModel end
struct MRICGCM3 <: ClimateModel end
struct MRIESM1 <: ClimateModel end
struct NorESM1M <: ClimateModel end
struct BCCCSM1 <: ClimateModel end
struct Inmcm4 <: ClimateModel end
struct BCCCSM2MR <: ClimateModel end
struct CNRMCM61 <: ClimateModel end
struct CNRMESM21 <: ClimateModel end
struct CanESM5 <: ClimateModel end
struct GFDLESM4 <: ClimateModel end
struct IPSLCM6ALR <: ClimateModel end
struct MIROCES2L <: ClimateModel end
struct MIROC6 <: ClimateModel end
struct MRIESM2 <: ClimateModel end
struct UKESM <: ClimateModel end
struct MPIESMHR <: ClimateModel end

"""
    CMIPphase

Abstract supertype for phases of the CMIP,
the Coupled Model Intercomparison Project.

Subtypes are `CMIP5` and `CMIP6`.
"""
abstract type CMIPphase end

"""
    CMIP5 <: CMIPphase

The Coupled Model Intercomparison Project, Phase 5.
"""
struct CMIP5 <: CMIPphase end

"""
    CMIP6 <: CMIPphase

The Coupled Model Intercomparison Project, Phase 6.
"""
struct CMIP6 <: CMIPphase end

"""
    ClimateScenario

Abstract supertype for scenarios used in [`CMIPphase`](@ref) models.
"""
abstract type ClimateScenario end

"""
    RepresentativeConcentrationPathway

Abstract supertype for Representative Concentration Pathways (RCPs) for [`CMIP5`](@ref).

Subtypes are: `RCP26`, `RCP45`, `RCP60`, `RCP85`
"""
abstract type RepresentativeConcentrationPathway <: ClimateScenario end

struct RCP26 <: RepresentativeConcentrationPathway end
struct RCP45 <: RepresentativeConcentrationPathway end
struct RCP60 <: RepresentativeConcentrationPathway end
struct RCP85 <: RepresentativeConcentrationPathway end

"""
    SharedSocioeconomicPathway

Abstract supertype for Shared Socio-economic Pathways (SSPs) for [`CMIP6`](@ref).

Subtypes are: `SSP126`, `SSP245`, SSP370`, SSP585`
"""
abstract type SharedSocioeconomicPathway <: ClimateScenario end

struct SSP126 <: SharedSocioeconomicPathway end
struct SSP245 <: SharedSocioeconomicPathway end
struct SSP370 <: SharedSocioeconomicPathway end
struct SSP585 <: SharedSocioeconomicPathway end

"""
    Future{<:RasterDataSet,<:CMIPphase,<:ClimateModel,<:ClimateScenario}

Future climate datasets specified with a dataset, phase, model, and scenario.

## Type Parameters

#### `RasterDataSet`

Currently [`BioClim`](@ref) and [`Climate`](@ref) are implemented
for the [`CHELSA`](@ref) data source.

#### `CMIPphase`

Can be either [`CMIP5`](@ref) or [`CMIP6`](@ref).

#### `ClimateModel`

Climate models can be chosen from:

`ACCESS1`, `BNUESM`, `CCSM4`, `CESM1BGC`, `CESM1CAM5`, `CMCCCMS`, `CMCCCM`,
`CNRMCM5`, `CSIROMk3`, `CanESM2`, `FGOALS`, `FIOESM`, `GFDLCM3`, `GFDLESM2G`,
`GFDLESM2M`, `GISSE2HCC`, `GISSE2H`, `GISSE2RCC`, `GISSE2R`, `HadGEM2AO`,
`HadGEM2CC`, `IPSLCM5ALR`, `IPSLCM5AMR`, `MIROCESMCHEM`, `MIROCESM`, `MIROC5`,
`MPIESMLR`, `MPIESMMR`, `MRICGCM3`, `MRIESM1`, `NorESM1M`, `BCCCSM1`, `Inmcm4`,
`BCCCSM2MR`, `CNRMCM61`, `CNRMESM21`, `CanESM5`, `MIROCES2L`, `MIROC6` for CMIP5;

`UKESM`, `MPIESMHR` `IPSLCM6ALR` `MRIESM2`, `GFDLESM4` for `CMIP6`.

#### `ClimateScenario`

CMIP5 Climate scenarios are all [`RepresentativeConcentrationPathway`](@ref)
and can be chosen from: `RCP26`, `RCP45`, `RCP60`, `RCP85`

CMIP6 Climate scenarios are all [`SharedSocioeconomicPathway`](@ref) and
can be chosen from: `SSP126`, `SSP245`, SSP370`, SSP585`

## Example

```jldoctest future
using RasterDataSources
dataset = Future{BioClim, CMIP5, BNUESM, RCP45}
# output
Future{BioClim, CMIP5, BNUESM, RCP45}
```
Currently `Future` is only implented for `CHELSA`

```jldoctest future
datasource = CHELSA{Future{BioClim, CMIP5, BNUESM, RCP45}}
```

"""
struct Future{D<:RasterDataSet,C<:CMIPphase,M<:ClimateModel,S<:ClimateScenario} end

_dataset(::Type{<:Future{D}}) where D = D
_phase(::Type{<:Future{<:Any,P}}) where P = P
_model(::Type{<:Future{<:Any,<:Any,M}}) where M = M
_scenario(::Type{<:Future{<:Any,<:Any,<:Any,S}}) where S = S

max_forks_repo_name: cesaraustralia/RasterDataSources


text:
@testset "@variable macro with Poly" begin
    @polyvar x y
    X = [1, y, x^2, x, y^2]

    m = Model()
    @variable(m, a)
    var_poly   = polynomial(a * x * y)
    var_poly_x = polynomial(a * x)
    aff_poly   = polynomial((a + 1) * x * y)
    aff_poly_x = polynomial((a + 1) * x)
    @test_throws ErrorException @variable m p Poly(X) unknown_kw=1
    @test_throws ErrorException @variable m p == 1 Poly(X)
    @test_throws ErrorException @variable m p Poly(X) start=1
    @test_throws ErrorException @variable m 0 <= p <= 1 Poly(X)
    @test_throws ErrorException @variable m 0 >= p >= 1 Poly(X)
    @test_throws ErrorException @variable m p <= 0 Poly(X)
    @test_throws ErrorException @variable m p >= 1 Poly(X)
    @test_throws ErrorException @variable m p == 1 Poly(X)
    @test_throws ErrorException @variable m p >= 0 Poly(X)
    @test_throws ErrorException @variable(m, p3[2:3] >= 0, Poly(X))

    function testvar(m, p, x, binary = false, integer = false,
                     vars = true, use_y = true)
        if vars
            PT = use_y ? typeof(var_poly) : typeof(var_poly_x)
        else
            PT = use_y ? typeof(aff_poly) : typeof(aff_poly_x)
        end
        @test isa(p, PT)
        @test monomials(p) == monovec(x)
        if vars
            @test all(α -> JuMP.is_binary(α) == binary, coefficients(p))
            @test all(α -> JuMP.is_integer(α) == integer, coefficients(p))
        else
            @test all(α -> JuMP.is_binary(first(keys(α.terms))) == binary, coefficients(p))
            @test all(α -> JuMP.is_integer(first(keys(α.terms))) == integer, coefficients(p))
        end
    end

    @testset "MonomialBasis" begin
        m = Model()
        @variable m p1[1:3] Poly(X)
        @test isa(p1, Vector{typeof(var_poly)})
        testvar(m, p1[1], X)
        @variable(m, p2, Poly(X), integer=true)
        testvar(m, p2, X, false, true)
        @variable(m, p3[2:3], Poly(X))
        @test isa(p3, JuMP.Containers.DenseAxisArray{typeof(var_poly),1,Tuple{UnitRange{Int}}})
        testvar(m, p3[2], X)
        @variable(m, p4[i=2:3,j=i:4], Poly(X), binary=true)
        testvar(m, p4[2,3], X, true)
    end

    @testset "ScaledMonomialBasis" begin
        m = Model()
        @variable(m, p1, Poly(MB.ScaledMonomialBasis([1, x, x^2])), Int)
        testvar(m, p1, monovec([1, x, x^2]), false, true, false, false)
    end

    @testset "FixedPolynomialBasis" begin
        m = Model()
        @variable(m, p1, Poly(MB.FixedPolynomialBasis([1 - x^2, x^2 + 2])), Bin)
        testvar(m, p1, monovec([x^2, 1]), true, false, false, false)
        @variable(m, p2[1:2], Poly(MB.FixedPolynomialBasis([1 - x^2, x^2 + 2])))
        testvar(m, p2[1], monovec([x^2, 1]), false, false, false, false)
        # Elements of the basis have type monomial
        @variable(m, p3[2:3], Poly(MB.FixedPolynomialBasis([x, x^2])))
        testvar(m, p3[2], monovec([x^2, x]), false, false, false, false)
        # Elements of the basis have type term
        @variable(m, p4[1:2], Poly(MB.FixedPolynomialBasis([1, x, x^2])), Int)
        testvar(m, p4[1], monovec([x^2, x, 1]), false, true, false, false)
        # Elements of the basis have type variable
        @variable(m, p5[-1:1], Poly(MB.FixedPolynomialBasis([x, y])), integer=true)
        testvar(m, p5[0], monovec([x, y]), false, true, false)
        @variable(m, p6[-1:1], Poly(MB.FixedPolynomialBasis([x])), integer=true)
        testvar(m, p6[0], monovec([x]), false, true, true, false)
    end

    @testset "ChebyshevBasis" begin
        m = Model()
        X = MB.maxdegree_basis(MB.ChebyshevBasis, [x, y], 2)
        @variable m p1[1:3] Poly(X)
        @test isa(p1, Vector{typeof(aff_poly)})
        testvar(m, p1[1], monomials([x,y], 0:2), false, false, false)
        @variable(m, p2, Poly(X), integer=true)
        testvar(m, p2,  monomials([x,y], 0:2), false, true, false)
        @variable(m, p3[2:3], Poly(X))
        @test isa(p3, JuMP.Containers.DenseAxisArray{typeof(aff_poly),1,Tuple{UnitRange{Int}}})
        testvar(m, p3[2], monomials([x,y], 0:2), false, false, false)
        @variable(m, p4[i=2:3,j=i:4], Poly(X), binary=true)
        testvar(m, p4[2,3], monomials([x,y], 0:2), true, false, false)
    end

end

@testset "JuMP.value function" begin
    m = Model()
    @variable m α
    @variable m β
    @polyvar x y
    p = α * x*y + β * x^2
    JuMP.fix(α, 2)
    JuMP.fix(β, 3)
    @test_broken JuMP.value(p) == 2x*y + 3x^2
end


text:
function tol_eq(x,y,eps)
    if abs(x-y)<eps
        return true
    else
        return false
    end
end


text:
# TODO ids?
# TODO docs
# TODO comments
# TODO multi word tokens
# TODO empty tokens
# TODO skip malformed trees?
"""
    Treebank

A corpus of sentences annotated as dependency trees on disk.
"""
struct Treebank{F}
    files::Vector{String}
    read_token::F
end

"""
    Treebank(file::String, read_token)

Read treebank from file `file`, calling `read_token` on each line to read tokens.
"""
Treebank(file::String, read_token) = Treebank([file], read_token)

"""
    Treebank(file)

Read a treebank from `files`, detecting the format from the filenames.
"""
Treebank(file::String) = Treebank([file])

"""
    Treebank(files)

Read a treebank from the files `files`, attempting to detect the format.
"""
function Treebank(files::Vector{String})
    if all(file -> endswith(file, ".conllu") || endswith(file, ".conll"), files)
        Treebank(files, from_conllu)
    else
        error("don't know how to read $files")
    end
end

function Base.iterate(tb::Treebank)
    r = TreeReader(first(tb.files), tb.read_token)
    return iterate(tb, (r, 1))
end
function Base.iterate(tb::Treebank, state)
    r, i = state
    tree = read_tree(r)
    if tree === nothing
        if i < length(tb.files)
            j = i + 1
            return iterate(tb, (TreeReader(tb.files[j], tb.read_token), j))
        end
    else
        return tree, state
    end
end

function Base.show(io::IO, tb::Treebank)
    fs = length(tb.files)
    len_str = length(tb.files) == 1 ? "1 file" : "$fs files"
    print(io, "Treebank ($len_str)")
end

Base.IteratorSize(treebank::Treebank) = Base.SizeUnknown()

struct TreeReader
    io::IO
    read_token
end

TreeReader(file::String, read_token) = TreeReader(open(file), read_token)

function read_tree(trees::TreeReader)
    io = trees.io
    if eof(io)
        close(io)
        return nothing
    end
    str = readuntil(io, "\n\n")
    return deptree(str, trees.read_token)
end


text:
using CSV, DataFrames, Dates, Test, BusinessDays
using AbnormalReturns

##

df_firm = CSV.File(joinpath("data", "daily_ret.csv")) |> DataFrame
df_mkt = CSV.File(joinpath("data", "mkt_ret.csv")) |> DataFrame
df_res = CSV.File(joinpath("data", "car_results.csv")) |> DataFrame
df_mkt[!, :mkt] = df_mkt.mktrf .+ df_mkt.rf
##

data = MarketData(
    df_mkt,
    df_firm
)

##

rr = quick_reg(data[18428, Date(2019, 4) .. Date(2019, 10)], @formula(ret ~ mktrf + hml))
@test coefnames(rr) == ["(Intercept)", "mktrf", "hml"]
@test responsename(rr) == "ret"
@test nobs(rr) == 126
@test all(isapprox.(coef(rr), [-.00125105, 1.40602071, 1.19924984]; atol=.00001))
@test isapprox(r2(rr), .42340085667)
@test isapprox(adjr2(rr), .41402526)
@test dof_residual(rr) == 123
@test islinear(rr)
@test alpha(rr) == rr.coef[1]
@test beta(rr) == rr.coef[2]

##

@test_throws KeyError quick_reg(data[1, Date(2020) .. Date(2021)], @formula(ret ~ mktrf))

rr = quick_reg(data[18428, Date(2019, 4) .. Date(2019, 10)], @formula(ret ~ mktrf + hml); minobs=1000)# not enough data case

@test ismissing(coefnames(rr))
@test ismissing(coef(rr))
@test nobs(rr) == 126
@test ismissing(r2(rr))
@test ismissing(adjr2(rr))
@test alpha(rr) |> ismissing
@test beta(rr) |> ismissing

##

# the SAS code that I tested this against appears to round results to 3 significant digits

# Since these are over specific periods, I specify those as functions to make this easier

event_start(x; data=data) = advancebdays(data.calendar, tobday(data.calendar, x), -10)
event_end(x; data=data) = advancebdays(data.calendar, tobday(data.calendar, x), 10)
est_end(x; data=data) = advancebdays(data.calendar, event_start(x), -16)
est_start(x; data=data) = advancebdays(data.calendar, est_end(x), -149)

rr_market = quick_reg(
    data[df_res.permno, est_start.(df_res.event_date) .. est_end.(df_res.event_date)],
    @formula(ret ~ mktrf)
)

@test isapprox(round.(alpha.(rr_market), digits=5), df_res.alpha_market_model_)
@test isapprox(round.(beta.(rr_market), digits=3), df_res.beta_market_model)
cars = car(data[df_res.permno, event_start.(df_res.event_date) .. event_end.(df_res.event_date)], rr_market)
@test isapprox(round.(cars, sigdigits=3), df_res.car_mm)
bhars = bhar(data[df_res.permno, event_start.(df_res.event_date) .. event_end.(df_res.event_date)], rr_market)
@test isapprox(round.(bhars, sigdigits=3), df_res.bhar_mm)
stds = std.(rr_market)
vars = var.(rr_market)
@test isapprox(round.(vars, digits=10), df_res.estimation_period_variance_market_model_)

##

rr_ff = quick_reg(
    data[df_res.permno, est_start.(df_res.event_date) .. est_end.(df_res.event_date)],
    @formula(ret ~ mktrf + smb + hml)
)

cars = car(data[df_res.permno, event_start.(df_res.event_date) .. event_end.(df_res.event_date)], rr_ff)
@test isapprox(round.(cars, sigdigits=3), df_res.car_ff)
bhars = bhar(data[df_res.permno, event_start.(df_res.event_date) .. event_end.(df_res.event_date)], rr_ff)
@test isapprox(round.(bhars, sigdigits=3), df_res.bhar_ff)
stds = std.(rr_ff)
vars = var.(rr_ff)
@test isapprox(round.(vars, digits=10), df_res.estimation_period_variance_ff_model_)

##

rr_ffm = quick_reg(
    data[df_res.permno, est_start.(df_res.event_date) .. est_end.(df_res.event_date)],
    @formula(ret ~ mktrf + smb + hml + umd)
)

cars = car(data[df_res.permno, event_start.(df_res.event_date) .. event_end.(df_res.event_date)], rr_ffm)
@test isapprox(round.(cars, sigdigits=3), df_res.car_ffm)
bhars = bhar(data[df_res.permno, event_start.(df_res.event_date) .. event_end.(df_res.event_date)], rr_ffm)
@test isapprox(round.(bhars, sigdigits=3), df_res.bhar_ffm)
stds = std.(rr_ffm)
vars = var.(rr_ffm)
@test isapprox(round.(vars, digits=10), df_res.estimation_period_variance_carhart_model_)

##


text:
function get_christoffel(x::TensorHead,TIT::TensorIndexType,metric::T ) where {T <: Tensor}
    """
    returns christoffel symbols as ``Γ_{ijk}``
    """
    @indices TIT i j k l
    g = metric

    simp(x) = contract_metric(x,TIT.metric)
    gg = diff(g(-i,-j),x(k))
    gg = simp(gg)

    h = (gg(-i,-j,-k) - gg(-j,-k,-i) + gg(-k,-i,-j))/2
    h = simp(h)

    return h(-i,-j,-k)
end

function get_riemann(x::TensorHead,h::S,TIT::TensorIndexType,metric::T) where {S <: Tensor, T <: Tensor}
    """
    compute Riemann curvature tensor as ``R_{ijkl}``
    """
    @indices TIT i j k l m n
    g = metric
    simp(x) = contract_metric(x,TIT.metric)
    dh = diff(h(-i,-j,-k),x(l))
    dh = simp(dh)

    hh = (TIT.metric(m,n)*h(-m,-i,-j))*h(-n,-k,-l)
    hh = simp(hh)

    Riemann = dh(-l,-i,-k,-j) - dh(-l,-j,-k,-i) + hh(-j,-k,-i,-l) - hh(-i,-k,-j,-l)
    return Riemann
end


text:
# file for functions related to parallel bookkeeping

# maximum number of matches a vertex can have (for 1, 2, and 3D meshes, 1D meshes
# are not), on any single remote process.  For Matches meshes, this determines
# the size of an array that is sent via MPI
global const Max_Vert_Matches = [0, 4, 8]

global const Max_Vert_Remotes = 1  # maximum number of remotes a vertex can have
                                   # on any single process

### Some Exposition about Pumi
# Pumi has 2 ways of describing relationships between meshentities of the
# same type: matches and remotes.  Matches are used to implement periodic
# boundary conditions and represent two entities that are logically the same
# entity but have different coordinates.  Remotes are used to implement
# ghosting of entities in parallel, and represent the same entity on two
# different processes (APF Parts).
# On top of the matches and remotes machinery is built the apf::Sharing
# machinery.  By default, the apf::Sharing object returns he union of
# the matches and the remotes for a given MeshEntity (although I think it is
# possible to construct a Sharing object that only returns one or the other)
# The functions below use the Sharing object so both matches and remotes are
# accounted for.  The only consideration is that matches, unlike remotes,
# could reside on the same Part as the original MeshEntity.

function getParallelInfo(mesh::PumiMeshDG)
# get information on number of shared entities
# TODO: consider making array of all mesh faces that are shared and iterating
#       over that subsequently.  Because a minority of faces should be shared
#       it might be worthwhile.

  myrank = mesh.myrank
  npeers = apf.countPeers(mesh.m_ptr, mesh.dim-1)  # get edge peers
  peer_nums = zeros(Cint, npeers)
  apf.getPeers(mesh.m_ptr, peer_nums)

  # count the number of edges shared with each peer
  partnums = zeros(Cint, 1)
  remotes = Array{Ptr{Void}}(1)
  counts = zeros(Int, npeers)  # hold the counts of the number of edges shared
                               # with each peer
  for i=1:mesh.numFace
    edge = mesh.faces[i]
    nshares = apf.countCopies(mesh.shr_ptr, edge)
    if nshares > 0
      @assert nshares == 1
      apf.getCopies(partnums, remotes)
      if partnums[1] == mesh.myrank
        continue
      end
      part_idx = getElIndex(peer_nums, partnums[1])
      if part_idx == 0  # this is a periodic neighbor not seen before
        push!(peer_nums, partnums[1])
        push!(counts, 1)
        npeers += 1
      else
        counts[part_idx] += 1
      end
    end
  end
  mesh.peer_face_counts = counts

  mesh.peer_parts = peer_nums
  mesh.npeers =  npeers
  mesh.send_reqs = Array{MPI.Request}(npeers)  # array of Requests for sends
  mesh.send_waited = Array{Bool}(npeers)
  mesh.recv_reqs = Array{MPI.Request}(npeers)  # array of Requests for receives
  mesh.recv_waited = Array{Bool}(npeers)
  mesh.send_stats = Array{MPI.Status}(npeers)
  mesh.recv_stats = Array{MPI.Status}(npeers)

  # for matching meshes, it is possible for one vertex to have several
  # Copies on a given remote process, so we send them all and sort it
  # out on the other side
  if apf.hasMatching(mesh.m_ptr)
    num_orientation_verts = 3*(Max_Vert_Matches[mesh.dim] + Max_Vert_Remotes)
  else
    num_orientation_verts = mesh.dim
  end

  # allocate memory
  edges_local = Array{Array{Ptr{Void}, 1}}(npeers)
  edges_remote = Array{Array{Ptr{Void}, 1}}(npeers)  # is this still used?
  mesh.bndries_local = bndries_local = Array{Array{Boundary,1}}(npeers)
  mesh.bndries_remote = bndries_remote = Array{Array{Boundary,1}}(npeers)
  mesh.shared_interfaces = shared_interfaces = Array{Array{Interface,1}}(npeers)
  orientations_send = Array{Array{Ptr{Void}, 2}}(npeers)
  orientations_recv = Array{Array{Ptr{Void}, 2}}(npeers)
  for i=1:npeers
    edges_local[i] = Array{Ptr{Void}}(counts[i])
    edges_remote[i] = Array{Ptr{Void}}(counts[i])
    bndries_local[i] = Array{Boundary}(counts[i])
    bndries_remote[i] = Array{Boundary}(counts[i])
    orientations_send[i] = Array{Ptr{Void}}(num_orientation_verts, counts[i])
    orientations_recv[i] = Array{Ptr{Void}}(num_orientation_verts, counts[i])
    fill!(orientations_send[i], C_NULL)
    fill!(orientations_recv[i], C_NULL)
  end

  # get all the (local pointers to) edges in a single pass,
  # even though they might not be needed
  curr_pos = ones(Int, npeers)  # hold the current position in each edge array
  for i=1:mesh.numFace
    edge_i = mesh.faces[i]
    nshares = apf.countCopies(mesh.shr_ptr, edge_i)
    if nshares > 0
      @assert nshares == 1
      apf.getCopies(partnums, remotes)
      peer_i = partnums[1]
      if peer_i == mesh.myrank
        continue
      end
      idx = getElIndex(peer_nums, peer_i)  # get the part boundary index
      edges_local[idx][curr_pos[idx]] = edge_i
      curr_pos[idx] += 1
    end
  end

  # get boundary info for the edges
  for i=1:npeers
    if mesh.myrank > peer_nums[i]
      getEdgeBoundaries(mesh, edges_local[i], bndries_local[i])
      sort!(bndries_local[i])
    end
  end

  # get the remote edge pointers if I am the higher number process
  # tell MPI that the pointers are really Ints
  dtype = MPI.mpitype(Int)
  MPI.mpitype_dict[Ptr{Void}] = dtype
  myrank = mesh.myrank
  down = Array{Ptr{Void}}(12)  # hold downward edges
  nrecvs = 0
  for i=1:npeers
    edges_i = edges_remote[i]  # array to hold the remote edge pointers
    if mesh.myrank > peer_nums[i]
      bndries_i = bndries_local[i]
      for j = 1:length(bndries_i)
        bndry_j = bndries_i[j]
        # get the edge
        el_j = mesh.elements[bndry_j.element]
        apf.getDownward(mesh.m_ptr, el_j, mesh.dim-1, down)
        edge_j = down[bndry_j.face]

        # get the remote edge pointer
        nremotes = apf.countCopies(mesh.shr_ptr, edge_j)
        apf.getCopies(partnums, remotes)
        @assert nremotes == 1
        edges_i[j] = remotes[1]
      end

      # now send the edges to the lower rank process
      mesh.send_reqs[i] = MPI.Isend(edges_i, peer_nums[i], 1, mesh.comm)
    else  # I am the receiving process
      mesh.recv_reqs[i] = MPI.Irecv!(edges_i, peer_nums[i], 1, mesh.comm)
      nrecvs += 1
    end
  end  # end loop over peers

  # get arrays of Requests
  recv_reqs_reduced = Array{MPI.Request}(nrecvs)
  recv_peers = Array{Int}(nrecvs)  # which process they came from
  pos = 1
  for i=1:npeers
    if mesh.myrank < peer_nums[i]
      recv_reqs_reduced[pos] = mesh.recv_reqs[i]
      recv_peers[pos] = i
      pos += 1
    end
  end

  # wait for all communication to finish
  #TODO: consider if it is possible to send only the relative orientation, not
  #      the vertices themselves, to the remote process
  for i=1:npeers
    peernum = peer_nums[i]
    if mesh.myrank > peer_nums[i]  # wait for send to finish
      mesh.send_stats[i] = MPI.Wait!(mesh.send_reqs[i])
      getBndryOrientations(mesh, peernum, bndries_local[i], orientations_send[i])
      # get edge orientations
    else  # wait for receive to finish
      j, stat = MPI.Waitany!(recv_reqs_reduced)
      mesh.recv_stats[j] = stat
      peeridx = recv_peers[j]
      peernum = peer_nums[peeridx]
      getEdgeBoundaries(mesh, edges_remote[peeridx], bndries_local[peeridx])
      getBndryOrientations(mesh, peernum, bndries_local[peeridx], orientations_send[peeridx])
      # get edge orientations
    end
  end

  # now send Boundary and orientation info
  send_reqs2 = Array{MPI.Request}(npeers)
  recv_reqs2 = Array{MPI.Request}(npeers)
  MPI.type_create(Boundary)
  MPI.type_create(EntityOrientation)
  for i=1:npeers
    mesh.send_reqs[i] = MPI.Isend(bndries_local[i], peer_nums[i], 1, mesh.comm)
    mesh.recv_reqs[i] = MPI.Irecv!(bndries_remote[i], peer_nums[i], 1, mesh.comm)
    send_reqs2[i] = MPI.Isend(orientations_send[i], peer_nums[i], 2, mesh.comm)
    recv_reqs2[i] = MPI.Irecv!(orientations_recv[i], peer_nums[i], 2, mesh.comm)
  end

  for i=1:npeers
    mesh.recv_stats[i] = MPI.Wait!(mesh.recv_reqs[i])
    MPI.Wait!(recv_reqs2[i])
  end

  # now create Interfaces from the two Boundary arrays
  peer_offsets = Array{Int}(npeers+1)  # record the starting element number
                                       # of the elements belonging to each peer
  mesh.local_element_lists = Array{Array{Int32, 1}}(mesh.npeers)
  curr_elnum = mesh.numEl + 1
  myrank = mesh.myrank

  for i=1:npeers
    peer_offsets[i] = curr_elnum
    curr_elnum, shared_interfaces[i] = numberBoundaryEls(mesh, curr_elnum,
                                 bndries_local[i], bndries_remote[i],
                                 orientations_recv[i])
    curr_elnum += 1

    mesh.local_element_lists[i] = getBoundaryElList(bndries_local[i])
  end
  peer_offsets[npeers+1] = curr_elnum
  mesh.shared_element_offsets = peer_offsets
  mesh.numGlobalEl = curr_elnum - 1
  mesh.numSharedEl = curr_elnum -1 - mesh.numEl  # number of non-local shared
                                                 # elnums
  mesh.local_element_counts = Array{Int}(mesh.npeers)
  mesh.remote_element_counts = Array{Int}(mesh.npeers)
  for i=1:mesh.npeers
    mesh.remote_element_counts[i] = peer_offsets[i+1] - peer_offsets[i]
    mesh.local_element_counts[i] = length(mesh.local_element_lists[i])
  end

  # create the dictonary that maps from locally owned element's numbers
  # to their adjacent non-local neighbors
  adj_dict, revadj = getLocalAdjacency(mesh)
  colordata = ColoringData(adj_dict, revadj)

  # delete the dangerous Ptr{Void} -> Int MPI type
  delete!(MPI.mpitype_dict, Ptr{Void})

  # wait for the sends to finish before exiting
  for i=1:npeers
    mesh.send_stats[i] = MPI.Wait!(mesh.send_reqs[i])
    MPI.Wait!(send_reqs2[i])
  end

  return colordata
end

function getEdgeBoundaries(mesh::PumiMeshDG, edges::Array{Ptr{Void}},
                           bndries::Array{Boundary})
# get the array of Boundaryies for an array of edges
# edges is the array of the edge MeshEnities
# bndries is the array to be populated with the Boundaryies

  faces = Array{Ptr{Void}}(1)

  for i=1:length(edges)
    edge_i = edges[i]

    numFace = apf.countAdjacent(mesh.m_ptr, edge_i, mesh.dim)  # should be count upward

    @assert( numFace == 1)

    apf.getAdjacent(faces)
    facenum = apf.getNumberJ(mesh.el_Nptr, faces[1], 0, 0) + 1
    edgenum = apf.getNumberJ(mesh.face_Nptr, edge_i, 0, 0) + 1
#    facenum = getFaceNumber2(faces[1]) + 1
#    edgenum = getEdgeNumber2(edge_i) + 1  # unneeded?
    edgenum_local = getFaceLocalNum(mesh, edgenum, facenum)

    # verify
    down, numdown = apf.getDownward(mesh.m_ptr, mesh.elements[facenum], mesh.dim-1)
    @assert down[edgenum_local] == edge_i

    bndries[i] = Boundary(facenum, edgenum_local)

    ncopies = apf.countCopies(mesh.shr_ptr, edge_i)
    if ncopies == 0
      println("face ", facenums, " has zero copies")
    end
  end

  return nothing
end

"""
  Gets all of the remote pointers (both matches and remotes) for the
  specified edges/faces on the specified peer process.

  **Inputs**

   * mesh
   * peer_num: the MPI rank of the peer process
   * bndries: array of Boundary objects that identify the faces

  **Inputs/Outputs**

   * orientations: an array of length n x length(bndries), where n is large
                   enough to store all the remotes/matches on the peer process.

  **Implementation Notes**

  This function asserts of n is not large enough
"""
function getBndryOrientations(mesh::PumiMeshDG, peer_num::Integer, bndries::AbstractArray{Boundary},
                             orientations::AbstractArray{Ptr{Void}, 2})
# gets the remote pointers to the vertices of the elements on the boundaries
# peer_num is the MPI rank of the peer process

  nfaces = length(bndries)
  myrank = mesh.myrank
  downward_verts = Array{Ptr{Void}}(4)
  downward_faces = Array{Ptr{Void}}(12)
  remote_partnums = Array{Cint}(400)  # equivalent of apf::up
  remote_ptrs = Array{Ptr{Void}}(400)
  vertmap = mesh.topo.face_verts

  for i=1:nfaces
    bndry_i = bndries[i]
    el_i = mesh.elements[bndry_i.element]
    facelocal_i = bndry_i.face
    apf.getDownward(mesh.m_ptr, el_i, 0, downward_verts)
    verts_startidx = 1

    # figure out if this face has a remote or a match (it can't have both)
    apf.getDownward(mesh.m_ptr, el_i, mesh.dim - 1, downward_faces)
    face_i = downward_faces[facelocal_i]

    # because we are dealing with faces, we don't need to check the part
    # number of the pair
    nremotes = apf.countRemotes(mesh.m_ptr, face_i)
    nmatches = apf.countMatches(mesh.m_ptr, face_i)

    is_remote = nremotes > 0
    @assert is_remote || nmatches > 0

    for j=1:(mesh.dim)  # dim = number of verts on a face

      vert_j = downward_verts[vertmap[j, facelocal_i]]
      # get the remote pointer for this vert
      # depending on whether the face has a remote or a match, send only
      # the remote or matched vertices
      # this helps resolve ambiguities on the receiving end when determining
      # orientation in some corner cases
      if is_remote
        nremotes = apf.countRemotes(mesh.m_ptr, vert_j)
        apf.getRemotes(remote_partnums, remote_ptrs)
      else
        nremotes = apf.countMatches(mesh.m_ptr, vert_j)
        apf.getMatches(remote_partnums, remote_ptrs)
      end

      @assert nremotes <= 400
      @assert nremotes >= 1

      # get the remote verts on the specified peer
      partnums_extract = sview(remote_partnums, 1:Int(nremotes))
      ptrs_extract = sview(remote_ptrs, 1:Int(nremotes))
      verts_j = sview(orientations, verts_startidx:size(orientations, 1), i)
      verts_startidx += getVertCopies(partnums_extract, ptrs_extract, peer_num, verts_j, mesh.f)

    end
  end

  return nothing
end

"""
  Given the remote parts numbers and the associated MeshEntity pointers, copies
  the pointers into the output array for only those MeshEntities on the
  specified peer process.

  **Inputs**

   * remote_partnums: array of part numbers (MPI ranks)
   * remote_ptrs: MeshEntity pointers (for each remote part)
   * peer_num: MPI rank of the process to get the MeshEntity pointers for
   * f: IO stream (unused?)

  **Inputs/Outputs**

   * vert_copies: array to put the selected MeshEntity pointers into

  **Outputs**

   * number of MeshEntity pointers copied to the output array
"""
function getVertCopies(remote_partnums::AbstractArray, remote_ptrs::AbstractArray, peer_num::Integer, vert_copies::AbstractArray, f=STDOUT)
# get all the vertices on peer peer_num and put them in the array
# returns the number of vertices inserted

  pos = 1
  for i=1:length(remote_partnums)
    if remote_partnums[i] == peer_num
      @assert pos <= length(vert_copies)
      vert_copies[pos] = remote_ptrs[i]
      pos += 1
    end
  end

  return pos - 1
end

function extractVertCopies(recv_verts::AbstractArray, local_verts::AbstractArray, facevertsR::AbstractArray)
# extract the vertices that are the same as the local_verts from recv_verts,
# preserving their order
# this is the inverse function to the way getVertCopies is used

  pos = 1
  for i=1:length(recv_verts)
    recv_vert = recv_verts[i]
    for j=1:length(local_verts)
      if recv_vert == local_verts[j]
        facevertsR[pos] = recv_vert
        pos += 1
      end
    end
  end

  @assert pos == 4  # check that we found all vertices

  return nothing
end

"""
  Assert that all entries of the array are unique

  Inputs:
    arr: array of values

"""
function assertUnique(arr::AbstractArray)

  n = length(arr)

  for i=1:n
    val_i = arr[i]
    for j=1:(i-1)
      @assert val_i != arr[j]
    end
  end

  return nothing
end

"""
  This function ensures that if an entry in matched_entities is present in
  other entities, then it occurs exactly once in matched_entities.  This ensures
  that vertices can be matched up to face orientation can be determined

  **Inputs**

   * matched_entities: array of entities
   * other_entities: array of other entities
"""
function assertUnique(matched_entities::AbstractArray,
                      other_entities::AbstractArray)

  for i=1:length(matched_entities)
    if (matched_entities[i] in other_entities)
      @assert countin(matched_entities[i], matched_entities) == 1
    end
  end

  return nothing
end

# this can be generalized once edge orientation is generalized
function numberBoundaryEls(mesh, startnum, bndries_local::Array{Boundary},
                           bndries_remote::Array{Boundary},
                           orientations_recv::AbstractArray{Ptr{Void}, 2},
                           f=STDOUT)
# create Interfaces out of the local + remote Boundary arrays
# also numbers the remote elements with numbers > numEl, storing them in
# the elementR field of the Interface

  ninterfaces = length(bndries_local)
  interfaces = Array{Interface}(ninterfaces)
  curr_elnum = startnum  # counter for
  new_elR = 0 # number of new element to create
  face_vertmap = mesh.topo.face_verts
  el_verts = Array{Ptr{Void}}(4)
  face_verts = Array{Ptr{Void}}(3)
  facevertsR = Array{Ptr{Void}}(3)
  for i=1:ninterfaces
    bndry_l = bndries_local[i]
    bndry_r = bndries_remote[i]

    # figure out elementR number
    old_iface_idx = isRepeated(bndries_remote, i)
    if old_iface_idx == 0
      new_elR = curr_elnum
      curr_elnum += 1
    else
      new_elR = interfaces[old_iface_idx].elementR
    end

    # figure out orientation
    if mesh.dim == 2
      orient = 1
    else
      # get the local face vertices
      el_ptr = mesh.elements[bndry_l.element]
      face_local = bndry_l.face
      apf.getDownward(mesh.m_ptr, el_ptr, 0, el_verts)
      for j=1:3
        face_verts[j] = el_verts[face_vertmap[j, face_local]]
      end
      faceverts_recv = sview(orientations_recv, :, i)
      assertUnique(faceverts_recv, face_verts)
      extractVertCopies(faceverts_recv, face_verts, facevertsR)
      orient = calcRelativeOrientation(face_verts, facevertsR)
    end

    interfaces[i] = Interface(bndry_l.element, new_elR,  bndry_l.face, bndry_r.face, UInt8(orient))
  end

  last_elnum = curr_elnum - 1
  return last_elnum, interfaces
end

function getBoundaryElList(bndries_local::Array{Boundary}, f=STDOUT)
# get the list of elemements on the boundary
  nfaces = length(bndries_local)
  elnums = Array{Int32}(nfaces)  # nfaces is the upper bound on
                                 # the number of elements
  pos = 1 # current position in elnums
  for i=1:nfaces
    old_iface_idx = isRepeated(bndries_local, i)
    if old_iface_idx == 0
      elnums[pos] = bndries_local[i].element
      pos += 1
    end
  end

  return elnums[1:(pos-1)]
end

function isRepeated(bndries::Array{Boundary}, idx)
# figures out if a specified element has been seen before, returning
# the index of the boundary containing it if so, and returning zero if not
  el = bndries[idx].element
  for i=1:(idx-1)
    if bndries[i].element == el
      return i
    end
  end

  return 0
end

# could be generalized with a few constants
function getLocalAdjacency(mesh::PumiMeshDG)

  # map from an element number of all the element numbers of the
  # non-local elements
  # an element can have a maximum of 2 non-local neighbor elements
  adj_dict = Dict{Int, Array{Int, 1}}()
  sizehint!(adj_dict, mesh.numSharedEl)

  # map from non-local element to local neighbors
  revadj = zeros(Int, mesh.numSharedEl, 5)

  for i=1:mesh.npeers
    ifaces_i = mesh.shared_interfaces[i]
    for j=1:length(ifaces_i)
      iface_j = ifaces_i[j]
      el_local = iface_j.elementL
      el_nonlocal = iface_j.elementR

      if !haskey(adj_dict, el_local)
        adj_dict[el_local] = Int[el_nonlocal, 0, 0, 0, 0]
      else
        old_tuple = adj_dict[el_local]
        @assert old_tuple[5] == 0
        idx = findfirst(old_tuple, 0)
        old_tuple[idx] = Int(el_nonlocal)
#        adj_dict[el_local] = (old_tuple[1], Int(el_nonlocal))
      end  # end if-else

      nonlocal_idx = el_nonlocal - mesh.numEl
      for k=1:5
        if revadj[nonlocal_idx, k] == 0
          revadj[nonlocal_idx, k] = el_local
          break
        end
      end
      #=
      first_neighbor = revadj[nonlocal_idx, 1]
      second_neighbor = revadj[nonlocal_idx, 2]
      if first_neighbor == 0
        revadj[nonlocal_idx, 1] = el_local
      elseif second_neighbor == 0
        revadj[nonlocal_idx, 2] = el_local
      else
        throw(ErrorException("Too many adjacent elements for peer $i interface $j"))
      end
      =#

    end # end loop over interfaces
  end  # end loop over peers

  return adj_dict, revadj
end

function getVertexParallelInfo(mesh::PumiMeshDG)
# get lists of shared vertices in a mutually agreed order
# this will enable sending data back and forth

  myrank = mesh.myrank
  npeers = apf.countPeers(mesh.m_ptr, 0)  # count vertex peers
  peer_nums = zeros(Cint, npeers)
  apf.getPeers(mesh.m_ptr, peer_nums)  # get the identifiers of the peers
  # storage for all shares of a vertex
  partnums = zeros(Cint, 400 + Max_Vert_Matches[mesh.dim])
  remotes = Array{Ptr{Void}}(400 + Max_Vert_Matches[mesh.dim])
  counts = zeros(Int, npeers)  # count number of verts shared with each peer

  for i=1:mesh.numVert
    vert = mesh.verts[i]
    nshares = apf.countCopies(mesh.shr_ptr, vert)
    if nshares > 0
      @assert nshares < (400 + Max_Vert_Matches[mesh.dim])
      apf.getCopies(partnums, remotes)

      for j=1:nshares
        if partnums[j] == myrank  # this is not shared with another part, ignore
          continue
        end

        part_idx = getElIndex(peer_nums, partnums[j])
        if part_idx == 0  # this part has not been seen before and it was not
                          # on the list from apf.getPeers(), so it must be
                          # a periodic share

          push!(peer_nums, partnums[j])
          push!(counts, 1)
          npeers += 1
        else  # part has been seen before
          counts[part_idx] += 1
        end
      end  # end loop j
    end  # end if nshares > 0
  end  # end loop i

  # now we know how many vertices are shared with whom
  # allocate arrays of the right size
  verts_local = Array{Array{Ptr{Void}, 1}}(npeers)  # local pointers to shared
                                                    # verts
  verts_remote = Array{Array{Ptr{Void}, 1}}(npeers) # remote pointers to shared
                                                    # verts
  for i=1:npeers
    verts_local[i] = Array{Ptr{Void}}(counts[i])
    #TODO: dont allocate unneeded arrays in verts_remote
    verts_remote[i] = Array{Ptr{Void}}(counts[i])
  end
  curr_pos = ones(Int, npeers)  # current position in each array inside
                                # verts_local
#  verts_remotes = Array{Array{Ptr{Void}, 1}}(npeers)
  for j=1:mesh.numVert
    vert_j = mesh.verts[j]
    nshares = apf.countCopies(mesh.shr_ptr, vert_j)
    if nshares > 0
      apf.getCopies(partnums, remotes)
      for k=1:nshares
        if partnums[k] == myrank
          continue
        end
        peer_k = partnums[k]
        vert_k = remotes[k]
        # if my rank is greater add the remote pointer to the list
        if myrank > peer_k
          peer_idx = getElIndex(peer_nums, peer_k)
          verts_local[peer_idx][curr_pos[peer_idx]] = vert_j
          verts_remote[peer_idx][curr_pos[peer_idx]] = vert_k
          curr_pos[peer_idx] += 1
        end
      end  # end loop k
    end  # end if nshares > 0
  end  # end loop j

  # hold the MPI.Requests for the send/recive
  # either process either sends or receives to each peer, so store all
  # send and receive requests in the same array
  all_reqs = Array{MPI.Request}(npeers)

  # tell MPI the Ptr{Voids} are really Ints
  dtype = MPI.mpitype(Int)
  MPI.mpitype_dict[Ptr{Void}] = dtype

  for i=1:npeers
    peernum = peer_nums[i]

    # send the remote verts pointers to lower process numbers
    # recieve the data into verts_local, which wasn't populated for
    # lower process numbers
    if myrank > peernum
      all_reqs[i] = MPI.Isend(verts_remote[i], peernum, 1, mesh.comm)
    else
      all_reqs[i] = MPI.Irecv!(verts_local[i], peernum, 1, mesh.comm)
    end
  end

  # wait for the receives to finish
  for i=1:npeers
    peernum = peer_nums[i]
    if myrank < peernum
      MPI.Wait!(all_reqs[i])
    end
  end

  vert_nums, rev_mapping = getVertReverseMapping(mesh, peer_nums, counts, verts_local)

  vshare = VertSharing(npeers, peer_nums, counts, vert_nums, rev_mapping)

  # wait for all sends to finish before exiting
  for i=1:npeers
    if myrank > peer_nums[i]
      MPI.Wait!(all_reqs[i])
    end
  end

  # delete the dangerous Ptr{Void} -> Int MPI type
  delete!(MPI.mpitype_dict, Ptr{Void})

  return vshare
end

# get the vertex numbers of all the shared verts (per peer)
# also get a dictionary that maps from local vert number to all the
# (part number, vert idx)s
function getVertReverseMapping(mesh::PumiMeshDG, peer_nums::Array{Cint, 1}, counts::Array{Int, 1}, verts_local::Array{Array{Ptr{Void}, 1}, 1})

  npeers = length(peer_nums)
  vert_nums = Array{Array{Int, 1}}(npeers)
  rev_mapping = Dict{Int, Pair{Array{Cint, 1}, Array{Int, 1}}}()

  if npeers == 0
    return vert_nums, rev_mapping
  end

  sizehint!(rev_mapping, maximum(counts))  # this is an underestimate, but there is no
                              # way to know the upper bound until after the
                              # dictionary is populated
  for i=1:npeers
    vert_nums[i] = Array{Int}(counts[i])
    vert_nums_i = vert_nums[i]
    verts_local_i = verts_local[i]
    for j=1:counts[i]
      vert_j = verts_local_i[j]
      vertnum_j = apf.getNumberJ(mesh.vert_Nptr, vert_j, 0, 0) + 1

      vert_nums_i[j] = vertnum_j

      if !haskey(rev_mapping, vertnum_j)
        rev_mapping[vertnum_j] = Pair(Cint[peer_nums[i]], Int[j])
      else
        oldpair = rev_mapping[vertnum_j]
        push!(oldpair.first, peer_nums[i])
        push!(oldpair.second, j)
      end
    end
  end

  return vert_nums, rev_mapping
end


text:
# Get residuals: (predicted - observed) data
lm_get_residual!(r::VecIO, μ::VecI, f::AbstractObjective) = get_residual!(r, μ, f)

# Get square sum of residuals
lm_get_squaresum(β::Real, r::VecI, n::Int) = 0.5 * β * dot(r, r, n)
lm_get_squaresum(Λ::MatI, r::VecI, n::Int) = 0.5 * dot(r, n, Λ, r, n)

function lm_get_squaresum(λ::VecI, r::VecI, n::Int)
    R = 0.0
    m = mod(n, 4)
    if m ≠ 0
        @inbounds for i in 1:m
            R += λ[i] * abs2(r[i])
        end
        n < 4 && return R
    end
    m += 1
    @inbounds for i in m:4:n
        R += λ[i] * abs2(r[i]) + λ[i+1] * abs2(r[i+1]) + λ[i+2] * abs2(r[i+2]) + λ[i+3] * abs2(r[i+3])
    end
    return 0.5 * R
end

# Get jacobian matrix of predicted model
lm_get_jacobian!(J::MatIO, μ::VecI, f::AbstractObjective) = get_jacobian!(J, μ, f)

# Get approximated hessian of least-square residuals
function lm_get_apprhess!(H::MatIO, β::Real, J::MatI)
    BLAS.syrk!('U', 'T', β, J, 0.0, H)
    @inbounds for j in axes(H, 1), i in 1:j-1
        H[j,i] = H[i,j]
    end
    return nothing
end

function lm_get_apprhess!(H::MatIO, λy::VecI, J::MatI, ΛJ::MatB)
    for j in axes(J, 2)
        @simd for i in eachindex(λy)
            @inbounds ΛJ[i,j] = λy[i] * J[i,j]
        end
    end
    BLAS.gemm!('T', 'N', 1.0, ΛJ, J, 0.0, H) # H ← JΛ * J
    return nothing
end

function lm_get_apprhess!(H::MatIO, Λy::MatI, J::MatI, ΛJ::MatB)
    BLAS.symm!('L', 'U', 1.0, Λy, J, 0.0, ΛJ) # ΛJ ← Λy * J
    BLAS.gemm!('T', 'N', 1.0, ΛJ, J, 0.0,  H) # H  ← JΛ * J
    return nothing
end

# Get 1st-order least-square gradient
lm_get_gradient!(d::VecIO, β::Real, J::MatI, r::VecI) = BLAS.gemv!('T', β, J, r, 0.0, d)
lm_get_gradient!(d::VecIO, ΛJ::MatI, r::VecI)         = BLAS.gemv!('T', 1.0, ΛJ, r, 0.0, d)

# Get 2nd-order directional gradient
function lm_get_gradient!(n::VecIO, r::VecI, J::MatI, d::VecI, ih1::Real, ih2::Real)
    BLAS.axpy!(-1.0, r, n)
    BLAS.gemv!('N', ih1, J, d, ih2, n)
    return nothing
end

function lm_get_cholesky!(A::MatIO)
    _, cholesky_state = LAPACK.potrf!('L', A)
    return cholesky_state
end

# Perform BLAS.trsv!
function lm_get_blastrsv!(d::VecIO, A::MatIO)
    BLAS.trsv!('L', 'N', 'N', A, d)
    BLAS.trsv!('L', 'T', 'N', A, d)
    return nothing
end

struct LevenbergMarquardtOptimizer
    @def prop Vector{Float64} rs rt μs μt dμ δμ
    @def prop Matrix{Float64} Js ΛJ Hs Hf
    @def prop Int nd ny

    function LevenbergMarquardtOptimizer(nd::Int, ny::Int)
        @def vars Vector{Float64}(undef, ny) rs rt
        @def vars Vector{Float64}(undef, nd) μs μt dμ δμ
        @def vars Matrix{Float64}(undef, ny, nd) Js ΛJ
        @def vars Matrix{Float64}(undef, nd, nd) Hs Hf # Hessian (solution, factorization)
        return new(rs, rt, μs, μt, dμ, δμ, Js, ΛJ, Hs, Hf, nd, ny)
    end
end

# *t := temp one of *
# *s := solution of *
function lm_get_optimize!(xs::VecIO, fn::AbstractObjective, μ0::VecI, βy::Real, τ::Real, h::Real, itmax::Int, lm::LevenbergMarquardtOptimizer)
    #### Allocations ####
    @get lm nd ny rs rt μs μt dμ δμ
    @get lm Js Hs Hf
    #### Initialization
    ih1 = inv(h)
    ih2 = ih1 * ih1
    toN = eachindex(1:nd)
    unsafe_copy!(μs, μ0)
    #### Compute precision matrix
    lm_get_residual!(rs, μs, fn)
    lm_get_jacobian!(Js, μs, fn)
    lm_get_apprhess!(Hs, βy, Js)
    # Controlling parameters of trust region
    a = τ * diagmax(Hs, nd); b = 2.0
    #### Compute least square cost
    Cnow = lm_get_squaresum(βy, rs, ny)
    #### Iteration
    it = 0
    while it < itmax
        it += 1
        # Damped Gauss-Newton approximation
        unsafe_copy!(Hf, Hs)
        λ = a * diagmax(Hf, nd)
        @simd for i in toN
            @inbounds Hf[i,i] += λ
        end
        cholesky_state = lm_get_cholesky!(Hf)
        # Levenberg-Marquardt step by 1st order linearization
        lm_get_gradient!(dμ, βy, Js, rs)
        lm_get_blastrsv!(dμ, Hf)
        # Finite-difference of 2nd order directional derivative
        @simd for i in toN
            @inbounds μt[i] = μs[i] + h * dμ[i]
        end
        lm_get_residual!(rt, μt, fn)
        lm_get_gradient!(rt, rs, Js, dμ, ih1, ih2)
        # Geodesic acceleration step by 2nd order linearization
        lm_get_gradient!(δμ, βy, Js, rt)
        lm_get_blastrsv!(δμ, Hf)
        gratio = 2.0 * nrm2(δμ) / nrm2(dμ)
        gratio > 0.75 && (a *= b; b *= 2.0; continue)
        # (Levenberg-Marquardt velocity) + (Geodesic acceleration)
        BLAS.axpy!(1.0, δμ, dμ) # dμ ← dμ + δμ
        # Compute predicted gain of the least square cost by linear approximation
        LinApprox = lm_get_squaresum(Hs, dμ, nd) + λ * dot(dμ, dμ, nd)
        # Compute gain ratio
        @simd for i in toN
            @inbounds μt[i] = μs[i] + dμ[i]
        end
        lm_get_residual!(rt, μt, fn)
        Cnew = lm_get_squaresum(βy, rt, ny)

        Δ = Cnow - Cnew; ρ = Δ / LinApprox
        if ρ > 0.0
            Cnow = Cnew
            unsafe_copy!(μs, μt)
            (maximum(δμ) < 1e-15 || Δ < 1e-10) && break # Check localized convergence
            lm_get_jacobian!(Js, μs, fn)
            lm_get_apprhess!(Hs, βy, Js)
            unsafe_copy!(rs, rt)
            a = ρ < 0.9367902323681495 ? a * (1.0 - cubic(2.0 * ρ - 1.0)) : a / 3.0
            b = 2.0
        else
            a > 1e10 && break # Check localized divergence
            a *= b; b *= 2.0
        end
    end
    unsafe_copy!(xs, μs); return xs
end

function lm_get_optimize!(xs::VecIO, fn::AbstractObjective, μ0::VecI, Λy::Union{VecI, MatI}, τ::Real, h::Real, itmax::Int, lm::LevenbergMarquardtOptimizer)
    #### Allocations ####
    @get lm nd ny rs rt μs μt dμ δμ
    @get lm Js ΛJ Hs Hf
    #### Initialization
    ih1 = inv(h)
    ih2 = ih1 * ih1
    toN = eachindex(1:nd)
    unsafe_copy!(μs, μ0)
    #### Compute precision matrix
    lm_get_residual!(rs, μs, fn)
    lm_get_jacobian!(Js, μs, fn)
    lm_get_apprhess!(Hs, Λy, Js, ΛJ) # ΛJ = Λy * Js
    # Controlling parameters of trust region
    a = τ * diagmax(Hs, nd); b = 2.0
    #### Compute least square cost
    Cnow = lm_get_squaresum(Λy, rs, ny)
    #### Iteration
    it = 0
    while it < itmax
        it += 1
        # Damped Gauss-Newton approximation
        unsafe_copy!(Hf, Hs)
        λ = a * diagmax(Hf, nd)
        @simd for i in toN
            @inbounds Hf[i,i] += λ
        end
        cholesky_state = lm_get_cholesky!(Hf)
        # Levenberg-Marquardt step by 1st order linearization
        lm_get_gradient!(dμ, ΛJ, rs)
        lm_get_blastrsv!(dμ, Hf)
        # Finite-difference of 2nd order directional derivative
        @simd for i in toN
            @inbounds μt[i] = μs[i] + h * dμ[i]
        end
        lm_get_residual!(rt, μt, fn)
        lm_get_gradient!(rt, rs, Js, dμ, ih1, ih2)
        # Geodesic acceleration step by 2nd order linearization
        lm_get_gradient!(δμ, ΛJ, rt)
        lm_get_blastrsv!(δμ, Hf)
        gratio = 2.0 * nrm2(δμ) / nrm2(dμ)
        gratio > 0.75 && (a *= b; b *= 2.0; continue)
        # (Levenberg-Marquardt velocity) + (Geodesic acceleration)
        BLAS.axpy!(1.0, δμ, dμ) # dμ ← dμ + δμ
        # Compute predicted gain of the least square cost by linear approximation
        LinApprox = lm_get_squaresum(Hs, dμ, nd) + λ * dot(dμ, dμ, nd)
        # Compute gain ratio
        @simd for i in toN
            @inbounds μt[i] = μs[i] + dμ[i]
        end
        lm_get_residual!(rt, μt, fn)
        Cnew = lm_get_squaresum(Λy, rt, ny)

        Δ = Cnow - Cnew; ρ = Δ / LinApprox
        if ρ > 0.0
            Cnow = Cnew
            unsafe_copy!(μs, μt)
            (maximum(δμ) < 1e-15 || Δ < 1e-10) && break # Check localized convergence
            lm_get_jacobian!(Js, μs, fn)
            lm_get_apprhess!(Hs, Λy, Js, ΛJ)
            unsafe_copy!(rs, rt)
            a = ρ < 0.9367902323681495 ? a * (1.0 - cubic(2.0 * ρ - 1.0)) : a / 3.0
            b = 2.0
        else
            a > 1e10 && break # Check localized divergence
            a *= b; b *= 2.0
        end
    end
    unsafe_copy!(xs, μs); return xs
end

optimize!(xs::VecIO, fn::AbstractObjective, μ0::VecI, Λy::Union{Real, VecI, MatI}, lm::LevenbergMarquardtOptimizer; τ::Real=1e-3, h::Real=0.1, itmax::Int=100) = lm_get_optimize!(xs, fn, μ0, Λy, τ, h, itmax, lm)


text:
import Base.show

@enum CONTROL_OP WAIT SYNC LOAD_REPEAT REPEAT GOTO CALL RETURN

immutable ControlFlow
	label::String
	op::CONTROL_OP
	target
	value
end

wait() = ControlFlow("WAIT", WAIT, 0, 0)
sync() = ControlFlow("SYNC", SYNC, 0, 0)
goto(target) = ControlFlow("GOTO", GOTO, target, 0)

function show(io::IO, cf::ControlFlow)
	print(io, cf.label)
	if cf.op == GOTO
		print(io, ": $(cf.target)")
	end
end


text:
"""
    Fermi.GaussianBasis

Module handling basis functions and basis set.
"""
module GaussianBasis

using Fermi
using Fermi.Geometry
using Fermi.Options
using Memoize
using Formatting

import Base: getindex, show

export BasisSet, BasisFunction

include("BasisParser.jl")

@doc raw"""
    Fermi.GaussianBasis.BasisFunction

Object representing a shell of Gaussian basis functions composed of ``N`` primitives:

```math
\chi_{l,m} = \sum_i^N C_i r^{l}e^{-\zeta_i r^2} Y_{l,m}(\theta,\phi)
```

# Fields

| Name    | Type | Description |
|:--------|:------|:-----------------------------------------------------------|
|`l`      |`Int32`            | Angular momentum number (e.g. 0, 1, 2 for S, P, D...)      |
|`coef`   |`Array{Float64,1}` | Array with coefficients (``C_i``) for each primitive       |
|`exp`    |`Array{Float64,1}` | Array with exponents (``\zeta_i``) for each primitive      |

# Examples

```julia
julia> bf = Fermi.GaussianBasis.BasisFunction(1, [1/√2, 1/√2], [5.0, 1.2])
P shell with 3 basis built from 2 primitive gaussians

χ₁₋₁ =    0.7071067812⋅Y₁₋₁⋅r¹⋅exp(-5.0⋅r²)
     +    0.7071067812⋅Y₁₋₁⋅r¹⋅exp(-1.2⋅r²)

χ₁₀  =    0.7071067812⋅Y₁₀⋅r¹⋅exp(-5.0⋅r²)
     +    0.7071067812⋅Y₁₀⋅r¹⋅exp(-1.2⋅r²)

χ₁₁  =    0.7071067812⋅Y₁₁⋅r¹⋅exp(-5.0⋅r²)
     +    0.7071067812⋅Y₁₁⋅r¹⋅exp(-1.2⋅r²)
```
"""
struct BasisFunction
    l::Cint
    coef::Array{Cdouble,1}
    exp::Array{Cdouble,1}
end

"""
    Fermi.GaussianBasis.BasisSet

Object holding a set of BasisFunction objects associated with each atom in a molecule.

# Fields

| Name        | Type                               |   Description |
|:------------|:-----------------------------------|:-----------------------------------------------------------|
|`molecule`   | `Molecule`                         | Fermi.Geometry.Molecule object  |
|`basis_name` | `String`                           | String holding the basis set name  |
|`basis`      | `Dict{Atom,Array{BasisFunction,1}}`| A dictionary that maps Atom objects to an Array of BasisFunction  |
|`natoms`     | `Int32`                            | Number of atoms in the BasisSet |
|`nbas`       | `Int32`                            | Number of basis functions (Note, not equal the number of BasisFunction objects) |
|`nshells`    | `Int32`                            | Number of shells, i.e. BasisFunction objects |
|`lc_atoms`   | `Array{Int32,1}`                   | Integer array maping data to libcint |
|`lc_bas`     | `Array{Int32,1}`                   | Integer array maping data to libcint |
|`lc_env`     | `Array{Float64,1}`                 | Float64 array maping data to libcint |

# Example

Build a basis set from default options
```julia
julia> bset = Fermi.GaussianBasis.BasisSet()
sto-3g Basis Set
Number of shells: 5
Number of basis:  7

O: 1s 2s 1p
H: 1s
H: 1s
```
The BasisSet object can be accessed as two-dimensional array.
```julia
julia> bset[2] # Show all basis for the second atom (H 1s)
1-element Vector{Fermi.GaussianBasis.BasisFunction}:
 S shell with 1 basis built from 3 primitive gaussians

χ₀₀  =    0.9817067283⋅Y₀₀⋅exp(-3.425250914⋅r²)
     +    0.9494640079⋅Y₀₀⋅exp(-0.6239137298⋅r²)
     +    0.2959064597⋅Y₀₀⋅exp(-0.168855404⋅r²)
julia> bset[1,2] # Show the second basis for the first atom (O 2s)
S shell with 1 basis built from 3 primitive gaussians

χ₀₀  =    0.8486970052⋅Y₀₀⋅exp(-5.033151319⋅r²)
     +    1.1352008076⋅Y₀₀⋅exp(-1.169596125⋅r²)
     +    0.8567529838⋅Y₀₀⋅exp(-0.38038896⋅r²)
```
You can also create your own crazy mix!
```julia
julia> @molecule {
    H 0.0 0.0 0.0
    H 0.0 0.0 0.7
}
julia> mol = Fermi.Geometry.Molecule()
# Lets create an S basis function for H
julia> s = Fermi.GaussianBasis.BasisFunction(0, [0.5215367271], [0.122])
# Now a P basis function
julia> p = Fermi.GaussianBasis.BasisFunction(1, [1.9584045349], [0.727])
# Now to create the basis set we need a dictionary maping atoms do basis functions
# Let's fetch the atom list from the molecule object
julia> Atoms = mol.atoms
# Now create the desired mapping as a dictionary. We want the first hydrogen
# to use only the S function, whereas the second will use both S and P
julia> shells = Dict(
    Atoms[1] => [s],
    Atoms[2] => [s,p]
)
# Now we can create the basis set!
julia> Fermi.GaussianBasis.BasisSet(mol, "My Crazy", shells)
My Crazy Basis Set
Number of shells: 3
Number of basis:  5

H: 1s
H: 1s 1p
```
"""
struct BasisSet
    molecule::Molecule
    basis_name::String
    basis::Dict{Atom,Array{BasisFunction,1}}
    natoms::Cint
    nbas::Cint
    nshells::Cint
    lc_atoms::Array{Cint,1}
    lc_bas::Array{Cint, 1}
    lc_env::Array{Cdouble,1}
end

function BasisSet()
    mol = Molecule()
    bname = Options.get("basis")
    BasisSet(mol, bname)
end

function BasisSet(mol::Molecule, basis_name::String)

    output("\n  ⇒ Preparing new Basis Set: {:s}", basis_name)
    shells = Dict{Atom, Array{BasisFunction,1}}()
    for A in mol.atoms
        shells[A] = read_basisset(basis_name, A.AtomicSymbol)
    end
    BasisSet(mol, basis_name, shells)
end

function BasisSet(mol::Molecule, basis_name::String, shells::Dict{Atom, Array{BasisFunction,1}})

    atoms = mol.atoms
    ATM_SLOTS = 6
    BAS_SLOTS = 8

    natm = length(atoms)
    nbas = 0
    nshells = 0
    nexps = 0
    nprims = 0

    for A in atoms
        basis = shells[A]
        for b in basis
            nshells += 1
            nbas += 2*b.l + 1
            nexps += length(b.exp)
            nprims += length(b.coef)
        end
    end

    lc_atm = zeros(Cint, natm*ATM_SLOTS)
    lc_bas = zeros(Cint, nshells*BAS_SLOTS)
    env = zeros(Cdouble, 3*natm+nexps+nprims)

    # Prepare the lc_atom input
    off = 0
    ib = 0
    for i = eachindex(atoms)
        A = atoms[i]
        # lc_atom has ATM_SLOTS (6) "spaces" for each atom
        # The first one (Z_INDEX) is the atomic number
        lc_atm[1 + ATM_SLOTS*(i-1)] = Cint(A.Z)
        # The second one is the env index address for xyz
        lc_atm[2 + ATM_SLOTS*(i-1)] = off
        env[off+1:off+3] .= A.xyz ./ Fermi.PhysicalConstants.bohr_to_angstrom
        off += 3
        # The remaining 4 slots are zero.

        # Prepare the lc_bas input
        for j = eachindex(shells[A])
            B = shells[A][j]
            Ne = length(B.exp)
            Nc = length(B.coef)
            # lc_bas has BAS_SLOTS for each basis set
            # The first one is the index of the atom starting from 0
            lc_bas[1 + BAS_SLOTS*ib] = i-1
            # The second one is the angular momentum
            lc_bas[2 + BAS_SLOTS*ib] = B.l
            # The third is the number of primitive functions
            lc_bas[3 + BAS_SLOTS*ib] = Nc
            # The fourth is the number of contracted functions
            lc_bas[4 + BAS_SLOTS*ib] = 1
            # The fifth is a κ parameter
            lc_bas[5 + BAS_SLOTS*ib] = 0
            # Sixth is the env index address for exponents
            lc_bas[6 + BAS_SLOTS*ib] = off
            env[off+1:off+Ne] .= B.exp
            off += Ne
            # Seventh is the env index address for contraction coeff
            lc_bas[7 + BAS_SLOTS*ib] = off
            env[off+1:off+Nc] .= B.coef
            off += Nc
            # Eigth, nothing
            ib += 1
        end
    end
    empty!(memoize_cache(read_basisset))
    output("")
    return BasisSet(mol, basis_name, shells, natm, nbas, nshells, lc_atm, lc_bas, env)
end

"""
    Fermi.GaussianBasis.gto_norm(n::Signed, a::AbstractFloat)

Function that returns the normalization factor for a gaussian basis function rⁿ⋅exp(-a⋅r²)
"""
function gto_norm(n::Signed, a::AbstractFloat)
   # normalization factor of function rⁿ exp(-ar²)
    s = 2^(2n+3) * factorial(n+1) * (2a)^(n+1.5) / (factorial(2n+2) * √π)
    return √s
end

"""
    Fermi.GaussianBasis.normalize_basisfunction!(B::BasisFunction)

Modify the BasisFunction object by normalizing each primitive.
"""
function normalize_basisfunction!(B::BasisFunction)
    for i = eachindex(B.coef)
        B.coef[i] *= gto_norm(B.l, B.exp[i])
    end
end

function getindex(B::BasisSet, N::Int)
    A = B.molecule.atoms[N]
    return B.basis[A]
end

function getindex(B::BasisSet, i::Int, j::Int)
    A = B.molecule.atoms[i]
    return B.basis[A][j]
end

function string_repr(B::BasisFunction)
    # Generate Unicode symbol for sub number
    l_sub = Char(0x2080 + B.l)

    # Unicode for superscript is a bit messier, so gotta use control flow
    l_sup = B.l == 1 ? Char(0x00B9) :
            B.l in [2,3] ? Char(0x00B1 + B.l) :
            Char(0x2070 + B.l)

    nbas = 2*B.l + 1
    mvals = collect(-B.l:B.l)
    nprim = length(B.exp)

    # Reverse Dict(symbol=>num) to get Symbols from B.l
    Lsymbol = Dict(value => key for (key, value) in Fermi.GaussianBasis.AMDict)[B.l]
    out = "$(Lsymbol) shell with $nbas basis built from $nprim primitive gaussians\n\n"
    for m in mvals
        # Add sub minus sign (0x208B) if necessary
        m_sub = m < 0 ? Char(0x208B)*Char(0x2080 - m) : Char(0x2080 + m)
        out *= format("{:<4s} = ","χ$(l_sub)$m_sub")
        for i in eachindex(B.coef)

            if i > 1
                out *= B.coef[i] > 0 ? "\n     + " : "\n     - "
            end

            #out *= "$(abs(B.coef[i]))⋅Y$(l_sub)$m_sub"
            out *= format("{:>15.10f}⋅Y$(l_sub)$m_sub", abs(B.coef[i]))

            if B.l != 0
                out *= "⋅r$l_sup"
            end

            out *= "⋅exp(-$(B.exp[i])⋅r²)"
        end
        out *="\n\n"
    end
    return strip(out)
end

function string_repr(B::BasisSet)
    out  =  "$(B.basis_name) Basis Set\n"
    out *= "Number of shells: $(B.nshells)\n"
    out *= "Number of basis:  $(B.nbas)\n\n"

    l_to_symbol = Dict(
        0 => "s",
        1 => "p",
        2 => "d",
        3 => "f",
        4 => "g",
        5 => "h",
        6 => "i",
    )

    for A in B.molecule.atoms
        # Count how many times s,p,d appears for numbering
        count = zeros(Int16, 7)
        out *= "$(A.AtomicSymbol): "
        for bfs in B.basis[A]
            L = bfs.l
            count[L+1] += 1
            out *= "$(count[L+1])$(l_to_symbol[L]) "
        end
        out *="\n"
    end

    return strip(out)
end

# Pretty printing
function show(io::IO, ::MIME"text/plain", X::T) where T<:Union{BasisFunction, BasisSet}
    print(io, string_repr(X))
end

end #module


text:
# test functional calculation

"""
  Tests for functional computations.  The input file loads a 2 element mesh,
  square, -1 to 1,
  with a uniform flow (ICRho1E2U3).  The functional calculation should be exact
  in this case, so we can test against an analytical value.  The entire
  boundary is on BC 1
"""
function test_functionals()

  @testset "Testing functionals" begin

  # test all functional derivatives
  opts = Dict{String, Any}(
    "physics" => "Euler",
    "operator_type" => "SBPOmega",
    "dimensions" => 2,
    "run_type" => 5,
    "jac_method" => 2,
    "jac_type" => 2,
    "order" => 1,
    "IC_name" => "ICIsentropicVortex",
    "use_DG" => true,
    "volume_integral_type" => 2,
    "Volume_flux_name" => "IRFlux",
    "face_integral_type" => 2,
    "FaceElementIntegral_name" => "ESLFFaceIntegral",
    "Flux_name" => "IRFlux",
    "numBC" => 3,
    "BC1" => [0],
    "BC1_name" => "isentropicVortexBC",  # outlet
    "BC2" => [2],
    "BC2_name" => "isentropicVortexBC", # inlet
    "BC3" => [1, 3],
    "BC3_name" => "noPenetrationBC",  # was noPenetrationBC
    "aoa" => 0.0,
    "smb_name" => "SRCMESHES/vortex_3x3_.smb",
    "dmg_name" => ".null",
    "itermax" => 20,
    "res_abstol" => 1e-9,
    "res_reltol" => 1e-9,
    "do_postproc" => true,
    "exact_soln_func" => "ICIsentropicVortex",
    "force_solution_complex" => true,
    "force_mesh_complex" => true,
    "need_adjoint" => true,
    )

  # make second sets of 2D objects with Diagonal E (to test LPS)
  opts2 = copy(opts)
  opts2["operator_type"] = "SBPDiagonalE"
  delete!(opts2, "face_integral_type")
  delete!(opts2, "FaceElementIntegral_name")
  opts2["Flux_name"] = "IRSLFFlux"
  opts2["use_lps"] = true

  mesh, sbp, eqn, opts = solvePDE(opts)
  mesh2, sbp2, eqn2, opts2 = solvePDE(opts2)
  mesh3, sbp3, eqn3, opts3 = solvePDE("input_vals_jac3d.jl")

  testEntropyDissFunctional(mesh, sbp, eqn, opts)
  testEntropyDissFunctional2(mesh2, sbp2, eqn2, opts2)

  # test derivative of all functionals

  funcs_diage = ["negboundaryentropydiss", "entropydissipation", "negentropydissipation", "totalentropydissipation", "negtotalentropydissipation"]
  funcs_skip_zero = ["totalentropydissipation", "negtotalentropydissipation"]
  for funcname in keys(EulerEquationMod.FunctionalDict)
    println("testing functional ", funcname)
    obj = createFunctional(mesh, sbp, eqn, opts, funcname, [1, 3])
    if !(funcname in funcs_skip_zero) &&
        typeof(obj) <: EulerEquationMod.EntropyPenaltyFunctional ||
        funcname == "solutiondeviation"

      test_functional_zero(mesh, sbp, eqn, opts, obj)
    end

    test_functional_deriv_q(mesh, sbp, eqn, opts, obj, shock=true)

    if funcname in funcs_diage
      obj = createFunctional(mesh2, sbp2, eqn2, opts2, funcname, [1, 3])
      test_functional_deriv_q(mesh2, sbp2, eqn2, opts2, obj)
    end
  end

  func1 = createFunctional(mesh, sbp, eqn, opts, "massflow", [1, 3])
  func2 = createFunctional(mesh, sbp, eqn, opts, "lift", [1, 3])
  test_compositefunctional(mesh, sbp, eqn, opts, func1, func2)

  func1 = createFunctional(mesh, sbp, eqn, opts, "entropydissipation", [0])
  func2 = createFunctional(mesh, sbp, eqn, opts, "entropyjump", [0])
  test_compositefunctional(mesh, sbp, eqn, opts, func1, func2, test_revm=true)

#=
  obj = createFunctional(mesh3, sbp3, eqn3, opts3, "entropydissipation", [1])
  test_functional_zero(mesh3, sbp3, eqn3, opts3, obj)
  test_functional_deriv_q(mesh3, sbp3, eqn3, opts3, obj)
=#
  functional_revm_names = ["boundaryentropydiss",
                           "negboundaryentropydiss",
                           "entropydissipation",
                           "lpsdissipation",
                           "scdissipation",
                           "entropydissipation2",
                           "totalentropydissipation",
                           "negentropydissipation",
                           "neglpsdissipation",
                           "negscdissipation",
                           "negentropydissipation2",
                           "negtotalentropydissipation",
                           "lift", "liftCoefficient",
                           "drag", "dragCoefficient",
                           "solutiondeviation"]

  for funcname in functional_revm_names
    println("testing revm of functional ", funcname)

    obj = createFunctional(mesh, sbp, eqn, opts, funcname, [1, 3])
    test_functional_deriv_m(mesh, sbp, eqn, opts, obj)

    if funcname in funcs_diage
      obj = createFunctional(mesh2, sbp2, eqn2, opts2, funcname, [1, 3])
      test_functional_deriv_m(mesh2, sbp2, eqn2, opts2, obj)
    end

    obj3 = createFunctional(mesh3, sbp3, eqn3, opts3, funcname, [1])
    test_functional_deriv_m(mesh3, sbp3, eqn3, opts3, obj3, shock=true)
  end

  end  # end testset

  return nothing
end

"""
  Test the entropy dissipation functional.  This function is defined over all
  faces, rather than a boundary, so it is a bit special
"""
function testEntropyDissFunctional(mesh, sbp, eqn, _opts)

  opts = copy(_opts)
  opts2 = read_input_file("input_vals_channel.jl")
  opts2["Flux_name"] = "RoeFlux"
  opts2["use_DG"] = true
  opts2["solve"] = false

  mesh2, sbp2, eqn2, opts2 = solvePDE(opts2)

  obj = createFunctional(mesh2, sbp2, eqn2, opts2, "entropyflux", [1])
  val = evalFunctional(mesh2, sbp2, eqn2, opts2, obj)
  # the functional is u_i * U dot n_i, so for a constant field around a closed
  # curve it is zero
  @test isapprox(val, 0.0) atol=1e-13

  # compute the functional the regular way
  func = createFunctional(mesh, sbp, eqn, opts, "entropydissipation", Int[])
  J1 = evalFunctional(mesh, sbp, eqn, opts, func)

  # compute the functional using evalResidual
  opts["addVolumeIntegrals"] = false
  opts["addBoundaryIntegrals"] = false
  eqn.face_element_integral_func = EulerEquationMod.ELFPenaltyFaceIntegral(mesh, eqn)
  array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)
  w_vec = copy(eqn.q_vec)
  EulerEquationMod.convertToIR(mesh, sbp, eqn, opts, w_vec)

  J2 = dot(w_vec, eqn.res_vec)

  @test isapprox(J1, J2) atol=1e-13

  return nothing
end

function testEntropyDissFunctional2(mesh, sbp, eqn, opts)

  # test that the two methods of computing the entropy dissipation give the
  # same answer
  # This depends on the eqn object using the IRSLF flux with diagonal E
  # operators

  func1 = createFunctional(mesh, sbp, eqn, opts, "entropydissipation", [1])
  func2 = createFunctional(mesh, sbp, eqn, opts, "entropydissipation2", [2])

  val1 = evalFunctional(mesh, sbp, eqn, opts, func1)
  val2 = evalFunctional(mesh, sbp, eqn, opts, func2)
  println("first functional value = ", val1)
  println("second functional value = ", val2)

  @test abs(val1 - val2) < 1e-13

  return nothing
end

"""
  Test that a constant state -> functional value = 0
"""
function test_functional_zero(mesh, sbp, eqn, opts, func)

  icfunc = EulerEquationMod.ICDict["ICRho1E2U3"]
  icfunc(mesh, sbp, eqn, opts, eqn.q_vec)
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

  f = evalFunctional(mesh, sbp, eqn, opts, func)

  @test abs(f) < 1e-13

  return nothing
end

"""
  Test the derivative of a boundary functional wrt to q

  The eqn object must have been created with Tsol = Complex128 for this to work
"""
function test_functional_deriv_q(mesh, sbp, eqn, opts, func; shock=false)

  println("testing functional deriv q for functional ", typeof(func))
  h = 1e-20
  pert = Complex128(0, h)

  # use a spatially varying solution
  icfunc = EulerEquationMod.ICDict["ICRho1E2U3"]
  icfunc(mesh, sbp, eqn, opts, eqn.q_vec)
  eqn.q_vec .+= 0.01*rand(length(eqn.q_vec))
  if shock
    for i=1:div(mesh.numEl, 2)
      eqn.q[1, 1, i] += 1
    end
  end
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

  q_dot = rand(size(eqn.q))
  q_bar = zeros(eqn.q)

  evalFunctionalDeriv_q(mesh, sbp, eqn, opts, func, q_bar)
  val = sum(q_bar .* q_dot)

  eqn.q_vec .+= pert*vec(q_dot)
  f = evalFunctional(mesh, sbp, eqn, opts, func)
  val2 = imag(f)/h
  eqn.q_vec .-= pert*vec(q_dot)

  println("val1 = ", val)
  println("val2 = ", val2)
  @test abs(val - val2) < 1e-12

  return nothing
end

function test_functional_deriv_m(mesh, sbp, eqn, opts, func; shock=false)

  h = 1e-20
  pert = Complex128(0, h)

  # use a spatially varying solution
  if mesh.dim == 3
    icfunc = EulerEquationMod.ICDict["ICExp"]
  else
    icfunc = EulerEquationMod.ICDict["ICIsentropicVortex"]
  end
  icfunc(mesh, sbp, eqn, opts, eqn.q_vec)
  eqn.q_vec .+= 0.1*rand(length(eqn.q_vec))
  if shock
    for i=1:div(mesh.numEl, 2)
      eqn.q[1, 1, i] += 1
    end
  end
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

  zeroBarArrays(mesh)
  func_bar = rand_realpart(mesh.numDof)

  dxidx_dot       = rand_realpart(size(mesh.dxidx))
  jac_dot         = rand_realpart(size(mesh.jac))
  nrm_bndry_dot   = rand_realpart(size(mesh.nrm_bndry))
  nrm_face_dot    = rand_realpart(size(mesh.nrm_face_bar))
  coords_bndry_dot = rand_realpart(size(mesh.coords_bndry))

  mesh.dxidx        .+= pert*dxidx_dot
  mesh.jac          .+= pert*jac_dot
  mesh.nrm_bndry    .+= pert*nrm_bndry_dot
  mesh.nrm_face     .+= pert*nrm_face_dot
  mesh.coords_bndry .+= pert*coords_bndry_dot

  val = evalFunctional(mesh, sbp, eqn, opts, func)
  println("functional value = ", val)
  val = imag(val/h)

  mesh.dxidx        .-= pert*dxidx_dot
  mesh.jac          .-= pert*jac_dot
  mesh.nrm_bndry    .-= pert*nrm_bndry_dot
  mesh.nrm_face     .-= pert*nrm_face_dot
  mesh.coords_bndry .-= pert*coords_bndry_dot

  evalFunctionalDeriv_m(mesh, sbp, eqn, opts, func)

  val2 = sum(mesh.dxidx_bar .* dxidx_dot)              +
         sum(mesh.jac_bar .* jac_dot)                  +
         sum(mesh.nrm_bndry_bar .* nrm_bndry_dot)      +
         sum(mesh.nrm_face_bar .* nrm_face_dot)        +
         sum(mesh.coords_bndry_bar .* coords_bndry_dot)

#  println("val = ", real(val))
#  println("val2 = ", real(val2))
#  println("max dxidx_bar = ", maximum(abs.(mesh.dxidx_bar)))
#  println("max jac_bar = ", maximum(abs.(mesh.jac_bar)))
#  println("max nrm_bndry_bar = ", maximum(abs.(mesh.nrm_bndry_bar)))
#  println("max nrm_face_bar = ", maximum(abs.(mesh.nrm_face_bar)))
#  println("max coords_bndry_bar = ", maximum(abs.(mesh.coords_bndry_bar)))
  @test abs(val - val2) < 1e-12

  # test val_bar != 1
  zeroBarArrays(mesh)
  evalFunctionalDeriv_m(mesh, sbp, eqn, opts, func, 2)

  val3 = sum(mesh.dxidx_bar .* dxidx_dot)              +
         sum(mesh.jac_bar .* jac_dot)                  +
         sum(mesh.nrm_bndry_bar .* nrm_bndry_dot)      +
         sum(mesh.nrm_face_bar .* nrm_face_dot)        +
         sum(mesh.coords_bndry_bar .* coords_bndry_dot)

  @test abs(val3 - 2*val2) < 1e-12

  # test accumulation
  evalFunctionalDeriv_m(mesh, sbp, eqn, opts, func, 2)

  val4 = sum(mesh.dxidx_bar .* dxidx_dot)              +
         sum(mesh.jac_bar .* jac_dot)                  +
         sum(mesh.nrm_bndry_bar .* nrm_bndry_dot)      +
         sum(mesh.nrm_face_bar .* nrm_face_dot)        +
         sum(mesh.coords_bndry_bar .* coords_bndry_dot)

  @test abs(val4 - 2*val3) < 1e-12

  return nothing
end

function test_compositefunctional(mesh, sbp, eqn, opts,
                      func1::AbstractFunctional, func2::AbstractFunctional;
                      test_revm=false)
  println("testing CompositeFunctional")

  icfunc = EulerEquationMod.ICDict["ICRho1E2U3"]
  icfunc(mesh, sbp, eqn, opts, eqn.q_vec)
  eqn.q_vec .+= 0.1*rand_realpart(length(eqn.q_vec))
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

  func3 = CompositeFunctional(func1, func2)
  J1 = evalFunctional(mesh, sbp, eqn, opts, func1)
  J2 = evalFunctional(mesh, sbp, eqn, opts, func2)
  J3 = evalFunctional(mesh, sbp, eqn, opts, func3)

  @test abs(J3 - (J1 + J2)) < 1e-13

  dJdq1 = zeros(eqn.q)
  dJdq2 = zeros(eqn.q)
  dJdq3 = zeros(eqn.q)

  evalFunctionalDeriv_q(mesh, sbp, eqn, opts, func1, dJdq1)
  evalFunctionalDeriv_q(mesh, sbp, eqn, opts, func2, dJdq2)
  evalFunctionalDeriv_q(mesh, sbp, eqn, opts, func3, dJdq3)

  @test maximum(abs.(dJdq3 - (dJdq1 + dJdq2))) < 1e-13

  if test_revm
    zeroBarArrays(mesh)
    evalFunctionalDeriv_m(mesh, sbp, eqn, opts, func1)
    evalFunctionalDeriv_m(mesh, sbp, eqn, opts, func2)

    dxidx_bar = copy(mesh.dxidx_bar)
    jac_bar = copy(mesh.jac_bar)
    nrm_face_bar = copy(mesh.nrm_face_bar)
    nrm_bndry_bar = copy(mesh.nrm_bndry_bar)
    coords_bndry_bar = copy(mesh.coords_bndry_bar)

    zeroBarArrays(mesh)
    evalFunctionalDeriv_m(mesh, sbp, eqn, opts, func3)

    @test maximum(abs.(mesh.dxidx_bar - dxidx_bar)) < 1e-13
    @test maximum(abs.(mesh.jac_bar - jac_bar)) < 1e-13
    @test maximum(abs.(mesh.nrm_face_bar - nrm_face_bar)) < 1e-13
    @test maximum(abs.(mesh.nrm_bndry_bar - nrm_bndry_bar)) < 1e-13
    @test maximum(abs.(mesh.coords_bndry_bar - coords_bndry_bar)) < 1e-13
  end
end

add_func1!(EulerTests, test_functionals, [TAG_FUNCTIONAL, TAG_SHORTTEST])