text:
function compute_average_speed_of_sound(settings, z, c_0, T_0)

    n_intermediate = 11
    dz = z ./ (n_intermediate-1)
    c_bar = c_0
    for k in 1:(n_intermediate-1)
        T_im = T_0 .+ settings.dT .- (z .- k .* dz) .* (-0.0065)
        c_im = sqrt.(1.4 * 287. .* T_im)
        c_bar = (k) / (k + 1) .* c_bar .+ c_im ./ (k + 1)
    end

    return c_bar
end

function compute_euler_transformation(n_vcr_a_1, n_vcr_a_2, n_vcr_a_3, theta_B, phi_B, psi_B)

    cth  = cos.(deg2rad.(theta_B))
    sth  = sin.(deg2rad.(theta_B))
    cphi = cos.(deg2rad.(phi_B))
    sphi = sin.(deg2rad.(phi_B))
    cpsi = cos.(deg2rad.(psi_B))
    spsi = sin.(deg2rad.(psi_B))
    n_vcr_s_1 = cth .* cpsi .* n_vcr_a_1 .+ cth .* spsi .* n_vcr_a_2 - sth .* n_vcr_a_3
    n_vcr_s_2 = (-spsi .* cphi .+ sphi .* sth .* cpsi) .* n_vcr_a_1 .+ ( cphi .* cpsi .+ sphi .* sth .* spsi) .* n_vcr_a_2 .+ sphi .* cth .* n_vcr_a_3
    n_vcr_s_3 = (spsi .* sphi .+ cphi .* sth .* cpsi) .* n_vcr_a_1 .+ ( -sphi .* cpsi .+ cphi .* sth .* spsi) .* n_vcr_a_2 .+ cphi .* cth .* n_vcr_a_3

    return n_vcr_s_1, n_vcr_s_2, n_vcr_s_3
end

function geometry(settings, x_obs::Array{Float64, 1}, x, y, z, alpha, gamma, t_s, c_0, T_0)

    # Size of inputs
    n_t = size(x)[1]

    # Initialize vectors
    T = eltype(x)
    r = zeros(T, n_t)
    beta = zeros(T, n_t)
    theta = zeros(T, n_t)
    phi = zeros(T, n_t)
    c_bar = zeros(T, n_t)
    t_o = zeros(T, n_t)

    # Compute body angles (psi_B, theta_B, phi_B): angle of body w.r.t. horizontal
    theta_B = alpha .+ gamma
    phi_B = zeros(T, n_t)
    psi_B = zeros(T, n_t)

    # Compute the relative observer-aircraft position vector i.e. difference between observer and ac coordinate
    r_1 =  x_obs[1] .- x
    r_2 =  x_obs[2] .- y
    r_3 = -x_obs[3] .+ z

    # Normalize the distance vector
    r = sqrt.(r_1 .^2 .+ r_2 .^2 .+ r_3 .^2)
    n_vcr_a_1 = r_1 ./ r
    n_vcr_a_2 = r_2 ./ r
    n_vcr_a_3 = r_3 ./ r

    # Define elevation angle
    beta = asind.(n_vcr_a_3)

    # Transformation direction cosines (Euler angles) to the source coordinate system (i.e. take position of the aircraft into account)
    n_vcr_s_1, n_vcr_s_2, n_vcr_s_3 = compute_euler_transformation(n_vcr_a_1, n_vcr_a_2, n_vcr_a_3, theta_B, phi_B, psi_B)

    # Define polar directivity angle
    theta = acosd.(n_vcr_s_1)

    # Define azimuthal directivity angle
    phi = -atand.(n_vcr_s_2, n_vcr_s_3)

    # Average speed of sound
    c_bar = compute_average_speed_of_sound(settings, z, c_0, T_0)

    # Define observed time
    t_o = t_s .+ r ./ c_bar

    return r, beta, theta, phi, c_bar, t_o
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
const RadialOrbitalArray{T,N,B<:Basis} = Mul{<:Any,<:Tuple{<:BasisOrRestricted{B},<:AbstractArray{T,N}}}
const AdjointRadialOrbitalArray{T,N,B<:Basis} = Mul{<:Any,<:Tuple{<:Adjoint{<:Any,<:AbstractArray{T,N}},<:AdjointBasisOrRestricted{B}}}

"""
    RadialOrbital

A radial orbital is represented using basis coupled with a vector of
expansion coefficients with respect to that basis; the basis
implemented as an `AbstractQuasimatrix`.
"""
const RadialOrbital{T,B<:Basis} = RadialOrbitalArray{T,1,B}
const AdjointRadialOrbital{T,B<:Basis} = AdjointRadialOrbitalArray{T,1,B}

Base.axes(v::RadialOrbital, p::Int) = p == 1 ? axes(v)[1] : 1

"""
    RadialOrbitals

A collection of radial orbitals is instead represented using a
matrix of such expansion coefficients, where each matrix column
corresponds to a single radial orbital.
"""
const RadialOrbitals{T,B<:Basis} = RadialOrbitalArray{T,2,B}
const AdjointRadialOrbitals{T,B<:Basis} = AdjointRadialOrbitalArray{T,2,B}

const RadialOperator{T,B<:Basis,M<:AbstractMatrix{T}} =
    Mul{<:Any,<:Tuple{<:BasisOrRestricted{B},M,<:AdjointBasisOrRestricted{B}}}

function Base.show(io::IO, ro::RadialOperator{T,B,M}) where {T,B,M}
    write(io, "RadialOperator(")
    show(io, M)
    write(io, ")")
end

Base.iszero(op::RadialOperator) = iszero(op.args[2])
Base.iszero(op::MulQuasiArray{<:Any,2,<:Tuple{
    <:BasisOrRestricted,<:AbstractMatrix,<:AdjointBasisOrRestricted}}) =
        iszero(op.applied)

matrix(o::RadialOperator) = o.args[2]

function Base.zero(o::RadialOperator)
    A,B,C = o.args
    applied(*, A, Diagonal(zeros(eltype(B),size(B,1))), C)
end

function Base.:(+)(a::RadialOperator{T,B}, b::RadialOperator{T,B}) where {T,B}
    R = first(a.args)
    applied(*, R, matrix(a) + matrix(b), R')
end

ext: jl


text:
module LuaCutscenesLuaCutsceneTrigger

using ..Ahorn, Maple

@mapdef Trigger "luaCutscenes/luaCutsceneTrigger" LuaCutsceneTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight, onlyOnce::Bool=true, filename::String="", arguments::String="", unskippable::Bool=false)

const placements = Ahorn.PlacementDict(
    "Lua Cutscenes (Lua Cutscene)" => Ahorn.EntityPlacement(
        LuaCutsceneTrigger,
        "rectangle"
    )
)

end

max_stars_repo_name: rhelmot/LuaCutscenes

max_issues_repo_name: rhelmot/LuaCutscenes

max_forks_repo_name: rhelmot/LuaCutscenes


text:
export FDTDEnv

abstract Acoustics
abstract AcousticEnvironment <:Acoustics
abstract FDTDAcousticEnvironment <: AcousticEnvironment

immutable FDTDEnv <: FDTDAcousticEnvironment
	Fs::Float64        # sampling frequency
	c::Float64         # speed of sound
	X::Float64         # spatial step
	T::Float64         # temp. step
        scheme::FDTDScheme #FDTD scheme
end

"""
Create acoustic environment for FDTD simulation

## Usage

* 'FDTD(X::Float64, scheme::FDTDScheme)'
   * Create an acoustic env. for a given spatial step `X`
* 'FDTD(scheme::FDTDScheme, Fs::Float64)':
   * Create an acoustic env. for a given sampling frequency `Fs`
"""
function FDTDEnv(X::Float64, scheme::FDTDScheme; c = 343.)

	Fs = 343/(X*scheme.λ)   # Sampling Frequency
	return FDTDEnv(Fs,c,X,1/Fs,scheme)

end

function FDTDEnv(scheme::FDTDScheme, Fs::Float64; c = 343.)

	X = 343/(Fs*scheme.λ)   # Sampling Frequency
	return FDTDEnv(Fs,c,X,1/Fs,scheme)

end

function Base.show(io::IO, env::FDTDEnv)
	println(io, "FDTD Acoustic Env.")
	println(io, "Sampling freq. : $(round(env.Fs/1000,1))   kHz")
	println(io, "Spatial step   : $(round(env.X*100,2)) cm")
	println(io, "Speed of sound : $(round(env.c,1)) m/s")
	println(io, "FDTD Scheme    : $(env.scheme.name)")
end

max_stars_repo_name: nantonel/JuFdtd

max_issues_repo_name: nantonel/JuFdtd

max_forks_repo_name: nantonel/JuFdtd


text:
# Convert full BATSRUS AMR output to VTK.

using FortranFiles, WriteVTK, Glob, LightXML

export convertTECtoVTU, convertIDLtoVTK, readhead, readtree, getConnectivity
export Batl

# Named indexes of iTree_IA
const status_   =  Int8(1) # used, unused, to be refined, to be coarsened, etc.
const level_    =  Int8(2) # grid AMR level for this node
const proc_     =  Int8(3) # processor index where the block is stored for active nodes
const block_    =  Int8(4) # local block index for this node
const minLevel_ =  Int8(5) # minimum AMR level allowed for this node
const maxLevel_ =  Int8(6) # maximum AMR level allowed for this node
const coord0_   =  Int8(6) # equal to coord1_-1
const coord1_   =  Int8(7) # coordinate of node in 1st dimension
const coord2_   =  Int8(8) # coordinate of node in 2nd dimension
const coord3_   =  Int8(9) # coordinate of node in 3rd dimension
const coordLast_=  Int8(9) # coord0_ + MaxDim (?)
const parent_   = Int8(10) # parent_ must be equal to child0_
const child0_   = Int8(10) #
const child1_   = Int8(child0_ + 1)
#const childLast_= child0_ + nChild # This can be replaced with end

# Possible values for the status variable
const unset_   = Int8(-100) # index for unset values (that are otherwise larger)
const unused_     = Int8(-1) # unused block (not a leaf)
const refine_     = Int8(-2) # parent block to be refined
const noCoarsen_  = Int8(-3) # block not to be coarsened
const coarsen_    = Int8(-4) # child block to be coarsened
const used_       = Int8( 1) # currently used block (leaf)
const refineNew_  = Int8( 2) # child block to be refined
const refined_    = Int8( 3) # refined child block
const coarsenNew_ = Int8( 4) # parent block to be coarsened
const coarsened_  = Int8( 5) # coarsened parent block

# Deepest AMR level relative to root nodes (limited by 32 bit integers)
const maxLevel = 30

# The maximum integer coordinate for a given level below root nodes
const maxCoord_I = Int32[2^i for i in 0:maxLevel]

"BATSRUS output standalone header information."
struct Head
   # Number of cells per block in each direction.
   # These values are set by the Config.pl script.
   # Set 1 for ignored directions!
   nI::Int32
   nJ::Int32
   nK::Int32
   # Maximum number of ghost cells set by Config.pl script.
   # Valid values are 0,1,2,3,4,5
   nG::Int32
   # Refinement ratios in the 3 dimensions. Either 1 or 2.
   # The values are set by the Config.pl script.
   iRatio::Int32
   jRatio::Int32
   kRatio::Int32
   # Number of dimensions in which grid adaptation is done
   nDimAmr::Int32
   # Number of children per node
   nChild::Int32
   # Number of root blocks
   nRoot_D::Vector{Int32}
   # Grid lower and upper edges
   CoordMin_D::Vector{Float64}
   CoordMax_D::Vector{Float64}
   # Indexes of AMR dimensions
   iDimAmr_D::Vector{Int32}
   isPeriodic_D::Vector{Bool}
   # Plotting resolution
   dxPlot_D::Vector{Float64}
end

struct Batl
   head::Head
   iTree_IA::Array{Int32,2}
   iRatio_D::Vector{Int32} # Array of refinement ratios
   nDim::Int8
end

"""
	convertTECtoVTU(head, data, connectivity, filename="out")

Convert unstructured Tecplot data to VTK. Note that if using voxel type data
in VTK, the connectivity sequence is different from Tecplot.
Note that the 3D connectivity sequence in Tecplot is the same with the
`hexahedron` type in VTK, but different with the `voxel` type.
The 2D connectivity sequence is the same as the `quad` type, but different with
the `pixel` type.
For example, in 3D the index conversion is:
```
# PLT to VTK voxel index_ = [1 2 4 3 5 6 8 7]
for i = 1:2
   connectivity = swaprows!(connectivity, 4*i-1, 4*i)
end
```
"""
function convertTECtoVTU(head, data, connectivity, filename="out")

   nVar = length(head.variables)
   points = @view data[1:head.nDim,:]
   cells = Vector{MeshCell{VTKCellType,Array{Int32,1}}}(undef,head.nCell)

   if head.nDim == 3
      @inbounds for i = 1:head.nCell
         cells[i] = MeshCell(VTKCellTypes.VTK_HEXAHEDRON, connectivity[:,i])
      end
   elseif head.nDim == 2
      @inbounds for i = 1:head.nCell
         cells[i] = MeshCell(VTKCellTypes.VTK_QUAD, connectivity[:,i])
      end
   end

   vtkfile = vtk_grid(filename, points, cells)

   for ivar = head.nDim+1:nVar
      if endswith(head.variables[ivar],"_x") # vector
         if head.nDim == 3
            var1 = @view data[ivar,:]
            var2 = @view data[ivar+1,:]
            var3 = @view data[ivar+2,:]
            namevar = replace(head.variables[ivar], "_x"=>"")
            vtk_point_data(vtkfile, (var1, var2, var3), namevar)
         elseif head.nDim == 2
            var1 = @view data[ivar,:]
            var2 = @view data[ivar+1,:]
            namevar = replace(head.variables[ivar], "_x"=>"")
            vtk_point_data(vtkfile, (var1, var2), namevar)
         end
      elseif endswith(head.variables[ivar],r"_y|_z")
         continue
      else
         var = @view data[ivar,:]
         vtk_point_data(vtkfile, var, head.variables[ivar])
      end
   end

   # Add meta data from Tecplot AUXDATA
   for i in 1:length(head.auxdata)
      vtkfile[head.auxdataname[i],VTKFieldData()] = head.auxdata[i]
   end

   outfiles = vtk_save(vtkfile)
end

"""
	convertIDLtoVTK(filename; dir=".", gridType=1, verbose=false)

Convert 3D BATSRUS *.out to VTK. If `gridType==1`, it converts to the
rectilinear grid; if `gridType==2`, it converts to the structured grid.
If `filename` does not end with "out", it tries to find the ".info" and ".tree"
file with the same name tag and generates 3D unstructured VTU file.
"""
function convertIDLtoVTK(filename::AbstractString; dir=".", gridType=1,
   verbose=false)

   if endswith(filename, ".out")
      data = readdata(filename, dir=dir)

      nVar = length(data.head.wnames)

      outname = filename[1:end-4]

      if gridType == 1 # rectilinear grid
         x = @view data.x[:,1,1,1]
         y = @view data.x[1,:,1,2]
         z = @view data.x[1,1,:,3]

         outfiles = vtk_grid(outname, x,y,z) do vtk
            for ivar = 1:nVar
               if data.head.wnames[ivar][end] == 'x' # vector
                  var1 = @view data.w[:,:,:,ivar]
                  var2 = @view data.w[:,:,:,ivar+1]
                  var3 = @view data.w[:,:,:,ivar+2]
                  namevar = data.head.wnames[ivar][1:end-1]
                  vtk_point_data(vtk, (var1, var2, var3), namevar)
               elseif data.head.wnames[ivar][end] in ('y','z')
                  continue
               else
                  var = @view data.w[:,:,:,ivar]
                  vtk_point_data(vtk, var, data.head.wnames[ivar])
               end
            end
         end
      elseif gridType == 2 # structured grid
         xyz = permutedims(data.x, [4,1,2,3])

         outfiles = vtk_grid(outname, xyz) do vtk
            for ivar = 1:nVar
               if data.head.wnames[ivar][end] == 'x' # vector
                  var1 = @view data.w[:,:,:,ivar]
                  var2 = @view data.w[:,:,:,ivar+1]
                  var3 = @view data.w[:,:,:,ivar+2]
                  namevar = data.head.wnames[ivar][1:end-1]
                  vtk_point_data(vtk, (var1, var2, var3), namevar)
               elseif data.head.wnames[ivar][end] in ('y','z')
                  continue
               else
                  var = @view data.w[:,:,:,ivar]
                  vtk_point_data(vtk, var, data.head.wnames[ivar])
               end
            end
         end
      elseif gridType == 3 # unstructured grid
         @error "No tree information for conversion!"
      end

   else
      # info, tree, and out files
      data = readdata(filename*".out", dir=dir)
      batl = Batl(readhead(filename*".info"), readtree(filename)...)
      connectivity = getConnectivity(batl)

      outname = filename

      nDim = batl.nDim
      nVar = length(data.head.wnames)
      nCell = size(connectivity, 2)
      if nDim == 3
         if batl.head.dxPlot_D[1] ≥ 0.0
            @error "Why are there duplicate points?"
         else
            points = data.x[:,1,1,:]'
         end
      elseif nDim == 2
         if batl.head.dxPlot_D[1] ≥ 0.0 # points are not sorted in postproc.f90
            @error "Why are there duplicate points? Ask!"
            points = data.x[:,1,:]'
         else # points are sorted in postproc.f90!
            @error "point original order cannot be retrieved!"
         end
      end
      cells = Vector{MeshCell{VTKCellType,Array{Int32,1}}}(undef,nCell)

      if nDim == 3
         @inbounds for i = 1:nCell
            cells[i] = MeshCell(VTKCellTypes.VTK_HEXAHEDRON, connectivity[:,i])
         end
      elseif nDim == 2
         @inbounds for i = 1:nCell
            cells[i] = MeshCell(VTKCellTypes.VTK_QUAD, connectivity[:,i])
         end
      end

      vtkfile = vtk_grid(filename, points, cells)

      for ivar = 1:nVar
         if endswith(data.head.wnames[ivar], "x") # vector
            if nDim == 3
               var1 = @view data.w[:,1,1,ivar]
               var2 = @view data.w[:,1,1,ivar+1]
               var3 = @view data.w[:,1,1,ivar+2]
               var = (var1, var2, var3)
               vtkfile[data.head.wnames[ivar][1:end-1], VTKPointData()] = var
            elseif nDim == 2 # not sure how VTK handles 2D vector!
               var1 = @view data.w[:,1,ivar]
               var2 = @view data.w[:,1,ivar+1]
               vtkfile[data.head.wnames[ivar][1:end-1], VTKPointData()] = (var1, var2)
            end
         elseif endswith(data.head.wnames[ivar], r"y|z")
            continue
         else
            if nDim == 3
               var = @view data.w[:,1,1,ivar]
            elseif nDim == 2
               var = @view data.w[:,1,ivar]
            end
            vtkfile[data.head.wnames[ivar], VTKPointData()] = var
         end
      end

      outfiles = vtk_save(vtkfile)
   end

   verbose && @info "$(filename) finished conversion."
   return
end

"Return matrix X with swapped rows i and j."
function swaprows!(X, i, j)
   m, n = size(X)
   if (1 ≤ i ≤ n) && (1 ≤ j ≤ n)
      @inbounds @simd for k = 1:n
         X[i,k],X[j,k] = X[j,k],X[i,k]
      end
      return X
   else
      throw(BoundsError())
   end
end

"Return header from info file. Currently only designed for 2D and 3D."
function readhead(filehead)

   nDim = 3
   nI, nJ, nK = 1, 1, 1
   nRoot_D = Int32[1,1,1]
   # Make sure that the thickness is unity in the ignored dimensions.
   CoordMin_D = fill(-0.5, 3)
   CoordMax_D = fill(0.5, 3)
   isPeriodic_D = fill(false, 3)
   dxPlot_D = fill(0.0,3)

   open(filehead) do f
      while !eof(f)
         ln = readline(f)
         if startswith(ln, "#NDIM")
            nDim = parse(Int32, split(readline(f))[1])
         elseif startswith(ln, "#GRIDBLOCKSIZE")
            nI = parse(Int32, split(readline(f))[1])
            nJ = parse(Int32, split(readline(f))[1])
            if nDim == 2
               nK = Int32(1)
            elseif nDim == 3
               nK = parse(Int32, split(readline(f))[1])
            end
         elseif startswith(ln, "#ROOTBLOCK")
            nRoot_D[1] = parse(Int32, split(readline(f))[1])
            nRoot_D[2] = parse(Int32, split(readline(f))[1])
            if nDim == 3
               nRoot_D[3] = parse(Int32, split(readline(f))[1])
            end
         elseif startswith(ln, "#GRIDGEOMETRYLIMIT")
            typeGeometry = split(readline(f))[1]
            CoordMin_D[1] = parse(Float64, split(readline(f))[1])
            CoordMax_D[1] = parse(Float64, split(readline(f))[1])
            CoordMin_D[2] = parse(Float64, split(readline(f))[1])
            CoordMax_D[2] = parse(Float64, split(readline(f))[1])
            if nDim == 3
               CoordMin_D[3] = parse(Float64, split(readline(f))[1])
               CoordMax_D[3] = parse(Float64, split(readline(f))[1])
            end
         elseif startswith(ln, "#PERIODIC")
            str = split(readline(f))[1]
            isPeriodic_D[1] = str == "T" ? true : false
            str = split(readline(f))[1]
            isPeriodic_D[2] = str == "T" ? true : false
            if nDim == 3
               str = split(readline(f))[1]
               isPeriodic_D[3] = str == "T" ? true : false
            end
         elseif startswith(ln, "#PLOTRESOLUTION")
            dxPlot_D[1] = parse(Float64, split(readline(f))[1])
            dxPlot_D[2] = parse(Float64, split(readline(f))[1])
            if nDim == 3
               dxPlot_D[3] = parse(Float64, split(readline(f))[1])
            end
         end
      end
   end

   nG = 2 # most outputs are 2nd order

   iRatio, jRatio, kRatio = min(2, nI), min(2, nJ), min(2, nK)

   nDimAmr = iRatio + jRatio + kRatio - 3

   nChild = 2^nDimAmr

   # Get size of domain (in generalized coordinates)
   DomainSize_D = CoordMax_D - CoordMin_D

   # Indexes of AMR dimensions.
   # The magic formulas should be correct from 1 to nDimAmr.
   iDimAmrTmp_D = Int32[ 1 + (2-iRatio)*(3-jRatio), 6-iRatio-jRatio, 3 ]
   iDimAmr_D = iDimAmrTmp_D[1:nDimAmr]

   Head(nI, nJ, nK, nG, iRatio, jRatio, kRatio, nDimAmr, nChild, nRoot_D,
      CoordMin_D, CoordMax_D, iDimAmr_D, isPeriodic_D, dxPlot_D)
end

"Return BATL tree structure."
function readtree(filetag)

   iRatio_D = Int32[1,1,1]
   nRoot_D = Int32[1,1,1]

   # Loading AMR tree
   f = FortranFile(filetag*".tree")

   nDim, nInfo, nNode = read(f, Int32, Int32, Int32)
   iRatio_D[1:nDim] = read(f, (Int32,nDim)) # Array of refinement ratios
   nRoot_D[1:nDim] = read(f, (Int32,nDim)) # The number of root nodes in all dimension
   iTree_IA = read(f, (Int32,(nInfo,nNode)))

   close(f)

   return iTree_IA, iRatio_D, nDim
end

"""
   find_grid_block(batl, xyz_D)

Return processor local block index that contains a point. Input location should
be given in Cartesian coordinates.
"""
function find_grid_block(batl::Batl, xyz_D)

   CoordMin_D = batl.head.CoordMin_D
   CoordMax_D = batl.head.CoordMax_D

   Coord_D = xyz_D

   # Calculate normalized coordinates for tree search
   CoordTree_D = (Coord_D - CoordMin_D) ./ (CoordMax_D - CoordMin_D)

   if any(CoordTree_D .< 0.0) || any(CoordTree_D .> 1.0)
      iBlock = unset_
      return iBlock
   end

   # Find node containing the point
   iNode = find_tree_node(batl, CoordTree_D)

   # Check if point was found
   if iNode > 0
      # Convert to block and processor indexes
      iBlock = batl.iTree_IA[block_,iNode]
   else
      iBlock = unset_
   end

   return iBlock
end

"""
   find_tree_node(batl, Coord_D)

Find the node that contains a point. The point coordinates should be given in
generalized coordinates normalized by the domain size.
"""
function find_tree_node(batl::Batl, Coord_D)

   nRoot_D = batl.head.nRoot_D
   nDimAmr = batl.head.nDimAmr
   iDimAmr_D = batl.head.iDimAmr_D
   iTree_IA = batl.iTree_IA

   # Scale coordinates so that 1 ≤ Coord_D ≤ nRoot_D+1
   Coord_D = @. 1.0 + nRoot_D*max(0.0, min(1.0, Coord_D))

   # Get root node index
   iRoot_D = min.(floor.(Int32, Coord_D), nRoot_D)

   # Root node indexes are ordered
   iNode = iRoot_D[1] + nRoot_D[1]*((iRoot_D[2]-1) + nRoot_D[2]*(iRoot_D[3]-1))

   if iTree_IA[status_,iNode] == used_
      return iNode
   end

   nLevelMax, iNodeMorton_I = order_tree(batl)

   # Get normalized coordinates within root node and scale it up
   # to the largest resolution: 0 <= iCoord_D <= maxCoord_I(nLevelMax)-1
   iCoord_D = min.(maxCoord_I[nLevelMax+1] - 1,
      floor.(Int32, (Coord_D[iDimAmr_D] - iRoot_D[iDimAmr_D])*maxCoord_I[nLevelMax+1]))

   # Go down the tree using bit information
   for iLevel = nLevelMax-1:-1:0
      # Get the binary bits based on the coordinates
      iBit_D = ibits.(iCoord_D, iLevel, 1)

      # Construct child index as iChild = Sum Bit_i*2^i
      iChild = sum(iBit_D.*maxCoord_I[1:nDimAmr]) + child1_
      iNode  = iTree_IA[iChild,iNode]

      if iTree_IA[status_,iNode] == used_
         return iNode
      end
   end

   # Did not find the point so set iNode as unset
   iNode = unset_

   return iNode
end

"Logical shifts as the Fortran instrinsic function."
function ibits(i, pos, len)
   # Treat it as 32 bits integer as in BATL.
   ds = digits(i, base=2, pad=32)[pos+1:pos+len]
   s = zero(eltype(i))
   for val in ds
      s = s * 2 + val
   end
   return s
end

"""
   order_tree(batl)

Return maximum AMR level in the used block and the Morton curve order.
Set iNodeMorton_I indirect index arrays according to
1. root node order
2. Morton ordering for each root node
"""
function order_tree(batl::Batl)

   nRoot_D = batl.head.nRoot_D
   iTree_IA = batl.iTree_IA

   nNode = prod(nRoot_D)
   iNode = 0
   iMorton = 0

   # Get used node counts
   nNodeUsed = count(==(used_), iTree_IA[status_,:])

   iNodeMorton_I = fill(Int32(unset_), nNodeUsed)

   for kRoot = 1:nRoot_D[3]
      for jRoot = 1:nRoot_D[2]
         for iRoot = 1:nRoot_D[1]
            # Root nodes are the first ones
            iNode += 1
            # All root nodes are handled as if they were first child
            iMorton = order_children!(batl, iNode, iMorton, iNodeMorton_I)
         end
      end
   end

   nNodeUsed = iMorton

   # Set min and max refinement levels
   nLevelMin = maxLevel
   nLevelMax = 0
   for iMorton = 1:nNodeUsed
      iNode = iNodeMorton_I[iMorton]
      iLevel = iTree_IA[level_,iNode]
      nLevelMin = min(iLevel, nLevelMin)
      nLevelMax = max(iLevel, nLevelMax)
   end

   return nLevelMax, iNodeMorton_I
end

"""
   order_children!(iNode)

Recursively apply Morton ordering for nodes below a root block.
Store result into iNodeMorton_I and iMortonNode_A using the iMorton index.
"""
function order_children!(batl::Batl, iNode, iMorton, iNodeMorton_I)

   iTree_IA = batl.iTree_IA

   if iTree_IA[status_, iNode] ≥ used_
      iMorton += 1
      iNodeMorton_I[iMorton] = iNode
   else
      for iChild = child1_:size(iTree_IA,1)
         iMorton = order_children!(batl, iTree_IA[iChild, iNode], iMorton, iNodeMorton_I)
      end
   end

   return iMorton
end

"Find neighbors for any node in the tree. Only for Cartesian now."
function find_neighbor_for_anynode(batl::Batl, iNode)

   nDim = batl.nDim
   nRoot_D = batl.head.nRoot_D
   iTree_IA = batl.iTree_IA
   iRatio_D = batl.iRatio_D

   IsLatitudeAxis = false
   IsSphericalAxis = false
   IsCylindricalAxis = false

   # Periodicity check is necessary for recovering the boundary conditions,
   # but it makes postprocessing harder. An easy solution is to turn it off.
   #isPeriodic_D = batl.head.isPeriodic_D
   isPeriodic_D = [false, false, false]

   DiLevelNei_III = fill(Int8(0), (3,3,3))
   iNodeNei_III = fill(Int32(-1), (4,4,4))

   # Get AMR level of the node
   iLevel = iTree_IA[level_,iNode]

   # Calculate scaling factor from integer index to 0<x,y,z<1 real coordinates
   Scale_D = 1.0 ./ nRoot_D
   for i = 1:3
      if iRatio_D[i] == 2
         Scale_D[i] /= maxCoord_I[iLevel+1]
      end
   end

   # Fill in self-referring info
   iNodeNei_III[2:3,2:3,2:3] .= Int32(iNode)
   DiLevelNei_III[2,2,2]     = 0

   # Loop through neighbors
   for k = 0:3
      Dk = round(Int8, (k - 1.5)/1.5)
      if nDim < 3
         if k != 1 continue end
         z = 0.3
      else
         z = (iTree_IA[coord3_, iNode] + 0.4*k - 1.1)*Scale_D[3]
         if z > 1.0 || z < 0.0
            if isPeriodic_D[3]
               z = mod(z, 1.0)
            elseif !IsLatitudeAxis
               iNodeNei_III[:,:,k+1] .= unset_
               DiLevelNei_III[:,:,Dk+2] .= unset_
               continue
            end
         end
      end
      # store z for spherical axis
      z0 = z
      for j = 0:3
         z = z0
         Dj = round(Int8, (j - 1.5)/1.5)
         if nDim < 2
            if j != 1 continue end
            y = 0.3
         else
            y = (iTree_IA[coord2_, iNode] + 0.4*j - 1.1)*Scale_D[2]
            if y > 1.0 || y < 0.0
               if isPeriodic_D[2]
                  y = mod(y, 1.0)
               elseif IsSphericalAxis
                  # Push back theta and go around half way in phi
                  y = max(0.0, min(1.0, y))
                  z = mod(z0+0.5, 1.0)
               else
                  iNodeNei_III[:,j+1,k+1] .= unset_
                  DiLevelNei_III[:,Dj+2,Dk+2] .= unset_
                  continue
               end
            end
            if z0 > 1.0 || z0 < 0.0
               # Push back latitude and go around half way in longitude
               z = max(0.0, min(1.0, z0))
               y = mod(y+0.5, 1.0)
            end
         end
         # store y for cylindrical axis case
         y0 = y
         for i = 0:3
            # Exclude inner points
            if 0<i<3 && 0<j<3 && 0<k<3 continue end

            Di = round(Int8, (i - 1.5)/1.5)

            # If neighbor is not finer, fill in the i=2 or j=2 or k=2 elements
            if DiLevelNei_III[Di+2,Dj+2,Dk+2] ≥ 0
               if i == 2
                  iNodeNei_III[i+1,j+1,k+1] = iNodeNei_III[2,j+1,k+1]
                  continue
               end
               if j == 2
                  iNodeNei_III[i+1,j+1,k+1] = iNodeNei_III[i+1,2,k+1]
                  continue
               end
               if k == 2
                  iNodeNei_III[i+1,j+1,k+1] = iNodeNei_III[i+1,j+1,2]
                  continue
               end
            end

            x = (iTree_IA[coord1_, iNode] + 0.4*i - 1.1)*Scale_D[1]
            y = y0
            if x > 1.0 || x < 0.0
               if isPeriodic_D[1]
                  x = mod(x, 1.0)
               elseif IsCylindricalAxis && x < 0.0
                  # Push back radius and go around half way in phi direction
                  x = 0.0
                  y = mod(y0+0.5, 1.0)
               else
                  iNodeNei_III[i+1,j+1,k+1] = unset_
                  DiLevelNei_III[Di+2,Dj+2,Dk+2] = unset_
                  continue
               end
            end

            jNode = find_tree_node(batl, [x, y, z])

            iNodeNei_III[i+1,j+1,k+1] = jNode
            DiLevelNei_III[Di+2,Dj+2,Dk+2] = iLevel - iTree_IA[level_, jNode]
         end
      end
   end

   return iNodeNei_III, DiLevelNei_III
end

"Return global block index for the node."
function nodeToGlobalBlock(batl::Batl, iNode::Int32, nBlock_P)

   iTree_IA = batl.iTree_IA

   iProc = iTree_IA[proc_, iNode]

   localNodes_B = findall(x->x==iProc, iTree_IA[proc_,:])
   localBlocks_B = iTree_IA[block_,localNodes_B]
   myRank = findfirst(x->x==iNode, localNodes_B)

   localBlocksSorted_B = sort(localBlocks_B)
   localOrder = findfirst(x->x==localBlocks_B[myRank], localBlocksSorted_B)

   globalBlock = localOrder + nBlock_P[iProc+1]

   return globalBlock
end

"Get cell connectivity list."
function getConnectivity(batl::Batl)

   iTree_IA = batl.iTree_IA
   nI, nJ, nK = batl.head.nI, batl.head.nJ, batl.head.nK
   nDim = batl.nDim
   if nDim == 3
      nConn = 8
   elseif nDim == 2
      nConn = 4
      @error "2D not working currently!"
   elseif nDim == 1
      @error "Detected 1D data, connectivity not implemented!"
   end

   iCell_G = Array{Int32,3}(undef, nI+2, nJ+2, nK+2)

   # It would be difficult to calculate this beforehand.
   # The current implementation is to compute this in two rounds, where the
   # first round only do the counting.
   nElem = 0
   nBlockBefore = 0

   # Count accumulated number of blocks for MPI ranks in ascending order
   # Local block indexes may have gaps in between!
   nProc = maximum(iTree_IA[proc_,:]) + 1
   nBlock_P = fill(Int32(0), nProc)
   if nProc > 1
      for iProc = 1:nProc-1
         nBlock_P[iProc+1] = nBlock_P[iProc] + count(==(iProc-1), iTree_IA[proc_,:])
      end
   end

   # Pre-allocate just to let Julia know this variable.
   connectivity = Array{Int32,2}(undef, nConn, nElem)

   for iRound = 1:2
      if iRound == 2
         connectivity = Array{Int32,2}(undef, nConn, nElem)
         iElem = 0
         nBlockBefore = 0 # Reset
      end
      for iProc = 0:nProc-1
         localNodes_B = findall(x->x==iProc, iTree_IA[proc_,:])
         # It needs to be sorted for filling the correct iCell_G indexes!
         localBlocks_B = iTree_IA[block_,localNodes_B]
         seq = sortperm(localBlocks_B)
         localNodes_B = localNodes_B[seq] # sorted according to local block index

         nBlock = length(localNodes_B) # number of blocks on this processor

         for iBlock = 1:nBlock
            # initial cell index
            iCell = nI*nJ*nK*(nBlockBefore + iBlock - 1)

            @inbounds for k = 1:nK, j = 1:nJ, i = 1:nI
               iCell += 1
               iCell_G[i+1,j+1,k+1] = iCell
            end

            iNodeNei_III, DiLevelNei_III = find_neighbor_for_anynode(batl, localNodes_B[iBlock]) # rename!

            fillCellNeighbors!(batl, iCell_G, DiLevelNei_III, iNodeNei_III, nBlock_P)

            # Check who is in charge of writing following the prescribed rules.
            if DiLevelNei_III[3,3,2] < 0
               iCell_G[end,end,:] .= 0
            end
            if DiLevelNei_III[3,2,3] < 0
               iCell_G[end,:,end] .= 0
            end
            if DiLevelNei_III[2,3,3] < 0
               iCell_G[:,end,end] .= 0
            end
            if nDim == 3 && DiLevelNei_III[3,3,3] < 0
               iCell_G[end,end,end] = 0
            end

            iMin = DiLevelNei_III[1,2,2] == 1 ? 0 : 1
            jMin = DiLevelNei_III[2,1,2] == 1 ? 0 : 1
            kMin = DiLevelNei_III[2,2,1] == 1 ? 0 : 1

            iMax = DiLevelNei_III[3,2,2] < 0 ? nI-1 : nI
            jMax = DiLevelNei_III[2,3,2] < 0 ? nJ-1 : nJ
            kMax = DiLevelNei_III[2,2,3] < 0 ? nK-1 : nK

            if nDim == 3
               for k = kMin:kMax, j = jMin:jMax, i = iMin:iMax
                  # Skip bricks that are not fully inside/usable
                  if any(iCell_G[i+1:i+2,j+1:j+2,k+1:k+2] .== 0)
                     continue
                  end
                  iRound == 1 ? nElem += 1 : iElem += 1
                  if iRound == 2
                     connectivity[:,iElem] = [
                        iCell_G[i+1,j+1,k+1],
                        iCell_G[i+2,j+1,k+1],
                        iCell_G[i+2,j+2,k+1],
                        iCell_G[i+1,j+2,k+1],
                        iCell_G[i+1,j+1,k+2],
                        iCell_G[i+2,j+1,k+2],
                        iCell_G[i+2,j+2,k+2],
                        iCell_G[i+1,j+2,k+2]]
                  end
               end
            elseif nDim == 2
               for j = jMin:jMax, i = iMin:iMax
                  connectivity = hcat(connectivity,
                     [
                     iCell_G[i+1,j+1,2],
                     iCell_G[i+2,j+1,2],
                     iCell_G[i+2,j+2,2],
                     iCell_G[i+1,j+2,2]]
                  )
               end
            end
         end
         nBlockBefore += nBlock # for next processor in ascending order
      end
   end

   # The current writing process will generate duplicate connectivity rows at
   # edge cells between coarse and fine blocks. WriteTecplot module function
   # has the same issue.
   #return unique(connectivity, dims=2)
   return connectivity
end

"Return sibling index (1-8) for the given block node matrix."
function getSibling(iNodeNei_III, iTree_IA)
   iMyNode = iNodeNei_III[2,2,2]
   iMotherNode = iTree_IA[parent_,iMyNode]
   motherNode = iTree_IA[:,iMotherNode]
   iSibling = findfirst(x->x==iMyNode, motherNode[child1_:end])
end

"""
    fillCellNeighbors!(batl, iCell_G, DiLevelNei_III, iNodeNei_III, nBlock_P)

Fill neighbor cell indexes for the given block. The faces, edges, and vertices
are ordered from left (-) to right (+) in x-y-z sequentially.
Vertices:        Edges: (10,11 ignored)
   7 ----- 8        . --4-- .
 - .     - .      7 .     8 .
5 ----- 6  .     . --3-- .  12
.  .    .  .     .  .    .  .
.  3 ----- 4     9  . --2-- .
. -     . -      . 5     . 6
1 ----- 2        . --1-- .
Only tested for 3D.
"""
function fillCellNeighbors!(batl, iCell_G, DiLevelNei_III, iNodeNei_III, nBlock_P)

   iTree_IA = batl.iTree_IA
   nI, nJ, nK = batl.head.nI, batl.head.nJ, batl.head.nK

   nIJ = nI*nJ
   nIJK = nI*nJ*nK

   ## Faces

   # -x face
   if DiLevelNei_III[1,2,2] == 1
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,2,2], nBlock_P)

      @inbounds for k = 1:nK, j = 1:nJ
         iSibling = getSibling(iNodeNei_III, iTree_IA)

         if iSibling == 1
            iCell_G[1,j+1,k+1] = nIJK*(neiBlock-1) +
               nI*(1 + (j-1)÷2) + nIJ*((k-1)÷2)
         elseif iSibling == 3
            iCell_G[1,j+1,k+1] = nIJK*(neiBlock-1) +
               nIJ/2 + nI*(1 + (j-1)÷2) + nIJ*((k-1)÷2)
         elseif iSibling == 5
            iCell_G[1,j+1,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nI*(1 + (j-1)÷2) + nIJ*((k-1)÷2)
         elseif iSibling == 7
            iCell_G[1,j+1,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nIJ/2 + nI*(1 + (j-1)÷2) + nIJ*((k-1)÷2)
         end
      end
   end

   # +x face
   if DiLevelNei_III[3,2,2] == 0
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,3,3], nBlock_P)

      @inbounds for k = 1:nK, j = 1:nJ
         iCell_G[end,j+1,k+1] = nIJK*(neiBlock-1) + 1 + nI*(j-1) + nIJ*(k-1)
      end
   elseif DiLevelNei_III[3,2,2] == 1
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,3,3], nBlock_P)

      @inbounds for k = 1:nK, j = 1:nJ
         iSibling = getSibling(iNodeNei_III, iTree_IA)

         if iSibling == 2
            iCell_G[end,j+1,k+1] = nIJK*(neiBlock-1) +
               1 + nI*((j-1)÷2) + nIJ*((k-1)÷2)
         elseif iSibling == 4
            iCell_G[end,j+1,k+1] = nIJK*(neiBlock-1) +
               nIJ/2 + 1 + nI*((j-1)÷2) + nIJ*((k-1)÷2)
         elseif iSibling == 6
            iCell_G[end,j+1,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + 1 + nI*((j-1)÷2) + nIJ*((k-1)÷2)
         elseif iSibling == 8
            iCell_G[end,j+1,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nIJ/2 + 1 + nI*((j-1)÷2) + nIJ*((k-1)÷2)
         end
      end
   end

   # -y face
   if DiLevelNei_III[2,1,2] == 1
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[2,1,2], nBlock_P)

      @inbounds for k = 1:nK, i = 1:nI
         iSibling = getSibling(iNodeNei_III, iTree_IA)

         if iSibling == 1
            iCell_G[i+1,1,k+1] = nIJK*(neiBlock-1) +
               nI*(nJ-1) + 1 + (i-1)÷2 + nIJ*((k-1)÷2)
         elseif iSibling == 2
            iCell_G[i+1,1,k+1] = nIJK*(neiBlock-1) +
               nI*(nJ-1) + nI/2 + 1 + (i-1)÷2 + nIJ*((k-1)÷2)
         elseif iSibling == 5
            iCell_G[i+1,1,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nI*(nJ-1) + 1 + (i-1)÷2 + nIJ*((k-1)÷2)
         elseif iSibling == 6
            iCell_G[i+1,1,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nI*(nJ-1) + nI/2 + 1 + (i-1)÷2 + nIJ*((k-1)÷2)
         end
      end
   end

   # +y face
   if DiLevelNei_III[2,3,2] == 0
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[3,4,3], nBlock_P)

      @inbounds for k = 1:nK, i = 1:nI
         iCell_G[i+1,end,k+1] = nIJK*(neiBlock-1) + i + nIJ*(k-1)
      end
   elseif DiLevelNei_III[2,3,2] == 1
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[3,4,3], nBlock_P)

      @inbounds for k = 1:nK, i = 1:nI
         iSibling = getSibling(iNodeNei_III, iTree_IA)

         if iSibling == 3
            iCell_G[i+1,end,k+1] = nIJK*(neiBlock-1) +
               1 + (i-1)÷2 + nIJ*((k-1)÷2)
         elseif iSibling == 4
            iCell_G[i+1,end,k+1] = nIJK*(neiBlock-1) +
               nI/2 + 1 + (i-1)÷2 + nIJ*((k-1)÷2)
         elseif iSibling == 7
            iCell_G[i+1,end,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + 1 + (i-1)÷2 + nIJ*((k-1)÷2)
         elseif iSibling == 8
            iCell_G[i+1,end,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nI/2 + 1 + (i-1)÷2 + nIJ*((k-1)÷2)
         end
      end
   end

   # -z face
   if DiLevelNei_III[2,2,1] == 1
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[2,2,1], nBlock_P)

      @inbounds for j = 1:nJ, i = 1:nI
         iSibling = getSibling(iNodeNei_III, iTree_IA)

         if iSibling == 1
            iCell_G[i+1,j+1,1] = nIJK*neiBlock -
               nIJ + 1 + (i-1)÷2 + nI*((j-1)÷2)
         elseif iSibling == 2
            iCell_G[i+1,j+1,1] = nIJK*neiBlock -
               nIJ + nI/2 + 1 + (i-1)÷2 + nI*((j-1)÷2)
         elseif iSibling == 3
            iCell_G[i+1,j+1,1] = nIJK*neiBlock -
               nIJ/2 + 1 + (i-1)÷2 + nI*((j-1)÷2)
         elseif iSibling == 4
            iCell_G[i+1,j+1,1] = nIJK*neiBlock -
               nIJ/2 + nI/2 + 1 + (i-1)÷2 + nI*((j-1)÷2)
         end
      end
   end

   # +z face
   if DiLevelNei_III[2,2,3] == 0
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[3,3,4], nBlock_P)

      @inbounds for j = 1:nJ, i = 1:nI
         iCell_G[i+1,j+1,end] = nIJK*(neiBlock-1) + i + nI*(j-1)
      end
   elseif DiLevelNei_III[2,2,3] == 1
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[3,3,4], nBlock_P)

      @inbounds for j = 1:nJ, i = 1:nI
         iSibling = getSibling(iNodeNei_III, iTree_IA)

         if iSibling == 5
            iCell_G[i+1,j+1,end] = nIJK*(neiBlock-1) +
               1 + (i-1)÷2 + nI*((j-1)÷2)
         elseif iSibling == 6
            iCell_G[i+1,j+1,end] = nIJK*(neiBlock-1) +
               nI/2 + 1 + (i-1)÷2 + nI*((j-1)÷2)
         elseif iSibling == 7
            iCell_G[i+1,j+1,end] = nIJK*(neiBlock-1) +
               nIJ/2 + 1 + (i-1)÷2 + nI*((j-1)÷2)
         elseif iSibling == 8
            iCell_G[i+1,j+1,end] = nIJK*(neiBlock-1) +
               nIJ/2 + nI/2 + 1 + (i-1)÷2 + nI*((j-1)÷2)
         end
      end
   end

   ## Edges, in total 12

   # edge 1
   if DiLevelNei_III[2,1,1] == 0
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[2,1,1], nBlock_P)

      @inbounds for i = 1:nI
         iCell_G[i+1,1,1] = nIJK*neiBlock - nI + i
      end
   elseif DiLevelNei_III[2,1,1] in (1,2)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[2,1,1], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      iAMR = 2*DiLevelNei_III[2,1,1]

      if iSibling == 1
         @inbounds for i = 1:nI
            iCell_G[i+1,1,1] = nIJK*neiBlock -
               nI + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 2
         @inbounds for i = 1:nI
            iCell_G[i+1,1,1] = nIJK*neiBlock -
               nI/2 + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 3
         @inbounds for i = 1:nI
            iCell_G[i+1,1,1] = nIJK*neiBlock -
               nIJ/2 - nI + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 4
         @inbounds for i = 1:nI
            iCell_G[i+1,1,1] = nIJK*neiBlock -
               nIJ/2 - nI/2 + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 5
         @inbounds for i = 1:nI
            iCell_G[i+1,1,1] = nIJK*(neiBlock-1) +
               nIJ*(nK/2-1) + nI*(nJ-1) + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 6
         @inbounds for i = 1:nI
            iCell_G[i+1,1,1] = nIJK*(neiBlock-1) +
               nIJ*(nK/2-1) + nI*(nJ-1) + nI/2 + 1 + fld(i-1, iAMR)
         end
      end
   end

   # edge 2
   if DiLevelNei_III[2,3,1] == 0
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[2,4,1], nBlock_P)

      @inbounds for i = 1:nI
         iCell_G[i+1,end,1] = nIJK*neiBlock - nIJ + i
      end
   elseif DiLevelNei_III[2,3,1] in (1,2)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[2,4,1], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      iAMR = 2*DiLevelNei_III[2,3,1]

      if iSibling == 1
         @inbounds for i = 1:nI
            iCell_G[i+1,end,1] = nIJK*neiBlock -
               nIJ/2 + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 2
         @inbounds for i = 1:nI
            iCell_G[i+1,end,1] = nIJK*neiBlock -
               nIJ/2 + nI/2 + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 3
         @inbounds for i = 1:nI
            iCell_G[i+1,end,1] = nIJK*neiBlock -
               nIJ + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 4
         @inbounds for i = 1:nI
            iCell_G[i+1,end,1] = nIJK*neiBlock -
               nIJ + nI/2 + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 7
         @inbounds for i = 1:nI
            iCell_G[i+1,end,1] = nIJK*(neiBlock-1) +
               nIJ*(nK/2-1) + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 8
         @inbounds for i = 1:nI
            iCell_G[i+1,end,1] = nIJK*(neiBlock-1) +
               nIJ*(nK/2-1) + nI/2 + 1 + fld(i-1, iAMR)
         end
      end
   end

   # edge 3
   if DiLevelNei_III[2,1,3] == 0
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[2,1,4], nBlock_P)

      @inbounds for i = 1:nI
         iCell_G[i+1,1,end] = nIJK*(neiBlock-1) + nI*(nJ-1) + i
      end
   elseif DiLevelNei_III[2,1,3] in (1,2)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[2,1,4], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      iAMR = 2*DiLevelNei_III[2,1,3]

      if iSibling == 1
         @inbounds for i = 1:nI
            iCell_G[i+1,1,end] = nIJK*(neiBlock-1) +
               nIJK/2 + nI*(nJ-1) + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 2
         @inbounds for i = 1:nI
            iCell_G[i+1,1,end] = nIJK*(neiBlock-1) +
               nIJK/2 + nI*(nJ-1) + nI/2 + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 5
         @inbounds for i = 1:nI
            iCell_G[i+1,1,end] = nIJK*(neiBlock-1) +
               nI*(nJ-1) + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 6
         @inbounds for i = 1:nI
            iCell_G[i+1,1,end] = nIJK*(neiBlock-1) +
               nI*(nJ-1) + nI/2 + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 7
         @inbounds for i = 1:nI
            iCell_G[i+1,1,end] = nIJK*(neiBlock-1) +
               nI*(nJ/2-1) + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 8
         @inbounds for i = 1:nI
            iCell_G[i+1,1,end] = nIJK*(neiBlock-1) +
               nI*(nJ/2-1) + nI/2 + 1 + fld(i-1, iAMR)
         end
      end
   end

   # edge 4
   if DiLevelNei_III[2,3,3] == 0
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[2,4,4], nBlock_P)
      @inbounds for i = 1:nI
         iCell_G[i+1,end,end] = nIJK*(neiBlock-1) + i
      end
   elseif DiLevelNei_III[2,3,3] in (1,2)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[2,4,4], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      iAMR = 2^DiLevelNei_III[2,3,3]

      if iSibling == 3
         @inbounds for i = 1:nI
            iCell_G[i+1,end,end] = nIJK*(neiBlock-1) +
               nIJK/2 + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 4
         @inbounds for i = 1:nI
            iCell_G[i+1,end,end] = nIJK*(neiBlock-1) +
               nIJK/2 + nI/2 + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 5
         @inbounds for i = 1:nI
            iCell_G[i+1,end,end] = nIJK*(neiBlock-1) +
               nIJ/2 + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 6
         @inbounds for i = 1:nI
            iCell_G[i+1,end,end] = nIJK*(neiBlock-1) +
               nIJ/2 + nI/2 + 1 + fld(i-1, iAMR)
         end
      elseif iSibling == 7
         @inbounds for i = 1:nI
            iCell_G[i+1,end,end] = nIJK*(neiBlock-1) +
               1 + fld(i-1, iAMR)
         end
      elseif iSibling == 8
         @inbounds for i = 1:nI
            iCell_G[i+1,end,end] = nIJK*(neiBlock-1) +
               nI/2 + 1 + fld(i-1, iAMR)
         end
      end
   end

   # edge 5
   if DiLevelNei_III[1,2,1] == 0
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,2,1], nBlock_P)

      @inbounds for j = 1:nJ
         iCell_G[1,j+1,1] = nIJK*neiBlock - nIJ + nI*j
      end
   elseif DiLevelNei_III[1,2,1] in (1,2)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,2,1], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      iAMR = 2^DiLevelNei_III[1,2,1]

      if iSibling == 1
         @inbounds for j = 1:nJ
            iCell_G[1,j+1,1] = nIJK*neiBlock -
               nIJ + nI*(1 + fld(j-1,iAMR))
         end
      elseif iSibling == 2
         @inbounds for j = 1:nJ
            iCell_G[1,j+1,1] = nIJK*neiBlock -
               nIJ + nI/2 + nI*fld(j-1,iAMR)
         end
      elseif iSibling == 3
         @inbounds for j = 1:nJ
            iCell_G[1,j+1,1] = nIJK*neiBlock -
               nIJ/2 + nI*(1 + fld(j-1,iAMR))
         end
      elseif iSibling == 4
         @inbounds for j = 1:nJ
            iCell_G[1,j+1,1] = nIJK*neiBlock -
               nIJ/2 + nI/2 + nI*fld(j-1,iAMR)
         end
      elseif iSibling == 5
         @inbounds for j = 1:nJ
            iCell_G[1,j+1,1] = nIJK*neiBlock -
               nIJ*(nK/2+1) + nI*(1 + fld(j-1,iAMR))
         end
      elseif iSibling == 7
         @inbounds for j = 1:nJ
            iCell_G[1,j+1,1] = nIJK*neiBlock -
               nIJ*(nK/2+1) + nIJ/2 + nI*(1 + fld(j-1,iAMR))
         end
      end
   end

   # edge 6
   if DiLevelNei_III[3,2,1] == 0
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,2,1], nBlock_P)

      @inbounds for j = 1:nJ
         iCell_G[end,j+1,1] = nIJK*neiBlock - nIJ + 1 + nI*(j-1)
      end
   elseif DiLevelNei_III[3,2,1] in (1,2)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,2,1], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      iAMR = 2^DiLevelNei_III[3,2,1]

      if iSibling == 1
         @inbounds for j = 1:nJ
            iCell_G[end,j+1,1] = nIJK*neiBlock -
               nIJ + nI/2 + 1 + nI*fld(j-1,iAMR)
         end
      elseif iSibling == 2
         @inbounds for j = 1:nJ
            iCell_G[end,j+1,1] = nIJK*neiBlock -
               nIJ + 1 + nI*fld(j-1,iAMR)
         end
      elseif iSibling == 3
         @inbounds for j = 1:nJ
            iCell_G[end,j+1,1] = nIJK*neiBlock -
               nIJ/2 + nI/2 + 1 + nI*fld(j-1,iAMR)
         end
      elseif iSibling == 4
         @inbounds for j = 1:nJ
            iCell_G[end,j+1,1] = nIJK*neiBlock -
               nIJ/2 + 1 + nI*fld(j-1,iAMR)
         end
      elseif iSibling == 6
         @inbounds for j = 1:nJ
            iCell_G[end,j+1,1] = nIJK*(neiBlock-1) +
               nIJ*(nK/2-1) + 1 + nI*fld(j-1,iAMR)
         end
      elseif iSibling == 8
         @inbounds for j = 1:nJ
            iCell_G[end,j+1,1] = nIJK*(neiBlock-1) +
               nIJ*(nK/2-1) + nIJ/2 + 1 + nI*fld(j-1,iAMR)
         end
      end
   end

   # edge 7
   if DiLevelNei_III[1,2,3] == 0
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,2,4], nBlock_P)

      @inbounds for j = 1:nJ
         iCell_G[1,j+1,end] = nIJK*(neiBlock-1) + nI*j
      end
   elseif DiLevelNei_III[1,2,3] in (1,2)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,2,4], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      iAMR = 2^DiLevelNei_III[1,2,3]

      if iSibling == 1
         @inbounds for j = 1:nJ
            iCell_G[1,j+1,end] = nIJK*(neiBlock-1) +
               nIJK/2 + nI*(1 + fld(j-1,iAMR))
         end
      elseif iSibling == 3
         @inbounds for j = 1:nJ
            iCell_G[1,j+1,end] = nIJK*(neiBlock-1) +
               nIJK/2 + nIJ/2 + nI*(1 + fld(j-1,iAMR))
         end
      elseif iSibling == 5
         @inbounds for j = 1:nJ
            iCell_G[1,j+1,end] = nIJK*(neiBlock-1) +
               nI*(1 + fld(j-1,iAMR))
         end
      elseif iSibling == 6
         @inbounds for j = 1:nJ
            iCell_G[1,j+1,end] = nIJK*(neiBlock-1) +
               nI*(0.5 + fld(j-1,iAMR))
         end
      elseif iSibling == 7
         @inbounds for j = 1:nJ
            iCell_G[1,j+1,end] = nIJK*(neiBlock-1) +
               nIJ/2 + nI*(1 + fld(j-1,iAMR))
         end
      elseif iSibling == 8
         @inbounds for j = 1:nJ
            iCell_G[1,j+1,end] = nIJK*(neiBlock-1) +
               nIJ/2 + nI*(0.5 + fld(j-1,iAMR))
         end
      end
   end

   # edge 8
   if DiLevelNei_III[3,2,3] == 0
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,2,4], nBlock_P)
      @inbounds for j = 1:nJ
         iCell_G[end,j+1,end] = nIJK*(neiBlock-1) + 1 + nI*(j-1)
      end
   elseif DiLevelNei_III[3,2,3] in (1,2)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,2,4], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      iAMR = 2^DiLevelNei_III[3,2,3]

      if iSibling == 2
         @inbounds for j = 1:nJ
            iCell_G[end,j+1,end] = nIJK*(neiBlock-1) +
               nIJK/2 + 1 + nI*fld(j-1,iAMR)
         end
      elseif iSibling == 4
         @inbounds for j = 1:nJ
            iCell_G[end,j+1,end] = nIJK*(neiBlock-1) +
               nIJK/2 + nIJ/2 + 1 + nI*fld(j-1,iAMR)
         end
      elseif iSibling == 5
         @inbounds for j = 1:nJ
            iCell_G[end,j+1,end] = nIJK*(neiBlock-1) +
               nI/2 + 1 + nI*fld(j-1,iAMR)
         end
      elseif iSibling == 6
         @inbounds for j = 1:nJ
            iCell_G[end,j+1,end] = nIJK*(neiBlock-1) +
               1 + nI*fld(j-1,iAMR)
         end
      elseif iSibling == 7
         @inbounds for j = 1:nJ
            iCell_G[end,j+1,end] = nIJK*(neiBlock-1) +
               nIJ/2 + nI/2 + 1 + nI*fld(j-1,iAMR)
         end
      elseif iSibling == 8
         @inbounds for j = 1:nJ
            iCell_G[end,j+1,end] = nIJK*(neiBlock-1) +
               nIJ/2 + 1 + nI*fld(j-1,iAMR)
         end
      end
   end

   # edge 9
   if DiLevelNei_III[1,1,2] == 0
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,1,2], nBlock_P)

      @inbounds for k = 1:nK
         iCell_G[1,1,k+1] = nIJK*(neiBlock-1) + nIJ*k
      end
   elseif DiLevelNei_III[1,1,2] in (1,2)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,1,2], nBlock_P)

      iAMR = 2^DiLevelNei_III[1,1,2]

      iSibling = getSibling(iNodeNei_III, iTree_IA)
      if iSibling == 1
         @inbounds for k = 1:nK
            iCell_G[1,1,k+1] = nIJK*(neiBlock-1) +
               nIJ*(1 + fld(k-1,iAMR))
         end
      elseif iSibling == 2
         @inbounds for k = 1:nK
            iCell_G[1,1,k+1] = nIJK*(neiBlock-1) +
               nI*(nJ-1) + nI/2 + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 3
         @inbounds for k = 1:nK
            iCell_G[1,1,k+1] = nIJK*(neiBlock-1) +
               nIJ/2 + nIJ*(1 + fld(k-1,iAMR))
         end
      elseif iSibling == 5
         @inbounds for k = 1:nK
            iCell_G[1,1,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nIJ*(1 + fld(k-1,iAMR))
         end
      elseif iSibling == 6
         @inbounds for k = 1:nK
            iCell_G[1,1,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nI*(nJ-1) + nI/2 + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 7
         @inbounds for k = 1:nK
            iCell_G[1,1,k+1] = nIJK*(neiBlock-1) +
               nIJ*nJ/2 + nIJ/2 + nIJ*fld(k-1,iAMR)
         end
      end
   end

   # edge 10
   if DiLevelNei_III[3,1,2] == 0
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,1,2], nBlock_P)

      @inbounds for k = 1:nK
         iCell_G[end,1,k+1] = nIJK*(neiBlock-1) + nI*(nJ-1) + 1 + nIJ*(k-1)
      end
   elseif DiLevelNei_III[3,1,2] in (1,2)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,1,2], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      iAMR = 2^DiLevelNei_III[3,1,2]

      if iSibling == 1
         @inbounds for k = 1:nK
            iCell_G[end,1,k+1] = nIJK*(neiBlock-1) +
               nI*(nJ-1) + nI/2 + 1 + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 2
         @inbounds for k = 1:nK
            iCell_G[end,1,k+1] = nIJK*(neiBlock-1) +
               nI*(nJ-1) + 1 + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 4
         @inbounds for k = 1:nK
            iCell_G[end,1,k+1] = nIJK*(neiBlock-1) +
               nI*(nJ/2-1) + 1 + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 5
         @inbounds for k = 1:nK
            iCell_G[end,1,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nI*(nJ-1) + nI/2 + 1 + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 6
         @inbounds for k = 1:nK
            iCell_G[end,1,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nI*(nJ-1) + 1 + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 8
         @inbounds for k = 1:nK
            iCell_G[end,1,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nI*(nJ/2-1) + 1 + nIJ*fld(k-1,iAMR)
         end
      end
   end

   # edge 11
   if DiLevelNei_III[1,3,2] == 0
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,4,2], nBlock_P)

      @inbounds for k = 1:nK
         iCell_G[1,end,k+1] = nIJK*(neiBlock-1) + nI + nIJ*(k-1)
      end
   elseif DiLevelNei_III[1,3,2] in (1,2)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,4,2], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      iAMR = 2^DiLevelNei_III[1,3,2]

      if iSibling == 1
         @inbounds for k = 1:nK
            iCell_G[1,end,k+1] = nIJK*(neiBlock-1) +
               nIJ/2 + nI + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 3
         @inbounds for k = 1:nK
            iCell_G[1,end,k+1] = nIJK*(neiBlock-1) +
               nI + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 4
         @inbounds for k = 1:nK
            iCell_G[1,end,k+1] = nIJK*(neiBlock-1) +
               nI/2 + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 5
         @inbounds for k = 1:nK
            iCell_G[1,end,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nIJ/2 + nI + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 7
         @inbounds for k = 1:nK
            iCell_G[1,end,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nI + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 8
         @inbounds for k = 1:nK
            iCell_G[1,end,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nI/2 + nIJ*fld(k-1,iAMR)
         end
      end
   end

   # edge 12
   if DiLevelNei_III[3,3,2] == 0
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,4,2], nBlock_P)
      @inbounds for k = 1:nK
         iCell_G[end,end,k+1] = nIJK*(neiBlock-1) + 1 + nIJ*(k-1)
      end
   elseif DiLevelNei_III[3,3,2] in (1,2)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,4,2], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      iAMR = 2^DiLevelNei_III[3,3,2]

      if iSibling == 2
         @inbounds for k = 1:nK
            iCell_G[end,end,k+1] = nIJK*(neiBlock-1) +
               nIJ/2 + 1 + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 3
         @inbounds for k = 1:nK
            iCell_G[end,end,k+1] = nIJK*(neiBlock-1) +
               nI/2 + 1 + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 4
         @inbounds for k = 1:nK
            iCell_G[end,end,k+1] = nIJK*(neiBlock-1) +
               1 + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 6
         @inbounds for k = 1:nK
            iCell_G[end,end,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nIJ/2 + 1 + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 7
         @inbounds for k = 1:nK
            iCell_G[end,end,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + nI/2 + 1 + nIJ*fld(k-1,iAMR)
         end
      elseif iSibling == 8
         @inbounds for k = 1:nK
            iCell_G[end,end,k+1] = nIJK*(neiBlock-1) +
               nIJK/2 + 1 + nIJ*fld(k-1,iAMR)
         end
      end
   end

   ## Node, in total 8

   # node 1
   if DiLevelNei_III[1,1,1] in (0,2,3)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,1,1], nBlock_P)

      iCell_G[1,1,1] = nIJK*neiBlock

   elseif DiLevelNei_III[1,1,1] == 1
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,1,1], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      if iSibling == 1
         iCell_G[1,1,1] = nIJK*neiBlock
      elseif iSibling == 2
         iCell_G[1,1,1] = nIJK*neiBlock - nI/2
      elseif iSibling == 3
         iCell_G[1,1,1] = nIJK*neiBlock - nIJ/2
      elseif iSibling == 4
         iCell_G[1,1,1] = nIJK*neiBlock - nIJ/2 - nI/2
      elseif iSibling == 5
         iCell_G[1,1,1] = nIJK*neiBlock - nIJK/2
      elseif iSibling == 6
         iCell_G[1,1,1] = nIJK*neiBlock - nIJK/2 - nI/2
      elseif iSibling == 7
         iCell_G[1,1,1] = nIJK*neiBlock - nIJK/2 - nIJ/2
      end
   end

   # node 2
   if DiLevelNei_III[3,1,1] in (0,2,3)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,1,1], nBlock_P)

      iCell_G[end,1,1] = nIJK*neiBlock - nI + 1

   elseif DiLevelNei_III[3,1,1] == 1
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,1,1], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      if iSibling == 1
         iCell_G[end,1,1] = nIJK*neiBlock - nI/2 + 1
      elseif iSibling == 2
         iCell_G[end,1,1] = nIJK*neiBlock - nI + 1
      elseif iSibling == 3
         iCell_G[end,1,1] = nIJK*neiBlock - nIJ/2 - nI/2 + 1
      elseif iSibling == 4
         iCell_G[end,1,1] = nIJK*neiBlock - nIJ/2 - nI + 1
      elseif iSibling == 5
         iCell_G[end,1,1] = nIJK*neiBlock - nIJK/2 - nI/2 + 1
      elseif iSibling == 6
         iCell_G[end,1,1] = nIJK*neiBlock - nIJK/2 - nI + 1
      elseif iSibling == 8
         iCell_G[end,1,1] = nIJK*neiBlock - nIJK/2 - nIJ/2 - nI + 1
      end
   end

   # node 3
   if DiLevelNei_III[1,3,1] in (0,2,3)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,4,1], nBlock_P)

      iCell_G[1,end,1] = nIJK*neiBlock - nIJ + nI

   elseif DiLevelNei_III[1,3,1] == 1
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,4,1], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      if iSibling == 1
         iCell_G[1,end,1] = nIJK*neiBlock - nIJ/2 + nI
      elseif iSibling == 2
         iCell_G[1,end,1] = nIJK*neiBlock - nIJ/2 + nI/2
      elseif iSibling == 3
         iCell_G[1,end,1] = nIJK*neiBlock - nIJ + nI
      elseif iSibling == 4
         iCell_G[1,end,1] = nIJK*neiBlock - nIJ + nI/2
      elseif iSibling == 5
         iCell_G[1,end,1] = nIJK*neiBlock - nIJK/2 - nIJ/2 + nI
      elseif iSibling == 7
         iCell_G[1,end,1] = nIJK*neiBlock - nIJK/2 - nIJ + nI
      elseif iSibling == 8
         iCell_G[1,end,1] = nIJK*neiBlock - nIJK/2 - nIJ + nI/2
      end
   end

   # node 4
   if DiLevelNei_III[3,3,1] in (0,2,3)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,4,1], nBlock_P)

      iCell_G[end,end,1] = nIJK*neiBlock - nIJ + 1

   elseif DiLevelNei_III[3,3,1] == 1
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,4,1], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      if iSibling == 1
         iCell_G[end,end,1] = nIJK*neiBlock - nIJ/2 + nI/2 + 1
      elseif iSibling == 2
         iCell_G[end,end,1] = nIJK*neiBlock - nIJ/2 + 1
      elseif iSibling == 3
         iCell_G[end,end,1] = nIJK*neiBlock - nIJ + nI/2 + 1
      elseif iSibling == 4
         iCell_G[end,end,1] = nIJK*neiBlock - nIJ + 1
      elseif iSibling == 6
         iCell_G[end,end,1] = nIJK*neiBlock - nIJK/2 - nIJ/2 + 1
      elseif iSibling == 7
         iCell_G[end,end,1] = nIJK*neiBlock - nIJK/2 - nIJ + nI/2 + 1
      elseif iSibling == 8
         iCell_G[end,end,1] = nIJK*neiBlock - nIJK/2 - nIJ + 1
      end
   end

   # node 5
   if DiLevelNei_III[1,1,3] in (0,2,3)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,1,4], nBlock_P)

      iCell_G[1,1,end] = nIJK*(neiBlock-1) + nIJ

   elseif DiLevelNei_III[1,1,3] == 1
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,1,4], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      if iSibling == 1
         iCell_G[1,1,end] = nIJK*(neiBlock-1) + nIJK/2 + nIJ
      elseif iSibling == 2
         iCell_G[1,1,end] = nIJK*(neiBlock-1) + nIJK/2 + nI*(nJ-1) + nI/2
      elseif iSibling == 3
         iCell_G[1,1,end] = nIJK*(neiBlock-1) + nIJK/2 + nIJ/2
      elseif iSibling == 5
         iCell_G[1,1,end] = nIJK*(neiBlock-1) + nIJ
      elseif iSibling == 6
         iCell_G[1,1,end] = nIJK*(neiBlock-1) + nI*(nJ-1) + nI/2
      elseif iSibling == 7
         iCell_G[1,1,end] = nIJK*(neiBlock-1) + nIJ/2
      elseif iSibling == 8
         iCell_G[1,1,end] = nIJK*(neiBlock-1) + nI*(nJ/2-1) + nI/2
      end
   end

   # node 6
   if DiLevelNei_III[3,1,3] in (0,2,3)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,1,4], nBlock_P)

      iCell_G[end,1,end] = nIJK*(neiBlock-1) + nI*(nJ-1) + 1

   elseif DiLevelNei_III[3,1,3] == 1
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,1,4], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      if iSibling == 1
         iCell_G[end,1,end] = nIJK*(neiBlock-1) + nIJK/2 + nI*(nJ-1) + nI/2 + 1
      elseif iSibling == 2
         iCell_G[end,1,end] = nIJK*(neiBlock-1) + nIJK/2 + nI*(nJ-1) + 1
      elseif iSibling == 4
         iCell_G[end,1,end] = nIJK*(neiBlock-1) + nIJK/2 + nI*(nJ/2-1) + 1
      elseif iSibling == 5
         iCell_G[end,1,end] = nIJK*(neiBlock-1) + nI*(nJ-1) + nI/2 + 1
      elseif iSibling == 6
         iCell_G[end,1,end] = nIJK*(neiBlock-1) + nI*(nJ-1) + 1
      elseif iSibling == 7
         iCell_G[end,1,end] = nIJK*(neiBlock-1) + nI*(nJ/2-1) + nI/2 + 1
      elseif iSibling == 8
         iCell_G[end,1,end] = nIJK*(neiBlock-1) + nI*(nJ/2-1) + 1
      end
   end

   # node 7
   if DiLevelNei_III[1,3,3] in (0,2,3)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,4,4], nBlock_P)

      iCell_G[1,end,end] = nIJK*(neiBlock-1) + nI

   elseif DiLevelNei_III[1,3,3] == 1
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[1,4,4], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      if iSibling == 1
         iCell_G[1,end,end] = nIJK*(neiBlock-1) + nIJK/2 + nIJ/2 + nI
      elseif iSibling == 3
         iCell_G[1,end,end] = nIJK*(neiBlock-1) + nIJK/2 + nI
      elseif iSibling == 4
         iCell_G[1,end,end] = nIJK*(neiBlock-1) + nIJK/2 + nI/2
      elseif iSibling == 5
         iCell_G[1,end,end] = nIJK*(neiBlock-1) + nIJ/2 + nI
      elseif iSibling == 6
         iCell_G[1,end,end] = nIJK*(neiBlock-1) + nIJ/2 + nI/2
      elseif iSibling == 7
         iCell_G[1,end,end] = nIJK*(neiBlock-1) + nI
      elseif iSibling == 8
         iCell_G[1,end,end] = nIJK*(neiBlock-1) + nI/2
      end
   end

   # node 8
   if DiLevelNei_III[3,3,3] in (0,2,3)
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,4,4], nBlock_P)

      iCell_G[end,end,end] = nIJK*(neiBlock-1) + 1

   elseif DiLevelNei_III[3,3,3] == 1
      neiBlock = nodeToGlobalBlock(batl, iNodeNei_III[4,4,4], nBlock_P)

      iSibling = getSibling(iNodeNei_III, iTree_IA)

      # Sibling 1 does not have to compute.
      if iSibling == 2
         iCell_G[end,end,end] = nIJK*(neiBlock-1) + nIJK/2 + nIJ/2 + 1
      elseif iSibling == 3
         iCell_G[end,end,end] = nIJK*(neiBlock-1) + nIJK/2 + nI/2 + 1
      elseif iSibling == 4
         iCell_G[end,end,end] = nIJK*(neiBlock-1) + nIJK/2 + 1
      elseif iSibling == 5
         iCell_G[end,end,end] = nIJK*(neiBlock-1) + nIJ/2 + nI/2 + 1
      elseif iSibling == 6
         iCell_G[end,end,end] = nIJK*(neiBlock-1) + nIJ/2 + 1
      elseif iSibling == 7
         iCell_G[end,end,end] = nIJK*(neiBlock-1) + nI/2 + 1
      elseif iSibling == 8
         iCell_G[end,end,end] = nIJK*(neiBlock-1) + 1
      end
   end
   nothing
end

"""
    create_pvd(filepattern)

Generate PVD file for a time series collection of VTK data.
# Example
```
create_pvd("*.vtu)
```
"""
function create_pvd(filepattern)

   filenames = glob(filepattern)

   # create an empty XML document
   xdoc = XMLDocument()

   # create & attach a root node
   xroot = create_root(xdoc, "VTKFile")

   type = "Collection"
   byte_order = "LittleEndian"
   compressor = "vtkZLibDataCompressor"

   set_attributes(xroot; type, byte_order, compressor)

   # create the first child
   xs1 = new_child(xroot, "Collection")

   for file in filenames
      i_end  = findfirst("_n",file)[1] - 1

      second = parse(Int32, file[i_end-1:i_end])
      minute = parse(Int32, file[i_end-3:i_end-2])
      timestep = 60*minute + second
      @show timestep

      e = new_child(xs1, "DataSet")

      set_attributes(e; timestep, group="", part="0", file)
   end

   name_index = findfirst("_t",filenames[1])[1] - 1

   save_file(xdoc, filenames[1][1:name_index]*".pvd")
   nothing
end


text:
module mT4quartercylnm
using FinEtools
using FinEtools.MeshExportModule: VTK
using Test
function test()
	Radius, Length, nperradius, nL = 3.0, 5.0, 3, 5
	fens,fes = T4quartercyln(Radius, Length, nperradius, nL)
	@test  count(fes) == 270
	File = "mesh.vtk"
	VTK.vtkexportmesh(File, fens, fes)
	rm(File)
	# @async run(`"paraview.exe" $File`)
	bfes = meshboundary(fes)
	@test  count(bfes) == 156
	File = "bmesh.vtk"
	VTK.vtkexportmesh(File, fens, fes)
	rm(File)
	# @async run(`"paraview.exe" $File`)
end
end
using .mT4quartercylnm
mT4quartercylnm.test()


text:
# MCState holds two samples at each Monte Carlo iteration, the current and the successive one

type MCState{S<:MCSample}
  successive::S # If proposed sample is accepted, then successive = proposed, otherwise successive = current
  current::S
  diagnostics::Dict
end

MCState{S<:MCSample}(p::S, c::S) = MCState(p, c, Dict())


text:
struct Layer
    transformations::Tuple
    data::Any
    mappings::Arguments
end

Layer(transformations::Tuple=()) = Layer(transformations, nothing, arguments())

data(df) = Layer((), columns(df), arguments())
mapping(args...; kwargs...) = Layer((), nothing, arguments(args...; kwargs...))

function Base.:*(l1::Layer, l2::Layer)
    t1, t2 = l1.transformations, l2.transformations
    d1, d2 = l1.data, l2.data
    m1, m2 = l1.mappings, l2.mappings
    transformations = (t1..., t2...)
    data = isnothing(d2) ? d1 : d2
    mappings = Arguments(
        vcat(m1.positional, m2.positional),
        merge(m1.named, m2.named)
    )
    return Layer(transformations, data, mappings)
end


text:
using AtomicGraphNets
using Test

@testset "layer_tests" begin
    include("layer_tests.jl")
end

@testset "model_tests" begin
    include("model_tests.jl")
end


text:
# ==============================================================================
# applique une perturbation sur la solution entiere faisant l'objet d'un cycle

function perturbSolution!(vg::Vector{tGenerateur}, k::Int64, c1::Array{Int,1}, c2::Array{Int,1}, d::tListDisplay)
    # nombre de variables maximum a considerer
    T = 20
    # nombre effectif de variables a flipper
    TT = rand(T/2:3*T/2)

    # liste des candidats (valeur, indice) et tri decroissant
    nbvar = length(vg[k].sInt.x)
    candidats=[( abs( vg[k].sPrj.x[i] - vg[k].sInt.x[i] ) , i ) for i=1:nbvar]
    sort!(candidats, rev=true, by = x -> x[1])
#    sort!(candidats,  by = x -> x[1])

#@show vg[k].sPrj.x
#@show vg[k].sInt.x
#@show candidats
#@show TT
#@show nbvar

    i = 1
    while (i<= nbvar) && (i<=TT)
        j=candidats[i][2]
#        @show candidats[i][2]
        if vg[k].sInt.x[j] == 0
            vg[k].sInt.x[j] = 1
            vg[k].sInt.y[1] = vg[k].sInt.y[1] + c1[j]
            vg[k].sInt.y[2] = vg[k].sInt.y[2] + c2[j]
        else
            vg[k].sInt.x[j] = 0
            vg[k].sInt.y[1] = vg[k].sInt.y[1] - c1[j]
            vg[k].sInt.y[2] = vg[k].sInt.y[2] - c2[j]
        end
        i+=1
    end
    @printf("  %2dC : [ %8.2f , %8.2f ] \n", k, vg[k].sInt.y[1], vg[k].sInt.y[2])

    # archive le point obtenu pour les besoins d'affichage
    if generateurVisualise == -1
        # archivage pour tous les generateurs
        push!(d.XPert,vg[k].sInt.y[1])
        push!(d.YPert,vg[k].sInt.y[2])
    elseif generateurVisualise == k
        # archivage seulement pour le generateur k
        push!(d.XPert,vg[k].sInt.y[1])
        push!(d.YPert,vg[k].sInt.y[2])
    end
#    @show vg[k].sInt.x
    return nothing
end

# ==============================================================================
# applique une perturbation sur la solution entiere faisant l'objet d'un cycle

function perturbSolution30!(vg::Vector{tGenerateur}, k::Int64, c1::Array{Int,1}, c2::Array{Int,1}, d::tListDisplay)

    # liste des candidats (valeur, indice) et tri decroissant
    nbvar = length(vg[k].sInt.x)
    idxTilde0, idxTilde1 = split01(vg[k].sInt.x)

#    candidats=[( abs( vg[k].sPrj.x[i] - vg[k].sInt.x[i] ) , i ) for i=1:nbvar if vg[k].sPrj.x[i]>0 && vg[k].sPrj.x[i]<1]
    candidats=[( vg[k].sPrj.x[i] , i ) for i=1:nbvar if vg[k].sPrj.x[i]>0 && vg[k].sPrj.x[i]<1]
#    sort!(candidats, rev=true, by = x -> x[1])

#@show vg[k].sPrj.x
#@show vg[k].sInt.x
#@show candidats
#@show nbvar

seq = randperm(length(candidats)) # melange les candidats afin d'avoir une composante variee
etat = vg[k].sInt.x[ candidats[seq[1]][2] ] # etat 0 ou 1 de la premiere variable candidate
for i = 1:length(candidats)
    j=candidats[seq[i]][2]
    #@show candidats[seq[i]][2]
    if etat == 0
        if vg[k].sInt.x[j] == 0
            vg[k].sInt.x[j] = 1
            vg[k].sInt.y[1] = vg[k].sInt.y[1] + c1[j]
            vg[k].sInt.y[2] = vg[k].sInt.y[2] + c2[j]
        end
    else
        if vg[k].sInt.x[j] == 1
            vg[k].sInt.x[j] = 0
            vg[k].sInt.y[1] = vg[k].sInt.y[1] - c1[j]
            vg[k].sInt.y[2] = vg[k].sInt.y[2] - c2[j]
        end
    end
    etat=(etat+1)%2
end

    @printf("  %2dC : [ %8.2f , %8.2f ] \n", k, vg[k].sInt.y[1], vg[k].sInt.y[2])

    # archive le point obtenu pour les besoins d'affichage
    if generateurVisualise == -1
        # archivage pour tous les generateurs
        push!(d.XPert,vg[k].sInt.y[1])
        push!(d.YPert,vg[k].sInt.y[2])
    elseif generateurVisualise == k
        # archivage seulement pour le generateur k
        push!(d.XPert,vg[k].sInt.y[1])
        push!(d.YPert,vg[k].sInt.y[2])
    end
#    @show vg[k].sInt.x

    return nothing
end

max_stars_repo_name: ducanh-le/GravityMachine

max_issues_repo_name: ducanh-le/GravityMachine

max_forks_repo_name: ducanh-le/GravityMachine


text:
press_data(T::Type, y) = [T.(yi) for yi in y]
compress_data(y) = press_data(Float16, y)
decompress_data(y) = press_data(Float64, y)


text:
## JL-ML SPGL Comparison
export scale_cs

"""
An example of using implicit JOLI operators to solve a system using GenSPGL

Use: x, r, g, info, SNR = ex_cs()
"""
function scale_cs(n,k,i)

    #Create a sparse vector
    #inds = randperm(n)
    #x0 = zeros(n)
    #keep = Int(floor(n/1000))
    #x0[inds[1:keep]] = randn(keep)
    #r_inds = inds[1:k]

    #= Load sparse vector x and restriction mask
    fid = matopen(Pkg.dir("GenSPGL")*"/src/Examples/scaletest/x0.mat")
    x0 = vec(read(fid, "x0"))
    r_inds = Int.(vec(read(fid, "r_inds")))
    n = length(vec(x0))
    close(fid)
    =#

    # Load sparse vector x and restriction mask
    fid = matopen(Pkg.dir("GenSPGL")*"/src/Examples/scaletest/x0.mat")
    x0 = vec(read(fid, "x0"))
    r_inds = Int.(vec(read(fid, "r_inds")))
    n = length(vec(x0))
    close(fid)

    # Sparsity Promoting Transform
    S = joDCT(n)

    # Restriction Operator
    R = joRestriction(n, r_inds)

    # Modelling op
    A = R*S

    # Create data
    b = A*x0

    # Set BLAS Threads
    BLAS.set_num_threads(8)

    # Solve
    opts = spgOptions(verbosity = 1, iterations = 100)
    t = [@elapsed spgl1(A, b, options = opts) for ii = 1:i]
    return t

end


text:
using FixedEffectModels, DataFrames, CSV, Test

df = CSV.read(joinpath(dirname(@__FILE__), "..", "dataset/Cigar.csv"))
df[:pState] = categorical(df[:State])
df[:pYear] = categorical(df[:Year])

model = @model Sales ~ NDI
result = reg(df, model)
show(result)
predict(result, df)
residuals(result, df)

model = @model Sales ~ CPI + (Price ~ Pimin)
result = reg(df, model)
show(result)
predict(result, df)
residuals(result, df)
response(result, df)
@test  nobs(result) == 1380
@test vcov(result)[1] ≈ 3.5384578251636785

model = @model Sales ~ Price  fe = pState
result = reg(df, model)
show(result)
model = @model Sales ~ CPI + (Price ~ Pimin) fe = pState
result = reg(df, model)
show(result)

model = @model Sales ~ Price
result = reg(df, model)
@test residuals(result, df)[1:10] ≈ [-39.2637, -37.48801, -34.38801, -36.09743, -36.97446, -43.15547, -41.22573, -40.83648, -34.52427, -28.91617] atol = 1e-4

model = @model Sales ~ Price fe = pState
result = reg(df, model, save = true)
@test result.augmentdf[:residuals][1:10] ≈ [-22.08499, -20.33318, -17.23318, -18.97645, -19.85547, -26.1161, -24.20627, -23.87674, -17.62624, -12.01018] atol = 1e-4

# test diffrent specification sf save

model = @model Sales ~ Price fe = pState
result = reg(df, model, save = true)
@test :residuals ∈ names(result.augmentdf)
@test :pState ∈ names(result.augmentdf)

model = @model Sales ~ Price fe = pState
result = reg(df, model, save = :residuals)
@test :residuals ∈ names(result.augmentdf)
@test :pState ∉ names(result.augmentdf)

model = @model Sales ~ Price fe = pState
result = reg(df, model, save = :fe)
@test :residuals ∉ names(result.augmentdf)
@test :pState ∈ names(result.augmentdf)


text:
######################################################################################################
# Merge from                                                                                #
# * Wouter Koolen: https://bitbucket.org/wmkoolen/tidnabbil/src/master/purex_games_paper/, and       #
# * Xuedong Shang: https://github.com/xuedong/LinBAI.jl                                              #                                      #
######################################################################################################
struct FTL
    L :: Array{Float64, 1};

    function FTL(K)
        new(zeros(K));
    end
end
function act(h::FTL)
    u = h.L .== minimum(h.L);
    u ./ sum(u);
end
function incur!(h::FTL, l)
    h.L .+= l;
end

mutable struct AdaHedge
    L :: Array{Float64, 1};
    Δ :: Float64;

    function AdaHedge(K)
        new(zeros(K), 0.01);
    end
end
function act(h::AdaHedge)
    η = log(length(h.L))/h.Δ;
    u = exp.(-η.*(h.L .- minimum(h.L)));
    u ./ sum(u);
end
function incur!(h::AdaHedge, l)
    K = length(l);
    u = act(h);
    η = log(length(h.L))/h.Δ;
    M(L) = η==0 ? sum(L)/K : minimum(L) - 1/η*log(sum(exp.(-η.*(L .- minimum(L))))/K);
    Mpre = M(h.L);
    h.L .+= l;
    Mpst = M(h.L);
    m = Mpst-Mpre;

    @assert u'l ≥ m - 1e-7 "hedge loss $(u'l), Δ=$(h.Δ), η=$η, mix loss $m, u=$u, L=$(h.L), l=$(l)";
    h.Δ  += u'l - m;
end

######################################################################################################
# Implementation by Xuedong Shang: https://github.com/xuedong/LinBAI.jl                              #                                                                     #
######################################################################################################
mutable struct LinBAIAdaHedge
    L :: Array{Float64, 1};
    Δ :: Float64;

    function LinBAIAdaHedge(K)
        new(zeros(K), 0.01);
    end
end
function act(h::LinBAIAdaHedge)
    η = log(length(h.L))/h.Δ;
    u = exp.(-η.*(h.L .- minimum(h.L)));
    u ./ sum(u);
end
function incur!(h::LinBAIAdaHedge, l)
    u = act(h);
    η = log(length(h.L))/h.Δ;
    h.L .+= l;
    m = minimum(l) - 1/η*log(u'exp.(-η.*(l .- minimum(l))));
    # println(η, " ", u, " ", u'l, " ", -1/η*log(u'exp.(-η.*l)));
    @assert isfinite(m) && u'l ≥ m - 1e-7 "hedge loss $(u'l), Δ=$(h.Δ), η=$η, mix loss $m, u=$u, L=$(h.L)";
    h.Δ  += u'l - m;
end

# Menard calls this "Lazy Mirror Descent" As perhaps the main feature
# here is the mixing step, we call it "Fixed Share".  This is the
# version with slowly decreasing O(1/sqrt(t)) learning and switching
# rates.  The switching is used to force exploration
mutable struct FixedShare
    w :: Array{Float64, 1};
    t; # time
    S; # loss range: S ≥ \max_k l_k - \min_k l_k
    function FixedShare(K; S=1)
        new(fill(1/K,K), 0, S);
    end
end
function act(h::FixedShare)
    h.w
end
function incur!(h::FixedShare, l)
    K = length(l);
    h.t += 1;
    # Exponential Weights update
    η = sqrt(log(K)/h.t)/h.S;   # decaying learning rate
    h.w .*= exp.(-η*(l .- minimum(l)));
    h.w ./= sum(h.w);
    # Fixed Share (to unifor prior)
    γ = 1/(4*sqrt(h.t)); # decaying switching rate
    h.w .= (1-γ).*h.w .+ γ/K;
end


text:
module Import

using ACE, JuLIP, JuLIP.Testing
using JuLIP: evaluate, evaluate_d
using ACE: PIBasisFcn, DAGEvaluator
using ACE.RPI: RnYlm1pBasis, RnYlmBasisFcn

function import_rbasis_v05(D, rtests = [])
   @assert D["__id__"] == "SHIPs_TransformedJacobi"
   trans = ACE.Transforms.PolyTransform(D["trans"])
   ru, rl = D["ru"], D["rl"]
   tu, tl = trans(ru), trans(rl)
   cutoff_id = D["cutoff"]["__id__"]
   if cutoff_id == "SHIPs_PolyCutoff1s"
      pl = 0
      pu = D["cutoff"]["P"]
   elseif cutoff_id == "SHIPs_PolyCutoff2s"
      pl = pu = D["cutoff"]["P"]
   else
      error("unknown cutoff %(cutoff_id)")
   end
   if tl > tu
      tl, tr = tu, tl
      pl, pr = pu, pl
   else
      tl, tr = tl, tu
      pl, pr = pl, pu
   end

   A = Float64.(D["jacobicoeffs"]["A"])
   B = Float64.(D["jacobicoeffs"]["B"])
   C = Float64.(D["jacobicoeffs"]["C"])
   nrm = Float64.(D["jacobicoeffs"]["nrm"])
   deg = D["deg"]
   α, β = D["a"], D["b"]

   An = zeros(length(A) + 1)
   Bn = zeros(length(An))
   Cn = zeros(length(An))
   An[1] = nrm[1]
   An[2] = -nrm[2] / nrm[1] * 0.5 * (α + β + 2)
   Bn[2] = nrm[2] / nrm[1] * ((α + 1) - 0.5 * (α + β + 2))
   for n = 3:length(An)
      An[n] = -nrm[n] / nrm[n-1] * A[n-1]
      Bn[n] = nrm[n] / nrm[n-1] * B[n-1]
      Cn[n] = nrm[n] / nrm[n-2] * C[n-1]
   end

   # transform to the correct interval
   a = (2 / (tr - tl))
   At = An * a
   Bt = Bn - An * ((tr + tl) / (tr - tl))
   Ct = Cn
   At[1] = An[1] * a^(pl+pr)

   Jt = ACE.OrthPolys.OrthPolyBasis(pl, tl, pr, tr, At, Bt, Ct,
                                      Float64[], Float64[])
   Jr = ACE.OrthPolys.TransformedPolys(Jt, trans, rl, ru)

   if !isempty(rtests) > 0
      r = rtests[1]["r"]
      Jr_test = rtests[1]["Pr"]
      Jr_new = evaluate(Jr, r)
      f = Jr_test[1] / Jr_new[1]
      Jr.J.A[1] *= f
   end
   return Jr
end

import_pipot_v05(fname::AbstractString) =
      import_pipot_v05(JuLIP.read_dict(fname))

function import_pipot_v05(D::Dict)
   species = D["zlist"]["list"]
   zlist = read_dict(D["zlist"])
   rawcoeffs = Float64.(D["coeffs_re"][1])

   # get the radial basis
   Pr = import_rbasis_v05(D["J"], D["rtests"])

   # get the 1-p and n-p basis specifications
   rawspec = D["aalists"][1]["ZKLM_list"]

   pispec = []
   pispec_i = []
   pispec_i_dict = Dict{Vector{Int}, Int}()
   spec1 = []
   coeffs = []
   for (idx, ZNLM) in enumerate(rawspec)
      izs = ZNLM[1]
      ns = ZNLM[2] .+ 1
      ls = ZNLM[3]
      ms = ZNLM[4]
      order = length(ns)
      zs = [AtomicNumber(species[iz]) for iz in izs]
      bs = RnYlmBasisFcn[]
      bs_i = Int[]
      for α = 1:order
         b = RnYlmBasisFcn(ns[α], ls[α], ms[α], zs[α])
         I1 = findall(isequal(b), spec1)
         if isempty(I1)
            push!(spec1, b)
            i1 = length(spec1)
         else
            @assert length(I1) == 1
            i1 = I1[1]
         end
         push!(bs, b)
         push!(bs_i, i1)
      end
      p = sortperm(bs_i)
      bs_i = bs_i[p]
      bb = PIBasisFcn(AtomicNumber(species[1]), bs[p])

      if haskey(pispec_i_dict, bs_i)
         In = pispec_i_dict[bs_i]
         coeffs[In] += rawcoeffs[idx]
      else
         push!(pispec, bb)
         push!(pispec_i, bs_i)
         push!(coeffs, rawcoeffs[idx])
         pispec_i_dict[bs_i] = length(pispec_i)
      end

   end

   # construct 1p basis from the specs
   basis1p = RnYlm1pBasis(Pr, zlist, identity.(spec1))
   # build a naive index allocation (we have only one species!!)
   AAindices = 1:length(pispec)
   # construct the inner PI basis from the specification
   inner = ACE.InnerPIBasis(spec1, pispec, AAindices, zlist.list[1])
   # ... and finally put it all together
   pibasis = PIBasis(basis1p, zlist, (inner,), DAGEvaluator())

   # now the inner basis has sorted the basis functins so we need to
   # fix the permutation of the coefficients
   perm = [ inner.b2iAA[b]  for b in pispec ]
   invperm = sortperm(perm)
   V = PIPotential(pibasis, (Float64.(coeffs[invperm]),))

   return V
end

end


text:
exclude(f, ext) = occursin(ext, f) || endswith(f, ".opj")

function replace_unicode(str)
    new_string = ""
    for letter in str
        if !isascii(letter)
            c = string(codepoint(letter), base=16, pad=4)
            new_string *= "\\x($c)"
        else
            new_string *= letter
        end
    end
    return new_string
end

function replace_powers(str)
    inv_one = r"(?<pre>[a-zA-Z]+) (?<unit>[a-zA-Z]+)\^-1"
    inv_pre = r"(?<pre>[a-zA-Z]+) (?<unit>[a-zA-Z]+)\^-(?<power>[2-9]+)"
    power_only_unit = r"(?<unit>[a-zA-Z]+)\^(?<power>-?[1-9]+)"
    str = replace(str, inv_one=>s"\g<pre>/\g<unit>")
    str = replace(str, inv_pre=>s"\g<pre>/\g<unit>\\+(\g<power>)")
    str = replace(str, power_only_unit => s"\g<unit>\\+(\g<power>)")

    return str
end

function to_origin(str)
    str = replace_powers(str)
    replace_unicode(str)
end

function comment_value(datafile)
    u = datafile.legend_units
    idx = datafile.idx

    value = filevalue(datafile)
    unicode_val = si_round(uparse(value*u), idx)
    return to_origin(unicode_val)
end


text:
using Swagger, Test

function as_taskfailedexception(ex)
    try
        task = @async throw(ex)
        wait(task)
    catch ex
        return ex
    end
end

function test_longpoll_exception_check()
    resp = Swagger.Downloads.Response("http", "http://localhost", 200, "no error", [])
    reqerr1 = Swagger.Downloads.RequestError("http://localhost", 500, "not timeout error", resp)
    reqerr2 = Swagger.Downloads.RequestError("http://localhost", 200, "Operation timed out after 300 milliseconds with 0 bytes received", resp) # timeout error

    @test Swagger.is_longpoll_timeout("not an exception") == false

    swaggerex1 = Swagger.ApiException(reqerr1)
    @test Swagger.is_longpoll_timeout(swaggerex1) == false
    @test Swagger.is_longpoll_timeout(as_taskfailedexception(swaggerex1)) == false

    swaggerex2 = Swagger.ApiException(reqerr2)
    @test Swagger.is_longpoll_timeout(swaggerex2)
    @test Swagger.is_longpoll_timeout(as_taskfailedexception(swaggerex2))

    @test Swagger.is_longpoll_timeout(CompositeException([swaggerex1, swaggerex2]))
    @test Swagger.is_longpoll_timeout(CompositeException([swaggerex1, as_taskfailedexception(swaggerex2)]))
    @test Swagger.is_longpoll_timeout(CompositeException([swaggerex1, as_taskfailedexception(swaggerex1)])) == false
end


text:
######################################################################
# FactCheck.jl
# A testing framework for Julia
# http://github.com/JuliaLang/FactCheck.jl
# MIT Licensed
######################################################################

module FactCheck

export @fact, @fact_throws, @pending,
       facts, context,
       # Assertion helpers
       not,
       exactly,
       roughly,
       anyof,
       less_than, less_than_or_equal,
       greater_than, greater_than_or_equal

######################################################################
# Success, Failure, Error <: Result
# Represents the result of a test. These are very similar to the types
# with the same names in Base.Test, except for the addition of the
# `ResultMetadata` type that is used to retain information about the test,
# such as its file, line number, description, etc.
abstract type Result end

mutable struct ResultMetadata
    line
    msg
    function ResultMetadata(;line=nothing, msg=nothing)
        new(line, msg)
    end
end

mutable struct Success <: Result
    expr::Expr
    fact_type::Symbol
    lhs  # What it was
    rhs  # What it should have been
    meta::ResultMetadata
end

mutable struct Failure <: Result
    expr::Expr
    fact_type::Symbol
    lhs  # What it was
    rhs  # What it should have been
    meta::ResultMetadata
end

mutable struct Error <: Result
    expr::Expr
    fact_type::Symbol
    err::Exception
    backtrace
    meta::ResultMetadata
end

struct Pending <: Result end

# Collection of all results across facts
allresults = Result[]
clear_results() = (global allresults; allresults = Result[])

# Formats a fact expression
function format_fact(ex::Expr)
    if ex.head == :(-->) || ex.head == :(=>)
        # :(fn(1) --> 2) to 'fn(1) --> 2'
        # :("1"*"1" --> "11") to '"1" * "1" --> "11"'
        # We handle non-expresion arguments differently,
        # otherwise, e.g. quote marks on strings disappear
        x, y = ex.args
        x_str = sprint(isa(x,Expr) || isa(x,Symbol) ? print : show, x)
        y_str = sprint(isa(y,Expr) || isa(y,Symbol) ? print : show, y)
        string(x_str, " --> ", y_str)
    else
        # Something else, that maybe didn't have a -->
        # such as @fact_throws. Punt and just stringify
        string(ex)
    end
end

# Define printing functions for the result types
function Base.show(io::IO, f::Failure)
    println(io, "\n<IT::>", f.meta.msg != nothing ? "$(f.meta.msg)" : format_fact(f.expr))
    println(io, "\n<FAILED::>Test Failed")
    if f.fact_type == :fact_throws
        # @fact_throws didn't get an error, or the right type of error
        if f.rhs != :fact_throws_noerror
            println(io, "\tExpected: ", f.rhs[1])
            println(io, "\tOccurred: ", f.rhs[2])
        end
    elseif f.fact_type == :fact
        # @fact didn't get the right result
        args = f.expr.args
        if length(args) >= 2 && _factcheck_function(args[2]) != nothing
            # Fancy helper fact
            fcFunc = _factcheck_function(args[2])
            if haskey(FACTCHECK_FUN_NAMES, fcFunc)
                println(io, replace_lf(string("\tExpected: ", sprint(show, f.lhs), " ", FACTCHECK_FUN_NAMES[fcFunc], " ", sprint(show, f.rhs))))
            else
                println(io, replace_lf(string("\tExpected: ", sprint(show, f.lhs), " --> ", fcFunc, "(", sprint(show, f.rhs), ")")))
            end
        else
            # Normal equality-test-style fact
            println(io, replace_lf(string("\tExpected: ", sprint(show, f.rhs), "\n", "\tOccurred: ", sprint(show, f.lhs))))
        end
    else
        error("Unknown fact type: ", f.fact_type)
    end
    println(io, "\n<COMPLETEDIN::>")
end

function Base.show(io::IO, e::Error)
    println(io, "\n<IT::>", e.meta.msg != nothing ? "$(e.meta.msg)" : format_fact(e.expr))
    println(io, "\n<ERROR::>Test Errored")
    println(io, replace_lf(sprint(showerror, e.err, e.backtrace)))
    println(io, "\n<COMPLETEDIN::>")
end

function Base.show(io::IO, s::Success)
    if s.rhs == :fact_throws_error
        println(io, "\n<IT::>", s.meta.msg != nothing ? "$(s.meta.msg)" : "Throws Error")
        println(io, "\n<PASSED::>Test Passed")
        println(io, "\n<COMPLETEDIN::>")
    else
        println(io, "\n<IT::>", s.meta.msg != nothing ? "$(s.meta.msg)" : format_fact(s.expr))
        println(io, "\n<PASSED::>Test Passed")
        println(io, "\n<COMPLETEDIN::>")
    end
end

function Base.show(io::IO, ::Pending)
    println(io, "\n<LOG::>Pending")
end

const SPECIAL_FACTCHECK_FUNCTIONS =
    Set([:not, :exactly, :roughly, :anyof,
         :less_than, :less_than_or_equal, :greater_than, :greater_than_or_equal])

const FACTCHECK_FUN_NAMES =
    Dict{Symbol,AbstractString}(
      :roughly => "≅",
      :less_than => "<",
      :less_than_or_equal => "≤",
      :greater_than => ">",
      :greater_than_or_equal => "≥")

isexpr(x) = isa(x, Expr)
iscallexpr(x) = isexpr(x) && x.head == :call
isdotexpr(x) = isexpr(x) && x.head == :.
isquoteexpr(x) = isexpr(x) && x.head == :quote
isparametersexpr(x) = isexpr(x) && x.head == :parameters

function _factcheck_function(assertion)
    iscallexpr(assertion) || return nothing

    # checking for lhs => roughly(rhs)
    if assertion.args[1] in SPECIAL_FACTCHECK_FUNCTIONS
        return assertion.args[1]
    end

    # checking for lhs => FactCheck.roughly(rhs)
    isdotexpr(assertion.args[1]) || return nothing
    dotexpr = assertion.args[1]
    length(dotexpr.args) >= 2 || return nothing
    if isquoteexpr(dotexpr.args[2])
        quoteexpr = dotexpr.args[2]
        if length(quoteexpr.args) >= 1 && quoteexpr.args[1] in SPECIAL_FACTCHECK_FUNCTIONS
            return quoteexpr.args[1]
        else
            return nothing
        end
    end

    # sometimes it shows up as a QuoteNode...
    if isa(dotexpr.args[2], QuoteNode) && dotexpr.args[2].value in SPECIAL_FACTCHECK_FUNCTIONS
        return dotexpr.args[2].value
    end
    nothing
end

######################################################################
# Core testing macros and functions

# @fact takes an assertion of the form LHS --> RHS, and replaces it
# with code to evaluate that fact (depending on the type of the RHS),
# and produce and record a result based on the outcome
macro fact(factex::Expr, args...)
    if factex.head != :(-->) && factex.head != :(=>)
        error("Incorrect usage of @fact: $factex")
    end
    # TODO: remove deprecated syntax
    if factex.head == :(=>)
        Base.warn_once("The `=>` syntax is deprecated, use `-->` instead")
    end
    # Extract the two sides of the fact
    lhs, initial_rhs = factex.args
    # If there is another argument to the macro, assume it is a
    # message and record it
    msg = length(args) > 0 ? args[1] : (:nothing)

    # rhs is the assertion, unless it's wrapped by a special FactCheck function
    rhs = initial_rhs
    if _factcheck_function(initial_rhs) != nothing
        rhs = initial_rhs.args[isparametersexpr(initial_rhs.args[2]) ? 3 : 2]
    end

    quote
        # Build a function (predicate) that, depending on the nature of
        # the RHS, either compares the sides or applies the RHS to the LHS
        predicate = function(lhs_value)
            rhs_value = $(esc(initial_rhs))
            if isa(rhs_value, Function)
                # The RHS is a function, so instead of testing for equality,
                # return the value of applying the RHS to the LHS
                (rhs_value(lhs_value), lhs_value, $(esc(rhs)))
            else
                # The RHS is a value, so test for equality
                (rhs_value == lhs_value, lhs_value, $(esc(rhs)))
            end
        end
        # Replace @fact with a call to the do_fact function that constructs
        # the test result object by evaluating the
        do_fact(() -> predicate($(esc(lhs))),
                $(Expr(:quote, factex)),
                :fact,
                ResultMetadata(line=getline(),
                               msg=$(esc(msg))))
    end
end

# `@fact_throws` is similar to `@fact`, except it only checks if
# the expression throws an error or not - there is no explict
# assertion to compare against.
macro fact_throws(args...)
    expr, extype, msg = nothing, nothing, nothing
    nargs = length(args)
    if nargs == 1
        if isa(args[1],Expr)
            expr = args[1]
        else
            throw(ArgumentError("invalid @fact_throws macro"))
        end
    elseif nargs == 2
        if (isa(args[1],Symbol) || isa(args[1],Expr)) && isa(args[2],Expr)
            extype, expr = args
        elseif isa(args[1],Expr)
            expr, msg = args
        else
            throw(ArgumentError("invalid @fact_throws macro"))
        end
    elseif nargs >= 3
        if (isa(args[1],Symbol) || isa(args[1], Expr)) && isa(args[2],Expr)
            extype, expr, msg = args
        else
            throw(ArgumentError("invalid @fact_throws macro"))
        end
    end

    quote
        do_fact(() -> try
                          $(esc(expr))
                          (false, "no exception was thrown", :fact_throws_noerror)
                      catch ex
                          $(if extype === nothing
                              :((true, "an exception was thrown", :fact_throws_error))
                            else
                              :(if isa(ex,$(esc(extype)))
                                  (true, "correct exception was thrown", :fact_throws_error)
                                else
                                  (false, "wrong exception was thrown",
                                    ($(esc(extype)), typeof(ex)))
                                end)
                            end)
                      end,
                $(Expr(:quote, expr)),
                :fact_throws,
                ResultMetadata(line=getline(),msg=$(esc(msg))))
    end
end

# `do_fact` constructs a Success, Failure, or Error depending on the
# outcome of a test and passes it off to the active test handler
# `FactCheck.handlers[end]`. It finally returns the test result.
function do_fact(thunk::Function, factex::Expr, fact_type::Symbol, meta::ResultMetadata)
    result = try
        res, val, rhs = thunk()
        res ? Success(factex, fact_type, val, rhs, meta) :
                Failure(factex, fact_type, val, rhs, meta)
    catch err
        Error(factex, fact_type, err, catch_backtrace(), meta)
    end

    !isempty(handlers) && handlers[end](result)
    result
end

# `@pending` is a no-op test - it doesn't do anything except record
# its existance in the final totals of tests "run"
macro pending(factex::Expr, args...)
    quote
        result = Pending()
        !isempty(handlers) && handlers[end](result)
        result
    end
end

######################################################################
# Grouping of tests
#
# `facts` describes a top-level test scope, which can contain
# `contexts` to group similar tests. Test results will be collected
# instead of throwing an exception immediately.

# A TestSuite collects the results of a series of tests, as well as
# some information about the tests such as their file and description.
mutable struct TestSuite
    filename
    desc
    successes::Vector{Success}
    failures::Vector{Failure}
    errors::Vector{Error}
    pending::Vector{Pending}
end
TestSuite(f, d) = TestSuite(f, d, Success[], Failure[], Error[], Pending[])

# The last handler function found in `handlers` will be passed
# test results.
const handlers = Function[]

# A list of test contexts. `contexts[end]` should be the
# inner-most context.
const contexts = AbstractString[]

# Constructs a function that handles Successes, Failures, and Errors,
# pushing them into a given TestSuite and printing Failures and Errors
# as they arrive (unless in compact mode, in which case we delay
# printing details until the end).
function make_handler(suite::TestSuite)
    function delayed_handler(r::Success)
        print(r)
    end
    function delayed_handler(r::Failure)
        print(r)
    end
    function delayed_handler(r::Error)
        print(r)
    end
    function delayed_handler(r::Pending)
        print(r)
    end
    delayed_handler
end

# facts
# Creates testing scope. It is responsible for setting up a testing
# environment, which means constructing a `TestSuite`, generating
# and registering test handlers, and reporting results.
function facts(f::Function, desc)
    println(string("\n<DESCRIBE::>", desc == nothing ? "Unnamed Facts" : "$(desc)"))
    handler = make_handler(TestSuite(nothing, desc))
    push!(handlers, handler)
    f()
    pop!(handlers)
    println("\n<COMPLETEDIN::>")
end
facts(f::Function) = facts(f, nothing)

# context
# Executes a battery of tests in some descriptive context, intended
# for use inside of `facts`. Displays the string in default mode.
function context(f::Function, desc::AbstractString)
    println(string("\n<DESCRIBE::>", desc == "" ? "Unnamed Context" : desc))
    push!(contexts, desc)
    try
        f()
    finally
        pop!(contexts)
        println("\n<COMPLETEDIN::>")
    end
end
context(f::Function) = context(f, "")

######################################################################

@noinline getline() = StackTraces.stacktrace()[2].line

replace_lf(s::AbstractString) = replace(s, "\n" => "<:LF:>")

############################################################
# Assertion helpers
include("helpers.jl")

end # module FactCheck