text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/UMAT.jl/blob/master/LICENSE

module UMAT

using Libdl
using Materials
using Parameters
using LinearAlgebra
using Tensors

const abaqus = [1 4 6; 0 2 5; 0 0 3]

pkg_dir = dirname(Base.find_package("UMAT"))
if Sys.iswindows()
    lib_dir = joinpath(pkg_dir,"..","deps","usr","bin")
else
    lib_dir = joinpath(pkg_dir,"..","deps","usr","lib")
end

"""
Variables updated by UMAT routine.
"""
@with_kw struct UmatVariableState <: AbstractMaterialState
    NTENS :: Int
    NSTATV :: Int = zero(Int)
    jacobian :: SymmetricTensor = zero(SymmetricTensor{4,convert(Int,NTENS/2),Float64})
    stress :: SymmetricTensor = zero(SymmetricTensor{2,convert(Int,NTENS/2),Float64})
    STATEV :: Array{Float64,1} = zeros(Float64, NSTATV)
    SSE :: Array{Float64,1} = zeros(Float64, 1)
    SPD :: Array{Float64,1} = zeros(Float64, 1)
    SCD :: Array{Float64,1} = zeros(Float64, 1)
    RPL :: Array{Float64,1} = zeros(Float64, 1)
    DDSDDT :: Array{Float64,1} = zeros(Float64, NTENS)
    DRPLDE :: Array{Float64,1} = zeros(Float64, NTENS)
    DRPLDT :: Array{Float64,1} = zeros(Float64, 1)
    PNEWDT :: Array{Float64,1} = ones(Float64, 1)
end

"""
Material parameters in order that is specific to chosen UMAT.
"""
@with_kw struct UmatParameterState <: AbstractMaterialState
    NPROPS :: Int = zero(Int)
    PROPS :: Array{Float64,1} = zeros(Float64, NPROPS)
end

"""
Variables passed in for information.
These drive evolution of the material state.
"""
@with_kw mutable struct UmatDriverState <: AbstractMaterialState
    NTENS :: Int
    strain :: SymmetricTensor = zero(SymmetricTensor{2,convert(Int,NTENS/2),Float64})
    time :: Float64 = zero(Float64)
    dtime :: Float64 = zero(Float64)
    TIME_ :: Array{Float64,1} = zeros(Float64, 2)
    TEMP :: Float64 = zero(Float64)
    PREDEF :: Float64 = zero(Float64)
end

"""
Other Abaqus UMAT variables passed in for information.
"""
@with_kw mutable struct UmatOtherState <: AbstractMaterialState
    CMNAME :: String = "U" # Ptr{Cuchar} # = 'U'
    NDI :: Int = 3
    NSHR :: Int = 3
    NTENS :: Int  # NDI + NSHR
    NSTATV :: Int
    NPROPS :: Int
    COORDS :: Array{Float64,1} = zeros(Float64, 3)
    DROT :: Array{Float64,2} = I + zeros(Float64, 3, 3)
    CELENT :: Float64 = zero(Float64)
    DFGRD0 :: Array{Float64,2} = I + zeros(Float64, 3, 3)
    DFGRD1 :: Array{Float64,2} = I + zeros(Float64, 3, 3)
    NOEL :: Int = zero(Int)
    NPT :: Int = zero(Int)
    LAYER :: Int = zero(Int)
    KSPT :: Int = zero(Int)
    JSTEP :: Array{Float64,1} = zeros(Int, 3)
    KSTEP :: Int = zero(Int)
    KINC :: Int = zero(Int)
end

"""
UMAT material structure.

`lib_path` is the path to the compiled shared library.
`behaviour` is the function name to call from the shared library
    generally this is compiler dependent, by default `umat_`.
    MFront Abaqus interface produces specific name, e.g. `ELASTICITY_3D`.
"""
@with_kw mutable struct UmatMaterial <: AbstractMaterial
    NTENS :: Int
    NSTATV :: Int = zero(Int)
    NPROPS :: Int = zero(Int)

    drivers :: UmatDriverState = UmatDriverState(NTENS=NTENS)
    ddrivers :: UmatDriverState = UmatDriverState(NTENS=NTENS)
    variables :: UmatVariableState = UmatVariableState(NTENS=NTENS, NSTATV=NSTATV)
    variables_new :: UmatVariableState = UmatVariableState(NTENS=NTENS, NSTATV=NSTATV)
    parameters :: UmatParameterState = UmatParameterState(NPROPS=NPROPS)
    dparameters :: UmatParameterState = UmatParameterState(NPROPS=NPROPS)

    umat_other :: UmatOtherState = UmatOtherState(NTENS=NTENS, NSTATV=NSTATV, NPROPS=NPROPS)

    lib_path :: String
    behaviour :: Symbol = :umat_

end

"""
Wrapper function to ccall the UMAT subroutine from the compiled shared library.
"""
function call_umat!(func_umat::Symbol, lib_path::String, STRESS_,STATEV,DDSDDE,
        SSE,SPD,SCD,RPL,DDSDDT,DRPLDE,DRPLDT,
        STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,NDI,NSHR,
        NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,DFGRD0,DFGRD1,
        NOEL,NPT,LAYER,KSPT,KSTEP,KINC)

    lib_umat = Libdl.dlopen(lib_path) # Open the library explicitly.
    sym_umat = Libdl.dlsym(lib_umat, func_umat) # Get a symbol for the umat function to call.

    ccall(sym_umat, Nothing,
         (Ref{Float64},Ref{Float64},Ref{Float64},Ref{Float64},Ref{Float64},
          Ref{Float64},Ref{Float64},Ref{Float64},Ref{Float64},Ref{Float64},
          Ref{Float64},Ref{Float64},Ref{Float64},Ref{Float64},Ref{Float64},
          Ref{Float64},Ref{Float64},Ref{Float64},Ptr{Cuchar},Ref{Int},Ref{Int},
          Ref{Int},Ref{Int},Ref{Float64},Ref{Int},Ref{Float64},Ref{Float64},
          Ref{Float64},Ref{Float64},Ref{Float64},Ref{Float64},
          Ref{Int},Ref{Int},Ref{Int},Ref{Int},Ref{Int},Ref{Int}),
         STRESS_,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,DRPLDE,DRPLDT,
         STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,NDI,NSHR,
         NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,DFGRD0,DFGRD1,
         NOEL,NPT,LAYER,KSPT,KSTEP,KINC)

    Libdl.dlclose(lib_umat) # Close the library explicitly.

    return Nothing
end

"""
Calls UMAT and updates MaterialVariableState writing the result to material.variables_new
"""
function Materials.integrate_material!(material::UmatMaterial)
    @unpack CMNAME,NDI,NSHR,NTENS,NSTATV,NPROPS,COORDS,DROT,CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC = material.umat_other
    @unpack PROPS = material.parameters
    strain, time, TEMP, PREDEF = material.drivers.strain, material.drivers.time, material.drivers.TEMP, material.drivers.PREDEF
    dstrain, dtime, DTEMP, DPRED = material.ddrivers.strain, material.drivers.dtime, material.ddrivers.TEMP, material.ddrivers.PREDEF
    @unpack stress,STATEV,jacobian,SSE,SPD,SCD,RPL,DDSDDT,DRPLDE,DRPLDT,PNEWDT = material.variables

    # TODO need to replace with toabaqus etc.
    STRAN = tovoigt(strain; order=abaqus, offdiagscale = 2.0)
    DSTRAN = tovoigt(dstrain; order=abaqus, offdiagscale = 2.0)
    STRESS_ = tovoigt(stress, order=abaqus)
    DDSDDE = tovoigt(jacobian, order=abaqus)
    TIME_ = [(time - floor(time)) time]
    DTIME_ = dtime
    call_umat!(material.behaviour, material.lib_path, STRESS_, STATEV, DDSDDE, SSE, SPD, SCD, RPL, DDSDDT, DRPLDE, DRPLDT,
        STRAN, DSTRAN, TIME_, DTIME_, TEMP, DTEMP, PREDEF, DPRED, CMNAME, NDI, NSHR,
        NTENS, NSTATV, PROPS, NPROPS, COORDS, DROT, PNEWDT, CELENT, DFGRD0, DFGRD1,
        NOEL, NPT, LAYER, KSPT, KSTEP, KINC)

    stress = fromvoigt(SymmetricTensor{2,convert(Int,NTENS/2)}, STRESS_, order=abaqus)
    jacobian = fromvoigt(SymmetricTensor{4,convert(Int,NTENS/2)}, DDSDDE, order=abaqus)

    variables_new = UmatVariableState(NSTATV=NSTATV,NTENS=NTENS,stress=stress,
                                      STATEV=STATEV,jacobian=jacobian,SSE=SSE,
                                      SPD=SPD,SCD=SCD,RPL=RPL,DDSDDT=DDSDDT,
                                      DRPLDE=DRPLDE,DRPLDT=DRPLDT,PNEWDT=PNEWDT)
    material.variables_new = variables_new
end

function Materials.reset_material!(material::UmatMaterial)
    material.ddrivers = UmatDriverState(NTENS=material.NTENS)
    material.dparameters = UmatParameterState(NPROPS=material.NPROPS)
    material.variables_new = UmatVariableState(NTENS=material.NTENS, NSTATV=material.NSTATV)
    return nothing
end

include("gurson_model.jl")
include("druckerprager_model.jl")

export UmatMaterial, UmatDriverState, UmatParameterState, UmatVariableState
export UmatOtherState, GursonMaterial, DruckerPragerMaterial, integrate_material!

# Re-export update_material! and uniaxial_increment! from Materials.jl for convenience
export update_material!, uniaxial_increment!
end # module


text:
"""Represent a technique for aligning two timeseries."""
abstract type Alignment end
struct LeftAlignment <: Alignment end
struct UnionAlignment <: Alignment end
struct IntersectAlignment <: Alignment end

"""
    LEFT

For inputs `(A, B, ...)`, tick whenever `A` ticks so long as all inputs are active.
"""
const LEFT = LeftAlignment()

"""
    UNION

For inputs `(A, B, ...)`, tick whenever any input ticks so long as all inputs are active.
"""
const UNION = UnionAlignment()

"""
    INTERSECT

For inputs `(A, B, ...)`, tick whenever all inputs tick simultaneously.
"""
const INTERSECT = IntersectAlignment()

"""The default alignment for operators when not specified."""
const DEFAULT_ALIGNMENT = UNION

"""
    UnaryNodeOp{T} <: NodeOp{T}

An abstract type representing a node op with a single parent.
"""
abstract type UnaryNodeOp{T} <: NodeOp{T} end

"""
    BinaryNodeOp{T,A<:Alignment} <: NodeOp{T}

An abstract type representing a node op with two parents, and using alignment `A`.
"""
abstract type BinaryNodeOp{T,A<:Alignment} <: NodeOp{T} end

"""
    NaryNodeOp{N,T,A<:Alignment} <: NodeOp{T}

An abstract type representing a node op with `N` parents, and using alignment `A`.

This type should be avoided for `N < 3`, since in these cases it would be more appropriate
to use either [`TimeDag.UnaryNodeOp`](@ref) or [`TimeDag.BinaryNodeOp`](@ref).
"""
abstract type NaryNodeOp{N,T,A<:Alignment} <: NodeOp{T} end

# A note on the design choice here, which is motivated by performance reasons & profiling.
#
# Options considered:
#   1. Return `(value::T, should_tick::Bool)` pair always
#   2. Take `out::Ref{T}` as a parameter, return `should_tick::Bool`
#   3. Return `Maybe{T}`.
#
# Option 1 was disfavoured, because if `should_tick` is false, one still needs to invent a
#   `value` of type `T`, which is possibly hard in general.
# Option 2 was found, from benchmarking, to be moderately slower than the other choices.
#
# Option 3 was found to retain the performance of option 1 with minimal overhead.
#
# For cases when we know the node will always tick (as indicated by `always_ticks(op)`) we
# just return a raw `T`, as the extra information will not be used.
"""
    operator!(op::UnaryNodeOp{T}, (state,), (time,) x) -> T / Maybe{T}
    operator!(op::BinaryNodeOp{T}, (state,), (time,) x, y) -> T / Maybe{T}
    operator!(op::NaryNodeOp{N,T}, (state,), (time,) x, y, z...) -> T / Maybe{T}

Perform the operation for this node.

When defining a method of this for a new op, follow these rules:
- `state` should be omitted iff [`TimeDag.stateless_operator`](@ref).
- `time` should be omitted iff [`TimeDag.time_agnostic`](@ref).
- All values `x, y, z...` should be omittted iff [`TimeDag.value_agnostic`](@ref).

For stateful operations, this operator should mutate `state` as required.

The return value `out` should be of type `T` iff [`TimeDag.always_ticks`](@ref) is true,
otherwise it should be of type [`TimeDag.Maybe`](@ref)`{T}`.

If `out <: Maybe{T}`, and has `!valid(out)`, this indicates that we do not wish to emit a
knot at this time, and it will be skipped. Otherwise, `value(out)` will be used as the
output value.
"""
function operator! end

"""
    always_ticks(node) -> Bool
    always_ticks(op) -> Bool

Returns true iff the return value from `operator` can be assumed to always be valid.

Note, that for sensible performance characteristics, this should be knowable from
`typeof(op)`
"""
always_ticks(node::Node) = always_ticks(node.op)
always_ticks(::NodeOp) = false

"""
    stateless_operator(node) -> Bool
    stateless_operator(op) -> Bool

Returns true iff `operator(op, ...)` would never look at or modify the evaluation state.

If this returns true, `create_operator_evaluation_state` will not be used.

Note that if an `op` has `stateless(op)` returning true, then it necessarily should also
return true here. The default implementation is to return `stateless(op)`, meaning that if
one is creating a node that is fully stateless, one need only define `stateless`.

For optimal performance, this should be knowable from the type of `op` alone.
"""
stateless_operator(node::Node) = stateless_operator(node.op)
stateless_operator(op::NodeOp) = stateless(op)

"""
    time_agnostic(node) -> Bool
    time_agnostic(op) -> Bool

Returns true iff `op` does not care about the time of the input knot(s).

For optimal performance, this should be knowable from the type of `op` alone.
"""
time_agnostic(node::Node) = time_agnostic(node.op)
time_agnostic(::NodeOp) = false

"""
    value_agnostic(node) -> Bool
    value_agnostic(op) -> Bool

Returns true iff `op` does not care about the value(s) of the input knot(s).

For optimal performance, this should be knowable from the type of `op` alone.
"""
value_agnostic(node::Node) = value_agnostic(node.op)
value_agnostic(::NodeOp) = false

"""
    create_operator_evaluation_state(parents, op::NodeOp) -> NodeEvaluationState

Create an empty evaluation state for the given node, when starting evaluation at the
specified time.

Note that this is state that will be passed to `operator`. The overall node may additionally
wrap this state with further state, if this is necessary for e.g. alignment.
"""
function create_operator_evaluation_state end

"""
    has_initial_values(op) -> Bool

If this returns true, it indicates that initial values for the `op`'s parents are specified.

See the documentation on [Initial values](@ref) for further information.
"""
has_initial_values(::BinaryNodeOp) = false
has_initial_values(::NaryNodeOp) = false

"""
    initial_left(op::BinaryNodeOp) -> L

Specify the initial value to use for the first parent of the given `op`.

Needs to be defined if [`has_initial_values`](@ref) returns true, and alignment is
[`UNION`](@ref). For other alignments it is not required.
"""
function initial_left end

"""
    initial_right(op::BinaryNodeOp) -> R

Specify the initial value to use for the second parent of the given `op`.

Needs to be defined if `has_initial_values(op)` returns true, and alignment is
[`UNION`](@ref) or [`LEFT`](@ref). For [`INTERSECT`](@ref) alignment it is not required.
"""
function initial_right end

"""
    initial_values(op::NaryNodeOp) -> Tuple

Specify the initial values to use for all parents of the given `op`.

Needs to be defined for nary ops for which [`has_initial_values`](@ref) returns true.
"""
function initial_values end

"""Convenience method to dispatch to reduced-argument `operator!` calls."""
function _operator!(op::NodeOp, state::NodeEvaluationState, time::DateTime, values...)
    return if stateless_operator(op) && time_agnostic(op)
        value_agnostic(op) ? operator!(op) : operator!(op, values...)
    elseif stateless_operator(op) && !time_agnostic(op)
        value_agnostic(op) ? operator!(op, time) : operator!(op, time, values...)
    elseif !stateless_operator(op) && time_agnostic(op)
        value_agnostic(op) ? operator!(op, state) : operator!(op, state, values...)
    else
        error("Error! We should have dispatched to a more specialised method.")
    end
end

function _can_propagate_constant(op::Union{UnaryNodeOp,BinaryNodeOp,NaryNodeOp})
    return always_ticks(op) && stateless_operator(op) && time_agnostic(op)
end
function _propagate_constant_value(op::UnaryNodeOp, parents::Tuple{Node})
    # NB, we know that time & state is ignored (due to _can_propagate_constant).
    return operator!(op, value(@inbounds(parents[1])))
end
function _propagate_constant_value(op::BinaryNodeOp, parents::Tuple{Node,Node})
    return operator!(op, value(@inbounds(parents[1])), value(@inbounds(parents[2])))
end
function _propagate_constant_value(op::NaryNodeOp{N}, parents::NTuple{N,Node}) where {N}
    return operator!(op, map(value, parents)...)
end

"""
    _create_operator_evaluation_state(parents, op) -> NodeEvaluationState

Internal function that will look at `stateless_operator`, and iff it is false call
`create_operator_evaluation_state`. Otherwise return an empty node state.
"""
function _create_operator_evaluation_state(parents, op)
    return if stateless_operator(op)
        EMPTY_NODE_STATE
    else
        create_operator_evaluation_state(parents, op)
    end
end

# An unary node has no alignment state, so any state comes purely from the operator.
function create_evaluation_state(parents::Tuple{Node}, op::UnaryNodeOp)
    return _create_operator_evaluation_state(parents, op)
end

function run_node!(
    node_op::UnaryNodeOp{T},
    state::NodeEvaluationState,
    ::DateTime,  # time_start
    ::DateTime,  # time_end
    input::Block{L},
) where {T,L}
    n = length(input)
    values = _allocate_values(T, n)
    return if always_ticks(node_op)
        # We can ignore the validity of the return value of the operator, since we have been
        # promised that it will always tick. Hence we can use a for loop too.
        @inbounds for i in 1:n
            time = input.times[i]
            values[i] = _operator!(node_op, state, time, input.values[i])
        end
        Block(unchecked, input.times, values)
    else
        times = _allocate_times(n)
        j = 1
        @inbounds for i in 1:n
            time = input.times[i]
            out = _operator!(node_op, state, time, input.values[i])
            if valid(out)
                values[j] = unsafe_value(out)
                times[j] = input.times[i]
                j += 1
            end
        end
        _trim!(times, j - 1)
        _trim!(values, j - 1)
        Block(unchecked, times, values)
    end
end

"""Apply, assuming `input_l` and `input_r` have identical alignment."""
function _apply_fast_align_binary!(
    op::BinaryNodeOp{T}, operator_state::NodeEvaluationState, input_l::Block, input_r::Block
) where {T}
    n = length(input_l)
    values = _allocate_values(T, n)
    return if always_ticks(op)
        # We shouldn't assume that it is valid to broadcast f over the inputs, so loop
        # manually.
        @inbounds for i in 1:n
            time = input_l.times[i]
            values[i] = _operator!(
                op, operator_state, time, input_l.values[i], input_r.values[i]
            )
        end
        Block(unchecked, input_l.times, values)
    else
        # FIXME Implement this branch!
        error("Not implemented!")
    end
end

abstract type UnionAlignmentState{L,R} <: NodeEvaluationState end

mutable struct UnionWithOpState{L,R,OperatorState} <: UnionAlignmentState{L,R}
    valid_l::Bool
    valid_r::Bool
    operator_state::OperatorState
    # These fields may initially be uninitialised.
    latest_l::L
    latest_r::R

    function UnionWithOpState{L,R}(operator_state::OperatorState) where {L,R,OperatorState}
        return new{L,R,OperatorState}(false, false, operator_state)
    end

    function UnionWithOpState{L,R}(
        operator_state::OperatorState, initial_l::L, initial_r::R
    ) where {L,R,OperatorState}
        return new{L,R,OperatorState}(true, true, operator_state, initial_l, initial_r)
    end
end

mutable struct UnionWithoutOpState{L,R} <: UnionAlignmentState{L,R}
    valid_l::Bool
    valid_r::Bool
    # These fields may initially be uninitialised.
    latest_l::L
    latest_r::R

    function UnionWithoutOpState{L,R}() where {L,R}
        return new{L,R}(false, false)
    end

    function UnionWithoutOpState{L,R}(initial_l::L, initial_r::R) where {L,R}
        return new{L,R}(true, true, initial_l, initial_r)
    end
end

operator_state(state::UnionWithOpState) = state.operator_state
operator_state(::UnionWithoutOpState) = EMPTY_NODE_STATE

function _set_l!(state::UnionAlignmentState{L}, x::L) where {L}
    state.latest_l = x
    state.valid_l = true
    return state
end

function _set_r!(state::UnionAlignmentState{L,R}, x::R) where {L,R}
    state.latest_r = x
    state.valid_r = true
    return state
end

function create_evaluation_state(
    parents::Tuple{Node,Node}, op::BinaryNodeOp{T,UnionAlignment}
) where {T}
    L = value_type(parents[1])
    R = value_type(parents[2])
    return if stateless_operator(op)
        if has_initial_values(op)
            UnionWithoutOpState{L,R}(initial_left(op), initial_right(op))
        else
            UnionWithoutOpState{L,R}()
        end
    else
        operator_state = create_operator_evaluation_state(parents, op)
        if has_initial_values(op)
            UnionWithOpState{L,R}(operator_state, initial_left(op), initial_right(op))
        else
            UnionWithOpState{L,R}(operator_state)
        end
    end
end

@inline function _maybe_add_knot!(
    node_op::NodeOp,
    operator_state::NodeEvaluationState,
    out_times::AbstractVector{DateTime},
    out_values::AbstractVector{T},
    j::Int,
    time::DateTime,
    in_values...,
) where {T}
    # Find the output value. For a given op this will either be of type T, or Maybe{T}, and
    # we can (at compile time) select the correct branch below based on `always_ticks(op)`.
    out = _operator!(node_op, operator_state, time, in_values...)

    return if always_ticks(node_op)
        # Output value is raw, and should always be used.
        @inbounds out_times[j] = time
        @inbounds out_values[j] = out
        j + 1
    else
        if valid(out)
            # Add the output only if the output is valid.
            @inbounds out_times[j] = time
            @inbounds out_values[j] = unsafe_value(out)
            j + 1
        else
            j
        end
    end
end

function run_node!(
    node_op::BinaryNodeOp{T,UnionAlignment},
    state::UnionAlignmentState,
    ::DateTime,  # time_start
    ::DateTime,  # time_end
    input_l::Block,
    input_r::Block,
) where {T}
    @inbounds if isempty(input_l) && isempty(input_r)
        # Nothing to do, since neither input has ticked.
        return Block{T}()
    elseif isempty(input_l) && !state.valid_l
        # Left is inactive and won't tick, so nothing gets emitted. But make sure we update
        # the state on the right.
        _set_r!(state, last(input_r.values))
        return Block{T}()
    elseif isempty(input_r) && !state.valid_r
        # Right is inactive and won't tick, so nothing gets emitted. But make sure we update
        # the state on the left.
        _set_l!(state, last(input_l.values))
        return Block{T}()
    end

    if _equal_times(input_l, input_r)
        # Times are indistinguishable
        # Update the alignment state.
        _set_l!(state, @inbounds last(input_l.values))
        _set_r!(state, @inbounds last(input_r.values))
        return _apply_fast_align_binary!(node_op, operator_state(state), input_l, input_r)
    end

    # Create our outputs as the maximum possible size.
    nl = length(input_l)
    nr = length(input_r)
    max_size = nl + nr
    times = _allocate_times(max_size)
    values = _allocate_values(T, max_size)

    # Store indices into the inputs. The index represents the next time point for
    # consideration for each series.
    il = 1
    ir = 1

    # Index into the output.
    j = 1

    # Loop until we exhaust inputs.
    @inbounds while (il <= nl || ir <= nr)
        # Store the next available time in the series, that is being pointed to by il & ir.
        next_time_l = il <= nl ? input_l.times[il] : DateTime(0)
        next_time_r = ir <= nr ? input_r.times[ir] : DateTime(0)

        new_time = if (il <= nl && next_time_l < next_time_r) || ir > nr
            # Left ticks next
            _set_l!(state, input_l.values[il])
            il += 1
            next_time_l
        elseif (ir <= nr && next_time_r < next_time_l) || il > nl
            # Right ticks next
            _set_r!(state, input_r.values[ir])
            ir += 1
            next_time_r
        else
            # A shared time point where neither x1 nor x2 have been exhausted.
            _set_l!(state, input_l.values[il])
            _set_r!(state, input_r.values[ir])
            il += 1
            ir += 1
            next_time_l
        end

        # We must only output a knot if both inputs are active.
        if !state.valid_l || !state.valid_r
            continue
        end

        # Compute and possibly output a knot.
        j = _maybe_add_knot!(
            node_op,
            operator_state(state),
            times,
            values,
            j,
            new_time,
            state.latest_l,
            state.latest_r,
        )
    end

    # Package the outputs into a block, resizing the outputs as necessary.
    _trim!(times, j - 1)
    _trim!(values, j - 1)
    return Block(unchecked, times, values)
end

function create_evaluation_state(
    parents::Tuple{Node,Node}, op::BinaryNodeOp{T,IntersectAlignment}
) where {T}
    # Intersect alignment doesn't require remembering any previous state, so just return
    # the operator state.
    return _create_operator_evaluation_state(parents, op)
end

function run_node!(
    node_op::BinaryNodeOp{T,IntersectAlignment},
    operator_state::NodeEvaluationState,
    ::DateTime,  # time_start
    ::DateTime,  # time_end
    input_l::Block,
    input_r::Block,
) where {T}
    if isempty(input_l) || isempty(input_r)
        # Output will be empty unless both inputs have ticked.
        return Block{T}()
    end

    if _equal_times(input_l, input_r)
        # Times are indistinguishable.
        return _apply_fast_align_binary!(node_op, operator_state, input_l, input_r)
    end

    # Create our outputs as the maximum possible size.
    nl = length(input_l)
    nr = length(input_r)
    max_size = min(nl, nr)
    times = _allocate_times(max_size)
    values = _allocate_values(T, max_size)

    # Store indices into the inputs. The index represents the next time point for
    # consideration for each series.
    il = 1
    ir = 1

    # Index into the output.
    j = 1

    # If we get to the end of either series, we know that we cannot add any more elements to
    # the output.
    @inbounds while (il <= nl && ir <= nr)
        # Obtain the *next available* times from each entity. We know that the current
        # state, and last emitted, time is strictly less than either of these.
        next_time_l = input_l.times[il]
        next_time_r = input_r.times[ir]

        if next_time_l < next_time_r
            # Left ticks next; consider the next knot.
            il += 1
        elseif next_time_r < next_time_l
            # Right ticks next; consider the next knot.
            ir += 1
        else  # next_time_l == next_time_r
            # Shared time point, so emit a knot.
            j = _maybe_add_knot!(
                node_op,
                operator_state,
                times,
                values,
                j,
                next_time_l,
                input_l.values[il],
                input_r.values[ir],
            )
            il += 1
            ir += 1
        end
    end

    # Package the outputs into a block, resizing the outputs as necessary.
    _trim!(times, j - 1)
    _trim!(values, j - 1)
    return Block(unchecked, times, values)
end

abstract type LeftAlignmentState{R} <: NodeEvaluationState end

mutable struct LeftWithOpState{R,OperatorState} <: LeftAlignmentState{R}
    valid_r::Bool
    operator_state::OperatorState
    latest_r::R

    function LeftWithOpState{R}(operator_state::OperatorState) where {R,OperatorState}
        return new{R,OperatorState}(false, operator_state)
    end

    function LeftWithOpState{R}(
        operator_state::OperatorState, initial_r::R
    ) where {R,OperatorState}
        return new{R,OperatorState}(true, operator_state, initial_r)
    end
end

mutable struct LeftWithoutOpState{R} <: LeftAlignmentState{R}
    valid_r::Bool
    latest_r::R

    LeftWithoutOpState{R}() where {R} = new{R}(false)
    LeftWithoutOpState{R}(initial_r::R) where {R} = new{R}(true, initial_r)
end

operator_state(::LeftWithoutOpState) = EMPTY_NODE_STATE
operator_state(state::LeftWithOpState) = state.operator_state

function _set_r!(state::LeftAlignmentState{R}, x::R) where {R}
    state.latest_r = x
    state.valid_r = true
    return state
end

function create_evaluation_state(
    parents::Tuple{Node,Node}, op::BinaryNodeOp{T,LeftAlignment}
) where {T}
    R = value_type(parents[2])
    return if stateless_operator(op)
        if has_initial_values(op)
            LeftWithoutOpState{R}(initial_right(op))
        else
            LeftWithoutOpState{R}()
        end
    else
        operator_state = create_operator_evaluation_state(parents, op)
        if has_initial_values(op)
            LeftWithOpState{R}(operator_state, initial_right(op))
        else
            LeftWithOpState{R}(operator_state)
        end
    end
end

function run_node!(
    node_op::BinaryNodeOp{T,LeftAlignment},
    state::LeftAlignmentState,
    ::DateTime,  # time_start
    ::DateTime,  # time_end
    input_l::Block,
    input_r::Block,
) where {T}
    have_initial_r = state.valid_r

    if isempty(input_l)
        # We will not tick, but update state if necessary.
        if !isempty(input_r)
            _set_r!(state, @inbounds last(input_r.values))
        end
        return Block{T}()
    elseif isempty(input_r) && !have_initial_r
        # We cannot tick, since we have no values on the right. No state to update either.
        return Block{T}()
    end

    if _equal_times(input_l, input_r)
        # Times are indistinguishable.
        return _apply_fast_align_binary!(node_op, operator_state(state), input_l, input_r)
    end

    # The most we can emit is one knot for every knot in input_l.
    nl = length(input_l)
    nr = length(input_r)
    times = _allocate_times(nl)
    values = _allocate_values(T, nl)

    # Start with 0, indicating that input_r hasn't started ticking yet.
    ir = 0

    # The index into the output.
    j = 1

    @inbounds for il in 1:nl
        # Consume r while it would leave us before the current time in l, or until we reach
        # the end of r.
        next_time_l = input_l.times[il]
        while (ir < nr && input_r.times[ir + 1] <= next_time_l)
            ir += 1
        end

        if ir > 0
            j = _maybe_add_knot!(
                node_op,
                operator_state(state),
                times,
                values,
                j,
                next_time_l,
                input_l.values[il],
                input_r.values[ir],
            )

        elseif have_initial_r
            # R hasn't ticked in this batch, but we have an initial value.
            j = _maybe_add_knot!(
                node_op,
                operator_state(state),
                times,
                values,
                j,
                next_time_l,
                input_l.values[il],
                state.latest_r,
            )
        end
    end

    # Update state
    if !isempty(input_r)
        _set_r!(state, @inbounds last(input_r.values))
    end

    # Package the outputs into a block, resizing the outputs as necessary.
    _trim!(times, j - 1)
    _trim!(values, j - 1)
    return Block(unchecked, times, values)
end

# Some thought went into how to store `latest` for the nary case.
#
# Ideally one would use a mutable heterogeneous collection (i.e. a "mutable tuple"), but
# such a thing doesn't exist in Julia. There are tricks one can use if your elements are
# isbits (see StaticArrays.MArray), but we want a solution that can work for *any* input
# Julia type.
#
# Options considered:
#   1. Nary alignment, with latest::Vector{Any} and accept that state *cannot* store latest
#       values efficiently (will be a lot of runtime type checking). This WILL be slow, as
#       it is in the inner loop of evaluation.
#   2. Nary alignment, and use a Tuple for storage.
#       But, because this isn't mutable, will have to use `Base.setindex` to update, which
#       could theoretically do a lot of allocating.
#   3. Use code generation to make Union2, Union3, Union4, ...
#       This bloats the code somewhat, but would be guaranteed to be fast.
#       However, there is then an arbitrary limit to the number of things that we can align.
#
# Benchmarking suggests that option 2 can actually be quite efficient in practice, since if
# we immediately assign the tuple, LLVM has the ability to re-use memory. We hence take this
# approach.
#
# In the future, we could generate a SMALL number of additional special-cases using option 3
# (if, for example, we found ourselves doing a lot of ternary alignment).
#
# Also, to simplify matters, we use a single alignment state. It is slightly less memory
# efficient, because we store a pointer for operator state even if it is unnecessary.

# `Types` will look something like Tuple{In1,In2,...}
mutable struct NaryAlignmentState{N,Types<:Tuple,OperatorState} <: NodeEvaluationState
    # Since we may not have a latest value, we use the partial initialisation trick inside
    # `Maybe` to avoid having to invent an unused placeholder for arbitrary types.
    latest::Tuple{Vararg{Maybe}}
    operator_state::OperatorState

    function NaryAlignmentState{N,Types}(
        operator_state::OperatorState
    ) where {N,Types,OperatorState}
        return new{N,Types,OperatorState}(
            Tuple(Maybe{T}() for T in Types.parameters), operator_state
        )
    end

    function NaryAlignmentState{N,Types}(
        operator_state::OperatorState, initial_values::Types
    ) where {N,Types,OperatorState}
        return new{N,Types,OperatorState}(
            Tuple(Maybe(v) for v in initial_values), operator_state
        )
    end
end

operator_state(state::NaryAlignmentState) = state.operator_state

Base.@propagate_inbounds function _set!(state::NaryAlignmentState, x, i::Integer)
    # TODO Consider making a mutable Maybe type, to avoid possible reconstruction of the
    #   tuple? It might just end up being slower though due to extra dereferencing.
    # The state tuple is not mutable, so we need to use this trick.
    state.latest = Base.setindex(state.latest, Maybe(x), i)
    return state
end

function create_evaluation_state(
    parents::NTuple{N,Node}, op::NaryNodeOp{N,T,UnionAlignment}
) where {N,T}
    # Work out the tuple of input types necessary for constructing the alignment state.
    Types = Tuple{map(value_type, parents)...}
    operator_state = _create_operator_evaluation_state(parents, op)
    return if has_initial_values(op)
        NaryAlignmentState{N,Types}(operator_state, initial_values(op))
    else
        NaryAlignmentState{N,Types}(operator_state)
    end
end

function create_evaluation_state(
    parents::NTuple{N,Node}, op::NaryNodeOp{N,T,IntersectAlignment}
) where {N,T}
    # Intersect alignment doesn't require remembering any previous state, so just return
    # the operator state.
    return _create_operator_evaluation_state(parents, op)
end

function create_evaluation_state(
    parents::NTuple{N,Node}, op::NaryNodeOp{N,T,LeftAlignment}
) where {N,T}
    # Work out the tuple of input types necessary for constructing the alignment state.
    # Note that we should ignore the left parent, so we only capture the rightmost M = N-1.
    _, parents_r... = parents
    M = N - 1
    Types = Tuple{map(value_type, parents_r)...}
    operator_state = _create_operator_evaluation_state(parents, op)
    return if has_initial_values(op)
        _, initial_values_r... = initial_values(op)
        NaryAlignmentState{M,Types}(operator_state, initial_values_r)
    else
        NaryAlignmentState{M,Types}(operator_state)
    end
end

"""Apply, assuming all `inputs` have identical alignment."""
function _apply_fast_align_nary!(
    op::NaryNodeOp{N,T}, operator_state::NodeEvaluationState, inputs::NTuple{N,Block}
) where {N,T}
    times = @inbounds first(inputs).times
    n = length(times)
    values = _allocate_values(T, n)
    return if always_ticks(op)
        @inbounds for i in 1:n
            time = times[i]
            values[i] = _operator!(
                op, operator_state, time, map(x -> x.values[i], inputs)...
            )
        end
        Block(unchecked, times, values)
    else
        # FIXME Implement this branch!
        error("Not implemented!")
    end
end

"""
    equivalence_classes(f, collection) -> Vector{Vector{eltype(collection)}}

Generate the set of equivalence classes for `collection` generated by the equivalence
relation `f : T × T → Bool`.

Note that behaviour is undefined if `f` is non-transitive.
"""
function equivalence_classes(
    f::Function, x::Union{AbstractVector{T},Tuple{Vararg{T}}}
) where {T}
    result = Vector{Vector{T}}()
    isempty(x) && return result

    @inbounds for el in x
        for class in result
            f(el, first(class)) || continue
            push!(class, el)
            # Jump to "found" label below. This means that we don't need to add a new
            # equivalence class.
            @goto found
        end

        # We end up here if we didn't find an existing equivalence class for `el`, so we
        # need to create a new one.
        push!(result, T[el])

        # Label used for escaping from the for loop above.
        @label found
    end
    return result
end

function run_node!(
    node_op::NaryNodeOp{N,T,UnionAlignment},
    state::NaryAlignmentState{N},
    ::DateTime,  # time_start
    ::DateTime,  # time_end
    inputs::Block...,
) where {N,T}
    @assert N == length(inputs)

    inputs_empty = map(isempty, inputs)
    @inbounds if all(inputs_empty)
        # Nothing to do, as no inputs have ticked.
        return Block{T}()
    elseif any(inputs_empty .& map(!valid, state.latest))
        # There is at least one input which has an empty input and does *not* have an
        # initial value.

        # We need to make sure we update any states which do have inputs.
        for (k, input) in enumerate(inputs)
            isempty(input) && continue
            _set!(state, last(input.values), k)
        end

        # The output will, however, be empty.
        return Block{T}()
    end

    input_classes = equivalence_classes(_equal_times, inputs)

    if length(input_classes) == 1
        # All times are indistinguishable.
        # Update alignment state and use fast alignment.
        @inbounds for (k, input) in enumerate(inputs)
            _set!(state, last(input.values), k)
        end
        return _apply_fast_align_nary!(node_op, operator_state(state), inputs)
    end

    # TODO `max_size` could be a massive overestimate.
    #   In practice it may be better to pick something smaller, and then increase the buffer
    #   size where necessary.

    # We use an optimisation here, whereby we use the `_equal_times` equivalence
    # relation to partition the inputs into sets which have equal times. We know that times
    # that are equal can only appear once, so this avoids some double counting.
    max_size = sum(x -> length(first(x)), input_classes)
    times = _allocate_times(max_size)
    values = _allocate_values(T, max_size)

    # Indices into the inputs. The index represents the next time point for
    # consideration for each series.
    is = MVector{N,Int64}(ones(Int64, N))

    # Index into the output.
    j = 1

    ns = map(length, inputs)
    @inbounds while true
        unfinisheds = (is .<= ns)
        any(unfinisheds) || break  # Loop until we exhaust inputs.

        # For each input, figure out the next time it would tick.
        next_times = ntuple(Val(N)) do k
            # If we're at the end of the series, just use a placeholder.
            !unfinisheds[k] && return typemax(DateTime)
            return inputs[k].times[is[k]]
        end

        # This is when we will next tick.
        new_time = minimum(next_times)

        # For every input that ticks at exactly this time, update to use the new value.
        for k in 1:N
            unfinisheds[k] || continue
            next_times[k] == new_time || continue
            # If we are here, this input should advance.
            _set!(state, inputs[k].values[is[k]], k)
            is[k] += 1
        end

        # We must only output a knot if all inputs are active.
        all(map(valid, state.latest)) || continue

        # Compute and possibly output a knot.
        j = _maybe_add_knot!(
            node_op,
            operator_state(state),
            times,
            values,
            j,
            new_time,
            map(unsafe_value, state.latest)...,
        )
    end

    # Package the outputs into a block, resizing the outputs as necessary.
    _trim!(times, j - 1)
    _trim!(values, j - 1)

    return Block(unchecked, times, values)
end

function run_node!(
    node_op::NaryNodeOp{N,T,IntersectAlignment},
    operator_state::NodeEvaluationState,
    ::DateTime,  # time_start
    ::DateTime,  # time_end
    inputs::Block...,
) where {N,T}
    @assert N == length(inputs)

    if any(isempty, inputs)
        # Nothing to do, as no inputs have ticked.
        return Block{T}()
    end

    input_classes = equivalence_classes(_equal_times, inputs)

    if length(input_classes) == 1
        # All times are indistinguishable, so use fast alignment.
        return _apply_fast_align_nary!(node_op, operator_state, inputs)
    end

    # Create our output as the maximum possible size.
    ns = map(length, inputs)
    max_size = minimum(ns)
    times = _allocate_times(max_size)
    values = _allocate_values(T, max_size)

    # Store indices into the inputs. The index represents the next time point for
    # consideration for each series.
    is = MVector{N,Int64}(ones(Int64, N))

    # Index into the output.
    j = 1

    # If we get to the end of any series, we know that we cannot add any more elements to
    # the output.
    @inbounds while all(is .<= ns)
        # For each input, figure out the next time it would tick.
        next_times = ntuple(k -> inputs[k].times[is[k]], Val(N))

        # This is when we *might* next tick, if all inputs tick at this time.
        new_time = minimum(next_times)

        tick_mask = next_times .== new_time
        if all(tick_mask)
            # All inputs tick simultaneously, this means we should compute and possibly
            # output a knot.
            input_values = map(zip(is, inputs)) do pair
                i, input = pair
                return input.values[i]
            end
            j = _maybe_add_knot!(
                node_op, operator_state, times, values, j, new_time, input_values...
            )
            # Advance all pointers.
            is .+= 1
        else
            # We aren't going to emit a knot, so advance all inputs that tick next.
            for k in 1:N
                next_times[k] == new_time || continue
                # If we are here, this input should advance.
                is[k] += 1
            end
        end
    end

    # Package the outputs into a block, resizing the outputs as necessary.
    _trim!(times, j - 1)
    _trim!(values, j - 1)
    return Block(unchecked, times, values)
end

function run_node!(
    node_op::NaryNodeOp{N,T,LeftAlignment},
    state::NaryAlignmentState{M},
    ::DateTime,  # time_start
    ::DateTime,  # time_end
    input_l::Block,
    inputs_r::Block...,
) where {N,T,M}
    @assert M == length(inputs_r)
    @assert M == N - 1

    wont_emit = @inbounds (
        isempty(input_l) ||  # Left input empty
        (
            # at least one right input inactive
            any(zip(inputs_r, state.latest)) do pair
                input_r, latest = pair
                return (isempty(input_r) && !valid(latest))
            end
        )
    )

    @inbounds if wont_emit
        # We will not be emitting a knot, but we should go through and update the state
        # where necessary
        for (kr, input_r) in enumerate(inputs_r)
            isempty(input_r) && continue
            _set!(state, last(input_r.values), kr)
        end
        return Block{T}()
    end

    all_times_equal = all(input_r -> _equal_times(input_l, input_r), inputs_r)
    if all_times_equal
        # All times are indistinguishable.
        # Update alignment state and use fast alignment.
        @inbounds for (k, input) in enumerate(inputs_r)
            _set!(state, last(input.values), k)
        end
        inputs = (input_l, inputs_r...)
        return _apply_fast_align_nary!(node_op, operator_state(state), inputs)
    end

    # The most we can emit is one knot for every knot in input_l.
    nl = length(input_l)
    times = _allocate_times(nl)
    values = _allocate_values(T, nl)

    # Indices into the right inputs. The index represents the next time point for
    # consideration for each series.
    # Use zero as the intitial index to indicate that the corresponding input hasn't
    # started ticking yet.
    irs = MVector{M,Int64}(zeros(Int64, M))

    # Index into the output.
    j = 1

    # Length of each 'right' input.
    nrs = map(length, inputs_r)

    @inbounds for il in 1:nl
        new_time = input_l.times[il]

        # Consume right inputs until we get to at or before `new_time`, or until we reach
        # the end of the input.
        for kr in 1:M
            nr = nrs[kr]
            nr < 1 && continue  # Nothing to do in this case.
            input_r = inputs_r[kr]
            advanced = false
            while (irs[kr] < nr && input_r.times[irs[kr] + 1] <= new_time)
                irs[kr] += 1
                advanced = true
            end

            # If this input has advanced, we should update the state.
            if advanced
                _set!(state, input_r.values[irs[kr]], kr)
            end
        end

        # Ascertain whether we can emit a knot. The answer is "yes" so long as all
        # right-states are valid.
        if all(valid, state.latest)
            j = _maybe_add_knot!(
                node_op,
                operator_state(state),
                times,
                values,
                j,
                new_time,
                input_l.values[il],
                map(unsafe_value, state.latest)...,
            )
        end
    end

    # We now need to update the state for the right-inputs, as we may not have exhausted
    # them all in the loop above.
    @inbounds for (k, input) in enumerate(inputs_r)
        _set!(state, last(input.values), k)
    end

    # Package the outputs into a block, resizing the outputs as necessary.
    _trim!(times, j - 1)
    _trim!(values, j - 1)
    return Block(unchecked, times, values)
end

ext: jl


text:
mutable struct nccpop_data
    n # number of variables
    m # number of constraints
    numeq # number of equality constraints
    supp # support data
    coe # coefficient data
    obj # "eigen" or "trace"
    basis # monomial bases
    ksupp # extending support at the k-th step
    sb # sizes of different blocks
    numb # numbers of different blocks
    blocks # block structure
    cl # numbers of blocks
    blocksize # sizes of blocks
end

"""
    opt,data = nctssos_first(pop::Vector{Polynomial{false, T}} where T<:Number, x::Vector{PolyVar{false}},
        d::Int; numeq=0, reducebasis=false, TS="block", obj="eigen", merge=false, md=3, solve=true, QUIET=false)

Compute the first step of the NCTSSOS hierarchy for constrained noncommutative polynomial optimization with
relaxation order `d`.
Return the optimum and other auxiliary data.

# Arguments
- `pop`: the vector of the objective function, inequality constraints, and equality constraints.
- `x`: the set of noncommuting variables.
- `d`: the relaxation order of the moment-SOHS hierarchy.
- `numeq`: the number of equality constraints.
"""

function nctssos_first(pop::Vector{Polynomial{false, T}} where T<:Number, x::Vector{PolyVar{false}},
    d::Int; numeq=0, reducebasis=false, TS="block", obj="eigen", merge=false, md=3, solve=true, QUIET=false)
    n,supp,coe = polys_info(pop, x)
    opt,data = nctssos_first(supp, coe, n, d, numeq=numeq, reducebasis=reducebasis, TS=TS, obj=obj, merge=merge, md=md, QUIET=QUIET, solve=solve)
    return opt,data
end

function polys_info(pop, x)
    n = length(x)
    m = length(pop)-1
    coe = Vector{Vector{Float64}}(undef, m+1)
    supp = Vector{Vector{Vector{UInt16}}}(undef, m+1)
    for k = 1:m+1
        mon = monomials(pop[k])
        coe[k] = coefficients(pop[k])
        supp[k] = [UInt16[] for i=1:length(mon)]
        for i = 1:length(mon)
            ind = mon[i].z .> 0
            vars = mon[i].vars[ind]
            exp = mon[i].z[ind]
            for j = 1:length(vars)
                l = ncbfind(x, n, vars[j], rev=true)
                append!(supp[k][i], l*ones(UInt16, exp[j]))
            end
        end
    end
    return n,supp,coe
end

function nctssos_first(supp::Vector{Vector{Vector{UInt16}}}, coe, n::Int64, d::Int64; numeq=0,
    reducebasis=false, TS="block", obj="eigen", merge=false, md=3, solve=true, QUIET=false)
    println("***************************NCTSSOS***************************")
    println("NCTSSOS is launching...")
    m = length(supp)-1
    dg = [maximum(length.(supp[i])) for i=2:m+1]
    if obj == "trace"
        supp[1],coe[1] = cyclic_canon(supp[1], coe[1])
    else
        supp[1],coe[1] = sym_canon(supp[1], coe[1])
    end
    basis = Vector{Vector{Vector{UInt16}}}(undef, m+1)
    basis[1] = get_ncbasis(n,d)
    ksupp = copy(supp[1])
    for i = 1:m
        basis[i+1] = get_ncbasis(n, d-Int(ceil(dg[i]/2)))
        if obj == "trace"
            append!(ksupp, [min(_cyclic_canon(word), _cyclic_canon(reverse(word))) for word in supp[i+1]])
        else
            append!(ksupp, _sym_canon.(supp[i+1]))
        end
    end
    if obj == "trace"
        append!(ksupp, [_cyclic_canon([basis[1][i][end:-1:1]; basis[1][i]]) for i=1:length(basis[1])])
    else
        append!(ksupp, [[basis[1][i][end:-1:1]; basis[1][i]] for i=1:length(basis[1])])
    end
    sort!(ksupp)
    unique!(ksupp)
    if TS != false && QUIET == false
        println("Starting to compute the block structure...")
    end
    time = @elapsed begin
    blocks,cl,blocksize,sb,numb,_ = get_nccblocks(m, ksupp, supp[2:end], basis, TS=TS, obj=obj, QUIET=QUIET, merge=merge, md=md)
    if reducebasis == true && obj == "eigen"
        gsupp = get_ncgsupp(m, supp, basis[2:end], blocks[2:end], cl[2:end], blocksize[2:end])
        psupp = copy(supp[1])
        push!(psupp, UInt16[])
        append!(psupp, gsupp)
        psupp = psupp[is_sym.(psupp)]
        basis[1],flag = reducebasis!(psupp, basis[1], blocks[1], cl[1], blocksize[1])
        if flag == 1
            ksupp = copy(supp[1])
            for i = 1:m
                append!(ksupp, _sym_canon.(supp[i+1]))
            end
            append!(ksupp, [[basis[1][i][end:-1:1]; basis[1][i]] for i=1:length(basis[1])])
            sort!(ksupp)
            unique!(ksupp)
            blocks,cl,blocksize,sb,numb,_ = get_nccblocks(m, ksupp, supp[2:end], basis, TS=TS, obj=obj, QUIET=QUIET, merge=merge, md=md)
        end
    end
    end
    if TS != false && QUIET == false
        mb = maximum(maximum.(sb))
        println("Obtained the block structure in $time seconds. The maximal size of blocks is $mb.")
    end
    opt,ksupp = ncblockcpop(m, supp, coe, basis, blocks, cl, blocksize, numeq=numeq, QUIET=QUIET, obj=obj, solve=solve)
    data = nccpop_data(n, m, numeq, supp, coe, obj, basis, ksupp, sb, numb, blocks, cl, blocksize)
    return opt,data
end

function nctssos_higher!(data::nccpop_data; TS="block", merge=false, md=3, solve=true, QUIET=false)
    m = data.m
    numeq = data.numeq
    supp = data.supp
    coe = data.coe
    obj = data.obj
    basis = data.basis
    ksupp = data.ksupp
    sb = data.sb
    numb = data.numb
    blocks = data.blocks
    cl = data.cl
    blocksize = data.blocksize
    if QUIET == false
        println("Starting to compute the block structure...")
    end
    time = @elapsed begin
    blocks, cl, blocksize, sb, numb, status = get_nccblocks(m, ksupp, supp[2:end], basis, blocks=blocks, cl=cl, blocksize=blocksize, sb=sb, numb=numb, TS=TS, obj=obj, QUIET=QUIET, merge=merge, md=md)
    end
    opt = nothing
    if status == 1
        if QUIET == false
            mb = maximum(maximum.(sb))
            println("Obtained the block structure in $time seconds. The maximal size of blocks is $mb.")
        end
        opt,ksupp = ncblockcpop(m, supp, coe, basis, blocks, cl, blocksize, numeq=numeq, QUIET=QUIET, obj=obj, solve=solve)
    end
    data.ksupp = ksupp
    data.sb = sb
    data.numb = numb
    data.blocks = blocks
    data.cl = cl
    data.blocksize = blocksize
    return opt,data
end

function get_ncgsupp(m, supp, gbasis, gblocks, gcl, gblocksize)
    gsupp = Vector{UInt16}[]
    for k = 1:m, i = 1:gcl[k], j = 1:gblocksize[k][i], r = j:gblocksize[k][i], s = 1:length(supp[k+1])
        @inbounds bi = [gbasis[k][gblocks[k][i][j]][end:-1:1]; supp[k+1][s]; gbasis[k][gblocks[k][i][r]]]
        push!(gsupp, bi)
    end
    return gsupp
end

function reducebasis!(psupp, basis, blocks, cl, blocksize)
    init = 0
    flag = 0
    check = 0
    while init == 0 || check > 0
        init = 1
        check = 0
        for i = 1:cl
            if blocksize[i] > 1
                for j = 1:blocksize[i], r = 1:blocksize[i]
                    if j != r
                        @inbounds bi = [basis[blocks[i][j]][end:-1:1]; basis[blocks[i][r]]]
                        if is_sym(bi)
                            push!(psupp, bi)
                        end
                    end
                end
            end
        end
        sort!(psupp)
        unique!(psupp)
        lpsupp = length(psupp)
        for i = 1:cl
            lo = blocksize[i]
            indexb = [k for k=1:lo]
            j = 1
            while lo >= j
                bi = [basis[blocks[i][indexb[j]]][end:-1:1]; basis[blocks[i][indexb[j]]]]
                Locb = ncbfind(psupp, lpsupp, bi)
                if Locb == 0
                   check = 1
                   flag = 1
                   deleteat!(indexb, j)
                   lo = lo-1
                else
                   j += 1
                end
            end
            blocks[i] = blocks[i][indexb]
            blocksize[i] = lo
        end
    end
    if flag == 1
       indexb = blocks[1]
       for i = 2:cl
           indexb = append!(indexb, blocks[i])
       end
       sort!(indexb)
       unique!(indexb)
       return basis[indexb],flag
    else
       return basis,flag
    end
end

function get_nccgraph(ksupp, supp, basis; obj="eigen")
    lb = length(basis)
    lksupp = length(ksupp)
    G = SimpleGraph(lb)
    for i = 1:lb, j = i+1:lb
        r = 1
        while r <= length(supp)
            bi = [basis[i][end:-1:1]; supp[r]; basis[j]]
            bi = reduce!(bi, obj=obj)
            if ncbfind(ksupp, lksupp, bi)!=0
               break
            else
                r += 1
            end
        end
        if r <= length(supp)
           add_edge!(G, i, j)
        end
    end
    return G
end

function get_nccblocks(m, ksupp, gsupp, basis; blocks=[], cl=[], blocksize=[], sb=[], numb=[], TS="block", obj="eigen", QUIET=true, merge=false, md=3)
    if isempty(blocks)
        blocks = Vector{Vector{Vector{UInt16}}}(undef, m+1)
        blocksize = Vector{Vector{UInt16}}(undef, m+1)
        cl = Vector{UInt16}(undef, m+1)
    end
    if TS == false
        for k = 1:m+1
            blocks[k] = [[i for i=1:length(basis[k])]]
            blocksize[k] = [length(basis[k])]
            cl[k] = 1
        end
        nsb = Int.(blocksize[1])
        nnumb = [1]
        status = 1
    else
        G = get_ncgraph(ksupp, basis[1], obj=obj)
        if TS == "block"
            blocks[1] = connected_components(G)
            blocksize[1] = length.(blocks[1])
            cl[1] = length(blocksize[1])
        else
            blocks[1],cl[1],blocksize[1] = chordal_cliques!(G, method=TS, minimize=false)
            if merge == true
                blocks[1],cl[1],blocksize[1] = clique_merge!(blocks[1], d=md, QUIET=true)
            end
        end
        nsb = sort(Int.(unique(blocksize[1])), rev=true)
        nnumb = [sum(blocksize[1].== i) for i in nsb]
        if isempty(sb) || nsb!=sb || nnumb!=numb
            status = 1
            if QUIET == false
                println("------------------------------------------------------")
                println("The sizes of PSD blocks:\n$nsb\n$nnumb")
                println("------------------------------------------------------")
            end
            for k = 1:m
                G = get_nccgraph(ksupp, gsupp[k], basis[k+1], obj=obj)
                if TS == "block"
                    blocks[k+1] = connected_components(G)
                    blocksize[k+1] = length.(blocks[k+1])
                    cl[k+1] = length(blocksize[k+1])
                else
                    blocks[k+1],cl[k+1],blocksize[k+1] = chordal_cliques!(G, method=TS, minimize=false)
                    if merge == true
                        blocks[k+1],cl[k+1],blocksize[k+1] = clique_merge!(blocks[k+1], d=md, QUIET=true)
                    end
                end
            end
        else
            status = 0
            if QUIET == false
                println("No higher NCTSSOS hierarchy!")
            end
        end
    end
    return blocks,cl,blocksize,nsb,nnumb,status
end

function ncblockcpop(m, supp, coe, basis, blocks, cl, blocksize; numeq=0, QUIET=true, obj="eigen", solve=true)
    ksupp = Vector{UInt16}[]
    for i = 1:cl[1], j = 1:blocksize[1][i], r = j:blocksize[1][i]
        @inbounds bi = [basis[1][blocks[1][i][j]][end:-1:1]; basis[1][blocks[1][i][r]]]
        @inbounds push!(ksupp, bi)
    end
    gsupp = get_ncgsupp(m, supp, basis[2:end], blocks[2:end], cl[2:end], blocksize[2:end])
    append!(ksupp, gsupp)
    ksupp = reduce!.(ksupp, obj=obj)
    sort!(ksupp)
    unique!(ksupp)
    lksupp = length(ksupp)
    if QUIET == false
        println("There are $lksupp affine constraints.")
    end
    objv = nothing
    if solve == true
        if QUIET==false
            println("Assembling the SDP...")
        end
        model = Model(optimizer_with_attributes(Mosek.Optimizer))
        set_optimizer_attribute(model, MOI.Silent(), QUIET)
        time = @elapsed begin
        cons = [AffExpr(0) for i=1:lksupp]
        pos = Vector{Union{VariableRef,Symmetric{VariableRef}}}(undef, cl[1])
        for i = 1:cl[1]
            bs = blocksize[1][i]
            if bs == 1
               @inbounds pos[i] = @variable(model, lower_bound=0)
               @inbounds bi = [basis[1][blocks[1][i][1]][end:-1:1]; basis[1][blocks[1][i][1]]]
               bi = reduce!(bi, obj=obj)
               Locb = ncbfind(ksupp, lksupp, bi)
               @inbounds add_to_expression!(cons[Locb], pos[i])
            else
               @inbounds pos[i] = @variable(model, [1:bs, 1:bs], PSD)
               for j = 1:bs, r = j:bs
                   @inbounds bi = [basis[1][blocks[1][i][j]][end:-1:1]; basis[1][blocks[1][i][r]]]
                   bi = reduce!(bi, obj=obj)
                   Locb = ncbfind(ksupp, lksupp, bi)
                   if j == r
                       @inbounds add_to_expression!(cons[Locb], pos[i][j,r])
                   else
                       @inbounds add_to_expression!(cons[Locb], 2, pos[i][j,r])
                   end
               end
            end
        end
        gpos = Vector{Vector{Union{VariableRef,Symmetric{VariableRef}}}}(undef, m)
        for k = 1:m
            gpos[k] = Vector{Union{VariableRef,Symmetric{VariableRef}}}(undef, cl[k+1])
            for i = 1:cl[k+1]
                bs = blocksize[k+1][i]
                if bs == 1
                    if k <= m-numeq
                        gpos[k][i] = @variable(model, lower_bound=0)
                    else
                        gpos[k][i] = @variable(model)
                    end
                    for s = 1:length(supp[k+1])
                        @inbounds bi = [basis[k+1][blocks[k+1][i][1]][end:-1:1]; supp[k+1][s]; basis[k+1][blocks[k+1][i][1]]]
                        bi = reduce!(bi, obj=obj)
                        Locb = ncbfind(ksupp, lksupp, bi)
                        @inbounds add_to_expression!(cons[Locb], coe[k+1][s], gpos[k][i])
                    end
                else
                    if k <= m-numeq
                       gpos[k][i] = @variable(model, [1:bs, 1:bs], PSD)
                    else
                       gpos[k][i] = @variable(model, [1:bs, 1:bs], Symmetric)
                    end
                    for j = 1:bs, r = j:bs, s = 1:length(supp[k+1])
                        @inbounds bi=[basis[k+1][blocks[k+1][i][j]][end:-1:1]; supp[k+1][s]; basis[k+1][blocks[k+1][i][r]]]
                        bi = reduce!(bi, obj=obj)
                        Locb = ncbfind(ksupp, lksupp, bi)
                        if j == r
                            @inbounds add_to_expression!(cons[Locb], coe[k+1][s], gpos[k][i][j,r])
                        else
                            @inbounds add_to_expression!(cons[Locb], 2*coe[k+1][s], gpos[k][i][j,r])
                        end
                    end
                end
            end
        end
        bc = zeros(lksupp)
        for i = 1:length(supp[1])
            Locb = ncbfind(ksupp, lksupp, supp[1][i])
            if Locb == 0
               @error "The monomial basis is not enough!"
               return nothing,nothing
            else
               bc[Locb] = coe[1][i]
           end
        end
        @variable(model, lower)
        cons[1] += lower
        @constraint(model, cons.==bc)
        @objective(model, Max, lower)
        end
        if QUIET == false
            println("SDP assembling time: $time seconds.")
            println("Solving the SDP...")
        end
        time = @elapsed begin
        optimize!(model)
        end
        if QUIET == false
            println("SDP solving time: $time seconds.")
        end
        status = termination_status(model)
        objv = objective_value(model)
        if status != MOI.OPTIMAL
           println("termination status: $status")
           status = primal_status(model)
           println("solution status: $status")
        end
        println("optimum = $objv")
    end
    return objv,ksupp
end


text:
Yao.measure(reg::StabilizerReg; nshot::Integer = 1) = [measure!(copy(reg)) for _ = 1:nshot]
function Yao.measure!(reg::StabilizerReg)
    n = Yao.nqubits(reg)
    if n < 64
        meas_res = zero(Int64)
    elseif n < 128
        meas_res = zero(Int128)
    else
        meas_res = zero(BigInt)
    end
    idn = ['I' for _ = 1:n]
    for l in 1:n
        op_str = copy(idn)
        op_str[l] = 'Z'
        op = QuantumClifford._P_str(string(op_str...))
        _, anti_ind, proj_res = project!(reg.st, op)
        if proj_res === nothing
            proj_res = rand((0x00, 0x02))
            reg.st.phases[anti_ind] = proj_res
        end
        proj_res == 0x02 && (meas_res = flip(meas_res, one(meas_res) << (l-1)))
    end
    return BitStr{n, typeof(meas_res)}(meas_res)
end

function Yao.expect(op::Add, reg::StabilizerReg)
    return sum(Yao.expect(o, reg) for o in op)
end
function Yao.expect(op::Scale, reg::StabilizerReg)
    c = op.alpha
    c === Val(-1) && (c = -1)
    c === Val(im) && (c = im)
    c === Val(-im) && (c = -im)
    return c * Yao.expect(op.content, reg)
end
function Yao.expect(op::PauliString, reg::StabilizerReg)
    Yao.nqubits(op) == Yao.nqubits(reg) || error("The number of qubits does not match!")
    n = Yao.nqubits(op)
    p = zero(PauliOperator, n)
    for i = 1:n
        op[i] === Yao.X && (p[i] = (true, false))
        op[i] === Yao.Y && (p[i] = (true, true))
        op[i] === Yao.Z && (p[i] = (false, true))
    end
    return QuantumClifford.expect(p, reg.st)
end


text:
"""

```julia
map_naive!(f,output,x,box)
```

Internal function or structure - interface may change.

# Extended help

Function that uses the naive pairwise mapping algorithm, for testing.

"""
function map_naive!(f, output, x, box::Box)
    @unpack unit_cell, cutoff_sq = box
    for i in 1:length(x) - 1
        xᵢ = x[i]
        for j in i + 1:length(x)
            xⱼ = wrap_relative_to(x[j], xᵢ, box)
            d2 = norm_sqr(xᵢ - xⱼ)
            if d2 <= cutoff_sq
                output = f(xᵢ, xⱼ, i, j, d2, output)
            end
        end
    end
    return output
end

function map_naive!(f, output, x, y, box::Box)
    @unpack unit_cell, cutoff_sq = box
    for i in 1:length(x)
        xᵢ = x[i]
        for j in 1:length(y)
            yⱼ = wrap_relative_to(y[j], xᵢ, box)
            d2 = norm_sqr(xᵢ - yⱼ)
            if d2 <= cutoff_sq
                output = f(xᵢ, yⱼ, i, j, d2, output)
            end
        end
    end
    return output
end

"""

```
view_celllist_particles(cl::CellList)
```

Internal function or structure - interface may change.

# Extended help

Auxiliary function to view the particles of a computing box, including images created
for computing purposes.

## Example

```julia-repl
julia> box = Box([ 100 50; 50 100 ],10);

julia> x = [ box.unit_cell_max .* rand(SVector{2,Float64}) for i in 1:1000 ];

julia> cl = CellList(x,box);

julia> p = CellListMap.view_celllist_particles(cl);

julia> using Plots

julia> scatter(Tuple.(p),label=nothing,xlims=(-10,180),ylims=(-10,180))

```

"""
function view_celllist_particles(cl::CellList{N,T}) where {N,T}
    x = SVector{N,T}[]
    for icell in 1:cl.n_cells_with_particles
        cell = cl.cells[icell]
        for ip in 1:cell.n_particles
            p = cell.particles[ip]
            push!(x, p.coordinates)
        end
    end
    return x
end

test_map2(box,cl;parallel=true) =
    map_pairwise!(
        (x, y, i, j, d2, s) -> begin
            s += 1
        end,
        0, box, cl, parallel=parallel
    )

function missing_pair(x1,x2)
    for x in eachrow(x2)
        found = false
        for y in eachrow(x1)
            (x[1] == y[1]) && (x[2] == y[2]) && (found = true)
            (x[1] == y[2]) && (x[2] == y[1]) && (found = true)
            if found
                break
            end
        end
        if !found
            return x
        end
    end
end

test_map(box,cl;parallel=true) = map_pairwise!((x, y, i, j, d2, s) -> s += 1/d2, 0., box, cl, parallel=parallel)
test_naive(x,box) = CellListMap.map_naive!((x, y, i, j, d2, s) -> s += 1/d2, 0., x, box)
function simple_test(x,box;parallel=true)
    cl = CellList(x,box,parallel=parallel)
    r_naive = test_naive(x,box)
    r_map = test_map(box,cl,parallel=parallel)
    println("naive = $r_naive")
    println("map   = $r_map")
    return "Test passed: $(r_naive ≈ r_map)"
end

function check_random_cells(
    N=3, M=10;
    show_progress=true,
    lcell=1,
    parallel=false,
    UnitCellType=OrthorhombicCell
)
    local x, box
    ntrial = 0
    show_progress && (p = Progress(10000, 0.5))
    while ntrial < 2000
        unit_cell_matrix = zeros(SMatrix{N,N,Float64})
        if UnitCellType == OrthorhombicCell
            for i in 1:N
                @set! unit_cell_matrix[i,i] = 1 + 10*rand()
            end
        else
            for i in 1:N, j in 1:N
                @set! unit_cell_matrix[i,j] = 1 + 10*rand()
            end
        end
        cutoff = 1 + rand()
        try
            if !check_unit_cell(unit_cell_matrix, cutoff, printerr=false)
                continue
            end
        catch
            @show unit_cell_matrix, cutoff
            return false
        end
        if UnitCellType == OrthorhombicCell
            box = Box([unit_cell_matrix[i,i] for i in 1:N], cutoff, lcell=lcell)
        else
            box = Box(unit_cell_matrix, cutoff, lcell=lcell)
        end
        if prod(box.nc) > 100000
            continue
        end
        show_progress && (next!(p))
        x = 10 * rand(SVector{N,Float64}, M)
        for i in eachindex(x)
            x[i] = x[i] .- 50
        end
        cl = CellList(x, box, parallel=parallel)
        test = test_map(box, cl, parallel=parallel)
        if test ≈ 0
            continue
        end
        ntrial += 1
        if !(test ≈ test_naive(x, box))
            show_progress && println("FOUND PROBLEMATIC SETUP.")
            return false, x, box
        else

        end
        if parallel
            aux = AuxThreaded(cl)
            cl = UpdateCellList!(x, box, cl, aux, parallel=parallel)
        else
            cl = UpdateCellList!(x, box, cl, parallel=parallel)
        end
        if !(test_map(box, cl, parallel=parallel) ≈ test_naive(x, box))
            show_progress && println("FOUND PROBLEMATIC SETUP.")
            return false, x, box
        end
    end
    show_progress && println(" ALL PASSED! ")
    return true, x, box
end

function drawbox(box::Box{UnitCellType,2}) where UnitCellType
    S = SVector{2,Float64}
    m = box.unit_cell.matrix
    x = [
        S(0., 0.),
        S(m[:,1]),
        S(m[:,1] + m[:,2]),
        S(m[:,2]),
        S(0., 0.)
    ]
    return x
end

function drawbox(box::Box{UnitCellType,3}) where UnitCellType
    S = SVector{3,Float64}
    m = box.unit_cell.matrix
    x = [
        S(0., 0., 0.),
        S(m[:,1]),
        S(m[:,1] + m[:,2]),
        S(m[:,2]),
        S(0., 0., 0.),
        S(m[:,1]),
        S(m[:,1] + m[:,3]),
        S(m[:,3]),
        S(0., 0., 0.),
        S(m[:,2]),
        S(m[:,2] + m[:,3]),
        S(m[:,2]),
        S(0., 0., 0.),
        S(m[:,1]),
        S(m[:,1] + m[:,2]),
        S(m[:,1] + m[:,2]) + m[:,3],
        S(m[:,1] + m[:,3]),
        S(m[:,3]),
        S(m[:,3]) + m[:,2],
        S(m[:,1] + m[:,2]) + m[:,3]
    ]
    return x
end

"""

```
draw_computing_cell(x,box::Box{UnitCellType,2}) where UnitCellType
```

```
draw_computing_cell(cl::CellList,box::Box{UnitCellType,2},x) where UnitCellType
```

This function creates a plot of the computing cell, in two dimensions.

"""
function draw_computing_cell(x, box::Box{UnitCellType,2};
    parallel=true,
    xticks=nothing,
    yticks=nothing,
) where UnitCellType
    cl = CellList(x, box, parallel=parallel)
    return draw_computing_cell(cl,box,x,xticks=xticks,yticks=yticks)
end
function draw_computing_cell(
    cl::CellList, box::Box{UnitCellType,2}, x;
    xticks = nothing,
    yticks = nothing,
) where UnitCellType
    box_points = drawbox(box)
    p = view_celllist_particles(cl)
    plt = Main.plot()
    Main.plot!(plt, Tuple.(box_points), label=:none)
    Main.scatter!(plt, Tuple.(p), label=:none, markeralpha=0.3)
    Main.scatter!(plt, Tuple.(wrap_to_first.(x, Ref(box))), label=:none)
    xmin = minimum(el[1] for el in p) - 3 * box.cell_size[1]
    xmax = maximum(el[1] for el in p) + 3 * box.cell_size[1]
    ymin = minimum(el[2] for el in p) - 3 * box.cell_size[2]
    ymax = maximum(el[2] for el in p) + 3 * box.cell_size[2]
    isnothing(xticks) && (xticks=(round.(digits=3, xmin:box.cell_size[1]:xmax)))
    isnothing(yticks) && (yticks=(round.(digits=3, ymin:box.cell_size[2]:ymax)))
    Main.plot!(plt,
        aspect_ratio=1,framestyle=:box,xrotation=60,
        xlims=(xmin, xmax),
        ylims=(ymin, ymax),
        xticks=xticks,
        yticks=yticks,
    )
    return plt
end

"""

```
draw_computing_cell(x,box::Box{UnitCellType,3}) where UnitCellType
```

This function creates a plot of the computing cell, in three dimensions.

"""
function draw_computing_cell(x, box::Box{UnitCellType,3}; parallel=true) where UnitCellType
    cl = CellList(x, box, parallel=parallel)
    box_points = drawbox(box)
    p = view_celllist_particles(cl)
    plt = Main.plot()
    Main.plot!(plt, Tuple.(box_points), label=:none)
    Main.scatter!(plt, Tuple.(p), label=:none, markeralpha=0.3)
    Main.scatter!(plt, Tuple.(wrap_to_first.(x, Ref(box))), label=:none, markeralpha=0.3)
    lims = Vector{Float64}[]
    for i in 1:3
        push!(lims, [ -2 * box.cell_size[i], box.nc[i] + 2 * box.cell_size[i] ])
    end
    Main.plot!(plt,
        aspect_ratio=1,framestyle=:box,xrotation=60,yrotation=-70,zrotation=0,
        xlims=lims[1],
        ylims=lims[2],
        zlims=lims[3],
        xticks=(round.(digits=3, lims[1][1]:box.cell_size[1]:lims[1][2])),
        yticks=(round.(digits=3, lims[2][1]:box.cell_size[2]:lims[2][2])),
        zticks=(round.(digits=3, lims[3][1]:box.cell_size[3]:lims[3][2])),
    )
    return plt
end

function compare_cells(cl1::CellList{N,T}, cl2::CellList) where {N,T}

    if cl1.n_cells_with_real_particles != cl2.n_cells_with_real_particles
        println("Error: n_cells_with_real_particles differ")
    end
    if cl1.n_particles != cl2.n_particles
        println("Error: n_particles differ")
    end
    for i in 1:cl1.n_cells_with_particles
        linear_index =
        j = findfirst(
            isequal(cl1.cell_indices[i]),
            @view(cl2.cell_indices[1:cl2.n_cells_with_particles])
        )
        if isnothing(j)
            println(" Warning: cells_with_particles: Could not find index $(cl1.cell_indices[i])")
            break
        end
    end
    for i in 1:cl1.n_cells_with_real_particles
        j = findfirst(
            isequal(cl1.cell_indices_real[i]),
            @view(cl2.cell_indices_real[1:cl2.n_cells_with_real_particles])
        )
        if isnothing(j)
            println(" Warning: cells_with_real_particles: Could not find index $(cl1.cell_indices_real[i])")
            break
        end
    end
    if length(cl1.projected_particles) != length(cl2.projected_particles)
        println("Error: length of project_particles differ.")
    end
    ncheck = 0
    nparticle_check = 0
    if length(cl1.cells) != length(cl2.cells)
        println("Error: lengths of cells lists differ")
    else
        differ = false
        for i in 1:cl1.n_cells_with_particles
            if differ
                break
            end
            ncheck += 1
            ci = cl1.cells[i]
            j = findfirst(c -> c.linear_index == ci.linear_index, cl2.cells)
            if isnothing(j)
                println("Error: linear index of $(ci.linear_index) not found.")
                differ = true
            else
                cj = cl2.cells[j]
                if ci.cartesian_index != cj.cartesian_index
                    println("Error: cartesian_index differ")
                    differ = true
                end
                if !(ci.center ≈ cj.center)
                    println("Error: centers differ")
                    differ = true
                end
                if ci.contains_real != cj.contains_real
                    println("Error: contains_real differ")
                    differ = true
                end
                if ci.n_particles != cj.n_particles
                    println("Error: n_particles differ")
                    differ = true
                end
                for p1 in ci.particles
                    j = findfirst(p -> p.index == p1.index, cj.particles)
                    if isnothing(j)
                        println("Error: could not find particle in cl2.")
                        @show p1
                        differ = true
                    else
                        nparticle_check += 1
                        p2 = cj.particles[j]
                        if p1 != p2
                            println("Error: particles differ. ")
                            @show p1
                            @show p2
                            differ = true
                            break
                        end
                    end
                    differ && break
                end
                if sum(p -> p.index, ci.particles[1:ci.n_particles]) !=
                   sum(p -> p.index, cj.particles[1:cj.n_particles])
                    println("Error: particles of cells differ")
                    differ = true
                end
                if differ
                    println("Error: on cell with linear_index: $(ci.linear_index)")
                end
            end
        end
    end
    println(" Checked $ncheck of $(cl1.n_cells_with_particles) cells.")
    println(" Checked $nparticle_check of $(cl1.n_particles) particles.")
end

function check_cl(cl,box;mark)
    for i in 1:cl.n_cells_with_particles-1
        cellᵢ = cl.cells[i]
        for j in i+1:cl.n_cells_with_particles
            cellⱼ = cl.cells[j]
            for ip in 1:cellᵢ.n_particles
                xᵢ = cellᵢ.particles[ip]
                for jp in 1:cellⱼ.n_particles
                    xⱼ = cellⱼ.particles[jp]
                    if norm(xᵢ.coordinates - xⱼ.coordinates) ≈ 0
                        println("duplicate found:")
                        @show mark
                        @show Threads.threadid()
                        @show i, j, cellᵢ.linear_index, cellⱼ.linear_index
                        @show cellᵢ.n_particles
                        @show cellⱼ.n_particles
                        @show xᵢ
                        @show xⱼ
                        @show particle_cell(xᵢ.coordinates,box)
                        @show particle_cell(xⱼ.coordinates,box)
                        @show cell_linear_index(box.nc,particle_cell(xᵢ.coordinates,box))
                        @show cell_linear_index(box.nc,particle_cell(xⱼ.coordinates,box))
                        error()
                    end
                end
            end
        end
    end
end

# Create a set of random particles with an atomic density simiilar to that
# of water at room temperature and pressure.
function xatomic(n=100_000;cutoff=12.,lcell=1,type=:vec)
    atomic_density_of_water = 0.1 # atoms/ Angstrom^3
    vol = n/atomic_density_of_water
    l = vol^(1/3)
    if type == :mat
        x = l*rand(3,n)
    elseif type == :vec
        x = [ l*rand(SVector{3,Float64}) for _ in 1:n ]
    end
    box = Box([l,l,l],cutoff,lcell=lcell)
    return x, box
end

# Create a set of  random particles with the typical density of galaxies
function xgalactic(n=8_000_000;cutoff=5.,lcell=1,type=:vec)
    galactic_density = 12 # galaxyes per cubic mega-parsecs (https://doi.org/10.1093/mnras/106.2.121)
    vol = n/galactic_density
    l = vol^(1/3) # mega-parsecs
    if type == :mat
        x = l*rand(3,n)
    elseif type == :vec
        x = [ l*rand(SVector{3,Float64}) for _ in 1:n ]
    end
    box = Box([l,l,l],cutoff,lcell=lcell)
    return x, box
end