text:
"""
`module Build` : this is a very poor man's version of `ase.build`. At the moment
only a subset of the `bulk` functionality is supported.
"""
module Build

import JuLIP
using ..Chemistry
using JuLIP: JVec, JMat, JVecF, JMatF, mat, vecs,
      Atoms, cell, cell_vecs, positions, momenta, masses, numbers, pbc,
      chemical_symbols, set_cell!, set_pbc!, update_data!,
      apply_defm!, calculator, set_calculator!

using LinearAlgebra: I, Diagonal, isdiag, norm, cross, det

import Base: union
import JuLIP: Atoms

export repeat, bulk, cluster, autocell!, append, rotation_matrix, rotate!

_auto_pbc(pbc::Tuple{Bool, Bool, Bool}) = JVec(pbc...)
_auto_pbc(pbc::Bool) = JVec(pbc, pbc, pbc)
_auto_pbc(pbc::AbstractVector) = JVec(pbc...)

_auto_cell(cell) = cell
_auto_cell(cell::AbstractMatrix) = JMat(cell)
_auto_cell(C::AbstractVector{T}) where {T <: Real} = (
   length(C) == 3 ? JMatF(C[1], 0.0, 0.0, 0.0, C[2], 0.0, 0.0, 0.0, C[3]) :
                    JMatF(C...) )  # (case 2 requires that length(C) == 0
_auto_cell(C::AbstractVector{T}) where {T <: AbstractVector} = JMatF([ C[1] C[2] C[3] ])
_auto_cell(C::AbstractVector) = _auto_cell([ c for c in C ])

"""
`autocell!(at::Atoms) -> Atoms`

generates cell vectors that contain all atom positions + a small buffer,
sets the cell of `at` accordingly (in-place) and returns the same atoms objet
with the new cell.
"""
autocell!(at::Atoms) = set_cell!(at, _autocell(positions(at)))

function _autocell(X::AbstractVector{JVec{T}}) where T
   ext = extrema(mat(X), dims = (2,))
   C = Diagonal([ e[2] - e[1] + 1.0  for e in ext ][:])
   return JMat{T}(C)
end

# if we have no clue about X just return it and see what happens
_auto_X(X) = X
# if the elements of X weren't inferred, try to infer before reading
# TODO: this might lead to a stack overflow!!
_auto_X(X::AbstractVector) = _auto_X( [x for x in X] )
# the cases where we know what to do ...
_auto_X(X::AbstractVector{T}) where {T <: AbstractVector} = [ JVecF(x) for x in X ]
_auto_X(X::AbstractVector{T}) where {T <: Real} = _auto_X(reshape(X, 3, :))
_auto_X(X::AbstractMatrix) = (@assert size(X)[1] == 3;
                              [ JVecF(X[:,n]) for n = 1:size(X,2) ])

_auto_M(M::AbstractVector{T}) where {T <: AbstractFloat} = Vector{T}(M)
_auto_M(M::AbstractVector) = Vector{Float64}(M)

_auto_Z(Z::AbstractVector) = Vector{AtomicNumber}(Z)

const PositionArray{T} = Union{AbstractMatrix{T}, AbstractVector{JVec{T}}}

Atoms(; kwargs...) = Atoms{Float64}(; kwargs...)

Atoms{T}(; X = JVec{T}[],
           P = JVec{T}[],
           M = T[],
           Z = AtomicNumber[],
           cell = zero(JMat{T}),
           pbc = JVec(false, false, false),
           calc = nothing ) where {T} =
      Atoms(X, P, M, Z, cell, pbc, calc)

Atoms(X, P, M, Z, cell, pbc, calc=nothing) =
      Atoms(_auto_X(X), _auto_X(P), _auto_M(M), _auto_Z(Z),
            _auto_cell(cell), _auto_pbc(pbc), calc)

Atoms(Z::Vector{AtomicNumber}, X::Vector{JVec{T}}; kwargs...) where {T} =
  Atoms{T}(; Z=Z, X=X, kwargs...)

Atoms( Z::Vector{<: Integer}, X::PositionArray{T}; kwargs...
       ) where {T <: AbstractFloat} =
   Atoms{T}(; Z=Z, X=X, kwargs...)

"""
simple way to construct an atoms object from just positions
"""
Atoms(s::Symbol, X::PositionArray{T}; kwargs...) where {T <: AbstractFloat} =
      Atoms{T}(; X = X,
                 M = fill(one(T), length(X)),
                 P = zeros(JVec{T}, length(X)),
                 Z = fill(atomic_number(s), length(X)),
                 cell = _autocell(X),
                 pbc = (false, false, false), kwargs... )

Atoms(s::Symbol, X::Matrix{Float64}) = Atoms(s, vecs(X))

# shallow copy constructor
Atoms(at::Atoms) = Atoms(at.X, at.P, at.M, at.Z, at.cell, at.pbc, at.calc)

import Base.copy

"""
   at2 = copy(at)

Return a copy of Atoms, referring to same arrays for positions, momenta, etc.
Use `deepcopy()` for a deep copy of the arrays as well.
"""
copy(at::Atoms) = Atoms(at)

"""
   rotation_matrix([x=X, y=Y, z=Z])

Construct a rotation matrix from two or more Miller indicies.

Example usage:

```
A = rotation_matrix(x=[1,1,1], y=[1,-1,0]) # fills in z from cross(x, y)
```
"""
function rotation_matrix(; x=nothing, y=nothing, z=nothing)
   @assert sum([x, y, z] .!== nothing) >= 2

   x !== nothing && (x /= norm(x))
   y !== nothing && (y /= norm(y))
   z !== nothing && (z /= norm(z))

   if x === nothing
      @assert isapprox(y' * z, 0, atol=1e-6)
      x = cross(y, z)
      x /= norm(x)
   end

   if y === nothing
      @assert isapprox(z' * x, 0, atol=1e-6)
      y = cross(z, x)
      y /= norm(y)
   end

   if z === nothing
      @assert isapprox.(x' * y, 0, atol=1e-6)
      z = cross(x, y)
      z /= norm(z)
   end

   A = hcat(x, y, z)
   (det(A) < 0.0) && (A = hcat(-x, y, z))
   return A
end

"""
Rotate atoms in place to align axes with Miller indices x, y, z, at least two of which must be given.
"""
rotate!(at::Atoms; x=nothing, y=nothing, z=nothing) = apply_defm!(at, rotation_matrix(x=x, y=y, z=z))

const _unit_cells = Dict(     # (positions, cell matrix, factor of a)
   :fcc => ( [ [0 0 0], ],
             [0 1 1; 1 0 1; 1 1 0],  0.5),
   :bcc => ( [ [0.0,0.0,0.0], ],
             [-1 1 1; 1 -1 1; 1 1 -1], 0.5),
   :diamond => ( [ [0.0, 0.0, 0.0], [0.5, 0.5, 0.5] ],
                 [0 1 1; 1 0 1; 1 1 0], 0.5)
)

const _cubic_cells = Dict(   # (positions, factor of a)
   :fcc => ( [ [0 0 0], [0 1 1], [1 0 1], [1 1 0] ], 0.5 ),
   :bcc => ( [ [0 0 0], [1 1 1] ], 0.5 ),
   :diamond => ( [ [0 0 0], [1 1 1], [0 2 2], [1 3 3], [2 0 2],
                   [3 1 3], [2 2 0], [3 3 1] ], 0.25 )
)

_simple_structures = [:fcc, :bcc, :diamond]

function _simple_bulk(sym::Symbol, cubic::Bool; a=nothing)
   if cubic
      X, scal = _cubic_cells[symmetry(sym)]
      C = Matrix(1.0I, 3, 3) / scal
   else
      X, C, scal = _unit_cells[symmetry(sym)]
   end
   a === nothing && (a = lattice_constant(sym))
   return [ JVecF(x) * a * scal  for x in X ], JMatF(C * a * scal)
end

function _bulk_hcp(sym::Symbol; a=nothing, c=nothing)
   D = Chemistry.refstate(sym)
   a === nothing && (a = D["a"])
   c === nothing && (c = a * D["c/a"])
   return [JVecF(0.0, 0.0, 0.0), JVecF(0.0, a / sqrt(3), c / 2)],
          JMatF( [a, -a / 2,  0.0, 0.0, a * sqrt(3)/2, 0.0, 0.0, 0.0, c] )
end

_convert_pbc(pbc::NTuple{3, Bool}) = pbc
_convert_pbc(pbc::Bool) = (pbc, pbc, pbc)

function bulk(sym::Symbol; T=Float64, cubic = false, pbc = (true,true,true), a=nothing, c=nothing, x=nothing, y=nothing, z=nothing)
   symm = symmetry(sym)
   if symm in _simple_structures
      X, C = _simple_bulk(sym, cubic; a=a)
   elseif symm == :hcp
      X, C = _bulk_hcp(sym; a=a, c=c)  # cubic parameter is irrelevant for hcp
   end
   m = atomic_mass(sym)
   Z = atomic_number(sym)
   nat = length(X)
   at = Atoms( convert(Vector{JVec{T}}, X),
                 fill(zero(JVec{T}), nat),
                 fill(T(m), nat),
                 fill(AtomicNumber(Z), nat),
                 convert(JMat{T}, C),
                 _convert_pbc(pbc)  )
   (x !== nothing || y !== nothing || z !== nothing) && rotate!(at, x=x, y=y, z=z)
   return at
end

"""
auxiliary function to convert a general cell to a cubic one; this is a bit of
a hack, so to not make a mess of things, this will only work in very restrictive
circumstances and otherwise throw an error. But it could be revisited.
"""
function _cubic_cell(atu::Atoms{T}) where {T}
   @assert length(atu) == 1
   ru = JuLIP.rmin(atu)
   at = bulk(chemical_symbol(atu.Z[1]), cubic=true)
   r = JuLIP.rmin(at)
   return apply_defm!(at, (ru/r) * one(JMat{T}))
end

"""
`cluster(args...; kwargs...) -> at::AbstractAtoms`

Produce a cluster of approximately radius R. The center
atom is always at index 1.

## Methods
```
cluster(atu::AbstractAtoms, R::Real)  # specify unit cell
cluster(sym::Symbol, R::Real)         # atu = bulk(sym; cubic=true)
```
Both methods assume that there is only a single species, and both require
the use of an orthorhombic unit cell (for now).

* `sym`: chemical symbol of the atomic species
* `R` : length-scale, meaning depends on shape of the cluster

## Keyword Arguments:
* `dims` : dimension into which the cluster is extended, typically
   `[1,2,3]` for 3D point defects and `[1,2]` for 2D dislocations, in the
   remaining dimension(s) the b.c. will be periodic.
* `shape` : shape of the cluster, at the moment only `:ball` is allowed
* `parity`: enforce an odd (where `parity=1`) or even (where `parity=0`)
   number of cells in each lattice direction. Use a `nothing` element to
   leave the number of cells unchanged (e.g. in periodic directions).
## todo
 * lift the restriction of single species
 * allow other shapes
"""
function cluster(atu::Atoms{T}, R::Real;
                 dims = findall(pbc(atu).==true),
                 shape = :ball, x0=nothing,
                 parity=nothing) where {T}
   sym = chemical_symbols(atu)[1]
   # check that the cell is orthorombic
   if !isdiag(cell(atu))
      if length(atu) > 1
         error("""`JuLIP.cluster` requires as argument either a cubic cell or a
                  one-atom cell.""")
      end
      atu = _cubic_cell(atu)
   end
   # # check that the first index is the centre
   # @assert norm(atu[1]) == 0.0
   # determine by how much to multiply in each direction
   Fu = cell(atu)'
   L = [ j ∈ dims ? 2 * (ceil(Int, R/Fu[j,j])+3) : 1    for j = 1:3]
   if parity !== nothing
      L[((L .% 2) .!= parity) .& (parity .!= nothing)] .+= 1
   end
   # multiply
   at = atu * L
   # find point closest to centre
   x̄ = sum( x[dims] for x in at.X ) / length(at.X)
   i0 = findmin( [norm(x[dims] - x̄) for x in at.X] )[2]
   if x0 === nothing
      x0 = at[i0]
   else
      x0 += at[i0]
   end
   # swap positions
   X = positions(at)
   X[1], X[i0] = X[i0], X[1]
   F = Diagonal([Fu[j,j]*L[j] for j = 1:3])
   # carve out a cluster with mini-buffer to account for round-off
   r = [ norm(x[dims] - x0[dims]) for x in X ]
   IR = findall( r .<= R+sqrt(eps(T)) )
   # generate new positions
   Xcluster = X[IR]
   # generate the cluster
   at_cluster = Atoms(sym, Xcluster)
   set_cell!(at_cluster, F')
   # take open boundary in all directions specified in dims, periodic in the rest
   set_pbc!(at_cluster, [ !(j ∈ dims) for j = 1:3 ])
   # return the cluster and the index of the center atom
   return at_cluster
end

cluster(sym::Symbol, R::Real; kwargs...) =
      cluster( bulk(sym, cubic=true), R; kwargs... )

"""
```
repeat(at::Atoms, n::NTuple{3}) -> Atoms
repeat(at::Atoms, n::Integer) = repeat(at, (n,n,n))
```

Takes an `Atoms` configuration / cell and repeats it n_j times
into the j-th cell-vector direction. For example,
```
at = repeat(bulk("C"), (3,2,4))
```
creates 3 x 2 x 4 unit cells of carbon.

The same can be achieved by `*`:
```
at = bulk("C") * (3, 2, 4)
```
"""
function Base.repeat(at::Atoms, n::NTuple{3})
   C0 = cell(at)
   c1, c2, c3 = cell_vecs(at)
   X0, P0, M0, Z0 = positions(at), momenta(at), masses(at), numbers(at)
   nrep = n[1] * n[2] * n[3]
   nat0 = length(at)
   X = Base.repeat(X0, outer=nrep)
   P = Base.repeat(P0, outer=nrep)
   M = Base.repeat(M0, outer=nrep)
   Z = Base.repeat(Z0, outer=nrep)
   i = 0
   for a in CartesianIndices( (1:n[1], 1:n[2], 1:n[3]) )
      b = c1 * (a[1]-1) + c2 * (a[2]-1) + c3 * (a[3]-1)
      X[i+1:i+nat0] = [b+x for x in X0]
      i += nat0
   end
   at1 = Atoms(X, P, M, Z, Diagonal(collect(n)) * cell(at), pbc(at))
   set_calculator!(at1, calculator(at))
   return at1
end

Base.repeat(at::Atoms, n::Integer) = repeat(at, (n,n,n))
Base.repeat(at::Atoms, n::AbstractArray) = repeat(at, (n...,))

import Base.*
*(at::Atoms, n) = repeat(at, n)
*(n, at::Atoms) = repeat(at, n)

"""
`deleteat!(at::Atoms, n) -> at`:

returns the same atoms object `at`, but with the atom(s) specified by `n`
removed.
"""
function Base.deleteat!(at::Atoms, n)
   deleteat!(at.X, n)
   deleteat!(at.P, n)
   deleteat!(at.M, n)
   deleteat!(at.Z, n)
   update_data!(at, Inf)
   JuLIP.reset_clamp!(at)
   return at
end

union(at1::Atoms, at2::Atoms) =
   Atoms( X = union(at1.X, at2.X),
          P = union(at1.P, at2.P),
          M = union(at1.M, at2.M),
          Z = union(at1.Z, at2.Z),
          cell = cell(at1),
          pbc = pbc(at1) )

append(at::Atoms, X::AbstractVector{<:JVec}) =
   Atoms( X = union(at.X, X),
          P = union(at.P, zeros(JVecF, length(X))),
          M = union(at.M, zeros(length(X))),
          Z = union(at.Z, zeros(Int, length(X))),
          cell = cell(at),
          pbc = pbc(at) )

end


text:
using Test
using ContinuousWavelet

@testset "ContinuousWavelet.jl" begin
    include("detrend.jl")
    include("surrogate.jl")
    include("cutofftime.jl")
    include("wtdelta.jl")
    include("coherence.jl")
    include("synthesis.jl")
    include("plots.jl")
end

ext: jl


text:
# vim: set ts=8 sw=4 sts=4 et:

using Unicode2LaTeX, Test

@testset "Unicode2LaTeX" begin

@testset "Conversion" begin

@test unicode2latex('a')   == "a"
@test unicode2latex('¥')   == "\\textyen"
@test unicode2latex("Foo") == "Foo"
@test unicode2latex("Föö") == "F\\\"o\\\"o"

end #testset

@testset "Generation" begin

@test Unicode2LaTeX.generate(nowarn=true) > 0

end #testset

@testset "LaTeX" begin

mktempdir() do d
    cd(d) do
        for s in values(Unicode2LaTeX.u2latex)
            # Bruteforce output encoding, can not be known in advance.
            t = false
            for e in ("T1", "X2", "T2B", "T2C", "OT2")
                p = open(`latex`, "w", devnull)
                println(p.in, """
                \\documentclass{article}

                \\usepackage[$e]{fontenc}
                % Additional text symbol support.
                \\usepackage{textcomp}

                \\begin{document}

                $s

                \\end{document}
                """)
                close(p.in)
                wait(p)
                t = p.exitcode == 0
                !t || break
            end
            @test t
        end
    end
end

end #testset

end #testset


text:
struct MPIData{A,B,C}
    comm::A
    world_sz::B
    rank::C
end

function MPIData(comm=default_comm())
    return MPIData(comm, MPI.Comm_size(comm), MPI.Comm_rank(comm))
end

function default_comm()
    # Initialize MPI if not already done
    MPI.REFCOUNT[] == -1 && MPI.Init()
    return MPI.COMM_WORLD
end

parallel_execution_cache(::ParallelMPI) = MPIData()
num_workers(tc::MPIData) = tc.world_sz

## stuff
function workers_mean!(target, data, par_type::MPIData)
    mean!(target, data)
    MPI.Allreduce!(MPI.IN_PLACE, target, MPI.SUM, par_type.comm)
    target ./= par_type.world_sz

    return target
end

function workers_mean!(target, par_type::MPIData)
    MPI.Allreduce!(MPI.IN_PLACE, target, MPI.SUM, par_type.comm)
    target ./= par_type.world_sz

    return target
end

function workers_mean(data, par_type::MPIData)
    data_m = mean(data)
    target = MPI.Allreduce(data_m, MPI.SUM, par_type.comm)

    return target/par_type.world_sz
end

function workers_sum!(target, par_type::MPIData)
    MPI.Allreduce!(target, MPI.SUM, par_type.comm)
    return target
end

function all_bcast_from_root!(data, par_type::MPIData)
    MPI.Bcast!(data, 0, par_type.comm)
    return data
end

function worker_local_seed(seed, par_cache::MPIData)
    rng = Random.MersenneTwister(seed)
    local_seeds = rand(rng, eltype(seed), par_cache.world_sz)
    MPI.Bcast!(local_seeds, 0, par_cache.comm)
    rank = MPI.Comm_rank(par_cache.comm) +1
    return local_seeds[rank]
end

function worker_block(iter, par_cache::MPIData)
    return _iterator_blocks(iter, MPI.Comm_rank(par_cache.comm) + 1, par_cache.world_sz)
end

my_block(stuff, par_cache::MPIData) = stuff[MPI.Comm_rank(par_cache.comm) + 1]

function gather_data!(data, intervals, par_cache::MPIData)
    loc_len = 3
end

function worker_allgatherv!(data, counts, par_cache::MPIData)
    MPI.Allgatherv!(MPI.IN_PLACE, data, counts, par_cache.comm)
    return data
end


text:
const ProgramParamMap = Dict(
   :delete_status                         => GL_DELETE_STATUS,
   :link_status                           => GL_LINK_STATUS,
   :validate_status                       => GL_VALIDATE_STATUS,
   :info_log_length                       => GL_INFO_LOG_LENGTH,
   :attached_shaders                      => GL_ATTACHED_SHADERS,
   :active_atomic_counter_buffers         => GL_ACTIVE_ATOMIC_COUNTER_BUFFERS,
   :active_attributes                     => GL_ACTIVE_ATTRIBUTES,
   :active_attribute_max_length           => GL_ACTIVE_ATTRIBUTE_MAX_LENGTH,
   :active_uniforms                       => GL_ACTIVE_UNIFORMS,
   :active_uniform_blocks                 => GL_ACTIVE_UNIFORM_BLOCKS,
   :active_uniform_block_max_name_length  => GL_ACTIVE_UNIFORM_BLOCK_MAX_NAME_LENGTH,
   :active_uniform_max_length             => GL_ACTIVE_UNIFORM_MAX_LENGTH,
   :program_binary_length                 => GL_PROGRAM_BINARY_LENGTH,
   :transform_feedback_buffer_mode        => GL_TRANSFORM_FEEDBACK_BUFFER_MODE,
   :transform_feedback_varyings           => GL_TRANSFORM_FEEDBACK_VARYINGS,
   :transform_feedback_varying_max_length => GL_TRANSFORM_FEEDBACK_VARYING_MAX_LENGTH,
   :geometry_vertices_out                 => GL_GEOMETRY_VERTICES_OUT,
   :geometry_input_type                   => GL_GEOMETRY_INPUT_TYPE,
   :geometry_output_type                  => GL_GEOMETRY_OUTPUT_TYPE,
)

const ShaderParamMap = Dict(
   :type                 => GL_SHADER_TYPE,
   :compile_status       => GL_COMPILE_STATUS,
   :delete_status        => GL_DELETE_STATUS,
   :info_log_length      => GL_INFO_LOG_LENGTH,
   :shader_source_length => GL_SHADER_SOURCE_LENGTH
)

const ShaderTypeMap = Dict(
   :vertex     => GL_VERTEX_SHADER,
   :vert       => GL_VERTEX_SHADER,
   :fragment   => GL_FRAGMENT_SHADER,
   :frag       => GL_FRAGMENT_SHADER,
   :geometry   => GL_GEOMETRY_SHADER,
   :geom       => GL_GEOMETRY_SHADER,
   :compute    => GL_COMPUTE_SHADER,
   :comp       => GL_COMPUTE_SHADER,
)

const BufferTypeMap = Dict(
   :array               => GL_ARRAY_BUFFER,
   :atomic_counter      => GL_ATOMIC_COUNTER_BUFFER,
   :copy_read           => GL_COPY_READ_BUFFER,
   :copy_write          => GL_COPY_WRITE_BUFFER,
   :dispatch_indirect   => GL_DISPATCH_INDIRECT_BUFFER,
   :draw_indirect       => GL_DRAW_INDIRECT_BUFFER,
   :element_array       => GL_ELEMENT_ARRAY_BUFFER,
   :pixel_pack          => GL_PIXEL_PACK_BUFFER,
   :pixel_unpack        => GL_PIXEL_UNPACK_BUFFER,
   :query               => GL_QUERY_BUFFER,
   :shader_storage      => GL_SHADER_STORAGE_BUFFER,
   :texture             => GL_TEXTURE_BUFFER,
   :transform_feedback  => GL_TRANSFORM_FEEDBACK_BUFFER,
   :uniform             => GL_UNIFORM_BUFFER,
)

const BufferUsageMap = Dict(
   :static_draw  => GL_STATIC_DRAW,
   :static_copy  => GL_STATIC_COPY,
   :static_read  => GL_STATIC_READ,
   :stream_draw  => GL_STREAM_DRAW,
   :stream_copy  => GL_STREAM_COPY,
   :stream_read  => GL_STREAM_READ,
   :dynamic_draw => GL_DYNAMIC_DRAW,
   :dynamic_copy => GL_DYNAMIC_COPY,
   :dynamic_read => GL_DYNAMIC_READ,
)

const DebugSource = Dict(
   GL_DEBUG_SOURCE_API => "API",
   GL_DEBUG_SOURCE_WINDOW_SYSTEM => "Window",
   GL_DEBUG_SOURCE_SHADER_COMPILER => "Shader Compiler",
   GL_DEBUG_SOURCE_THIRD_PARTY => "Third Party",
   GL_DEBUG_SOURCE_APPLICATION => "Application",
   GL_DEBUG_SOURCE_OTHER => "Unknown"
)

const DebugType = Dict(
   GL_DEBUG_TYPE_ERROR => "Error",
   GL_DEBUG_TYPE_DEPRECATED_BEHAVIOR => "Deprecated Behavior",
   GL_DEBUG_TYPE_UNDEFINED_BEHAVIOR => "Undefined Behavior",
   GL_DEBUG_TYPE_PORTABILITY => "Portability",
   GL_DEBUG_TYPE_PERFORMANCE => "Performance",
   GL_DEBUG_TYPE_OTHER => "Other",
   GL_DEBUG_TYPE_MARKER => "Marker",
)

const DebugSeverity = Dict(
   GL_DEBUG_SEVERITY_HIGH => "high",
   GL_DEBUG_SEVERITY_MEDIUM => "medium",
   GL_DEBUG_SEVERITY_LOW => "low",
   GL_DEBUG_SEVERITY_NOTIFICATION => "notification",
)

max_stars_repo_name: jorge-brito/Glutils

max_issues_repo_name: jorge-brito/Glutils

max_forks_repo_name: jorge-brito/Glutils


text:
module PerceptualHashes

using Images, Colors, FixedPointNumbers

if VERSION > v"0.6"
  using FFTW
end

import Base.convert

"""

  `hist_equalise(image, levels=256)`

  Equalise a greyscale `image` by applying a histogram of `levels` levels to it.
  Equalisation 'sharpens' an image by redistributing pixel values across the
  range of possible outputs. Returns the equalised `ImageMeta`, also greyscale.
"""
function hist_equalise(image::ImageMeta, levels=256)
  edges, counts = imhist(image, levels)
  cdf = cumsum(counts)
  cdf_min = cdf[1]
  denum = (width(image) * height(image)) - cdf_min

  h = Dict()
  for i in 1:length(edges)
    h[edges[i]] = (round( ((cdf[i] - cdf_min) / denum) * (levels-2) ) + 1)/levels
  end

  imnew = map( x -> h[edges[max(findlast(y -> y <= Float64(gray(x)), edges),edges[1])]], data(greyscale(image)))
  return ImageMeta(imnew, colorspace="Gray", spatialorder="x y")
end

function greyscale(image::ImageMeta)
  return convert(ImageMeta{Gray{Normed{UInt8, 8}}}, image)
end

function imapply(func, image)
  return ImageMeta(map(func, data(image)), colorspace=image.properties["colorspace"], spatialorder=image.properties["spatialorder"])
end

"""
 convert(::Type{UInt32}, hash::BitArray{1})

Converts a bitarray into an UInt. I feel like I shouldn't have to write this,
but `reinterpret` doesn't seem to work.
"""
function convert(::Type{UInt32}, hsh::BitArray{1})
  return uintof(0x00000000, hsh)
end

"""
 convert(::Type{UInt32}, hash::BitArray{1})

Converts a bitarray into an UInt. I feel like I shouldn't have to write this,
but `reinterpret` doesn't seem to work.
"""
function convert(::Type{UInt64}, hsh::BitArray{1})
  return uintof(0x0000000000000000, hsh)
end

function uintof(i, hsh)
  for hbit in hsh
    i = i << 1
    if hbit
      i = i + 0x1
    end
  end
  return i
end

"""
  `hdist(hashone, hashtwo)`

  Returns the Hamming distance between two hashes, expressed as `BitArrays`.
"""
function hdist(hashone, hashtwo)
  return sum(hashone $ hashtwo)
end

"""
  `bmv_hash(image, imgsize=256, blocksize=16)`

  Implements the Block Mean Value hashing system proposed by Yang et al (2006),
  discarding the index encryption step in line with Zauner (2010). Returns a
  hash represented as a `BitArray`, with a default length of 256 bits
  (32 bytes).
"""
function bmv_hash(image::ImageMeta, imgsize=256, blocksize=16)
  greyed = greyscale(image)
  prepped = Images.imresize(greyed, (imgsize, imgsize))
  raw = map(x -> gray(x), data(prepped))

  blockmeans = []
  #Move a blocksize-square window over the image, taking the mean.
  for i in 0:blocksize:(imgsize-blocksize), j in 0:blocksize:(imgsize-blocksize)
    block = raw[i+1:i+blocksize,j+1:j+blocksize]
    blockmean = mean(block)
    push!(blockmeans, blockmean)
  end

  #Identifiy the median block and code all blocks as >= or not.
  medblock = median(blockmeans)
  hash = BitArray(map(bm-> bm >= medblock, blockmeans))
  return hash
end

function dct_matrix(size)
  c = Matrix(size,size)
  coeff = sqrt(2/size)
  div = 2*size
  for n in 0:(size-1)
    for m in 0:(size-1)
      c[n+1, m+1] =  coeff * cos( (((2*m) + 1) * n * pi) / div )
    end
  end
  return c
end

"""
  `dct_hash(image)`

  Implements a perceptual hash based on the discrete cosine transform, as
  described by Zauner (2010) and implemented by the [pHash library](http://www.phash.org/docs/design.html)
  While in some sense functional, this implementation seems to be flawed in
  a manner not yet resolved by the author. Returns a 64-bit hash as a
  `BitArray`.

  You may want to use `convert(UInt64, BitArray)`.
"""
function dct_hash(image::ImageMeta, imgsize=32)
  greyed = greyscale(image)

  mfilter = centered(imaverage((7,7)))
  smoothed = imfilter(greyed, mfilter)

  prepped = Images.imresize(smoothed, (imgsize, imgsize))

  singular = map(x -> gray(x), data(prepped))
  dcted = dct(singular)

  #This is the implementation given by pHash, but it seems stupid to me when
  #compared to the option of taking the low-frequency components from each of
  #the 8x8 grid points usually used in compression.
  #It also produces different values to that of their library, for no reason
  #I can understand from comparing their source to mine.
  pixvals = []
  for i in 2:9
    for j in 2:9
      push!(pixvals, dcted[j,i])
    end
  end

  med = median(pixvals)
  hash = BitArray(map(x -> x >= med, pixvals))
  return hash
end

function mh_kern(alpha, level)
  sigma = 4 * (alpha^level)
  size = 2*sigma+1
  mul = 1/(alpha^level)
  kern = Matrix(size,size)
  for x in 1:size, y in 1:size
    xpos = mul*(x-sigma)
    ypos = mul*(y-sigma)
    a = xpos^2 + ypos^2
    kern[x,y] = (2-a)*exp(-a/2)
  end
  return kern
end

function mh_hash(image::ImageMeta, alpha=2, lvl=1)
  greyed = greyscale(image)
  blurred = imfilter_gaussian(convert(ImageMeta{RGB}, greyed), ones(ndims(greyed)))
  resized = Images.imresize(greyscale(blurred), (512, 512))
  equalised = hist_equalise(resized)

  mhkern = mh_kern(alpha, level)
#  imfilter_LoG(img, mhkern) <- not this, need to implement 'correlate'
  #weird hash bullshit
end

"""
  `resolve_images(dirname, tolerance=0, method=bmv)`

  Given a directory name, produces a dictionary mapping all contained images to
  those other images similar within a tolerance of `tolerance`.
"""
function resolve_images(dirname, tolerance=0, method=bmv_hash)
  fnames = readdir(dirname)
  if length(fnames) < 2
   return Dict()
  end
  fnames = map(x-> string(dirname,"/",x), fnames)
  hashes = map(x -> perceptual_hash(x, method), fnames)
  joint = zip(fnames, hashes)
  imgmap = Dict()
  for pair in combinations(collect(joint), 2)
    imgone = pair[1]
    imgtwo = pair[2]
    if ! (imgone[1] in keys(imgmap))
      imgmap[imgone[1]] = []
    end
    if ! (imgtwo[1] in keys(imgmap))
      imgmap[imgtwo[1]] = []
    end
    if hdist(imgone[2], imgtwo[2]) <= tolerance
      push!(imgmap[imgone[1]], imgtwo[1])
      push!(imgmap[imgtwo[1]], imgone[1])
    end
  end
  return imgmap
end

"""
  `perceptual_hash(filename, method=bmv_hash)`

  Create a perceptual hash of an image file found at `filename` using `method`.

  See also: `bmv`
"""
function perceptual_hash(filename, method=bmv_hash)
  imfile = ImageMeta(load(filename))
  return method(imfile)
end

export perceptual_hash, bmv_hash, dct_hash, hdist, resolve_images

end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__,".."))
using DCAUtils
using Statistics

using Distances

function get_seq_len(fasta)
    # ATTENTION: ASSUMES NO LINE BREAKS IN SEQUENCES

    seq_len = 0
    open(fasta) do fid
        for (ind, line) in enumerate(eachline(fid))
            if ind == 1
                @assert line[1]  == '>'
                continue
            end
            if ind == 2
                seq_len = length(strip(line))
                continue
            end
            if ind == 3
                @assert line[1] == '>'
                break
            end
        end
    end

    return seq_len
end

function create_hamming_distances(data_dir)

    dist = Hamming()

    train_fastas = filter(s->endswith(s, "a2m.train"), readdir(data_dir))

    seq_len = readlines

    for train_fasta in train_fastas

        println("Doing $train_fasta")

        test_fasta = replace(train_fasta, "a2m.train" => "a2m.test")

        train_fasta_path = joinpath(data_dir, train_fasta)
        test_fasta_path = joinpath(data_dir, test_fasta)

        if !isfile(train_fasta_path)
            error("cannot find $train_fasta_path")
        end

        if !isfile(test_fasta_path)
            error("cannot find $test_fasta_path")
        end

        # get number of test sequences for consisteny checks later
        ntest = parse(Int, read(`grep -c ">" $test_fasta_path`, String))

        seq_len = get_seq_len(train_fasta_path)

        train_seqs = read_fasta_alignment(train_fasta_path, 1.0)
        test_seqs = read_fasta_alignment(test_fasta_path, 1.0)

        hmd = pairwise(dist, train_seqs, test_seqs, dims=2)

        mean_hmds = vec(mean(hmd, dims=1))
        min_hmds = vec(minimum(hmd, dims=1))

        mean_min_hmds = mean(min_hmds)
        println("mean minimum hamming distance = $mean_min_hmds")

        @assert length(mean_hmds) == ntest

        # assert that there are no zeros
        @assert minimum(min_hmds) > 0

        out_path = replace(train_fasta_path, ".train"=>".train_test_hamming")

        @assert out_path != train_fasta_path

        open(out_path, "w") do fid
            for i in 1:length(mean_hmds)
                s = string.([mean_hmds[i], mean_hmds[i]/seq_len, min_hmds[i], min_hmds[i]/seq_len])
                print(fid, join(s, " "))
                print(fid, '\n')
              end
        end

    end

end

data_dir=ARGS[1]

if !isdir(data_dir)
    error("$data_dir does not exist")
else
    create_hamming_distances(data_dir)
end

max_stars_repo_name: christophfeinauer/PairwiseDistillations

max_issues_repo_name: christophfeinauer/PairwiseDistillations

max_forks_repo_name: christophfeinauer/PairwiseDistillations


text:
abstract type AbstractSyntheticSatellite{T} <: AbstractEmitter{T} end

struct SyntheticSatellite{
    T <: AbstractFloat,
    CS <: ConstantDopplerSatellite{<: AbstractGNSS, T},
    S <: Union{Nothing, <:AbstractVector}
} <: AbstractSyntheticSatellite{T}
    sat::CS
    steer_vec::S
end

function SyntheticSatellite(sat::ConstantDopplerSatellite{<:AbstractGNSS, T}) where T <: AbstractFloat
    SyntheticSatellite(sat, nothing)
end

function gen_signal!(
    synthetic_sat::SyntheticSatellite,
    sample_frequency,
    intermediate_frequency,
    n0,
    num_samples::Integer,
    rng
)
    gen_signal!(
        synthetic_sat.sat,
        sample_frequency,
        intermediate_frequency,
        n0,
        num_samples,
        rng
    )
end

function propagate(
    synthetic_sat::SyntheticSatellite,
    num_samples,
    intermediate_frequency,
    sample_frequency,
    n0,
    rng
)
    SyntheticSatellite(
        propagate(
            synthetic_sat.sat,
            num_samples,
            intermediate_frequency,
            sample_frequency,
            n0,
            rng
        ),
        synthetic_sat.steer_vec
    )
end

function get_steer_vec(
    manifold::AbstractManifold{N},
    emitter::SyntheticSatellite{T, CS, Nothing},
    attitude
) where {N, T <: AbstractFloat, CS <: ConstantDopplerSatellite}
    ones(SVector{N, T})
end

function get_steer_vec(
    manifold::IdealManifold{1},
    emitter::SyntheticSatellite{T, CS, Nothing},
    attitude
) where {T <: AbstractFloat, CS <: ConstantDopplerSatellite}
    one(T)
end

function get_steer_vec(
    manifold::AbstractManifold,
    emitter::SyntheticSatellite,
    attitude
)
    emitter.steer_vec
end

@inline get_doa(ss::SyntheticSatellite) = get_doa(ss.sat)
@inline get_existence(ss::SyntheticSatellite) = get_existence(ss.sat)
@inline get_carrier_doppler(ss::SyntheticSatellite) = get_carrier_doppler(ss.sat)
@inline get_code_doppler(ss::SyntheticSatellite) = get_code_doppler(ss.sat)
@inline get_carrier_phase(ss::SyntheticSatellite) = get_carrier_phase(ss.sat)
@inline get_code_phase(ss::SyntheticSatellite) = get_code_phase(ss.sat)
@inline get_prn(ss::SyntheticSatellite) = get_prn(ss.sat)
@inline get_amplitude(ss::SyntheticSatellite, n0) = get_amplitude(ss.sat, n0)
@inline get_gnss_system(ss::SyntheticSatellite) = get_gnss_system(ss.sat)
@inline get_carrier_to_noise_density_ratio(ss::SyntheticSatellite) =
    get_carrier_to_noise_density_ratio(ss.sat)


text:
using Multilane
using MCTS
using RobustMCTS
using POMCP
using GenerativeModels
using POMDPToolbox
using JLD
using POMDPs
using DataFrames
using Plots

N=500

nb_lanes = 4
desired_lane_reward = 10.
rmodel = NoCrashRewardModel(desired_lane_reward*10., desired_lane_reward,2.5,nb_lanes)

pp = PhysicalParam(nb_lanes,lane_length=100.)

normal_behavior = IDMMOBILBehavior("normal", pp.v_med, pp.l_car, 1)
behaviors=IDMMOBILBehavior[
                    normal_behavior,
                    IDMMOBILBehavior("cautious", pp.v_slow+0.5, pp.l_car, 2),
                    IDMMOBILBehavior("aggressive", pp.v_fast, pp.l_car, 3)
                    ]

_discount = 1.0
nb_cars = 10
dmodel = NoCrashIDMMOBILModel(nb_cars, pp,
                              behaviors=behaviors,
                              behavior_probabilities=WeightVec([0.2, 0.4, 0.4]),
                              lane_terminate=false)

base_problem = NoCrashMDP(dmodel, rmodel, _discount)

problems = Dict{String, Any}()
key_rng = MersenneTwister(123)

is_rng = MersenneTwister(123)
init_phys_states = [MLPhysicalState(initial_state(base_problem, is_rng)) for i in 1:N]
initial_states = Dict{String, Any}()
state_lists = Dict{String, Any}()

params = Dict{Symbol, Any}()
params[:lambda] = Float64[0.1, 1., 10., 17.78, 31.62, 56.23, 100., 1000.]
# params[:lambda] = Float64[0.1, 1.]
params[:p_normal] = Float64[0.2, 1.0]

initial_relevant = [:p_normal]

param_table = DataFrame()
for p in params
    if isempty(param_table)
        param_table = DataFrame(Dict(p))
    else
        param_table = join(param_table, DataFrame(Dict(p)), kind = :cross)
    end
end

param_table[:problem_key] = DataArray(String, nrow(param_table))
param_table[:state_list_key] = DataArray(String, nrow(param_table))

for row in eachrow(param_table)
    key = randstring(key_rng)
    problem = deepcopy(base_problem)
    # lambda
    problem.rmodel.cost_dangerous_brake = row[:lambda]*problem.rmodel.reward_in_desired_lane
    # p_normal
    nb_behaviors = length(problem.dmodel.behaviors)
    probs = ones(nb_behaviors) * (1.0 - row[:p_normal])/(nb_behaviors-1)
    probs[1] = row[:p_normal]
    problem.dmodel.behavior_probabilities = WeightVec(probs)

    problems[key] = problem
    row[:problem_key] = key
end

for g in groupby(param_table, initial_relevant)
    key = randstring(key_rng)
    pk = g[1,:problem_key]
    p = problems[pk]
    these_states = assign_keys([initial_state(p, init_phys_states[i], is_rng) for i in 1:N])
    merge!(initial_states, these_states)
    g[:,:state_list_key] = key
    state_lists[key] = collect(keys(these_states))
end

objects = Dict{String, Any}()
objects["problems"] = problems
objects["param_table"] = param_table
objects["initial_states"] = initial_states
objects["state_lists"] = state_lists

# Just Figure 1 with 5 simulations

dpws = DPWSolver(depth=20,
                 n_iterations=300,
                 exploration_constant=50.0,
                 k_state=4.0,
                 alpha_state=1/8,
                 rollout_solver=SimpleSolver())

rsolver = RobustMCTSSolver(
    depth=20,
    c=50.0,
    n_iterations=1500,
    k_state=4.0,
    alpha_state=1/8,
    k_nature=4.0,
    alpha_nature=1/8,
    rollout_solver=SimpleSolver(),
    rollout_nature=StochasticBehaviorNoCrashMDP(collect(values(problems))[1]))

pomcps = POMCPDPWSolver(
    eps=0.001,
    max_depth=20,
    c=50.0,
    tree_queries=1500,
    k_observation=4.0,
    alpha_observation=1/8,
    rollout_solver=SimpleSolver())

# for 20% and 100 percent, we need the upper bound, the heuristic and random
point_solvers = Dict{String, Solver}(
    "random" => RandomSolver(),
    "heuristic" => SimpleSolver()
)

objects["solvers"] = deepcopy(point_solvers)
objects["solvers"]["upper_bound"] = dpws
objects["solvers"]["pomcp"] = MLPOMDPSolver(pomcps, BehaviorRootUpdaterStub(0.1))
objects["solvers"]["robust"] = RobustMLSolver(rsolver)

results = objects

problems_100 = param_table[param_table[:p_normal].==1.0, :problem_key]
problems_20 = param_table[param_table[:p_normal].==0.2, :problem_key]
solver_pairs = [("upper_bound", problems_100),
                ("pomcp", problems_20),
                ("robust", problems_20)]

for p_normal in unique(param_table[:p_normal])
    list = param_table[param_table[:p_normal].==p_normal, :state_list_key][1]
    inits = state_lists[list]
    problems = param_table[param_table[:p_normal].==p_normal, :problem_key]

    r = evaluate([problems[1]], inits, collect(keys(point_solvers)), objects, N=N, desc="$p_normal Point: ")
    merge_results!(results, r)

    r = evaluate(problems, inits, ["upper_bound"], objects, N=N, desc="$p_normal Upper: ")
    merge_results!(results, r)

    for sp in solver_pairs
        r = evaluate(problems, inits, [sp[1]], objects, soln_problem_keys=sp[2], N=N, desc="$p_normal $(sp[1]): ")
        merge_results!(results, r)
    end
end

#=
# Calculate single points
point_solvers = Dict{String, Solver}(
    "random" => RandomSolver(),
    "heuristic" => SimpleSolver()
)

point_results = evaluate(representative_problem, initial_states, point_solvers, parallel=true, N=N)

dpws = DPWSolver(depth=20,
                 n_iterations=300,
                 exploration_constant=50.0,
                 k_state=4.0,
                 alpha_state=1/8,
                 rollout_solver=SimpleSolver())

rsolver = RobustMCTSSolver(
    depth=20,
    c=50.0,
    n_iterations=1500,
    k_state=4.0,
    alpha_state=1/8,
    k_nature=4.0,
    alpha_nature=1/8,
    rollout_solver=SimpleSolver(),
    rollout_nature=StochasticBehaviorNoCrashMDP(collect(values(problems))[1]))

pomcps = POMCPDPWSolver(
    eps=0.001,
    max_depth=20,
    c=50.0,
    tree_queries=1500,
    k_observation=4.0,
    alpha_observation=1/8,
    rollout_solver=SimpleSolver()
)

curve_solvers = Dict{String, Solver}(
    "upper_bound" => dpws,
    "robust" => RobustMLSolver(rsolver),
    "single_behavior" =>
        SingleBehaviorSolver(dpws,
                 IDMMOBILBehavior("normal", 30.0, 10.0, 1)),
    "pomcp" => MLPOMDPSolver(pomcps, BehaviorRootUpdaterStub(0.1))
)

curve_results = evaluate(problems, initial_states, curve_solvers, parallel=true, N=N)

results = merge_results!(curve_results, point_results)
=#

println(results["stats"])

filename = string("results_", Dates.format(Dates.now(),"u_d_HH_MM"), ".jld")
results["histories"] = nothing
save(filename, results)
println("results saved to $filename")

#=
stats = results["stats"]
mean_performance = by(stats, :solver_key) do df
    by(df, :lambda) do df
        DataFrame(steps_in_lane=mean(df[:steps_in_lane]),
                  nb_brakes=mean(df[:nb_brakes]),
                  time_to_lane=mean(df[:time_to_lane]),
                  steps=mean(df[:steps]),
                  )
    end
end

unicodeplots()
for g in groupby(mean_performance, :solver_key)
    if size(g,1) > 1
        # plot!(g, :steps_in_lane, :nb_brakes, group=:solver_key)
        plot!(g, :time_to_lane, :nb_brakes, group=:solver_key)
    else
        # scatter!(g, :steps_in_lane, :nb_brakes, group=:solver_key)
        scatter!(g, :time_to_lane, :nb_brakes, group=:solver_key)
    end
end
gui()
=#


text:
using ArgParseLite
using Test

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


text:
export update_halo!

import MPI
@static if ENABLE_CUDA
    using CUDA
end
using Base.Threads

"""
    update_halo!(A)
    update_halo!(A...)

Update the halo of the given GPU/CPU-array(s).

# Typical use cases:
    update_halo!(A)        # Update the halo of the array A.
    update_halo!(A, B, C)  # Update the halos of the arrays A, B and C.

!!! note "Performance note"
    Group subsequent calls to `update_halo!` in a single call for better performance (enables additional pipelining).
    Consider activating CUDA-aware MPI (see [`ImplicitGlobalGrid`](@ref)).
"""
function update_halo!(A::GGArray...)
    check_initialized();
    check_fields(A...);
    _update_halo!(A...);  # Asignment of A to fields in the internal function _update_halo!() as vararg A can consist of multiple fields; A will be used for a single field in the following (The args of update_halo! must however be "A..." for maximal simplicity and elegance for the user).
    return nothing
end

function _update_halo!(fields::GGArray...)
    allocate_bufs(fields...);
    if any_array(fields...) allocate_tasks(fields...); end
    @enable_if_cuda if any_cuarray(fields...) allocate_custreams(fields...); end

    for dim = 1:NDIMS_MPI  # NOTE: this works for 1D-3D (e.g. if nx>1, ny>1 and nz=1, then for d=3, there will be no neighbors, i.e. nothing will be done as desired...).
        for ns = 1:NNEIGHBORS_PER_DIM,  i = 1:length(fields)
            if has_neighbor(ns, dim) iwrite_sendbufs!(ns, dim, fields[i], i); end
        end
        # Send / receive if the neighbors are other processes (usual case).
        reqs = fill(MPI.REQUEST_NULL, length(fields), NNEIGHBORS_PER_DIM, 2);
        if all(neighbors(dim) .!= me())  # Note: handling of send/recv to itself requires special configurations for some MPI implementations (e.g. self BTL must be activated with OpenMPI); so we handle this case without MPI to avoid this complication.
            for nr = NNEIGHBORS_PER_DIM:-1:1,  i = 1:length(fields) # Note: if there were indeed more than 2 neighbors per dimension; then one would need to make sure which neigbour would communicate with which.
                if has_neighbor(nr, dim) reqs[i,nr,1] = irecv_halo!(nr, dim, fields[i], i); end
            end
            for ns = 1:NNEIGHBORS_PER_DIM,  i = 1:length(fields)
                if has_neighbor(ns, dim)
                    wait_iwrite(ns, fields[i], i);  # Right before starting to send, make sure that the data of neighbor ns and field i has finished writing to the sendbuffer.
                    reqs[i,ns,2] = isend_halo(ns, dim, fields[i], i);
                end
            end
        # Copy if I am my own neighbors (when periodic boundary and only one process in this dimension).
        elseif all(neighbors(dim) .== me())
            for ns = 1:NNEIGHBORS_PER_DIM,  i = 1:length(fields)
                wait_iwrite(ns, fields[i], i);  # Right before starting to send, make sure that the data of neighbor ns and field i has finished writing to the sendbuffer.
                sendrecv_halo_local(ns, dim, fields[i], i);
                nr = NNEIGHBORS_PER_DIM - ns + 1;
                iread_recvbufs!(nr, dim, fields[i], i);
            end
        else
            error("Incoherent neighbors in dimension $dim: either all neighbors must equal to me, or none.")
        end
        for nr = NNEIGHBORS_PER_DIM:-1:1,  i = 1:length(fields)  # Note: if there were indeed more than 2 neighbors per dimension; then one would need to make sure which neigbour would communicate with which.
            if (reqs[i,nr,1]!=MPI.REQUEST_NULL) MPI.Wait!(reqs[i,nr,1]); end
            if (has_neighbor(nr, dim) && neighbor(nr, dim)!=me()) iread_recvbufs!(nr, dim, fields[i], i); end  # Note: if neighbor(nr,dim) != me() is done directly in the sendrecv_halo_local loop above for better performance (thanks to pipelining)
        end
        for nr = NNEIGHBORS_PER_DIM:-1:1,  i = 1:length(fields) # Note: if there were indeed more than 2 neighbors per dimension; then one would need to make sure which neigbour would communicate with which.
            if has_neighbor(nr, dim) wait_iread(nr, fields[i], i); end
        end
        for ns = 1:NNEIGHBORS_PER_DIM
            if (any(reqs[:,ns,2].!=[MPI.REQUEST_NULL])) MPI.Waitall!(reqs[:,ns,2]); end
        end
    end
end

##---------------------------
## FUNCTIONS FOR SYNTAX SUGAR

halosize(dim::Integer, A::GGArray) = (ndims(A)>1) ? size(A)[1:ndims(A).!=dim] : (1,);

##---------------------------------------
## FUNCTIONS RELATED TO BUFFER ALLOCATION

let
    global free_update_halo_buffers, allocate_bufs, sendbuf, recvbuf, sendbuf_flat, recvbuf_flat, cusendbuf, curecvbuf, cusendbuf_flat, curecvbuf_flat
    sendbufs_raw = nothing
    recvbufs_raw = nothing
    cusendbufs_raw = nothing
    curecvbufs_raw = nothing
    cusendbufs_raw_h = nothing
    curecvbufs_raw_h = nothing

    function free_update_halo_buffers()
        @enable_if_cuda if any(cudaaware_MPI()) free_cubufs(cusendbufs_raw) end
        @enable_if_cuda if any(cudaaware_MPI()) free_cubufs(curecvbufs_raw) end
        @enable_if_cuda if none(cudaaware_MPI()) unregister_cubufs(cusendbufs_raw_h) end
        @enable_if_cuda if none(cudaaware_MPI()) unregister_cubufs(curecvbufs_raw_h) end
        sendbufs_raw = nothing;
        recvbufs_raw = nothing;
        cusendbufs_raw = nothing
        curecvbufs_raw = nothing
        cusendbufs_raw_h = nothing
        curecvbufs_raw_h = nothing
        GC.gc();
    end

    @static if ENABLE_CUDA
        function free_cubufs(bufs)
            if (bufs !== nothing)
                for i = 1:length(bufs)
                    for n = 1:length(bufs[i])
                        if is_cuarray(bufs[i][n]) CUDA.unsafe_free!(bufs[i][n]); bufs[i][n] = []; end
                    end
                end
            end
        end
        function unregister_cubufs(bufs)
            if (bufs !== nothing)
                for i = 1:length(bufs)
                    for n = 1:length(bufs[i])
                        if (isa(bufs[i][n],CUDA.Mem.HostBuffer)) CUDA.Mem.unregister(bufs[i][n]); bufs[i][n] = []; end
                    end
                end
            end
        end
    end

    # Allocate for each field two send and recv buffers (one for the left and one for the right neighbour of a dimension). The required length of the buffer is given by the maximal number of halo elements in any of the dimensions. Note that buffers are not allocated separately for each dimension, as the updates are performed one dimension at a time (required for correctness).
    function allocate_bufs(fields::GGArray{T}...) where T <: GGNumber
        if (isnothing(sendbufs_raw) || isnothing(recvbufs_raw))
            free_update_halo_buffers();
            init_bufs_arrays();
            @enable_if_cuda init_cubufs_arrays();
        end
        init_bufs(T, fields...);
        @enable_if_cuda init_cubufs(T, fields...);
        for i = 1:length(fields)
            A = fields[i];
            for n = 1:NNEIGHBORS_PER_DIM # Ensure that the buffers are interpreted to contain elements of the same type as the array.
                reinterpret_bufs(T, i, n);
                @enable_if_cuda reinterpret_cubufs(T, i, n);
            end
            max_halo_elems = (ndims(A) > 1) ? prod(sort([size(A)...])[2:end]) : 1;
            if (length(sendbufs_raw[i][1]) < max_halo_elems)
                for n = 1:NNEIGHBORS_PER_DIM
                    reallocate_bufs(T, i, n, max_halo_elems);
                    @enable_if_cuda if (is_cuarray(A) && none(cudaaware_MPI())) reregister_cubufs(T, i, n); end  # Host memory is page-locked (and mapped to device memory) to ensure optimal access performance (from kernel or with 3-D memcopy).
                end
                GC.gc(); # Too small buffers had been replaced with larger ones; free the now unused memory.
            end
            if (!isnothing(cusendbufs_raw) && length(cusendbufs_raw[i][1]) < max_halo_elems)
                for n = 1:NNEIGHBORS_PER_DIM
                    @enable_if_cuda if (is_cuarray(A) &&  any(cudaaware_MPI())) reallocate_cubufs(T, i, n, max_halo_elems); GC.gc(); end # Too small buffers had been replaced with larger ones; free the unused memory immediately.
                end
            end
        end
    end

    function init_bufs_arrays()
        sendbufs_raw = Array{Array{Any,1},1}();
        recvbufs_raw = Array{Array{Any,1},1}();
    end

    function init_bufs(T::DataType, fields::GGArray...)
        while (length(sendbufs_raw) < length(fields)) push!(sendbufs_raw, [zeros(T,0), zeros(T,0)]); end
        while (length(recvbufs_raw) < length(fields)) push!(recvbufs_raw, [zeros(T,0), zeros(T,0)]); end
    end

    function reinterpret_bufs(T::DataType, i::Integer, n::Integer)
        if (eltype(sendbufs_raw[i][n]) != T) sendbufs_raw[i][n] = reinterpret(T, sendbufs_raw[i][n]); end
        if (eltype(recvbufs_raw[i][n]) != T) recvbufs_raw[i][n] = reinterpret(T, recvbufs_raw[i][n]); end
    end

    function reallocate_bufs(T::DataType, i::Integer, n::Integer, max_halo_elems::Integer)
        sendbufs_raw[i][n] = zeros(T, Int(ceil(max_halo_elems/GG_ALLOC_GRANULARITY))*GG_ALLOC_GRANULARITY); # Ensure that the amount of allocated memory is a multiple of 4*sizeof(T) (sizeof(Float64)/sizeof(Float16) = 4). So, we can always correctly reinterpret the raw buffers even if next time sizeof(T) is greater.
        recvbufs_raw[i][n] = zeros(T, Int(ceil(max_halo_elems/GG_ALLOC_GRANULARITY))*GG_ALLOC_GRANULARITY);
    end

    @static if ENABLE_CUDA
        function init_cubufs_arrays()
            cusendbufs_raw = Array{Array{Any,1},1}();
            curecvbufs_raw = Array{Array{Any,1},1}();
            cusendbufs_raw_h = Array{Array{Any,1},1}();
            curecvbufs_raw_h = Array{Array{Any,1},1}();
        end

        function init_cubufs(T::DataType, fields::GGArray...)
            while (length(cusendbufs_raw) < length(fields)) push!(cusendbufs_raw, [CuArray{T}(undef,0), CuArray{T}(undef,0)]); end
            while (length(curecvbufs_raw) < length(fields)) push!(curecvbufs_raw, [CuArray{T}(undef,0), CuArray{T}(undef,0)]); end
            while (length(cusendbufs_raw_h) < length(fields)) push!(cusendbufs_raw_h, [[], []]); end
            while (length(curecvbufs_raw_h) < length(fields)) push!(curecvbufs_raw_h, [[], []]); end
        end

        function reinterpret_cubufs(T::DataType, i::Integer, n::Integer)
            if (eltype(cusendbufs_raw[i][n]) != T) cusendbufs_raw[i][n] = reinterpret(T, cusendbufs_raw[i][n]); end
            if (eltype(curecvbufs_raw[i][n]) != T) curecvbufs_raw[i][n] = reinterpret(T, curecvbufs_raw[i][n]); end
        end

        function reallocate_cubufs(T::DataType, i::Integer, n::Integer, max_halo_elems::Integer)
            cusendbufs_raw[i][n] = cuzeros(T, Int(ceil(max_halo_elems/GG_ALLOC_GRANULARITY))*GG_ALLOC_GRANULARITY); # Ensure that the amount of allocated memory is a multiple of 4*sizeof(T) (sizeof(Float64)/sizeof(Float16) = 4). So, we can always correctly reinterpret the raw buffers even if next time sizeof(T) is greater.
            curecvbufs_raw[i][n] = cuzeros(T, Int(ceil(max_halo_elems/GG_ALLOC_GRANULARITY))*GG_ALLOC_GRANULARITY);
        end

        function reregister_cubufs(T::DataType, i::Integer, n::Integer)
            if (isa(cusendbufs_raw_h[i][n],CUDA.Mem.HostBuffer)) CUDA.Mem.unregister(cusendbufs_raw_h[i][n]); cusendbufs_raw_h[i][n] = []; end # It is always initialized registered... if (cusendbufs_raw_h[i][n].bytesize > 32*sizeof(T))
            if (isa(curecvbufs_raw_h[i][n],CUDA.Mem.HostBuffer)) CUDA.Mem.unregister(curecvbufs_raw_h[i][n]); cusendbufs_raw_h[i][n] = []; end # It is always initialized registered... if (curecvbufs_raw_h[i][n].bytesize > 32*sizeof(T))
            cusendbufs_raw[i][n], cusendbufs_raw_h[i][n] = register(sendbufs_raw[i][n]);
            curecvbufs_raw[i][n], curecvbufs_raw_h[i][n] = register(recvbufs_raw[i][n]);
        end
    end

    function sendbuf_flat(n::Integer, dim::Integer, i::Integer, A::GGArray{T}) where T <: GGNumber
        return view(sendbufs_raw[i][n]::AbstractVector{T},1:prod(halosize(dim,A)));
    end

    function recvbuf_flat(n::Integer, dim::Integer, i::Integer, A::GGArray{T}) where T <: GGNumber
        return view(recvbufs_raw[i][n]::AbstractVector{T},1:prod(halosize(dim,A)));
    end

    function sendbuf(n::Integer, dim::Integer, i::Integer, A::GGArray)
        return reshape(sendbuf_flat(n,dim,i,A), halosize(dim,A));
    end

    function recvbuf(n::Integer, dim::Integer, i::Integer, A::GGArray)
        return reshape(recvbuf_flat(n,dim,i,A), halosize(dim,A));
    end

    @static if ENABLE_CUDA
        function cusendbuf_flat(n::Integer, dim::Integer, i::Integer, A::CuArray{T}) where T <: GGNumber
            return view(cusendbufs_raw[i][n]::CuVector{T},1:prod(halosize(dim,A)));
        end

        function curecvbuf_flat(n::Integer, dim::Integer, i::Integer, A::CuArray{T}) where T <: GGNumber
            return view(curecvbufs_raw[i][n]::CuVector{T},1:prod(halosize(dim,A)));
        end

        #TODO: see if I should remove T here and in other cases for CuArray or Array (but then it does not verify that CuArray is of type GGNumber) or if I should instead change GGArray to GGArrayUnion and create: GGArray = Array{T} where T <: GGNumber  and  GGCuArray = CuArray{T} where T <: GGNumber; This is however more difficult to read and understand for others.
        function cusendbuf(n::Integer, dim::Integer, i::Integer, A::CuArray{T}) where T <: GGNumber
            return reshape(cusendbuf_flat(n,dim,i,A), halosize(dim,A));
        end

        function curecvbuf(n::Integer, dim::Integer, i::Integer, A::CuArray{T}) where T <: GGNumber
            return reshape(curecvbuf_flat(n,dim,i,A), halosize(dim,A));
        end
    end

    # Make sendbufs_raw and recvbufs_raw accessible for unit testing.
    global get_sendbufs_raw, get_recvbufs_raw, get_cusendbufs_raw, get_curecvbufs_raw
    get_sendbufs_raw()   = deepcopy(sendbufs_raw)
    get_recvbufs_raw()   = deepcopy(recvbufs_raw)
    get_cusendbufs_raw() = deepcopy(cusendbufs_raw)
    get_curecvbufs_raw() = deepcopy(curecvbufs_raw)
end

##----------------------------------------------
## FUNCTIONS TO WRITE AND READ SEND/RECV BUFFERS

# NOTE: the tasks and custreams are stored here in a let clause to have them survive the end of a call to update_boundaries. This avoids the creation of new tasks and cuda streams every time. Besides, that this could be relevant for performance, it is important for debugging the overlapping the communication with computation (if at every call new stream/task objects are created this becomes very messy and hard to analyse).

function allocate_tasks(fields::GGArray...)
    allocate_tasks_iwrite(fields...);
    allocate_tasks_iread(fields...);
end

let
    global iwrite_sendbufs!, allocate_tasks_iwrite, wait_iwrite

    tasks = Array{Task}(undef, NNEIGHBORS_PER_DIM, 0);

    wait_iwrite(n::Integer, A::Array{T}, i::Integer) where T <: GGNumber = (schedule(tasks[n,i]); wait(tasks[n,i]);) # The argument A is used for multiple dispatch. #NOTE: The current implementation only starts a task when it is waited for, in order to make sure that only one task is run at a time and that they are run in the desired order (best for performance as the tasks are mapped only to one thread via context switching).

    function allocate_tasks_iwrite(fields::GGArray...)
        if length(fields) > size(tasks,2)  # Note: for simplicity, we create a tasks for every field even if it is not an Array
            tasks = [tasks Array{Task}(undef, NNEIGHBORS_PER_DIM, length(fields)-size(tasks,2))];  # Create (additional) emtpy tasks.
        end
    end

    function iwrite_sendbufs!(n::Integer, dim::Integer, A::Array{T}, i::Integer) where T <: GGNumber  # Function to be called if A is a CPU array.
        tasks[n,i] = @task begin
            if ol(dim,A) >= 2  # There is only a halo and thus a halo update if the overlap is at least 2...
                write_h2h!(sendbuf(n,dim,i,A), A, sendranges(n,dim,A), dim);
            end
        end
    end

    # Make tasks accessible for unit testing.
    global get_tasks_iwrite
    get_tasks_iwrite() = deepcopy(tasks)
end

let
    global iread_recvbufs!, allocate_tasks_iread, wait_iread

    tasks = Array{Task}(undef, NNEIGHBORS_PER_DIM, 0);

    wait_iread(n::Integer, A::Array{T}, i::Integer) where T <: GGNumber = (schedule(tasks[n,i]); wait(tasks[n,i]);) # The argument A is used for multiple dispatch. #NOTE: The current implementation only starts a task when it is waited for, in order to make sure that only one task is run at a time and that they are run in the desired order (best for performance currently as the tasks are mapped only to one thread via context switching).

    function allocate_tasks_iread(fields::GGArray...)
        if length(fields) > size(tasks,2)  # Note: for simplicity, we create a tasks for every field even if it is not an Array
            tasks = [tasks Array{Task}(undef, NNEIGHBORS_PER_DIM, length(fields)-size(tasks,2))];  # Create (additional) emtpy tasks.
        end
    end

    function iread_recvbufs!(n::Integer, dim::Integer, A::Array{T}, i::Integer) where T <: GGNumber  # Function to be called if A is a CPU array.
        tasks[n,i] = @task begin
            if ol(dim,A) >= 2  # There is only a halo and thus a halo update if the overlap is at least 2...
                read_h2h!(recvbuf(n,dim,i,A), A, recvranges(n,dim,A), dim);
            end
        end
    end

    # Make tasks accessible for unit testing.
    global get_tasks_iread
    get_tasks_iread() = deepcopy(tasks)
end

@static if ENABLE_CUDA
    function allocate_custreams(fields::GGArray...)
        allocate_custreams_iwrite(fields...);
        allocate_custreams_iread(fields...);
    end

    let
        global iwrite_sendbufs!, allocate_custreams_iwrite, wait_iwrite

        custreams = Array{CuStream}(undef, NNEIGHBORS_PER_DIM, 0)

        wait_iwrite(n::Integer, A::CuArray{T}, i::Integer) where T <: GGNumber = synchronize(custreams[n,i]); # The argument A is used for multiple dispatch.

        function allocate_custreams_iwrite(fields::GGArray...)
    	    if length(fields) > size(custreams,2)  # Note: for simplicity, we create a stream for every field even if it is not a CuArray
                custreams = [custreams [CuStream(; flags=CUDA.STREAM_NON_BLOCKING, priority=CUDA.priority_range()[end]) for n=1:NNEIGHBORS_PER_DIM, i=1:(length(fields)-size(custreams,2))]];  # Create (additional) maximum priority nonblocking streams to enable overlap with computation kernels.
            end
        end

        function iwrite_sendbufs!(n::Integer, dim::Integer, A::CuArray{T}, i::Integer) where T <: GGNumber  # Function to be called if A is a GPU array.
            if ol(dim,A) >= 2  # There is only a halo and thus a halo update if the overlap is at least 2...
                if dim == 1 || cudaaware_MPI(dim) # Use a custom copy kernel for the first dimension to obtain a good copy performance (the CUDA 3-D memcopy does not perform well for this extremely strided case).
                    ranges = sendranges(n, dim, A);
                    nthreads = (dim==1) ? (1, 32, 1) : (32, 1, 1);
                    halosize = [r[end] - r[1] + 1 for r in ranges];
                    nblocks  = Tuple(ceil.(Int, halosize./nthreads));
                    @cuda blocks=nblocks threads=nthreads stream=custreams[n,i] write_d2x!(cusendbuf(n,dim,i,A), A, ranges[1], ranges[2], ranges[3], dim);
                else
                    write_d2h_async!(sendbuf_flat(n,dim,i,A), A, sendranges(n,dim,A), dim, custreams[n,i]);
                end
            end
        end
    end

    let
        global iread_recvbufs!, allocate_custreams_iread, wait_iread

        custreams = Array{CuStream}(undef, NNEIGHBORS_PER_DIM, 0)

        wait_iread(n::Integer, A::CuArray{T}, i::Integer) where T <: GGNumber = synchronize(custreams[n,i]); # The argument A is used for multiple dispatch.

        function allocate_custreams_iread(fields::GGArray...)
    	    if length(fields) > size(custreams,2)  # Note: for simplicity, we create a stream for every field even if it is not a CuArray
                custreams = [custreams [CuStream(; flags=CUDA.STREAM_NON_BLOCKING, priority=CUDA.priority_range()[end]) for n=1:NNEIGHBORS_PER_DIM, i=1:(length(fields)-size(custreams,2))]];  # Create (additional) maximum priority nonblocking streams to enable overlap with computation kernels.
            end
        end

        function iread_recvbufs!(n::Integer, dim::Integer, A::CuArray{T}, i::Integer) where T <: GGNumber # Function to be called if A is a GPU array.
            if ol(dim,A) >= 2  # There is only a halo and thus a halo update if the overlap is at least 2...
                if dim == 1 || cudaaware_MPI(dim)  # Use a custom copy kernel for the first dimension to obtain a good copy performance (the CUDA 3-D memcopy does not perform well for this extremely strided case).
                    ranges = recvranges(n, dim, A);
                    nthreads = (dim==1) ? (1, 32, 1) : (32, 1, 1);
                    halosize = [r[end] - r[1] + 1 for r in ranges];
                    nblocks  = Tuple(ceil.(Int, halosize./nthreads));
                    @cuda blocks=nblocks threads=nthreads stream=custreams[n,i] read_x2d!(curecvbuf(n,dim,i,A), A, ranges[1], ranges[2], ranges[3], dim);
                else
                    read_h2d_async!(recvbuf_flat(n,dim,i,A), A, recvranges(n,dim,A), dim, custreams[n,i]);
                end
            end
        end
    end
end

# Return the ranges from A to be sent. It will always return ranges for the dimensions x,y and z even if the A is 1D or 2D (for 2D, the 3rd range is 1:1; for 1D, the 2nd and 3rd range are 1:1).
function sendranges(n::Integer, dim::Integer, A::GGArray)
    if (ol(dim, A) < 2) error("Incoherent arguments: ol(A,dim)<2."); end
    if     (n==2) ixyz_dim = size(A, dim) - (ol(dim, A) - 1);
    elseif (n==1) ixyz_dim = 1            + (ol(dim, A) - 1);
    end
    sendranges      = [1:size(A,1), 1:size(A,2), 1:size(A,3)];  # Initialize with the ranges of A.
    sendranges[dim] = ixyz_dim:ixyz_dim;
    return sendranges
end

# Return the ranges from A to be received. It will always return ranges for the dimensions x,y and z even if the A is 1D or 2D (for 2D, the 3rd range is 1:1; for 1D, the 2nd and 3rd range are 1:1).
function recvranges(n::Integer, dim::Integer, A::GGArray)
    if (ol(dim, A) < 2) error("Incoherent arguments: ol(A,dim)<2."); end
    if     (n==2) ixyz_dim = size(A, dim);
    elseif (n==1) ixyz_dim = 1;
    end
    recvranges      = [1:size(A,1), 1:size(A,2), 1:size(A,3)];  # Initialize with the ranges of A.
    recvranges[dim] = ixyz_dim:ixyz_dim;
    return recvranges
end

# Write to the send buffer on the host from the array on the host (h2h). Note: it works for 1D-3D, as sendranges contains always 3 ranges independently of the number of dimensions of A (see function sendranges).
function write_h2h!(sendbuf::AbstractArray{T}, A::Array{T}, sendranges::Array{UnitRange{T2},1}, dim::Integer) where T <: GGNumber where T2 <: Integer
    ix = (length(sendranges[1])==1) ? sendranges[1][1] : sendranges[1];
    iy = (length(sendranges[2])==1) ? sendranges[2][1] : sendranges[2];
    iz = (length(sendranges[3])==1) ? sendranges[3][1] : sendranges[3];
    if     (dim == 1 && length(ix)==1     && iy == 1:size(A,2) && iz == 1:size(A,3)) memcopy!(view(sendbuf,:), view(view(A,ix, :, :),:));
    elseif (dim == 1 && length(ix)==1     && iy == 1:size(A,2) && length(iz)==1    ) memcopy!(view(sendbuf,:), view(view(A,ix, :,iz),:));
    elseif (dim == 2 && ix == 1:size(A,1) && length(iy)==1     && iz == 1:size(A,3)) memcopy!(view(sendbuf,:), view(view(A, :,iy, :),:));
    elseif (dim == 2 && ix == 1:size(A,1) && length(iy)==1     && length(iz)==1    ) memcopy!(view(sendbuf,:), view(view(A, :,iy,iz),:));
    elseif (dim == 3 && ix == 1:size(A,1) && iy == 1:size(A,2)                     ) memcopy!(view(sendbuf,:), view(view(A, :, :,iz),:));
    elseif (dim == 1 || dim == 2 || dim == 3)                                        memcopy!(view(sendbuf,:), view(view(A,sendranges...),:)); # This general case is slower than the three optimised cases above (the result would be the same, of course).
    end
end

# Read from the receive buffer on the host and store on the array on the host (h2h). Note: it works for 1D-3D, as recvranges contains always 3 ranges independently of the number of dimensions of A (see function recvranges).
function read_h2h!(recvbuf::AbstractArray{T}, A::Array{T}, recvranges::Array{UnitRange{T2},1}, dim::Integer) where T <: GGNumber where T2 <: Integer
    ix = (length(recvranges[1])==1) ? recvranges[1][1] : recvranges[1];
    iy = (length(recvranges[2])==1) ? recvranges[2][1] : recvranges[2];
    iz = (length(recvranges[3])==1) ? recvranges[3][1] : recvranges[3];
    if     (dim == 1 && length(ix)==1     && iy == 1:size(A,2) && iz == 1:size(A,3)) memcopy!(view(view(A,ix, :, :), :), view(recvbuf,:));
    elseif (dim == 1 && length(ix)==1     && iy == 1:size(A,2) && length(iz)==1    ) memcopy!(view(view(A,ix, :,iz), :), view(recvbuf,:));
    elseif (dim == 2 && ix == 1:size(A,1) && length(iy)==1     && iz == 1:size(A,3)) memcopy!(view(view(A, :,iy, :), :), view(recvbuf,:));
    elseif (dim == 2 && ix == 1:size(A,1) && length(iy)==1     && length(iz)==1    ) memcopy!(view(view(A, :,iy,iz), :), view(recvbuf,:));
    elseif (dim == 3 && ix == 1:size(A,1) && iy == 1:size(A,2)                     ) memcopy!(view(view(A, :, :,iz), :), view(recvbuf,:));
    elseif (dim == 1 || dim == 2 || dim == 3)                                        memcopy!(view(view(A,recvranges...),:), view(recvbuf,:)); # This general case is slower than the three optimised cases above (the result would be the same, of course).
    end
end

@static if ENABLE_CUDA
    # Write to the send buffer on the host or device from the array on the device (d2x).
    function write_d2x!(cusendbuf::CuDeviceArray{T}, A::CuDeviceArray{T}, sendrangex::UnitRange{Int64}, sendrangey::UnitRange{Int64}, sendrangez::UnitRange{Int64},  dim::Integer) where T <: GGNumber
        ix = (blockIdx().x-1) * blockDim().x + threadIdx().x + sendrangex[1] - 1
        iy = (blockIdx().y-1) * blockDim().y + threadIdx().y + sendrangey[1] - 1
        iz = (blockIdx().z-1) * blockDim().z + threadIdx().z + sendrangez[1] - 1
        if !(ix in sendrangex && iy in sendrangey && iz in sendrangez) return nothing; end
        if     (dim == 1) cusendbuf[iy,iz] = A[ix,iy,iz];
        elseif (dim == 2) cusendbuf[ix,iz] = A[ix,iy,iz];
        elseif (dim == 3) cusendbuf[ix,iy] = A[ix,iy,iz];
        end
        return nothing
    end

    # Read from the receive buffer on the host or device and store on the array on the device (x2d).
    function read_x2d!(curecvbuf::CuDeviceArray{T}, A::CuDeviceArray{T}, recvrangex::UnitRange{Int64}, recvrangey::UnitRange{Int64}, recvrangez::UnitRange{Int64}, dim::Integer) where T <: GGNumber
        ix = (blockIdx().x-1) * blockDim().x + threadIdx().x + recvrangex[1] - 1
        iy = (blockIdx().y-1) * blockDim().y + threadIdx().y + recvrangey[1] - 1
        iz = (blockIdx().z-1) * blockDim().z + threadIdx().z + recvrangez[1] - 1
        if !(ix in recvrangex && iy in recvrangey && iz in recvrangez) return nothing; end
        if     (dim == 1) A[ix,iy,iz] = curecvbuf[iy,iz];
        elseif (dim == 2) A[ix,iy,iz] = curecvbuf[ix,iz];
        elseif (dim == 3) A[ix,iy,iz] = curecvbuf[ix,iy];
        end
        return nothing
    end

    # Write to the send buffer on the host from the array on the device (d2h).
    function write_d2h_async!(sendbuf::AbstractArray{T}, A::CuArray{T}, sendranges::Array{UnitRange{T2},1}, dim::Integer, custream::CuStream) where T <: GGNumber where T2 <: Integer
        Mem.unsafe_copy3d!(
            pointer(sendbuf), Mem.Host, pointer(A), Mem.Device,
            length(sendranges[1]), length(sendranges[2]), length(sendranges[3]);
            srcPos=(sendranges[1][1], sendranges[2][1], sendranges[3][1]),
            srcPitch=sizeof(T)*size(A,1), srcHeight=size(A,2),
            dstPitch=sizeof(T)*length(sendranges[1]), dstHeight=length(sendranges[2]),
            async=true, stream=custream
        )
    end

    # Read from the receive buffer on the host and store on the array on the device (h2d).
    function read_h2d_async!(recvbuf::AbstractArray{T}, A::CuArray{T}, recvranges::Array{UnitRange{T2},1}, dim::Integer, custream::CuStream) where T <: GGNumber where T2 <: Integer
        Mem.unsafe_copy3d!(
            pointer(A), Mem.Device, pointer(recvbuf), Mem.Host,
            length(recvranges[1]), length(recvranges[2]), length(recvranges[3]);
            dstPos=(recvranges[1][1], recvranges[2][1], recvranges[3][1]),
            srcPitch=sizeof(T)*length(recvranges[1]), srcHeight=length(recvranges[2]),
            dstPitch=sizeof(T)*size(A,1), dstHeight=size(A,2),
            async=true, stream=custream
        )
    end
end

##------------------------------
## FUNCTIONS TO SEND/RECV FIELDS

function irecv_halo!(n::Integer, dim::Integer, A::GGArray, i::Integer; tag::Integer=0)
    req = MPI.REQUEST_NULL;
    if ol(dim,A) >= 2  # There is only a halo and thus a halo update if the overlap is at least 2...
        if cudaaware_MPI(dim) && is_cuarray(A)
            @enable_if_cuda req = MPI.Irecv!(curecvbuf_flat(n,dim,i,A), neighbor(n,dim), tag, comm());
        else
            req = MPI.Irecv!(recvbuf_flat(n,dim,i,A), neighbor(n,dim), tag, comm());
        end
    end
    return req
end

function isend_halo(n::Integer, dim::Integer, A::GGArray, i::Integer; tag::Integer=0)
    req = MPI.REQUEST_NULL;
    if ol(dim,A) >= 2  # There is only a halo and thus a halo update if the overlap is at least 2...
        if cudaaware_MPI(dim) && is_cuarray(A)
            @enable_if_cuda req = MPI.Isend(cusendbuf_flat(n,dim,i,A), neighbor(n,dim), tag, comm());
        else
            req = MPI.Isend(sendbuf_flat(n,dim,i,A), neighbor(n,dim), tag, comm());
        end
    end
    return req
end

function sendrecv_halo_local(n::Integer, dim::Integer, A::GGArray, i::Integer)
    if ol(dim,A) >= 2  # There is only a halo and thus a halo update if the overlap is at least 2...
        if cudaaware_MPI(dim) && is_cuarray(A)
            if n == 1
                @enable_if_cuda cumemcopy!(curecvbuf_flat(2,dim,i,A), cusendbuf_flat(1,dim,i,A));
            elseif n == 2
                @enable_if_cuda cumemcopy!(curecvbuf_flat(1,dim,i,A), cusendbuf_flat(2,dim,i,A));
            end
        else
            if n == 1
                memcopy!(recvbuf_flat(2,dim,i,A), sendbuf_flat(1,dim,i,A));
            elseif n == 2
                memcopy!(recvbuf_flat(1,dim,i,A), sendbuf_flat(2,dim,i,A));
            end
        end
    end
end

function memcopy!(dst::AbstractArray{T}, src::AbstractArray{T}) where T <: GGNumber
	if nthreads() > 1
		@threads for ix = 1:length(dst)  # NOTE: Set the number of threads e.g. as: export JULIA_NUM_THREADS=12
		    @inbounds dst[ix] = src[ix]  # NOTE: We fix here exceptionally the use of @inbounds as this copy between two flat vectors (which must have the right length) is considered safe.
		end
	else
		@inbounds dst .= src
	end
end

@static if ENABLE_CUDA
    function cumemcopy!(dst::CuArray{T}, src::CuArray{T}) where T <: GGNumber
    	@inbounds dst .= src
    end
end

##-------------------------------------------
## FUNCTIONS FOR CHECKING THE INPUT ARGUMENTS

function check_fields(fields::GGArray...)
    # Raise an error if any of the given fields does not have a halo.
    no_halo = Int[];
    for i = 1:length(fields)
        A = fields[i];
        if all([ol(dim, A) < 2 for dim = 1:ndims(A)]) # There is no halo if the overlap is less than 2...
            push!(no_halo, i);
        end
    end
    if length(no_halo) > 1
        error("The fields at positions $(join(no_halo,", "," and ")) have no halo; remove them from the call.")
    elseif length(no_halo) > 0
        error("The field at position $(no_halo[1]) has no halo; remove it from the call.")
    end

    # Raise an error if any of the given fields contains any duplicates.
    duplicates = [[i,j] for i=1:length(fields) for j=i+1:length(fields) if fields[i]===fields[j]];
    if length(duplicates) > 2
        error("The pairs of fields with the positions $(join(duplicates,", "," and ")) are the same; remove any duplicates from the call.")
    elseif length(duplicates) > 0
        error("The field at position $(duplicates[1][2]) is a duplicate of the one at the position $(duplicates[1][1]); remove the duplicate from the call.")
    end

    # Raise an error if not all fields are of the same datatype (restriction comes from buffer handling).
    different_types = [i for i=2:length(fields) if typeof(fields[i])!=typeof(fields[1])];
    if length(different_types) > 1
        error("The fields at positions $(join(different_types,", "," and ")) are of different type than the first field; make sure that in a same call all fields are of the same type.")
    elseif length(different_types) == 1
        error("The field at position $(different_types[1]) is of different type than the first field; make sure that in a same call all fields are of the same type.")
    end
end


text:
function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{St, D},
    ::Type{S},
) where {St <: PSY.Storage, D <: AbstractStorageFormulation, S <: PM.AbstractPowerModel}
    devices = get_available_components(St, sys)

    if !validate_available_devices(St, devices)
        return
    end

    # Variables
    add_variables!(optimization_container, ActivePowerInVariable, devices, D())
    add_variables!(optimization_container, ActivePowerOutVariable, devices, D())
    add_variables!(optimization_container, ReactivePowerVariable, devices, D())
    add_variables!(optimization_container, EnergyVariable, devices, D())

    # Initial Conditions
    initial_conditions!(optimization_container, devices, D())

    # Constraints
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    energy_capacity_constraints!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(optimization_container, devices, model, get_feedforward(model))

    # Energy Balanace limits
    add_constraints!(
        optimization_container,
        EnergyBalanceConstraint,
        EnergyVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )

    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{St, D},
    ::Type{S},
) where {
    St <: PSY.Storage,
    D <: AbstractStorageFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(St, sys)

    if !validate_available_devices(St, devices)
        return
    end

    # Variables
    add_variables!(optimization_container, ActivePowerInVariable, devices, D())
    add_variables!(optimization_container, ActivePowerOutVariable, devices, D())
    add_variables!(optimization_container, EnergyVariable, devices, D())

    # Initial Conditions
    initial_conditions!(optimization_container, devices, D())

    # Constraints
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    energy_capacity_constraints!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(optimization_container, devices, model, get_feedforward(model))

    # Energy Balanace limits
    add_constraints!(
        optimization_container,
        EnergyBalanceConstraint,
        EnergyVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )

    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{St, BookKeepingwReservation},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractPowerModel}
    devices = get_available_components(St, sys)

    if !validate_available_devices(St, devices)
        return
    end

    # Variables
    add_variables!(
        optimization_container,
        ActivePowerInVariable,
        devices,
        BookKeepingwReservation(),
    )
    add_variables!(
        optimization_container,
        ActivePowerOutVariable,
        devices,
        BookKeepingwReservation(),
    )
    add_variables!(
        optimization_container,
        ReactivePowerVariable,
        devices,
        BookKeepingwReservation(),
    )
    add_variables!(
        optimization_container,
        EnergyVariable,
        devices,
        BookKeepingwReservation(),
    )
    add_variables!(
        optimization_container,
        ReserveVariable,
        devices,
        BookKeepingwReservation(),
    )

    # Initial Conditions
    initial_conditions!(optimization_container, devices, BookKeepingwReservation())

    # Constraints
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    energy_capacity_constraints!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(optimization_container, devices, model, get_feedforward(model))

    # Energy Balanace limits
    add_constraints!(
        optimization_container,
        EnergyBalanceConstraint,
        EnergyVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )

    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{St, BookKeepingwReservation},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(St, sys)

    if !validate_available_devices(St, devices)
        return
    end

    # Variables
    add_variables!(
        optimization_container,
        ActivePowerInVariable,
        devices,
        BookKeepingwReservation(),
    )
    add_variables!(
        optimization_container,
        ActivePowerOutVariable,
        devices,
        BookKeepingwReservation(),
    )
    add_variables!(
        optimization_container,
        EnergyVariable,
        devices,
        BookKeepingwReservation(),
    )
    add_variables!(
        optimization_container,
        ReserveVariable,
        devices,
        BookKeepingwReservation(),
    )

    # Initial Conditions
    initial_conditions!(optimization_container, devices, BookKeepingwReservation())

    # Constraints
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    energy_capacity_constraints!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(optimization_container, devices, model, get_feedforward(model))

    # Energy Balanace limits
    add_constraints!(
        optimization_container,
        EnergyBalanceConstraint,
        EnergyVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )

    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{St, EnergyTarget},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractPowerModel}
    devices = get_available_components(St, sys)

    if !validate_available_devices(St, devices)
        return
    end

    # Variables
    add_variables!(optimization_container, ActivePowerInVariable, devices, EnergyTarget())
    add_variables!(optimization_container, ActivePowerOutVariable, devices, EnergyTarget())
    add_variables!(optimization_container, ReactivePowerVariable, devices, EnergyTarget())
    add_variables!(optimization_container, EnergyVariable, devices, EnergyTarget())
    add_variables!(optimization_container, EnergyShortageVariable, devices, EnergyTarget())
    add_variables!(optimization_container, EnergySurplusVariable, devices, EnergyTarget())
    add_variables!(optimization_container, ReserveVariable, devices, EnergyTarget())
    # Initial Conditions
    initial_conditions!(optimization_container, devices, EnergyTarget())

    # Constraints
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    energy_capacity_constraints!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(optimization_container, devices, model, get_feedforward(model))

    # Energy Balanace limits
    add_constraints!(
        optimization_container,
        EnergyBalanceConstraint,
        EnergyVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    energy_target_constraint!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )

    # Cost Function
    cost_function!(optimization_container, devices, model, S, get_feedforward(model))

    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{St, EnergyTarget},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(St, sys)

    if !validate_available_devices(St, devices)
        return
    end

    # Variables
    add_variables!(optimization_container, ActivePowerInVariable, devices, EnergyTarget())
    add_variables!(optimization_container, ActivePowerOutVariable, devices, EnergyTarget())
    add_variables!(optimization_container, EnergyVariable, devices, EnergyTarget())
    add_variables!(optimization_container, EnergyShortageVariable, devices, EnergyTarget())
    add_variables!(optimization_container, EnergySurplusVariable, devices, EnergyTarget())
    add_variables!(optimization_container, ReserveVariable, devices, EnergyTarget())
    # Initial Conditions
    initial_conditions!(optimization_container, devices, EnergyTarget())

    # Constraints
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    energy_capacity_constraints!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(optimization_container, devices, model, get_feedforward(model))

    # Energy Balanace limits
    add_constraints!(
        optimization_container,
        EnergyBalanceConstraint,
        EnergyVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    energy_target_constraint!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )

    # Cost Function
    cost_function!(optimization_container, devices, model, S, get_feedforward(model))

    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{St, BatteryAncillaryServices},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractPowerModel}
    devices = get_available_components(St, sys)

    if !validate_available_devices(St, devices)
        return
    end

    # Variables
    add_variables!(
        optimization_container,
        ActivePowerInVariable,
        devices,
        BatteryAncillaryServices(),
    )
    add_variables!(
        optimization_container,
        ActivePowerOutVariable,
        devices,
        BatteryAncillaryServices(),
    )
    add_variables!(
        optimization_container,
        ReactivePowerVariable,
        devices,
        BatteryAncillaryServices(),
    )
    add_variables!(
        optimization_container,
        EnergyVariable,
        devices,
        BatteryAncillaryServices(),
    )
    add_variables!(
        optimization_container,
        ReserveVariable,
        devices,
        BatteryAncillaryServices(),
    )
    # Initial Conditions
    initial_conditions!(optimization_container, devices, BatteryAncillaryServices())

    # Constraints
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    energy_capacity_constraints!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(optimization_container, devices, model, get_feedforward(model))

    # Energy Balanace limits
    add_constraints!(
        optimization_container,
        EnergyBalanceConstraint,
        EnergyVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    reserve_contribution_constraint!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )

    return
end

function construct_device!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    model::DeviceModel{St, BatteryAncillaryServices},
    ::Type{S},
) where {St <: PSY.Storage, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(St, sys)

    if !validate_available_devices(St, devices)
        return
    end

    # Variables
    add_variables!(
        optimization_container,
        ActivePowerInVariable,
        devices,
        BatteryAncillaryServices(),
    )
    add_variables!(
        optimization_container,
        ActivePowerOutVariable,
        devices,
        BatteryAncillaryServices(),
    )
    add_variables!(
        optimization_container,
        EnergyVariable,
        devices,
        BatteryAncillaryServices(),
    )
    add_variables!(
        optimization_container,
        ReserveVariable,
        devices,
        BatteryAncillaryServices(),
    )
    # Initial Conditions
    initial_conditions!(optimization_container, devices, BatteryAncillaryServices())

    # Constraints
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        optimization_container,
        RangeConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    energy_capacity_constraints!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(optimization_container, devices, model, get_feedforward(model))

    # Energy Balanace limits
    add_constraints!(
        optimization_container,
        EnergyBalanceConstraint,
        EnergyVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    reserve_contribution_constraint!(
        optimization_container,
        devices,
        model,
        S,
        get_feedforward(model),
    )

    return
end


text:
;; cycle-placement -- toggle a window's placement between multiple candidate
;; positions.
;;
;; Copyright 2014 Cory Echols <cory@sixtimesnine.org>

(define-structure cycle-placement
	(export
		cycle-window-placement-on-current-head
		move-resize-frame-by-percent
		move-resize-frame-by-percent-on-current-head
	)

	(open rep
		rep.system
		sawfish.wm.state.maximize
		sawfish.wm.windows
		sawfish.wm.misc
		sawfish.wm.util.rects
	)

	;; Standard reduce thing
	(defun reduce (f seq start)
		(cond
			((equal (length seq) 0) start)
			(t (reduce f (cdr seq) (funcall f start (car seq))))
		)
	)

	;; Find dock windows on the current window's head.
	(defun docks-on-current-head (win)
		(filter-windows (lambda (w)
			(and
				(equal (current-head win) (current-head w))
				(window-get w 'dock-type))))
	)

	;; Return a new size for window w, so that it will not overlap any dock
	;; windows within the rectangle defined by (x y width height)
	(defun bump-window-position-for-docks (win x y width height)
		(let* ((win-rect (list x y (+ x width) (+ y height)))
			   (largest-rect (reduce
							   find-largest-subrect
							   (rectangles-from-windows (docks-on-current-head (input-focus)))
							   win-rect)
			   )
			   (left (nth 0 largest-rect))
			   (top (nth 1 largest-rect))
			   (width (- (nth 2 largest-rect) left))
			   (height (- (nth 3 largest-rect) top))
			  )
			(list left top width height)
		)
	)

	(defun max-by-weight-func-helper (f seq max-weight max-elem)
		(cond
			((equal (length seq) 0) max-elem)
			((> max-weight (funcall f (car seq))) (max-by-weight-func-helper f (cdr seq) max-weight max-elem))
			(t (max-by-weight-func-helper f (cdr seq) (funcall f (car seq)) (car seq)))
		)
	)

	;; Apply f to each elem in seq.  Return the elem in seq that results in the
	;; max value for f.
	(defun max-by-weight-func (f seq)
		(cond
			((equal (length seq) 0) seq)
			(t (max-by-weight-func-helper f (cdr seq) (funcall f (car seq)) (car seq)))
		)
	)

	;; Find the largest (by area) rectangle that would result from
	;; removing r1 from r2
	(defun find-largest-subrect (r1 r2)
		(max-by-weight-func rectangle-area (rect-minus r1 r2))
	)

	;; Move and resize the frame to the supplied position and size
	(defun move-resize-frame-to (window x y width height)
		(move-window-to window x y)
		(resize-frame-to window width height)
	)

	;; Resize the window's frame to the provided width and height.
	(defun resize-frame-to (window width height)
		(let* ((dim-difference (window-frame-difference window))
			   (w-diff (car dim-difference))
			   (h-diff (cdr dim-difference))
			   (new-width (- width w-diff))
			   (new-height (- height h-diff)))
			(unmaximize-window window)
			(resize-window-to window new-width new-height))
	)

	;; Calculate how much frame decorations add to a window's size.
	(defun window-frame-difference (window)
		(let* ((win-dim (window-dimensions window))
			  (win-w (car win-dim))
			  (win-h (cdr win-dim))
			  (frame-dim (window-frame-dimensions window))
			  (frame-w (car frame-dim))
			  (frame-h (cdr frame-dim)))
			(cons (- frame-w win-w) (- frame-h win-h)))
	)

	;; Move and resize a frame, where supplied coordinates and size are
	;; a percentage of the current screen dimensions.
	;; Also avoid any docks or panels within (x y width height)
	(defun move-resize-frame-by-percent (window x y width height)
		(let* ((screen-x (car (screen-dimensions)))
			   (screen-y (cdr (screen-dimensions)))
			   (new-x (exact-percent x screen-x))
			   (new-y (exact-percent y screen-y))
			   (new-width (exact-percent width screen-x))
			   (new-height (exact-percent height screen-y))
			   (bumped-position (bump-window-position-for-docks window new-x new-y new-width new-height))
			  )
			(move-resize-frame-to
			  window
			  (nth 0 bumped-position)
			  (nth 1 bumped-position)
			  (nth 2 bumped-position)
			  (nth 3 bumped-position)
			))
	)

	;; Move and resize a frame, where suppplied coordinates and size are a
	;; percentage of the current monitor dimensions
	(defun move-resize-frame-by-percent-on-current-head (window x y width height)
		(let* ((head-x (car (current-head-dimensions window)))
			   (head-y (cdr (current-head-dimensions window)))
			   (head-offset-x (car (current-head-offset window)))
			   (head-offset-y (cdr (current-head-offset window)))
			   (new-x (+ head-offset-x (exact-percent x head-x)))
			   (new-y (+ head-offset-y (exact-percent y head-y)))
			   (new-width (exact-percent width head-x))
			   (new-height (exact-percent height head-y))
			   (bumped-position (bump-window-position-for-docks window new-x new-y new-width new-height))
			  )
			(move-resize-frame-to
			  window
			  (nth 0 bumped-position)
			  (nth 1 bumped-position)
			  (nth 2 bumped-position)
			  (nth 3 bumped-position)
			)
		)
	)

	;; Return 'percent' of 'quantity' as an exact integer.
	;; librep has muliple numeric types, but it looks like sawfish wants
	;; integers for everything.  There's probably a better way to do this.
	(defun exact-percent (percent quantity)
		(inexact->exact (exact->inexact (* ( / percent 100) quantity) ))
	)

	;; Cycle window through the list of candidate positions
	;; _window_ is a window, as would be returned by sawfish functions such as
	;; `(input-focus)`.
    ;;
	;; _list-of-candidates_ is a list-of-lists.  Each of the contained lists is
	;; a candidate position.  The candidate position describes where the window
	;; should be placed as `(x-percent y-percent width-percent height-percent)`.
	;; The first invocation of the function will place the window according to
	;; the first entry in `list-of-candidates`.  Successive invocations will
	;; select from the following entries in the list.
	(defun cycle-window-placement-on-current-head (window candidate-positions)
		(let* ((last-cycled-position (window-get window 'cte-last-cycled-pos 'none))
			   (last-cycled-member (member last-cycled-position candidate-positions))
			   (next-cycled-position (cond
										((< 1 (length last-cycled-member))
											(nth 1 last-cycled-member))
										(t
											(nth 0 candidate-positions)))))
			(window-put window 'cte-last-cycled-pos next-cycled-position)
			(move-resize-frame-by-percent-on-current-head
				window
				(nth 0 next-cycled-position)
				(nth 1 next-cycled-position)
				(nth 2 next-cycled-position)
				(nth 3 next-cycled-position)
			)
		)
	)
)

max_stars_repo_name: ctechols/cycle-window-placement

max_issues_repo_name: ctechols/cycle-window-placement

max_forks_repo_name: ctechols/cycle-window-placement