text:
"""
    abstract type CellField <: GridapType end

This is a plain array of fields plus some metadata.
It is interpreted as a `Field` object on each cell of a computational
mesh.
"""
abstract type CellField <: GridapType end

"""
    get_array(cf::CellField)
"""
function get_array(cf::CellField)
  @abstractmethod
end

"""
    get_cell_map(cf::CellField)
"""
function get_cell_map(cf::CellField)
  @abstractmethod
end

"""
    get_cell_axes(cf::CellField)
"""
function get_cell_axes(cf::CellField)
  @abstractmethod
end

function get_memo(cf::CellField)
  @abstractmethod
end

"""
This trait returns `Val{true}()` when the `CellField` is defined in a
reference finite element space, and `Val{false}()` when it is defined in the
physical space
"""
RefStyle(::Type{<:CellField}) = @notimplemented

# We use duck typing here for all types marked with the RefStyle
RefStyle(::Type) = @notimplemented
RefStyle(a) = RefStyle(typeof(a))
is_in_ref_space(::Type{T}) where T = get_val_parameter(RefStyle(T))
is_in_ref_space(::T) where T = is_in_ref_space(T)
is_in_physical_space(::Type{T}) where T = !is_in_ref_space(T)
is_in_physical_space(a::T) where T = !is_in_ref_space(T)

"""
This trait provides information about the field size at a single
evaluation point.

For physical fields `MetaSizeStyle(T) == Val( () )`
For test bases `MetaSizeStyle(T) == Val( (:,) )`
For trial bases `MetaSizeStyle(T) = Val( (1,:) )`
For fields representing elemental matrices `MetaSizeStyle(T) = Val( (:,:) )`
"""
MetaSizeStyle(::Type{<:CellField}) = @notimplemented

# We use duck typing here for all types marked with the MetaSizeStyle trait
MetaSizeStyle(::Type) = @notimplemented
MetaSizeStyle(::T) where T = MetaSizeStyle(T)
get_metasize(::Type{T}) where T = get_val_parameter(MetaSizeStyle(T))
get_metasize(::T) where T = get_val_parameter(MetaSizeStyle(T))
is_test(::Type{T}) where T = get_metasize(T) == (:,)
is_test(::T) where T = is_test(T)
is_trial(::Type{T}) where T = get_metasize(T) == (1,:)
is_trial(::T) where T = is_trial(T)
is_basis(::Type{T}) where T =  is_test(T) || is_trial(T)
is_basis(::T) where T = is_basis(T)

# Tester

"""
    test_cell_field(
      cf::CellField,
      x::AbstractArray,
      b::AbstractArray,
      pred=(==);
      grad=nothing)
"""
function test_cell_field(cf::CellField,x::AbstractArray,b::AbstractArray,pred=(==);grad=nothing)
  cell_map = get_cell_map(cf)
  @test isa(cell_map,AbstractArray)
  a = evaluate(cf,x)
  test_array(a,b,pred)
  if grad != nothing
    g = evaluate(gradient(cf),x)
    test_array(g,grad,pred)
  end
  rs = RefStyle(cf)
  @test isa(get_val_parameter(rs),Bool)
  _cf = change_ref_style(cf)
  @test get_array(_cf) === get_array(cf)
  @test is_in_ref_space(cf) == !is_in_ref_space(_cf)
  @test is_in_physical_space(cf) == !is_in_physical_space(_cf)
  @test get_metasize(_cf) == get_metasize(cf)
  @test get_cell_axes(_cf) == get_cell_axes(cf)
  @test isa(get_cell_axes(cf),AbstractArray{<:Tuple})
  @test isa(get_metasize(cf),Tuple)
  @test isa(MetaSizeStyle(cf),Val)
  @test isa(get_memo(cf),Dict)
end

# Concrete implementation

"""
struct GenericCellField{R,S} <: CellField
  # Private fields
end

"""
struct GenericCellField{R,S} <: CellField
  array::AbstractArray
  cell_map::AbstractArray
  ref_trait::Val{R}
  cell_axes::AbstractArray
  metasize::Val{S}
  memo::Dict
  function GenericCellField(
    array::AbstractArray,
    cell_map::AbstractArray,
    ref_trait::Val{R},
    cell_axes::AbstractArray,
    metasize::Val{S}) where {R,S}
    memo = Dict()
    new{R,S}(array,cell_map,ref_trait,cell_axes,metasize,memo)
  end
end

function GenericCellField(array::AbstractArray,cell_map::AbstractArray)
  GenericCellField(array,cell_map,Val{true}(),Fill((),length(cell_map)))
end

function GenericCellField(array::AbstractArray,cell_map::AbstractArray,ref_trait::Val)
  GenericCellField(array,cell_map,ref_trait,Fill((),length(cell_map)))
end

function GenericCellField(
  array::AbstractArray,cell_map::AbstractArray,ref_trait::Val,cell_axes::AbstractArray)
  GenericCellField(array,cell_map,ref_trait,cell_axes,Val(()))
end

function get_array(cf::GenericCellField)
  cf.array
end

function get_cell_map(cf::GenericCellField)
  cf.cell_map
end

function get_memo(cf::GenericCellField)
  cf.memo
end

function get_cell_axes(cf::GenericCellField)
  cf.cell_axes
end

function RefStyle(::Type{<:GenericCellField{R}}) where {R}
  Val{R}()
end

function MetaSizeStyle(::Type{<:GenericCellField{R,S}}) where {R,S}
  Val{S}()
end

# The rest of the file is some default API using the CellField interface

# Preserving and manipulating the metadata

"""
    similar_object(cf::CellField,array::AbstractArray,cell_axes::AbstractArray,msize_style::Val)
"""
function similar_object(cf::CellField,array::AbstractArray,cell_axes::AbstractArray,msize_style::Val)
  cm = get_cell_map(cf)
  similar_object(cf,array,cm,cell_axes,msize_style)
end

function similar_object(
  cf::CellField,array::AbstractArray,cell_map::AbstractArray,cell_axes::AbstractArray,msize_style::Val)
  GenericCellField(array,cell_map,RefStyle(cf),cell_axes,msize_style)
end

"""
   change_ref_style(cf::CellField)

Return an object with the same array and metadata as in `cf`, except for `RefStyle` which is changed.
"""
function change_ref_style(cf::CellField)
  ref_sty = RefStyle(cf)
  bool = !get_val_parameter(ref_sty)
  new_sty = Val{bool}()
  ar = get_array(cf)
  cm = get_cell_map(cf)
  GenericCellField(ar,cm,new_sty,get_cell_axes(cf),MetaSizeStyle(cf))
end

# Move between ref and phys spaces by using the underling cell_map

to_ref_space(a::CellField) = _to_ref_space(a,RefStyle(a))
_to_ref_space(a,::Val{true}) = a
function _to_ref_space(a,::Val{false})
  cell_map = get_cell_map(a)
  array = compose(  get_array(a), cell_map  )
  no = similar_object(a,array,get_cell_axes(a),MetaSizeStyle(a))
  change_ref_style(no)
end

to_physical_space(a::CellField) = _to_physical_space(a,RefStyle(a))
_to_physical_space(a,::Val{true}) = @notimplemented # and probably not doable in some cases
_to_physical_space(a,::Val{false}) = a

# Assumption : x ALWAIS defined in the reference space
# In the future we can also add the RefStyle to x by defining CellPoints
"""
    evaluate(cf::CellField,x)
"""
function evaluate(cf::CellField,x::AbstractArray)
  key = (:evaluate,objectid(x))
  memo = get_memo(cf)
  if !haskey(memo,key)
    memo[key] = _evaluate(cf,x,RefStyle(cf))
  end
  memo[key]
end

function _evaluate(cf::CellField,x::AbstractArray,::Val{true})
  evaluate_field_array(get_array(cf),x)
end

function _evaluate(cf::CellField,x::AbstractArray,::Val{false})
  cm = get_cell_map(cf)
  _x = evaluate(cm,x)
  evaluate_field_array(get_array(cf),_x)
end

"""
    length(cf::CellField)
"""
function Base.length(cf::CellField)
  a = get_array(cf)
  length(a)
end

function Arrays.reindex(cf::CellField,a::AbstractVector)
  array = reindex(get_array(cf),a)
  cell_axes = reindex(get_cell_axes(cf),a)
  cell_map = reindex(get_cell_map(cf),a)
  similar_object(cf,array,cell_map,cell_axes,MetaSizeStyle(cf))
end

# Bases-related
function trialize_cell_basis(test::CellField)
  _trialize_cell_basis(test,MetaSizeStyle(test))
end

function _trialize_cell_basis(test,metasize)
  @unreachable
end

function _trialize_cell_basis(trial,metasize::Val{(1,:)})
  trial
end

function _trialize_cell_basis(test::CellField,metasize::Val{(:,)})
  array = trialize_array_of_bases(get_array(test))
  axs = apply(Fields._add_singleton_block,get_cell_axes(test))
  similar_object(test,array,axs,Val((1,:)))
end

function Fields.lincomb(a::CellField,b::AbstractArray)
  @notimplementedif !is_test(a)
  array = lincomb(get_array(a),b)
  axs = Fill((),length(b))
  similar_object(a,array,axs,Val(()))
end

# Diff ops

"""
    gradient(cf::CellField)
"""
function gradient(cf::CellField)
  key = :gradient
  memo = get_memo(cf)
  if !haskey(memo,key)
    a = get_array(cf)
    ag = field_array_gradient(a)
    memo[key] = similar_object(cf,ag,get_cell_axes(cf),MetaSizeStyle(cf))
  end
  memo[key]
end

# Operations

function operate(op,cf::CellField)
  key = objectid(op)
  memo = get_memo(cf)
  if !haskey(memo,key)
    a = get_array(cf)
    b = operate_arrays_of_fields(Fields._UnimplementedField,op,a)
    memo[key] = similar_object(cf,b,get_cell_axes(cf),MetaSizeStyle(cf))
  end
  memo[key]
end

function operate(op,cf1::CellField,cf2::CellField)
  @assert length(cf1) == length(cf2)
  @assert RefStyle(cf1) == RefStyle(cf2)
  a1 = get_array(cf1)
  a2 = get_array(cf2)
  b = operate_arrays_of_fields(Fields._UnimplementedField,op,a1,a2)
  axs = apply(field_operation_axes,get_cell_axes(cf1),get_cell_axes(cf2))
  metasize = field_operation_metasize(get_metasize(cf1),get_metasize(cf2))
  similar_object(cf1,b,axs,Val(metasize))
end

function operate(op,cf1::CellField,object)
  cm = get_cell_map(cf1)
  cf2 = convert_to_cell_field(object,cm,RefStyle(cf1))
  operate(op,cf1,cf2)
end

function operate(op,object,cf2::CellField)
  cm = get_cell_map(cf2)
  cf1 = convert_to_cell_field(object,cm,RefStyle(cf2))
  operate(op,cf1,cf2)
end

function operate(op,args::CellField...)
  a1 = first(args)
  @assert all( map( a->length(a) == length(a1), args) )
  @assert all( map( a->RefStyle(a)==RefStyle(a1), args) )
  arrs = map(get_array,args)
  m = operate_arrays_of_fields(Fields._UnimplementedField,op,arrs...)
  axs = apply(field_operation_axes,map(get_cell_axes,args)...)
  metasize = field_operation_metasize(map(get_metasize,args)...)
  similar_object(a1,m,axs,Val(metasize))
end

function operate(op,a1::CellField,args...)
  cm = get_cell_map(a1)
  rs = RefStyle(cm)
  cfs = map(obj->convert_to_cell_field(obj,cm,rs),args)
  operate(op,a1,cfs...)
end

# Conversions

"""
    convert_to_cell_field(object,cell_map,ref_style)
"""
function convert_to_cell_field(object,cell_map,ref_style::Val)
  @abstractmethod
end

function convert_to_cell_field(object,cell_map)
  ref_style = Val{true}()
  convert_to_cell_field(object,cell_map,ref_style)
end

function convert_to_cell_field(object::CellField,cell_map)
  @assert length(object) == length(cell_map)
  object
end

# pre-defined conversions

function convert_to_cell_field(object::CellField,cell_map,ref_style::Val)
  @assert RefStyle(object) == ref_style
  object
end

function convert_to_cell_field(object::AbstractArray,cell_map,ref_style::Val)
  @assert length(object) == length(cell_map)
  GenericCellField(object,cell_map,ref_style)
end

function convert_to_cell_field(object::Function,cell_map,ref_style::Val{true})
  b = compose(object,cell_map)
  GenericCellField(b,cell_map,Val{true}())
end

function convert_to_cell_field(fun::Function,cell_map,ref_style::Val{false})
  field = function_field(fun)
  cell_field = Fill(field,length(cell_map))
  GenericCellField(cell_field,cell_map,Val{false}())
end

function convert_to_cell_field(object::Number,cell_map,ref_style::Val)
  a = Fill(object,length(cell_map))
  GenericCellField(a,cell_map,ref_style)
end

# Skeleton related

"""
    struct SkeletonCellField <: GridapType
      left::CellField
      right::CellField
    end

Supports the same differential and algebraic operations than [`CellField`](@ref)
"""
struct SkeletonCellField <: GridapType
  left::CellField
  right::CellField
end

get_inward(a::SkeletonCellField) = a.left

get_outward(a::SkeletonCellField) = a.right

function Base.getproperty(x::SkeletonCellField, sym::Symbol)
  if sym in (:inward,:⁺)
    get_inward(x)
  elseif sym in (:outward,:⁻)
    get_outward(x)
  else
    getfield(x, sym)
  end
end

"""
    get_cell_map(a::SkeletonCellField)
"""
function get_cell_map(a::SkeletonCellField)
  get_cell_map(a.left)
end

"""
    jump(sf::SkeletonCellField)
"""
function jump(sf::SkeletonCellField)
  sf.⁺ - sf.⁻
end

"""
    mean(sf::SkeletonCellField)
"""
function mean(sf::SkeletonCellField)
  operate(_mean,sf.⁺,sf.⁻)
end

_mean(x,y) = 0.5*x + 0.5*y

function gradient(cf::SkeletonCellField)
  left = gradient(cf.left)
  right = gradient(cf.right)
  SkeletonCellField(left,right)
end

function operate(op,cf::SkeletonCellField)
  left = operate(op,cf.left)
  right = operate(op,cf.right)
  SkeletonCellField(left,right)
end

function operate(op,cf1::SkeletonCellField,cf2::SkeletonCellField)
  left = operate(op,cf1.left,cf2.left)
  right = operate(op,cf1.right,cf2.right)
  SkeletonCellField(left,right)
end

function operate(op,cf1::SkeletonCellField,cf2::CellField)
  left = operate(op,cf1.left,cf2)
  right = operate(op,cf1.right,cf2)
  SkeletonCellField(left,right)
end

function operate(op,cf1::CellField,cf2::SkeletonCellField)
  left = operate(op,cf1,cf2.left)
  right = operate(op,cf1,cf2.right)
  SkeletonCellField(left,right)
end

function operate(op,cf1::SkeletonCellField,object)
  cm = get_cell_map(cf1)
  cf2 = convert_to_cell_field(object,cm,RefStyle(cf1.left))
  operate(op,cf1,cf2)
end

function operate(op,object,cf2::SkeletonCellField)
  cm = get_cell_map(cf2)
  cf1 = convert_to_cell_field(object,cm,RefStyle(cf2.left))
  operate(op,cf1,cf2)
end

function merge_cell_dofs_at_skeleton(idsL,idsR,axesL,axesR)
  blocks = (idsL,idsR)
  blockids = [(1,),(2,)]
  axs = create_array_of_blocked_axes(axesL,axesR)
  VectorOfBlockArrayCoo(blocks,blockids,axs)
end

function merge_cell_fields_at_skeleton(cfL,cfR)
  @assert is_basis(cfL) == is_basis(cfR)
  if !is_basis(cfL)
    SkeletonCellField(cfL,cfR)
  else
    ax1 = get_cell_axes(cfL)
    ax2 = get_cell_axes(cfR)
    arrL = insert_array_of_bases_in_block(1,get_array(cfL),ax1,ax2)
    cfSL = similar_object(cfL,arrL,arrL.axes,MetaSizeStyle(cfL))
    arrR = insert_array_of_bases_in_block(2,get_array(cfR),ax1,ax2)
    cfSR = similar_object(cfR,arrR,arrR.axes,MetaSizeStyle(cfR))
    SkeletonCellField(cfSL,cfSR)
  end
end


text:
using TupleVectors
using AbstractMCMC

using AbstractDifferentiation

struct IsNested end

struct IsMCMC end

struct HasDeriv end
struct NoDeriv end

export sample

function AbstractMCMC.sample(post::Posterior, sampler::S, args...; init_params=nothing, kwargs...) where {S}
    θ0 = init_params
    if isnothing(init_params)
        θ0 = rand(post.prior)
    end
    return _sample(samplertype(S), post, sampler, args...; init_params=θ0, kwargs...)
end

function _sample(::IsMCMC, post, sampler, args...; init_params, kwargs...)
    tpost = asflat(post)
    θ0 = inverse(tpost.transform, init_params)
    return sample(tpost, sampler, args...; init_params=θ0, kwargs...)
end

function _sample(::IsNested, post, sampler, args...; init_params, kwargs...)
    tpost = ascube(post)
    θ0 = inverse(tpost.transform, init_params)
    return sample(tpost, sampler, args...; init_params=θ0, kwargs...)
end

ext: jl


text:
using ImmutableArrays
function scalematrix{T}(scale::Vector3{T})
    result      = eye(T, 4, 4)
    result[1,1] = scale[1]
    result[2,2] = scale[2]
    result[3,3] = scale[3]

    return Matrix4x4(result)
end
immutable AABB{T}
  min::Vector3{T}
  max::Vector3{T}
end
function Base.max{T, NDIM}(x::Array{Vector3{T},NDIM})
    reduce(x) do v0, v1
        Vector3(v0[1] > v1[1] ? v0[1] : v1[1],
            v0[2] > v1[2] ? v0[2] : v1[2],
            v0[3] > v1[3] ? v0[3] : v1[3])
    end
end

function maxper(v0::Vector3, v1::Vector3)
	return Vector3(max(v0[1], v1[1]),
            max(v0[2], v1[2]),
            max(v0[3], v1[3]))
end
function minper(v0::Vector3, v1::Vector3)
	return Vector3(min(v0[1], v1[1]),
            min(v0[2], v1[2]),
            min(v0[3], v1[3]))
end

function Base.min{T, NDIM}(x::Array{Vector3{T},NDIM}) = reduce(minper, x)
function Base.max{T, NDIM}(x::Array{Vector3{T},NDIM}) = reduce(maxper, x)

function Base.min{T, NDIM}(x::Array{Vector4{T},NDIM})
    reduce(x) do v0, v1
        Vector4(v0[1] < v1[1] ? v0[1] : v1[1],
            v0[2] < v1[2] ? v0[2] : v1[2],
            v0[3] < v1[3] ? v0[3] : v1[3],
            v0[4] < v1[4] ? v0[4] : v1[4])
    end
end
const Vec3 = Vector3{Float32}
const Vec2 = Vector2{Float32}
const Vec4 = Vector4{Float32}
const Mat4 = Matrix4x4{Float32}
function genquad{T}(downleft::Vector3{T}, width::Vector3{T}, height::Vector3{T})
    v = Vector3{T}[
        downleft,
        downleft + height,
        downleft + width + height,
        downleft + width
    ]
    uv = Vector2{T}[
        Vector2{T}(0, 1),
        Vector2{T}(0, 0),
        Vector2{T}(1, 0),
        Vector2{T}(1, 1)
    ]
    indexes = Uint32[0,1,2,2,3,0]

    normal = unit(cross(width, height))
    (v, uv, Vector3{T}[normal for i=1:4], indexes)
end

COLOR_QUAD = genquad(Vec3(0, 0, 0), Vec3(1, 0, 0), Vec3(0, 1, 0))
color = Dict(
  :middle                   => Vec2(0.5, 0.5),

  :swatchsize               => 0.1f0,
  :border_size              => 0.02f0,

  :antialiasing_value       => 0.01f0,
  :model                    => scalematrix(Vec3(200,200,1))
  )
function color_chooser_boundingbox(obj)
  middle      = obj[:middle]
  swatchsize  = obj[:swatchsize]
  border_size = obj[:border_size]
  model       = obj[:model]
  verts       = COLOR_QUAD[1]
  verts = map(verts) do x
    Vec3(model*Vec4(x...,0f0))
  end
  AABB(min(verts), max(verts))
end

abstract AbstractFixedVector{T,C}
immutable GLGlyph{T} <: AbstractFixedVector{T, 4}
  glyph::T
  line::T
  row::T
  style_group::T
end

function GLGlyph(glyph::Integer, line::Integer, row::Integer, style_group::Integer)
  if !isascii(char(glyph))
    glyph = char('1')
  end
  GLGlyph{Uint16}(uint16(glyph), uint16(line), uint16(row), uint16(style_group))
end
function GLGlyph(glyph::Char, line::Integer, row::Integer, style_group::Integer)
  if !isascii(glyph)
    glyph = char('1')
  end
  GLGlyph{Uint16}(uint16(glyph), uint16(line), uint16(row), uint16(style_group))
end

GLGlyph() = GLGlyph(' ', typemax(Uint16), typemax(Uint16), 0)

Base.length{T}(::GLGlyph{T})                   = 4
Base.length{T}(::Type{GLGlyph{T}})             = 4
Base.eltype{T}(::GLGlyph{T})                   = T
Base.eltype{T}(::Type{GLGlyph{T}})             = T
Base.size{T}(::GLGlyph{T})                     = (4,)

Base.start{T}(::GLGlyph{T})                    = 1
Base.next{T}(x::GLGlyph{T}, state::Integer)    = (getfield(x, state), state+1)
Base.done{T}(x::GLGlyph{T}, state::Integer)    = state > 4

function setindex1D!{T <: AbstractFixedVector, ElType}(a::Union(Matrix{T}, Vector{T}), x::ElType, i::Integer, accessor::Integer)
    if length(a) < i
        error("Out of Bounds. 1D index: ", i, " Matrix: ", typeof(a), " length: ", length(a), " size: ", size(a))
    end
    cardinality = length(T)
    if length(accessor) > cardinality
        error("Out of Bounds. 1D index: ", i, " Matrix: ", typeof(a), " length: ", length(a), " size: ", size(a))
    end

  ptr = convert(Ptr{eltype(T)}, pointer(a))
  unsafe_store!(ptr, convert(eltype(T), x), ((i-1)*cardinality)+accessor)
end
function setindex1D!{T <: AbstractFixedVector, ElType}(a::Union(Matrix{T}, Vector{T}), x::Vector{ElType}, i::Integer, accessor::UnitRange)
   if length(a) < i
     error("Out of Bounds. 1D index: ", i, " Matrix: ", typeof(a), " length: ", length(a), " size: ", size(a))
   end
   cardinality = length(T)
   if length(accessor) > cardinality
     error("Out of Bounds. 1D index: ", i, " Matrix: ", typeof(a), " length: ", length(a), " size: ", size(a))
   end
   eltp  = eltype(T)
   x     = convert(Vector{eltp}, x)
   ptr   = convert(Ptr{eltp}, pointer(a))
   unsafe_copy!(ptr + (sizeof(eltp)*((i-1)*(cardinality)+first(accessor-1))), pointer(x), length(accessor))
end

import Base: (+)

function (+){T}(a::Array{GLGlyph{T}, 1}, b::GLGlyph{T})
  for i=1:length(a)
    a[i] = a[i] + b
  end
end
function (+){T}(a::GLGlyph{T}, b::GLGlyph{T})
  GLGlyph{T}(a.glyph + b.glyph, a.line + b.line, a.row + b.row, a.style_group + b.style_group)
end
Base.utf16(glypharray::Array{GLGlyph{Uint16}}) = utf16(Uint16[c.glyph for c in glypharray])
Base.utf8(glypharray::Array{GLGlyph{Uint16}})  = utf8(Uint8[uint8(c.glyph) for c in glypharray])

function escape_regex(x::String)
  result = ""
  for elem in x
      if elem in regex_literals
          result *= string('\\')
      end
      result *= string(elem)
  end
  result
end
regreduce(arr, prefix="(", suffix=")") = Regex(reduce((v0, x) -> v0*"|"*prefix*escape_regex(x)*suffix, prefix*escape_regex(arr[1])*suffix, arr[2:end]))

function update_glyphpositions!{T}(text_array::AbstractArray{GLGlyph{T}}, start=1, stop=length(text_array))
  line = text_array[start].line
  row  = text_array[start].row
  for i=1:stop
    glyph = text_array[i].glyph
    setindex1D!(text_array, T[line, row], i, 2:3)
    if glyph == '\n'
      row = zero(T)
      line += one(T)
    else
      row += one(T)
    end
  end
end

function makedisplayable(text::String)
  result = map(collect(text)) do x
    str = string(x)
    if !is_valid_utf8(str)
      return utf8([one(Uint8)]) # replace with something that yields a missing symbol
    elseif str == "\r"
      return "\n"
    else
      return str == "\t" ? utf8(" "^tab) : utf8(str) # also replace tabs
    end
  end
  join(result)
end

function toglypharray(text::String, tab=3)
  #@assert is_valid_utf16(text) # future support for utf16
  text = makedisplayable(text)
  #Allocate some more memory, to reduce growing the texture residing on VRAM
  texturesize = div(length(text),     1024)+1 # a texture size of 1024 should be supported on every GPU
  text_array  = Array(GLGlyph{Uint16}, 1024, texturesize)
  setindex1D!(text_array, 1, 1, 2) # set first line
  setindex1D!(text_array, 0, 1, 3) # set first row
  #Set text
  for (i, elem) in enumerate(text)
    setindex1D!(text_array, uint16(char(elem)), i, 1)
    setindex1D!(text_array, 0, i, 4)
  end
  update_glyphpositions!(text_array) # calculate glyph positions
  text_array
end

text = Dict(
  :scalingfactor    => Vec2(1.0, 1.5), #Multiplicator for advance, newline
  :model            => eye(Mat4),
  :newline          => -Vec3(0, 24, 0),
  :advance          => Vec3(12, 0, 0),
  :text 			=> toglypharray("alkjsdölkjdsöfljk\n"^15),
  :textlength 		=> length("alkjsdölkjdsöfljk\n"^15)
)

function text_boundingbox(obj)
  glypharray  = obj[:text]
  advance  	  = obj[:advance]
  newline  	  = obj[:newline]

  maxv = Vector3(typemin(Float32))
  minv = Vector3(typemax(Float32))
  for elem in glypharray[1:obj[:textlength]]
  	println(elem.row)
  	currentpos = elem.row*advance + elem.line*newline
  	maxv = maxper(maxv, currentpos)
  	minv = minper(minv, currentpos)
  end
  AABB(minv, maxv)
end

println(text_boundingbox(text))


text:
# Compute densities of states

# IDOS (integrated density of states)
# N(ε) = sum_n f_n = sum_n f((εn-ε)/temperature)
# DOS (density of states)
# D(ε) = N'(ε)
# LDOS (local density of states)
# LD = sum_n f_n ψn

@doc raw"""
    compute_nos(ε, basis, eigenvalues; smearing=basis.model.smearing,
                temperature=basis.model.temperature)

The number of Kohn-Sham states in a temperature window of width `temperature` around the
energy `ε` contributing to the DOS at temperature `T`.

This quantity is not a physical quantity, but rather a dimensionless approximate measure
for how well properties near the Fermi surface are sampled with the passed `smearing`
and temperature `T`. It increases with both `T` and better sampling of the BZ with
``k``-points. A value ``\gg 1`` indicates a good sampling of properties near the
Fermi surface.
"""
function compute_nos(ε, basis, eigenvalues; smearing=basis.model.smearing,
                     temperature=basis.model.temperature)
    N = zeros(typeof(ε), basis.model.n_spin_components)
    if (temperature == 0) || smearing isa Smearing.None
        error("compute_nos only supports finite temperature")
    end

    # Note the differences to the DOS and LDOS functions: We are not counting states
    # per BZ volume (like in DOS), but absolute number of states. Therefore n_symeqk
    # is used instead of kweigths. We count the states inside a temperature window
    # `temperature` centred about εik. For this the states are weighted by the distribution
    # -f'((εik - ε)/temperature).
    #
    # To explicitly show the similarity with DOS and the temperature dependence we employ
    # -f'((εik - ε)/temperature) = temperature * ( d/dε f_τ(εik - ε') )|_{ε' = ε}
    for σ in 1:basis.model.n_spin_components, ik = krange_spin(basis, σ)
        n_symeqk = length(basis.ksymops[ik])  # Number of symmetry-equivalent k-points
        for (iband, εnk) in enumerate(eigenvalues[ik])
            enred = (εnk - ε) / temperature
            N[σ] -= n_symeqk * Smearing.occupation_derivative(smearing, enred)
        end
    end
    N = mpi_sum(N, basis.comm_kpts)
end

"""
Total density of states at energy ε
"""
function compute_dos(ε, basis, eigenvalues; smearing=basis.model.smearing,
                     temperature=basis.model.temperature)
    if (temperature == 0) || smearing isa Smearing.None
        error("compute_dos only supports finite temperature")
    end
    filled_occ = filled_occupation(basis.model)

    D = zeros(typeof(ε), basis.model.n_spin_components)
    for σ in 1:basis.model.n_spin_components, ik = krange_spin(basis, σ)
        for (iband, εnk) in enumerate(eigenvalues[ik])
            enred = (εnk - ε) / temperature
            D[σ] -= (filled_occ * basis.kweights[ik] / temperature
                     * Smearing.occupation_derivative(smearing, enred))
        end
    end
    D = mpi_sum(D, basis.comm_kpts)
end

"""
Local density of states, in real space
"""
function compute_ldos(ε, basis, eigenvalues, ψ; smearing=basis.model.smearing,
                      temperature=basis.model.temperature)
    if (temperature == 0) || smearing isa Smearing.None
        error("compute_ldos only supports finite temperature")
    end
    filled_occ = filled_occupation(basis.model)

    weights = deepcopy(eigenvalues)
    for (ik, εk) in enumerate(eigenvalues)
        for (iband, εnk) in enumerate(εk)
            enred = (εnk - ε) / temperature
            weights[ik][iband] = (-filled_occ / temperature
                                  * Smearing.occupation_derivative(smearing, enred))
        end
    end

    # Use compute_density routine to compute LDOS, using just the modified
    # weights (as "occupations") at each k-point. Note, that this automatically puts in the
    # required symmetrization with respect to kpoints and BZ symmetry
    compute_density(basis, ψ, weights)
end

"""
Plot the density of states over a reasonable range
"""
function plot_dos end


text:
module EquationsOfState

abstract type EquationOfStateOfSolidsParameters{T} end
const Parameters = EquationOfStateOfSolidsParameters

abstract type EquationOfState{T<:Parameters} end

end


text:
"A string whose data is stored in a Blob."
struct BlobString <: AbstractString
    data::Blob{UInt8}
    len::Int64 # in bytes
end

Base.pointer(blob::BlobString) = pointer(blob, 1)
function Base.pointer(blob::BlobString, i::Integer)
    # TODO(jamii) would prefer to boundscheck on load, but this will do for now
    getindex(blob.data + (blob.len - 1))
    pointer(blob.data + (i-1))
end

function Base.unsafe_copyto!(blob::BlobString, string::Union{BlobString, String})
    @assert blob.len >= sizeof(string)
    unsafe_copyto!(pointer(blob), pointer(string), sizeof(string))
end

function Base.String(blob::BlobString)
    unsafe_string(pointer(blob), blob.len)
end

# string interface - this is largely copied from Base

Base.sizeof(s::BlobString) = s.len

Base.ncodeunits(s::BlobString) = sizeof(s)
Base.codeunit(s::BlobString) = UInt8

@inline function Base.codeunit(s::BlobString, i::Integer)
    @boundscheck checkbounds(s, i)
    GC.@preserve s unsafe_load(pointer(s, i))
end

## comparison ##

function Base.cmp(a::BlobString, b::BlobString)
    al, bl = sizeof(a), sizeof(b)
    c = ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, UInt),
    a, b, min(al,bl))
    return c < 0 ? -1 : c > 0 ? +1 : cmp(al,bl)
end

function Base.:(==)(a::BlobString, b::BlobString)
    al = sizeof(a)
    al == sizeof(b) && 0 == ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, UInt), a, b, al)
end

## thisind, nextind ##

Base.thisind(s::BlobString, i::Int) = Base._thisind_str(s, i)

Base.nextind(s::BlobString, i::Int) = Base._nextind_str(s, i)

## checking UTF-8 & ACSII validity ##

byte_string_classify(s::BlobString) =
ccall(:u8_isvalid, Int32, (Ptr{UInt8}, Int), s, sizeof(s))
# 0: neither valid ASCII nor UTF-8
# 1: valid ASCII
# 2: valid UTF-8

## required core functionality ##

Base.@propagate_inbounds function Base.iterate(s::BlobString, i::Int=firstindex(s))
    i > ncodeunits(s) && return nothing
    b = codeunit(s, i)
    u = UInt32(b) << 24
    Base.between(b, 0x80, 0xf7) || return reinterpret(Char, u), i+1
    return next_continued(s, i, u)
end

function next_continued(s::BlobString, i::Int, u::UInt32)
    u < 0xc0000000 && (i += 1; @goto ret)
    n = ncodeunits(s)
    # first continuation byte
    (i += 1) > n && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 16
    # second continuation byte
    ((i += 1) > n) | (u < 0xe0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 8
    # third continuation byte
    ((i += 1) > n) | (u < 0xf0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b); i += 1
    @label ret
    return reinterpret(Char, u), i
end

Base.@propagate_inbounds function Base.getindex(s::BlobString, i::Int)
    b = codeunit(s, i)
    u = UInt32(b) << 24
    Base.between(b, 0x80, 0xf7) || return reinterpret(Char, u)
    return getindex_continued(s, i, u)
end

function getindex_continued(s::BlobString, i::Int, u::UInt32)
    if u < 0xc0000000
        # called from `getindex` which checks bounds
        @inbounds isvalid(s, i) && @goto ret
        Base.string_index_err(s, i)
    end
    n = ncodeunits(s)

    (i += 1) > n && @goto ret
    @inbounds b = codeunit(s, i) # cont byte 1
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 16

    ((i += 1) > n) | (u < 0xe0000000) && @goto ret
    @inbounds b = codeunit(s, i) # cont byte 2
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 8

    ((i += 1) > n) | (u < 0xf0000000) && @goto ret
    @inbounds b = codeunit(s, i) # cont byte 3
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b)
    @label ret
    return reinterpret(Char, u)
end

Base.getindex(s::BlobString, r::UnitRange{<:Integer}) = s[Int(first(r)):Int(last(r))]

function Base.getindex(s::BlobString, r::UnitRange{Int})
    isempty(r) && return ""
    i, j = first(r), last(r)
    @boundscheck begin
        checkbounds(s, r)
        @inbounds isvalid(s, i) || Base.string_index_err(s, i)
        @inbounds isvalid(s, j) || Base.string_index_err(s, j)
    end
    j = nextind(s, j) - 1
    # n = j - i + 1
    # ss = _string_n(n)
    # p = pointer(ss)
    # for k = 1:n
    #     unsafe_store!(p, codeunit(s, i + k - 1), k)
    # end
    # return ss
    BlobString(s.data + (i-1), j-i+1)
end

function Base.length(s::BlobString, i::Int, j::Int)
    @boundscheck begin
        0 < i ≤ ncodeunits(s)+1 || throw(BoundsError(s, i))
        0 ≤ j < ncodeunits(s)+1 || throw(BoundsError(s, j))
    end
    j < i && return 0
    @inbounds i, k = thisind(s, i), i
    c = j - i + (i == k)
    _length(s, i, j, c)
end

Base.length(s::BlobString) = _length(s, 1, ncodeunits(s), ncodeunits(s))

@inline function _length(s::BlobString, i::Int, n::Int, c::Int)
    i < n || return c
    @inbounds b = codeunit(s, i)
    @inbounds while true
        while true
            (i += 1) ≤ n || return c
            0xc0 ≤ b ≤ 0xf7 && break
            b = codeunit(s, i)
        end
        l = b
        b = codeunit(s, i) # cont byte 1
        c -= (x = b & 0xc0 == 0x80)
        x & (l ≥ 0xe0) || continue

        (i += 1) ≤ n || return c
        b = codeunit(s, i) # cont byte 2
        c -= (x = b & 0xc0 == 0x80)
        x & (l ≥ 0xf0) || continue

        (i += 1) ≤ n || return c
        b = codeunit(s, i) # cont byte 3
        c -= (b & 0xc0 == 0x80)
    end
end

## overload methods for efficiency ##

Base.isvalid(s::BlobString, i::Int) = checkbounds(Bool, s, i) && thisind(s, i) == i


text:
#include("should/not/appear1.jl")
#include("shouldnt/appear2.jl")

include("some/path/file1.jl");
include("another/path/file2.jl")


text:
module Lava

include("LavaCore.jl")
using .LavaCore
include("features.jl")
using .features

const lava = LavaCore
export lava
export features
export VkExt
export Utils
export Camera

end


text:
module AbstractYieldMacros

export @yield, @yieldfrom

using ContextualMacros

ContextualMacros.@contextualmacro @yield
ContextualMacros.@contextualmacro @yieldfrom

end


text:
function run(count::Int32)
    start = time()
    x = zero(Int32)
    while x < count
        x = ccall((:plusone, "./newplus/libnewplus.so"), Int32, (Cint,), x)
    end
    return round(Int32, (time() - start) * 1000)
end

function start()
    if isempty(ARGS)
        println("First arg (0 - 2000000000) is required.")
        return nothing
    end
    count = parse(Int32, ARGS[1])
    if count <= 0 || count > 2*10^9
        println("Must be a positive number not exceeding 2 billion.")
        return nothing
    end
    duration = run(count)
    println(duration)
end

start()

max_stars_repo_name: minond/ffi-overhead

max_issues_repo_name: minond/ffi-overhead

max_forks_repo_name: minond/ffi-overhead


text:
"""
    struct EffectiveChargeDensity{T, N, S, AT} <: AbstractArray{T, N}

Effective charge density needed to calculate the [`ElectricPotential`](@ref).
The effective charge density is the charge density (in C/m³) times the volume of the voxel of the respective
grid point (in m³). Thus, the unit of the effective charge density is Coulomb (C).

## Parametric types
* `T`: Element type of `data`.
* `N`: Dimension of the `grid` and `data` array.
* `S`: Coordinate system (`Cartesian` or `Cylindrical`).
* `AT`: Axes type.

## Fields
* `data::Array{T, N}`: Array containing the values of the effective charge density at the discrete points of the `grid`.
* `grid::Grid{T, N, S, AT}`: [`Grid`](@ref) defining the discrete points at which the electric potential is determined.
"""
struct EffectiveChargeDensity{T, N, S, AT} <: AbstractArray{T, N}
    data::Array{T, N}
    grid::Grid{T, N, S, AT}
end

@inline size(ρ::EffectiveChargeDensity{T, N, S}) where {T, N, S} = size(ρ.data)
@inline length(ρ::EffectiveChargeDensity{T, N, S}) where {T, N, S} = length(ρ.data)
@inline getindex(ρ::EffectiveChargeDensity{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = getindex(ρ.data, I...)
@inline getindex(ρ::EffectiveChargeDensity{T, N, S}, i::Int) where {T, N, S} = getindex(ρ.data, i)
@inline getindex(ρ::EffectiveChargeDensity{T, N, S}, s::Symbol) where {T, N, S} = getindex(ρ.grid, s)

function NamedTuple(ρ::EffectiveChargeDensity{T}) where {T <: SSDFloat}
    return (
        grid = NamedTuple(ρ.grid),
        values = ρ.data * internal_voltage_unit,
    )
end
Base.convert(T::Type{NamedTuple}, x::EffectiveChargeDensity) = T(x)

function EffectiveChargeDensity(nt::NamedTuple)
    grid = Grid(nt.grid)
    T = typeof(ustrip(nt.values[1]))
    S = get_coordinate_system(grid)
    N = get_number_of_dimensions(grid)
    EffectiveChargeDensity{T, N, S, typeof(grid.axes)}( ustrip.(uconvert.(internal_voltage_unit, nt.values)), grid)
end
Base.convert(T::Type{EffectiveChargeDensity}, x::NamedTuple) = T(x)


text:
struct TimeVaryingHMM{F,T} <: AbstractHMM{F}
    a::Vector{T}
    A::AbstractArray{T,3}
    B::Matrix{Distribution{F}}
    TimeVaryingHMM{F,T}(a, A, B) where {F,T} = assert_timevaryinghmm(a, A, B) && new(a, A, B)
end

TimeVaryingHMM(
    a::AbstractVector{T},
    A::AbstractArray{T,3},
    B::AbstractMatrix{<:Distribution{F}},
) where {F,T} = TimeVaryingHMM{F,T}(a, A, B)

TimeVaryingHMM(A::AbstractArray{T,3}, B::AbstractMatrix{<:Distribution{F}}) where {F,T} =
    TimeVaryingHMM{F,T}(ones(size(A, 1)) ./ size(A, 1), A, B)

function assert_timevaryinghmm(a::AbstractVector, A::AbstractArray{T,3} where T, B::AbstractMatrix{<:Distribution})
    @argcheck isprobvec(a)
    @argcheck istransmats(A)
    @argcheck all(length.(B) .== length(B[1])) ArgumentError("All distributions must have the same dimensions")
    @argcheck size(A,3) + 1 == size(B,2) ArgumentError("Number of transition rates must match length of chain")
    @argcheck length(a) == size(A,1) == size(B,1)
    return true
end

istransmats(A::AbstractArray{T,3}) where {T} = all(i->istransmat(@view A[:,:,i]), 1:size(A,3))

==(h1::TimeVaryingHMM, h2::TimeVaryingHMM) = (h1.a == h2.a) && (h1.A == h2.A) && (h1.B == h2.B)

function rand(
    rng::AbstractRNG,
    hmm::TimeVaryingHMM;
    init = rand(rng, Categorical(hmm.a)),
    seq = false,
)
    T = size(hmm.B, 2)
    z = Vector{Int}(undef, T)
    (T >= 1) && (z[1] = init)
    for t = 2:T
        z[t] = rand(rng, Categorical(hmm.A[z[t-1],:,t-1]))
    end
    y = rand(rng, hmm, z)
    seq ? (z, y) : y
end

function rand(rng::AbstractRNG, hmm::TimeVaryingHMM{Univariate}, z::AbstractVector{<:Integer})
    y = Vector{Float64}(undef, length(z))
    for t in eachindex(z)
        y[t] = rand(rng, hmm.B[z[t],t])
    end
    y
end

function rand(
    rng::AbstractRNG,
    hmm::TimeVaryingHMM{Multivariate},
    z::AbstractVector{<:Integer},
)
    y = Matrix{Float64}(undef, length(z), size(hmm,2))
    for t in eachindex(z)
        y[t, :] = rand(rng, hmm.B[z[t],t])
    end
    y
end

"""
    size(hmm, [dim]) -> Int | Tuple

Return the number of states in `hmm`, the dimension of the observations and the length of the chain.
"""
size(hmm::TimeVaryingHMM, dim = :) = (size(hmm.B, 1), length(hmm.B[1]), size(hmm.B, 2))[dim]

copy(hmm::TimeVaryingHMM) = TimeVaryingHMM(copy(hmm.a), copy(hmm.A), copy(hmm.B))

function nparams(hmm::TimeVaryingHMM)
    (length(hmm.a) - 1) + (size(hmm.A,1) * size(hmm.A,2) - size(hmm.A,1)) * size(hmm.A,3)
end

###

function loglikelihoods!(LL::AbstractMatrix, hmm::TimeVaryingHMM{Univariate}, observations)
    T, K = size(observations, 1), size(hmm, 1)
    @argcheck size(LL) == (T, K)
    @argcheck T == size(hmm, 3)
    @inbounds for i in OneTo(K), t in OneTo(T)
        LL[t,i] = logpdf(hmm.B[i,t], observations[t])
    end
end

function loglikelihoods!(LL::AbstractMatrix, hmm::TimeVaryingHMM{Multivariate}, observations)
    T, K = size(observations, 1), size(hmm, 1)
    @argcheck size(LL) == (T, K)
    @argcheck T == size(hmm, 3)
    @inbounds for i in OneTo(K), t in OneTo(T)
        LL[t,i] = logpdf(hmm.B[i,t], view(observations, t, :))
    end
end

####

function update_A!(
    A::AbstractArray{T,3} where {T},
    ξ::AbstractArray,
    α::AbstractMatrix,
    β::AbstractMatrix,
    LL::AbstractMatrix,
)
    @argcheck size(α, 1) == size(β, 1) == size(LL, 1) == size(ξ, 1) == size(A,3) + 1
    @argcheck size(α, 2) ==
              size(β, 2) ==
              size(LL, 2) ==
              size(A, 1) ==
              size(A, 2) ==
              size(ξ, 2) ==
              size(ξ, 3)

    T, K = size(LL)

    @inbounds for t in OneTo(T - 1)
        m = vec_maximum(view(LL, t + 1, :))
        c = 0.0

        for i in OneTo(K), j in OneTo(K)
            ξ[t, i, j] = α[t, i] * A[i, j, t] * exp(LL[t+1, j] - m) * β[t+1, j]
            c += ξ[t, i, j]
        end

        for i in OneTo(K), j in OneTo(K)
            ξ[t, i, j] /= c
        end
    end

    fill!(A, 0.0)

    @inbounds for t in OneTo(T - 1)
        for i in OneTo(K)
            c = 0.0

            for j in OneTo(K)
                A[i, j, t] += ξ[t, i, j]
                c += A[i, j, t]
            end

            for j in OneTo(K)
                A[i, j, t] /= c
            end
        end
    end
end

# # In-place update of the observations distributions.
# function update_B!(B::AbstractMatrix, γ::AbstractMatrix, observations, estimator)
#     @argcheck size(γ, 1) == size(observations, 1)
#     @argcheck size(γ, 2) == size(B, 1)
#     K = length(B)
#     for i in OneTo(K)
#         if sum(γ[:, i]) > 0
#             B[i] = estimator(typeof(B[i]), permutedims(observations), γ[:, i])
#         end
#     end
# end

function fit_mle!(
    hmm::TimeVaryingHMM,
    observations;
    display = :none,
    maxiter = 100,
    tol = 1e-3,
    robust = false,
    estimator = fit_mle,
    fit_dists = false
)
    @argcheck display in [:none, :iter, :final]
    @argcheck maxiter >= 0
    @argcheck fit_dists === false "Not yet supported"

    T, K = size(observations, 1), size(hmm, 1)
    @argcheck T == size(hmm.B, 2)
    history = EMHistory(false, 0, [])

    # Allocate memory for in-place updates
    c = zeros(T)
    α = zeros(T, K)
    β = zeros(T, K)
    γ = zeros(T, K)
    ξ = zeros(T, K, K)
    LL = zeros(T, K)

    loglikelihoods!(LL, hmm, observations)
    robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

    forwardlog!(α, c, hmm.a, hmm.A, LL)
    backwardlog!(β, c, hmm.a, hmm.A, LL)
    posteriors!(γ, α, β)

    logtot = sum(c)
    (display == :iter) && println("Iteration 0: logtot = $logtot")

    for it = 1:maxiter
        update_a!(hmm.a, α, β)
        update_A!(hmm.A, ξ, α, β, LL)
        fit_dists && update_B!(hmm.B, γ, observations, estimator)

        # Ensure the "connected-ness" of the states,
        # this prevents case where there is no transitions
        # between two extremely likely observations.
        robust && (hmm.A .+= eps())

        @check isprobvec(hmm.a)
        @check istransmats(hmm.A)

        if fit_dists
            loglikelihoods!(LL, hmm, observations)
            robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
        end

        forwardlog!(α, c, hmm.a, hmm.A, LL)
        backwardlog!(β, c, hmm.a, hmm.A, LL)
        posteriors!(γ, α, β)

        logtotp = sum(c)
        (display == :iter) && println("Iteration $it: logtot = $logtotp")

        push!(history.logtots, logtotp)
        history.iterations += 1

        if abs(logtotp - logtot) < tol
            (display in [:iter, :final]) &&
                println("EM converged in $it iterations, logtot = $logtotp")
            history.converged = true
            break
        end

        logtot = logtotp
    end

    if !history.converged
        if display in [:iter, :final]
            println("EM has not converged after $(history.iterations) iterations, logtot = $logtot")
        end
    end

    history
end

# In-place forward pass, where α and c are allocated beforehand.
function forwardlog!(
    α::AbstractMatrix,
    c::AbstractVector,
    a::AbstractVector,
    A::AbstractArray{T,3} where {T},
    LL::AbstractMatrix,
)
    @argcheck size(α, 1) == size(LL, 1) == size(c, 1)
    @argcheck size(α, 2) == size(LL, 2) == size(a, 1) == size(A, 1) == size(A, 2)

    T, K = size(LL)
    @argcheck T == size(A,3) + 1

    fill!(α, 0.0)
    fill!(c, 0.0)

    m = vec_maximum(view(LL, 1, :))

    for j in OneTo(K)
        α[1, j] = a[j] * exp(LL[1, j] - m)
        c[1] += α[1, j]
    end

    for j in OneTo(K)
        α[1, j] /= c[1]
    end

    c[1] = log(c[1]) + m

    @inbounds for t = 2:T
        m = vec_maximum(view(LL, t, :))

        for j in OneTo(K)
            for i in OneTo(K)
                α[t,j] += α[t-1, i] * A[i,j,t-1]
            end
            α[t, j] *= exp(LL[t, j] - m)
            c[t] += α[t, j]
        end

        for j in OneTo(K)
            α[t, j] /= c[t]
        end

        c[t] = log(c[t]) + m
    end
end

# In-place backward pass, where β and c are allocated beforehand.
function backwardlog!(
    β::AbstractMatrix,
    c::AbstractVector,
    a::AbstractVector,
    A::AbstractArray{T,3} where {T},
    LL::AbstractMatrix,
)
    @argcheck size(β, 1) == size(LL, 1) == size(c, 1)
    @argcheck size(β, 2) == size(LL, 2) == size(a, 1) == size(A, 1) == size(A, 2)

    T, K = size(LL)
    @argcheck T == size(A,3) + 1
    L = zeros(K)
    (T == 0) && return

    fill!(β, 0.0)
    fill!(c, 0.0)

    for j in OneTo(K)
        β[end, j] = 1.0
    end

    @inbounds for t = T-1:-1:1
        m = vec_maximum(view(LL, t + 1, :))

        for i in OneTo(K)
            L[i] = exp(LL[t+1, i] - m)
        end

        for j in OneTo(K)
            for i in OneTo(K)
                β[t,j] += β[t+1, i] * A[j,i,t] * L[i]
            end
            c[t+1] += β[t, j]
        end

        for j in OneTo(K)
            β[t, j] /= c[t+1]
        end

        c[t+1] = log(c[t+1]) + m
    end

    m = vec_maximum(view(LL, 1, :))

    for j in OneTo(K)
        c[1] += a[j] * exp(LL[1, j] - m) * β[1, j]
    end

    c[1] = log(c[1]) + m
end

function forward(a::AbstractVector, A::AbstractArray{T,3} where {T}, LL::AbstractMatrix)
    m = Matrix{Float64}(undef, size(LL))
    c = Vector{Float64}(undef, size(LL, 1))
    forwardlog!(m, c, a, A, LL)
    m, sum(c)
end

function backward(a::AbstractVector, A::AbstractArray{T,3} where {T}, LL::AbstractMatrix)
    m = Matrix{Float64}(undef, size(LL))
    c = Vector{Float64}(undef, size(LL, 1))
    backwardlog!(m, c, a, A, LL)
    m, sum(c)
end

function posteriors(a::AbstractVector, A::AbstractArray{T,3} where {T}, LL::AbstractMatrix; kwargs...)
    α, _ = forward(a, A, LL; kwargs...)
    β, _ = backward(a, A, LL; kwargs...)
    posteriors(α, β)
end

function viterbilog!(
    T1::AbstractMatrix,
    T2::AbstractMatrix,
    z::AbstractVector,
    a::AbstractVector,
    A::AbstractArray{T,3} where {T},
    LL::AbstractMatrix,
)
    T, K = size(LL)
    @argcheck T == size(A,3) + 1
    (T == 0) && return

    fill!(T1, 0.0)
    fill!(T2, 0)

    al = log.(a)
    Al = log.(A)

    for i in OneTo(K)
        T1[1, i] = al[i] + LL[1, i]
    end

    @inbounds for t = 2:T
        for j in OneTo(K)
            amax = 0
            vmax = -Inf

            for i in OneTo(K)
                v = T1[t-1, i] + Al[i, j, t-1]
                if v > vmax
                    amax = i
                    vmax = v
                end
            end

            T1[t, j] = vmax + LL[t, j]
            T2[t, j] = amax
        end
    end

    z[T] = argmax(T1[T, :])
    for t = T-1:-1:1
        z[t] = T2[t+1, z[t+1]]
    end
end

function viterbi(a::AbstractVector, A::AbstractArray{T,3} where {T}, LL::AbstractMatrix)
    T1 = Matrix{Float64}(undef, size(LL))
    T2 = Matrix{Int}(undef, size(LL))
    z = Vector{Int}(undef, size(LL, 1))
    viterbilog!(T1, T2, z, a, A, LL)
    z
end


text:
# codon_set.jl
# ============
#
# A tiny immutable set type for DNA and RNA codons.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/NaturalSelection.jl/blob/master/LICENSE.md

# Type definition
# ---------------

"""
CodonSet is a primitive type that represents a set of all possible 64 codons.

It is useful to compactly represent such a set in a single 64 bit word, when
one is working with variations between biological sequences.
"""
primitive type CodonSet{T<:NucleicAcid} 64 end

# Conversion and construction
# ---------------------------

CodonSet{T}() where T <: NucleicAcid = reinterpret(CodonSet{T}, UInt64(0))

# Base operators
# --------------

@inline function Base.:|(x::CodonSet{T}, y::UInt64) where T<:NucleicAcid
    return reinterpret(CodonSet{T}, reinterpret(UInt64, x) | y)
end

@inline function Base.:|(x::CodonSet{T}, y::Kmer{T,3}) where T<:NucleicAcid
    return x | UInt64(1) << UInt64(y)
end

@inline function Base.:&(x::CodonSet{T}, y::CodonSet{T}) where T <: NucleicAcid
    return reinterpret(CodonSet{T}, reinterpret(UInt64, x) & reinterpret(UInt64, y))
end

# Iteration interface
# -------------------

@inline function getnext(x::CodonSet, start::Integer)
    return trailing_zeros(reinterpret(UInt64, x) >> start) + start
end

Base.iteratorsize(::Type{CodonSet}) = HasLength()
Base.iteratoreltype(::Type{CodonSet}) = HasEltype()
Base.eltype(::Type{CodonSet{T}}) where T<:NucleicAcid = Kmer{T,3}
Base.length(x::CodonSet) = count_ones(reinterpret(UInt64, x))
Base.start(x::CodonSet) = getnext(x, UInt64(0))

@inline function Base.next(x::CodonSet{T}, state::UInt64) where T<:NucleicAcid
    return Kmer{T,3}(state), getnext(x, state + UInt64(1))
end

@inline function Base.done(x::CodonSet, state::UInt64)
    return state > UInt64(63)
end


text:
using Test, OMEinsum
using OMEinsum: match_rule, Sum, Tr, DefaultRule, Diag, Duplicate,
                Permutedims, nopermute,
                Identity, SimpleBinaryRule

@testset "match rule" begin
    ixs = ((1,2), (2,3))
    iy = (1,3)
    @test match_rule(ixs, iy) == SimpleBinaryRule(ein"ij,jk->ik")
    @test match_rule(EinCode(ixs, iy)) == SimpleBinaryRule(ein"ij,jk->ik")

    ixs = ((1,2), (2,3), (3,4))
    iy = (1,4)
    @test match_rule(ixs, iy) == DefaultRule()

    ixs = ((1,1),)
    iy = ()
    @test match_rule(ixs, iy) == Tr()

    ixs = ((1,1,2),)
    iy = (2,)
    @test match_rule(ixs, iy) == DefaultRule()

    ixs = ((1,2),)
    iy = (1,)
    @test match_rule(ixs, iy) == Sum()

    ixs = ((1,2,3),)
    iy = (2,1,3)
    @test match_rule(ixs, iy) != Sum()

    ixs = ((1,2),(1,2),(1,2))
    iy = (1,2)
    @test match_rule(ixs, iy) == DefaultRule()

    ixs = ((1,2,3),)
    iy = (1,3,2)
    @test match_rule(ixs, iy) == Permutedims()

    ixs = ((1,2),)
    iy = (1,1)
    @test match_rule(ixs, iy) == DefaultRule()

    ix = (1,2,3)
    iy = (1,2)
    @test nopermute(ix,iy)

    ix = (1,2,3)
    iy = (2,1)
    @test !nopermute(ix,iy)

    # 3-argument match_rule
    @test match_rule(((1,1),), ()) == Tr()
    @test match_rule(((1,2),), ()) == Sum()
    @test match_rule(((1,1),), (1,)) == Diag()

    @test match_rule(((1,2),(2,3),(3,4),(4,5)), (1,5)) == DefaultRule()
    @test match_rule(((1,2),(2,3),(5,4),(3,4)), (1,5)) == DefaultRule()
    @test match_rule(((1,2),(2,3),(3,4),(3,5)), (1,5)) == DefaultRule()

    @test match_rule(((1,2,3,4),), (1,2)) == Sum()
    @test match_rule(((1,2,3,4),), (1,2)) == Sum()
    @test match_rule(((1,2,3,4),), (1,2,3,4)) == Identity()
    @test match_rule(((1,2,3,4),), (2,1)) == Sum()
    @test match_rule(((1,2,3,4),), (2,1,1)) == DefaultRule()

    @test match_rule(((1,2),), (2,1)) == Permutedims()
    @test match_rule(((1,2,3),), (2,1)) == Sum()
    @test match_rule(((3,3),), (1,3)) == DefaultRule()

    @test match_rule(((1,2),(1,2)), (1,2)) == DefaultRule()
    @test match_rule(((1,2),(2,1)), (1,2)) == DefaultRule()
    @test match_rule(((1,2),(1,2)), (1,2,1)) == DefaultRule()

    @test match_rule(((1,1,2),), (2,)) == DefaultRule()
    @test match_rule(((1,1,2,3),), (2,3)) == DefaultRule()
    @test match_rule(((1,1,2,3),), (3,2)) == DefaultRule()
    @test match_rule(((1,1,2),), (1,2,)) == Diag()

    @test match_rule(((1,2),(2,3)), (1,3)) == SimpleBinaryRule(ein"ij,jk->ik")
    @test match_rule(((1,2),(2,3)), (3,1)) == SimpleBinaryRule(ein"ij,jk->ki")
    @test match_rule(((2,1),(2,3)), (3,1)) == SimpleBinaryRule(ein"ji,jk->ki")
    @test match_rule(((2,1),(3,2)), (3,1)) == SimpleBinaryRule(ein"ji,kj->ki")
    @test match_rule(((1,2,3),(2,3,4)), (1,4)) == DefaultRule()

    @test match_rule(((1,2,3),), (1,2,3)) == Identity()
    @test match_rule(((1,2,3),), (1,3,2)) == Permutedims()

    @test match_rule(((1,2,1),), (1,2)) == Diag()
    @test match_rule(((1,2),), (1,2,1)) == Duplicate()
end

@testset "isbatchmul" begin
    @test match_rule(((1,2), (2,3)), (1,3)) == SimpleBinaryRule(ein"ij,jk->ik") # matmul
    @test match_rule(((1,2,3), (2,3)), (1,3)) == SimpleBinaryRule(ein"ijl,jl->il") # matvec-batched
    @test match_rule(((7,1,2,3), (2,4,3,7)), (1,4,3)) == DefaultRule() # matmat-batched
    @test match_rule(((3,), (3,)), (3,)) == SimpleBinaryRule(ein"l,l->l")  # pure batch
    @test match_rule(((3,1), (3,)), (3,1)) == DefaultRule()  # batched vector-vector
    @test match_rule(((2,2), (2,3)), (1,3)) == DefaultRule()
    @test match_rule(((2,1), (2,3)), (1,1)) == DefaultRule()
    @test match_rule(((1,2,3), (2,4,3)), (1,3)) == DefaultRule()
end

@testset "match_rule eye candies" begin
    @test match_rule(ein"ij,jk,kl->il") == DefaultRule()
end


text:
using ApproxFunOrthogonalPolynomials, ApproxFunBase, LinearAlgebra, Test
    import ApproxFunBase: testspace, recA, recB, recC
    import ApproxFunOrthogonalPolynomials: forwardrecurrence

@testset "Chebyshev" begin
    @testset "Forward recurrence" begin
        @test forwardrecurrence(Float64,Chebyshev(),0:9,1.0) == ones(10)
    end

    @testset "ChebyshevInterval" begin
        @test Fun(x->2,10)(.1) ≈ 2
        @test Fun(x->2)(.1) ≈ 2
        @test Fun(Chebyshev,Float64[]).([0.,1.]) ≈ [0.,0.]
        @test Fun(Chebyshev,[])(0.) ≈ 0.
        @test Fun(x->4,Chebyshev(),1).coefficients == [4.0]
        @test Fun(x->4).coefficients == [4.0]
        @test Fun(4).coefficients == [4.0]

        @test Fun(x->4).coefficients == [4.0]
        @test Fun(4).coefficients == [4.0]

        f = Fun(ChebyshevInterval(), [1])
        @test f(0.1) == 1

        ef = Fun(exp)
        @test ef(0.1) ≈ exp(0.1)

        for d in (ChebyshevInterval(),Interval(1.,2.),Segment(1.0+im,2.0+2im))
            testspace(Chebyshev(d))
        end

        ef = Fun(exp,ChebyshevInterval())

        @test ef == -(-ef)
        @test ef == (ef-1) + 1

        ef = Fun(exp)

        cf = Fun(cos)

        ecf = Fun(x->cos(x).*exp(x))
        eocf = Fun(x->cos(x)./exp(x))

        @test ef(.5) ≈ exp(.5)
        @test ecf(.123456) ≈ cos(.123456).*exp(.123456)
    end

    @testset "Algebra" begin
        ef = Fun(exp,ChebyshevInterval())

        @test ef == -(-ef)
        @test ef == (ef-1) + 1

        @test ef / 3 == (3 \ ef)

        cf = Fun(cos)

        ecf = Fun(x->cos(x).*exp(x))
        eocf = Fun(x->cos(x)./exp(x))

        @test ef(.5) ≈ exp(.5)
        @test ecf(.123456) ≈ cos(.123456).*exp(.123456)

        r=2 .* rand(100) .- 1

        @test maximum(abs,ef.(r)-exp.(r))<200eps()
        @test maximum(abs,ecf.(r).-cos.(r).*exp.(r))<200eps()

        @test (cf .* ef)(0.1) ≈ ecf(0.1) ≈ cos(0.1)*exp(0.1)
        @test domain(cf.*ef) ≈ domain(ecf)
        @test domain(cf.*ef) == domain(ecf)

        @test norm((ecf-cf.*ef).coefficients)<200eps()
        @test maximum(abs,(eocf-cf./ef).coefficients)<1000eps()
        @test norm(((ef/3).*(3/ef)-1).coefficients)<1000eps()
    end

    @testset "Diff and cumsum" begin
        ef = Fun(exp)
        cf = Fun(cos)
        @test norm((ef - ef').coefficients)<10E-11

        @test norm((ef - cumsum(ef)').coefficients) < 20eps()
        @test norm((cf - cumsum(cf)').coefficients) < 20eps()

        @test sum(ef)  ≈ 2.3504023872876028
        @test norm(ef)  ≈ 1.90443178083307
    end

    @testset "other domains" begin
        ef = Fun(exp,1..2)
        cf = Fun(cos,1..2)

        ecf = Fun(x->cos(x).*exp(x),1..2)
        eocf = Fun(x->cos(x)./exp(x),1..2)

        x=1.5
        @test ef(x) ≈ exp(x)

        r=rand(100) .+ 1
        @test maximum(abs,ef.(r)-exp.(r))<400eps()
        @test maximum(abs,ecf.(r).-cos.(r).*exp.(r))<100eps()
    end

    @testset "Other interval" begin
        ef = Fun(exp,1..2)
        cf = Fun(cos,1..2)

        ecf = Fun(x->cos(x).*exp(x),1..2)
        eocf = Fun(x->cos(x)./exp(x),1..2)

        r=rand(100) .+ 1
        x=1.5

        @test ef(x) ≈ exp(x)

        @test maximum(abs,ef.(r)-exp.(r))<400eps()
        @test maximum(abs,ecf.(r).-cos.(r).*exp.(r))<100eps()
        @test norm((ecf-cf.*ef).coefficients)<500eps()
        @test maximum(abs,(eocf-cf./ef).coefficients)<1000eps()
        @test norm(((ef/3).*(3/ef)-1).coefficients)<1000eps()

        ## Diff and cumsum
        @test norm((ef - ef').coefficients)<10E-11
        @test norm((ef - cumsum(ef)').coefficients) < 10eps()
        @test norm((cf - cumsum(cf)').coefficients) < 10eps()

        @test sum(ef) ≈ 4.670774270471604
        @test norm(ef) ≈ 4.858451087240335
    end

    @testset "Roots" begin
        f=Fun(x->sin(10(x-.1)))
        @test norm(f.(roots(f)))< 1000eps()

        @test_throws ArgumentError roots(Fun(zero))
        @test_throws ArgumentError roots(Fun(Chebyshev(),Float64[]))
    end

    @testset "Aliasing" begin
        f=Fun(x->cos(50acos(x)))
        @test norm(f.coefficients-Matrix(I,ncoefficients(f),ncoefficients(f))[:,51])<100eps()
    end

    @testset "Integer values" begin
        @test Fun(x->2,10)(.1) ≈ 2
        @test Fun(x->2)(.1) ≈ 2

        @test Fun(Chebyshev,Float64[]).([0.,1.]) ≈ [0.,0.]
        @test Fun(Chebyshev,[])(0.) ≈ 0.
    end

    @testset "large intervals #121" begin
        x = Fun(identity,0..10)
        f = sin(x^2)
        g = cos(x)
        @test exp(x)(0.1) ≈ exp(0.1)
        @test f(.1) ≈ sin(.1^2)

        x = Fun(identity,0..100)
        f = sin(x^2)
        @test ≈(f(.1),sin(.1^2);atol=1E-12)
    end

    @testset "Reverse" begin
        f=Fun(exp)
        @test ApproxFunBase.default_Fun(f, Chebyshev(Segment(1 , -1)), ncoefficients(f))(0.1) ≈ exp(0.1)
        @test Fun(f,Chebyshev(Segment(1,-1)))(0.1) ≈ f(0.1)
    end

    @testset "minimum/maximum" begin
        x=Fun()
        @test minimum(x) == -1
        @test maximum(x) == 1
    end

    @testset "Do not overresolve #7" begin
        @test ncoefficients(Fun(x->sin(400*pi*x),-1..1)) ≤ 1400
    end

    @testset "Bug from Trogdon" begin
        δ = .03 # should be less than 0.03
        @test 0.0 ∈ Domain(1-8*sqrt(δ)..1+8*sqrt(δ))
        @test 0.00001 ∈ Domain(1-8*sqrt(δ)..1+8*sqrt(δ))

        ϕfun = Fun(x -> 1/sqrt(2*pi*δ)*exp(-abs2.(x-1)/(2*δ)), 1-8sqrt(δ)..1+8sqrt(δ))
        ϕfun(0.00001) ≈ 1/sqrt(2*pi*δ)*exp(-abs2.(0.00001-1)/(2*δ))

        iϕfun = 1-cumsum(ϕfun)
        @test iϕfun(0.00001) ≈ 1
    end

    @testset "Large scaling" begin
        w = Fun(x -> 1e5/(x*x+1), 283.72074879785936 .. 335.0101119042838)
        @test w(leftendpoint(domain(w))) ≈ 1e5/(leftendpoint(domain(w))^2+1)
    end

    @testset "supremum norm" begin
        x = Fun()
        f = 1/(1 + 25*(x^2))
        @test norm(f, Inf) ≈ 1.0
    end

    @testset "Jacobi" begin
        S=Chebyshev()
        @test S.a==-0.5
        @test S.b==-0.5
    end
end


text:
struct TaskNotPresentError{T,U} <: Exception
    t :: T
    task :: U
end
function Base.showerror(io::IO, err::TaskNotPresentError)
    print(io, "could not find the task $(err.task) in the list $(err.t)")
end

"""
    AbstractConstrainedProduct{T, N, Q}

Supertype of [`ProductSplit`](@ref) and [`ProductSection`](@ref).
"""
abstract type AbstractConstrainedProduct{T, N, Q} end
Base.eltype(::AbstractConstrainedProduct{T}) where {T} = T

_niterators(::AbstractConstrainedProduct{<:Any, N}) where {N} = N

const IncreasingAbstractConstrainedProduct{T, N} =
    AbstractConstrainedProduct{T, N, <:NTuple{N, AbstractUnitRange}}

"""
    ProductSection{T, N, Q<:NTuple{N,AbstractRange}}

Iterator that loops over a specified section of the
outer product of ranges in. If the ranges are strictly increasing, the
iteration will be in reverse - lexicographic order.
Given `N` ranges, each element returned by the iterator will be
a tuple of length `N` with one element from each range.

See also: [`ProductSplit`](@ref)
"""
struct ProductSection{T, N, Q <: NTuple{N,AbstractRange}} <: AbstractConstrainedProduct{T, N, Q}
    iterators :: Q
    togglelevels :: NTuple{N, Int}
    firstind :: Int
    lastind :: Int

    function ProductSection(iterators::Tuple{Vararg{AbstractRange, N}}, togglelevels::NTuple{N, Int},
        firstind::Int, lastind::Int) where {N}

        # Ensure that all the iterators are strictly increasing
        all(x->step(x)>0, iterators) ||
        throw(ArgumentError("all the ranges need to be strictly increasing"))

        T = Tuple{map(eltype, iterators)...}

        new{T, N, typeof(iterators)}(iterators, togglelevels, firstind, lastind)
    end
end

function _cumprod(len::Tuple)
    (0, _cumprod(first(len), Base.tail(len))...)
end

_cumprod(::Integer,::Tuple{}) = ()
function _cumprod(n::Integer, tl::Tuple)
    (n, _cumprod(n*first(tl), Base.tail(tl))...)
end

function takedrop(ps::ProductSection)
    drop = ps.firstind - 1
    take = ps.lastind - ps.firstind + 1
    Iterators.take(Iterators.drop(Iterators.product(ps.iterators...), drop), take)
end

"""
    ProductSection(iterators::Tuple{Vararg{AbstractRange}}, inds::AbstractUnitRange)

Construct a `ProductSection` iterator that represents a 1D view of the outer product
of the ranges provided in `iterators`, with the range of indices in the view being
specified by `inds`.

# Examples
```jldoctest
julia> p = ParallelUtilities.ProductSection((1:3, 4:6), 5:8);

julia> collect(p)
4-element $(Vector{Tuple{Int, Int}}):
 (2, 5)
 (3, 5)
 (1, 6)
 (2, 6)

julia> collect(p) == collect(Iterators.product(1:3, 4:6))[5:8]
true
```
"""
function ProductSection(iterators::Tuple{Vararg{AbstractRange}}, inds::AbstractUnitRange)
    firstind, lastind = first(inds), last(inds)

    len = map(length, iterators)
    Nel = prod(len)
    1 <= firstind || throw(
        ArgumentError("the range of indices must start from a number ≥ 1"))
    lastind <= Nel || throw(
        ArgumentError("the maximum index must be less than or equal to the total number of elements = $Nel"))
    togglelevels = _cumprod(len)
    ProductSection(iterators, togglelevels, firstind, lastind)
end
ProductSection(::Tuple{}, ::AbstractUnitRange) = throw(ArgumentError("need at least one iterator"))

"""
    ProductSplit{T, N, Q<:NTuple{N,AbstractRange}}

Iterator that loops over a section of the outer product of ranges.
If the ranges are strictly increasing, the iteration is in reverse - lexicographic order.
Given `N` ranges, each element returned by the iterator will be
a tuple of length `N` with one element from each range.

See also: [`ProductSection`](@ref)
"""
struct ProductSplit{T, N, Q<:NTuple{N, AbstractRange}} <: AbstractConstrainedProduct{T, N, Q}
    ps :: ProductSection{T, N, Q}
    np :: Int
    p :: Int

    function ProductSplit(ps::ProductSection{T, N, Q}, np::Integer, p::Integer) where {T, N, Q}
        1 <= p <= np || throw(ArgumentError("processor rank out of range"))
        new{T, N, Q}(ps, np, p)
    end
end

function nelementsdroptake(len, np, p)
    d, r = divrem(len, np)
    drop = d*(p - 1) + min(r, p - 1)
    lastind = d*p + min(r, p)
    take = lastind - drop
    drop, take
end

"""
    ProductSplit(iterators::Tuple{Vararg{AbstractRange}}, np::Integer, p::Integer)

Construct a `ProductSplit` iterator that represents the outer product
of the iterators split over `np` workers, with this instance reprsenting
the values on the `p`-th worker.

!!! note
    `p` here refers to the rank of the worker, and is unrelated to the worker ID obtained by
    executing `myid()` on that worker.

# Examples
```jldoctest
julia> ParallelUtilities.ProductSplit((1:2, 4:5), 2, 1) |> collect
2-element $(Vector{Tuple{Int, Int}}):
 (1, 4)
 (2, 4)

julia> ParallelUtilities.ProductSplit((1:2, 4:5), 2, 2) |> collect
2-element $(Vector{Tuple{Int, Int}}):
 (1, 5)
 (2, 5)
```
"""
function ProductSplit(iterators::Tuple{Vararg{AbstractRange}}, np::Integer, p::Integer)
    # d, r = divrem(prod(length, iterators), np)
    # firstind = d*(p - 1) + min(r, p - 1) + 1
    # lastind = d*p + min(r, p)
    drop, take = nelementsdroptake(prod(length, iterators), np, p)
    firstind = drop + 1
    lastind = drop + take
    ProductSplit(ProductSection(iterators, firstind:lastind), np, p)
end
ProductSplit(::Tuple{}, ::Integer, ::Integer) = throw(ArgumentError("Need at least one iterator"))

takedrop(ps::ProductSplit) = takedrop(ProductSection(ps))

workerrank(ps::ProductSplit) = ps.p
Distributed.nworkers(ps::ProductSplit) = ps.np

ProductSection(ps::ProductSection) = ps
ProductSection(ps::ProductSplit) = ps.ps

getiterators(ps::AbstractConstrainedProduct) = ProductSection(ps).iterators
togglelevels(ps::AbstractConstrainedProduct) = ProductSection(ps).togglelevels

function Base.summary(io::IO, ps::AbstractConstrainedProduct)
    print(io, length(ps), "-element ", string(nameof(typeof(ps))))
end
function Base.show(io::IO, ps::AbstractConstrainedProduct)
    summary(io, ps)
    if !isempty(ps)
        print(io, " [", repr(first(ps)) * ", ... , " * repr(last(ps)), "]")
    end
end

Base.isempty(ps::AbstractConstrainedProduct) = (firstindexglobal(ps) > lastindexglobal(ps))

function Base.first(ps::AbstractConstrainedProduct)
    isempty(ps) && throw(ArgumentError("collection must be non - empty"))
    _first(getiterators(ps), childindex(ps, firstindexglobal(ps))...)
end

function _first(t::Tuple, ind::Integer, rest::Integer...)
    (1 <= ind <= length(first(t))) || throw(BoundsError(first(t), ind))
    (first(t)[ind], _first(Base.tail(t), rest...)...)
end
_first(::Tuple{}) = ()

function Base.last(ps::AbstractConstrainedProduct)
    isempty(ps) && throw(ArgumentError("collection must be non - empty"))
    _last(getiterators(ps), childindex(ps, lastindexglobal(ps))...)
end

function _last(t::Tuple, ind::Integer, rest::Integer...)
    (1 <= ind <= length(first(t))) || throw(BoundsError(first(t), ind))
    (first(t)[ind], _last(Base.tail(t), rest...)...)
end
_last(::Tuple{}) = ()

Base.length(ps::AbstractConstrainedProduct) = lastindex(ps)

Base.firstindex(ps::AbstractConstrainedProduct) = 1
Base.lastindex(ps::AbstractConstrainedProduct) = lastindexglobal(ps) - firstindexglobal(ps) + 1

firstindexglobal(ps::AbstractConstrainedProduct) = ProductSection(ps).firstind
lastindexglobal(ps::AbstractConstrainedProduct) = ProductSection(ps).lastind

# SplittablesBase interface
function SplittablesBase.halve(ps::AbstractConstrainedProduct)
    iter = getiterators(ps)
    firstind = firstindexglobal(ps)
    lastind = lastindexglobal(ps)
    nleft = length(ps) ÷ 2
    firstindleft = firstind
    lastindleft = firstind + nleft - 1
    firstindright = lastindleft + 1
    lastindright = lastind
    tl = togglelevels(ps)
    ProductSection(iter, tl, firstindleft, lastindleft),
    ProductSection(iter, tl, firstindright, lastindright)
end

"""
    childindex(ps::AbstractConstrainedProduct, ind)

Return a tuple containing the indices of the individual `AbstractRange`s
corresponding to the element that is present at index `ind` in the
outer product of the ranges.

!!! note
    The index `ind` corresponds to the outer product of the ranges, and not to `ps`.

# Examples
```jldoctest
julia> iters = (1:5, 2:4, 1:3);

julia> ps = ParallelUtilities.ProductSplit(iters, 7, 1);

julia> ind = 6;

julia> cinds = ParallelUtilities.childindex(ps, ind)
(1, 2, 1)

julia> v = collect(Iterators.product(iters...));

julia> getindex.(iters, cinds) == v[ind]
true
```

See also: [`childindexshifted`](@ref)
"""
function childindex(ps::AbstractConstrainedProduct, ind)
    tl = reverse(Base.tail(togglelevels(ps)))
    reverse(childindex(tl, ind))
end

function childindex(tl::Tuple, ind)
    t = first(tl)
    k = div(ind - 1, t)
    (k + 1, childindex(Base.tail(tl), ind - k*t)...)
end

# First iterator gets the final remainder
childindex(::Tuple{}, ind) = (ind,)

"""
    childindexshifted(ps::AbstractConstrainedProduct, ind)

Return a tuple containing the indices in the individual iterators
given an index of `ps`.

If the iterators `(r1, r2, ...)` are used to generate
`ps`, then return `(i1, i2, ...)` such that `ps[ind] == (r1[i1], r2[i2], ...)`.

# Examples
```jldoctest
julia> iters = (1:5, 2:4, 1:3);

julia> ps = ParallelUtilities.ProductSplit(iters, 7, 3);

julia> psind = 4;

julia> cinds = ParallelUtilities.childindexshifted(ps, psind)
(3, 1, 2)

julia> getindex.(iters, cinds) == ps[psind]
true
```

See also: [`childindex`](@ref)
"""
function childindexshifted(ps::AbstractConstrainedProduct, ind)
    childindex(ps, (ind - 1) + firstindexglobal(ps))
end

function Base.getindex(ps::AbstractConstrainedProduct, ind)
    1 <= ind <= length(ps) || throw(BoundsError(ps, ind))
    _getindex(ps, childindexshifted(ps, ind)...)
end
# This needs to be a separate function to deal with the case of a single child iterator, in which case
# it's not clear if the single index is for the ProductSplit or the child iterator

# This method asserts that the number of indices is correct
function _getindex(ps::AbstractConstrainedProduct{<:Any, N}, inds::Vararg{Integer, N}) where {N}
    _getindex(getiterators(ps), inds...)
end

function _getindex(t::Tuple, ind::Integer, rest::Integer...)
    (1 <= ind <= length(first(t))) || throw(BoundsError(first(t), ind))
    (first(t)[ind], _getindex(Base.tail(t), rest...)...)
end
_getindex(::Tuple{}, ::Integer...) = ()

function Base.iterate(ps::AbstractConstrainedProduct, state...)
    iterate(takedrop(ps), state...)
end

function _firstlastalongdim(ps::AbstractConstrainedProduct, dims,
    firstindchild::Tuple = childindex(ps, firstindexglobal(ps)),
    lastindchild::Tuple = childindex(ps, lastindexglobal(ps)))

    iter = getiterators(ps)[dims]

    fic = firstindchild[dims]
    lic = lastindchild[dims]

    first_iter = iter[fic]
    last_iter = iter[lic]

    (first_iter, last_iter)
end

function _checkrollover(ps::AbstractConstrainedProduct, dims,
    firstindchild::Tuple = childindex(ps, firstindexglobal(ps)),
    lastindchild::Tuple = childindex(ps, lastindexglobal(ps)))

    _checkrollover(getiterators(ps), dims, firstindchild, lastindchild)
end

function _checkrollover(t::Tuple, dims, firstindchild::Tuple, lastindchild::Tuple)
    if dims > 0
        return _checkrollover(Base.tail(t), dims - 1, Base.tail(firstindchild), Base.tail(lastindchild))
    end

    !_checknorollover(reverse(t), reverse(firstindchild), reverse(lastindchild))
end

function _checknorollover(t, firstindchild, lastindchild)
    iter = first(t)
    first_iter = iter[first(firstindchild)]
    last_iter = iter[first(lastindchild)]

    (last_iter == first_iter) &
        _checknorollover(Base.tail(t), Base.tail(firstindchild), Base.tail(lastindchild))
end
_checknorollover(::Tuple{}, ::Tuple{}, ::Tuple{}) = true

function _nrollovers(ps::AbstractConstrainedProduct, dims::Integer)
    dims == _niterators(ps) && return 0
    nelements(ps; dims = dims + 1) - 1
end

"""
    nelements(ps::AbstractConstrainedProduct{T, N, <:NTuple{N,AbstractUnitRange}}; dims::Integer) where {T,N}

Compute the number of unique values in the section of the `dims`-th range contained in `ps`.

The function is defined currently only for iterator products of `AbstractUnitRange`s.

# Examples
```jldoctest
julia> ps = ParallelUtilities.ProductSplit((1:5, 2:4, 1:3), 7, 3);

julia> collect(ps)
7-element $(Vector{Tuple{Int, Int, Int}}):
 (5, 4, 1)
 (1, 2, 2)
 (2, 2, 2)
 (3, 2, 2)
 (4, 2, 2)
 (5, 2, 2)
 (1, 3, 2)

julia> ParallelUtilities.nelements(ps, dims = 1)
5

julia> ParallelUtilities.nelements(ps, dims = 2)
3

julia> ParallelUtilities.nelements(ps, dims = 3)
2
```
"""
function nelements(ps::IncreasingAbstractConstrainedProduct; dims::Integer)
    1 <= dims <= _niterators(ps) || throw(ArgumentError("1 ⩽ dims ⩽ N=$(_niterators(ps)) not satisfied for dims=$dims"))

    iter = getiterators(ps)[dims]

    if _nrollovers(ps, dims) == 0
        st = first(ps)[dims]
        en = last(ps)[dims]
        stind = findfirst(isequal(st), iter)
        enind = findfirst(isequal(en), iter)
        nel = length(stind:enind)
    elseif _nrollovers(ps, dims) > 1
        nel = length(iter)
    else
        st = first(ps)[dims]
        en = last(ps)[dims]
        stind = findfirst(isequal(st), iter)
        enind = findfirst(isequal(en), iter)
        if stind > enind
            # some elements are missed out
            nel = length(stind:length(iter)) + length(1:enind)
        else
            nel = length(iter)
        end
    end
    return nel
end

"""
    maximumelement(ps::AbstractConstrainedProduct; dims::Integer)

Compute the maximum value of the section of the range number `dims` contained in `ps`.

# Examples
```jldoctest
julia> ps = ParallelUtilities.ProductSplit((1:2, 4:5), 2, 1);

julia> collect(ps)
2-element $(Vector{Tuple{Int, Int}}):
 (1, 4)
 (2, 4)

julia> ParallelUtilities.maximumelement(ps, dims = 1)
2

julia> ParallelUtilities.maximumelement(ps, dims = 2)
4
```
"""
function maximumelement(ps::IncreasingAbstractConstrainedProduct; dims::Integer)
    isempty(ps) && throw(ArgumentError("collection must be non - empty"))

    firstindchild = childindex(ps, firstindexglobal(ps))
    lastindchild = childindex(ps, lastindexglobal(ps))

    _, last_iter = _firstlastalongdim(ps, dims, firstindchild, lastindchild)

    v = last_iter

    # The last index will not roll over so this can be handled easily
    if dims == _niterators(ps)
        return v
    end

    if _checkrollover(ps, dims, firstindchild, lastindchild)
        iter = getiterators(ps)[dims]
        v = maximum(iter)
    end

    return v
end

function maximumelement(ps::IncreasingAbstractConstrainedProduct{<:Any, 1})
    isempty(ps) && throw(ArgumentError("range must be non - empty"))
    lastindchild = childindex(ps, lastindexglobal(ps))
    lic_dim = lastindchild[1]
    iter = getiterators(ps)[1]
    iter[lic_dim]
end

"""
    minimumelement(ps::AbstractConstrainedProduct; dims::Integer)

Compute the minimum value of the section of the range number `dims` contained in `ps`.

# Examples
```jldoctest
julia> ps = ParallelUtilities.ProductSplit((1:2, 4:5), 2, 1);

julia> collect(ps)
2-element $(Vector{Tuple{Int, Int}}):
 (1, 4)
 (2, 4)

julia> ParallelUtilities.minimumelement(ps, dims = 1)
1

julia> ParallelUtilities.minimumelement(ps, dims = 2)
4
```
"""
function minimumelement(ps::IncreasingAbstractConstrainedProduct; dims::Integer)
    isempty(ps) && throw(ArgumentError("collection must be non - empty"))

    firstindchild = childindex(ps, firstindexglobal(ps))
    lastindchild = childindex(ps, lastindexglobal(ps))

    first_iter, last_iter = _firstlastalongdim(ps, dims, firstindchild, lastindchild)

    v = first_iter

    # The last index will not roll over so this can be handled easily
    if dims == _niterators(ps)
        return v
    end

    if _checkrollover(ps, dims, firstindchild, lastindchild)
        iter = getiterators(ps)[dims]
        v = minimum(iter)
    end

    return v
end

function minimumelement(ps::IncreasingAbstractConstrainedProduct{<:Any, 1})
    isempty(ps) && throw(ArgumentError("range must be non - empty"))
    firstindchild = childindex(ps, firstindexglobal(ps))
    fic_dim = firstindchild[1]
    iter = getiterators(ps)[1]
    iter[fic_dim]
end

"""
    extremaelement(ps::AbstractConstrainedProduct; dims::Integer)

Compute the `extrema` of the section of the range number `dims` contained in `ps`.

# Examples
```jldoctest
julia> ps = ParallelUtilities.ProductSplit((1:2, 4:5), 2, 1);

julia> collect(ps)
2-element $(Vector{Tuple{Int, Int}}):
 (1, 4)
 (2, 4)

julia> ParallelUtilities.extremaelement(ps, dims = 1)
(1, 2)

julia> ParallelUtilities.extremaelement(ps, dims = 2)
(4, 4)
```
"""
function extremaelement(ps::IncreasingAbstractConstrainedProduct; dims::Integer)
    isempty(ps) && throw(ArgumentError("collection must be non - empty"))

    firstindchild = childindex(ps, firstindexglobal(ps))
    lastindchild = childindex(ps, lastindexglobal(ps))

    first_iter, last_iter = _firstlastalongdim(ps, dims, firstindchild, lastindchild)

    v = (first_iter, last_iter)
    # The last index will not roll over so this can be handled easily
    if dims == _niterators(ps)
        return v
    end

    if _checkrollover(ps, dims, firstindchild, lastindchild)
        iter = getiterators(ps)[dims]
        v = extrema(iter)
    end

    return v
end

function extremaelement(ps::IncreasingAbstractConstrainedProduct{<:Any, 1})
    isempty(ps) && throw(ArgumentError("collection must be non - empty"))
    firstindchild = childindex(ps, firstindexglobal(ps))
    lastindchild = childindex(ps, lastindexglobal(ps))
    fic_dim = firstindchild[1]
    lic_dim = lastindchild[1]
    iter = getiterators(ps)[1]

    (iter[fic_dim], iter[lic_dim])
end

for (f, g) in [(:maximumelement, :maximum), (:minimumelement, :minimum), (:extremaelement, :extrema)]
    @eval $f(ps::AbstractConstrainedProduct{<:Any, 1}) = $g(first, takedrop(ps))
    @eval $f(ps::AbstractConstrainedProduct; dims::Integer) = $g(x -> x[dims], takedrop(ps))
end

"""
    extremadims(ps::AbstractConstrainedProduct)

Compute the extrema of the sections of all the ranges contained in `ps`.
Functionally this is equivalent to

```julia
map(i -> extrema(ps, dims = i), 1:_niterators(ps))
```

but it is implemented more efficiently.

Returns a `Tuple` containing the `(min, max)` pairs along each
dimension, such that the `i`-th index of the result contains the `extrema` along the section of the `i`-th range
contained locally.

# Examples
```jldoctest
julia> ps = ParallelUtilities.ProductSplit((1:2, 4:5), 2, 1);

julia> collect(ps)
2-element $(Vector{Tuple{Int, Int}}):
 (1, 4)
 (2, 4)

julia> ParallelUtilities.extremadims(ps)
((1, 2), (4, 4))
```
"""
function extremadims(ps::AbstractConstrainedProduct)
    _extremadims(ps, 1, getiterators(ps))
end

function _extremadims(ps::AbstractConstrainedProduct, dims::Integer, iterators::Tuple)
    (extremaelement(ps; dims = dims), _extremadims(ps, dims + 1, Base.tail(iterators))...)
end
_extremadims(::AbstractConstrainedProduct, ::Integer, ::Tuple{}) = ()

"""
    extrema_commonlastdim(ps::AbstractConstrainedProduct{T, N, <:NTuple{N,AbstractUnitRange}}) where {T,N}

Return the reverse - lexicographic extrema of values taken from
ranges contained in `ps`, where the pairs of ranges are constructed
by concatenating the ranges along each dimension with the last one.

For two ranges this simply returns `([first(ps)], [last(ps)])`.

# Examples
```jldoctest
julia> ps = ParallelUtilities.ProductSplit((1:3, 4:7, 2:7), 10, 2);

julia> collect(ps)
8-element $(Vector{Tuple{Int, Int, Int}}):
 (3, 6, 2)
 (1, 7, 2)
 (2, 7, 2)
 (3, 7, 2)
 (1, 4, 3)
 (2, 4, 3)
 (3, 4, 3)
 (1, 5, 3)

julia> ParallelUtilities.extrema_commonlastdim(ps)
$((Tuple{Int,Int}[(1, 2), (6, 2)], Tuple{Int,Int}[(3, 3), (5, 3)]))
```
"""
function extrema_commonlastdim(ps::IncreasingAbstractConstrainedProduct)
    isempty(ps) && return nothing

    m = extremadims(ps)
    lastvar_min, lastvar_max = last(m)

    val_first = first(ps)
    val_last = last(ps)
    min_vals = collect(Base.front(val_first))
    max_vals = collect(Base.front(val_last))

    for val in ps
        val_rev = reverse(val)
        lastvar = first(val_rev)
        (lastvar_min < lastvar < lastvar_max) && continue

        for (ind, vi) in enumerate(Base.tail(val_rev))
            if lastvar == lastvar_min
                min_vals[_niterators(ps) - ind] = min(min_vals[_niterators(ps) - ind], vi)
            end
            if lastvar == lastvar_max
                max_vals[_niterators(ps) - ind] = max(max_vals[_niterators(ps) - ind], vi)
            end
        end
    end

    [(m, lastvar_min) for m in min_vals], [(m, lastvar_max) for m in max_vals]
end

_infullrange(val::T, ps::AbstractConstrainedProduct{T}) where {T} = _infullrange(val, getiterators(ps))

function _infullrange(val, t::Tuple)
    first(val) in first(t) && _infullrange(Base.tail(val), Base.tail(t))
end
_infullrange(::Tuple{}, ::Tuple{}) = true

"""
    indexinproduct(iterators::NTuple{N, AbstractRange}, val::NTuple{N, Any}) where {N}

Return the index of `val` in the outer product of `iterators`.
Return nothing if `val` is not present.

# Examples
```jldoctest
julia> iterators = (1:4, 1:3, 3:5);

julia> val = (2, 2, 4);

julia> ind = ParallelUtilities.indexinproduct(iterators, val)
18

julia> collect(Iterators.product(iterators...))[ind] == val
true
```
"""
function indexinproduct(iterators::NTuple{N, AbstractRange}, val::Tuple{Vararg{Any, N}}) where {N}
    all(map(in, val, iterators)) || return nothing

    ax = map(x -> 1:length(x), iterators)
    individual_inds = map((it, val) -> findfirst(isequal(val), it), iterators, val)

    LinearIndices(ax)[individual_inds...]
end

indexinproduct(::Tuple{}, ::Tuple{}) = throw(ArgumentError("need at least one iterator"))

function Base.in(val::T, ps::AbstractConstrainedProduct{T}) where {T}
    _infullrange(val, ps) || return false

    ind = indexinproduct(getiterators(ps), val)
    firstindexglobal(ps) <= ind <= lastindexglobal(ps)
end

function Base.in(val::T, ps::IncreasingAbstractConstrainedProduct{T}) where {T}
    _infullrange(val, ps) || return false
    ReverseLexicographicTuple(first(ps)) <= ReverseLexicographicTuple(val) <= ReverseLexicographicTuple(last(ps))
end

# This struct is just a wrapper to flip the tuples before comparing
struct ReverseLexicographicTuple{T<:Tuple}
    t :: T
end

Base.isless(a::ReverseLexicographicTuple{T}, b::ReverseLexicographicTuple{T}) where {T} = reverse(a.t) < reverse(b.t)
Base.isequal(a::ReverseLexicographicTuple, b::ReverseLexicographicTuple) = a.t == b.t

"""
    whichproc(iterators::Tuple{Vararg{AbstractRange}}, val::Tuple, np::Integer)

Return the processor rank that will contain `val` if the outer
product of the ranges contained in `iterators` is split evenly
across `np` processors.

# Examples
```jldoctest
julia> iters = (1:4, 2:3);

julia> np = 2;

julia> ParallelUtilities.ProductSplit(iters, np, 2) |> collect
4-element $(Vector{Tuple{Int, Int}}):
 (1, 3)
 (2, 3)
 (3, 3)
 (4, 3)

julia> ParallelUtilities.whichproc(iters, (2, 3), np)
2
```
"""
function whichproc(iterators::Tuple{AbstractRange, Vararg{AbstractRange}}, val, np::Integer)
    _infullrange(val, iterators) || return nothing
    np >= 1 || throw(ArgumentError("np must be >= 1"))
    np  == 1 && return 1

    # We may carry out a binary search as the iterators are sorted
    left, right = 1, np

    val_t = ReverseLexicographicTuple(val)

    while left < right
        mid = div(left + right, 2)
        ps = ProductSplit(iterators, np, mid)

        # If np is greater than the number of ntasks then it's possible
        # that ps is empty. In this case the value must be somewhere in
        # the previous workers. Otherwise each worker has some tasks and
        # these are sorted, so carry out a binary search

        if isempty(ps) || val_t < ReverseLexicographicTuple(first(ps))
            right = mid - 1
        elseif val_t > ReverseLexicographicTuple(last(ps))
            left = mid + 1
        else
            return mid
        end
    end

    return left
end

whichproc(ps::ProductSplit, val) = whichproc(getiterators(ps), val, ps.np)

# This function tells us the range of processors that would be involved
# if we are to compute the tasks contained in the list ps on np_new processors.
# The total list of tasks is contained in iterators, and might differ from
# getiterators(ps) (eg if ps contains a subsection of the parameter set)
"""
    procrange_recast(iterators::Tuple{Vararg{AbstractRange}}, ps, np_new::Integer)

Return the range of processor ranks that would contain the values in `ps` if
the outer produce of the ranges in `iterators` is split across `np_new`
workers.

The values contained in `ps` should be a subsection of the outer product of
the ranges in `iterators`.

# Examples
```jldoctest
julia> iters = (1:10, 4:6, 1:4);

julia> ps = ParallelUtilities.ProductSplit(iters, 5, 2);

julia> ParallelUtilities.procrange_recast(iters, ps, 10)
3:4
```
"""
function procrange_recast(iterators::Tuple{AbstractRange, Vararg{AbstractRange}}, ps::AbstractConstrainedProduct, np_new::Integer)
    isempty(ps) && return nothing

    procid_start = whichproc(iterators, first(ps), np_new)
    if procid_start === nothing
        throw(TaskNotPresentError(iterators, first(ps)))
    end
    if length(ps) == 1
        procid_end = procid_start
    else
        procid_end = whichproc(iterators, last(ps), np_new)
        if procid_end === nothing
            throw(TaskNotPresentError(iterators, last(ps)))
        end
    end

    return procid_start:procid_end
end

"""
    procrange_recast(ps::AbstractConstrainedProduct, np_new::Integer)

Return the range of processor ranks that would contain the values in `ps` if the
iterators used to construct `ps` were split across `np_new` processes.

# Examples
```jldoctest
julia> iters = (1:10, 4:6, 1:4);

julia> ps = ParallelUtilities.ProductSplit(iters, 5, 2); # split across 5 processes initially

julia> ParallelUtilities.procrange_recast(ps, 10) # If `iters` were spread across 10 processes
3:4
```
"""
function procrange_recast(ps::AbstractConstrainedProduct, np_new::Integer)
    procrange_recast(getiterators(ps), ps, np_new)
end

"""
    localindex(ps::AbstractConstrainedProduct{T}, val::T) where {T}

Return the index of `val` in `ps`. Return `nothing` if the value
is not found.

# Examples
```jldoctest
julia> ps = ParallelUtilities.ProductSplit((1:3, 4:5:20), 3, 2);

julia> collect(ps)
4-element $(Vector{Tuple{Int, Int}}):
 (2, 9)
 (3, 9)
 (1, 14)
 (2, 14)

julia> ParallelUtilities.localindex(ps, (3, 9))
2
```
"""
function localindex(ps::AbstractConstrainedProduct{T}, val::T) where {T}
    (isempty(ps) || val ∉ ps) && return nothing

    indflat = indexinproduct(getiterators(ps), val)
    indflat - firstindexglobal(ps) + 1
end

"""
    whichproc_localindex(iterators::Tuple{Vararg{AbstractRange}}, val::Tuple, np::Integer)

Return `(rank, ind)`, where `rank` is the
rank of the worker that `val` will reside on if the outer product
of the ranges in `iterators` is spread over `np` workers, and `ind` is
the index of `val` in the local section on that worker.

# Examples
```jldoctest
julia> iters = (1:4, 2:8);

julia> np = 10;

julia> ParallelUtilities.whichproc_localindex(iters, (2, 4), np)
(4, 1)

julia> ParallelUtilities.ProductSplit(iters, np, 4) |> collect
3-element $(Vector{Tuple{Int, Int}}):
 (2, 4)
 (3, 4)
 (4, 4)
```
"""
function whichproc_localindex(iterators::Tuple{Vararg{AbstractRange}}, val::Tuple, np::Integer)
    procid = whichproc(iterators, val, np)
    procid === nothing && return nothing
    index = localindex(ProductSplit(iterators, np, procid), val)
    index === nothing && return nothing
    return procid, index
end

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

"""
    dropleading(ps::AbstractConstrainedProduct{T, N, NTuple{N,AbstractUnitRange}}) where {T,N}

Return a `ProductSection` leaving out the first iterator contained in `ps`.
The range of values of the remaining iterators in the
resulting `ProductSection` will be the same as in `ps`.

# Examples
```jldoctest
julia> ps = ParallelUtilities.ProductSplit((1:5, 2:4, 1:3), 7, 3);

julia> collect(ps)
7-element $(Vector{Tuple{Int, Int, Int}}):
 (5, 4, 1)
 (1, 2, 2)
 (2, 2, 2)
 (3, 2, 2)
 (4, 2, 2)
 (5, 2, 2)
 (1, 3, 2)

julia> ParallelUtilities.dropleading(ps) |> collect
3-element $(Vector{Tuple{Int, Int}}):
 (4, 1)
 (2, 2)
 (3, 2)
```
"""
function dropleading(ps::IncreasingAbstractConstrainedProduct)
    isempty(ps) && throw(ArgumentError("need at least one iterator"))
    iterators = Base.tail(getiterators(ps))
    first_element = Base.tail(first(ps))
    last_element = Base.tail(last(ps))
    firstind = indexinproduct(iterators, first_element)
    lastind = indexinproduct(iterators, last_element)
    ProductSection(iterators, firstind:lastind)
end


text:
function varnetinjectiterate!(netinjection::A, variable::JumpVariable, time_periods::Int64, devices::Array{T}) where {A <: JumpExpressionMatrix, T <: PowerSystems.Generator}

    for t in 1:time_periods, d in devices

        isassigned(netinjection,  d.bus.number,t) ? JuMP.add_to_expression!(netinjection[d.bus.number,t], variable[d.name, t]) : netinjection[d.bus.number,t] = variable[d.name, t];

    end

end

function varnetinjectiterate!(netinjection::A, variable::JumpVariable, time_periods::Int64, devices::Array{T}) where {A <: JumpExpressionMatrix, T <: PowerSystems.ElectricLoad}

    for t in 1:time_periods, d in devices

        isassigned(netinjection,  d.bus.number,t) ? JuMP.add_to_expression!(netinjection[d.bus.number,t], -1*variable[d.name, t]) : netinjection[d.bus.number,t] = -1*variable[d.name, t];

    end

end

function varnetinjectiterate!(netinjection::A, variable_in::JumpVariable, variable_out::JumpVariable, time_range::UnitRange{Int64}, devices::Array{T}) where {A <: JumpExpressionMatrix, T <: PowerSystems.Storage}

        for t in time_range, d in devices

            isassigned(netinjection,  d.bus.number,t) ? JuMP.add_to_expression!(netinjection[d.bus.number,t], variable_in[d.name,t] - variable_out[d.name,t]) : netinjection[d.bus.number,t] = variable_in[d.name,t] - variable_out[d.name,t];

        end

end


text:
using CUDA
include("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\kernelEvolutions.jl")
include("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\utils\\CUDAGpuUtils.jl")
using ..BasicPreds, ..CUDAGpuUtils,Cthulhu,BenchmarkTools , CUDA

goldBoolGPU= CuArray(trues(16,16,16));

segmBoolGPU= CuArray(trues(16,16,16));

fn = CuArray([0])

function kernelFunct(goldBoolGPU::CuDeviceArray{Bool, 3, 1}, segmBoolGPU::CuDeviceArray{Bool, 3, 1},fn)

    mainAfDiv, zz = fldmod1(blockIdx().x,8)

    i= (mainAfDiv-1) * blockDimX() + threadIdxX()

    j = (blockIdx().y-1) * blockDimY() + threadIdxY()

    z = (blockIdx().z-1) * blockDimZ() + threadIdxZ()

    CUDA.@cuprint """ i $(i)  j $(j) z $(z)   blockIdx().x $(blockIdx().x)  blockDimX()  $(blockDimX())   threadIdxX() $(threadIdxX())     blockIdx().y $(blockIdx().y)  blockDimY()  $(blockDimY())   threadIdxY() $(threadIdxY())  blockIdx().z $(blockIdx().z)  blockDimZ()  $(blockDimZ())   threadIdxZ() $(threadIdxZ()) \n """

    sync_threads()

    # if (goldBoolGPU[i,j,z] & !segmBoolGPU[i,j,z] )

    #    CUDA.@atomic fn[]+=1

    #     end

    return

    end

@cuda threads=(2, 2,2) blocks=(2,2,2)  kernelFunct(goldBoolGPU,segmBoolGPU,fn)

using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils,..ResultListUtils, ..Housdorff

using CUDA
dataBdim= (32,24,32)
fp = CUDA.zeros(UInt16,1)
sumInBits = (dataBdim[1]+2)+(dataBdim[2]+2)+(dataBdim[3]+2)+dataBdim[1]+dataBdim[2]+dataBdim[3]
shmemSum = cld(sumInBits,8)#in bytes
function testKernelA(dataBdim,fp)
    resShmem =  @cuDynamicSharedMem(Bool,((dataBdim[1]+2),(dataBdim[2]+2),(dataBdim[3]+2)))
       sourceShmem =  @cuDynamicSharedMem(Bool,(dataBdim[1],dataBdim[2],dataBdim[3]))
    #
    for i in 1:(dataBdim[1]+2),j in 1:(dataBdim[2]+2), n in 1:(dataBdim[3]+2)
        resShmem[i,j,n]=false
    end

    for i in 1:(dataBdim[1]),j in 1:(dataBdim[2]), n in 1:(dataBdim[3])
        sourceShmem[i,j,n]=false
    end
    fp[1]=1
return
end
@cuda threads=(32,5) blocks=(2) shmem=shmemSum  testKernelA(dataBdim,fp)
fp[1]

using CUDA
dataBdim= (32,24,32)
fp = CUDA.zeros(UInt16,1)
sumInBits = (dataBdim[1]+2)*(dataBdim[2]+2)*(dataBdim[3]+2)+dataBdim[1]*dataBdim[2]*dataBdim[3]
shmemSum = cld(sumInBits,8)#in bytes
function testKernelA(dataBdim,fp)
    resShmem =  @cuStaticSharedMem(Bool,(34,26,34))
    #sourceShmem =  @cuStaticSharedMem(Bool,(32,24,32))
    #
    for i in 1:(dataBdim[1]+2),j in 1:(dataBdim[2]+2), n in 1:(dataBdim[3]+2)
        resShmem[i,j,n]=false
    end

    # for i in 1:(dataBdim[1]),j in 1:(dataBdim[2]), n in 1:(dataBdim[3])
    #     sourceShmem[i,j,n]=false
    # end
    fp[1]=1
return
end
@cuda threads=(32,5) blocks=(2) shmem=shmemSum  testKernelA(dataBdim,fp)
fp[1]


text:
# select the two seeds maximizing separation, O(N)
function _splitseeds_linear(node::Node)
    @assert length(node) > 1

    seed1 = seed2 = 0
    max_separation = -Inf
    for n in 1:ndims(node)
        argmax_low = argmin_high = 1
        min_low = max_low = mbr(node[1]).low[n]
        min_high = max_high = mbr(node[1]).high[n]

        for i in 2:length(node)
            @inbounds lowi = mbr(node[i]).low[n]
            @inbounds highi = mbr(node[i]).high[n]
            if lowi > max_low
                max_low = lowi
                argmax_low = i
            elseif lowi < min_low
                min_low = lowi
            end
            if highi < min_high
                min_high = highi
                argmin_high = i
            elseif highi > max_high
                max_high = highi
            end
        end

        sepn = max_low - min_high
        if max_high > min_low
            sepn /= (max_high - min_low)
        end

        if sepn > max_separation
            max_separation = sepn
            seed1 = argmin_high
            seed2 = argmax_low
        end
    end

    @assert seed1 > 0 && seed2 > 0
    if seed1 == seed2 # avoid collisions
        if seed2 == 1
            seed2 += 1
        else
            seed2 -= 1
        end
    end
    return seed1, seed2
end

# select the two seeds maximizing inefficiency, O(N²)
function _splitseeds_quadratic(node::Node)
    @assert length(node) > 1

    seed1 = 0
    seed2 = 0
    inefficiency = -Inf
    # for each pair of Regions (account for overflow Region too!)
    for i in eachindex(node)
        @inbounds mbri = mbr(node[i])
        areai = area(mbri)

        for j in (i+1):length(node)
            # get the combined MBR of those two entries.
            @inbounds mbrj = mbr(node[j])
            @inbounds r = combine(mbri, mbrj)

            # find the inefficiency of grouping these entries together.
            d = area(r) - areai - area(mbrj)
            if d > inefficiency
                inefficiency = d;
                seed1 = i
                seed2 = j
            end
        end
    end
    return seed1, seed2
end

# pick the two seeds (children indices) for splitting the node into two
function _splitseeds(node::Node, tree::RTree)
    length(node) > 1 ||
        throw(SpatialIndexError("Cannot split the node with less than 2 children"))

    if variant(tree) == RTreeLinear || variant(tree) == RTreeStar
        return _splitseeds_linear(node)
    elseif variant(tree) == RTreeQuadratic
        return _splitseeds_quadratic(node)
    else
        throw(SpatialIndexException("RTree variant not supported"))
    end
end

# R-tree split node (select the node by minimal enlargment and area)
function _split!_rtree(node::Node, tree::RTree)
    nnodes_min = floor(Int, capacity(node, tree) * tree.fill_factor)
    #@debug "_split!_rtree(): lev=$(level(node)) len=$(length(node)) nmin=$nnodes_min"

    # initialize each group with the seed entries.
    seed1, seed2 = _splitseeds(node, tree)
    n1 = _attach!(acquire(tree, typeof(node), level(node)), node[seed1], tree)
    n2 = _attach!(acquire(tree, typeof(node), level(node)), node[seed2], tree)

    # use this mask for marking visited entries.
    available = trues(length(node))
    available[seed1] = available[seed2] = false
    navail = length(node) - 2
    while navail > 0
        if length(n1) + navail == nnodes_min || length(n2) + navail == nnodes_min
            # all remaining entries must be assigned to n1 or n2 to comply with minimum load requirement.
            sink = length(n1) + navail == nnodes_min ? n1 : n2
            i = findfirst(available)
            while i !== nothing
                _attach!(sink, node[i], tree)
                available[i] = false
                i = findnext(available, i+1)
            end
            navail = 0
        else
            # For all remaining entries compute the difference of the cost of grouping an
            # entry in either group. When done, choose the entry that yielded the maximum
            # difference. In case of linear split, select any entry (e.g. the first one.)
            area1 = area(mbr(n1))
            area2 = area(mbr(n2))

            max_cost_diff = -Inf
            sel_cost1 = sel_cost2 = NaN
            sel_child_ix = 0
            i = findfirst(available)
            while i !== nothing
                child_mbr = mbr(node[i])
                enl1 = area(combine(mbr(n1), child_mbr)) - area1
                enl2 = area(combine(mbr(n2), child_mbr)) - area2
                enl_delta = abs(enl1 - enl2)

                if enl_diff > max_enl_diff
                    max_cost_diff = enl_diff
                    sel_enl1 = enl1
                    sel_enl2 = enl2
                    sel_child_ix = i
                    (variant(tree) != RTreeQuadratic) && break
                end
                i = findnext(available, i+1)
            end

            # select the node where we should add the new entry.
            sink = sel_enl1 < sel_enl2 || (sel_enl1 == sel_enl2 &&
                  (area1 < area2 || (area1 == area2 &&
                   length(n1) <= length(n2)))) ? n1 : n2
            available[sel_child_ix] = false
            navail -= 1
            _attach!(sink, node[sel_child_ix], tree)
        end
    end
    @assert !any(available)

    return n1, n2
end

# get low/high of a specific MBR dimension
# accessor functions to speedup sortperm!()
_mbr_low(node, dim::Integer) = @inbounds mbr(node).low[dim]
_mbr_high(node, dim::Integer) = @inbounds mbr(node).high[dim]

# R*-star node split
function _split!_rstar(node::Node, tree::RTree)
    nsplit = floor(Int, length(node) * tree.split_factor)
    nsplit_distr = length(node) - 2 * nsplit + 1
    @assert nsplit_distr > 0

    #@debug "_split!_rstar(): lev=$(level(node)) len=$(length(node)) nsplit=$nsplit nsplit_distr=$nsplit_distr"

    # find split_dim that minimizes the perimeter in all splits
    min_perim = Inf
    split_dim = 0
    loworder = zeros(Int, length(node))
    highorder = similar(loworder)
    use_low = false
    for dim in 1:ndims(node)
        sortperm!(loworder, children(node), by=Base.Fix2(_mbr_low, dim))
        sortperm!(highorder, children(node), by=Base.Fix2(_mbr_high, dim))

        # calculate the sum of perimiters for all splits
        low_perim = 0.0
        high_perim = 0.0
        for i in nsplit:(nsplit + nsplit_distr)
            @inbounds br_low1 = mapreduce(j -> @inbounds(mbr(node[j])), combine, view(loworder, 1:i))
            @inbounds br_high1 = mapreduce(j -> @inbounds(mbr(node[j])), combine, view(highorder, 1:i))
            @inbounds br_low2 = mapreduce(j -> @inbounds(mbr(node[j])), combine, view(loworder, (i+1):length(node)))
            @inbounds br_high2 = mapreduce(j -> @inbounds(mbr(node[j])), combine, view(highorder, (i+1):length(node)))

            low_perim += perimeter(br_low1) + perimeter(br_low2)
            high_perim += perimeter(br_high1) + perimeter(br_high2)
        end

        perim = min(low_perim, high_perim)
        if perim < min_perim
            min_perim = perim
            split_dim = dim
            use_low = low_perim < high_perim
        end
    end

    # final sorting
    selorder = use_low ?
        sortperm!(loworder, children(node), by=Base.Fix2(_mbr_low, split_dim)) :
        sortperm!(highorder, children(node), by=Base.Fix2(_mbr_high, split_dim))

    # find the best split point (minimizes split overlap and area)
    min_overlap = Inf
    best_area_sum = Inf
    split_ix = 0
    for i in nsplit:(nsplit + nsplit_distr)
        @inbounds br1 = mapreduce(i -> mbr(node[selorder[i]]), combine, 1:i)
        @inbounds br2 = mapreduce(i -> mbr(node[selorder[i]]), combine, (i+1):length(node))

        overlap = overlap_area(br1, br2)
        area_sum = area(br1) + area(br2)
        if overlap < min_overlap || (overlap == min_overlap &&
           area_sum < best_area_sum)
            split_ix = i
            min_overlap = overlap
            best_area_sum = area_sum
        end
    end

    n1 = acquire(tree, typeof(node), level(node))
    for i in 1:split_ix
        _attach!(n1, node[selorder[i]], tree)
    end
    n2 = acquire(tree, typeof(node), level(node))
    for i in (split_ix+1):length(selorder)
        _attach!(n2, node[selorder[i]], tree)
    end
    return n1, n2
end

# split the node into two
function _split!(node::Node, tree::RTree)
    #@debug "_split!(): lev=$(level(node)) len=$(length(node))"

    if variant(tree) == RTreeLinear || variant(tree) == RTreeQuadratic
        return _split!_rtree(node, tree)
    elseif variant(tree) == RTreeStar
        return _split!_rstar(node, tree)
    else
        throw(SpatialIndexError("RTree variant not supported"))
    end
end