text:
#
# Utilities
#

"""
```
tau(sim::Clock, check::Symbol, x::Union{Number,Symbol}; m::Module=Main)
tau(check::Symbol, x::Union{Number,Symbol}; m::Module=Main)
```
Compare the current simulation time against a number or a variable.

# Arguments
- `sim::Clock`: clock variable, if not given, it is 𝐶.
- `check::Symbol`: a comparison operator as a symbol like `:>`,
- `x::Union{Number,Symbol}`: a number or a symbolic variable like `:a`,
- `m::Module=Main`: the evaluation scope, if a symbolic variable is given.

# Examples
```jldoctest
julia> using Simulate

julia> tau(:>=, 1)
false
julia> tau(:<, 1)
true
julia> a = 1
1
julia> tau(:<=, :a)
true
```
"""
tau(sim::Clock, check::Symbol, x::Union{Number,Symbol}, m::Module=Main) =
        x isa Number ? eval(check)(tau(sim), x) : eval(check)(tau(sim), Core.eval(m, x))
tau(check::Symbol, x::Union{Number,Symbol}, m::Module=Main) = tau(𝐶, check, x, m)

"""
```
val(a::Union{Number, Symbol}, check::Symbol, x::Union{Number, Symbol}, m::Module=Main)
```
Compare two variables or numbers.

# Examples
```jldoctest
julia> using Simulate

julia> m = @__MODULE__;  # necessary for doctest

julia> val(1, :<=, 2)
true
julia> a = 1
1
julia> val(:a, :<=, 2, m)
true
```
"""
function val(a::Union{Number, Symbol}, check::Symbol, x::Union{Number, Symbol},
             m::Module=Main)
    a = a isa Number ? a : Core.eval(m, a)
    x = x isa Number ? x : Core.eval(m, x)
    eval(check)(a, x)
end

"""
```
@SF(f::Symbol, arg...)
@SF f arg...
```
create a `SimFunction` from arguments f, arg...

# Arguments
- `f::Symbol`: a function given as a symbol, e.g. `:f` if f() is your function,
- `arg...`: further arguments to your function

!!! note
    1. keyword arguments don't work with this macro, use SF instead.
    2. if you give @SF as argument(s) to a function, you must enclose it/them
        in parentheses ( @SF ... ) or ( (@SF ...), (@SF ...) )

# Examples
```@jldoctest
julia> using Simulate

julia> @SF :sin pi
SimFunction(sin, (π,), Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}())
julia> a = 1
1
julia> incra() = global a += 1             # create a simple increment function
incra (generic function with 1 method)
julia> event!((@SF :incra), after, 3)      # schedule an increment after 3 time units
3.0
julia> a
1
julia> run!(𝐶, 5)
"run! finished with 1 clock events, 0 sample steps, simulation time: 5.0"
julia> a
2
julia> event!((@SF :incra), (@tau :>= 8))  # schedule a conditional increment
5.0
julia> run!(𝐶, 5)
"run! finished with 0 clock events, 500 sample steps, simulation time: 10.0"
julia> a
3
julia> event!(((@SF :incra), (@SF :incra)), ((@tau :>= 12), (@val :a :<= 3)))
10.0
julia> run!(𝐶, 5)
"run! finished with 0 clock events, 500 sample steps, simulation time: 15.0"
julia> a
5
```
"""
macro SF(f::QuoteNode, arg...)
    return :( SimFunction( Core.eval(@__MODULE__, $f), $(arg...) ) )
end

"""
```
@SP(id, f::Symbol, input::Channel, output::Channel, arg...)
@SP id f input output arg...
```
create a `SimProcess` from arguments f, arg...

!!! note
    keyword arguments don't work with this macro, use SP instead.
"""
macro SP(id, f::Symbol, input::Channel, output::Channel, arg...)
    return :( SimProcess($id, $f, $input, $output, $(arg...) ) )
end

"""
```
@tau(sim::Clock)
@tau sim
@tau()
@tau
```
return the current simulation time.

# Arguments
- `sim::Clock`: if no clock argument is given, it returns 𝐶's time.
"""
macro tau(sim)
    return :( tau($sim) )
end
macro tau()
    return :( tau(𝐶) )
end

"""
```
@tau(sim::Clock, check::Symbol, val)
@tau sim check val
@tau(check::Symbol, val)
@tau check val
```
create a `SimFunction` comparing current simulation time with a given value or
variable.

# Arguments
- `sim::Clock`: if no clock is given, it compares with 𝐶's time,
- `check::Symbol`: the check operator must be given as a symbol e.g. `:<`,
- `val::Union{Number, QuoteNode}`: a value or a symbolic variable,

!!! note
    If you give @tau as argument(s) to a function, you must enclose it/them
    in parentheses ( @tau ... ) or ( (@tau ...), (@tau ...) )!

# Examples
```jldoctest
julia> using Simulate

julia> reset!(𝐶)
"clock reset to t₀=0.0, sampling rate Δt=0.0."
julia> s = @tau :≥ 100
SimFunction(Simulate.tau, (:≥, 100, Main), Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}())
julia> Simulate.simExec(s)
false
julia> Simulate.simExec(@tau < 100)                ### wrong !!
ERROR: syntax: "<" is not a unary operator
julia> Simulate.simExec(@tau :< 100)               ### correct
true
julia> a = 1
1
julia> Simulate.simExec(@tau :< :a)
true
julia> event!(SF(()->global a+=1), (@tau :>= 3))   ### create a conditional event
0.0
julia> a
1
julia> run!(𝐶, 5)                                  ### run
"run! finished with 0 clock events, 500 sample steps, simulation time: 5.0"
julia> a
2
```
"""
macro tau(sim, check::Symbol, val::Union{Number, QuoteNode})
    return :( SimFunction(tau, $sim, $check, $val, @__MODULE__) )
end
macro tau(check::QuoteNode, val::Union{Number, QuoteNode})
    return :( SimFunction(tau, $check, $val, @__MODULE__) )
end

"""
```
@val(a, check::QuoteNode, x)
@val a check b
```
Create a Simfunction comparing two values a and b or two symbolic variables
:a and :b. The comparison operator must be given symbolically, e.g. `:≤`.

# Arguments
- `a, b::: a number, expression or symbol
- `check::QuoteNode`: a comparison operator as a symbol like `:≤`
- `, m::Module=Main`: a module scope for evaluation of given symbolic variables

!!! note
    If you give @val as argument(s) to a function, you must enclose it/them
    in parentheses ( @val ... ) or e.g. ( (@tau ...), (@val ...) )!

# Examples
```jldoctest
julia> using Simulate

julia> reset!(𝐶)
"clock reset to t₀=0.0, sampling rate Δt=0.0."
julia> @val 1 :≤ 2
SimFunction(Simulate.val, (1, :≤, 2, Main), Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}())
julia> Simulate.simExec(@val 1 :≤ 2)
true
julia> a = 1
1
julia> Simulate.simExec(@val :a :≤ 2)
true
julia> event!(SF(()->global a+=1), ((@tau :>= 3), (@val :a :<= 3))) ### a conditional event
0.0
julia> run!(𝐶, 5)
"run! finished with 0 clock events, 500 sample steps, simulation time: 5.0"
julia> a
2
```
"""
macro val(a, check::QuoteNode, x)
    return :( SimFunction(val, $a, $check, $x, @__MODULE__) )
end


text:
function ZipfDistribution(arg0::jint, arg1::jdouble)
    return ZipfDistribution((jint, jdouble), arg0, arg1)
end

function cumulative_probability(obj::ZipfDistribution, arg0::jint)
    return jcall(obj, "cumulativeProbability", jdouble, (jint,), arg0)
end

function get_exponent(obj::ZipfDistribution)
    return jcall(obj, "getExponent", jdouble, ())
end

function get_number_of_elements(obj::ZipfDistribution)
    return jcall(obj, "getNumberOfElements", jint, ())
end

function get_numerical_mean(obj::ZipfDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::ZipfDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_support_lower_bound(obj::ZipfDistribution)
    return jcall(obj, "getSupportLowerBound", jint, ())
end

function get_support_upper_bound(obj::ZipfDistribution)
    return jcall(obj, "getSupportUpperBound", jint, ())
end

function is_support_connected(obj::ZipfDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end

function log_probability(obj::ZipfDistribution, arg0::jint)
    return jcall(obj, "logProbability", jdouble, (jint,), arg0)
end

function probability(obj::ZipfDistribution, arg0::jint)
    return jcall(obj, "probability", jdouble, (jint,), arg0)
end

ext: jl


text:
"""Функции для численного дифференцирования с использованием центральных разностей"""
module NumDiff
"""
Нахождение частной производной первого порядка
от векторозначной функции многих переменных

``derivative(func, [t,] x, k; p, h, order) -> df``

Для вычисления производных используется метод центральных разностей.

# Аргументы

* `func::Function`: функция в виде func(t, x) или func(x).
* `t::Float64`: обциональный скалярный аргумент функции (время).
* `x::Vector{Float64}`: векторный аргумент функции `x`.
* `k::Int`: номер переменной, по которой необходимо взять производную.
* `h::Float64`: шаг приращения, который необходимо использовать.
* `order::Int`: порядок точности {2, 4, 6, 8}.

# Возвращаемые значения

`df::Vector{Float64}`: вектор частных производных.
"""
function derivative(func::Function, t::Float64, x::Vector{Float64}; k=1, h=0.1, order=4)
  # Используются центральные конечные разности коэффициенты этих разностей задаются массивом
  # Заполняем массив коэффициентов
  a = [[-0.5, 0.0, 0.5], # Порядка 2
       [1.0/12.0, -2.0/3.0, 0.0, 2.0/3.0, -1.0/12.0], # Порядка 4
       [-1.0/60.0, 3.0/20.0, -3.0/4.0, 0.0, 3.0/4.0, -3.0/20.0, 1.0/60.0], # Порядка 6
       [1.0/280.0, -4.0/105.0, 1.0/5.0, -4.0/5.0, 0.0, 4.0/5.0, -1.0/5.0, 4.0/105.0, -1.0/280.0]] # Порядка 8

  # Проверяем по возвращаемому значению,
  # является ли функция векторной или скалярной
  if isa(func(t, x), Vector{Float64})
    df = zeros(func(t, x))
  elseif isa(func(t, x), Float64)
    df = 0.0
  end

  Ih = h * eye(Float64, length(x))

  if order == 2
    for j = 0:order
      df = df + a[1][j+1] * func(t, x + (j - 1) * Ih[:, k])
    end
  elseif order == 4
    for j = 0:order
      df = df + a[2][j+1] * func(t, x + (j - 2) * Ih[:, k])
    end
  elseif order == 6
    for j = 0:order
      df = df + a[3][j+1] * func(t, x + (j - 3) * Ih[:, k])
    end
  elseif order == 8
    for j = 0:order
      df = df + a[4][j+1] * func(t, x + (j - 4) * Ih[:, k])
    end
  else
    println(STDERR, "Variable order = {2, 4, 6, 8}! Using order=4")
    for j = 0:order
      df = df + a[2][j+1] * func(t, x + (j - 2) * Ih[:, k])
    end
  end

  return df / h
end

# Частные производные от функции без параметра t
function derivative(func::Function, x::Vector{Float64}; k=1, h=0.1, order=4)
  # Заполняем массив коэффициентов
  a = [[-0.5, 0.0, 0.5], # Порядка 2
       [1.0/12.0, -2.0/3.0, 0.0, 2.0/3.0, -1.0/12.0], # Порядка 4
       [-1.0/60.0, 3.0/20.0, -3.0/4.0, 0.0, 3.0/4.0, -3.0/20.0, 1.0/60.0], # Порядка 6
       [1.0/280.0, -4.0/105.0, 1.0/5.0, -4.0/5.0, 0.0, 4.0/5.0, -1.0/5.0, 4.0/105.0, -1.0/280.0]] # Порядка 8

  # Проверяем по возвращаемому значению,
  # является ли функция векторной или скалярной
  if isa(func(x), Vector{Float64})
    df = zeros(func(x))
  elseif isa(func(x), Float64)
    df = 0.0
  end

  Ih = h * eye(Float64, length(x))

  if order == 2
    for j = 0:order
      df = df + a[1][j+1] * func(x + (j - 1) * Ih[:, k])
    end
  elseif order == 4
    for j = 0:order
      df = df + a[2][j+1] * func(x + (j - 2) * Ih[:, k])
    end
  elseif order == 6
    for j = 0:order
      df = df + a[3][j+1] * func(x + (j - 3) * Ih[:, k])
    end
  elseif order == 8
    for j = 0:order
      df = df + a[4][j+1] * func(x + (j - 4) * Ih[:, k])
    end
  else
    println(STDERR, "Variable order = {2, 4, 6, 8}! Using order=4")
    for j = 0:order
      df = df + a[2][j+1] * func(x + (j - 2) * Ih[:, k])
    end
  end

  return df / h
end

# Частные производные от функции одного аргумента без параметра t
function derivative(func::Function, x::Float64; h=0.1, order=4)
  # Заполняем массив коэффициентов
  a = [[-0.5, 0.0, 0.5], # Порядка 2
       [1.0/12.0, -2.0/3.0, 0.0, 2.0/3.0, -1.0/12.0], # Порядка 4
       [-1.0/60.0, 3.0/20.0, -3.0/4.0, 0.0, 3.0/4.0, -3.0/20.0, 1.0/60.0], # Порядка 6
       [1.0/280.0, -4.0/105.0, 1.0/5.0, -4.0/5.0, 0.0, 4.0/5.0, -1.0/5.0, 4.0/105.0, -1.0/280.0]] # Порядка 8

  # Проверяем по возвращаемому значению,
  # является ли функция векторной или скалярной
  if isa(func(x), Vector{Float64})
    df = zeros(func(x))
  elseif isa(func(x), Float64)
    df = 0.0
  end

  if order == 2
    for j = 0:order
      df = df + a[1][j+1] * func(x + (j - 1) * h)
    end
  elseif order == 4
    for j = 0:order
      df = df + a[2][j+1] * func(x + (j - 2) * h)
    end
  elseif order == 6
    for j = 0:order
      df = df + a[3][j+1] * func(x + (j - 3) * h)
    end
  elseif order == 8
    for j = 0:order
      df = df + a[4][j+1] * func(x + (j - 4) * h)
    end
  else
    println(STDERR, "Variable order = {2, 4, 6, 8}! Using order=4")
    for j = 0:order
      df = df + a[2][j+1] * func(x + (j - 2) * h)
    end
  end

  return df / h
end

"""
Нахождение матрицы Якоби от векторозначной функции многих переменных

``jacobian(func, t, x; h=0.1, order=4) -> J``

# Аргументы

* `func::Function`: функция в виде func(t, x, p).
* `t::Float64`: скалярный аргумент функции (время).
* `x::Vector{Float64}`: векторный аргумент функции `x`.
* `p::Tuple`: параметры функции (даже если их нет).
* `h::Float64`: шаг приращения, который необходимо использовать.
* `order::Int`: порядок точности {2, 4, 6, 8}.

# Возвращаемые значения

* `J::Matrix{Float64}(n, n)`: матрица Якоби.
"""
function jacobian(func::Function, t::Float64, x::Vector{Float64}; h=0.1, order=4)
  n = length(x)
  J = Matrix{Float64}(n, n)

  for i in 1:n
    J[:, i] = derivative(func, t, x; k=i, h=h, order=order)
  end
  return J
end
end # module


text:
module Functionals

export J_T_ss, J_T_sm, J_T_re
export gate_functional, make_gate_chi
export make_gradient, make_chi

import ..WeightedObjective

using LinearAlgebra
using Zygote: Zygote
using FiniteDifferences: FiniteDifferences

@doc raw"""
Average complex overlap of the target states with forward-propagated states.

```julia
f_tau(ϕ, objectives; τ=nothing)
```

calculates

```math
f_τ = \frac{1}{N} \sum_{k=1}^{N} w_k τ_k
```

with

```math
τ_k = ⟨ϕ_k^\tgt|ϕ_k(T)⟩
```

in Hilbert space, or

```math
τ_k = \tr[ρ̂_k^{\tgt\,\dagger} ρ̂_k(T)]
```

in Liouville space, where ``|ϕ_k⟩`` or ``ρ̂_k`` are the elements
of `ϕ`, and ``|ϕ_k^\tgt⟩`` or ``ρ̂_k^\tgt`` are the
target states from the `target_state` field of the `objectives`. If `τ` is
given as a keyword argument, it must contain the values `τ_k` according to the
above definition. Otherwise, the ``τ_k`` values will be calculated internally.

``N`` is the number of objectives, and ``w_k`` is an optional weight for each
objective. For any objective that has a `weight` attribute (cf.
[`WeightedObjective`](@ref)), the ``w_k`` is taken from that attribute;
otherwise, ``w_k = 1``. The weights, if present, are not automatically
normalized, they are assumed to have values such that the resulting ``f_τ``
lies in the unit circle of the complex plane. Usually, this means that the
weights should sum to ``N``.
"""
function f_tau(ϕ, objectives; τ=nothing)
    N = length(objectives)
    if τ === nothing
        τ = [dot(objectives[k].target_state, ϕ[k]) for k = 1:N]
    end
    f::ComplexF64 = 0
    for k = 1:N
        obj = objectives[k]
        w = isa(obj, WeightedObjective) ? obj.weight : 1.0
        f += w * τ[k]
    end
    return f / N
end

@doc raw"""State-to-state phase-insensitive fidelity.

```julia
F_ss(ϕ, objectives; τ=nothing)
```

calculates

```math
F_{\text{ss}} = \frac{1}{N} \sum_{k=1}^{N} w_k |τ_k|^2 \quad\in [0, 1]
```

with ``N``, ``w_k`` and ``τ_k`` as in [`f_tau`](@ref).
"""
function F_ss(ϕ, objectives; τ=nothing)
    N = length(objectives)
    if τ === nothing
        τ = [dot(objectives[k].target_state, ϕ[k]) for k = 1:N]
    end
    F::ComplexF64 = f_tau(ϕ, objectives; τ=abs.(τ) .^ 2)
    @assert imag(F) < 1e-10
    return real(F)
end

@doc raw"""State-to-state phase-insensitive functional.

```julia
J_T_ss(ϕ, objectives; τ=nothing)
```

calculates

```math
J_{T,\text{ss}} = 1 - F_{\text{ss}} \in [0, 1].
```

All arguments are passed to [`F_ss`](@ref).
"""
function J_T_ss(ϕ, objectives; τ=nothing)
    return 1.0 - F_ss(ϕ, objectives; τ=τ)
end

@doc raw"""Krotov-states χ for functional [`J_T_ss`](@ref).

```julia
chi_ss!(χ, ϕ, objectives; τ=nothing)
```

sets the elements of `χ` according to

```math
|χ_k⟩
= -\frac{∂ J_{T,\text{ss}}}{∂ ⟨ϕ_k(T)|}
= \frac{1}{N} w_k τ_k |ϕ^{\tgt}_k⟩
```

with ``|ϕ^{\tgt}_k⟩``, ``τ_k`` and ``w_k`` as defined in [`f_tau`](@ref).

Note: this function can be obtained with `make_chi(J_T_ss, objectives)`.
"""
function chi_ss!(χ, ϕ, objectives; τ=nothing)
    N = length(objectives)
    if τ === nothing
        τ = [dot(objectives[k].target_state, ϕ[k]) for k = 1:N]
    end
    for k = 1:N
        obj = objectives[k]
        ϕₖ_tgt = obj.target_state
        copyto!(χ[k], ϕₖ_tgt)
        w = isa(obj, WeightedObjective) ? obj.weight : 1.0
        lmul!((τ[k] * w) / N, χ[k])
    end
end

@doc raw"""Square-modulus fidelity.

```julia
F_sm(ϕ, objectives; τ=nothing)
```

calculates

```math
F_{\text{sm}} = |f_τ|^2  \quad\in [0, 1].
```

All arguments are passed to [`f_tau`](@ref) to evaluate ``f_τ``.
"""
function F_sm(ϕ, objectives; τ=nothing)
    return abs(f_tau(ϕ, objectives; τ=τ))^2
end

@doc raw"""Square-modulus functional.

```julia
J_T_sm(ϕ, objectives; τ=nothing)
```

calculates

```math
J_{T,\text{sm}} = 1 - F_{\text{sm}} \quad\in [0, 1].
```

All arguments are passed to [`f_tau`](@ref) while evaluating ``F_{\text{sm}}``
in [`F_sm`](@ref).
"""
function J_T_sm(ϕ, objectives; τ=nothing)
    return 1.0 - F_sm(ϕ, objectives; τ=τ)
end

@doc raw"""Gradient for [`J_T_sm`](@ref).

```julia
grad_J_T_sm!(∇J_T, τ, ∇τ)
```

analytically sets the elements of the gradient `∇J_T` according to

```math
\frac{∂ J_{T,\text{sm}}(\{τ_k\})}{\partial ϵ_{ln}}
= \frac{1}{N^2} \sum_{k=1}^N \sum_{k'=1}^N\left[
        \frac{\partial τ_{k'}^*}{∂ϵ_{ln}} τ_k +
        τ_{k'}^* \frac{\partial τ_k^*}{∂ϵ_{ln}}
   \right]
= -\frac{2}{N} \Re \sum_{k=1}^N \sum_{k'=1}^N
  τ_{k'}^* \frac{∂τ_k}{\partial ϵ_{ln}}
```

with all quantities as defined in [`make_gradient`](@ref).

Note: this function can be obtained with
`make_gradient(J_T_sm, objectives, via=:tau)`.
"""
function grad_J_T_sm!(∇J_T, τ, ∇τ)
    N = length(τ) # number of objectives
    L, N_T = size(∇τ[1])  # number of controls/time intervals
    ∇J_T′ = reshape(∇J_T, L, N_T)  # writing to ∇J_T′ modifies ∇J_T
    for l = 1:L
        for n = 1:N_T
            ∇J_T′[l, n] = real(sum([conj(τ[k′]) * ∇τ[k][l, n] for k′ = 1:N for k = 1:N]))
        end
    end
    lmul!(-2 / N^2, ∇J_T)
    return ∇J_T
end

@doc raw"""Krotov-states χ for functional [`J_T_sm`](@ref).

```julia
chi_sm!(χ, ϕ, objectives; τ=nothing)
```

sets the elements of `χ` according to

```math
|χ_k⟩
= -\frac{\partial J_{T,\text{sm}}}{\partial ⟨ϕ_k(T)|}
= \frac{1}{N^2} w_k \sum_{j}^{N} w_j τ_j |ϕ_k^{\tgt}⟩
```

with ``|ϕ^{\tgt}_k⟩``, ``τ_j`` and ``w_k`` as defined in [`f_tau`](@ref).

Note: this function can be obtained with `make_chi(J_T_sm, objectives)`.
"""
function chi_sm!(χ, ϕ, objectives; τ=nothing)

    N = length(objectives)
    if τ === nothing
        τ = [dot(objectives[k].target_state, ϕ[k]) for k = 1:N]
    end

    w = ones(N)
    for k = 1:N
        obj = objectives[k]
        w[k] = isa(obj, WeightedObjective) ? obj.weight : 1.0
    end

    for k = 1:N
        obj = objectives[k]
        ϕₖ_tgt = obj.target_state
        copyto!(χ[k], ϕₖ_tgt)
        lmul!(w[k] * sum(w .* τ) / N^2, χ[k])
    end

end

@doc raw"""Real-part fidelity.

```julia
F_re(ϕ, objectives; τ=nothing)
```

calculates

```math
F_{\text{re}} = \Re[f_{τ}] \quad\in \begin{cases}
    [-1, 1] & \text{in Hilbert space} \\
    [0, 1] & \text{in Liouville space.}
\end{cases}
```

All arguments are passed to [`f_tau`](@ref) to evaluate ``f_τ``.
"""
function F_re(ϕ, objectives; τ=nothing)
    return real(f_tau(ϕ, objectives; τ=τ))
end

@doc raw"""Real-part functional.

```julia
J_T_re(ϕ, objectives; τ=nothing)
```

calculates

```math
J_{T,\text{re}} = 1 - F_{\text{re}} \quad\in \begin{cases}
    [0, 2] & \text{in Hilbert space} \\
    [0, 1] & \text{in Liouville space.}
\end{cases}
```

All arguments are passed to [`f_tau`](@ref) while evaluating ``F_{\text{re}}``
in [`F_re`](@ref).
"""
function J_T_re(ϕ, objectives; τ=nothing)
    return 1.0 - F_re(ϕ, objectives; τ=τ)
end

@doc raw"""Krotov-states χ for functional [`J_T_re`](@ref).

```julia
chi_re!(χ, ϕ, objectives; τ=nothing)
```

sets the elements of `χ` according to

```math
|χ_k⟩
= -\frac{∂ J_{T,\text{re}}}{∂ ⟨ϕ_k(T)|}
= \frac{1}{2N} w_k |ϕ^{\tgt}_k⟩
```

with ``|ϕ^{\tgt}_k⟩`` and ``w_k`` as defined in [`f_tau`](@ref).

Note: this function can be obtained with `make_chi(J_T_re, objectives)`.
"""
function chi_re!(χ, ϕ, objectives; τ=nothing)
    N = length(objectives)
    if τ === nothing
        τ = [dot(objectives[k].target_state, ϕ[k]) for k = 1:N]
    end
    for k = 1:N
        obj = objectives[k]
        ϕₖ_tgt = obj.target_state
        copyto!(χ[k], ϕₖ_tgt)
        w = isa(obj, WeightedObjective) ? obj.weight : 1.0
        lmul!(w / (2N), χ[k])
    end
end

"""Convert a functional from acting on a gate to acting on propagated states.

```
J_T = gate_functional(J_T_U; kwargs...)
```

constructs a functional `J_T` that meets the requirements for
[`make_gradient`](@ref) and [`make_chi`](@ref). That is, the output `J_T` takes
positional positional arguments `ϕ` and `objectives`. The input functional
`J_T_U` is assumed to have the signature `J_T_U(U; kwargs...)` where `U` is a
matrix with elements ``U_{ij} = ⟨Ψ_i|ϕ_j⟩``, where ``|Ψ_i⟩`` is the
`initial_state` of the i'th `objectives` (assumed to be the i'th canonical
basis state) and ``|ϕ_j⟩`` is the result of forward-propagating ``|Ψ_j⟩``. That
is, `U` is the projection of the time evolution operator into the subspace
defined by the basis in the `initial_states` of the  `objectives`.
"""
function gate_functional(J_T_U; kwargs...)

    function J_T(ϕ, objectives; τ=nothing)
        N = length(objectives)
        U = [(objectives[i].initial_state ⋅ ϕ[j]) for i = 1:N, j = 1:N]
        return J_T_U(U; kwargs...)
    end

    return J_T

end

@doc raw"""
Return a function to evaluate ``|χ_k⟩ = -∂J_T(Û)/∂⟨ϕ_k|`` via the chain rule.

```julia
chi! = make_gate_chi(J_T_U, objectives; use_finite_differences=false, kwargs...)
```

returns a function equivalent to

```julia
chi! = make_chi(gate_functional(J_T_U; kwargs...), objectives)
```

```math
\begin{split}
    |χ_k⟩
    &= -\frac{∂}{∂⟨ϕ_k|} J_T \\
    &= - \frac{1}{2} \sum_i (∇_U J_T)_{ik} \frac{∂ U_{ik}}{∂⟨ϕ_k|} \\
    &= - \frac{1}{2} \sum_i (∇_U J_T)_{ik} |Ψ_i⟩
\end{split}
```

where ``|Ψ_i⟩`` is the basis state stored as the `initial_state` of the i'th
`objective`, see [`gate_functional`](@ref).

The gradient ``∇_U J_T`` is obtained via automatic differentiation, or via
finite differences if `use_finite_differences=true`.

Compared to the more general [`make_chi`](@ref), `make_gate_chi` will generally
have a slightly smaller numerical overhead, as it pushes the use of automatic
differentiation down by one level.

With `use_finite_differences=true`, this routine serves to test and debug
gradients for gate functionals obtained by automatic differentiation.
"""
function make_gate_chi(J_T_U, objectives; use_finite_differences=false, kwargs...)

    N = length(objectives)
    basis = [obj.initial_state for obj in objectives]

    function zygote_gate_chi!(χ, ϕ, objectives; τ=nothing)
        function _J_T(U)
            -J_T_U(U; kwargs...)
        end
        U = [basis[i] ⋅ ϕ[j] for i = 1:N, j = 1:N]
        if use_finite_differences
            fdm = FiniteDifferences.central_fdm(5, 1)
            ∇J = FiniteDifferences.grad(fdm, gate -> _J_T(gate), U)[1]
        else
            ∇J = Zygote.gradient(gate -> _J_T(gate), U)[1]
        end
        for k = 1:N
            χ[k] .= 0.5 * sum([∇J[i, k] * basis[i] for i = 1:N])
        end
    end

    return zygote_gate_chi!

end

@doc raw"""Gradient for an arbitrary functional evaluated via χ-states.

```julia
grad_J_T_via_chi!(∇J_T, τ, ∇τ)
```

sets the (vectorized) elements of the gradient `∇J_T` to the gradient
``∂J_T/∂ϵ_{ln}`` for an arbitrary functional ``J_T=J_T(\{|ϕ_k(T)⟩\})``, under
the assumption that

```math
\begin{aligned}
    τ_k &= ⟨χ_k|ϕ_k(T)⟩ \quad \text{with} \quad |χ_k⟩ &= -∂J_T/∂⟨ϕ_k(T)|
    \quad \text{and} \\
    ∇τ_{kln} &= ∂τ_k/∂ϵ_{ln}\,,
\end{aligned}
```

where ``|ϕ_k(T)⟩`` is a state resulting from the forward propagation of some
initial state ``|ϕ_k⟩`` under the pulse values ``ϵ_{ln}`` where ``l`` numbers
the controls and ``n`` numbers the time slices. The ``τ_k`` are the elements of
`τ` and ``∇τ_{kln}`` corresponds to `∇τ[k][l, n]`.

In this case,

```math
(∇J_T)_{ln} = ∂J_T/∂ϵ_{ln} = -2 \Re \sum_k ∇τ_{kln}\,,
```

see [`make_gradient`](@ref).

Note that the definition of the ``|χ_k⟩`` matches exactly the definition
of the boundary condition for the backward propagation in Krotov's method, see
[`make_chi`](@ref). Specifically, there is a minus sign in front of the
derivative, compensated by the minus sign in the factor ``(-2)`` of the final
``(∇J_T)_{ln}``.
"""
function grad_J_T_via_chi!(∇J_T, τ, ∇τ)
    N = length(τ) # number of objectives
    L, N_T = size(∇τ[1])  # number of controls/time intervals
    ∇J_T′ = reshape(∇J_T, L, N_T)  # writing to ∇J_T′ modifies ∇J_T
    for l = 1:L
        for n = 1:N_T
            ∇J_T′[l, n] = real(sum([∇τ[k][l, n] for k = 1:N]))
        end
    end
    lmul!(-2, ∇J_T)
    return ∇J_T
end

@doc raw"""Return a function that evaluates the gradient ``∇J_T``.

```julia
grad_func! = make_gradient(J_T, objectives; via=:tau, force_zygote=false)
```

creates a function `gradfunc!(∇J_T, τ, ∇τ)` that takes a vector `τ` of values
``τ_k`` and a vector `∇τ` of gradients ``∇τ_k`` where the ``(ln)``'th element
of ``∇τ_k`` is ``∂τ_k/∂ϵ_{ln}``, and writes the element ``(ln)`` of the
gradient `∇J_T` as ``∂J_T/∂ϵ_{ln}``. The definition of ``τ_k`` depends on
`via`, see below. The ``ϵ_{ln}`` are the values of the control control field
discretized to the midpoints of a time grid. The index ``l`` numbers the
control and ``n`` numbers the time slice. The gradient (like the controls) are
assumed to be vectorized. That is, ``∇J_T`` is a vector of values with a
double-index ``(ln)``.

The passed `J_T` parameter corresponding to the functional ``J_T`` must be a
function that takes a vector of forward-propagates states `ϕ` and a vector of
objectives as positional parameters.  It must also accept a vector `τ` as a
keyword argument, which contains the overlaps of the states in `ϕ` and the
`target_state` fields of the `objectives`. If the `objectices` do not define
`target_states`, or if the ``τ``-values are not available, `J_T` must accept
`τ=nothing`. See [`J_T_sm`](@ref) for an example.

## Gradient via τ

For `via=:tau` (default), we define

```math
τ_k ≡ ⟨ϕ_k^\tgt|ϕ_k(T)⟩
```

as the overlap of ``|ϕ_k(T)⟩`` resulting from the forward propagation of the
`initial_state` ``|ϕ_k⟩``  of the k'th objective under the pulse values
``ϵ_{ln}``, and ``|ϕ_k^\tgt⟩`` as the `target_state` of the k'th objective.

We then understand ``J_T`` as a function of the ``τ_k``, and evaluate the
elements of ``∇J_T`` via the chain rule:

```math
(∇J_T)_{ln} ≡ \frac{∂J_T(\{τ_k\})}{∂ϵ_{ln}}
= 2\Re\sum_k
    \frac{∂J_T}{∂τ_k}
    \frac{∂τ_k}{∂ϵ_{ln}}\,.
```

Since the ``τ_k`` are complex numbers,

```math
\frac{∂J_T}{∂τ_k} = \frac{1}{2}\left(
    \frac{∂ J_T}{∂ \Re[τ_k]}
    - i \frac{∂ J_T}{∂ \Im[τ_k]}
\right)
```

is defined as the [Wirtinger
derivative](https://www.ekinakyurek.me/complex-derivatives-wirtinger/), and

```math
\frac{∂τ_k}{∂ϵ_{ln}} = \frac{∂\Re[τ_k]}{∂ϵ_{ln}} + i \frac{∂\Im[τ_k]}{∂ϵ_{ln}}
```

is simply the derivative of a complex number with respect to the real-valued
``ϵ_{ln}``.

Thus, the returned `grad_func!` effectively encodes the outer derivative
``∂J_T/∂τ_k``. For functionals where that derivative is known analytically,
the analytic expression is used, e.g., [`J_T_sm`](@ref) →
[`grad_J_T_sm!`](@ref).

Otherwise, or if `force_zygote=true`, the outer derivative is
determined directly from `J_T`, via automatic differentiation (using
[Zygote](https://fluxml.ai/Zygote.jl)).

## Gradient via χ

For `via=:chi`, the functional ``J_T`` is understood directly as a function of
the forward-propagated states ``|ϕ_k⟩`` instead of a function of overlaps with
the target states. This is useful in particular if the `objectives` do not
define objectives and/or the functional `J_T` cannot be expressed in terms of
overlaps.

Again we apply a chain rule to calculate the elements of the gradient ``∇J_T``:

```math
\begin{split}
(∇J_T)_{ln} &≡ \frac{∂J_T(\{|ϕ_k(T)⟩\})}{∂ϵ_{ln}}\\
&= 2\Re\sum_k
    \frac{∂J_T}{∂|ϕ_k(T)⟩}
    \frac{∂|ϕ_k(T)⟩}{∂ϵ_{ln}} \\
&= -2 \Re \sum_k \frac{∂}{∂ϵ_{ln}} ⟨χ_k(T)|ϕ_k(T)⟩\,,
\end{split}
```

with

```math
|χ_k⟩
= -\frac{J_T}{⟨ϕ_k(T)|}
= -\frac{1}{2}\left(
    \left\vert \frac{∂J_T}{∂\Re[ϕ_k]} \right\rangle
    + i \left\vert \frac{∂J_T}{∂\Im[ϕ_k]} \right\rangle
    \right)
```

as a matrix-calculus extension of the Wirtinger derivative. This definition of
``|χ_k⟩`` (note the minus sign!) matches the definition of the boundary
condition in Krotov's method, and for a given functional `J_T`, the
states ``|χ_k⟩`` can be obtained with [`make_chi`](@ref).

We define

```math
τ_k ≡ ⟨χ_k(T)|ϕ_k(T)⟩
```

and associate the `gradfunc!` argument `∇τ[k][l, n]` with

```math
(∇τ_k)_{ln} = \frac{∂τ_k}{∂ ϵ_{ln}}
```

so that structurally, ``(∇J_T)_{ln}`` is the same as for `via=:tau`, just that
``τ_k`` is now defined with respect to the boundary condition state ``|χ_k⟩``
instead of the target state ``|ϕ_k^\tgt⟩``.

The returned `grad_func!` that encodes the above equations is
[`grad_J_T_via_chi!`](@ref). This is independent of `J_T`, since the dependency
on the functional `J_T` is entirely encoded in the states ``|χ_k(T)⟩``, and
thus the gradient `∇τ`. Also, `force_zygote=true` has no effect for `via=:chi`.
Instead, `force_zygote` should be passed to the underlying [`make_chi`](@ref).

!!! tip

    In order to extend `make_gradient` with an analytic implementation for a
    new `J_T` function, define a new method like so:

    ```julia
    make_gradient(::typeof(J_T_sm), objectives, via::Val{:tau}) = grad_J_T_sm!
    ```

    which links `make_gradient` for [`J_T_sm`](@ref) to [`grad_J_T_sm!`](@ref).
"""
function make_gradient(J_T, objectives; via::Symbol=:tau, force_zygote=false)
    if force_zygote && (via == :tau)
        return make_zygote_gradient(J_T, objectives)
    else
        return make_gradient(J_T, objectives, Val(via))
    end
end

function make_zygote_gradient(J_T, objectives)

    function zygote_gradfunc!(∇J_T, τ, ∇τ)
        ∇J = Zygote.gradient(τ -> J_T(nothing, objectives; τ=τ), τ)[1]
        ∇J_T .= vec(sum(real(∇J) .* real.(∇τ) + imag(∇J) .* imag.(∇τ)))
    end

    return zygote_gradfunc!
end

make_gradient(J_T, objectives, via::Val{:chi}) = grad_J_T_via_chi!
make_gradient(J_T, objectives, via::Val{:tau}) = make_zygote_gradient(J_T, objectives)
make_gradient(::typeof(J_T_sm), objectives, via::Val{:tau}) = grad_J_T_sm!

@doc raw"""Return a function that evaluates ``|χ_k⟩ = -∂J_T/∂⟨ϕ_k|``.

```julia
chi! = make_chi(J_T, objectives; force_zygote=false)
```

creates a function `chi!(χ, ϕ, objectives; τ=nothing)` that sets
``|χ_k⟩ = -∂J_T/∂⟨ϕ_k|``. This is the state used as the boundary condition for
the backward propagation propagation in Krotov's method, as well as GRAPE if
[`grad_J_T_via_chi!`](@ref) is used. It is defined as a [Wirtinger
derivative](https://www.ekinakyurek.me/complex-derivatives-wirtinger/),
see [`make_gradient`](@ref).

The function `J_T` must take a vector of states `ϕ`
and a vector of `objectives` as positional parameters, and a vector `τ` as a
keyword argument, see e.g. [`J_T_sm`](@ref). If all objectives define a
`target_state`, then `τ` will be the overlap of the states `ϕ` with those
target states. The functional `J_T` may or may not use those overlaps.
Likewise, the resulting `chi!` may or may not use the keyword parameter `τ`.

For functionals where ``-∂J_T/∂⟨ϕ_k|`` is known analytically, that analytic
derivative will be returned, e.g.,

* [`J_T_sm`](@ref) → [`chi_sm!`](@ref),
* [`J_T_re`](@ref) → [`chi_re!`](@ref),
* [`J_T_ss`](@ref) → [`chi_ss!`](@ref).

Otherwise, or if `force_zygote=true`, automatic differentiation via Zygote is
used to calculate the derivative directly from `J_T`.

!!! tip

    In order to extend `make_chi` with an analytic implementation for a new
    `J_T` function, define a new method `make_analytic_chi` like so:

    ```julia
    make_analytic_chi(::typeof(J_T_sm), objectives) = chi_sm!
    ```

    which links `make_chi` for [`J_T_sm`](@ref) to [`chi_sm!`](@ref).
"""
function make_chi(J_T, objectives; force_zygote=false)
    if force_zygote
        return make_zygote_chi(J_T, objectives)
    else
        return make_analytic_chi(J_T, objectives)
    end
end

function make_zygote_chi(J_T, objectives)

    N = length(objectives)

    function zygote_chi!(χ, ϕ, objectives; τ=nothing)
        function _J_T(Ψ...)
            -J_T(Ψ, objectives)
        end
        for (k, ∇J) ∈ enumerate(Zygote.gradient(_J_T, ϕ...))
            copyto!(χ[k], ∇J)
            lmul!(0.5, χ[k])
        end
    end

    return zygote_chi!

end

make_analytic_chi(J_T, objectives) = make_zygote_chi(J_T, objectives)
# Well, the above fallback to `make_zygote_chi` is clearly not "analytic", but
# the intent of this name is to allow users to define new analytic
# implementation, cf. the explanation in the doc of `make_chi`.
make_analytic_chi(::typeof(J_T_sm), objectives) = chi_sm!
make_analytic_chi(::typeof(J_T_re), objectives) = chi_re!
make_analytic_chi(::typeof(J_T_ss), objectives) = chi_ss!

end


text:
export TemperatureSensor, getTemperatureSensors, getTemperatureSensor, numChannels, getTemperatures, getTemperature

abstract type TemperatureSensor <: Device end

include("DummyTemperatureSensor.jl")
include("ArduinoTemperatureSensor.jl")
#include("FOTemp.jl")

Base.close(t::TemperatureSensor) = nothing

@mustimplement numChannels(sensor::TemperatureSensor)
@mustimplement getTemperatures(sensor::TemperatureSensor)
@mustimplement getTemperature(sensor::TemperatureSensor, channel::Int)

getTemperatureSensors(scanner::MPIScanner) = getDevices(scanner, TemperatureSensor)
getTemperatureSensor(scanner::MPIScanner) = getDevice(scanner, TemperatureSensor)


text:
module SweepOperator
export sweep!

using LinearAlgebra
import LinearAlgebra: BlasFloat, checksquare

const AMat = AbstractMatrix
const AVec = AbstractVector

"""
    sweep!(A, k ; inv=false)
    sweep!(A, ks; inv=false)

Perform the sweep operation (or inverse sweep if `inv=true`) on matrix `A` on element `k`
(or each element in `ks`).  Only the upper triangle is read/swept.

# Example:

    x = randn(100, 10)
    xtx = x'x
    sweep!(xtx, 1)
    sweep!(xtx, 1, true)
"""
function sweep!(A::AMat, k::Integer, inv::Bool = false)
    sweep_with_buffer!(Vector{eltype(A)}(undef, size(A, 2)), A, k, inv)
end

function sweep_with_buffer!(akk::AVec{T}, A::AMat{T}, k::Integer, inv::Bool = false) where
        {T<:BlasFloat}
    # ensure @inbounds is safe
    p = checksquare(A)
    p == length(akk) || throw(DimensionError("incorrect buffer size"))
    @inbounds d = one(T) / A[k, k]  # pivot
    # get column A[:, k] (hack because only upper triangle is available)
    for j in 1:k
        @inbounds akk[j] = A[j, k]
    end
    for j in (k+1):p
        @inbounds akk[j] = A[k, j]
    end
    BLAS.syrk!('U', 'N', -d, akk, one(T), A)  # everything not in col/row k
    rmul!(akk, d * (-one(T)) ^ inv)
    for i in 1:(k-1)  # col k
        @inbounds A[i, k] = akk[i]
    end
    for j in (k+1):p  # row k
        @inbounds A[k, j] = akk[j]
    end
    @inbounds A[k, k] = -d  # pivot element
    A
end

function sweep!(A::AMat{T}, ks::AVec{I}, inv::Bool = false) where {T<:BlasFloat, I<:Integer}
    akk = zeros(T, size(A, 1))
    for k in ks
        sweep_with_buffer!(akk, A, k, inv)
    end
    A
end
function sweep_with_buffer!(akk::AVec{T}, A::AMat{T}, ks::AVec{I}, inv::Bool = false) where
        {T<:BlasFloat, I<:Integer}
    for k in ks
        sweep_with_buffer!(akk, A, k, inv)
    end
    A
end

end # module


text:
"""
Similar to `filelist`, `filelistall(dir)` returns all files in the directory (**as well as its subdirectory**) `dir`.
"""
function filelistall(dir)
    wd = walkdir(dir)
    allpaths = String[]
    for (root, dirs, files) in wd
        for file in files
            path_i = joinpath(root, file)
            push!(allpaths, path_i)
        end
    end
    return allpaths
end

"""
Similar to `filelist`, `filelistall(expr::Regex, dir)` returns files whose name matches `expr` in the directory (**as well as its subdirectory**) `dir`.
"""
function filelistall(expr::Regex, dir)
    allpaths = filelistall(dir)
    desired_ind = occursin.(expr, basename.(allpaths))
    return allpaths[desired_ind]
end


text:
module JungleHelperCobweb

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/Cobweb" Cobweb(x::Integer, y::Integer, sprite::String="")

const placements = Ahorn.PlacementDict(
    "Cobweb (Jungle Helper)" => Ahorn.EntityPlacement(
        Cobweb
    )
)

sprite = "JungleHelper/Cobweb/idle00"

function Ahorn.selection(entity::Cobweb)
    x, y = Ahorn.position(entity)
    return Ahorn.getSpriteRectangle(sprite, x, y)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Cobweb, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
cipher(input::AbstractString) = map(x->isalpha(x)?Char(219-Int(x)):x, lowercase(filter(isalnum, input)))
encode(input::AbstractString) = join(matchall(r"(.{1,5})", cipher(input)), ' ')
decode(input::AbstractString) = cipher(input)


text:
using HorizonSideRobots

function s_peregorodki(r)
    sum = 0
    Z = 0
    goinSWcorner(r)
    while !isborder(r,Nord)
        Z = search(r)
        if sum < Z
            sum = Z
        end
    end
    return sum
end

function goinSWcorner(r)
    sign = true
    while !isborder(r,Sud) || sign
        if !isborder(r,Sud)
            move!(r,Sud)
        end
        if isborder(r,Sud)
            sign = obhod(r,2)
        end
    end
    while !isborder(r,West)
        move!(r,West)
    end
end

function obhod2(r,side)
    a = 0
    b = 0
    x = 0
    sign = true
    while isborder(r,HorizonSide(side)) && !isborder(r,HorizonSide((side + 1) % 4))
        move!(r,HorizonSide((side + 1) % 4))
        x += 1
        a += 1
    end
    if isborder(r,HorizonSide(side)) && isborder(r,HorizonSide((side + 1) % 4))
        sign = false
    end
    if  !isborder(r,HorizonSide(side))
        move!(r,HorizonSide(side))
    end
    while !isborder(r,HorizonSide(side)) && isborder(r,HorizonSide((side + 3) % 4))
        move!(r,HorizonSide(side))
        b += 1
    end
    for i = 1:x
        move!(r,HorizonSide((side+3) % 4))
    end
    return sign, (a)*(b)
end

function obhod(r,side)
    x = 0
    sign = true
    while isborder(r,HorizonSide(side)) && !isborder(r,HorizonSide((side + 3) % 4))
        move!(r,HorizonSide((side + 3) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(side)) && isborder(r,HorizonSide((side + 3) % 4))
        sign = false
    end
    if  !isborder(r,HorizonSide(side))
        move!(r,HorizonSide(side))
    end
    while !isborder(r,HorizonSide(side)) && isborder(r,HorizonSide((side + 1) % 4))
        move!(r,HorizonSide(side))
    end
    for i = 1:x
        move!(r,HorizonSide((side+1) % 4))
    end
    return sign
end

function search(r)
    sign = true
    sum = 0
    Z = 0
    if !isborder(r,Nord)
        sign = true
        move!(r,Nord)
        while !isborder(r,Ost) || sign
            if !isborder(r,Ost)
                move!(r,Ost)
            else
                sign, Z = obhod2(r,3)
                if sum < Z
                    sum = Z
                end
            end
        end
    end
    sign = true
    while !isborder(r, West) || sign
        if !isborder(r, West)
            move!(r, West)
        else
            sign = obhod(r,1)
        end
    end
    return sum
end

s_peregorodki(r)

max_stars_repo_name: Cheburkova/mirea-progs

max_issues_repo_name: Cheburkova/mirea-progs

max_forks_repo_name: Cheburkova/mirea-progs


text:
using CUDA.CUTENSOR
using CUDA
using LinearAlgebra

# using host memory with CUTENSOR doesn't work on Windows
can_pin = !Sys.iswindows()

eltypes = ((Float16, Float16),
            #(Float16, Float32),
            (Float32, Float32),
            #(Float32, Float64),
            (Float64, Float64),
            #(ComplexF16, ComplexF16),
            #(ComplexF16, ComplexF32),
            (ComplexF32, ComplexF32),
            #(ComplexF32, ComplexF64),
            (ComplexF64, ComplexF64))
@testset for N=2:5
    @testset for (eltyA, eltyC) in eltypes
        # setup
        dmax = 2^div(18,N)
        dims = rand(2:dmax, N)
        p = randperm(N)
        indsA = collect(('a':'z')[1:N])
        indsC = indsA[p]
        dimsA = dims
        dimsC = dims[p]
        A = rand(eltyA, dimsA...)
        can_pin && Mem.pin(A)
        dA = CuArray(A)
        dC = similar(dA, eltyC, dimsC...)

        # simple case
        dC = CUTENSOR.permutation!(one(eltyA), dA, indsA, dC, indsC)
        C  = collect(dC)
        @test C == permutedims(A, p) # exact equality
        if can_pin
            Csimple = zeros(eltyC, dimsC...)
            Mem.pin(Csimple)
            Csimple = CUTENSOR.permutation!(one(eltyA), A, indsA, Csimple, indsC)
            synchronize()
            @test Csimple == permutedims(A, p) # exact equality
        end

        # with scalar
        α  = rand(eltyA)
        dC = CUTENSOR.permutation!(α, dA, indsA, dC, indsC)
        C  = collect(dC)
        @test C ≈ α * permutedims(A, p) # approximate, floating point rounding
        if can_pin
            Cscalar = zeros(eltyC, dimsC...)
            Mem.pin(Cscalar)
            Cscalar = CUTENSOR.permutation!(α, A, indsA, Cscalar, indsC)
            synchronize()
            @test Cscalar ≈ α * permutedims(A, p) # approximate, floating point rounding
        end
    end
end


text:
module IOUtils

include("misc.jl")
include("boxes.jl")
include("indents.jl")
include("displays.jl")
include("columns.jl")

export redirect_output, print_boxed, indent, @display,
    horizontal_line, columns, display_matrix

end # module


text:
"""
`daqaddinput(dev, ...)`

Add channels that should be acquired.
"""
daqaddinput(dev::AbstractDAQ)=error("Not implemented for AbstractDAQ")

"""
`daqacquire(dev)`

Start a synchronous data acquisition run.
"""
daqacquire(dev::AbstractDAQ) = error("Not implemented for AbstractDAQ")
daqacquire!(dev::AbstractDAQ) = error("Not implemented for AbstractDAQ")

"""
`daqstart(dev, ...)`

Initiate a data acquisition run asyncrhonously.
"""
daqstart(dev::AbstractDAQ, usethread=false) =
    error("Not implemented for AbstractDAQ")

"""
`daqread(dev)`

Wait to finish the data acquisition run and return the data.
"""
daqread(dev::AbstractDAQ) = error("Not implemented for AbstractDAQ")
daqread!(dev::AbstractDAQ) = error("Not implemented for AbstractDAQ")

"""
`daqstop(dev)`

Stop asynchronous data acquisition.
"""
daqstop(dev::AbstractDAQ) = error("Not implemented for AbstractDAQ")

"""
`daqreference(dev)`

Use a measurement point as a reference. Specific channels can be specified.

"""
daqreference(dev::AbstractDAQ) =error("Not implemented for AbstractDAQ")

"""
`daqconfig(dev; rate, nsamples, time, avg=1)`

Generic configuration of data acquisition. Different devices might
have other capabilities and different terminologies. To use the device specific
parameters and terminology, use function [`daqconfigdev`](@ref).

In this generic interface, the following keyword parameters are allowed:

 * `rate` or `dt` (only one of them)
    - `rate` Sampling rate in Hz
    - `dt` Sampling time in s
 * `nsamples` or `time` (only one of them)
    - `nsamples` Number of samples to be read. 0 usually means continous reading
    - `time` sampling time in seconds
 * `avg` Number of samples that should be read and averaged for each output.
 * `trigger` An integer specifying the trigger type 0 - internal trigger, other values depend on the specific device.
"""
daqconfig(dev::AbstractDAQ; kw...) =
    error("Not implemented for AbstractDAQ")

"""
`daqconfigdev(dev; kw...)`

Device configuration.

Does the samething as [`daqconfig`](@ref) but uses the devices terminology and exact
parameters.
"""
daqconfigdev(dev::AbstractDAQ; kw...) =
    error("Not implemented for AbstractDAQ")

"""
`daqzero(dev)`

Perform a zero calibration of the DAQ device. The exact nature of this zero calibration.
"""
daqzero(dev::AbstractDAQ) =
    error("Not implemented for AbstractDAQ")

"""
`samplesread(dev)`

Return the number of samples read since the beginning of data aquisition.
"""
samplesread(dev::AbstractDAQ) =
    error("Not implemented for AbstractDAQ")

"""
`isreading(dev)`

Returns `true` if data acquisition is ongoing, `false` otherwise.
"""
isreading(dev::AbstractDAQ) =
    error("Not implemented for AbstractDAQ")

"""
`isdaqfinished(dev)`

Returns true if the device has completed its operations.
"""
isdaqfinished(dev::AbstractDAQ) =
    error("Not implemented for AbstractDAQ")
"""
`issamplesavailable(dev)`

Are samples available for reading?
"""
issamplesavailable(dev::AbstractDAQ) =
    error("Not implemented for AbstractDAQ")

"""
`numchannels(dev)`

Number of channels available or configured in the DAQ device.
"""
numchannels(dev::AbstractDAQ) =
    error("Not implemented for AbstractDAQ")

"""
`daqchannels(dev)`

Returns the DAQ channels available or configured in the DAQ device.
"""
daqchannels(dev::AbstractDAQ) =
    error("Not implemented for AbstractDAQ")


text:
#########################################################################
# PyVector: no-copy wrapping of a Julia object around a Python sequence

"""
    PyVector(o::PyObject)

This returns a PyVector object, which is a wrapper around an arbitrary Python list or sequence object.

Alternatively, `PyVector` can be used as the return type for a `pycall` that returns a sequence object (including tuples).
"""
mutable struct PyVector{T} <: AbstractVector{T}
    o::PyObject
    function PyVector{T}(o::PyObject) where T
        if ispynull(o)
            throw(ArgumentError("cannot make PyVector from NULL PyObject"))
        end
        new{T}(o)
    end
end

PyVector(o::PyObject) = PyVector{PyAny}(o)
PyObject(a::PyVector) = a.o
convert(::Type{PyVector}, o::PyObject) = PyVector(o)
convert(::Type{PyVector{T}}, o::PyObject) where {T} = PyVector{T}(o)
unsafe_convert(::Type{PyPtr}, a::PyVector) = PyPtr(a.o)
PyVector(a::PyVector) = a
PyVector(a::AbstractVector{T}) where {T} = PyVector{T}(array2py(a))

# when a PyVector is copied it is converted into an ordinary Julia Vector
similar(a::PyVector, T, dims::Dims) = Array{T}(dims)
similar(a::PyVector{T}) where {T} = similar(a, pyany_toany(T), size(a))
similar(a::PyVector{T}, dims::Dims) where {T} = similar(a, pyany_toany(T), dims)
similar(a::PyVector{T}, dims::Int...) where {T} = similar(a, pyany_toany(T), dims)
eltype(::PyVector{T}) where {T} = pyany_toany(T)
eltype(::Type{PyVector{T}}) where {T} = pyany_toany(T)

size(a::PyVector) = (length(a.o),)

getindex(a::PyVector) = getindex(a, 1)
getindex(a::PyVector{T}, i::Integer) where {T} = convert(T, PyObject(@pycheckn ccall((@pysym :PySequence_GetItem), PyPtr, (PyPtr, Int), a, i-1)))

setindex!(a::PyVector, v) = setindex!(a, v, 1)
function setindex!(a::PyVector, v, i::Integer)
    @pycheckz ccall((@pysym :PySequence_SetItem), Cint, (PyPtr, Int, PyPtr), a, i-1, PyObject(v))
    v
end

summary(a::PyVector{T}) where {T} = string(Base.dims2string(size(a)), " ",
                                          string(pyany_toany(T)), " PyVector")

splice!(a::PyVector, i::Integer) = splice!(a.o, i)
function splice!(a::PyVector{T}, indices::AbstractVector{I}) where {T,I<:Integer}
    v = pyany_toany(T)[a[i] for i in indices]
    for i in sort(indices, rev=true)
        @pycheckz ccall((@pysym :PySequence_DelItem), Cint, (PyPtr, Int), a, i-1)
    end
    v
end
pop!(a::PyVector) = pop!(a.o)
popfirst!(a::PyVector) = popfirst!(a.o)
empty!(a::PyVector) = empty!(a.o)

# only works for List subtypes:
push!(a::PyVector, item) = push!(a.o, item)
insert!(a::PyVector, i::Integer, item) = insert!(a.o, i, item)
pushfirst!(a::PyVector, item) = pushfirst!(a.o, item)
prepend!(a::PyVector, items) = prepend!(a.o, items)
append!(a::PyVector{T}, items) where {T} = PyVector{T}(append!(a.o, items))


text:
struct DominanceSum{Tv, Ti} <: AbstractMatrix{Tv}
    m::Int
    n::Int
    N::Int
    b::Int
    b′::Int
    H::Int
    pos::Vector{Ti}
    qos::Vector{Ti}
    byt::Vector{Ti}
    cnt::Array{Ti, 3}
    wgt::Array{Tv, 2}
    scn::Array{Tv, 3}
end

mutable struct SparseStepwiseDominanceSum{Tv, Ti} <: AbstractMatrix{Tv}
    m::Int
    n::Int
    N::Int
    i::Ti
    j::Ti
    pos::Vector{Ti}
    idx::Vector{Ti}
    val::Vector{Tv}
    Δ::Vector{Ti}
    c::Tv
end

Base.size(arg::DominanceSum) = (arg.m + 1, arg.n + 1)

dominancesum(args...; kwargs...) = dominancesum(NoHint(), args...; kwargs...)
dominancesum(::AbstractHint, args...; kwargs...) = @assert false
function dominancesum(hint::AbstractHint, A::SparseMatrixCSC{Tv, Ti}; kwargs...) where {Tv, Ti}
    (m, n) = size(A)
    N = nnz(A)
    pos = copy(A.colptr)
    idx = copy(A.rowval)
    val = A.nzval
    return dominancesum!(hint, m, n, N, pos, idx, val; kwargs...)
end

dominancesum!(args...; kwargs...) = dominancesum!(NoHint(), args...; kwargs...)
dominancesum!(::AbstractHint, args...; kwargs...) = @assert false
function dominancesum!(hint::AbstractHint, m, n, N, pos, idx, val; H = nothing, b = nothing, b′ = nothing, kwargs...)
    if H === b === b′ === nothing
        return DominanceSum(hint, m, n, N, pos, idx, val; b = 4, kwargs...)
    else
        return DominanceSum(hint, m, n, N, pos, idx, val; H = H, b = b, b′ = b′, kwargs...)
    end
end
function dominancesum!(hint::SparseHint, m, n, N, pos, idx, val; kwargs...)
    return DominanceSum(hint, m, n, N, pos, idx, val; kwargs...)
end
function dominancesum!(hint::StepHint, m, n, N, pos, idx, val; kwargs...)
    SparseStepwiseDominanceCount(hint, m, n, N, pos, idx, val; kwargs...)
end

DominanceSum(hint::AbstractHint, m, n, N, pos::Vector{Ti}, idx::Vector{Ti}, val::Vector{Tv}; kwargs...) where {Tv, Ti} =
    DominanceSum{Tv, Ti}(hint, m, n, N, pos, idx, val; kwargs...)
function DominanceSum{Tv, Ti}(hint::AbstractHint, m, n, N, pos::Vector{Ti}, idx::Vector{Ti}, val::Vector{Tv}; b = nothing, H = nothing, b′ = nothing) where {Tv, Ti}
    @inbounds begin
        if b === nothing #b = branching factor of tree
            if H === nothing
                b = cld(cllog2(m + 1), 3)
            else
                b = cld(cllog2(m + 1), H)
            end
        end

        if H === nothing #H = height of tree
            H = cld(cllog2(m + 1), b)
        end

        if b′ === nothing #b′ = cache period
            b′ = b + cllog2(H)
        end

        qos = zeros(Ti, m + 2)
        qos[1] = 1
        qos[end] = N + 1
        bkt = undefs(Ti, 1 << b + 1) #bkt = scratch space for cumulative count
        cnt = undefs(Ti, 1 << b, N >> b′ + 1, H - 1) #cnt = cached cumulative counts
        wgt = undefs(Tv, N, H) #wgt = permuted copies of weights
        for q = 1:N
            wgt[q, end] = val[q]
        end
        pre = undefs(Tv, 1 << b + 1) #pre = scratch space for cumulative sum
        scn = undefs(Tv, 1 << b + 1, N >> b′ + 1, H) #scn = cached cumulative sums

        byt = undefs(Ti, N)

        for h = H : -1 : 2 #h = level in tree
            for i′ = 1 : 1 << (h * b) : m + 1 #i′ = "quotient" of sorts
                zero!(bkt)
                for q = qos[i′] : qos[min(i′ + (1 << (h * b)), end)] - 1 #q = position in current level
                    i = idx[q]
                    d = i >> ((h - 1) * b) & ((1 << b) - 1) + 1 #d = a "remainder" of sorts
                    bkt[d + 1] += 1
                end
                bkt[1] = qos[i′]
                for d = 1:(1 << b)
                    bkt[d + 1] = bkt[d] + bkt[d + 1]
                end
                for q = qos[i′] : qos[min(i′ + (1 << (h * b)), end)] - 1
                    i = idx[q]
                    d = i >> ((h - 1) * b) & (1 << b - 1) + 1
                    q′ = bkt[d]
                    byt[q′] = (idx[q′] & ~(1 << ((h - 1) * b) - 1)) | (i & (1 << ((h - 1) * b) - 1))
                    wgt[q′, h - 1] = wgt[q, h]
                    bkt[d] = q′ + 1
                end
                for d = 1 : 1 << b
                    qos[min(i′ + d << ((h - 1) * b), end)] = bkt[d]
                end
            end
            for d = 1 : 1 << b
                cnt[d, 1, h - 1] = 0
            end
            scn[1, 1, h] = 0
            for d = 1 : 1 << b
                scn[d + 1, 1, h] = 0
            end
            zero!(bkt)
            zero!(pre)
            for q = 1 : N
                i = idx[q]
                d = i >> ((h - 1) * b) & ((1 << b) - 1) + 1 #d = a "remainder" of sorts
                bkt[d] += 1
                pre[d] += wgt[q, h]
                if q & (1 << b′ - 1) == 0 #Cache bkt at the end of each 1 << b′ block
                    Q = (q >> b′) + 1
                    scn[1, Q, h] = 0
                    for d = 1 : 1 << b
                        cnt[d, Q, h - 1] = bkt[d]
                        scn[d + 1, Q, h] = pre[d] + scn[d, Q, h]
                    end
                end
            end

            idx, byt = byt, idx
        end

        for i′ = 1 : 1 << b : m + 1 #i′ = "quotient" of sorts
            zero!(bkt)
            for q = qos[i′] : qos[min(i′ + (1 << b), end)] - 1 #q = position in current level
                i = idx[q]
                d = i & ((1 << b) - 1) + 1 #d = a "remainder" of sorts
                bkt[d + 1] += 1
            end
            bkt[1] = qos[i′]
            for d = 1:(1 << b)
                bkt[d + 1] = bkt[d] + bkt[d + 1]
            end
            for q = qos[i′] : qos[min(i′ + (1 << b), end)] - 1
                i = idx[q]
                d = i & (1 << b - 1) + 1
                bkt[d] += 1
            end
            for d = 1 : 1 << b
                qos[min(i′ + d, end)] = bkt[d]
            end
        end
        scn[1, 1, 1] = 0
        for d = 1 : 1 << b
            scn[d + 1, 1, 1] = 0
        end
        zero!(pre)
        for q = 1 : N
            i = idx[q]
            d = i & ((1 << b) - 1) + 1 #d = a "remainder" of sorts
            pre[d] += wgt[q, 1]
            if q & (1 << b′ - 1) == 0 #Cache bkt at the end of each 1 << b′ block
                Q = (q >> b′) + 1
                scn[1, Q, 1] = 0
                for d = 1 : 1 << b
                    scn[d + 1, Q, 1] = pre[d] + scn[d, Q, 1]
                end
            end
        end

        byt = idx

        return DominanceSum{Tv, Ti}(m, n, N, b, b′, H, pos, qos, byt, cnt, wgt, scn)
    end
end

Base.getindex(arg::DominanceSum{Tv, Ti}, i::Integer, j::Integer) where {Tv, Ti} = arg(i, j)
function (arg::DominanceSum{Tv, Ti})(i::Integer, j::Integer) where {Tv, Ti}
    @inbounds begin
        b = arg.b
        b′ = arg.b′
        H = arg.H
        pos = arg.pos
        qos = arg.qos
        byt = arg.byt
        cnt = arg.cnt
        wgt = arg.wgt
        scn = arg.scn

        Δq = pos[j] - 1
        i = Ti(i - 1)::Ti
        s = Tv(0)
        for h = H : -1 : 2
            i′ = i & ~(1 << (h * b) - 1) + 1
            q₁ = qos[i′] - 1
            q₂ = q₁ + Δq
            d = i >> ((h - 1) * b) & ((1 << b) - 1) + 1
            Q₁ = (q₁ >> b′) + 1
            Q₂ = (q₂ >> b′) + 1
            s += scn[d, Q₂, h] - scn[d, Q₁, h]
            Δq = cnt[d, Q₂, h - 1] - cnt[d, Q₁, h - 1]
            for q = (Q₁ - 1) << b′ + 1 : q₁
                d′ = byt[q] >> ((h - 1) * b) & ((1 << b) - 1) + 1
                if d′ < d
                    s -= wgt[q, h]
                end
                Δq -= d′ == d
            end
            for q = (Q₂ - 1) << b′ + 1 : q₂
                d′ = byt[q] >> ((h - 1) * b) & ((1 << b) - 1) + 1
                if d′ < d
                    s += wgt[q, h]
                end
                Δq += d′ == d
            end
        end
        i′ = i & ~(1 << b - 1) + 1
        q₁ = qos[i′] - 1
        q₂ = q₁ + Δq

        d = i & ((1 << b) - 1) + 1
        Q₁ = (q₁ >> b′) + 1
        Q₂ = (q₂ >> b′) + 1
        s += scn[d + 1, Q₂, 1] - scn[d + 1, Q₁, 1]
        for q = (Q₁ - 1) << b′ + 1 : q₁
            d′ = byt[q] & ((1 << b) - 1) + 1
            if d′ <= d
                s -= wgt[q, 1]
            end
        end
        for q = (Q₂ - 1) << b′ + 1 : q₂
            d′ = byt[q] & ((1 << b) - 1) + 1
            if d′ <= d
                s += wgt[q, 1]
            end
        end

        return s
    end
end

SparseStepwiseDominanceSum(hint::AbstractHint, m, n, N, pos::Vector{Ti}, idx::Vector{Ti}, val::Vector{Tv}; kwargs...) where {Tv, Ti} =
    SparseStepwiseDominanceSum{Tv, Ti}(hint, m, n, N, pos, idx, val; kwargs...)
function SparseStepwiseDominanceSum{Tv, Ti}(hint::AbstractHint, m, n, N, pos::Vector{Ti}, idx::Vector{Ti}, val::Vector{Tv}; kwargs...) where {Tv, Ti}
    @inbounds begin
        i = j = 0
        Δ = zeros(Ti, m)
        c = Tv(0)
        return SparseStepwiseDominanceSum(m, n, N, i, j, pos, idx, val, Δ, c)
    end
end

Base.getindex(arg::SparseStepwiseDominanceSum{Tv, Ti}, i::Integer, j::Integer) where {Tv, Ti} = arg(i, j)
function (arg::SparseStepwiseDominanceSum{Tv, Ti})(i::Integer, j::Integer) where {Tv, Ti}
    @inbounds begin
        i -= 1
        j -= 1
        c = arg.c
        Δ = arg.Δ
        pos = arg.pos
        idx = arg.idx
        val = arg.val
        arg_i = arg.i
        arg_j = arg.j
        #reset case
        if (arg.m + pos[j + 1]) < pos[arg_j + 1] - pos[j + 1]
            arg_j = zero(Ti)
            c = zero(Tv)
            zero!(Δ)
        end
        for q = pos[j + 1] : pos[arg_j + 1] - 1
            Δ[idx[q]] -= val[q]
            if idx[q] <= i
                c -= val[q]
            end
        end
        for q = pos[arg_j + 1] : pos[j + 1] - 1
            Δ[idx[q]] += val[q]
            if idx[q] <= i
                c += val[q]
            end
        end
        for q = i + 1 : arg_i
            c -= Δ[q]
        end
        for q = arg_i + 1 : i
            c += Δ[q]
        end
        arg.i = i
        arg.j = j
        arg.c = c
        return c
    end
end

@propagate_inbounds function (stp::Step{Cnt})(i::Same, j::Same) where {Tv, Ti, Cnt <: SparseStepwiseDominanceSum{Tv, Ti}}
    begin
        arg = stp.ocl
        return arg.c
    end
end

@propagate_inbounds function (stp::Step{Cnt})(_i::Same, _j::Next) where {Tv, Ti, Cnt <: SparseStepwiseDominanceSum{Tv, Ti}}
    begin
        arg = stp.ocl
        i = destep(_i)
        j = destep(_j)
        i -= 1
        j -= 1
        c = arg.c
        Δ = arg.Δ
        pos = arg.pos
        idx = arg.idx
        for q = pos[j] : pos[j + 1] - 1
            Δ[idx[q]] += val[q]
            if idx[q] <= i
                c += val[q]
            end
        end
        arg.j = j
        arg.c = c
        return c
    end
end

@propagate_inbounds function (stp::Step{Cnt})(_i::Same, _j::Prev) where {Tv, Ti, Cnt <: SparseStepwiseDominanceSum{Tv, Ti}}
    begin
        arg = stp.ocl
        i = destep(_i)
        j = destep(_j)
        i -= 1
        j -= 1
        c = arg.c
        Δ = arg.Δ
        pos = arg.pos
        idx = arg.idx
        val = arg.val
        for q = pos[j + 1] : pos[j + 2] - 1
            Δ[idx[q]] -= val[q]
            if idx[q] <= i
                c -= val[q]
            end
        end
        arg.j = j
        arg.c = c
        return c
    end
end

@propagate_inbounds function (stp::Step{Cnt})(_i::Next, _j::Same) where {Tv, Ti, Cnt <: SparseStepwiseDominanceSum{Tv, Ti}}
    begin
        arg = stp.ocl
        i = destep(_i)
        j = destep(_j)
        i -= 1
        j -= 1
        c = arg.c
        Δ = arg.Δ
        c += Δ[i]
        arg.i = i
        arg.c = c
        return c
    end
end

@propagate_inbounds function (stp::Step{Cnt})(_i::Prev, _j::Same) where {Tv, Ti, Cnt <: SparseStepwiseDominanceSum{Tv, Ti}}
    begin
        arg = stp.ocl
        i = destep(_i)
        j = destep(_j)
        i -= 1
        j -= 1
        c = arg.c
        Δ = arg.Δ
        c -= Δ[i + 1]
        arg.i = i
        arg.c = c
        return c
    end
end

struct DominanceCount{Ti} <: AbstractMatrix{Ti}
    m::Int
    n::Int
    N::Int
    b::Int
    b′::Int
    H::Int
    pos::Vector{Ti}
    qos::Vector{Ti}
    byt::Vector{Ti}
    cnt::Array{Ti, 3}
end

Base.size(arg::DominanceCount) = (arg.m + 1, arg.n + 1)

struct BinaryDominanceCount{Ti, Tb} <: AbstractMatrix{Ti}
    m::Int
    n::Int
    N::Int
    H::Int
    pos::Vector{Ti}
    qos::Vector{Ti}
    byt::Array{Tb, 2}
    cnt::Array{Ti, 2}
end

Base.size(arg::BinaryDominanceCount) = (arg.m + 1, arg.n + 1)

mutable struct SparseStepwiseDominanceCount{Ti} <: AbstractMatrix{Ti}
    m::Int
    n::Int
    N::Int
    i::Ti
    j::Ti
    pos::Vector{Ti}
    idx::Vector{Ti}
    Δ::Vector{Ti}
    c::Ti
end

Base.size(arg::SparseStepwiseDominanceCount) = (arg.m + 1, arg.n + 1)

dominancecount(args...; kwargs...) = dominancecount(NoHint(), args...; kwargs...)
dominancecount(::AbstractHint, args...; kwargs...) = @assert false
function dominancecount(hint::AbstractHint, A::SparseMatrixCSC{Tv, Ti}; kwargs...) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)

        return dominancecount!(hint, m, n, nnz(A), copy(A.colptr), copy(A.rowval); kwargs...)
    end
end

dominancecount!(args...; kwargs...) = dominancecount!(NoHint(), args...; kwargs...)
dominancecount!(::AbstractHint, args...; kwargs...) = @assert false
function dominancecount!(hint::AbstractHint, m, n, N, pos::Vector{Ti}, idx::Vector{Ti}; kwargs...) where {Ti}
    BinaryDominanceCount(hint, m, n, N, pos, idx; kwargs...)
end
function dominancecount!(hint::SparseHint, m, n, N, pos::Vector{Ti}, idx::Vector{Ti}; kwargs...) where {Ti}
    DominanceCount(hint, m, n, N, pos, idx; kwargs...)
end
function dominancecount!(hint::StepHint, m, n, N, pos::Vector{Ti}, idx::Vector{Ti}; kwargs...) where {Ti}
    SparseStepwiseDominanceCount(hint, m, n, N, pos, idx; kwargs...)
end

DominanceCount(hint::AbstractHint, m, n, N, pos::Vector{Ti}, idx::Vector{Ti}; kwargs...) where {Ti} =
    DominanceCount{Ti}(hint, m, n, N, pos, idx; kwargs...)
function DominanceCount{Ti}(hint::AbstractHint, m, n, N, pos::Vector{Ti}, idx::Vector{Ti}; b = nothing, H = nothing, b′ = nothing, kwargs...) where {Ti}
    @inbounds begin
        if b === nothing #b = branching factor of tree
            if H === nothing
                b = cld(cllog2(m + 1), 3)
            else
                b = cld(cllog2(m + 1), H)
            end
        end

        if H === nothing #H = height of tree
            H = cld(cllog2(m + 1), b)
        end

        if b′ === nothing #b′ = cache period
            b′ = b + cllog2(H)
        end

        qos = zeros(Ti, m + 2)
        qos[1] = 1
        qos[end] = N + 1
        bkt = undefs(Ti, 1 << b + 1) #bkt = scratch space for cumulative count
        cnt = undefs(Ti, 1 << b + 1, N >> b′ + 1, H) #cnt = cached cumulative counts

        byt = undefs(Ti, N)

        for h = H : -1 : 1 #h = level in tree
            for i′ = 1 : 1 << (h * b) : m + 1 #i′ = "quotient" of sorts
                zero!(bkt)
                for q = qos[i′] : qos[min(i′ + (1 << (h * b)), end)] - 1 #q = position in current level
                    i = idx[q]
                    d = i >> ((h - 1) * b) & ((1 << b) - 1) + 1 #d = a "remainder" of sorts
                    bkt[d + 1] += 1
                end
                bkt[1] = qos[i′]
                for d = 1:(1 << b)
                    bkt[d + 1] = bkt[d] + bkt[d + 1]
                end
                for q = qos[i′] : qos[min(i′ + (1 << (h * b)), end)] - 1
                    i = idx[q]
                    d = i >> ((h - 1) * b) & (1 << b - 1) + 1
                    q′ = bkt[d]
                    byt[q′] = (idx[q′] & ~(1 << ((h - 1) * b) - 1)) | (i & (1 << ((h - 1) * b) - 1))
                    bkt[d] = q′ + 1
                end
                for d = 1 : 1 << b
                    qos[min(i′ + d << ((h - 1) * b), end)] = bkt[d]
                end
            end
            cnt[1, 1, h] = 0
            for d = 1 : 1 << b
                cnt[d + 1, 1, h] = 0
            end
            zero!(bkt)
            for q = 1 : N
                i = idx[q]
                d = i >> ((h - 1) * b) & ((1 << b) - 1) + 1 #d = a "remainder" of sorts
                bkt[d] += 1
                if q & (1 << b′ - 1) == 0 #Cache bkt at the end of each 1 << b′ block
                    Q = (q >> b′) + 1
                    cnt[1, Q, h] = 0
                    for d = 1 : 1 << b
                        cnt[d + 1, Q, h] = bkt[d] + cnt[d, Q, h]
                    end
                end
            end

            idx, byt = byt, idx
        end
        byt = idx
        return DominanceCount{Ti}(m, n, N, b, b′, H, pos, qos, byt, cnt)
    end
end

Base.getindex(arg::DominanceCount{Ti}, i::Integer, j::Integer) where {Ti} = arg(i, j)
function (arg::DominanceCount{Ti})(i::Integer, j::Integer) where {Ti}
    @inbounds begin
        b = arg.b
        b′ = arg.b′
        H = arg.H
        pos = arg.pos
        qos = arg.qos
        byt = arg.byt
        cnt = arg.cnt

        Δq = pos[j] - 1
        i = Ti(i - 1)::Ti
        s = Ti(0)
        for h = H : -1 : 2
            i′ = i & ~(1 << (h * b) - 1) + 1
            q₁ = qos[i′] - 1
            q₂ = q₁ + Δq
            d = i >> ((h - 1) * b) & ((1 << b) - 1) + 1
            Q₁ = (q₁ >> b′) + 1
            Q₂ = (q₂ >> b′) + 1
            s += cnt[d, Q₂, h] - cnt[d, Q₁, h]
            Δq = (cnt[d + 1, Q₂, h] - cnt[d, Q₂, h]) - (cnt[d + 1, Q₁, h] - cnt[d, Q₁, h])
            msk = ((1 << b) - 1) << ((h - 1) * b)
            cmp = (d - 1) << ((h - 1) * b)
            for q = (Q₁ - 1) << b′ + 1 : q₁
                d′ = byt[q] & msk
                s -= d′ < cmp
                Δq -= d′ == cmp
                #d′ = byt[q] >> ((h - 1) * b) & ((1 << b) - 1) + 1
                #s -= d′ < d
                #Δq -= d′ == d
            end
            for q = (Q₂ - 1) << b′ + 1 : q₂
                d′ = byt[q] & msk
                s += d′ < cmp
                Δq += d′ == cmp
                #d′ = byt[q] >> ((h - 1) * b) & ((1 << b) - 1) + 1
                #s += d′ < d
                #Δq += d′ == d
            end
        end

        i′ = i & ~(1 << b - 1) + 1
        q₁ = qos[i′] - 1
        q₂ = q₁ + Δq

        d = i & ((1 << b) - 1) + 1
        Q₁ = (q₁ >> b′) + 1
        Q₂ = (q₂ >> b′) + 1
        s += cnt[d + 1, Q₂, 1] - cnt[d + 1, Q₁, 1]
        msk = (1 << b) - 1
        cmp = d - 1
        for q = (Q₁ - 1) << b′ + 1 : q₁
            d′ = byt[q] & msk
            s -= d′ <= cmp
            #d′ = byt[q] & ((1 << b) - 1) + 1
            #s -= d′ <= d
        end
        for q = (Q₂ - 1) << b′ + 1 : q₂
            d′ = byt[q] & msk
            s += d′ <= cmp
            #d′ = byt[q] & ((1 << b) - 1) + 1
            #s += d′ <= d
        end

        return s
    end
end

BinaryDominanceCount(hint::AbstractHint, m, n, N, pos, idx; kwargs...) =
    BinaryDominanceCount{UInt}(hint, m, n, N, pos, idx; kwargs...)
BinaryDominanceCount{Tb}(hint::AbstractHint, m, n, N, pos::Vector{Ti}, idx::Vector{Ti}; kwargs...) where {Tb, Ti} =
    BinaryDominanceCount{Tb, Ti}(hint, m, n, N, pos, idx; kwargs...)
function BinaryDominanceCount{Tb, Ti}(hint::AbstractHint, m, n, N, pos::Vector{Ti}, idx::Vector{Ti}; kwargs...) where {Tb, Ti}
    @inbounds begin
        #b = branching factor of tree = 1

        #H = height of tree
        H = cllog2(m + 1)

        qos = zeros(Ti, m + 2)
        qos[1] = 1
        qos[end] = N + 1
        cnt = zeros(Ti, 1 + cld(N, nbits(Tb)), H) #cnt = cached cumulative counts

        byt = zeros(Tb, 1 + cld(N, nbits(Tb)), H)
        idx′ = undefs(Ti, N)

        for h = H : -1 : 1 #h = level in tree
            _cnt = 0
            for i′ = 1 : 1 << h : m + 1 #i′ = "quotient" of sorts
                bkt_1 = 0
                bkt_2 = 0
                for q = qos[i′] : qos[min(i′ + (1 << h), end)] - 1 #q = position in current level
                    i = idx[q]
                    d = i >> (h - 1) & 1 #d = a "remainder" of sorts
                    Q = ((q - 1) >> log2nbits(Tb)) + 1
                    byt[Q, h] |= d << ((q - 1) & (nbits(Tb) - 1))
                    cnt[Q + 1, h] = _cnt += d
                    bkt_2 += 1 - d
                end
                bkt_1 = qos[i′]
                bkt_2 += bkt_1
                for q = qos[i′] : qos[min(i′ + (1 << h), end)] - 1
                    i = idx[q]
                    d = i >> (h - 1) & 1
                    q′ = ifelse(d == 0, bkt_1, bkt_2)
                    idx′[q′] = i
                    bkt_1 += 1 - d
                    bkt_2 += d
                end
                qos[min(i′ + 1 << (h - 1), end)] = bkt_1
            end

            idx, idx′ = idx′, idx
        end
        return BinaryDominanceCount{Ti, Tb}(m, n, N, H, pos, qos, byt, cnt)
    end
end

#we have removed b′ because increasing it doesn't significantly reduce the storage cost or preprocessing cost of the structure.

Base.getindex(arg::BinaryDominanceCount{Ti, Tb}, i::Integer, j::Integer) where {Ti, Tb} = arg(i, j)
function (arg::BinaryDominanceCount{Ti, Tb})(i::Integer, j::Integer) where {Ti, Tb}
    @inbounds begin
        H = arg.H
        pos = arg.pos
        qos = arg.qos
        byt = arg.byt
        cnt = arg.cnt

        Δq = pos[j] - 1
        i = Ti(i - 1)::Ti
        s = Ti(0)
        for h = H : -1 : 1
            i′ = i & ~(1 << h - 1) + 1
            q₁ = qos[i′] - 1
            q₂ = q₁ + Δq
            d = i >> (h - 1) & Ti(1)

            Q₁ = (q₁ >> log2nbits(Tb)) + 1
            Q₂ = (q₂ >> log2nbits(Tb)) + 1
            bkt_2 = cnt[Q₂, h] - cnt[Q₁, h]
            bkt_2 += count_ones(byt[Q₂, h] & ((one(Tb) << (q₂ & (nbits(Tb) - 1))) - 1))
            bkt_2 -= count_ones(byt[Q₁, h] & ((one(Tb) << (q₁ & (nbits(Tb) - 1))) - 1))
            bkt_1 = Δq - bkt_2
            s += ifelse(d == 0, 0, bkt_1)
            Δq = ifelse(d == 0, bkt_1, bkt_2)
        end

        return s + Δq
    end
end

SparseStepwiseDominanceCount(hint::AbstractHint, m, n, N, pos::Vector{Ti}, idx::Vector{Ti}; kwargs...) where {Ti} =
    SparseStepwiseDominanceCount{Ti}(hint, m, n, N, pos, idx; kwargs...)
function SparseStepwiseDominanceCount{Ti}(hint::AbstractHint, m, n, N, pos::Vector{Ti}, idx::Vector{Ti}; kwargs...) where {Ti}
    @inbounds begin
        i = j = 0
        Δ = zeros(Ti, m)
        c = Ti(0)
        return SparseStepwiseDominanceCount(m, n, N, i, j, pos, idx, Δ, c)
    end
end

Base.getindex(arg::SparseStepwiseDominanceCount{Ti}, i::Integer, j::Integer) where {Ti} = arg(i, j)
function (arg::SparseStepwiseDominanceCount{Ti})(i::Integer, j::Integer) where {Ti}
    @inbounds begin
        i -= 1
        j -= 1
        c = arg.c
        Δ = arg.Δ
        pos = arg.pos
        idx = arg.idx
        arg_i = arg.i
        arg_j = arg.j
        #reset case
        if (arg.m + pos[j + 1]) < pos[arg_j + 1] - pos[j + 1]
            arg_j = zero(Ti)
            c = zero(Ti)
            zero!(Δ)
        end
        for q = pos[j + 1] : pos[arg_j + 1] - 1
            Δ[idx[q]] -= 1
            c -= idx[q] <= arg_i
        end
        for q = pos[arg_j + 1] : pos[j + 1] - 1
            Δ[idx[q]] += 1
            c += idx[q] <= arg_i
        end
        for q = i + 1 : arg_i
            c -= Δ[q]
        end
        for q = arg_i + 1 : i
            c += Δ[q]
        end
        arg.i = i
        arg.j = j
        arg.c = c
        return c
    end
end

@propagate_inbounds function (stp::Step{Cnt})(i::Same, j::Same) where {Ti, Cnt <: SparseStepwiseDominanceCount{Ti}}
    begin
        arg = stp.ocl
        return arg.c
    end
end

@propagate_inbounds function (stp::Step{Cnt})(_i::Same, _j::Next) where {Ti, Cnt <: SparseStepwiseDominanceCount{Ti}}
    begin
        arg = stp.ocl
        i = destep(_i)
        j = destep(_j)
        i -= 1
        j -= 1
        c = arg.c
        Δ = arg.Δ
        pos = arg.pos
        idx = arg.idx
        for q = pos[j] : pos[j + 1] - 1
            Δ[idx[q]] += 1
            c += idx[q] <= i
        end
        arg.j = j
        arg.c = c
        return c
    end
end

@propagate_inbounds function (stp::Step{Cnt})(_i::Same, _j::Prev) where {Ti, Cnt <: SparseStepwiseDominanceCount{Ti}}
    begin
        arg = stp.ocl
        i = destep(_i)
        j = destep(_j)
        i -= 1
        j -= 1
        c = arg.c
        Δ = arg.Δ
        pos = arg.pos
        idx = arg.idx
        for q = pos[j + 1] : pos[j + 2] - 1
            Δ[idx[q]] -= 1
            c -= idx[q] <= i
        end
        arg.j = j
        arg.c = c
        return c
    end
end

@propagate_inbounds function (stp::Step{Cnt})(_i::Next, _j::Same) where {Ti, Cnt <: SparseStepwiseDominanceCount{Ti}}
    begin
        arg = stp.ocl
        i = destep(_i)
        j = destep(_j)
        i -= 1
        j -= 1
        c = arg.c
        Δ = arg.Δ
        c += Δ[i]
        arg.i = i
        arg.c = c
        return c
    end
end

@propagate_inbounds function (stp::Step{Cnt})(_i::Prev, _j::Same) where {Ti, Cnt <: SparseStepwiseDominanceCount{Ti}}
    begin
        arg = stp.ocl
        i = destep(_i)
        j = destep(_j)
        i -= 1
        j -= 1
        c = arg.c
        Δ = arg.Δ
        c -= Δ[i + 1]
        arg.i = i
        arg.c = c
        return c
    end
end

struct RookSum{Ti, Tv} <: AbstractMatrix{Tv}
    N::Int
    b::Int
    b′::Int
    H::Int
    byt::Vector{Ti}
    cnt::Array{Ti, 3}
    wgt::Array{Tv, 2}
    scn::Array{Tv, 3}
end

Base.size(arg::RookSum) = (arg.N + 1, arg.N + 1)

rooksum!(args...; kwargs...) = rooksum!(NoHint(), args...; kwargs...)
rooksum!(::AbstractHint, args...; kwargs...) = @assert false
function rooksum!(hint::AbstractHint, N, idx, val; H = nothing, b = nothing, b′ = nothing, kwargs...)
    if H === b === b′ === nothing
        return RookSum(hint, N, idx, val; b = 4, kwargs...)
    else
        return RookSum(hint, N, idx, val; H = H, b = b, b′ = b′, kwargs...)
    end
end
function rooksum!(hint::SparseHint, N, idx, val; kwargs...)
    return RookSum(hint, N, idx, val; kwargs...)
end

RookSum(hint::AbstractHint, N, idx::Vector{Ti}, val::Vector{Tv}; kwargs...) where {Ti, Tv} =
    RookSum{Ti, Tv}(hint, N, idx, val; kwargs...)
function RookSum{Ti, Tv}(hint::AbstractHint, N, idx::Vector{Ti}, val::Vector{Tv}; b = nothing, H = nothing, b′ = nothing, kwargs...) where {Ti, Tv}
    @inbounds begin
        @debug begin
            @assert length(idx) == N
            hash = falses(N)
            for i in idx
                @assert 1 <= i <= N
                @assert !hash[i]
                hash[i] = true
            end
        end

        if b === nothing #b = branching factor of tree
            if H === nothing
                b = cld(cllog2(N + 1), 3)
            else
                b = cld(cllog2(N + 1), H)
            end
        end

        if H === nothing #H = height of tree
            H = cld(cllog2(N + 1), b)
        end

        if b′ === nothing #b′ = cache period
            b′ = b + cllog2(H)
        end

        bkt = undefs(Ti, 1 << b) #bkt = scratch space for cumulative count
        cnt = undefs(Ti, 1 << b, N >> b′ + 1, H) #cnt = cached cumulative counts
        pre = undefs(Tv, 1 << b + 1) #pre = scratch space for cumulative sum
        scn = undefs(Tv, 1 << b + 1, N >> b′ + 1, H) #scn = cached cumulative sums
        wgt = undefs(Tv, N, H) #wgt = permuted copies of weights
        for q = 1:N
            wgt[q, end] = val[q]
        end

        byt = undefs(Ti, N)

        for h = H : -1 : 2 #h = level in tree
            for i′ = 1 : 1 << (h * b) : N + 1 #i′ = "quotient" of sorts
                for d = 1:(1 << b)
                    bkt[d] = max(i′ + (d - 1) << ((h - 1) * b) - 1, 1)
                end
                for q = max(i′ - 1, 1) : min(i′ + (1 << (h * b)) - 1, N + 1) - 1
                    i = idx[q]
                    d = i >> ((h - 1) * b) & (1 << b - 1) + 1
                    q′ = bkt[d]
                    byt[q′] = (idx[q′] & ~(1 << ((h - 1) * b) - 1)) | (i & (1 << ((h - 1) * b) - 1))
                    wgt[q′, h - 1] = wgt[q, h]
                    bkt[d] = q′ + 1
                end
            end

            for d = 1 : 1 << b
                cnt[d, 1, h - 1] = 0
            end
            scn[1, 1, h] = 0
            for d = 1 : 1 << b
                scn[d + 1, 1, h] = 0
            end
            zero!(bkt)
            zero!(pre)
            for q = 1 : N
                i = idx[q]
                d = i >> ((h - 1) * b) & ((1 << b) - 1) + 1 #d = a "remainder" of sorts
                bkt[d] += 1
                pre[d] += wgt[q, h]
                if q & (1 << b′ - 1) == 0 #Cache bkt at the end of each 1 << b′ block
                    Q = (q >> b′) + 1
                    scn[1, Q, h] = 0
                    for d = 1 : 1 << b
                        cnt[d, Q, h - 1] = bkt[d]
                        scn[d + 1, Q, h] = pre[d] + scn[d, Q, h]
                    end
                end
            end
            idx, byt = byt, idx
        end

        scn[1, 1, 1] = 0
        for d = 1 : 1 << b
            scn[d + 1, 1, 1] = 0
        end
        zero!(pre)
        for q = 1 : N
            i = idx[q]
            d = i & ((1 << b) - 1) + 1 #d = a "remainder" of sorts
            pre[d] += wgt[q, 1]
            if q & (1 << b′ - 1) == 0 #Cache bkt at the end of each 1 << b′ block
                Q = (q >> b′) + 1
                scn[1, Q, 1] = 0
                for d = 1 : 1 << b
                    scn[d + 1, Q, 1] = pre[d] + scn[d, Q, 1]
                end
            end
        end

        byt = idx
        return RookSum{Ti, Tv}(N, b, b′, H, byt, cnt, wgt, scn)
    end
end

Base.getindex(arg::RookSum{Ti, Tv}, i::Integer, j::Integer) where {Ti, Tv} = arg(i, j)
function (arg::RookSum{Ti, Tv})(i::Integer, j::Integer) where {Ti, Tv}
    @inbounds begin
        N = arg.N
        b = arg.b
        b′ = arg.b′
        H = arg.H
        byt = arg.byt
        cnt = arg.cnt
        wgt = arg.wgt
        scn = arg.scn

        Δq = j - 1
        i = Ti(i - 1)::Ti
        s = Tv(0)
        for h = H : -1 : 2
            i′ = i & ~(1 << (h * b) - 1) + 1
            q₁ = max(i′ - 1, 1) - 1
            q₂ = q₁ + Δq
            d = i >> ((h - 1) * b) & ((1 << b) - 1) + 1
            Q₁ = (q₁ >> b′) + 1
            Q₂ = (q₂ >> b′) + 1
            s += scn[d, Q₂, h] - scn[d, Q₁, h]
            Δq = cnt[d, Q₂, h - 1] - cnt[d, Q₁, h - 1]
            for q = (Q₁ - 1) << b′ + 1 : q₁
                d′ = byt[q] >> ((h - 1) * b) & ((1 << b) - 1) + 1
                if d′ < d
                    s -= wgt[q, h]
                end
                Δq -= d′ == d
            end
            for q = (Q₂ - 1) << b′ + 1 : q₂
                d′ = byt[q] >> ((h - 1) * b) & ((1 << b) - 1) + 1
                if d′ < d
                    s += wgt[q, h]
                end
                Δq += d′ == d
            end
        end

        i′ = i & ~(1 << b - 1) + 1
        q₁ = max(i′ - 1, 1) - 1
        q₂ = q₁ + Δq

        d = i & ((1 << b) - 1) + 1
        Q₁ = (q₁ >> b′) + 1
        Q₂ = (q₂ >> b′) + 1
        s += scn[d + 1, Q₂, 1] - scn[d + 1, Q₁, 1]
        for q = (Q₁ - 1) << b′ + 1 : q₁
            d′ = byt[q] & ((1 << b) - 1) + 1
            if d′ <= d
                s -= wgt[q, 1]
            end
        end
        for q = (Q₂ - 1) << b′ + 1 : q₂
            d′ = byt[q] & ((1 << b) - 1) + 1
            if d′ <= d
                s += wgt[q, 1]
            end
        end

        return s
    end
end

struct RookCount{Ti} <: AbstractMatrix{Ti}
    N::Int
    b::Int
    b′::Int
    H::Int
    byt::Vector{Ti}
    cnt::Array{Ti, 3}
end

Base.size(arg::RookCount) = (arg.N + 1, arg.N + 1)

struct BinaryRookCount{Tb, Ti} <: AbstractMatrix{Ti}
    N::Int
    H::Int
    byt::Array{Tb, 2}
    cnt::Array{Ti, 2}
end

Base.size(arg::BinaryRookCount) = (arg.N + 1, arg.N + 1)

mutable struct SparseRookStepCounter{Ti} <: AbstractMatrix{Ti}
    N::Int
    i::Ti
    j::Ti
    idx::Vector{Ti}
    Δ::Vector{Ti}
    c::Ti
end

Base.size(arg::SparseRookStepCounter) = (arg.N + 1, arg.N + 1)

rookcount!(args...; kwargs...) = rookcount!(NoHint(), args...; kwargs...)
rookcount!(::AbstractHint, args...; kwargs...) = @assert false
function rookcount!(hint::AbstractHint, N, idx; kwargs...)
    return BinaryRookCount(hint, N, idx; kwargs...)
end
function rookcount!(hint::SparseHint, N, idx; kwargs...)
    RookCount(hint, N, idx; kwargs...)
end

RookCount(hint::AbstractHint, N, idx::Vector{Ti}; kwargs...) where {Ti} =
    RookCount{Ti}(hint, N, idx; kwargs...)
function RookCount{Ti}(hint::AbstractHint, N, idx::Vector{Ti}; b = nothing, H = nothing, b′ = nothing, kwargs...) where {Ti}
    @inbounds begin
        @debug begin
            @assert length(idx) == N
            hash = falses(N)
            for i in idx
                @assert 1 <= i <= N
                @assert !hash[i]
                hash[i] = true
            end
        end

        if b === nothing #b = branching factor of tree
            if H === nothing
                b = cld(cllog2(N + 1), 3)
            else
                b = cld(cllog2(N + 1), H)
            end
        end

        if H === nothing #H = height of tree
            H = cld(cllog2(N + 1), b)
        end

        if b′ === nothing #b′ = cache period
            b′ = b + cllog2(H)
        end

        bkt = undefs(Ti, 1 << b + 1) #bkt = scratch space for cumulative count
        cnt = undefs(Ti, 1 << b + 1, N >> b′ + 1, H) #cnt = cached cumulative counts

        byt = undefs(Ti, N)

        for h = H : -1 : 1 #h = level in tree
            for i′ = 1 : 1 << (h * b) : N + 1 #i′ = "quotient" of sorts
                for d = 1:(1 << b) + 1
                    bkt[d] = max(i′ + (d - 1) << ((h - 1) * b) - 1, 1)
                end
                for q = max(i′ - 1, 1) : min(i′ + (1 << (h * b)) - 1, N + 1) - 1
                    i = idx[q]
                    d = i >> ((h - 1) * b) & (1 << b - 1) + 1
                    q′ = bkt[d]
                    byt[q′] = (idx[q′] & ~(1 << ((h - 1) * b) - 1)) | (i & (1 << ((h - 1) * b) - 1))
                    bkt[d] = q′ + 1
                end
            end
            cnt[1, 1, h] = 0
            for d = 1 : 1 << b
                cnt[d + 1, 1, h] = 0
            end
            zero!(bkt)
            for q = 1 : N
                i = idx[q]
                d = i >> ((h - 1) * b) & ((1 << b) - 1) + 1 #d = a "remainder" of sorts
                bkt[d] += 1
                if q & (1 << b′ - 1) == 0 #Cache bkt at the end of each 1 << b′ block
                    Q = (q >> b′) + 1
                    cnt[1, Q, h] = 0
                    for d = 1 : 1 << b
                        cnt[d + 1, Q, h] = bkt[d] + cnt[d, Q, h]
                    end
                end
            end

            idx, byt = byt, idx
        end
        byt = idx
        return RookCount{Ti}(N, b, b′, H, byt, cnt)
    end
end

Base.getindex(arg::RookCount{Ti}, i::Integer, j::Integer) where {Ti} = arg(i, j)
function (arg::RookCount{Ti})(i::Integer, j::Integer) where {Ti}
    @inbounds begin
        N = arg.N
        b = arg.b
        b′ = arg.b′
        H = arg.H
        byt = arg.byt
        cnt = arg.cnt

        Δq = j - 1
        i = Ti(i - 1)::Ti
        s = Ti(0)
        for h = H : -1 : 2
            i′ = i & ~(1 << (h * b) - 1) + 1
            q₁ = max(i′ - 1, 1) - 1
            q₂ = q₁ + Δq
            d = i >> ((h - 1) * b) & ((1 << b) - 1) + 1
            Q₁ = (q₁ >> b′) + 1
            Q₂ = (q₂ >> b′) + 1
            s += cnt[d, Q₂, h] - cnt[d, Q₁, h]
            Δq = (cnt[d + 1, Q₂, h] - cnt[d, Q₂, h]) - (cnt[d + 1, Q₁, h] - cnt[d, Q₁, h])
            msk = ((1 << b) - 1) << ((h - 1) * b)
            cmp = (d - 1) << ((h - 1) * b)
            for q = (Q₁ - 1) << b′ + 1 : q₁
                d′ = byt[q] & msk
                s -= d′ < cmp
                Δq -= d′ == cmp
                #d′ = byt[q] >> ((h - 1) * b) & ((1 << b) - 1) + 1
                #s -= d′ < d
                #Δq -= d′ == d
            end
            for q = (Q₂ - 1) << b′ + 1 : q₂
                d′ = byt[q] & msk
                s += d′ < cmp
                Δq += d′ == cmp
                #d′ = byt[q] >> ((h - 1) * b) & ((1 << b) - 1) + 1
                #s += d′ < d
                #Δq += d′ == d
            end
        end

        i′ = i & ~(1 << b - 1) + 1
        q₁ = max(i′ - 1, 1) - 1
        q₂ = q₁ + Δq

        d = i & ((1 << b) - 1) + 1
        Q₁ = (q₁ >> b′) + 1
        Q₂ = (q₂ >> b′) + 1
        s += cnt[d + 1, Q₂, 1] - cnt[d + 1, Q₁, 1]
        msk = (1 << b) - 1
        cmp = d - 1
        for q = (Q₁ - 1) << b′ + 1 : q₁
            d′ = byt[q] & msk
            s -= d′ <= cmp
            #d′ = byt[q] & ((1 << b) - 1) + 1
            #s -= d′ <= d
        end
        for q = (Q₂ - 1) << b′ + 1 : q₂
            d′ = byt[q] & msk
            s += d′ <= cmp
            #d′ = byt[q] & ((1 << b) - 1) + 1
            #s += d′ <= d
        end
        return s
    end
end

BinaryRookCount(hint::AbstractHint, N, idx; kwargs...) =
    BinaryRookCount{UInt}(hint, N, idx; kwargs...)
BinaryRookCount{Tb}(hint::AbstractHint, N, idx::Vector{Ti}; kwargs...) where {Tb, Ti} =
    BinaryRookCount{Tb, Ti}(hint, N, idx; kwargs...)
function BinaryRookCount{Tb, Ti}(hint::AbstractHint, N, idx::Vector{Ti}; kwargs...) where {Tb, Ti}
    @inbounds begin
        #b = branching factor of tree = 1

        #H = height of tree
        H = cllog2(N + 1)

        cnt = zeros(Ti, 1 + cld(N, nbits(Tb)), H) #cnt = cached cumulative counts

        byt = zeros(Tb, 1 + cld(N, nbits(Tb)), H)
        idx′ = undefs(Ti, N)

        for h = H : -1 : 1 #h = level in tree
            _cnt = 0
            for i′ = 1 : 1 << h : N + 1 #i′ = "quotient" of sorts
                bkt_1 = max(i′ - 1, 1)
                bkt_2 = max(i′ + 1 << (h - 1) - 1, 1)
                for q = max(i′ - 1, 1) : min(i′ + (1 << h) - 1, N + 1) - 1
                    i = idx[q]
                    d = i >> (h - 1) & 1
                    q′ = ifelse(d == 0, bkt_1, bkt_2)
                    idx′[q′] = i
                    bkt_1 += 1 - d
                    bkt_2 += d

                    Q = ((q - 1) >> log2nbits(Tb)) + 1
                    byt[Q, h] |= d << ((q - 1) & (nbits(Tb) - 1))
                    cnt[Q + 1, h] = _cnt += d
                end
            end

            idx, idx′ = idx′, idx
        end
        return BinaryRookCount{Tb, Ti}(N, H, byt, cnt)
    end
end

Base.getindex(arg::BinaryRookCount{Tb, Ti}, i::Integer, j::Integer) where {Ti, Tb} = arg(i, j)
function (arg::BinaryRookCount{Tb, Ti})(i::Integer, j::Integer) where {Ti, Tb}
    @inbounds begin
        H = arg.H
        byt = arg.byt
        cnt = arg.cnt

        Δq = j - 1
        i = Ti(i - 1)::Ti
        s = Ti(0)
        for h = H : -1 : 1
            i′ = i & ~(1 << h - 1) + 1
            q₁ = max(i′ - 1, 1) - 1
            q₂ = q₁ + Δq
            d = i >> (h - 1) & Ti(1)

            Q₁ = (q₁ >> log2nbits(Tb)) + 1
            Q₂ = (q₂ >> log2nbits(Tb)) + 1
            bkt_2 = cnt[Q₂, h] - cnt[Q₁, h]
            bkt_2 += count_ones(byt[Q₂, h] & ((one(Tb) << (q₂ & (nbits(Tb) - 1))) - 1))
            bkt_2 -= count_ones(byt[Q₁, h] & ((one(Tb) << (q₁ & (nbits(Tb) - 1))) - 1))
            bkt_1 = Δq - bkt_2
            s += ifelse(d == 0, 0, bkt_1)
            Δq = ifelse(d == 0, bkt_1, bkt_2)
        end

        return s + Δq
    end
end


text:
@testset "Lattices" begin
  include("Lattices/Lattices.jl")
  include("Lattices/Absolute.jl")
end