text:
# Utility functions for processing VLSV data.

"""
    getcell(meta, location) -> UInt

Return cell ID containing the given spatial `location` in meter, excluding domain
boundaries. Only accept 3D location.
"""
function getcell(meta::MetaVLSV, loc)
   (;coordmin, coordmax, dcoord, ncells, cellid, maxamr) = meta

   foreach( (i,comp) -> coordmin[i] < loc[i] < coordmax[i] ? nothing :
      error("$comp coordinate out of bound!"), 1:3, 'x':'z')

   # Get cell indices
   indices = @inbounds ntuple(i -> round(UInt, (loc[i] - coordmin[i]) ÷ dcoord[i]), Val(3))
   # Get cell id
   cid = @inbounds indices[1] + indices[2]*ncells[1] + indices[3]*ncells[1]*ncells[2] + 1

   ncells_lowerlevel = UInt(0)
   ncell = prod(ncells)

   @inbounds for ilevel = 0:maxamr
      cid in cellid && break

      ncells_lowerlevel += 2^(3*ilevel)*ncell

      ratio = 2^(ilevel+1)

      indices = ntuple(i -> floor(UInt, (loc[i] - coordmin[i]) / dcoord[i] * ratio), Val(3))

      cid = ncells_lowerlevel + indices[1] +
         ratio*ncells[1]*indices[2] + ratio^2*ncells[1]*ncells[2]*indices[3] + 1
   end

   cid
end

"""
    getlevel(meta, cid) -> Int

Return the AMR level of a given cell ID `cid`.
!!! warning
    This function does not check if the VLSV file of `meta` actually contains `cid`; it may
    be shadowed by refined children.
"""
function getlevel(meta::MetaVLSV, cid::Integer)
   ncell = prod(meta.ncells)
   ilevel = 0
   c = Int(cid) - ncell
   while c > 0
      ilevel += 1
      c -= 2^(3*ilevel)*ncell
   end
   ilevel
end

"""
    getparent(meta, cid) -> Int

Return the parent cell ID of given child `cid`.
"""
function getparent(meta::MetaVLSV, cid::Integer)
   @inbounds xcell, ycell = meta.ncells[1], meta.ncells[2]
   ncell = prod(meta.ncells)

   mylvl = getlevel(meta, cid)
   parentlvl = mylvl - 1

   if parentlvl < 0
      throw(ArgumentError("Cell ID $cid has no parent cell!"))
   else
      # get the first cellid on my level
      cid1st = get1stcell(mylvl, ncell) + 1
      # get row and column sequence on my level (starting with 0)
      xcell <<= mylvl
      ycell <<= mylvl

      myseq = cid - cid1st
      ix = myseq % xcell
      iz = myseq ÷ (xcell*ycell)
      iy = (myseq - iz*xcell*ycell) ÷ xcell
      # indexes on the parent level
      ixparent = ix ÷ 2
      iyparent = iy ÷ 2
      izparent = iz ÷ 2

      # get the first cellid on parent level
      cid1st -= ncell*8^parentlvl
      # get parent cellid (may not exist!!!)
      parentid = cid1st + izparent*xcell*ycell÷4 + iyparent*xcell÷2 + ixparent
   end
   parentid
end

"""
    getchildren(meta, cid) -> Vector{Int}

Return direct children of `cid`.
"""
function getchildren(meta::MetaVLSV, cid::Integer)
   xcell, ycell, zcell = meta.ncells
   ncell = prod(meta.ncells)

   mylvl = getlevel(meta, cid)

   # get the first cell ID on the my level
   cid1st = 1
   for i = 0:mylvl-1
      cid1st += ncell * 8^i
   end
   # get my row and column sequence on my level (starting with 0)
   xcell <<= mylvl
   ycell <<= mylvl

   myseq = cid - cid1st
   ix = myseq % xcell
   iz = myseq ÷ (xcell*ycell)
   iy = (myseq - iz*xcell*ycell) ÷ xcell

   # get the children sequences on the finer level
   ix *= 2
   iy *= 2
   iz *= 2

   nchildren = 2^ndims(meta)
   cid = zeros(Int, nchildren)
   # get the first cell ID on the finer level
   cid1st += ncell*8^mylvl
   ix_, iy_ = (ix, ix+1), (iy, iy+1)
   iz_ = zcell != 1 ? (iz, iz+1) : iz
   for (n,i) in enumerate(Iterators.product(ix_, iy_, iz_))
      @inbounds cid[n] = cid1st + i[3]*xcell*ycell*4 + i[2]*xcell*2 + i[1]
   end
   (cid)
end

"""
    getsiblings(meta, cid) -> Vector{Int}

Return sibling cells of a given `cid`, including itself.
"""
function getsiblings(meta::MetaVLSV, cid::Integer)
   xcell, ycell, zcell = meta.ncells
   ncell = prod(meta.ncells)

   mylvl = getlevel(meta, cid)

   mylvl == 0 && throw(ArgumentError("CellID $cid is not a child cell!"))

   xcell = xcell << mylvl
   ycell = ycell << mylvl

   # 1st cellid on my level
   cid1st = get1stcell(mylvl, ncell) + 1

   # xyz sequences on my level (starting with 0)
   myseq = cid - cid1st
   ix = myseq % xcell
   iz = myseq ÷ (xcell*ycell)
   iy = (myseq - iz*xcell*ycell) ÷ xcell

   ix1 = iseven(ix) ? ix + 1 : ix - 1
   iy1 = iseven(iy) ? iy + 1 : iy - 1
   iz1 = iseven(iz) ? iz + 1 : iz - 1
   # reorder
   ix, ix1 = minmax(ix, ix1)
   iy, iy1 = minmax(iy, iy1)
   iz, iz1 = minmax(iz, iz1)

   nsiblings = 2^ndims(meta)
   cid = zeros(Int, nsiblings)
   ix_, iy_ = (ix, ix1), (iy, iy1)
   iz_ = zcell != 1 ? (iz, iz1) : iz
   for (n,i) in enumerate(Iterators.product(ix_, iy_, iz_))
      @inbounds cid[n] = cid1st + i[3]*xcell*ycell + i[2]*xcell + i[1]
   end
   (cid)
end

"""
    isparent(meta, cid) -> Bool

Check if `cid` is a parent cell.
"""
function isparent(meta::MetaVLSV, cid::Integer)
   ncell_accum = get1stcell(meta.maxamr, prod(meta.ncells))

   cid ∉ meta.cellid && 0 < cid ≤ ncell_accum
end

"""
    getcellcoordinates(meta, cid) -> NTuple{3, Float64}

Return a given cell's spatial coordinates.
"""
function getcellcoordinates(meta::MetaVLSV, cid::Integer)
   (;ncells, coordmin, coordmax) = meta
   cid -= 1 # for easy divisions

   ncells_refmax = collect(ncells)
   reflevel = 0
   subtraction = prod(ncells) * (2^reflevel)^3
   # sizes on the finest level
   while cid ≥ subtraction
      cid -= subtraction
      reflevel += 1
      subtraction *= 8
      ncells_refmax[1] *= 2
      ncells_refmax[2] *= 2
      ncells_refmax[3] *= 2
   end

   indices = @inbounds (
      cid % ncells_refmax[1],
      cid ÷ ncells_refmax[1] % ncells_refmax[2],
      cid ÷ (ncells_refmax[1] * ncells_refmax[2]) )

   coords = @inbounds ntuple( i ->
      coordmin[i] + (indices[i] + 0.5) * (coordmax[i] - coordmin[i]) / ncells_refmax[i],
      Val(3))

   coords
end

"""
    getvcellcoordinates(meta, vcellids; species="proton")

Return velocity cells' coordinates of `species` and `vcellids`.
"""
function getvcellcoordinates(meta::MetaVLSV, vcellids; species="proton")
   (;vblocks, vblock_size, dv, vmin) = meta.meshes[species]

   bsize = prod(vblock_size)
   blockid = @. vcellids ÷ bsize
   # Get block coordinates
   blockInd = [(
      bid % vblocks[1],
      bid ÷ vblocks[1] % vblocks[2],
      bid ÷ (vblocks[1] * vblocks[2]) )
      for bid in blockid]
   blockCoord = [(
      bInd[1] * dv[1] * vblock_size[1] + vmin[1],
      bInd[2] * dv[2] * vblock_size[2] + vmin[2],
      bInd[3] * dv[3] * vblock_size[3] + vmin[3] )
      for bInd in blockInd]

   # Get cell indices
   vcellblockids = @. vcellids % bsize
   cellidxyz = [(
      cid % vblock_size[1],
      cid ÷ vblock_size[1] % vblock_size[2],
      cid ÷ (vblock_size[1] * vblock_size[2]) )
      for cid in vcellblockids]

   # Get cell coordinates
   cellCoords = [SVector(0.0f0, 0.0f0, 0.0f0) for _ in vcellblockids]
   @inbounds @simd for i in eachindex(vcellblockids)
      cellCoords[i] = [blockCoord[i][j] + (cellidxyz[i][j] + 0.5) * dv[j] for j in 1:3]
   end

   cellCoords
end

"""
    getdensity(meta, VDF; species="proton")
    getdensity(meta, vcellf; species="proton")
    getdensity(vmesh::VMeshInfo, vcellf)

Get density from `VDF` of `species` associated with `meta`, n = ∫ f(r,v) dV. Alternatively,
one can directly pass `vcellids` as original indices of nonzero VDFs and `vcellf` as their
corresponding values.
"""
function getdensity(meta::MetaVLSV, VDF::Array{T};
   species="proton") where T <: AbstractFloat

   (;dv) = meta.meshes[species]
   n = sum(VDF) * convert(T, prod(dv))
end

function getdensity(vmesh::VMeshInfo, vcellf::Vector{T}) where T <: AbstractFloat
   n = sum(vcellf) * convert(T, prod(vmesh.dv))
end

getdensity(meta::MetaVLSV, vcellf; species="proton") =
   getdensity(meta.meshes[species], vcellf)

"""
    getvelocity(meta, VDF; species="proton")
    getvelocity(meta, vcellids, vcellf; species="proton")
    getvelocity(vmesh::VMeshInfo, vcellids, vcellf)

Get bulk velocity from `VDF` of `species`, u = ∫ v * f(r,v) dV / n. Alternatively, one can
directly pass `vcellids`, `vcellf`, as in [`getdensity`](@ref).
"""
function getvelocity(meta::MetaVLSV, VDF::Array{T};
   species="proton") where T <: AbstractFloat

   (;dv, vmin) = meta.meshes[species]
   u = zeros(T, 3)

   @inbounds for k in axes(VDF,3), j in axes(VDF,2), i in axes(VDF,1)
      vx = vmin[1] + (i - 0.5f0)*dv[1]
      vy = vmin[2] + (j - 0.5f0)*dv[2]
      vz = vmin[3] + (k - 0.5f0)*dv[3]
      u[1] += vx*VDF[i,j,k]
      u[2] += vy*VDF[i,j,k]
      u[3] += vz*VDF[i,j,k]
   end

   n = sum(VDF)

   u ./= n
   SVector{3}(u)
end

function getvelocity(vmesh::VMeshInfo, vcellids::Vector{UInt32}, vcellf::Vector{T}) where
   T <: AbstractFloat

   (;vblock_size, vblocks, dv, vmin) = vmesh
   vsize = @inbounds ntuple(i -> vblock_size[i] * vblocks[i], Val(3))
   slicez = vsize[1]*vsize[2]
   blocksize = prod(vblock_size)
   sliceBz = vblocks[1]*vblocks[2]
   sliceCz = vblock_size[1]*vblock_size[2]

   u = zeros(T, 3)

   @inbounds @simd for ic in eachindex(vcellids)
      id = findindex(vcellids[ic], vblocks, vblock_size, blocksize, vsize, sliceBz, sliceCz)
      i = id % vsize[1]
      j = id % slicez ÷ vsize[1]
      k = id ÷ slicez

      vx = vmin[1] + (i + 0.5f0)*dv[1]
      vy = vmin[2] + (j + 0.5f0)*dv[2]
      vz = vmin[3] + (k + 0.5f0)*dv[3]
      u[1] += vx*vcellf[ic]
      u[2] += vy*vcellf[ic]
      u[3] += vz*vcellf[ic]
   end

   n = sum(vcellf)

   u ./= n
   SVector{3}(u)
end

getvelocity(meta::MetaVLSV, vcellids, vcellf; species="proton") =
   getvelocity(meta.meshes[species], vcellids, vcellf)

"""
    getpressure(meta, VDF; species="proton")
    getpressure(meta, vcellids, vcellf; species="proton")
    getpressure(vmesh::VMeshInfo, vcellids, vcellf)

Get pressure tensor (6 components) of `species` from `VDF` associated with `meta`,
pᵢⱼ = m/3 * ∫ (v - u)ᵢ(v - u)ⱼ * f(r,v) dV. Alternatively, one can directly pass `vcellids`,
`vcellf`, as in [`getdensity`](@ref).
"""
function getpressure(meta::MetaVLSV, VDF::Array{T};
   species="proton") where T <: AbstractFloat

   (;dv, vmin) = meta.meshes[species]
   p = zeros(T, 6)

   u = getvelocity(meta, VDF; species)

   @inbounds for k in axes(VDF,3), j in axes(VDF,2), i in axes(VDF,1)
      vx = vmin[1] + (i - 0.5f0)*dv[1]
      vy = vmin[2] + (j - 0.5f0)*dv[2]
      vz = vmin[3] + (k - 0.5f0)*dv[3]

      p[1] += (vx - u[1])*(vx - u[1])*VDF[i,j,k]
      p[2] += (vy - u[2])*(vy - u[2])*VDF[i,j,k]
      p[3] += (vz - u[3])*(vz - u[3])*VDF[i,j,k]
      p[4] += (vy - u[2])*(vz - u[3])*VDF[i,j,k]
      p[5] += (vx - u[1])*(vz - u[3])*VDF[i,j,k]
      p[6] += (vx - u[1])*(vy - u[2])*VDF[i,j,k]
   end

   factor = mᵢ * convert(T, prod(dv))
   p .*= factor
   SVector{6}(p)
end

function getpressure(vmesh::VMeshInfo, vcellids::Vector{UInt32}, vcellf::Vector{T}) where
   T <: AbstractFloat

   (;vblock_size, vblocks, dv, vmin) = vmesh
   vsize = @inbounds ntuple(i -> vblock_size[i] * vblocks[i], Val(3))
   slicez = vsize[1]*vsize[2]
   blocksize = prod(vblock_size)
   sliceBz = vblocks[1]*vblocks[2]
   sliceCz = vblock_size[1]*vblock_size[2]

   u = getvelocity(vmesh, vcellids, vcellf)

   p = zeros(T, 6)

   @inbounds @simd for ic in eachindex(vcellids)
      id = findindex(vcellids[ic], vblocks, vblock_size, blocksize, vsize, sliceBz, sliceCz)
      i = id % vsize[1]
      j = id % slicez ÷ vsize[1]
      k = id ÷ slicez

      vx = vmin[1] + (i + 0.5f0)*dv[1]
      vy = vmin[2] + (j + 0.5f0)*dv[2]
      vz = vmin[3] + (k + 0.5f0)*dv[3]
      p[1] += (vx - u[1])*(vx - u[1])*vcellf[ic]
      p[2] += (vy - u[2])*(vy - u[2])*vcellf[ic]
      p[3] += (vz - u[3])*(vz - u[3])*vcellf[ic]
      p[4] += (vy - u[2])*(vz - u[3])*vcellf[ic]
      p[5] += (vx - u[1])*(vz - u[3])*vcellf[ic]
      p[6] += (vx - u[1])*(vy - u[2])*vcellf[ic]
   end

   factor = mᵢ * convert(T, prod(dv))
   p .*= factor
   SVector{6}(p)
end

getpressure(meta::MetaVLSV, vcellids, vcellf; species="proton") =
   getpressure(meta.meshes[species], vcellids, vcellf)

"Get the original vcell index without blocks from raw vcell index `i` (0-based)."
@inline function findindex(i, vblocks, vblock_size, blocksize, vsize, sliceBz, sliceCz)
   iB = i ÷ blocksize
   iBx = iB % vblocks[1]
   iBy = iB % sliceBz ÷ vblocks[1]
   iBz = iB ÷ sliceBz
   iCellInBlock = i % blocksize
   iCx = iCellInBlock % vblock_size[1]
   iCy = iCellInBlock % sliceCz ÷ vblock_size[1]
   iCz = iCellInBlock ÷ sliceCz
   iBCx = iBx*vblock_size[1] + iCx
   iBCy = iBy*vblock_size[2] + iCy
   iBCz = iBz*vblock_size[3] + iCz
   iOrigin = iBCz*vsize[1]*vsize[2] + iBCy*vsize[1] + iBCx
end

"""
    reorder(vmesh::VMeshInfo, vcellids) -> vcellids_origin

Reorder vblock-organized VDF indexes into x-->y-->z indexes. `vcellids` are raw indices
of nonzero VDFs ordered by blocks.
"""
function reorder(vmesh::VMeshInfo, vcellids)
   (;vblock_size, vblocks) = vmesh
   blocksize = prod(vblock_size)
   sliceBz = vblocks[1]*vblocks[2]
   vsize = @inbounds ntuple(i -> vblock_size[i] * vblocks[i], Val(3))
   sliceCz = vblock_size[1]*vblock_size[2]

   vcellids_origin = similar(vcellids)
   # IDs are 0-based
   @inbounds @simd for i in eachindex(vcellids)
      vcellids_origin[i] = 1 +
         findindex(vcellids[i], vblocks, vblock_size, blocksize, vsize, sliceBz, sliceCz)
   end

   vcellids_origin
end

"""
    reconstruct(vmesh::VMeshInfo, vcellids, vcellf)

Reconstruct the full VDFs in 3D. `vcellids` are raw indices of nonzero VDFs ordered by
blocks, and `vcellf` are the corresponding values in each cell.
"""
function reconstruct(vmesh::VMeshInfo, vcellids, vcellf)
   (;vblock_size, vblocks) = vmesh
   blocksize = prod(vblock_size)
   sliceBz = vblocks[1]*vblocks[2]
   vsize = @inbounds ntuple(i -> vblock_size[i] * vblocks[i], Val(3))
   sliceCz = vblock_size[1]*vblock_size[2]
   # Reconstruct the full velocity space
   VDF = zeros(Float32, vsize)
   # Raw IDs are 0-based
   @inbounds @simd for i in eachindex(vcellids)
      j = 1 +
         findindex(vcellids[i], vblocks, vblock_size, blocksize, vsize, sliceBz, sliceCz)
      VDF[j] = vcellf[i]
   end

   VDF
end

"""
    getmaxwellianity(meta, VDF; species="proton")
    getmaxwellianity(meta, vcellids, vcellf; species="proton")

Obtain the Maxwellian similarity factor -log(1/(2n) * ∫ |f - g| dv), where `f` is the VDF
from Vlasiator and `g` is the analytical Maxwellian distribution that generates the same
density as `f`. The value ranges from [0, +∞], with 0 meaning not Maxwellian-distributed at
all, and +∞ a perfect Maxwellian distribution.
Alternatively, one can pass original `vcellids` and `vcellf` directly.
"""
function getmaxwellianity(meta, VDF; species="proton")
   (;dv, vmin) = meta.meshes[species]

   n = getdensity(meta, VDF)
   u = getvelocity(meta, VDF)
   P = getpressure(meta, VDF)
   p = (P[1] + P[2] + P[3]) / 3
   T = p / (n *kB) # temperature from scalar pressure

   ϵₘ = zero(eltype(VDF))
   vth2Inv = mᵢ / (2kB*T)

   @inbounds for k in axes(VDF,3), j in axes(VDF,2), i in axes(VDF,1)
      vx = vmin[1] + (i - 0.5f0)*dv[1]
      vy = vmin[2] + (j - 0.5f0)*dv[2]
      vz = vmin[3] + (k - 0.5f0)*dv[3]
      dv2 = (vx - u[1])^2 + (vy - u[2])^2 + (vz - u[3])^2

      g = n * sqrt(vth2Inv/π) * (vth2Inv / π) * exp(-vth2Inv*dv2)
      ϵₘ += abs(VDF[i,j,k] - g)
   end

   ϵₘ = -log(0.5 / n * convert(eltype(VDF), prod(dv)) * ϵₘ)
end

function getmaxwellianity(meta, vcellids, vcellf; species="proton")
   (;vblock_size, vblocks, dv, vmin) = meta.meshes[species]
   vsize = @inbounds ntuple(i -> vblock_size[i] * vblocks[i], Val(3))
   slicez = vsize[1]*vsize[2]
   blocksize = prod(vblock_size)
   sliceBz = vblocks[1]*vblocks[2]
   sliceCz = vblock_size[1]*vblock_size[2]

   n = getdensity(meta, vcellf)
   u = getvelocity(meta, vcellids, vcellf)
   P = getpressure(meta, vcellids, vcellf)
   p = (P[1] + P[2] + P[3]) / 3
   T = p / (n *kB) # temperature from scalar pressure

   ϵₘ = zero(eltype(vcellf))
   vth2Inv = mᵢ / (2kB*T)

   @inbounds @simd for ic in eachindex(vcellids)
      id = findindex(vcellids[ic], vblocks, vblock_size, blocksize, vsize, sliceBz, sliceCz)
      i = id % vsize[1]
      j = id % slicez ÷ vsize[1]
      k = id ÷ slicez

      vx = vmin[1] + (i + 0.5f0)*dv[1]
      vy = vmin[2] + (j + 0.5f0)*dv[2]
      vz = vmin[3] + (k + 0.5f0)*dv[3]
      dv2 = (vx - u[1])^2 + (vy - u[2])^2 + (vz - u[3])^2

      g = n * sqrt(vth2Inv/π) * (vth2Inv / π) * exp(-vth2Inv*dv2)
      ϵₘ += abs(vcellf[ic] - g)
   end

   ϵₘ = -log(0.5 / n * convert(eltype(vcellf), prod(dv)) * ϵₘ)
end

function isInsideDomain(meta::MetaVLSV, point)
   (;coordmin, coordmax) = meta

   if coordmin[1] < point[1] ≤ coordmax[1] &&
      coordmin[2] < point[2] ≤ coordmax[2] &&
      coordmin[3] < point[3] ≤ coordmax[3]
      return true
   else
      return false
   end
end

"""
    getcellinline(meta, point1, point2) -> cellids, distances, coords

Returns cell IDs, distances and coordinates for every cell in a line between two given
points `point1` and `point2`. TODO: preallocation?
"""
function getcellinline(meta::MetaVLSV, point1::Vector{T}, point2::Vector{T}) where T
   (;coordmin, coordmax, ncells) = meta

   if !isInsideDomain(meta, point1)
      throw(DomainError(point1, "point location outside simulation domain!"))
   elseif !isInsideDomain(meta, point2)
      throw(DomainError(point2, "point location outside simulation domain!"))
   end

   cell_lengths = @inbounds ntuple(i -> (coordmax[i] - coordmin[i]) / ncells[i], Val(3))

   distances = [zero(T)]
   cellids = [getcell(meta, point1)]
   coords = [SVector{3}(point1)]
   ϵ = eps(T)
   unit_vector = @. (point2 - point1) / $norm(point2 - point1 + ϵ)
   p = coords[1]
   coef_min = zeros(T, 3)
   coef_max = zeros(T, 3)

   @inbounds while true
      cid = getcell(meta, p)
      amrlvl = getlevel(meta, cid)

      # Get the max and min cell boundaries
      min_bounds = getcellcoordinates(meta, cid) .- 0.5.*cell_lengths.*0.5.^amrlvl
      max_bounds = min_bounds .+ cell_lengths

      # Check which face we hit first
      @. coef_min = (min_bounds - p) / unit_vector
      @. coef_max = (max_bounds - p) / unit_vector

      # Negative coefficients indicates the opposite direction
      for i = 1:3
         if unit_vector[i] == 0.0
            coef_min[i] = Inf
            coef_max[i] = Inf
         end
         if coef_min[i] ≤ 0  coef_min[i] = Inf end
         if coef_max[i] ≤ 0  coef_max[i] = Inf end
      end

      # Find the minimum distance from a boundary times a factor
      d = min(minimum(coef_min), minimum(coef_max)) * 1.00001

      coordnew = SVector(
         p[1] + d*unit_vector[1],
         p[2] + d*unit_vector[2],
         p[3] + d*unit_vector[3])

      dot(point2 .- coordnew, unit_vector) ≥ 0 || break

      cellidnew = getcell(meta, coordnew)

      push!(cellids, cellidnew)
      push!(coords, coordnew)
      push!(distances, norm(coordnew .- point1))

      p = coordnew
   end

   cellids, distances, coords
end

"""
    getslicecell(meta, sliceoffset, idim, minCoord, maxCoord) -> idlist, indexlist

Find the cell IDs `idlist` which are needed to plot a 2d cut through of a 3d mesh, in a
direction `idim` at `sliceoffset`, and the `indexlist`, which is a mapping from original
order to the cut plane and can be used to select data onto the plane.
"""
function getslicecell(meta::MetaVLSV, sliceoffset, idim, minCoord, maxCoord)
   idim ∉ (1,2,3) && @error "Unknown slice direction $idim"
   (;ncells, maxamr, cellid, cellindex) = meta

   nsize = ncells[idim]
   sliceratio = sliceoffset / (maxCoord - minCoord)
   0.0 ≤ sliceratio ≤ 1.0 || error("slice plane index out of bound!")

   # Find the ids
   nlen = 0
   ncell = prod(ncells)
   # number of cells up to each refinement level
   nStart = (vcat(0, accumulate(+, (ncell*8^ilvl for ilvl = 0:maxamr))))

   indexlist = Int[]
   idlist = UInt[]

   cellidsorted = cellid[cellindex]

   @inbounds for ilvl = 0:maxamr
      nLow, nHigh = nStart[ilvl+1], nStart[ilvl+2]
      ids = cellidsorted[nLow .< cellidsorted .≤ nHigh]
      ix, iy, iz = getindexes(ilvl, ncells[1], ncells[2], nLow, ids)

      coords =
         if idim == 1
            ix
         elseif idim == 2
            iy
         else # 3
            iz
         end

      # Find the cut plane index for each refinement level (0-based)
      depth = floor(Int, sliceratio*nsize*2^ilvl)
      # Find the needed elements to create the cut and save the results
      elements = coords .== depth
      append!(indexlist, (nlen+1:nlen+length(ids))[elements])
      append!(idlist, ids[elements])

      nlen += length(ids)
   end

   idlist, indexlist
end

"""
    refineslice(meta, idlist, data, normal) -> Array

Generate scalar data on the finest refinement level given cellids `idlist` and variable
`data` on the slice perpendicular to `normal`.
"""
function refineslice(meta::MetaVLSV, idlist, data, normal)
   (;ncells, maxamr) = meta

   dims = let ratio = 2^maxamr
      if normal == :x
         i1, i2 = 2, 3
      elseif normal == :y
         i1, i2 = 1, 3
      elseif normal == :z
         i1, i2 = 1, 2
      end
      (ncells[i1]*ratio, ncells[i2]*ratio)
   end

   dpoints = zeros(eltype(data), dims...)

   # Create the plot grid
   ncell = prod(ncells)
   nHigh, nLow = ncell, 0

   @inbounds for i = 0:maxamr
      ids = idlist[nLow .< idlist .≤ nHigh]
      d = data[nLow .< idlist .≤ nHigh]

      ix, iy, iz = getindexes(i, ncells[1], ncells[2], nLow, ids)

      # Get the correct coordinate values and the widths for the plot
      a, b =
         if normal == :x
            iy, iz
         elseif normal == :y
            ix, iz
         elseif normal == :z
            ix, iy
         end

      # Insert the data values into dpoints
      refineRatio = 2^(maxamr - i)
      iRange = 0:refineRatio-1
      X, Y = ndgrid(iRange, iRange)

      coords = [(0, 0) for _ in a, _ in 1:2^(2*(maxamr-i))]

      @inbounds for ir = 1:2^(2*(maxamr-i)), ic in eachindex(a, b)
         @fastmath coords[ic,ir] = (muladd(a[ic], refineRatio, 1+X[ir]),
                                    muladd(b[ic], refineRatio, 1+Y[ir]) )
      end

      for ir = 1:2^(2*(maxamr-i)), ic in eachindex(d)
         dpoints[ coords[ic,ir]... ] = d[ic]
      end

      nLow = nHigh
      nHigh += ncell*8^(i+1)
   end

   dpoints
end

"Compute every cell id's x, y and z indexes on the given refinement level (0-based)."
@inline function getindexes(ilevel, xcells, ycells, nCellUptoLowerLvl, ids)
   slicesize = xcells*ycells*4^ilevel

   iz = @. (ids - nCellUptoLowerLvl - 1) ÷ slicesize
   iy = similar(iz)
   ix = similar(iz)
   @inbounds for i in eachindex(ids, iz)
      # number of ids up to the coordinate z in the refinement level ilevel
      idUpToZ = muladd(iz[i], slicesize, nCellUptoLowerLvl)
      iy[i] = (ids[i] - idUpToZ - 1) ÷ (xcells*2^ilevel)
      ix[i] = ids[i] - idUpToZ - iy[i]*xcells*2^ilevel - 1
   end
   ix, iy, iz
end

@inline function getindexes(ilvl, xcells, ycells, nCellUptoLowerLvl, id::Integer)
   slicesize = xcells*ycells*4^ilvl
   iz = (id - nCellUptoLowerLvl - 1) ÷ slicesize
   idUpToZ = muladd(iz, slicesize, nCellUptoLowerLvl)
   iy = (id - idUpToZ - 1) ÷ (xcells*2^ilvl)
   ix = id - idUpToZ - iy*xcells*2^ilvl - 1
   ix, iy, iz
end

"""
    getnearestcellwithvdf(meta, id) -> UInt

Find the nearest spatial cell with VDF saved of a given cell `id` associated with `meta`.
"""
function getnearestcellwithvdf(meta::MetaVLSV, id)
   cells = getcellwithvdf(meta)
   isempty(cells) && throw(ArgumentError("No distribution saved in $(meta.name)"))
   coords = [(0.0f0, 0.0f0, 0.0f0) for _ in cells]
   @inbounds for i in eachindex(cells)
      coords[i] = getcellcoordinates(meta, cells[i])
   end
   coords_orig = getcellcoordinates(meta, id)
   d2 = [sum((c .- coords_orig).^2) for c in coords]
   cells[argmin(d2)]
end

"""
    getcellwithvdf(meta) -> cellids

Get all the cell IDs with VDF saved associated with `meta`.
"""
function getcellwithvdf(meta::MetaVLSV)
   fid, footer = meta.fid, meta.footer
   cellsWithVDF = readmesh(fid, footer, "SpatialGrid", "CELLSWITHBLOCKS")::Vector{UInt}
   nblock_C = readmesh(fid, footer, "SpatialGrid", "BLOCKSPERCELL")::Vector{UInt32}

   innerBCCells = findall(==(0), nblock_C)

   deleteat!(cellsWithVDF, innerBCCells)
   cellsWithVDF
end

"Return the first cellid - 1 on `mylevel` given `ncells` on this level."
function get1stcell(mylevel, ncells)
   cid1st = 0
   for i = 0:mylevel-1
      cid1st += ncells*8^i
   end
   cid1st
end

fillmesh(meta::MetaVLSV, vars::AbstractString) = fillmesh(meta, [vars])

"""
    fillmesh(meta::MetaVLSV, vars; verbose=false) -> celldata, vtkGhostType

Fill the DCCRG mesh with quantity of `vars` on all refinement levels.
# Return arguments
- `celldata::Vector{Vector{Array}}`: data for each variable on each AMR level.
- `vtkGhostType::Array{UInt8}`: cell status (to be completed!).
"""
function fillmesh(meta::MetaVLSV, vars; verbose=false)
   (;maxamr, fid, footer, ncells, cellid, cellindex) = meta

   nvarvg = findall(!startswith("fg_"), vars)
   nv = length(vars)
   T = Vector{DataType}(undef, nv)
   offset = zeros(Int, nv)
   arraysize = zeros(Int, nv)
   vsize  = zeros(Int, nv)
   @inbounds for i = 1:nv
      T[i], offset[i], arraysize[i], _, vsize[i] =
         getObjInfo(footer, vars[i], "VARIABLE", "name")
   end

   Tout = copy(T)
   for i in eachindex(T)
      if T[i] == Float64 Tout[i] = Float32 end
   end

   @inbounds celldata =
      [[zeros(Tout[iv], vsize[iv], ncells[1] << i, ncells[2] << i, ncells[3] << i)
      for i = 0:maxamr] for iv in 1:nv]

   @inbounds vtkGhostType =
      [zeros(UInt8, ncells[1] << i, ncells[2] << i, ncells[3] << i) for i = 0:maxamr]

   if maxamr == 0
      @inbounds for iv = 1:nv
         celldata[iv][1][:] = readvariable(meta, vars[iv])
      end
      return celldata, vtkGhostType
   end

   # Find the ids
   ncell = prod(ncells)
   nLow, nHigh = 0, ncell
   cellidsorted = cellid[cellindex]

   @inbounds for ilvl = 0:maxamr
      verbose && @info "scanning AMR level $ilvl..."

      idfirst_ = searchsortedfirst(cellidsorted, nLow+1)
      idlast_  = searchsortedlast(cellidsorted, nHigh)

      ids = cellidsorted[idfirst_:idlast_]

      # indicate the condition of non-existing cells
      idrefined = setdiff(nLow+1:nHigh, ids)

      @simd for id in idrefined
         ix, iy, iz = getindexes(ilvl, ncells[1], ncells[2], nLow, id) .+ 1
         vtkGhostType[ilvl+1][ix,iy,iz] = 8
      end

      rOffsetsRaw = indexin(ids, cellid)

      if ilvl != maxamr
         for iv in nvarvg
            verbose && @info "reading variable $(vars[iv])..."
            a = mmap(fid, Vector{UInt8}, sizeof(T[iv])*vsize[iv]*arraysize[iv], offset[iv])
            dataRaw = reshape(reinterpret(T[iv], a), vsize[iv], arraysize[iv])
            data = @view dataRaw[:,rOffsetsRaw]

            fillcell!(ilvl, ids, ncells, maxamr, nLow, celldata[iv], data)
         end
      else # max refinement level
         for (iv, var) = enumerate(vars)
            verbose && @info "reading variable $var..."
            if startswith(var, "fg_")
               celldata[iv][end][:] = readvariable(meta, var)
            else
               a = mmap(fid, Vector{UInt8}, sizeof(T[iv])*vsize[iv]*arraysize[iv],
                  offset[iv])
               dataRaw = reshape(reinterpret(T[iv], a), vsize[iv], arraysize[iv])
               data = @view dataRaw[:,rOffsetsRaw]

               fillcell!(ids, ncells, maxamr, nLow, celldata[iv][end], data)
            end
         end
      end
      nLow = nHigh
      nHigh += ncell*8^(ilvl+1)
   end

   celldata, vtkGhostType
end

function fillcell!(ilvl, ids, ncells, maxamr, nLow, dataout, datain)
   @inbounds for ilvlup = ilvl:maxamr
      r = 2^(ilvlup-ilvl) # ratio on refined level
      for c in eachindex(ids)
         ixr, iyr, izr = getindexes(ilvl, ncells[1], ncells[2], nLow, ids[c]) .* r
         for k = 1:r, j = 1:r, i = 1:r
            _fillcelldata!(dataout[ilvlup+1], datain, ixr+i, iyr+j, izr+k, c)
         end
      end
   end
end

function fillcell!(ids, ncells, maxamr, nLow, dataout, datain)
   @inbounds for i in eachindex(ids)
      ix, iy, iz = getindexes(maxamr, ncells[1], ncells[2], nLow, ids[i]) .+ 1
      _fillcelldata!(dataout, datain, ix, iy, iz, i)
   end
end

@inline function _fillcelldata!(dataout, datain, i, j, k, index)
   @inbounds @simd for icomp in axes(datain,1)
      dataout[icomp,i,j,k] = datain[icomp,index]
   end
end

"""
    write_vtk(meta::MetaVLSV; kwargs...)
    write_vtk(file; kwargs...)

Convert VLSV file to VTK format.
# Keyword arguments
- `vars=[""]`: select which variables to convert.
- `ascii=false`: output stored in ASCII or compressed binary format.
- `maxamronly=false`: generate image files on the highest refinement level only.
- `verbose=false`: display logs during conversion.
"""
function write_vtk(meta::MetaVLSV; vars=[""], ascii=false, maxamronly=false, verbose=false)
   (;ncells, maxamr, dcoord, coordmin) = meta

   append = ascii ? false : true

   filedata = Vector{String}(undef, maxamr+1)
   @inbounds for i in 1:maxamr+1
      filedata[i] = meta.name[1:end-5]*"_$i.vti"
   end

   if isempty(vars[1])
      vars = meta.variable
      cellid_ = findfirst(==("CellID"), vars)
      if !isnothing(cellid_) deleteat!(vars, cellid_) end
   end

   data, vtkGhostType = fillmesh(meta, vars; verbose)

   if maxamronly
      save_image(meta, meta.name[1:end-4]*"vti", vars, data, vtkGhostType[end], maxamr,
         append)
   else
      # Generate image file on each refinement level
      @inbounds for i in eachindex(vtkGhostType, filedata)
         fdata, ghost = filedata[i], vtkGhostType[i]
         save_image(meta, fdata, vars, data, ghost, i-1, append)
      end

      # Generate vthb file
      filemeta = meta.name[1:end-4]*"vthb"
      doc = XMLDocument()
      elm = ElementNode("VTKFile")
      setroot!(doc, elm)
      link!(elm, AttributeNode("type", "vtkOverlappingAMR"))
      link!(elm, AttributeNode("version", "1.1"))
      link!(elm, AttributeNode("byte_order", "LittleEndian")) # x86
      link!(elm, AttributeNode("header_type", "UInt64"))

      xamr = addelement!(elm, "vtkOverlappingAMR")

      origin = @sprintf "%f %f %f" coordmin[1] coordmin[2] coordmin[3]
      link!(xamr, AttributeNode("origin", origin))
      link!(xamr, AttributeNode("grid_description", "XYZ"))

      @inbounds for i = 0:maxamr
         xBlock = addelement!(xamr, "Block")
         link!(xBlock, AttributeNode("level", string(i)))
         spacing_str = @sprintf "%f %f %f" dcoord[1]/2^i dcoord[2]/2^i dcoord[3]/2^i
         link!(xBlock, AttributeNode("spacing", spacing_str))
         xDataSet = addelement!(xBlock, "DataSet")
         link!(xDataSet, AttributeNode("index", "0"))
         amr_box = (0, ncells[1]*2^i-1, 0, ncells[2]*2^i-1, 0, ncells[3]*2^i-1)
         box_str = @sprintf("%d %d %d %d %d %d", amr_box[1], amr_box[2], amr_box[3],
            amr_box[4], amr_box[5], amr_box[6])
         link!(xDataSet, AttributeNode("amr_box", box_str))
         link!(xDataSet, AttributeNode("file", filedata[i+1]))
      end

      write(filemeta, doc)
   end

   return
end

write_vtk(file; kwargs...) = write_vtk(load(file); kwargs...)

"""
    save_image(meta::MetaVLSV, file, vars, data, vtkGhostType, level,
       ascii=false, append=true)

Save `data` of name `vars` at AMR `level` into VTK image file of name `file`.
# Arguments
- `file::String`: output file name.
- `vars::Vector{String}`: variable names to be saved.
- `data::Vector{Vector}`: data for all the variables on each refinement level.
- `vtkGhostType::Array{UInt8}`: array for visibility control.
- `level::Int`: refinement level (0-based).
- `ascii=false`: save output in ASCII or binary format.
- `append=true`: determines whether to append data at the end of file or do in-block writing.
"""
function save_image(meta::MetaVLSV, file, vars, data, vtkGhostType, level, ascii=false,
   append=true)
   (;coordmin, dcoord, ncells) = meta
   origin = (coordmin[1], coordmin[2], coordmin[3])
   ratio = 2^level
   spacing = (dcoord[1] / ratio, dcoord[2] / ratio, dcoord[3] / ratio)

   vtk = vtk_grid(file, ncells[1]*ratio+1, ncells[2]*ratio+1, ncells[3]*ratio+1;
      origin, spacing, append, ascii)

   @inbounds for (iv, var) in enumerate(vars)
      vtk[var, VTKCellData()] = data[iv][level+1]
   end

   vtk["vtkGhostType", VTKCellData()] = vtkGhostType

   vtk_save(vtk)
end

"""
    write_vlsv(filein, fileout, newvars::Vector{Tuple{Vector, String, VarInfo}};
       force=false)

Generate a new VLSV `fileout` based on `filein`, with `newvars` added.
`force=true` overwrites the existing `fileout`.
"""
function write_vlsv(filein::AbstractString, fileout::AbstractString,
   newvars::Vector{Tuple{VecOrMat, String, VarInfo}}; force=false)
   if isfile(fileout) && !force
      error("Output target $fileout exists!")
   end

   fid = open(filein)
   endian_offset = 8 # First 8 bytes indicate big-endian or else
   seek(fid, endian_offset)
   # Obtain the offset of the XML footer
   offset = read(fid, UInt64)
   # Store all non-footer part as raw data
   raw_data = zeros(UInt8, offset)

   seekstart(fid)
   readbytes!(fid, raw_data, offset)
   # Read input VLSV file footer
   doc = read(fid, String) |> parsexml
   footer = doc |> root
   close(fid)
   # Get new variables' offsets
   offsets = accumulate(+,
      [offset, [sizeof(newvars[i][1]) for i in eachindex(newvars)[1:end-1]]...])
   # Create new children for footer
   for i in eachindex(newvars, offsets)
      elm = addelement!(footer, "VARIABLE", string(offsets[i]))

      a1 = AttributeNode("arraysize", string(length(newvars[i][1])))
      a2 = AttributeNode("datasize", string(sizeof(eltype(newvars[i][1]))))
      a3 =
         if eltype(newvars[i][1]) <: Signed
            AttributeNode("datatype", "int")
         elseif eltype(newvars[i][1]) <: AbstractFloat
            AttributeNode("datatype", "float")
         elseif eltype(newvars[i][1]) <: Unsigned
            AttributeNode("datatype", "uint")
         end
      a4 = AttributeNode("mesh", "SpatialGrid")
      a5 = AttributeNode("name", newvars[i][2])
      a6 = AttributeNode("unit", newvars[i][3].unit)
      a7 = AttributeNode("unitConversion", newvars[i][3].unitConversion)
      a8 = AttributeNode("unitLaTeX", newvars[i][3].unitLaTeX)
      a9 = AttributeNode("variableLaTeX", newvars[i][3].variableLaTeX)
      a10 =
         if ndims(newvars[i][1]) == 1
            AttributeNode("vectorsize", "1")
         else
            AttributeNode("vectorsize", string(size(newvars[i][1], 1)))
         end

      for attributenode in (a1, a2, a3, a4, a5, a6, a7, a8, a9, a10)
         link!(elm, attributenode)
      end
   end
   # Write to fileout
   open(fileout, "w") do io
      write(io, @view raw_data[1:8]) # endianness
      # Compute footer offset
      totalnewsize = 0
      for var in newvars
         totalnewsize += sizeof(var[1])
      end
      write(io, offset+totalnewsize) # record new footer offset
      write(io, @view raw_data[17:end]) # copy original data
      for var in newvars
         write(io, var[1])
      end
      write(io, string(footer), '\n')
   end
   return
end

"""
    issame(file1, file2, tol=1e-4; verbose=false) -> Bool

Check if two VLSV files `file1` and `file2` are approximately identical, under relative
tolerance `tol`.
"""
function issame(f1, f2, tol::AbstractFloat=1e-4; verbose=false)
   # 1st sanity check: minimal filesize difference
   if abs(filesize(f1) - filesize(f2)) / filesize(f2) > 1e-2
      verbose && println("The sizes of files are already quite different!")
      return false
   end

   meta1 = load(f1)
   meta2 = load(f2)
   varnames = meta1.variable
   strskip = r"CellID|rank|blocks"
   deleteat!(varnames, findall(endswith(strskip), varnames))

   isIdentical = true
   for vname in varnames
      v1 = readvariable(meta1, vname)
      v2 = readvariable(meta2, vname)

      s1, s2 = sum(v1), sum(v2)
      if abs(s1 - s2) > tol * abs(s1) && abs(s1 - s2) > tol * abs(s2)
         isIdentical = false
         verbose && println("$vname is quite different!")
         break
      end
   end
   verbose && isIdentical && println("$f1 and $f2 are identical under tolerance $tol.")
   close(meta1.fid)
   close(meta2.fid)
   return isIdentical
end


text:
using DataFrames, FreqTables, MixedModels, Random, StatsModels
#using Test

const LMM = LinearMixedModel;
const GLMM = GeneralizedLinearMixedModel;
const HC = HelmertCoding();
const contrasts = Dict(:A=>HC,:B=>HC,:C=>HC);

const rng = Random.MersenneTwister(4321234);

const LETTERS = ['A':'Z';];
const letters = ['a':'z';];
const Letters = vcat(LETTERS, letters);
const plusminus = ['+','-'];
const dat = (
    A = repeat(plusminus, outer=4096),
    B = repeat(plusminus, inner=2, outer=2048),
    C = repeat(plusminus, inner=4, outer=1024),
    G = repeat(Letters[1:32], inner=8, outer=32),
    H = repeat(letters[1:16], inner=16, outer=32),
    I = repeat(Letters[1:32], inner=256),
    Y = ones(8192),
);
show(DataFrame(dat))
println()
println(freqtable(dat, :H, :G))  # G is nested within H
println(freqtable(dat, :I, :G))  # G and I are crossed
println(freqtable(dat, :A, :G))
println(freqtable(dat, :B, :G))
println(freqtable(dat, :C, :G))
println(freqtable(dat, :A, :I))
println(freqtable(dat, :B, :I))
println(freqtable(dat, :C, :I))
m1 = LMM(@formula(Y ~ 1 + A*B*C + (1|G) + (1|I)), dat, contrasts=contrasts);
β = rand(rng, 8);
show(refit!(simulate!(rng, m1, β=β, σ=0.4, θ=[0.8, 0.6])))
println()
@show(β);
println()
show(m1.optsum);

m2 = LMM(@formula(Y ~ 1 + A*B*C + (1+A+B+C|G) + (1|I)), dat, contrasts=contrasts);
show(refit!(m2, copy(response(m1))))
m3 = LMM(@formula(Y ~ 1 + A*B*C + (1+A+B+C|G) + (1+A+B+C|I)), dat, contrasts=contrasts);
show(refit!(m3, copy(response(m1))))
m4 = LMM(@formula(Y ~ 1 + A*B*C + (1+A*B*C|G) + (1+A*B*C|I)), dat, contrasts=contrasts);
show(refit!(m4, copy(response(m1))))

logistic(η) = inv(1 + exp(-η))
"""
    bresp(rng::AbstractRNG, v::Vector{<:AbstractFloat})

Return a random sample from `[Bernoulli(logistic(x)) for x in v]`
"""
bresp(rng::AbstractRNG, v::AbstractVector{<:AbstractFloat}) =
    rand(rng, length(v)) .< logistic.(v)

copyto!(dat.Y, bresp(rng, response(m1)))
g1 = GLMM(@formula(Y ~ 1 + A*B*C + (1|G) + (1|H)), dat, Bernoulli(), contrasts=contrasts)
fit!(g1, fast=true);
show(deviance(g1))
#fit!(g1)

ext: jl


text:
using KitBase, Plots
using KitBase.JLD2
using KitBase.ProgressMeter: @showprogress
pyplot()
cd(@__DIR__)

function flux_sgks!(
    fw::X,
    wL::Y,
    wR::Y,
    inK::Real,
    γ::Real,
    μᵣ::Real,
    ω::Real,
    dt::Real,
    dxL::Real,
    dxR::Real,
    dy::Real,
    swL = zero(wL)::Y,
    swR = zero(wR)::Y,
) where {X,Y}

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    Mu1, Mv1, Mw1, MuL1, MuR1 = gauss_moments(primL, inK)
    Mu2, Mv2, Mw2, MuL2, MuR2 = gauss_moments(primR, inK)

    w =
        primL[1] .* moments_conserve(MuL1, Mv1, Mw1, 0, 0, 0) .+
        primR[1] .* moments_conserve(MuR2, Mv2, Mw2, 0, 0, 0)
    prim = conserve_prim(w, γ)
    tau =
        vhs_collision_time(prim, μᵣ, ω) +
        2 * dt * abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end])

    Mu, Mv, Mw, MuL, MuR = gauss_moments(prim, inK)
    sw0L = (w .- (wL .- swL .* dxL)) ./ dxL
    sw0R = ((wR .+ swR .* dxR) .- w) ./ dxR
    gaL = pdf_slope(prim, sw0L, inK)
    gaR = pdf_slope(prim, sw0R, inK)
    sw =
        -prim[1] .* (
            moments_conserve_slope(gaL, MuL, Mv, Mw, 1, 0, 0) .+
            moments_conserve_slope(gaR, MuR, Mv, Mw, 1, 0, 0)
        )
    gaT = pdf_slope(prim, sw, inK)

    # time-integration constants
    Mt = zeros(5)
    Mt[4] = tau * (1.0 - exp(-dt / tau))
    Mt[5] = -tau * dt * exp(-dt / tau) + tau * Mt[4]
    Mt[1] = dt - Mt[4]
    Mt[2] = -tau * Mt[1] + Mt[5]
    Mt[3] = 0.5 * dt^2 - tau * Mt[1]

    # flux related to central distribution
    Muv = moments_conserve(Mu, Mv, Mw, 1, 0, 0)
    MauL = moments_conserve_slope(gaL, MuL, Mv, Mw, 2, 0, 0)
    MauR = moments_conserve_slope(gaR, MuR, Mv, Mw, 2, 0, 0)
    MauT = moments_conserve_slope(gaT, Mu, Mv, Mw, 1, 0, 0)

    fw .= Mt[1] .* prim[1] .* Muv

    # flux related to upwind distribution
    MuvL = moments_conserve(MuL1, Mv1, Mw1, 1, 0, 0)
    #MauL = moments_conserve_slope(faL, MuL1, Mv1, Mw1, 2, 0, 0)
    #MauLT = moments_conserve_slope(faTL, MuL1, Mv1, Mw1, 1, 0, 0)

    MuvR = moments_conserve(MuR2, Mv2, Mw2, 1, 0, 0)
    #MauR = moments_conserve_slope(faR, MuR2, Mv2, Mw2, 2, 0, 0)
    #MauRT = moments_conserve_slope(faTR, MuR2, Mv2, Mw2, 1, 0, 0)

    #@. fw +=
    #    Mt[4] * primL[1] * MuvL - (Mt[5] + tau * Mt[4]) * primL[1] * MauL -
    #    tau * Mt[4] * primL[1] * MauLT + Mt[4] * primR[1] * MuvR -
    #    (Mt[5] + tau * Mt[4]) * primR[1] * MauR - tau * Mt[4] * primR[1] * MauRT
    @. fw += Mt[4] * primL[1] * MuvL + Mt[4] * primR[1] * MuvR

    fw .*= dy

    return nothing

end

function ev!(KS, ctr, a1face, a2face, dt)
    nx, ny, dx, dy = KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
    idx0, idx1 = 1, nx+1
    idy0, idy1 = 1, ny+1

    # r direction
    @inbounds Threads.@threads for j = 1:ny
        w = local_frame(ctr[1, j].w, a1face[1, j].n[1], a1face[1, j].n[2])
        flux_boundary_maxwell!(
            a1face[1, j].fw,
            [1.0, 0.0, 0.0, 1.0],
            w,
            KS.gas.K,
            KS.gas.γ,
            dt,
            a1face[1, j].len,
            1,
        )
        a1face[1, j].fw .=
                global_frame(a1face[1, j].fw, a1face[1, j].n[1], a1face[1, j].n[2])
    end
    @inbounds Threads.@threads for j = 2:ny
        for i = idx0:idx1
            wL = local_frame(ctr[i-1, j].w, a1face[i, j].n[1], a1face[i, j].n[2])
            wR = local_frame(ctr[i, j].w, a1face[i, j].n[1], a1face[i, j].n[2])

            flux_sgks!(
                a1face[i, j].fw,
                wL,
                wR,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                dt,
                ks.ps.dr[i-1, j] / 2,
                ks.ps.dr[i, j] / 2,
                a1face[i, j].len,
            )

            a1face[i, j].fw .=
                global_frame(a1face[i, j].fw, a1face[i, j].n[1], a1face[i, j].n[2])
        end
    end

    # θ direction
    @inbounds Threads.@threads for j = idy0:idy1
        for i = 1:nx
            wL = local_frame(ctr[i, j-1].w, a2face[i, j].n[1], a2face[i, j].n[2])
            wR = local_frame(ctr[i, j].w, a2face[i, j].n[1], a2face[i, j].n[2])

            flux_sgks!(
                a2face[i, j].fw,
                wL,
                wR,
                KS.gas.K,
                KS.gas.γ,
                KS.gas.μᵣ,
                KS.gas.ω,
                dt,
                ks.ps.darc[i, j-1] / 2,
                ks.ps.darc[i, j] / 2,
                a2face[i, j].len,
            )

            a2face[i, j].fw .=
                global_frame(a2face[i, j].fw, a2face[i, j].n[1], a2face[i, j].n[2])
        end
    end

    return nothing
end

begin
    set = Setup(
        case = "cylinder",
        space = "2d0f0v",
        boundary = ["isothermal", "extra", "mirror", "mirror"],
        limiter = "minmod",
        cfl = 0.1,
        maxTime = 15.0, # time
    )
    ps = CSpace2D(1.0, 6.0, 60, 0.0, π, 50, 1, 1)
    vs = nothing
    gas = Gas(Kn = 1e-2, Ma = 5.0, K = 1.0)

    prim0 = [1.0, 0.0, 0.0, 1.0]
    prim1 = [1.0, gas.Ma * sound_speed(1.0, gas.γ), 0.0, 1.0]
    fw = (args...) -> prim_conserve(prim1, gas.γ)
    bc = function(x, y)
        if abs(x^2 + y^2 - 1) < 1e-3
            return prim0
        else
            return prim1
        end
    end
    ib = IB(fw, bc)

    ks = SolverSet(set, ps, vs, gas, ib)
    ctr, a1face, a2face = init_fvm(ks, ks.ps)
end

t = 0.0
dt = timestep(ks, ctr, 0.0)
nt = ks.set.maxTime ÷ dt |> Int
res = zeros(4)

@showprogress for iter = 1:nt
    ev!(ks, ctr, a1face, a2face, dt)
    update!(ks, ctr, a1face, a2face, dt, res)

    for j = ks.ps.nθ÷2+1:ks.ps.nθ
        ctr[ks.ps.nr+1, j].w .= ks.ib.fw(6, 0)
        ctr[ks.ps.nr+1, j].prim .= conserve_prim(ctr[ks.ps.nr+1, j].w, ks.gas.γ)
    end
    for j = 0:ks.ps.nθ
        KitBase.bc_isothermal!(ctr[0, j], ctr[1, j], ks.gas.γ)
    end

    global t += dt
    if iter % 1000 == 0
        println("residuals: $(res)")
    end

    if maximum(res) < 1e-6
        break
    end
end
@save "kn2ns.jld2" ctr a1face a2face

begin
    sol = zeros(ks.ps.nr, ks.ps.nθ, 4)
    for i in axes(sol, 1), j in axes(sol, 2)
        sol[i, j, :] .= ctr[i, j].prim
        sol[i, j, end] = 1 / sol[i, j, end]
    end

    contourf(
        ps.x[1:ks.ps.nr, 1:ks.ps.nθ],
        ps.y[1:ks.ps.nr, 1:ks.ps.nθ],
        sol[:, :, 4],
        ratio = 1,
        xlabel = "x",
        ylabel = "y",
    )
end

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
struct GomokuEnv <: GameEnv
  N::Int
  n_in_row::Int
  action_space::Int
  planes::Int
  max_action_space::Int
  EMPTY_BOARD::Array{Int8, 2}

  function GomokuEnv(board_size::Int = 15, connect_row = 5, planes::Int = 17)
    N = board_size
    n_in_row = connect_row
    action_space = N ^ 2
    EMPTY_BOARD = zeros(Int8, N, N)
    @assert planes % 2 == 1
    planes = (planes - 1) ÷ 2

    new(N, n_in_row, action_space, planes, 361 , EMPTY_BOARD)
  end
end

function Base.show(io::IO, env::GomokuEnv)
  N, to_win = env.N, env.n_in_row
  println("Gomok(board size: $N, connect to win: $to_win")
end

include("board.jl")
include("coords.jl")


text:
import DataStructures: MutableBinaryHeap, MutableBinaryMinHeap
using Random

"""
Next reaction by Hazards
Also called Anderson's method.
"""
struct TransitionRecord
	exponential_interval::Float64
	heap_handle::Int64
end

mutable struct NextReactionHazards
	firing_queue::MutableBinaryHeap{NRTransition,Base.Order.ForwardOrdering}
	transition_state::Dict{Any,TransitionRecord}
	init::Bool
end

"""
Construct a Next Reaction sampler.
"""
function NextReactionHazards()
    heap = MutableBinaryMinHeap{NRTransition}()
    @debug("SampleSemiMarkov.NextReactionHazards type ", typeof(heap))
    state = Dict{Any,TransitionRecord}()
    NextReactionHazards(heap, state, true)
end

# Finds the next one without removing it from the queue.
function Next(propagator::NextReactionHazards, system, rng)
	if propagator.init
		Hazards(system, rng) do clock, now, updated, rng2
			Enable(propagator, clock, now, updated, rng2)
	    end
	    propagator.init = false
	end

	NotFound = NRTransition(nothing, Inf)
	if !isempty(propagator.firing_queue)
		least = top(propagator.firing_queue)
	else
		least = NotFound
	end
	@debug("SampleSemiMarkov.next queue length ",
			length(propagator.firing_queue), " least ", least)
	(least.time, least.key)
end

"""
Returns an observer of intensities to decide what to
do when they change.
"""
function Observer(propagator::NextReactionHazards)
	function nrobserve(clock, time, updated, rng)
		if updated == :Disabled || updated == :Fired
			Disable(propagator, clock, time, updated, rng)
		else
			Enable(propagator, clock, time, updated, rng)
		end
	end
end

function unit_hazard_interval(rng::MersenneTwister)
	-log(rand(rng))
end

# Enable or modify a hazard.
function Enable(propagator::NextReactionHazards, clock,
		now, updated, rng)
	key = clock
	clock_started = haskey(propagator.transition_state, key)
	if clock_started
		record = propagator.transition_state[key]
		when_fire = Putative(clock.intensity, now, record.exponential_interval)

		@assert(when_fire >= now)
		if record.heap_handle >= 0
			@debug("SampleSemiMarkov.enable keyu ", key, " interval ",
				record.exponential_interval, " when ", when_fire,
				" dist ", clock)
			update!(propagator.firing_queue, record.heap_handle,
				NRTransition(key, when_fire))
		else
			record.heap_handle = push!(propagator.firing_queue,
				NRTransition(key, when_fire))
			@debug("SampleSemiMarkov.enable keyp ", key, " interval ",
				record.exponential_interval, " when ", when_fire,
				" dist ", clock)
		end
	else
		firing_time, interval = MeasuredSample(clock.intensity, now, rng)
		@assert(firing_time >= now)
        handle = push!(propagator.firing_queue, NRTransition(key, firing_time))
        @debug("SampleSemiMarkov.enable Adding key ", key, " interval ",
        	interval, " when ", firing_time, " dist ", clock)
		record = TransitionRecord(interval, handle)
		propagator.transition_state[key] = record
	end
    @debug("SampleSemiMarkov.enable exit")
end

# Remove a transition from the queue because it was disabled.
function Disable(propagator::NextReactionHazards, key, now,
        updated, rng)
	record = propagator.transition_state[key]
	# We store distributions in order to calculate remaining hazard
	# which will happen AFTER the state has changed.
	update!(propagator.firing_queue, record.heap_handle,
		NRTransition(key, -1.))
	todelete = pop!(propagator.firing_queue)
	@assert(todelete.key == key && todelete.time == -1)
    if updated == :Disabled
    	record.heap_handle = -1 # This is the official sign it was disabled.
    elseif updated == :Fired
        # Deleting the key is slower for small, finite systems,
        # but it makes infinite (meaning long-running) systems possible.
        delete!(propagator.transition_state, key)
    else
        assert(updated == :Disabled || updated == :Fired)
    end
end

function print_next_reaction_hazards(propagator::NextReactionHazards)
    arr = Any[]
    for x in keys(propagator.transition_state)
        push!(arr, x)
    end
    sort!(arr)
    for trans in arr
        rec = propagator.transition_state[trans]
        if rec.distribution !== nothing
	        p=parameters(rec.distribution)
    	end
    end
end


text:
module PooksoftAlphaVantageDataStore

# include -
include("Include.jl")

# export functions -
# low-level functions -
export log_api_call
export execute_sts_intraday_api_call
export execute_sts_daily_api_call
export execute_sts_adjusted_daily_api_call
export execute_sts_weekly_api_call
export execute_sts_adjusted_weekly_api_call
export execute_sts_monthly_api_call
export execute_sts_adjusted_monthly_api_call
export execute_sts_quote_api_call
export execute_sts_search_api_call

export execute_simple_moving_average_api_call
export execute_exponential_moving_average_api_call
export execute_relative_strength_index_api_call
export execute_company_overview_api_call
export execute_company_earnings_api_call
export execute_company_income_statement_api_call

# high-level functions -
export build_api_user_model
export build_datastore_apicall_model
export execute_api_call

# read/write -
export write_data_to_filestore
export read_data_from_filestore

# export types -
export PSUserModel
export PSDataStoreAPICallModel

end # module


text:
@testset "Test NLP Stopping constrained" begin
  include("rosenbrock.jl")
  #Warning: see https://github.com/JuliaSmoothOptimizers/NLPModels.jl/blob/master/src/autodiff_model.jl
  #for the proper way of defining an ADNLPModel
  x0 = ones(6)
  c(x) = [sum(x)]
  meta = NLPModelMeta(
    6,
    x0 = x0,
    lvar = fill(-10.0, size(x0)),
    uvar = fill(10.0, size(x0)),
    ncon = 1,
    y0 = [0.0],
    lcon = [-Inf],
    ucon = [6.0],
  )
  nlp2 =
    ADNLPModel(meta, Counters(), ADNLPModels.ForwardDiffAD(6, 1, rosenbrock, x0), rosenbrock, c)

  nlp_at_x_c = NLPAtX(x0, NaN * ones(nlp2.meta.ncon))
  stop_nlp_c = NLPStopping(nlp2, nlp_at_x_c)

  a = zeros(6)
  fill_in!(stop_nlp_c, a)

  @test cons(nlp2, a) == stop_nlp_c.current_state.cx
  @test jac(nlp2, a) == stop_nlp_c.current_state.Jx

  @test stop!(stop_nlp_c) == false
  # we make sure the counter of stop works properly
  @test stop_nlp_c.meta.nb_of_stop == 1

  sol = ones(6)
  fill_in!(stop_nlp_c, sol)

  @test stop!(stop_nlp_c) == true
  @test stop_nlp_c.current_state.current_score == KKT(stop_nlp_c.pb, stop_nlp_c.current_state)

  stop_nlp_default = NLPStopping(nlp2, atol = 1.0)
  fill_in!(stop_nlp_default, sol)
  @test stop_nlp_default.meta.atol == 1.0
  @test stop!(stop_nlp_default) == true

  stop_nlp_default.meta.optimality_check = unconstrained_check
  update_and_stop!(stop_nlp_default, cx = nothing)
  @test stop_nlp_default.current_state.cx != nothing

  @test stop_nlp_default.pb.meta.minimize
  reinit!(stop_nlp_default, rstate = true)
  update!(stop_nlp_default, x = zeros(6), cx = [Inf])
  stop!(stop_nlp_default)
  @test status(stop_nlp_default, list=true) == [:Unknown]

  reinit!(stop_nlp_default, rstate = true)
  update!(stop_nlp_default, fx = Inf)
  stop!(stop_nlp_default)
  @test status(stop_nlp_default, list=true) == [:Unknown]

  reinit!(stop_nlp_default, rstate = true)
  update!(stop_nlp_default, fx = -Inf)
  stop!(stop_nlp_default)
  @test status(stop_nlp_default, list=true) == [:UnboundedPb]

  #Keywords in the stop! call
  nlp_at_x_kargs = NLPAtX(x0, NaN * ones(nlp2.meta.ncon))
  stop_nlp_kargs = NLPStopping(
    nlp2,
    nlp_at_x_c,
    optimality_check = (x, y; test = 1.0, kwargs...) -> Stopping.KKT(x, y, kwargs...) + test,
  )
  fill_in!(stop_nlp_kargs, sol)
  @test stop!(stop_nlp_kargs) == false
  @test stop!(stop_nlp_kargs, test = 0.0) == true
end


text:
mutable struct WordClusters
    vocab::Vector{String}
    clusters::Vector{Int}
    vocab_hash::Dict{String, Int}
    function WordClusters(vocab, clusters)
        vocab_hash = Dict{String, Int}()
        for (i, word) in enumerate(vocab)
            vocab_hash[word] = i
        end
        new(vocab, clusters, vocab_hash)
    end
end

function show(io::IO, wc::WordClusters)
    print(io, "WordClusters $(length(wc.vocab)) words, $(length(unique(wc.clusters))) clusters")
end

"""
    vocabulary(wc)

Return all the vocabulary of the WordClusters `wc`.
"""
vocabulary(wc::WordClusters) = wc.vocab

"""
    in_vocabulary(wc, word)

For the WordCluters `wc`, return `true` if `word` is part of the
vocabulary of `wc` and `false` otherwise.
"""
in_vocabulary(wc::WordClusters, word::AbstractString) = word in wc.vocab

"""
    index(wc, word)

Return the index of `word` from the WordCluaters `wc`.
"""
index(wc::WordClusters, word) = wc.vocab_hash[word]

"""
    get_cluster(wc, word)

Return the cluster number for a word in the vocabulary.
"""
get_cluster(wc::WordClusters, word) = (idx = index(wc, word); wc.clusters[idx])

"""
    clusters(wc)

Return all the clusters from the WordClusters `wc`.
"""
clusters(wc::WordClusters) = sort(unique(wc.clusters))

"""
    get_words(wc, cluster)

For the WordCluster `wc`, return all the words from a given cluster
number `cluster`.
"""
function get_words(wc::WordClusters, cluster::Int)
    inds = findall(in(cluster), wc.clusters)
    return wc.vocab[inds]
end

"""
    wordclusters(fname)

Generate a WordClusters type object from the text file `fname`.
"""
function wordclusters(fname::AbstractString)
    vocab = String[]
    clusters = Int[]
    open(fname) do f
        entries = split(strip(readline(f)), ' ')
        while length(entries) == 2
            push!(vocab, entries[1])
            push!(clusters, parse(Int, entries[2]))
            entries = split(strip(readline(f)), ' ')
        end
    end
    return WordClusters(vocab, clusters)
end


text:
module DecisionMakingProblems

using Distributions
using Parameters
using Random
using LinearAlgebra
using GridInterpolations
using Parameters
using Statistics
using Printf

export
    MDP, HexWorld, StraightLineHexWorld, TwentyFortyEight, CartPole, MountainCar, LQR, CollisionAvoidance,
    POMDP, DiscretePOMDP, CryingBaby, MachineReplacement, Catch,
    SimpleGame, PrisonersDilemma, RockPaperScissors, Travelers,
    MG, PredatorPreyHexWorld, CirclePredatorPreyHexWorld,
    POMG, MultiCaregiverCryingBaby,
    DecPOMDP, CollaborativePredatorPreyHexWorld, SimpleCollaborativePredatorPreyHexWorld, CircleCollaborativePredatorPreyHexWorld

import Base: <, ==, rand, vec

include("support_code.jl")

include("mdp/mdp.jl")
include("mdp/discrete_mdp.jl")
include("mdp/2048.jl")
include("mdp/hexworld.jl")
include("mdp/simple_lqr.jl")
include("mdp/cart_pole.jl")
include("mdp/mountain_car.jl")
include("mdp/collision_avoidance.jl")

include("pomdp/pomdp.jl")
include("pomdp/discrete_pomdp.jl")
include("pomdp/crying_baby.jl")
include("pomdp/machine_replacement.jl")
include("pomdp/catch.jl")

include("simple_game/simplegame.jl")
include("simple_game/prisoners_dilemma.jl")
include("simple_game/travelers.jl")
include("simple_game/rock_paper_scissors.jl")

include("markov_game/mg.jl")
include("markov_game/predator_prey.jl")

include("pomg/pomg.jl")
include("pomg/multicaregiver.jl")

include("decpomdp/decpomdp.jl")
include("decpomdp/collab_predator_prey.jl")

end # module


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BisectFractionPartitioner(normal, fraction=0.5, maxiter=10)

A method for partitioning spatial data into two half spaces
defined by a `normal` direction and a `fraction` of points.
The partition is returned within `maxiter` bisection iterations.
"""
struct BisectFractionPartitioner{T,N} <: AbstractPartitioner
  normal::SVector{N,T}
  fraction::Float64
  maxiter::Int

  function BisectFractionPartitioner{T,N}(normal, fraction, maxiter) where {N,T}
    new(normalize(normal), fraction, maxiter)
  end
end

BisectFractionPartitioner(normal::SVector{N,T}, fraction=0.5, maxiter=10) where {T,N} =
  BisectFractionPartitioner{T,N}(normal, fraction, maxiter)

BisectFractionPartitioner(normal::NTuple{N,T}, fraction=0.5, maxiter=10) where {T,N} =
  BisectFractionPartitioner(SVector(normal), fraction, maxiter)

function partition(object, partitioner::BisectFractionPartitioner)
  bbox = boundbox(object)
  n = partitioner.normal
  f = partitioner.fraction
  c = center(bbox)
  d = diagonal(bbox)

  # maximum number of bisections
  maxiter = partitioner.maxiter

  iter = 0; p = 0
  a = c - d/2 * n
  b = c + d/2 * n
  while iter < maxiter
    m = (a + b) / 2

    p = partition(object, BisectPointPartitioner(n, m))
    g = nelms(p[1]) / nelms(object)

    g ≈ f && break
    g > f && (b = m)
    g < f && (a = m)

    iter += 1
  end

  p
end


text:
using PyPlot

# 1000 samples with ND gaussians
Ndims = [1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30]'
ts = [0.279603	0.637024	1.028485	1.003426	1.250441	1.374436	1.414937	1.588211	1.967722	2.024116	2.55868	2.589394	2.618068	2.919893	2.95502	3.258155	3.547371	3.728454	3.939705	4.212554	3.871222	4.072321	4.340777	4.284822	4.755864	4.718708	5.135071	5.260214	6.393843	6.390321]'
js = [2	3	5	5	6	7	7	8	9	9	10	11	11	12	12	14	14	14	14	15	14	15	13	14	15	15	15	15	17	17]'

slowdown = ts ./ts[1]

fig = figure(figsize= [10,10])
plt.plot(Ndims, js, linewidth = 2, color = "blue")

title("1000 samples of N-D gaussians", fontsize = 22)
xlabel("Num of dimension", fontsize = 22)
ylabel("Number of iterations", fontsize = 22)

fig = figure(figsize= [10,10])
plt.plot(Ndims, slowdown, linewidth = 2, color = "blue")

title("1000 samples of N-D gaussians", fontsize = 22)
xlabel("Num of dimension", fontsize = 22)
ylabel("slowdown", fontsize = 22)


text:
using DataBlockStorage
using Test

@testset "DataBlockStorage.jl" begin
    # Write your tests here.
end


text:
struct Node
    children::Dict{Char,Node}
    posi::Array{Int}
end

function parse_string(sentence::String)::Array{Tuple{String,Int}}
	ret::Array{Tuple{String,Int}, 1} = []

	done = false

	i::Int = 1
	while !done
		if sentence[i] != ' '
			stop = i + 1
			while stop <= length(sentence) && sentence[stop] != ' '
				stop += 1
			end
			push!(ret, (sentence[i:stop - 1], i))
			i = stop + 1
		else
			i += 1
		end

		if i > length(sentence)
			done = true
		end
	end

	return ret
end

function build(list_of_words::Array{Tuple{String,Int}})::Node
	top::Node = Node(Dict([]), [])

	for i in list_of_words
		cur_node = top

		for j in i[1]
			if !haskey(cur_node.children, j)
				push!(cur_node.children, j => Node(Dict([]), []))
			end

			cur_node = cur_node.children[j]
		end

		push!(cur_node.posi, i[2])
	end

	return top
end

function positions(word::String, node::Node, index::Int=1)::Array{Int}
	indices::Array{Int, 1} = []

	if index > length(word)
		return indices
	elseif index == length(word)
		if word[index] != '?'
			if haskey(node.children, word[index])
				append!(indices, node.children[word[index]].posi)
			end
		else
			for i in node.children
				append!(indices, i[2].posi)
			end
		end
		return indices
	end

	if word[index] == '?'
		for i in node.children
			append!(indices, positions(word, i[2], index + 1))
		end
	elseif haskey(node.children, word[index])
		append!(indices, positions(word, node.children[word[index]], index + 1))
	end

	return indices
end


text:
MeshCat.LinearMap(q::Quaternion) = MeshCat.LinearMap(rotation_matrix(q))
MeshCat.js_quaternion(q::Quaternion) = [q.v1, q.v2, q.v3, q.s]

RotX(θ) = Quaternion(cos(θ/2), sin(θ/2), 0, 0)
RotY(θ) = Quaternion(cos(θ/2), 0, sin(θ/2), 0)
RotZ(θ) = Quaternion(cos(θ/2), 0, 0, sin(θ/2))

quateltype(x) = eltype(x) # TODO not super elegant
quateltype(::Quaternion{T}) where T = T

vector(q::Quaternion) = SA[q.s, q.v1, q.v2, q.v3]
vector(q::AbstractVector) = q

function Lmat(q::Quaternion)
    SA[
        q.s  -q.v1 -q.v2 -q.v3;
        q.v1  q.s  -q.v3  q.v2;
        q.v2  q.v3  q.s  -q.v1;
        q.v3 -q.v2  q.v1  q.s;
    ]
end

function Rmat(q::Quaternion)
    SA[
        q.s  -q.v1 -q.v2 -q.v3;
        q.v1  q.s   q.v3 -q.v2;
        q.v2 -q.v3  q.s   q.v1;
        q.v3  q.v2 -q.v1  q.s;
    ]
end

Lᵀmat(q) = Lmat(q)'
Rᵀmat(q) = Rmat(q)'

function Tmat(::Type{T}=Float64) where T
    SA{T}[
        1  0  0  0;
        0 -1  0  0;
        0  0 -1  0;
        0  0  0 -1;
    ]
end

function Vmat(::Type{T}=Float64) where T
    SA{T}[
        0 1 0 0;
        0 0 1 0;
        0 0 0 1;
    ]
end

function Vᵀmat(::Type{T}=Float64) where T
    SA{T}[
        0 0 0;
        1 0 0;
        0 1 0;
        0 0 1;
    ]
end

Vmat(q::Quaternion) = SA[q.v1, q.v2, q.v3]

function VLmat(q::Quaternion)
    SA[
        q.v1  q.s -q.v3  q.v2;
        q.v2  q.v3  q.s -q.v1;
        q.v3 -q.v2  q.v1  q.s;
    ]
end

function VLᵀmat(q::Quaternion)
    SA[
        -q.v1  q.s   q.v3 -q.v2;
        -q.v2 -q.v3  q.s   q.v1;
        -q.v3  q.v2 -q.v1  q.s;
    ]
end

function VRmat(q::Quaternion)
    SA[
        q.v1  q.s  q.v3 -q.v2;
        q.v2 -q.v3  q.s  q.v1;
        q.v3  q.v2 -q.v1  q.s;
    ]
end

function VRᵀmat(q::Quaternion)
    SA[
        -q.v1  q.s -q.v3  q.v2;
        -q.v2  q.v3  q.s -q.v1;
        -q.v3 -q.v2  q.v1  q.s;
    ]
end

function LVᵀmat(q::Quaternion)
    SA[
        -q.v1 -q.v2 -q.v3;
         q.s -q.v3  q.v2;
         q.v3  q.s -q.v1;
        -q.v2  q.v1  q.s;
    ]
end

function LᵀVᵀmat(q::Quaternion)
    SA[
         q.v1  q.v2  q.v3;
         q.s  q.v3 -q.v2;
        -q.v3  q.s  q.v1;
         q.v2 -q.v1  q.s;
    ]
end

function RVᵀmat(q::Quaternion)
    SA[
        -q.v1 -q.v2 -q.v3;
         q.s  q.v3 -q.v2;
        -q.v3  q.s  q.v1;
         q.v2 -q.v1  q.s;
    ]
end

function RᵀVᵀmat(q::Quaternion)
    SA[
         q.v1  q.v2  q.v3;
         q.s -q.v3  q.v2;
         q.v3  q.s -q.v1;
        -q.v2  q.v1  q.s;
    ]
end

################################################################################
# Matrix-Vector Product Jacobian
################################################################################
function ∂VLmat∂q(p::AbstractVector) # 𝞉(VLmat(q)*p)/∂q
	SA[
    	0     p[1]  p[2]  p[3];
    	p[1]  0     p[3] -p[2];
    	p[2] -p[3]  0     p[1];
    	p[3]  p[2] -p[1]  0;
    ]
end

function ∂LVᵀmat∂q(p::AbstractVector) # 𝞉(∂LVᵀmat∂q(q)*p)/∂q
	SA[
    	0    -p[1] -p[2] -p[3];
    	p[1]  0     p[3] -p[2];
    	p[2] -p[3]  0     p[1];
    	p[3]  p[2] -p[1]  0;
    ]
end

function ∂VLᵀmat∂q(p::AbstractVector) # 𝞉(VLᵀmat(q)*p)/∂q
	SA[
		p[2] -p[1] -p[4]  p[3];
		p[3]  p[4] -p[1] -p[2];
		p[4] -p[3]  p[2] -p[1];
    ]
end

function ∂LᵀVᵀmat∂q(p::AbstractVector) # 𝞉(LᵀVᵀmat(q)*p)/∂q
	SA[
    	0     p[1]  p[2]  p[3];
    	p[1]  0    -p[3]  p[2];
    	p[2]  p[3]  0    -p[1];
    	p[3] -p[2]  p[1]  0;
    ]
end

function ∂VRmat∂q(p::AbstractVector) # 𝞉(VRmat(q)*p)/∂q
	SA[
		p[2]  p[1] -p[4]  p[3];
		p[3]  p[4]  p[1] -p[2];
		p[4] -p[3]  p[2]  p[1];
    ]
end

function ∂RᵀVᵀmat∂q(p::AbstractVector) # 𝞉(RᵀVᵀmat(q)*p)/∂q
	SA[
    	p[2]  p[1]  p[4] -p[3];
    	p[3] -p[4]  p[1]  p[2];
    	p[4]  p[3] -p[2]  p[1];
    ]
end

function ∂VRᵀmat∂q(p::AbstractVector) # 𝞉(RᵀVᵀmat(q)*p)/∂q
	SA[
    	p[2] -p[1]  p[4] -p[3];
    	p[3] -p[4] -p[1]  p[2];
    	p[4]  p[3] -p[2] -p[1];
    ]
end

function ∂Rᵀmat∂q(p::AbstractVector) # 𝞉(Rᵀmat(q)*p)/∂q
	SA[
    	p[1]  p[2]  p[3]  p[4];
    	p[2] -p[1]  p[4] -p[3];
    	p[3] -p[4] -p[1]  p[2];
    	p[4]  p[3] -p[2] -p[1];
    ]
end

function ∂Lmat∂q(p::AbstractVector) # 𝞉(Lmat(q)*p)/∂q
	SA[
    	p[1] -p[2] -p[3] -p[4];
    	p[2]  p[1]  p[4] -p[3];
    	p[3] -p[4]  p[1]  p[2];
    	p[4]  p[3] -p[2]  p[1];
    ]
end

function skew(p)
    SA[
    	 0    -p[3]  p[2];
    	 p[3]  0    -p[1];
    	-p[2]  p[1]  0;
    ]
end

function ∂skew∂p(λ) # 𝞉(skew(p)*λ)/∂p
	skew(-λ)
end


text:
# SPDX-License-Identifier: X11
# 2020-11-09 (Timed)
# Almost GCD (200pt)

function main()
  readline()
  A = parse.(Int, split(readline()))

  max,gcd = -1,-1
  for i ∈ 2:maximum(A)
    if (g = gcdness(i, A)) > gcd
      max = i
      gcd = g
    end
  end

  println(max)
end

function gcdness(k::Integer, A::Array{<:Integer})::Integer
  cnt = 0
  for a ∈ A
    a % k == 0 && (cnt += 1)
  end
  return cnt
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
mutable struct SpanningTree{V, E} <: AbstractGraph{V, E}
    graph::DirectedGraph{V, E}
    root::V
    edges::Vector{E}
    inedges::Vector{E}
    outedges::Vector{Vector{E}}
    edge_tree_indices::Vector{Int64} # mapping from edge_index(edge) to index into edges(tree)
end

# AbstractGraph interface
vertices(tree::SpanningTree) = vertices(tree.graph)
edges(tree::SpanningTree) = tree.edges
source(edge, tree::SpanningTree{V, E}) where {V, E} = source(edge, tree.graph) # note: doesn't check that edge is in spanning tree!
target(edge, tree::SpanningTree{V, E}) where {V, E} = target(edge, tree.graph) # note: doesn't check that edge is in spanning tree!
in_edges(vertex::V, tree::SpanningTree{V, E}) where {V, E} = (tree.inedges[vertex_index(vertex)],)
out_edges(vertex::V, tree::SpanningTree{V, E}) where {V, E} = tree.outedges[vertex_index(vertex)]

root(tree::SpanningTree) = tree.root
edge_to_parent(vertex::V, tree::SpanningTree{V, E}) where {V, E} = tree.inedges[vertex_index(vertex)]
edges_to_children(vertex::V, tree::SpanningTree{V, E}) where {V, E} = out_edges(vertex, tree)
tree_index(edge, tree::SpanningTree{V, E}) where {V, E} = tree.edge_tree_indices[edge_index(edge)]
tree_index(vertex::V, tree::SpanningTree{V, E}) where {V, E} = vertex === root(tree) ? 1 : tree_index(edge_to_parent(vertex, tree), tree) + 1

function SpanningTree(g::DirectedGraph{V, E}, root::V, edges::AbstractVector{E}) where {V, E}
    n = num_vertices(g)
    length(edges) == n - 1 || error("Expected n - 1 edges.")
    inedges = Vector{E}(undef, n)
    outedges = [E[] for i = 1 : n]
    edge_tree_indices = Int64[]
    treevertices = V[]
    for (i, edge) in enumerate(edges)
        parent = source(edge, g)
        child = target(edge, g)
        isempty(treevertices) && push!(treevertices, parent)
        @assert parent ∈ treevertices
        inedges[vertex_index(child)] = edge
        push!(outedges[vertex_index(parent)], edge)
        push!(treevertices, child)
        resize!(edge_tree_indices, max(edge_index(edge), length(edge_tree_indices)))
        edge_tree_indices[edge_index(edge)] = i
    end
    SpanningTree(g, root, edges, inedges, outedges, edge_tree_indices)
end

function SpanningTree(g::DirectedGraph{V, E}, root::V, flipped_edge_map::Union{AbstractDict, Nothing} = nothing;
        next_edge = first #= breadth first =#) where {V, E}
    tree_edges = E[]
    tree_vertices = [root]
    frontier = E[]
    append!(frontier, out_edges(root, g))
    append!(frontier, in_edges(root, g))
    while !isempty(frontier)
        # select a new edge
        e = next_edge(frontier)

        # remove edges from frontier
        flip = source(e, g) ∉ tree_vertices
        child = flip ? source(e, g) : target(e, g)
        filter!(x -> x ∉ in_edges(child, g), frontier)
        filter!(x -> x ∉ out_edges(child, g), frontier)

        # flip current edge if necessary
        if flip
            rewire!(g, e, target(e, g), source(e, g))
            newedge = flip_direction(e)
            replace_edge!(g, e, newedge)
            if flipped_edge_map !== nothing
                flipped_edge_map[e] = newedge
            end
            e = newedge
        end

        # update tree
        push!(tree_edges, e)
        push!(tree_vertices, child)

        # add new edges to frontier
        append!(frontier, x for x in out_edges(child, g) if target(x, g) ∉ tree_vertices)
        append!(frontier, x for x in in_edges(child, g) if source(x, g) ∉ tree_vertices)
    end
    length(tree_vertices) == num_vertices(g) || error("Graph is not connected.")
    SpanningTree(g, root, tree_edges)
end

# adds an edge and vertex to both the tree and the underlying graph
function add_edge!(tree::SpanningTree{V, E}, source::V, target::V, edge::E) where {V, E}
    @assert target ∉ vertices(tree)
    add_edge!(tree.graph, source, target, edge)

    push!(tree.edges, edge)
    push!(tree.inedges, edge)
    push!(tree.outedges, E[])
    push!(out_edges(source, tree), edge)
    resize!(tree.edge_tree_indices, max(edge_index(edge), length(tree.edge_tree_indices)))
    tree.edge_tree_indices[edge_index(edge)] = num_edges(tree)
    tree
end

function replace_edge!(tree::SpanningTree{V, E}, old_edge::E, new_edge::E) where {V, E}
    @assert old_edge ∈ edges(tree)
    src = source(old_edge, tree)
    dest = target(old_edge, tree)
    tree.edges[tree_index(old_edge, tree)] = new_edge
    tree.inedges[vertex_index(dest)] = new_edge
    out_edge_index = findfirst(isequal(old_edge), out_edges(src, tree))
    out_edges(src, tree)[out_edge_index] = new_edge
    replace_edge!(tree.graph, old_edge, new_edge)
    tree
end

function Base.show(io::IO, tree::SpanningTree, vertex = root(tree), level::Int64 = 0)
    for i = 1 : level print(io, "  ") end
    print(io, "Vertex: ")
    show(IOContext(io, :compact => true), vertex)
    if vertex == root(tree)
        print(io, " (root)")
    else
        print(io, ", ")
        print(io, "Edge: ")
        show(IOContext(io, :compact => true), edge_to_parent(vertex, tree))
    end
    for edge in edges_to_children(vertex, tree)
        print(io, "\n")
        child = target(edge, tree)
        show(io, tree, child, level + 1)
    end
end

function ancestors(vertex::V, tree::SpanningTree{V, E}) where {V, E}
    ret = [vertex]
    while vertex != root(tree)
        vertex = source(edge_to_parent(vertex, tree), tree)
        push!(ret, vertex)
    end
    ret
end

function lowest_common_ancestor(v1::V, v2::V, tree::SpanningTree{V, E}) where {V, E}
    while v1 != v2
        if tree_index(v1, tree) > tree_index(v2, tree)
            v1 = source(edge_to_parent(v1, tree), tree)
        else
            v2 = source(edge_to_parent(v2, tree), tree)
        end
    end
    v1
end


text:
module Bonobo

using DataStructures
using NamedTupleTools

"""
    AbstractNode

The abstract type for a tree node. Your own type for `Node` given to [`initialize`](@ref) needs to subtype it.
The default if you don't provide your own is [`DefaultNode`](@ref).
"""
abstract type AbstractNode end

"""
    AbstractSolution{Node<:AbstractNode, Value}

The abstract type for a `Solution` object. The default is [`DefaultSolution`](@ref).
It is parameterized by `Node` and `Value` where `Value` is the value which describes the full solution i.e the value for every variable.
"""
abstract type AbstractSolution{Node<:AbstractNode, Value} end

"""
    BnBNodeInfo

Holds the necessary information of every node.
This needs to be added by every `AbstractNode` as `std::BnBNodeInfo`

```julia
id :: Int
lb :: Float64
ub :: Float64
```
"""
mutable struct BnBNodeInfo
    id :: Int
    lb :: Float64
    ub :: Float64
end

"""
    DefaultNode <: AbstractNode

The default structure for saving node information.
Currently this includes only the necessary `std::BnBNodeInfo` which needs to be part of every [`AbstractNode`](@ref).
"""
mutable struct DefaultNode <: AbstractNode
    std :: BnBNodeInfo
end

"""
    DefaultSolution{Node<:AbstractNode,Value} <: AbstractSolution{Node, Value}

The default struct to save a solution of the branch and bound run.
It holds
```julia
objective :: Float64
solution  :: Value
node      :: Node
```
Both the `Value` and the `Node` type are determined by the [`initialize`](@ref) method.

`solution` holds the information to obtain the solution for example the values of all variables.
"""
mutable struct DefaultSolution{Node<:AbstractNode,Value} <: AbstractSolution{Node, Value}
    objective :: Float64
    solution  :: Value
    node      :: Node
end

"""
    AbstractTraverseStrategy

The abstract type for a traverse strategy.
If you implement a new traverse strategy this must be the supertype.

If you want to implement your own strategy the [`get_next_node`](@ref) function needs a new method
which dispatches on the `traverse_strategy` argument.
"""
abstract type AbstractTraverseStrategy end

"""
    AbstractBranchStrategy

The abstract type for a branching strategy.
If you implement a new branching strategy, this must be the supertype.

If you want to implement your own strategy, you must implement a new method for [`get_branching_variable`](@ref)
which dispatches on the `branch_strategy` argument.
"""
abstract type AbstractBranchStrategy end

"""
    BFS <: AbstractTraverseStrategy

The BFS traverse strategy always picks the node with the lowest bound first.
If there is a tie then the smallest node id is used as a tie breaker.
"""
struct BFS <: AbstractTraverseStrategy end

"""
    FIRST <: AbstractBranchStrategy

The `FIRST` strategy always picks the first variable which isn't fixed yet and can be branched on.
"""
struct FIRST <: AbstractBranchStrategy end

"""
    MOST_INFEASIBLE <: AbstractBranchStrategy

The `MOST_INFEASIBLE` strategy always picks the variable which is furthest away from being "fixed" and can be branched on.
"""
struct MOST_INFEASIBLE <: AbstractBranchStrategy end

mutable struct Options
    traverse_strategy   :: AbstractTraverseStrategy
    branch_strategy     :: AbstractBranchStrategy
    atol                :: Float64
    rtol                :: Float64
end

"""
    BnBTree{Node<:AbstractNode,Root,Value,Solution<:AbstractSolution{Node,Value}}

Holds all the information of the branch and bound tree.

```
incumbent::Float64 - The best objective value found so far. Is stores as problem is a minimization problem
incumbent_solution::Solution - The currently best solution object
lb::Float64        - The highest current lower bound
solutions::Vector{Solution} - A list of solutions
node_queue::PriorityQueue{Int,Tuple{Float64, Int}} - A priority queue with key being the node id and the priority consists of the node lower bound and the node id.
nodes::Dict{Int, Node}  - A dictionary of all nodes with key being the node id and value the actual node.
root::Root      - The root node see [`set_root!`](@ref)
branching_indices::Vector{Int} - The indices to be able to branch on used for [`get_branching_variable`](@ref)
num_nodes::Int  - The number of nodes created in total
sense::Symbol   - The objective sense: `:Max` or `:Min`.
options::Options  - All options for the branch and bound tree. See [`Options`](@ref).
```
"""
mutable struct BnBTree{Node<:AbstractNode,Root,Value,Solution<:AbstractSolution{Node,Value}}
    incumbent::Float64
    incumbent_solution::Union{Nothing,Solution}
    lb::Float64
    solutions::Vector{Solution}
    node_queue::PriorityQueue{Int,Tuple{Float64, Int}}
    nodes::Dict{Int, Node}
    root::Root
    branching_indices::Vector{Int}
    num_nodes::Int
    sense::Symbol
    options::Options
end

Base.broadcastable(x::BnBTree) = Ref(x)

include("util.jl")
include("node.jl")
include("branching.jl")

"""
    initialize(; kwargs...)

Initialize the branch and bound framework with the the following arguments.
Later it can be dispatched on `BnBTree{Node, Root, Solution}` for various methods.

# Keyword arguments
- `traverse_strategy` [`BFS`] currently the only supported traverse strategy is [`BFS`](@ref). Should be an [`AbstractTraverseStrategy`](@ref)
- `branch_strategy` [`FIRST`] currently the only supported branching strategies are [`FIRST`](@ref) and [`MOST_INFEASIBLE`](@ref). Should be an [`AbstractBranchStrategy`](@ref)
- `atol` [1e-6] the absolute tolerance to check whether a value is discrete
- `rtol` [1e-6] the relative tolerance to check whether a value is discrete
- `Node` [`DefaultNode`](@ref) can be special structure which is used to store all information about a node.
    - needs to have `AbstractNode` as the super type
    - needs to have `std :: BnBNodeInfo` as a field (see [`BnBNodeInfo`](@ref))
- `Solution` [`DefaultSolution`](@ref) stores the node and several other information about a solution
- `root` [`nothing`] the information about the root problem. The type can be used for dispatching on types
- `sense` [`:Min`] can be `:Min` or `:Max` depending on the objective sense
- `Value` [`Vector{Float64}`] the type of a solution

Return a [`BnBTree`](@ref) object which is the input for [`optimize!`](@ref).
"""
function initialize(;
    traverse_strategy = BFS(),
    branch_strategy = FIRST(),
    atol = 1e-6,
    rtol = 1e-6,
    Node = DefaultNode,
    Value = Vector{Float64},
    Solution = DefaultSolution{Node,Value},
    root = nothing,
    sense = :Min,
)
    return BnBTree{Node,typeof(root),Value,Solution}(
        Inf,
        nothing,
        -Inf,
        Vector{Solution}(),
        PriorityQueue{Int,Tuple{Float64, Int}}(),
        Dict{Int,Node}(),
        root,
        get_branching_indices(root),
        0,
        sense,
        Options(traverse_strategy, branch_strategy, atol, rtol)
    )
end

"""
    get_branching_indices(root)

Return a vector of variables to branch on from the current root object.
"""
function get_branching_indices end

"""
    optimize!(tree::BnBTree; callback=(args...; kwargs...)->())

Optimize the problem using a branch and bound approach.

The steps, repeated until terminated is true, are the following:
```julia
# 1. get the next open node depending on the traverse strategy
node = get_next_node(tree, tree.options.traverse_strategy)
# 2. evaluate the current node and return the lower and upper bound
# if the problem is infeasible both values should be set to NaN
lb, ub = evaluate_node!(tree, node)
# 3. update the upper and lower bound of the node struct
set_node_bound!(tree.sense, node, lb, ub)

# 4. update the best solution
updated = update_best_solution!(tree, node)
updated && bound!(tree, node.id)

# 5. remove the current node
close_node!(tree, node)
# 6. compute the node children and adds them to the tree
# internally calls get_branching_variable and branch_on_variable!
branch!(tree, node)
```

A `callback` function can be provided which will be called whenever a node is closed.
It always has the arguments `tree` and `node` and is called after the `node` is closed.
Additionally the callback function **must** accept additional keyword arguments (`kwargs`)
which are set in the following ways:
1. If the node is infeasible the kwarg `node_infeasible` is set to `true`.
2. If the node has a higher lower bound than the incumbent the kwarg `worse_than_incumbent` is set to `true`.
"""
function optimize!(tree::BnBTree; callback=(args...; kwargs...)->())
    while !terminated(tree)
        node = get_next_node(tree, tree.options.traverse_strategy)
        lb, ub = evaluate_node!(tree, node)
        # if the problem was infeasible we simply close the node and continue
        if isnan(lb) && isnan(ub)
            close_node!(tree, node)
            callback(tree, node; node_infeasible=true)
            continue
        end

        set_node_bound!(tree.sense, node, lb, ub)
        tree.node_queue[node.id] = (node.lb, node.id)
        _ , prio = peek(tree.node_queue)
        @assert tree.lb <= prio[1]
        tree.lb = prio[1]

        # if the evaluated lower bound is worse than the best incumbent -> close and continue
        if node.lb >= tree.incumbent
            close_node!(tree, node)
            callback(tree, node; worse_than_incumbent=true)
            continue
        end

        updated = update_best_solution!(tree, node)
        if updated
            bound!(tree, node.id)
            if isapprox(tree.incumbent, tree.lb; atol=tree.options.atol, rtol=tree.options.rtol)
                break
            end
        end

        close_node!(tree, node)
        branch!(tree, node)
        callback(tree, node)
    end
    sort_solutions!(tree.solutions, tree.sense)
end

"""
    sort_solutions!(solutions::Vector{<:AbstractSolution}, sense::Symbol)

Sort the solutions vector by objective value such that the best solution is at index 1.
"""
function sort_solutions!(solutions::Vector{<:AbstractSolution}, sense::Symbol)
    sort!(solutions; rev=:sense == :Max, by=s->s.objective)
end

"""
    terminated(tree::BnBTree)

Return true when the branch and bound loop in [`optimize!`](@ref) should be terminated.
Default behavior is to terminate the loop only when no nodes exist in the priority queue.
"""
terminated(tree::BnBTree) = isempty(tree.nodes)

"""
    set_node_bound!(objective_sense::Symbol, node::AbstractNode, lb, ub)

Set the bounds of the `node` object to the lower and upper bound given.
Internally everything is stored as a minimization problem. Therefore the objective_sense `:Min`/`:Max` is needed.
"""
function set_node_bound!(objective_sense::Symbol, node::AbstractNode, lb, ub)
    if isnan(ub)
        ub = Inf
    end
    if objective_sense == :Min
        node.lb = max(lb, node.lb)
        node.ub = ub
    else
        node.lb = max(-lb, node.lb)
        node.ub = -ub
    end
end

"""
    bound!(tree::BnBTree, current_node_id::Int)

Close all nodes which have a lower bound higher or equal to the incumbent
"""
function bound!(tree::BnBTree, current_node_id::Int)
    for (_, node) in tree.nodes
        if node.id != current_node_id && node.lb >= tree.incumbent
            close_node!(tree, node)
        end
    end
end

"""
    close_node!(tree::BnBTree, node::AbstractNode)

Delete the node from the nodes dictionary and the priority queue.
"""
function close_node!(tree::BnBTree, node::AbstractNode)
    delete!(tree.nodes, node.id)
    delete!(tree.node_queue, node.id)
end

"""
    update_best_solution!(tree::BnBTree, node::AbstractNode)

Update the best solution when we found a better incumbent.
Calls [`add_new_solution!`] if this is the case, returns whether a solution was added.
"""
function update_best_solution!(tree::BnBTree, node::AbstractNode)
    isinf(node.ub) && return false
    node.ub >= tree.incumbent && return false

    tree.incumbent = node.ub

    add_new_solution!(tree, node)
    return true
end

"""
    add_new_solution!(tree::BnBTree{N,R,V,S}, node::AbstractNode) where {N,R,V,S<:DefaultSolution{N,V}}

Currently it changes the general solution itself by calling [`get_relaxed_values`](@ref) which needs to be implemented by you.

This function needs to be implemented by you if you have a different type of Solution object than [`DefaultSolution`](@ref).
"""
function add_new_solution!(tree::BnBTree{N,R,V,S}, node::AbstractNode) where {N,R,V,S<:DefaultSolution{N,V}}
    sol = DefaultSolution(node.ub, get_relaxed_values(tree, node), node)
    push!(tree.solutions, sol)
    if tree.incumbent_solution === nothing || sol.objective < tree.incumbent_solution.objective
        tree.incumbent_solution = sol
    end
end

"""
    get_relaxed_values(tree::BnBTree, node::AbstractNode)

Get the values of the current node. This is always called only after [`evaluate_node!`](@ref) is called.
It is used to store a `Solution` object.
Return the type of `Value` given to the [`initialize`](@ref) method.
"""
function get_relaxed_values end

"""
    get_solution(tree::BnBTree; result=1)

Return the solution values of the problem.
See [`get_objective_value`](@ref) to obtain the objective value instead.
"""
function get_solution(tree::BnBTree{N,R,V,S}; result=1) where {N,R,V,S<:DefaultSolution{N,V}}
    return tree.solutions[result].solution
end

"""
    get_objective_value(tree::BnBTree; result=1)

Return the objective value
"""
function get_objective_value(tree::BnBTree{N,R,V,S}; result=1) where {N,R,V,S<:DefaultSolution{N,V}}
    if tree.sense == :Max
        return -tree.solutions[result].objective
    else
        return tree.solutions[result].objective
    end
end

"""
    get_num_solutions(tree::BnBTree)

Return the number of solutions available.
"""
get_num_solutions(tree::BnBTree) = length(tree.solutions)

export BnBTree, BnBNodeInfo, AbstractNode, AbstractSolution

export AbstractTraverseStrategy, AbstractBranchStrategy

end


text:
module Subgroups

"""
Generate a sub-group of the multiplicative group of integers modulo p
from powers of x.
"""
function generate(x::Integer, p::Integer)
    Set(BigInt[x^i%p for i=1:p])
end

end

max_stars_repo_name: drichardson/examples

max_issues_repo_name: drichardson/examples

max_forks_repo_name: drichardson/examples


text:
using DiffEqSensitivity,OrdinaryDiffEq, ParameterizedFunctions,
      RecursiveArrayTools, DiffEqBase, ForwardDiff, Calculus
using Test
using DiffEqSensitivity: SensitivityAlg

f = @ode_def LotkaVolterra begin
  dx = a*x - b*x*y
  dy = -c*y + x*y
end a b c

function LotkaVolt!(du, u, p, t)
  du[1] = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = -p[3]*u[2] + u[1]*u[2]
end

p = [1.5,1.0,3.0]
prob = ODELocalSensitivityProblem(f,[1.0;1.0],(0.0,10.0),p)
probInpl = ODELocalSensitivityProblem(LotkaVolt!,[1.0;1.0],(0.0,10.0),p)
probnoad = ODELocalSensitivityProblem(LotkaVolt!,[1.0;1.0],(0.0,10.0),p,
                                      SensitivityAlg(autodiff=false))
sol = solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14)
@test_broken solInpl = solve(probInpl,KenCarp4(),abstol=1e-14,reltol=1e-14)
@test_broken solInpl2 = solve(probInpl,Rodas4(autodiff=false),abstol=1e-14,reltol=1e-14)
solInpl = solve(probInpl,KenCarp4(autodiff=false),abstol=1e-14,reltol=1e-14)
solInpl2 = solve(probInpl,Rodas4(autodiff=false),abstol=1e-14,reltol=1e-14)
solnoad = solve(probnoad,KenCarp4(autodiff=false),abstol=1e-14,reltol=1e-14)

x = sol[1:sol.prob.f.numindvar,:]

@test sol(5.0) ≈ solnoad(5.0)
@test sol(5.0) ≈ solInpl(5.0)
@test solInpl(5.0) ≈ solInpl2(5.0)

# Get the sensitivities

da = sol[sol.prob.f.numindvar+1:sol.prob.f.numindvar*2,:]
db = sol[sol.prob.f.numindvar*2+1:sol.prob.f.numindvar*3,:]
dc = sol[sol.prob.f.numindvar*3+1:sol.prob.f.numindvar*4,:]

sense_res1 = [da[:,end] db[:,end] dc[:,end]]

prob = ODELocalSensitivityProblem(f.f,[1.0;1.0],(0.0,10.0),p,SensitivityAlg(autojacvec=true))
sol = solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14)
x = sol[1:sol.prob.f.numindvar,:]

# Get the sensitivities

da = sol[sol.prob.f.numindvar+1:sol.prob.f.numindvar*2,:]
db = sol[sol.prob.f.numindvar*2+1:sol.prob.f.numindvar*3,:]
dc = sol[sol.prob.f.numindvar*3+1:sol.prob.f.numindvar*4,:]

sense_res2 = [da[:,end] db[:,end] dc[:,end]]

function test_f(p)
  prob = ODEProblem(f,eltype(p).([1.0,1.0]),(0.0,10.0),p)
  solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14,save_everystep=false)[end]
end

p = [1.5,1.0,3.0]
fd_res = ForwardDiff.jacobian(test_f,p)
calc_res = Calculus.finite_difference_jacobian(test_f,p)

@test sense_res1 ≈ sense_res2 ≈ fd_res
@test sense_res1 ≈ sense_res2 ≈ calc_res

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

# Now do from a plain parameterized function

function f2(du,u,p,t)
  du[1] = p[1] * u[1] - p[2] * u[1]*u[2]
  du[2] = -p[3] * u[2] + u[1]*u[2]
end
p = [1.5,1.0,3.0]
prob = ODELocalSensitivityProblem(f2,[1.0;1.0],(0.0,10.0),p)
sol = solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14)
res = sol[1:sol.prob.f.numindvar,:]

# Get the sensitivities

da = sol[sol.prob.f.numindvar+1:sol.prob.f.numindvar*2,:]
db = sol[sol.prob.f.numindvar*2+1:sol.prob.f.numindvar*3,:]
dc = sol[sol.prob.f.numindvar*3+1:sol.prob.f.numindvar*4,:]

sense_res = [da[:,end] db[:,end] dc[:,end]]

p = [1.5,1.0,3.0]
fd_res = ForwardDiff.jacobian(test_f,p)
calc_res = Calculus.finite_difference_jacobian(test_f,p)

@test sense_res ≈ fd_res
@test sense_res ≈ calc_res

## Check extraction

x,dp = extract_local_sensitivities(sol)
x == res
dp[1] == da

x,dp = extract_local_sensitivities(sol,length(sol.t))
sense_res2 = hcat(dp...)
sense_res == sense_res2

extract_local_sensitivities(sol,sol.t[3]) == extract_local_sensitivities(sol,3)

tmp = similar(sol[1])
extract_local_sensitivities(tmp,sol,sol.t[3]) == extract_local_sensitivities(sol,3)


text:
# SolarIrradiance.jl
# Copyright (c) 2021 Pierre Galloo-Beauvais, Pierre Haessig
# This code is subject to the terms of the MIT license (see LICENSE.txt)

module GHIProcess

export GHI_download, GHI_read, sliceday
export forwardfill, backwardfill, centeredfill, interpolation
export GHIavg

using CSV
using DataFrames
using Dates
using HTTP

"""
    GHI_download(year, lat, lon, basename="GHI")

Download of GHI hourly data from PVGIS at location `lat`, `lon` (°) for `year`.

Data is saved in a CSV file with the `basename` prefix (default: `"GHI"`).

⚠ `lat`, `lon` must be expressed with maximum 3 significant digits,
otherwise the PVGIS API refuses to generate the CSV containing the data

See also: [`GHI_read`](@ref) to load the data.
"""
function GHI_download(year, lat, lon; basename="GHI")
    baseurl = "https://re.jrc.ec.europa.eu/api/seriescalc?"
    angle = 0
    aspect = 0
    url = string(baseurl, "lat=",lat, "&lon=",lon, "&angle=", angle, "&aspect=", aspect,
                 "&startyear=", year, "&endyear=", year,
                 "&outputformat=csv&browser=1")
    fname = string(basename, "_", lat, "_", lon, "_SA_", angle, "deg_", aspect, "deg_",
                   year, "_", year, ".csv")

    io = open(fname, "w")
    r = HTTP.request("GET", url, response_stream=io)
    close(io)
    print("PVGIS GHI data downloaded and saved as $fname\n")
    return fname
end

"""
    GHI_read(source, year)

Load the CSV `source` containing PVGIS GHI hourly data for `year`.

Return a Vector of length 8760 or 8784 depending on whether `year` is a leap year.

See also: [`GHI_download`](@ref) to download the data.

# Example:

Download 2012 data for Rennes, France:

```jldoctest
julia> year = 2012; lon = -1.678; lat = 48.117;

julia> fname = GHI_download(year, lat, lon; basename="GHI")
PVGIS GHI data downloaded and saved as GHI_48.117_-1.678_SA_0deg_0deg_2012_2012.csv
"GHI_48.117_-1.678_SA_0deg_0deg_2012_2012.csv"

julia> GHI_year = GHI_read(fname, year)
8784-element Vector{Float64}:
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
 70.0
 97.0
  ⋮
 85.0
 49.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
```
"""
function GHI_read(source, year)
    nhours = 24 * (365 + Dates.isleapyear(year))
    df = CSV.read(source, DataFrame; header=9, select=[2], limit=nhours)
    GHI_year = df[!,1]
    return GHI_year
end

"""
    sliceday(GHI_year,n)

Renvoie une liste de 24 éléments, correspondant à la GHI d'un jour n d'une année.

La sélection de l'année se fait via l'utilisation préalable de la fonction GHIannee(annee,phi_deg,lambda_deg), qui permet de récupérer la GHI pour toute une année. Il faut donc l'utiliser en amont, et stocker le résultat dans une liste, nommée ici "GHI_year".

#Exemple :

julia> GHI_year = GHIannee(2012,48.117,-1.678)
julia> sliceday(GHI_year,1)

24-element Vector{Any}:
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
 70.0
 97.0
 55.0
 68.0
 83.0
 58.0
 15.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0

"""
function sliceday(GHI_year,n)
    GHI_day = GHI_year[1+(n-1)*24:n*24]
    return GHI_day
end

"""
    forwardfill(lesgjour,heure)

La fonction forwardfill renvoie la valeur correspondant au GHI (moyen sur une minute) à une certaine heure, en réalisant une interpolation type bloqueur d'ordre 0 aux données d'entrées. La valeur de GHI à l'heure demandée correspond à la valeur de la mesure précédant directement cette heure.

La fonction prend en entrée une liste correspondant aux données de GHI d'une journée, sous forme d'une liste de 24 éléments, et une heure, qui peut être un flottant.

Elle doit donc être précédée de l'utilisation de la fonction GHIannee et de la fonction sliceday.

#Exemple :

julia>GHI_year=GHIannee(2012,48.117,-1.678)
julia>lesgjour=sliceday(GHI_year,1)
julia>forwardfill(lesgjour,12.5)

68.0

"""
function forwardfill(lesgjour,heure)

    if heure>=0 && heure<1+1/6
        return(lesgjour[1])
    elseif heure>=23+1/6 && heure<=24
        return(lesgjour[24])
    end

    for k in [1:22...]
        if heure>=k+1/6 && heure<k+1+1/6
            return(lesgjour[k+1])
        end
    end
end

"""
    backwardfill(lesgjour,heure)

La fonction backwardfill renvoie la valeur correspondant au GHI (moyen sur une minute) à une certaine heure, en réalisant une interpolation type bloqueur d'ordre 0 inversé aux données d'entrées. La valeur de GHI à l'heure demandée correspond à la valeur de la mesure suivant directement cette heure.

La fonction prend en entrée une liste correspondant aux données de GHI d'une journée, sous forme d'une liste de 24 éléments, et une heure, qui peut être un flottant.

Elle doit donc être précédée de l'utilisation de la fonction GHIannee et de la fonction sliceday.

#Exemple :

julia>GHI_year=GHIannee(2012,48.117,-1.678)
julia>lesgjour=sliceday(GHI_year,1)
julia>backwardfill(lesgjour,12.5)

83.0

"""
function backwardfill(lesgjour,heure)

    if heure>=0 && heure<=1/6
        return(lesgjour[1])
    elseif heure>=22+1/6 && heure<=24
        return(lesgjour[24])
    end

    for k in [0:21...]
        if heure>k+1/6 && heure<=k+1+1/6
            return(lesgjour[k+2])
        end
    end
end

"""
    centeredfill(lesgjour,heure)

La fonction centeredfill renvoie la valeur correspondant au GHI (moyen sur une minute) à une certaine heure, en réalisant une interpolation type bloqueur centré sur la mesure. La valeur de GHI à l'heure demandée correspond à la valeur de la mesure la plus proche de cette heure.

La fonction prend en entrée une liste correspondant aux données de GHI d'une journée, sous forme d'une liste de 24 éléments, et une heure, qui peut être un flottant.

Elle doit donc être précédée de l'utilisation de la fonction GHIannee et de la fonction sliceday.

#Exemple :

julia>GHI_year=GHIannee(2012,48.117,-1.678)
julia>lesgjour=sliceday(GHI_year,1)
julia>centeredfill(lesgjour,12.5)

68.0

julia>GHI_year=GHIannee(2012,48.117,-1.678)
julia>lesgjour=sliceday(GHI_year,1)
julia>centeredfill(lesgjour,12.8)

83.0

"""
function centeredfill(lesgjour,heure)

    if heure>=0 && heure<4/6
        return(lesgjour[1])
    elseif heure>=23+4/6 && heure<=24
        return(lesgjour[24])
    end

    for k in [0:22...]
        if heure>=k+4/6 && heure<k+1+4/6
            return(lesgjour[k+2])
        end
    end
end

"""
    interpolation(lesgjour,heure)

La fonction interpolation renvoie la valeur correspondant au GHI (moyen sur une minute) à une certaine heure, en réalisant une interpolation linéaire sur les données d'entrées.

La fonction prend en entrée une liste correspondant aux données de GHI d'une journée, sous forme d'une liste de 24 éléments, et une heure, qui peut être un flottant.

Elle doit donc être précédée de l'utilisation de la fonction GHIannee et de la fonction sliceday.

#Exemple :

julia>GHI_year=GHIannee(2012,48.117,-1.678)
julia>lesgjour=sliceday(GHI_year,1)
julia>interpolation(lesgjour,12.5)

77.50000000000001

"""
function interpolation(lesgjour,heure)
    if heure>=0 && heure<=1/6
        return(heure*lesgjour[1])
    elseif heure>23+1/6 && heure<=24
        return(lesgjour[24]+(heure-23-1/6)*(lesgjour[24]-lesgjour[23]))
    end

    for k in [0:22...]
        if heure>k+1/6 && heure<=k+1+1/6
            return(lesgjour[k+1]+(heure-k-1/6)*(lesgjour[k+2]-lesgjour[k+1]))
        end
    end
end

"""
    GHIavg(GHI_year,n,tcini,dt,fonctionGHI)

Cette fonction évalue la valeur moyenne du GHI sur un intervalle de temps donné, un jour donné d'une année donnée, avec le choix d'une interpolation sur les données du jour.

Cette fonction prend en entrée les données de GHI sur une année, obtenues via la fonction GHIannee.

tcini est l'heure initiale, le début de l'intervalle sur lequel on veut obtenir la valeur moyenne de GHI. Peut être un flottant. Ne peut pas dépasser 24.

dt est la durée de l'intervalle sur lequel on veut obtenir la valeur moyenne de GHI, en heures. Peut être un flottant. Ne peut pas dépasser 24.

fonctionGHI est une chaîne de caractères correspondant à la signature de la fonction d'interpolation choisie pour reconstruire les données de GHI. On peut donc utiliser les chaînes suivantes : "forwardfill", "backwardfill", "centeredfill", "interpolation".

# Exemple :

julia> GHI_year=GHIannee(2012,48.117,-1.678)
julia> GHIavg(GHI_year,1,12.5,0.5,interpolation)

73.0

julia> GHI_year=GHIannee(2012,48.117,-1.678)
julia> GHIavg(GHI_year,1,12.5,0.5,forwardfill)

68.0

"""
function GHIavg(GHI_year,n,tcini,dt,fonctionGHI)
    lesgjour=sliceday(GHI_year,n)

    G=0
    compteur=60*dt
    for k in range(tcini,length=convert(Int64,60*dt))
        G+=fonctionGHI(lesgjour,tcini)
    end
    return(G/compteur)
end

end # module