text:
const NOpsType = Int#Union{Int,Vector{Int}}

struct UpperBoundedInteger{N, T<: Base.BitInteger} <: Integer
  i::T
  @inline UpperBoundedInteger{N}(i::T) where {N,T<:Base.BitInteger} = new{N,T}(i)
end

@inline UpperBoundedInteger(i::T, ::StaticInt{N}) where {N,T<:Base.BitInteger} = UpperBoundedInteger{N}(i)
@inline UpperBoundedInteger(::StaticInt{M}, ::StaticInt{N}) where {N,M} = StaticInt{M}()
@inline UpperBoundedInteger{N}(::StaticInt{M}) where {N,M} = StaticInt{M}()
@inline Base.:(%)(a::UpperBoundedInteger, ::Type{T}) where {T<:Base.BitInteger} = a.i % T
Base.promote_rule(::Type{T}, ::Type{UpperBoundedInteger{N,S}}) where {N,T<:Base.BitInteger,S} = promote_rule(T,S)
Base.promote_rule(::Type{UpperBoundedInteger{N,S}}, ::Type{T}) where {N,T<:Base.BitInteger,S} = promote_rule(S,T)
Base.convert(::Type{T}, i::UpperBoundedInteger) where {T<:Number} = convert(T, i.i)
Base.convert(::Type{UpperBoundedInteger{N,T}}, i::UpperBoundedInteger{N,T}) where {N,T<:Base.BitInteger} = i
Base.convert(::Type{Any}, i::UpperBoundedInteger) = i
upper_bound(_) = typemax(Int)
upper_bound(::Type{CloseOpen{T,UpperBoundedInteger{N,S}}}) where {T,N,S} = N - 1

@inline Base.last(r::CloseOpen{<:Integer,<:UpperBoundedInteger}) = getfield(getfield(r,:upper),:i) - One()
@inline ArrayInterface.static_last(r::CloseOpen{<:Integer,<:UpperBoundedInteger}) = getfield(getfield(r,:upper),:i) - One()
@inline Base.length(r::CloseOpen{<:Integer,<:UpperBoundedInteger}) = getfield(getfield(r,:upper),:i) - getfield(r,:start)
@inline Base.length(r::CloseOpen{Zero,<:UpperBoundedInteger}) = getfield(getfield(r,:upper),:i)

function Loop(ls::LoopSet, ex::Expr, sym::Symbol, f, s, l, ub::Int)
    if (f !== nothing) && (s !== nothing) && (l !== nothing)
        return static_loop(sym, f, s, l)
    end
    ssym = String(sym)
    rangesym = gensym(ssym * "_loop");
    lensym = gensym(ssym * "_looplen")
    pushpreamble!(ls, Expr(:(=), rangesym, ex))
    pushpreamble!(ls, Expr(:(=), lensym, Expr(:call, GlobalRef(ArrayInterface,:static_length), rangesym)))
    F = if f === nothing
        start = gensym(ssym*"_loopstart")
        pushpreamble!(ls, Expr(:(=), start, Expr(:call, %, Expr(:call, lv(:first), rangesym), Int)))
        MaybeKnown(start, 1)
    else
        MaybeKnown(f)
    end
    S = if s === nothing
        step = gensym(ssym*"_loopstep")
        pushpreamble!(ls, Expr(:(=), step, Expr(:call, %, Expr(:call, lv(:step), rangesym), Int)))
        MaybeKnown(step, 1)
    else
        MaybeKnown(s)
    end
    L = if l === nothing
        stop = gensym(ssym*"_loopstop")
        pushpreamble!(ls, Expr(:(=), stop, Expr(:call, %, Expr(:call, lv(:last), rangesym), Int)))
        MaybeKnown(stop, min(ub, 1024))
    else
        MaybeKnown(l)
    end
    loopiteratesatleastonce!(ls, Loop(sym, F, L, S, rangesym, lensym))
end
function Loop(ls::LoopSet, ex::Expr, sym::Symbol, ::Type{R}) where {R<:AbstractRange}
  f = ArrayInterface.known_first(R)
  s = ArrayInterface.known_step(R)
  l = ArrayInterface.known_last(R)
  ub = upper_bound(R)
  Loop(ls, ex, sym, f, s, l, ub)
end

function static_loop(sym::Symbol, L::Int, S::Int, U::Int)
    Loop(sym, MaybeKnown(L,0), MaybeKnown(U,0), MaybeKnown(S,0), Symbol(""), Symbol(""))
end
function Loop(::LoopSet, ::Expr, sym::Symbol, ::Type{OptionallyStaticUnitRange{Static{L}, Static{U}}}) where {L,U}
    static_loop(sym, L, 1, U)
end
function Loop(::LoopSet, ::Expr, sym::Symbol, ::Type{ArrayInterface.OptionallyStaticStepRange{StaticInt{L}, StaticInt{S}, StaticInt{U}}}) where {L,S,U}
    static_loop(sym, L, S, U)
end
function Loop(::LoopSet, ::Expr, sym::Symbol, ::Type{CloseOpen{Static{L}, Static{U}}}) where {L,U}
    static_loop(sym, L, 1, U - 1)
end

extract_loop(l) = Expr(:call, GlobalRef(Core,:getfield), Symbol("#loop#bounds#"), l, false)

function add_loops!(ls::LoopSet, LPSYM, LB)
    n = max(length(LPSYM), length(LB))
    for i = 1:n
        sym, l = LPSYM[i], LB[i]
        if l<:CartesianIndices
            add_loops!(ls, i, sym, l)
        else
            add_loop!(ls, Loop(ls, extract_loop(i), sym, l)::Loop)
            push!(ls.loopsymbol_offsets, ls.loopsymbol_offsets[end]+1)
        end
    end
end
function add_loops!(ls::LoopSet, i::Int, sym::Symbol, @nospecialize(l::Type{<:CartesianIndices}))
    N, T = l.parameters
    ssym = String(sym)
    for k = N:-1:1
        axisexpr = :(getfield(getfield(getfield(var"#loop#bounds#", $i, false), :indices), $k, false))
        add_loop!(ls, Loop(ls, axisexpr, Symbol(ssym*'#'*string(k)*'#'), T.parameters[k])::Loop)
    end
    push!(ls.loopsymbol_offsets, ls.loopsymbol_offsets[end]+N)
end

function ArrayReferenceMeta(
    ls::LoopSet, @nospecialize(ar::ArrayRefStruct), arraysymbolinds::Vector{Symbol},
    opsymbols::Vector{Symbol}, nopsv::Vector{NOpsType}, expandedv::Vector{Bool}
)
    # unpack the `ArrayRefStruct`
    # we don't want to specialize on it, as it is typed on symbols.
    index_types = (ar.index_types)::UInt128
    indices = (ar.indices)::UInt128
    offsets = (ar.offsets)::UInt128
    strides = (ar.strides)::UInt128
    arrayar, ptrar = array_and_ptr(ar)::Tuple{Symbol,Symbol}
    # ptrar = ptr(ar)::Symbol
    # TODO, check if this matters at all. The compiler still knows it is an `::ArrayRefStruct`, just not `arrayar` or `ptrar`?
    ArrayReferenceMeta(
        ls, index_types, indices, offsets, strides,
        arrayar, ptrar, arraysymbolinds, opsymbols, nopsv, expandedv
    )
end
function ArrayReferenceMeta(
    ls::LoopSet, index_types::UInt128, indices::UInt128, offsets::UInt128, strides::UInt128,
    arrayar::Symbol, ptrar::Symbol, arraysymbolinds::Vector{Symbol},
    opsymbols::Vector{Symbol}, nopsv::Vector{NOpsType}, expandedv::Vector{Bool}
)
    ni = filled_8byte_chunks(index_types)
    index_vec = Symbol[]
    offset_vec = Int8[]
    stride_vec = Int8[]
    loopedindex = Bool[]
    while index_types != zero(UInt128)
        ind = indices % UInt8
        offsetᵢ = offsets % Int8
        strideᵢ = strides % Int8
        if index_types == LoopIndex
            if ind == zero(Int8) # CONSTANTZEROINDEX
                pushfirst!(index_vec, CONSTANTZEROINDEX)
                pushfirst!(offset_vec, offsetᵢ)
                pushfirst!(stride_vec, strideᵢ)
                pushfirst!(loopedindex, true)
            else
                for inda in ls.loopsymbol_offsets[ind]+1:ls.loopsymbol_offsets[ind+1]
                    pushfirst!(index_vec, ls.loopsymbols[inda])
                    pushfirst!(offset_vec, offsetᵢ)
                    pushfirst!(stride_vec, strideᵢ)
                    pushfirst!(loopedindex, true)
                end
            end
        else#if index_types == ComputedIndex
            @assert index_types == ComputedIndex
            opsym = opsymbols[ind]
            if expandedv[ind]
                nops = nopsv[ind]
                for j ∈ 0:nops-1
                    pushfirst!(index_vec, expandedopname(opsym, j))
                    pushfirst!(offset_vec, offsetᵢ)
                    pushfirst!(stride_vec, strideᵢ)
                    pushfirst!(loopedindex, false)
                end
            else
                pushfirst!(index_vec, opsym)
                pushfirst!(offset_vec, offsetᵢ)
                pushfirst!(stride_vec, strideᵢ)
                pushfirst!(loopedindex, false)
            end
        # else
            # @assert index_types == SymbolicIndex
            # pushfirst!(index_vec, arraysymbolinds[ind])
            # pushfirst!(offset_vec, offset)
            # pushfirst!(loopedindex, false)
        end
        index_types >>>= 8
        indices >>>= 8
        offsets >>>= 8
        strides >>>= 8
        ni -= 1
    end
    ArrayReferenceMeta(
        ArrayReference(arrayar, index_vec, offset_vec, stride_vec),
        loopedindex, ptrar
    )
end

extract_varg(i) = :(getfield(var"#vargs#", $i, false))
# _extract(::Type{Static{N}}) where {N} = N
extract_gsp!(sptrs::Expr, name::Symbol) = (push!(sptrs.args, name); nothing)
tupleranks(R::NTuple{8,Int}) = ntuple(n -> sum(R[n] .≥ R), Val{8}())
function rank_to_sortperm((R,N)::Tuple{NTuple{8,Int},Int})::Vector{Int}
  r = tupleranks(R)
  sp = Vector{Int}(undef, N)
  @inbounds for n = 1:N
    sp[r[n]] = n
  end
  sp
end
function loop_indexes_bit!(ls::LoopSet, ar::ArrayReferenceMeta)
    li = ar.loopedindex;
    ind = first(getindices(ar))
    ind === DISCONTIGUOUS && return
    first(li) || throw(LoopError("The contiguous index of a `BitArray` shouldn't be a complex function.", )ind)
    ls.loopindexesbit[getloopid(ls,ind)] = true
    nothing
end
function add_mref!(
    sptrs::Expr, ls::LoopSet, ar::ArrayReferenceMeta, @nospecialize(_::Type{Ptr{T}}),
    C::Int, B::Int, sp::Vector{Int}, name::Symbol
) where {T}
  @assert B ≤ 0 "Batched arrays not supported yet."
  _add_mref!(sptrs, ls, ar, typetosym(T), C, B, sp, name)
  sizeof(T)
end
typetosym(::Type{T}) where {T<:NativeTypes} = (VectorizationBase.JULIA_TYPES[T])::Symbol
typetosym(T) = T
function _add_mref!(sptrs::Expr, ls::LoopSet, ar::ArrayReferenceMeta, T_sym, C::Int, B::Int, sp::Vector{Int}, name::Symbol)
    # maybe no change needed? -- optimize common case
  li = ar.loopedindex;
  if sp == eachindex(sp) || isone(length(li))
    # don't set `bit` to true if our vector width is ≥ 8
    ((T_sym === :Bit) && (ls.vector_width < 8)) && loop_indexes_bit!(ls, ar)
    return extract_gsp!(sptrs, name)
  end
  permute_mref!(ar, C, sp)
  # must now sort array's inds, and stack pointer's
  tmpsp = gensym(name)
  extract_gsp!(sptrs, tmpsp)
  strd_tup = Expr(:tuple)
  offsets_tup = Expr(:tuple)
  gf = GlobalRef(Core,:getfield)
  offsets = gensym(:offsets); strides = gensym(:strides)
  pushpreamble!(ls, Expr(:(=), offsets, Expr(:call, gf, tmpsp, QuoteNode(:offsets))))
  pushpreamble!(ls, Expr(:(=), strides, Expr(:call, gf, tmpsp, QuoteNode(:strd))))
  for (i, p) ∈ enumerate(sp)
    push!(strd_tup.args, Expr(:call, gf, strides, p, false))
    push!(offsets_tup.args, Expr(:call, gf, offsets, p, false))
  end
  #TODO: fix for `T === Bit`.
  column_major = Expr(:tuple)
  for n ∈ eachindex(sp)
    push!(column_major.args, n)
  end
  sptype = Expr(:curly, lv(:StridedPointer), T_sym, length(sp), (C == -1 ? -1 : 1), B, column_major)
  sptr = Expr(:call, sptype, Expr(:call, :pointer, tmpsp), strd_tup, offsets_tup)
  pushpreamble!(ls, Expr(:(=), name, sptr))
  nothing
end
function permute_mref!(ar::ArrayReferenceMeta, C::Int, sp::Vector{Int})
  sp == eachindex(sp) && return nothing
  li = ar.loopedindex; lic = copy(li);
  inds = getindices(ar); indsc = copy(inds);
  offsets = ar.ref.offsets; offsetsc = copy(offsets);
  strides = ar.ref.strides; stridesc = copy(strides);
  for (i, p) ∈ enumerate(sp)
    li[i] = lic[p]
    inds[i] = indsc[p]
    offsets[i] = offsetsc[p]
    strides[i] = stridesc[p]
  end
  C == -1 && makediscontiguous!(getindices(ar))
  return nothing
end
function add_mref!(
    sptrs::Expr, ::LoopSet, ::ArrayReferenceMeta, @nospecialize(_::Type{VectorizationBase.FastRange{T,F,S,O}}),
    ::Int, ::Int, sp::Vector{Int}, name::Symbol
) where {T,F,S,O}
  extract_gsp!(sptrs, name)
  sizeof(T)
end
function create_mrefs!(
    ls::LoopSet, arf::Vector{ArrayRefStruct}, as::Vector{Symbol}, os::Vector{Symbol},
    nopsv::Vector{NOpsType}, expanded::Vector{Bool}, ::Type{Tuple{}}
)
    length(arf) == 0 || throw(ArgumentError("Length of array ref vector should be 0 if there are no stridedpointers."))
    Vector{ArrayReferenceMeta}(undef, length(arf)), Int[]
end
function stabilize_grouped_stridedpointer_type(C, B, R)
  N = (length(C))::Int
  Cv = Vector{Int}(undef, N)
  Bv = Vector{Int}(undef, N)
  Rv = Vector{Tuple{NTuple{8,Int},Int}}(undef, N)
  for n ∈ 1:N
    Cv[n] = C[n]
    Bv[n] = B[n]
    Rₙ = R[n]
    let L::Int = length(Rₙ)
      Rv[n] = (ntuple(i -> i > L ? typemax(Int) : (Rₙ[i])::Int, Val(8))::NTuple{8,Int}, L)
    end
  end
  Cv, Bv, Rv
end
function create_mrefs!(
  ls::LoopSet, arf::Vector{ArrayRefStruct}, as::Vector{Symbol}, os::Vector{Symbol},
  nopsv::Vector{NOpsType}, expanded::Vector{Bool}, @nospecialize(_::Type{VectorizationBase.GroupedStridedPointers{P,C,B,R,I,X,O}})
) where {P,C,B,R,I,X,O}

  Cv,Bv,Rv = stabilize_grouped_stridedpointer_type(C, B, R)
  _create_mrefs!(ls, arf, as, os, nopsv, expanded, P.parameters, Cv, Bv, Rv)
end
function _create_mrefs!(
  ls::LoopSet, arf::Vector{ArrayRefStruct}, as::Vector{Symbol}, os::Vector{Symbol},
  nopsv::Vector{NOpsType}, expanded::Vector{Bool}, P::Core.SimpleVector, C::Vector{Int}, B::Vector{Int}, R::Vector{Tuple{NTuple{8,Int},Int}}
)
    mrefs::Vector{ArrayReferenceMeta} = Vector{ArrayReferenceMeta}(undef, length(arf))
    elementbytes::Vector{Int} = Vector{Int}(undef, length(arf))
    sptrs = Expr(:tuple)
    # pushpreamble!(ls, Expr(:(=), sptrs, :(VectorizationBase.stridedpointers(getfield(vargs, 1, false)))))
    pushpreamble!(ls, Expr(:(=), sptrs, :(VectorizationBase.stridedpointers(getfield(var"#vargs#", 1, false)))))
    j = 0
    rank_to_sps = Vector{Tuple{Int,Vector{Int}}}(undef, length(arf))
    for i ∈ eachindex(arf)
        ar = ArrayReferenceMeta(ls, arf[i], as, os, nopsv, expanded)
        duplicate = false
        vptrar = vptr(ar)
        for k ∈ 1:i-1
            if vptr(mrefs[k]) === vptrar
                duplicate = true
                # if isassigned(rank_to_sps, k)
                Cₖ, sp = rank_to_sps[k]
                permute_mref!(ar, Cₖ, sp)
                elementbytes[i] = elementbytes[k]
                # end
                break
            end
        end
        if !duplicate
            j += 1
            sp = rank_to_sortperm(R[j])::Vector{Int}
            rank_to_sps[i] = (C[j],sp)
            elementbytes[i] = add_mref!(sptrs, ls, ar, P[j], C[j], B[j], sp, vptr(ar))
        end
        mrefs[i] = ar
    end
    mrefs, elementbytes
end

function num_parameters(AM)
    num_param::Int = AM[1]
    # num_param += length(AM[2].parameters)
    num_param + length(AM[3])
end
function gen_array_syminds(AM)
    Symbol[Symbol("##arraysymbolind##"*i*'#') for i ∈ 1:(AM[1])::Int]
end
function process_metadata!(ls::LoopSet, AM, extractind::Int)
    opoffsets = ls.operation_offsets
    expandbyoffset!(ls.outer_reductions, AM[2], opoffsets)
    for (i,si) ∈ enumerate(AM[3])
        sii = si::Int
        opid = opoffsets[sii]+1
        if instruction(operations(ls)[opid]) ≠ DROPPEDCONSTANT
            s = gensym(:symlicm)
            push!(ls.preamble_symsym, (opid, s))
            pushpreamble!(ls, Expr(:(=), s, extract_varg((extractind += 1))))
        end
    end
    expandbyoffset!(ls.preamble_symint, AM[4], opoffsets)
    expandbyoffset!(ls.preamble_symfloat, AM[5], opoffsets)
    expandbyoffset!(ls.preamble_zeros, AM[6], opoffsets)
    expandbyoffset!(ls.preamble_funcofeltypes, AM[7], opoffsets)
    return extractind
end
function expandbyoffset!(indexpand::Vector{T}, inds, offsets::Vector{Int}, expand::Bool = true) where {T <: Union{Int,Tuple{Int,<:Any}}}
    for _ind ∈ inds
        ind = T === Int ? _ind : first(_ind)
        base = offsets[ind] + 1
        for inda ∈ base:(expand ? offsets[ind+1] : base)
            if T === Int
                push!(indexpand, inda)
            elseif T === Tuple{Int,Tuple{Int,Int32,Bool}}
                li = last(_ind)
                push!(indexpand, (inda,(li[1],li[2],li[3])))
            else
                push!(indexpand, (inda,last(_ind)))
            end
        end
    end
    indexpand
end
expandbyoffset(inds::Vector{Int}, offsets::Vector{Int}, expand::Bool) = expandbyoffset!(Int[], inds, offsets, expand)
function loopindex!(idxs::Vector{Int}, ls::LoopSet, u::Unsigned, shift::Unsigned)
  mask = (one(shift) << shift) - one(shift) # mask to zero out all but shift-bits
  while u != zero(u)
    pushfirst!(idxs, ( u % typeof(shift) ) & mask)
    u >>= shift
  end
  idxs
end
loopindex(ls::LoopSet, u::Unsigned, shift::Unsigned) = reverse!(loopindex!(Int[], ls, u, shift))
function loopindexoffset(ls::LoopSet, u::Unsigned, li::Bool, expand::Bool = false)
    if li
        shift = 0x04
        offsets = ls.loopsymbol_offsets
    else
        shift = 0x08
        offsets = ls.operation_offsets
    end
    idxs = loopindex(ls, u, shift)
    expandbyoffset(idxs, offsets, expand)
end
function parents_symvec(ls::LoopSet, u::Unsigned, expand, offset)
    idxs = loopindexoffset(ls, u, true, expand)   # TODO DRY  (undesirable that this gets hard-coded in multiple places)
    return Symbol[getloopsym(ls, i + offset) for i ∈ idxs]
end
loopdependencies(ls::LoopSet, os::OperationStruct, expand = false, offset = 0) = parents_symvec(ls, os.loopdeps, expand, offset)
reduceddependencies(ls::LoopSet, os::OperationStruct, expand = false, offset = 0) = parents_symvec(ls, os.reduceddeps, expand, offset)
childdependencies(ls::LoopSet, os::OperationStruct, expand = false, offset = 0) = parents_symvec(ls, os.childdeps, expand, offset)

# parents(ls::LoopSet, u::UInt128) = loopindexoffset(ls, u, false)
function parents(ls::LoopSet, u₀::UInt128, u₁::UInt128)
  idxs = Int[]
  u₁ == zero(u₁) || loopindex!(idxs, ls, u₁, 0x0010)
  loopindex!(idxs, ls, u₀, 0x0010)
  reverse!(idxs)
end
parents(ls::LoopSet, os::OperationStruct) = parents(ls, os.parents₀, os.parents₁)

expandedopname(opsymbol::Symbol, offset::Integer) = Symbol(String(opsymbol)*'#'*string(offset+1)*'#')
function calcnops(ls::LoopSet, os::OperationStruct)
    optyp = optype(os)
    if (optyp != loopvalue) && (optyp != compute)
        return 1
    end
    offsets = ls.loopsymbol_offsets
    idxs = loopindex(ls, os.loopdeps, 0x04)  # TODO DRY
    iszero(length(idxs)) && return 1
    return maximum(i->offsets[i+1]-offsets[i], idxs)
end
function isexpanded(ls::LoopSet, ops::Vector{OperationStruct}, nopsv::Vector{NOpsType}, i::Int)
    nops = nopsv[i]
    # nops isa Vector{Int} only if accesses_memory(os), which means isexpanded must be false
    (nops === 1 || isa(nops, Vector{Int})) && return false
    os = ops[i]
    optyp = optype(os)
    if optyp == compute
        any(j -> isexpanded(ls, ops, nopsv, j), parents(ls, os))
    elseif optyp == loopvalue
        true
    else
        false
    end
end
function mref_elbytes(os::OperationStruct, mrefs::Vector{ArrayReferenceMeta}, elementbytes::Vector{Int})
  if isload(os) | isstore(os)
    mrefs[os.array], elementbytes[os.array]
  else
    NOTAREFERENCE, 4
  end
end
function add_op!(
    ls::LoopSet, instr::Instruction, ops::Vector{OperationStruct}, nopsv::Vector{NOpsType}, expandedv::Vector{Bool}, i::Int,
    mrefs::Vector{ArrayReferenceMeta}, opsymbol, elementbytes::Vector{Int}
)
    os = ops[i]
    mref, elbytes = mref_elbytes(os, mrefs, elementbytes)
    # opsymbol = (isconstant(os) && instr != LOOPCONSTANT) ? instr.instr : opsymbol
    # If it's a CartesianIndex add or subtract, we may have to add multiple operations
    expanded = expandedv[i]# isexpanded(ls, ops, nopsv, i)
    opoffsets = ls.operation_offsets
    # offsets = ls.loopsymbol_offsets
    optyp = optype(os)
    if !expanded
        op = Operation(
            length(operations(ls)), opsymbol, elbytes, instr,
            optyp, loopdependencies(ls, os, true), reduceddependencies(ls, os, true),
            Operation[], mref, childdependencies(ls, os, true)
        )
        push!(ls.operations, op)
        push!(opoffsets, opoffsets[end] + 1)
        return
    end
    nops = (nopsv[i])::Int # if it were a vector, it would have to have been expanded
    # if expanded, optyp must be either loopvalue, or compute (with loopvalues in its ancestry, not cutoff by loads)
    for offset = 0:nops-1
        sym = nops === 1 ? opsymbol : expandedopname(opsymbol, offset)
        op = Operation(
            length(operations(ls)), sym, elbytes, instr, optyp,
            loopdependencies(ls, os, false, offset), reduceddependencies(ls, os, false, offset),
            Operation[], mref, childdependencies(ls, os, false, offset)
        )
        push!(ls.operations, op)
    end
    push!(opoffsets, opoffsets[end] + nops)
    nothing
end
function add_parents_to_op!(ls::LoopSet, op::Operation, up₀::UInt128, up₁::UInt128, k::Int, Δ::Int)
    vparents = parents(op)
    ops = operations(ls)
    offsets = ls.operation_offsets
    if isone(Δ) # not expanded
        @assert isone(k)
        for i ∈ parents(ls, up₀, up₁)
            # FIXME; children also filled in cacheunrolled
            for j ∈ offsets[i]+1:offsets[i+1] # if parents are expanded, add them all
                opp = ops[j]
                pushfirst!(vparents, opp)
                push!(children(opp), op)
            end
        end
    else#if isexpanded
        # Do we want to require that all Δidxs are equal?
        # Because `CartesianIndex((2,3)) - 1` results in a methoderorr, I think this is reasonable for now
        # FIXME; children also filled in cacheunrolled
        for i ∈ parents(ls, up₀, up₁)
            opp = ops[offsets[i]+k]
            pushfirst!(vparents, opp)
            push!(children(opp), op)
        end
    end
end
function add_parents_to_ops!(ls::LoopSet, ops::Vector{OperationStruct}, constoffset)
    offsets = ls.operation_offsets
    for i in 1:length(offsets)-1
        pos = offsets[i]
        Δ = offsets[i+1]-pos
        for k ∈ 1:Δ
            op = ls.operations[pos+k]
            if isconstant(op)
                instr = instruction(op)
                if !skip_constant(instr)
                    constoffset += 1
                    pushpreamble!(ls, Expr(:(=), instr.instr, extract_varg(constoffset)))
                end
            elseif !isloopvalue(op)
                add_parents_to_op!(ls, op, ops[i].parents₀, ops[i].parents₁, k, Δ)
            end
        end
    end
    constoffset
end
function add_ops!(
    ls::LoopSet, instr::Vector{Instruction}, ops::Vector{OperationStruct}, mrefs::Vector{ArrayReferenceMeta}, elementbytes::Vector{Int},
    opsymbols::Vector{Symbol}, constoffset::Int, nopsv::Vector{NOpsType}, expandedv::Vector{Bool}
)
    # @show ls.loopsymbols ls.loopsymbol_offsets
    for i ∈ eachindex(ops)
        os = ops[i]
        opsymbol = opsymbols[os.symid]
        add_op!(ls, instr[i], ops, nopsv, expandedv, i, mrefs, opsymbol, elementbytes)
    end
    add_parents_to_ops!(ls, ops, constoffset)
    # for op ∈ operations(ls)
    #     if isstore(op) && isreduction(op) && iszero(length(loopdependencies(op)))
    #         addreduct_to_outer_reductions!(ls, op)
    #     end
    # end
    # for op in operations(ls)
        # @show op
    # end
end

# elbytes(::VectorizationBase.AbstractPointer{T}) where {T} = sizeof(T)::Int
# typeeltype(::Type{P}) where {T,P<:VectorizationBase.AbstractStridedPointer{T}} = T
typeeltype(::Type{Ptr{T}}) where {T} = T
# typeeltype(::Type{Core.LLVMPtr{T,0}}) where {T} = T
typeeltype(::Type{VectorizationBase.FastRange{T,F,S,O}}) where {T,F,S,O} = T
typeeltype(::Type{T}) where {T<:Real} = T
# typeeltype(::Any) = Int8

function add_array_symbols!(ls::LoopSet, arraysymbolinds::Vector{Symbol}, offset::Int)
    for as ∈ arraysymbolinds
        pushpreamble!(ls, Expr(:(=), as, extract_varg((offset+=1))))
    end
    return offset
end
function extract_external_functions!(ls::LoopSet, offset::Int, vargs)
    for op ∈ operations(ls)
        if iscompute(op)
            instr = instruction(op)
            if instr.mod != :LoopVectorization
                offset += 1
                instr_new = get(FUNCTIONSYMBOLS, vargs[offset], instr)
                if instr_new === instr
                    extractf = Expr(:call, GlobalRef(Core,:getfield), Symbol("#vargs#"), offset, false)
                    pushpreamble!(ls, Expr(:(=), instr.instr, extractf))
                else
                    op.instruction = instr_new
                end
            end
        end
    end
    offset
end
outer_reduct_init_typename(op::Operation) = Symbol(mangledvar(op), "#or#init#type#")
function extract_outerreduct_types!(ls::LoopSet, offset::Int, vargs)
  # for op
  for or ∈ ls.outer_reductions
    extractt = Expr(:call, GlobalRef(Core,:getfield), Symbol("#vargs#"), (offset+=1), false)
    pushpreamble!(ls, Expr(:(=), outer_reduct_init_typename(operations(ls)[or]), extractt))
  end
  offset
end
function sizeofeltypes(v)::Int
    num_arrays = length(v)::Int
    if num_arrays == 0
        return 8
    end
    T = typeeltype(v[1])
    sz = if (VectorizationBase.simd_integer_register_size() != VectorizationBase.register_size()) && T <: Integer # hack
        (VectorizationBase.register_size() ÷ VectorizationBase.simd_integer_register_size()) * sizeof(T)
    else
        sz = sizeof(T)
    end
    for i ∈ 2:num_arrays
        Ttemp = typeeltype(v[i])
        szᵢ = if (VectorizationBase.simd_integer_register_size() != VectorizationBase.register_size()) && T <: Integer # hack
            (VectorizationBase.register_size() ÷ VectorizationBase.simd_integer_register_size()) * sizeof(T)
        else
            sizeof(Ttemp)
        end
        # if !VectorizationBase.SIMD_NATIVE_INTEGERS && Ttemp <: Integer # hack
        #     return VectorizationBase.register_size()
        # end
        # T = promote_type(T, Ttemp)
        sz = max(szᵢ, sz)
    end
    sz
    # sizeof(T)
end

function avx_loopset!(
    ls::LoopSet, instr::Vector{Instruction}, ops::Vector{OperationStruct}, arf::Vector{ArrayRefStruct},
    AM::Vector{Any}, LPSYM::Vector{Any}, LB::Core.SimpleVector, vargs::Core.SimpleVector
)
    pushpreamble!(ls, :((var"#loop#bounds#", var"#vargs#") = var"#lv#tuple#args#"))
    add_loops!(ls, LPSYM, LB)
    resize!(ls.loop_order, ls.loopsymbol_offsets[end])
    arraysymbolinds = gen_array_syminds(AM)
    opsymbols = [gensym("op") for _ ∈ eachindex(ops)]
    nopsv = NOpsType[calcnops(ls, op) for op in ops]
    expandedv = [isexpanded(ls, ops, nopsv, i) for i ∈ eachindex(ops)]

    resize!(ls.loopindexesbit, length(ls.loops)); fill!(ls.loopindexesbit, false);
    mrefs, elementbytes = create_mrefs!(ls, arf, arraysymbolinds, opsymbols, nopsv, expandedv, vargs[1])
    for mref ∈ mrefs
        push!(ls.includedactualarrays, vptr(mref))
    end
    # extra args extraction
    extractind = add_ops!(ls, instr, ops, mrefs, elementbytes, opsymbols, 1, nopsv, expandedv)
    extractind = process_metadata!(ls, AM, extractind)
    extractind = add_array_symbols!(ls, arraysymbolinds, extractind)
    extractind = extract_external_functions!(ls, extractind, vargs)
    extractind = extract_outerreduct_types!(ls, extractind, vargs)
    ls
end
function avx_body(ls::LoopSet, UNROLL::Tuple{Bool,Int8,Int8,Int8,Bool,Int,Int,Int,Int,Int,Int,Int,UInt})
    inline, u₁, u₂, v, isbroadcast, W, rs, rc, cls, l1, l2, l3, nt = UNROLL
    q = (iszero(u₁) & iszero(v)) ? lower_and_split_loops(ls, inline % Int) : lower(ls, u₁ % Int, u₂ % Int, v % Int, inline % Int)
    ls.isbroadcast = isbroadcast
    iszero(length(ls.outer_reductions)) ? push!(q.args, nothing) : push!(q.args, loopset_return_value(ls, Val(true)))
    q
end

function _turbo_loopset_debug(::Val{UNROLL}, ::Val{OPS}, ::Val{ARF}, ::Val{AM}, ::Val{LPSYM}, _vargs::Tuple{LB,V}) where {UNROLL, OPS, ARF, AM, LPSYM, LB, V}
    # @show OPS ARF AM LPSYM _vargs
    _turbo_loopset(OPS, ARF, AM, LPSYM, _vargs[1].parameters, V.parameters, UNROLL)
end
function tovector(@nospecialize(t))
    v = Vector{Any}(undef, length(t))
    for i ∈ eachindex(v)
        tᵢ = t[i]
        if tᵢ isa Tuple # reduce specialization?
            v[i] = tovector(tᵢ)
        else
            v[i] = tᵢ
        end
    end
    v
end
function _turbo_loopset(
  @nospecialize(OPSsv), @nospecialize(ARFsv), @nospecialize(AMsv), @nospecialize(LPSYMsv), LBsv::Core.SimpleVector, vargs::Core.SimpleVector,
  UNROLL::Tuple{Bool,Int8,Int8,Int8,Bool,Int,Int,Int,Int,Int,Int,Int,UInt}
)
    nops = length(OPSsv) ÷ 3
    instr = Instruction[Instruction(OPSsv[3i+1], OPSsv[3i+2]) for i ∈ 0:nops-1]
    ops = OperationStruct[ OPSsv[3i] for i ∈ 1:nops ]
    ls = LoopSet(:LoopVectorization)
    inline, u₁, u₂, v, isbroadcast, W, rs, rc, cls, l1, l2, l3, nt = UNROLL
    set_hw!(ls, rs, rc, cls, l1, l2, l3); ls.vector_width = W; ls.isbroadcast = isbroadcast
    arsv = Vector{ArrayRefStruct}(undef, length(ARFsv))
    for i ∈ eachindex(arsv)
        arsv[i] = ARFsv[i]
    end
    avx_loopset!(ls, instr, ops, arsv, tovector(AMsv), tovector(LPSYMsv), LBsv, vargs)
end

@static if VERSION ≥ v"1.7.0-DEV.421"
  using Base: @aggressive_constprop
else
  macro aggressive_constprop(ex); esc(ex); end
end

"""
    _turbo_!(unroll, ops, arf, am, lpsym, lb, vargs...)

Execute an `@turbo` block. The block's code is represented via the arguments:
- `unroll` is `Val((u₁,u₂))` and specifies the loop unrolling factor(s).
  These values may be supplied manually via the `unroll` keyword
  of [`@turbo`](@ref).
- `ops` is `Tuple{mod1, sym1, op1, mod2, sym2, op2...}` encoding the operations of the loop.
  `mod` and `sym` encode the module and symbol of the called function; `op` is an [`OperationStruct`](@ref)
  encoding the details of the operation.
- `arf` is `Tuple{arf1, arf2...}`, where each `arfi` is an [`ArrayRefStruct`](@ref) encoding
  an array reference.
- `am` contains miscellaneous data about the LoopSet (see `process_metadata!`)
- `lpsym` is `Tuple{:i,:j,...}`, a Tuple of the "loop symbols", i.e. the item variable `i` in `for i ∈ iter`
- `lb` is `Tuple{RngTypei,RngTypej,...}`, a Tuple encoding syntactically-knowable information about
  the iterators corresponding to `lpsym`. For example, in `for i ∈ 1:n`, the `1:n` would be encoded with
  `StaticLowerUnitRange(1)` because the lower bound of the iterator can be determined to be 1.
- `vargs...` holds the encoded pointers of all the arrays (see `VectorizationBase`'s various pointer types).
"""
@aggressive_constprop @generated function _turbo_!(
    ::Val{var"#UNROLL#"}, ::Val{var"#OPS#"}, ::Val{var"#ARF#"}, ::Val{var"#AM#"}, ::Val{var"#LPSYM#"}, ::Val{Tuple{var"#LB#",var"#V#"}}, var"#flattened#var#arguments#"::Vararg{Any,var"#num#vargs#"}
) where {var"#UNROLL#", var"#OPS#", var"#ARF#", var"#AM#", var"#LPSYM#", var"#LB#", var"#V#", var"#num#vargs#"}
  # 1 + 1 # Irrelevant line you can comment out/in to force recompilation...
  ls = _turbo_loopset(var"#OPS#", var"#ARF#", var"#AM#", var"#LPSYM#", var"#LB#".parameters, var"#V#".parameters, var"#UNROLL#")
  pushfirst!(ls.preamble.args, :(var"#lv#tuple#args#" = reassemble_tuple(Tuple{var"#LB#",var"#V#"}, var"#flattened#var#arguments#")))
  # return @show avx_body(ls, var"#UNROLL#")
  if last(var"#UNROLL#") > 1
    inline, u₁, u₂, v, isbroadcast, W, rs, rc, cls, l1, l2, l3, nt = var"#UNROLL#"
    # wrap in `var"#OPS#", var"#ARF#", var"#AM#", var"#LPSYM#"` in `Expr` to homogenize types
    avx_threads_expr(
      ls, (inline, u₁, u₂, v, isbroadcast, W, rs, rc, cls, l1, l2, l3, one(UInt)), nt,
      :(Val{$(var"#OPS#")}()), :(Val{$(var"#ARF#")}()), :(Val{$(var"#AM#")}()), :(Val{$(var"#LPSYM#")}())
    )
  else
    # Main.BODY[] = avx_body(ls, var"#UNROLL#")
    # @show avx_body(ls, var"#UNROLL#")
    avx_body(ls, var"#UNROLL#")
  end
  # @show var"#UNROLL#", var"#OPS#", var"#ARF#", var"#AM#", var"#LPSYM#", var"#LB#"
end


text:
"""
    Interval{T<:Real} <: Domain{T}

An Interval of type 'T' between 'a' and 'b' represents all the values of type 'T' between 'a' and 'b'.
"""
struct Interval{T<:Real} <: Domain{T}
  l::T
  u::T
  Interval{T}(l,u) where {T<:Real} = if u > l new(l, u) else new(u,l) end
end

Interval(x::Real) = Interval{typeof(x)}(x,x)
Interval(v::Vector{T}) where T <: Real = Interval(v[1],v[2])
Interval(x::T,y::T) where T <: Real = Interval{T}(x,y)
Interval(x::T1,y::T2) where {T1<:Real, T2<:Real} = Interval{promote_type(T1,T2)}(promote(x,y)...)

## Conversions and Promotion
## =========================

# A concrete number can be coerced into an interval with no width
function convert(::Type{Interval{T1}}, x::Interval{T2}) where {T1<:Real, T2<:Real}
  T = promote_type(T1,T2)
  Interval{T}(convert(T,x.l),convert(T,x.u))
end
convert(::Type{Interval}, c::Real) = Interval{typeof(c)}(c,c)
convert(::Type{Interval{T}}, c::T) where {T<:Real} = Interval{T}(c,c)
convert(::Type{Interval{T1}}, c::T2) where {T1<:Real, T2<:Real} = Interval{T1}(c,c)

promote_rule(::Type{Interval{T1}}, ::Type{T2}) where {T1<:Real, T2<:Real} = Interval{T1}
promote_rule(::Type{Interval{T1}}, ::Type{Interval{T2}}) where {T1<:Real, T2<:Real} = Interval{promote_type(T1,T2)}

## Domain operations
## =================
ndims(i::Interval) = 1
subsumes(x::Interval, y::Interval) = y.l >= x.l && y.u <= x.u
issubset(x::Interval, y::Interval) = x.l >= y.l && x.u <= y.u
⊑ = issubset

isintersect(x::Interval, y::Interval) = y.l <= x.u && x.l <= y.u
domaineq(x::Interval, y::Interval) = x.u == y.u && x.l == y.l

"""
    intersect(x::Interval{T}, y::Interval{T}) where T

Construct interval which is intersection of two intervals
"""
intersect(x::Interval{T}, y::Interval{T}) where T = Interval(max(x.l, y.l), min(x.u, y.u))
intersect(a::Interval{T}, b::Interval{S}) where {T,S} = intersect(promote(a,b)...)
⊓ = intersect

## Union/Join
union(a::Interval, b::Interval) = ⊔(a, b)
function ⊔(a::Interval, b::Interval)
  l = min(a.l,b.l)
  u = max(a.u, b.u)
  Interval(l,u)
end

⊔(a::Interval, b::Real) = ⊔(promote(a,b)...)
⊔(b::Real, a::Interval) = ⊔(promote(b,a)...)
⊔(a::Interval) = a
⊔(a::Vector{Interval}) = reduce(⊔,a)

isequal(x::Interval,y::Interval) = domaineq(x,y)
isrelational(::Type{Interval}) = false

isabstract(c::Type{T}, a::Type{Interval{T}}) where {T<:Real} = true

## Interval Arithmetic and Inequalities
## ====================================

# ==, != return values in AbstractBool
function ==(x::Interval, y::Interval)
  if x.u == y.u == x.l == y.l t
  elseif isintersect(x,y) tf
  else f end
end

!=(x::Interval,y::Interval) = !(==(x,y))

==(x::Interval,y::Real) = ==(promote(x,y)...)
==(y::Real,x::Interval) = ==(promote(y,x)...)

!=(x::Interval, y::Real) = !==(x,y)
!=(y::Real, x::Interval) = !==(y,x)

>(x::Interval, y::Interval) = if x.l > y.u t elseif x.u <= y.l f else tf end
>(x::Interval, y::Real) = if x.l > y t elseif x.u <= y f else tf end
>(y::Real, x::Interval) =  if y > x.u t elseif y <= x.l f else tf end

<(x::Interval, y::Interval) = y > x
<(x::Interval, y::Real) = y > x
<(y::Real, x::Interval) = x > y

<=(x::Interval, y::Interval) = !(x > y)
>=(x::Interval, y::Interval) = !(x < y)
<=(x::Interval, y::Real) = !(x > y)
<=(y::Real, x::Interval) = !(y > x)

>=(x::Interval, y::Real) = !(x < y)
>=(y::Real, x::Interval) = !(y < x)

+(x::Interval, y::Interval) = Interval(x.l + y.l, x.u + y.u)
-(x::Interval, y::Interval) = Interval(x.l - y.u, x.u - y.l)
+(x::Interval, y::Real) = Interval(x.l + y, x.u + y)
+(y::Real, x::Interval) = x + y
+(x::Interval) = x
-(x::Interval, y::Real) = Interval(x.l - y, x.u - y)
-(y::Real, x::Interval) = Interval(y - x.l, y - x.u)
-(x::Interval{T}) where {T} = zero{T} - x
*(x::Interval, y::Real) = Interval(x.l * y, x.u * y)
*(y::Real, x::Interval) = x * y

sqrt(x::Interval) = Interval(sqrt(x.l), sqrt(x.u))

# CODEREVIEW: Generalise to even powers
function sqr(x::Interval)
  a,b,c,d = x.l * x.l, x.l * x.u, x.u * x.l, x.u * x.u
  Interval(max(min(a,b,c,d),0),max(a,b,c,d,0))
end

function *(x::Interval, y::Interval)
  a,b,c,d = x.l * y.l, x.l * y.u, x.u * y.l, x.u * y.u
  Interval(min(a,b,c,d),max(a,b,c,d))
end

# is c inside the interval
# CODREVIEW: TESTME
in(c::Real, y::Interval) = y.l <= c <= y.u

# CODREVIEW: TESTME
inv(x::Interval) = Interval(1/x.u,1/x.l)

# Ratz Interval Division
# CODREVIEW: TESTME
function /(x::Interval, y::Interval)
  a,b,c,d = x.l,x.u,y.l,y.u
  if !(0 ∈ y)
    x * inv(y)
  elseif (0 ∈ x)
    Interval(-Inf,Inf)
  elseif b < 0 && c < d == 0
    Interval(b/c,Inf)
  elseif b < 0 && c < 0 < d
    Interval(-Inf,Inf)
  elseif b < 0 && 0 == c < d
    Interval(-Inf,b/d)
  elseif 0 < a && c < d == 0
    Interval(-Inf,a/c)
  elseif 0 < a && c < 0 < d
    Interval(-Inf,Inf)
  elseif 0 < a && 0 == c < d
    Interval(a/d, Inf)
  else
    Inf
  end
end

/(c::Real, x::Interval) = convert(Interval,c) / x
/(x::Interval, c::Real) = x / convert(Interval,c)
## Rationals
//(x::Interval, y::Interval) = x / y
//(x::Interval, c::Real) = x / c
//(c::Real, x::Interval) = c / x

## Functions on Interval type
## ==========================

unit(::Type{Interval{T}}) where T = Interval{T}(zero(T), one(T))
unit(::Interval{T}) where T = Interval{T}(zero(T), one(T))

## It's all Ones and Zeros
zero(::Type{Interval}) = Interval(0.0,0.0)
one(::Type{Interval{T}}) where T = Interval(one(T))
one(::Interval{T}) where T = Interval(one(T))
zero(::Type{Interval{T}}) where T = Interval(zero(T))
zero(::Interval{T}) where T = Interval(zero(T))

## Functions on interval abstraction itself
## =======================================
reflect(x::Interval) = Interval(-x.l,-x.u)
makepos(x::Interval) = Interval(max(x.l,0), max(x.u,0))
mid(x::Interval) = (x.u - x.l) / 2 + x.l

## Non primitive functions
## =======================
function abs(x::Interval)
  if x.l >= 0.0 && x.u >= 0.0 x
  elseif x.u >= 0.0 Interval(0,max(abs(x.l), abs(x.u)))
  else makepos(reflect(x))
  end
end

round(x::Interval) = Interval(round(x.l), round(x.u))

function isinf(x::Interval)
  if isinf(x.l) || isinf(x.u)
    x.u == x.l ? T : t
  else
    F
  end
end

function isapprox(x::Interval, y::Interval; epsilon::Real = 1E-5)
  ifelse(isinf(x) | isinf(y), x == y, abs(x - y) <= epsilon)
end

isapprox(x::Interval, y::Real) = isapprox(promote(x,y)...)
isapprox(x::Real, y::Interval) = isapprox(promote(x,y)...)

## Vector Interop
## ==============
l(v::Vector{T}) where {T<:Real} = v[1]
u(v::Vector{T}) where {T<:Real} = v[2]
l(x::Interval) = x.l
u(x::Interval) = x.u
pair(::Type{Interval{T}},low,up) where T = Interval(low,up)
pair(::Type{Vector{Float64}},low,up) = [low,up]
Pair = Union{Vector{Float64},Interval}

## Splitting
## =========
function split_box(i::P, split_point::Float64) where {P<:Pair}
  @assert l(i) <= split_point <= u(i) "Split point must be within interval"
  # @assert l(i) != u(i) "Can't split a single point interval into disjoint sets"

  if l(i) == u(i)  #Degenrate case
    P[pair(P, l(i), u(i)), pair(P, l(i), u(i))]
  elseif split_point < u(i)
    P[pair(P, l(i), split_point), pair(P, nextfloat(split_point), u(i))]
  else
    P[pair(P, l(i), prevfloat(split_point)), pair(P, split_point, u(i))]
  end
end

mid_split(i::Interval) = split_box(i,mid(i))

# Split along the middle n times
function mid_split(i::Interval, n::Int64)
  A = [i]
  for i = 1:n
    res = Interval[]
    for a in A
      splitted = mid_split(a)
      push!(res,splitted[1],splitted[2])
    end
    A = res
  end
  A
end

## Sampling
## ========
rand(x::Interval{T}) where {T<:AbstractFloat}  =  x.l + (x.u - x.l) * rand(T)
rand(x::Interval{T}) where {T<:Integer} = rand(UnitRange(x.l,x.u))
rand(x::Interval{T},n::Int) where {T<:Real} = T[rand(x) for i = 1:n]

## Print
## =====
string(x::Interval) = "[$(x.l) $(x.u)]"
print(io::IO, x::Interval) = print(io, string(x))
show(io::IO, x::Interval) = print(io, string(x))

ext: jl


text:
abstract type Equation end

struct DiagnosticEquation <: Equation
    var::Symbol
    terms::Tuple
end

struct PrognosticEquation{T<:Function} <: Equation
    var::Symbol
    inc::Symbol
    next::Symbol
    terms::Tuple
    integrator::T

    function PrognosticEquation(var::Symbol, terms::Tuple, integrator::T) where {T<:Function}
        new{T}(var, var, np1_symbol(var), terms, integrator)
    end
end

is_diagnostic(eq) = false
is_diagnostic(eq::DiagnosticEquation) = true
is_prognostic(eq) = false
is_prognostic(eq::PrognosticEquation) = true


text:
myfunc(x) = log(sum(x))

@testset "others" begin
    test_compare(myfunc; x=rand(2))
end


text:
    # This proposal samples from prior and updates exogenous variables
# Using inverse transform

# function stdproposal(qω, x::T, ω) where T
#   display(ω)
#   @show T
#   if x in keys(ω.data)
#     nothing
#   else
#     x_ = rand(qω, x.class)
#     @show x => x_
#   end
#   @show T
#   @assert false
# end

stdproposal(qω, x, ω) = nothing

# function stdproposal(qω, x::Member{<:Distribution}, ω)
#   @show x
#   display(ω)
#   @show ω.data[x]
#   if x in keys(ω)
#     @show "here"
#     @assert false
#   else
#     @show "not here"
#     @assert false
#   end
# end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# push!(LOAD_PATH, abspath("src"))
# push!(LOAD_PATH, abspath("../MPSCircuit.jl/src"))
using Yao, MSQR, MPSCircuit

# Basic parameters.
nBitT = 4
vBit = 2
rBit = 1
depth = 2
nMeasure = 2000
regTar = rand_state(nBitT)
MPSGen = MPSC(("DC",depth),nBitT,vBit,rBit)
circuit = MScircuit(nBitT, vBit, rBit, MPSGen.mpsBlocks)

# Setup Training environments
## Gradient Optimization using ADAM algorithm(https://arxiv.org/abs/1412.6980v8).
method1 = ("ADAM", 0.05)
nTrain1 = 20
## Gradient Optimization using fixed step size(default size=0.1).
method2 = ("default", 0.05)
nTrain2 = 20

# Training Program.
## ADAM based SGD method.
### MSQR training.
MSQRtrain!(regTar, deepcopy(circuit), nTrain1, nMeasure=nMeasure, GDmethod = method1, show=true)
### SWAP Test + SGD training.
SWAPtrain!(regTar, deepcopy(MPSGen.cExtend), nTrain1, nMeasure=nMeasure, GDmethod = method1, show=true)

## default SGD method.
### MSQR training.
MSQRtrain!(regTar, deepcopy(circuit), nTrain2, nMeasure=nMeasure, GDmethod = method2, show=true)
### SWAP Test + SGD training.
SWAPtrain!(regTar, deepcopy(MPSGen.cExtend), nTrain2, nMeasure=nMeasure, GDmethod = method2, show=true)


text:
"""
Extends the getindex method to access symbolic variables from the Community structure.

# Arguments
 - **a** (Community) Community structure
 - **var** (Symbol) Symbol to be extracted

# Returns
Float if symbol is global parameter
Array{Float} otherwise
"""
function Base.getproperty(a::Community,var::Symbol)

    if var in fieldnames(Community)
        return getfield(a,var)
    else
        if var in a.declaredSymbols_["Local"]
            pos = findfirst(a.declaredSymbols_["Local"].==var)
            return @views a.local_[:,pos]
        elseif var in a.declaredSymbols_["Identity"]
            pos = findfirst(a.declaredSymbols_["Identity"].==var)
            return @views a.identity_[:,pos]
        elseif var in a.declaredSymbols_["Global"]
            pos = findfirst(a.declaredSymbols_["Global"].==var)
            return @views a.global_[pos]
        elseif var in a.declaredSymbols_["GlobalArray"]
            pos = findfirst(a.declaredSymbols_["GlobalArray"].==var)
            return a.globalArray_[pos]
        elseif var in a.declaredSymbols_["Medium"]
            pos = findfirst(a.declaredSymbols_["Medium"].==var)
            if a.dims == 1
                return @views a.medium_[:,pos]
            elseif a.dims == 2
                return @views a.medium_[:,:,pos]
            elseif a.dims == 3
                println(size(a.medium_))
                return @views a.medium_[:,:,:,pos]
            end
        else
            error("Parameter ", var, " not found in the community.")
        end
    end
end

"""
Extends the setindex method to assign symbolic variables from the Community structure.

# Arguments
 - **a** (Community) Community structure
 - **v** (Float, Array{Float}) Parameters to be associated with the parameter. If Array, it has to be of the same length as the number of agents in the Community.
 - **var** (Symbol) Symbol to be associated

# Returns
Float if symbol is global parameter
Array{Float} otherwise
"""
function Base.setproperty!(a::Community,var::Symbol,v::Array{<:Number})

    if var in a.declaredSymbols_["Local"]
        if size(v) != (a.N,)
            error("Trying to assign array with shape ", shape(v), ". It should be of size (N,)")
        end
        pos = findfirst(a.declaredSymbols_["Local"].==var)
        vec = a.local_
        vec[:,pos] = v
        setfield!(a,:local_,vec)
    elseif var in a.declaredSymbols_["Identity"]
        if size(v) != (a.N,)
            error("Trying to assign array with shape ", shape(v), ". It should be of size (N,)")
        end
        pos = findfirst(a.declaredSymbols_["Identity"].==var)
        vec = a.identity_
        vec[:,pos] = v
        setfield!(a,:identity_,vec)
    elseif var in a.declaredSymbols_["GlobalArray"]
        pos = findfirst(a.declaredSymbols_["GlobalArray"].==var)
        a.globalArray_[pos] = v
    elseif var in a.declaredSymbols_["Medium"]
        if size(v) != size(a.var)[1:end-1]
            error("Trying to assign array with shape ", shape(v), ". It should be of size ", size(a.var)[1:end-1])
        end
        pos = findfirst(a.declaredSymbols_["Identity"].==var)
        vec = a.medium_
        if a.dims == 1
            vec[:,pos] = v
        elseif a.dims == 2
            vec[:,:,pos] = v
        elseif a.dims == 3
            vec[:,:,:,pos] = v
        end
        setfield!(a,:medium_,vec)
    else
        error("Parameter ", var, " not fount in the community.")
    end

end

function Base.setproperty!(a::Community,var::Symbol,v::Number)

    if var in a.declaredSymbols_["Local"]
        error("Local parameter ", var, " cannot be assigned with a scalar. Use .= instead.")
    elseif var in a.declaredSymbols_["Identity"]
        error("Identity parameter ", var, " cannot be assigned with a scalar. Use .= instead.")
    elseif var in a.declaredSymbols_["Global"]
        pos = findfirst(a.declaredSymbols_["Global"].==var)
        a.global_[pos] = v
    elseif var in a.declaredSymbols_["GlobalArray"]
        error("GlobalArray ", var, " cannot be assigned with a scalar.")
    elseif var == :N
        setfield!(a,:N,v)
    elseif var == :t
        setfield!(a,:t,v)
    else
        error("Parameter ", var," not fount in the community.")
    end

end

function Base.push!(a::CommunityInTime,c::Community)

    push!(a.com,c)

    return
end

function Base.length(a::CommunityInTime)

    return length(a.com)
end

function Base.getproperty(a::CommunityInTime,var::Symbol)

    if var in fieldnames(CommunityInTime)
        return getfield(a,var)
    else

        if var == :t
            return [i.t for i in a.com]
        elseif var == :N
            return [i.N for i in a.com]
        end

        if var in a.com[1].declaredSymbols_["GlobalArray"]

            out = []
            for i in 1:length(a)
                push!(out,getproperty(a.com[i],var))
            end

            return out

        else

            nMax = maximum(a.com[end].id)
            #Create array
            out = fill(NaN,length(a),nMax)
            #Fill array
            for i in 1:length(a)
                a.com[i].id
                out[i,a[i].id] .= getproperty(a.com[i],var)
            end

            return out

        end

    end

    return
end

function Base.getindex(a::CommunityInTime,var::Int)

    return a.com[var]

end

Base.first(a::CommunityInTime) = 1
Base.lastindex(a::CommunityInTime) = length(a)

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
# mpg123 wrappers of required functions to read mp3 files

# mpg123 encodings
const MPG123_ENC_8           = 0x00f                                            # 0000 0000 0000 1111 Some 8 bit  integer encoding.
const MPG123_ENC_16          = 0x040                                            # 0000 0000 0100 0000 Some 16 bit integer encoding.
const MPG123_ENC_24          = 0x4000                                           # 0100 0000 0000 0000 Some 24 bit integer encoding.
const MPG123_ENC_32          = 0x100                                            # 0000 0001 0000 0000 Some 32 bit integer encoding.
const MPG123_ENC_SIGNED      = 0x080                                            # 0000 0000 1000 0000 Some signed integer encoding.
const MPG123_ENC_FLOAT       = 0xe00                                            # 0000 1110 0000 0000 Some float encoding.
const MPG123_ENC_SIGNED_16   = MPG123_ENC_16 | MPG123_ENC_SIGNED | 0x10         # 0000 0000 1101 0000 signed 16 bit
const MPG123_ENC_UNSIGNED_16 = MPG123_ENC_16 | 0x20                             # 0000 0000 0110 0000 unsigned 16 bit
const MPG123_ENC_UNSIGNED_8  = 0x01                                             # 0000 0000 0000 0001 unsigned 8 bit
const MPG123_ENC_SIGNED_8    = MPG123_ENC_SIGNED | 0x02                         # 0000 0000 1000 0010 signed 8 bit
const MPG123_ENC_ULAW_8      = 0x04                                             # 0000 0000 0000 0100 ulaw 8 bit
const MPG123_ENC_ALAW_8      = 0x08                                             # 0000 0000 0000 0100 alaw 8 bit
const MPG123_ENC_SIGNED_32   = MPG123_ENC_32 | MPG123_ENC_SIGNED | 0x1000       # 0001 0001 1000 0000 signed 32 bit
const MPG123_ENC_UNSIGNED_32 = MPG123_ENC_32 | 0x2000                           # 0010 0001 0000 0000 unsigned 32 bit
const MPG123_ENC_SIGNED_24   = MPG123_ENC_24 | MPG123_ENC_SIGNED | 0x1000       # 0101 0000 1000 0000 signed 24 bit
const MPG123_ENC_UNSIGNED_24 = MPG123_ENC_24 | 0x2000                           # 0110 0000 0000 0000 unsigned 24 bit
const MPG123_ENC_FLOAT_32    = 0x200                                            # 0000 0010 0000 0000 32bit float
const MPG123_ENC_FLOAT_64    = 0x400                                            # 0000 0100 0000 0000 64bit float

# Any possibly known encoding from the list above.
const MPG123_ENC_ANY = ( MPG123_ENC_SIGNED_16  | MPG123_ENC_UNSIGNED_16
	                   | MPG123_ENC_UNSIGNED_8 | MPG123_ENC_SIGNED_8
	                   | MPG123_ENC_ULAW_8     | MPG123_ENC_ALAW_8
	                   | MPG123_ENC_SIGNED_32  | MPG123_ENC_UNSIGNED_32
	                   | MPG123_ENC_SIGNED_24  | MPG123_ENC_UNSIGNED_24
	                   | MPG123_ENC_FLOAT_32   | MPG123_ENC_FLOAT_64    )

"""represents the C pointer mpg123_handle*. used by all mpg123 functions"""
typealias MPG123 Ptr{Void}

const MPG123_DONE               = -12
const MPG123_NEW_FORMAT         = -11
const MPG123_NEED_MORE          = -10
const MPG123_ERR                = -1
const MPG123_OK                 = 0

"""return a string that explains given error code"""
function mpg123_plain_strerror(err)
    str = ccall((:mpg123_plain_strerror, libmpg123), Ptr{Cchar}, (Cint,), err)
    bytestring(str)
end

"""initialize mpg123 library"""
function mpg123_init()
    err = ccall((:mpg123_init, libmpg123), Cint, ())
    if err != MPG123_OK
        error("Could not initialize mpg123: ", mpg123_plain_strerror(err))
    end
end

"""create new mpg123 handle"""
function mpg123_new()
    err = Ref{Cint}(0)
    mpg123 = ccall((:mpg123_new, libmpg123), MPG123,
                   (Ptr{Cchar}, Ref{Cint}),
                   C_NULL, err)

    if err.x != MPG123_OK
        error("Could not create mpg123 handle: ", mpg123_plain_strerror(err.x))
    end

    mpg123
end

"""open an mp3 file at fiven path"""
function mpg123_open(mpg123::MPG123, path::AbstractString)
    err = ccall((:mpg123_open, libmpg123), Cint,
                (MPG123, Ptr{Cchar}),
                mpg123, path)

    if err != MPG123_OK
        mpg123_delete(mpg123)
        error("Could not open $path: ", mpg123_plain_strerror(err))
    end

    err
end

"""close a file that is opened by given handle"""
function mpg123_close(mpg123::MPG123)
    err = ccall((:mpg123_close, libmpg123), Cint, (MPG123,), mpg123)

    if err != MPG123_OK
        warn("Could not close mpg123 $mpg123: ", mpg123_plain_strerror(err))
    end

    err
end

"""delete mpg123 handle"""
function mpg123_delete(mpg123::MPG123)
    ccall((:mpg123_delete, libmpg123), Cint, (MPG123,), mpg123)
end

"""return birtate, number of channels and encoding of the mp3 file"""
function mpg123_getformat(mpg123::MPG123)
    bitrate = Ref{Clong}(0)
    nchannels = Ref{Cint}(0)
    encoding = Ref{Cint}(0)
    err = ccall((:mpg123_getformat, libmpg123), Cint,
                (MPG123, Ref{Clong}, Ref{Cint}, Ref{Cint}),
                mpg123, bitrate, nchannels, encoding)

    if err != MPG123_OK
        error("Could not read format: ", mpg123_plain_strerror(err))
    end

    bitrate.x, nchannels.x, encoding.x
end

"""return the appropriate block size for handling this mpg123 handle"""
function mpg123_outblock(mpg123::MPG123)
    ccall((:mpg123_outblock, libmpg123), Csize_t, (MPG123,), mpg123)
end

"""return the number of samples in the file"""
function mpg123_length(mpg123::MPG123)
    length = ccall((:mpg123_length, libmpg123), Coff_t, (MPG123,), mpg123)
    if length == MPG123_ERR
        error("Could not determine the frame length")
    end
    convert(Int64, length)
end

"""return how many bytes a sample (in a channel) uses"""
function mpg123_encsize(encoding::Cint)
    ccall((:mpg123_encsize, libmpg123), Cint, (Cint,), encoding)
end

"""
read audio samples from the mpg123 handle

# Arguments
* `mpg123::MPG123`: the mpg123 handle
* `out::Array{T}`: Array with appropriate data type, to store the samples
* `size::Integer`: the amount to read, in bytes. nchannels * encsize * nsamples
"""
function mpg123_read!{T}(mpg123::MPG123, out::Array{T}, size::Integer)
    done = Ref{Csize_t}(0)
    err = ccall((:mpg123_read, libmpg123), Cint,
                (MPG123, Ptr{T}, Csize_t, Ref{Csize_t}),
                mpg123, out, size, done)

    if err != MPG123_OK && err != MPG123_DONE
        error("Error while reading $mpg123: ", mpg123_plain_strerror(err))
    end

    Int(done.x)
end


text:
# VE pointer types

export VEPtr, PtrOrVEPtr, VEArrayPtr, VE_NULL, VERef, RefOrVERef

#
# VE device pointer
#

# FIXME: should be called VEDevicePtr...

"""
    VEPtr{T}

A memory address that refers to data of type `T` that is accessible from the VE. A `VEPtr`
is ABI compatible with regular `Ptr` objects, e.g. it can be used to `ccall` a function that
expects a `Ptr` to VE memory, but it prevents erroneous conversions between the two.
"""
VEPtr

primitive type VEPtr{T} 64 end

# constructor
VEPtr{T}(x::Union{Int,UInt,VEPtr}) where {T} = Base.bitcast(VEPtr{T}, x)

const VE_NULL = VEPtr{Cvoid}(0)

## getters

Base.eltype(::Type{<:VEPtr{T}}) where {T} = T

## conversions

# to and from integers
## pointer to integer
Base.convert(::Type{T}, x::VEPtr) where {T<:Integer} = T(UInt(x))
## integer to pointer
Base.convert(::Type{VEPtr{T}}, x::Union{Int,UInt}) where {T} = VEPtr{T}(x)
Int(x::VEPtr)  = Base.bitcast(Int, x)
UInt(x::VEPtr) = Base.bitcast(UInt, x)

# between regular and VE pointers
Base.convert(::Type{Ptr{T}}, p::VEPtr) where {T} =
    throw(ArgumentError("cannot convert from device to host pointer"))
#Base.unsafe_convert(::Type{Ptr{T}}, x::VEPtr) where {T} = convert(Ptr{T}, x)
Base.convert(::Type{VEPtr{T}}, p::Ptr{T}) where {T} = # Base.bitcast(Ptr{T}, p)
    throw(ArgumentError("cannot convert from host to device pointer"))
#Base.unsafe_convert(::Type{VEPtr{T}}, x::Ptr{T}) where {T} = convert(VEPtr{T}, x)

# between VE pointers
Base.convert(::Type{VEPtr{T}}, p::VEPtr) where {T} = Base.bitcast(VEPtr{T}, p)

# defer conversions to unsafe_convert
Base.cconvert(::Type{<:VEPtr}, x) = x

# fallback for unsafe_convert
Base.unsafe_convert(::Type{P}, x::VEPtr) where {P<:VEPtr} = convert(P, x)

## limited pointer arithmetic & comparison

Base.isequal(x::VEPtr, y::VEPtr) = (x === y)
Base.isless(x::VEPtr{T}, y::VEPtr{T}) where {T} = x < y

Base.:(==)(x::VEPtr, y::VEPtr) = UInt(x) == UInt(y)
Base.:(<)(x::VEPtr,  y::VEPtr) = UInt(x) < UInt(y)
Base.:(-)(x::VEPtr,  y::VEPtr) = UInt(x) - UInt(y)

Base.:(+)(x::VEPtr, y::Integer) = oftype(x, Base.add_ptr(UInt(x), (y % UInt) % UInt))
Base.:(-)(x::VEPtr, y::Integer) = oftype(x, Base.sub_ptr(UInt(x), (y % UInt) % UInt))
Base.:(+)(x::Integer, y::VEPtr) = y + x

#
# Host or device pointer
#

"""
    PtrOrVEPtr{T}

A special pointer type, ABI-compatible with both `Ptr` and `VEPtr`, for use in `ccall`
expressions to convert values to either a VE or a CPU type (in that order). This is
required for VEDA APIs which accept pointers that either point to host or device memory.
"""
PtrOrVEPtr

primitive type PtrOrVEPtr{T} 64 end

function Base.cconvert(::Type{PtrOrVEPtr{T}}, val) where {T}
    # `cconvert` is always implemented for both `Ptr` and `VEPtr`, so pick the first result
    # that has done an actual conversion

    ve_val = Base.cconvert(VEPtr{T}, val)
    if ve_val !== val
        return ve_val
    end

    cpu_val = Base.cconvert(Ptr{T}, val)
    if cpu_val !== val
        return cpu_val
    end

    return val
end

function Base.unsafe_convert(::Type{PtrOrVEPtr{T}}, val) where {T}
    ptr = if Core.Compiler.return_type(Base.unsafe_convert,
                                       Tuple{Type{Ptr{T}}, typeof(val)}) !== Union{}
        Base.unsafe_convert(Ptr{T}, val)
    elseif Core.Compiler.return_type(Base.unsafe_convert,
                                     Tuple{Type{VEPtr{T}}, typeof(val)}) !== Union{}
        Base.unsafe_convert(VEPtr{T}, val)
    else
        throw(ArgumentError("cannot convert to either a host or device pointer"))
    end

    return Base.bitcast(PtrOrVEPtr{T}, ptr)
end

#
# Device reference objects
#

primitive type VERef{T} 64 end

# general methods for VERef{T} type
Base.eltype(x::Type{<:VERef{T}}) where {T} = @isdefined(T) ? T : Any

Base.convert(::Type{VERef{T}}, x::VERef{T}) where {T} = x

# conversion for the actual ccall
Base.unsafe_convert(::Type{VERef{T}}, x::VERef{T}) where {T} = Base.bitcast(VERef{T}, Base.unsafe_convert(VEPtr{T}, x))
Base.unsafe_convert(::Type{VERef{T}}, x) where {T} = Base.bitcast(VERef{T}, Base.unsafe_convert(VEPtr{T}, x))

Base.unsafe_convert(::Type{VERef{T}}, x::Base.RefValue{T}) where {T} =
    Base.bitcast(VERef{T}, Base.unsafe_convert(Ptr{T}, x))

# VERef from literal pointer
Base.convert(::Type{VERef{T}}, x::VEPtr{T}) where {T} = x

# indirect constructors using VERef
Base.convert(::Type{VERef{T}}, x) where {T} = VERef{T}(x)

# dereferencing
@inline Base.getindex(p::VERef{T}) where {T} = Base.bitcast(Ref{T}, p)[]
@inline function Base.setindex(p::VERef{T}, x::T) where {T}
    r = Base.bitcast(Ref{T}, p)
    r[] = x
    x
end

#
# VE array pointer
#

primitive type VEArrayPtr{T} 64 end

# constructor
VEArrayPtr{T}(x::Union{Int,UInt,VEArrayPtr}) where {T} = Base.bitcast(VEArrayPtr{T}, x)

## getters

Base.eltype(::Type{<:VEArrayPtr{T}}) where {T} = T

## conversions

# to and from integers
## pointer to integer
Base.convert(::Type{T}, x::VEArrayPtr) where {T<:Integer} = T(UInt(x))
## integer to pointer
Base.convert(::Type{VEArrayPtr{T}}, x::Union{Int,UInt}) where {T} = VEArrayPtr{T}(x)
Int(x::VEArrayPtr)  = Base.bitcast(Int, x)
UInt(x::VEArrayPtr) = Base.bitcast(UInt, x)

# between regular and CUDA pointers
Base.convert(::Type{<:Ptr}, p::VEArrayPtr) =
    throw(ArgumentError("cannot convert a GPU array pointer to a CPU pointer"))

# between CUDA array pointers
Base.convert(::Type{VEArrayPtr{T}}, p::VEArrayPtr) where {T} = Base.bitcast(VEArrayPtr{T}, p)

# defer conversions to unsafe_convert
Base.cconvert(::Type{<:VEArrayPtr}, x) = x

# fallback for unsafe_convert
Base.unsafe_convert(::Type{P}, x::VEArrayPtr) where {P<:VEArrayPtr} = convert(P, x)

## limited pointer arithmetic & comparison

Base.isequal(x::VEArrayPtr, y::VEArrayPtr) = (x === y)
Base.isless(x::VEArrayPtr{T}, y::VEArrayPtr{T}) where {T} = x < y

Base.:(==)(x::VEArrayPtr, y::VEArrayPtr) = UInt(x) == UInt(y)
Base.:(<)(x::VEArrayPtr,  y::VEArrayPtr) = UInt(x) < UInt(y)
Base.:(-)(x::VEArrayPtr,  y::VEArrayPtr) = UInt(x) - UInt(y)

Base.:(+)(x::VEArrayPtr, y::Integer) = oftype(x, Base.add_ptr(UInt(x), (y % UInt) % UInt))
Base.:(-)(x::VEArrayPtr, y::Integer) = oftype(x, Base.sub_ptr(UInt(x), (y % UInt) % UInt))
Base.:(+)(x::Integer, y::VEArrayPtr) = y + x

#
## VERef object backed by an array at index i
#

struct VERefArray{T,A<:AbstractArray{T}} <: Ref{T}
    x::A
    i::Int
    VERefArray{T,A}(x,i) where {T,A<:AbstractArray{T}} = new(x,i)
end
VERefArray{T}(x::AbstractArray{T}, i::Int=1) where {T} = VERefArray{T,typeof(x)}(x, i)
VERefArray(x::AbstractArray{T}, i::Int=1) where {T} = VERefArray{T}(x, i)
Base.convert(::Type{VERef{T}}, x::AbstractArray{T}) where {T} = VERefArray(x, 1)

function Base.unsafe_convert(P::Type{VEPtr{T}}, b::VERefArray{T}) where T
    return pointer(b.x, b.i)
end
function Base.unsafe_convert(P::Type{VEPtr{Any}}, b::VERefArray{Any})
    return convert(P, pointer(b.x, b.i))
end
Base.unsafe_convert(::Type{VEPtr{Cvoid}}, b::VERefArray{T}) where {T} =
    convert(VEPtr{Cvoid}, Base.unsafe_convert(VEPtr{T}, b))

#
## Union with all VERef 'subtypes'
#

const VERefs{T} = Union{VEPtr{T}, VERefArray{T}}

## RefOrVERef

primitive type RefOrVERef{T} 64 end

Base.convert(::Type{RefOrVERef{T}}, x::Union{RefOrVERef{T}, Ref{T}, VERef{T}, VERefs{T}}) where {T} = x

# prefer conversion to CPU ref: this is generally cheaper
Base.convert(::Type{RefOrVERef{T}}, x) where {T} = Ref{T}(x)
Base.unsafe_convert(::Type{RefOrVERef{T}}, x::Ref{T}) where {T} =
    Base.bitcast(RefOrVERef{T}, Base.unsafe_convert(Ptr{T}, x))
Base.unsafe_convert(::Type{RefOrVERef{T}}, x) where {T} =
    Base.bitcast(RefOrVERef{T}, Base.unsafe_convert(Ptr{T}, x))

# support conversion from VE device ref
Base.unsafe_convert(::Type{RefOrVERef{T}}, x::VERefs{T}) where {T} =
    Base.bitcast(RefOrVERef{T}, Base.unsafe_convert(VEPtr{T}, x))

# support conversion from arrays
Base.convert(::Type{RefOrVERef{T}}, x::Array{T}) where {T} = convert(Ref{T}, x)
Base.convert(::Type{RefOrVERef{T}}, x::AbstractArray{T}) where {T} = convert(VERef{T}, x)
Base.unsafe_convert(P::Type{RefOrVERef{T}}, b::VERefArray{T}) where T =
    Base.bitcast(RefOrVERef{T}, Base.unsafe_convert(VERef{T}, b))


text:
@testset "Linear programs ($T)" for T in [Float32, Float64]

    using ProximalOperators
    using ProximalAlgorithms
    using LinearAlgebra

    # Solving LP with AFBA
    #
    #   minimize    c'x             -> f = <c,.>
    #   subject to  Ax = b          -> h = ind{b}, L = A
    #               x >= 0          -> g = ind{>=0}
    #
    # Dual LP
    #
    #   maximize    b'y
    #   subject to  A'y <= c
    #
    # Optimality conditions
    #
    #   x >= 0              [primal feasibility 1]
    #   Ax = b              [primal feasibility 2]
    #   A'y <= c            [dual feasibility]
    #   X'(c - A'y) = 0     [complementarity slackness]
    #

    function assert_lp_solution(c, A, b, x, y, tol)
        # Check and print solution quality measures (for some reason the
        # returned dual iterate is the negative of the dual LP variable y)

        nonneg = -minimum(min.(0.0, x))
        @test nonneg <= tol

        primal_feasibility = norm(A*x - b)
        @test primal_feasibility <= tol

        dual_feasibility = maximum(max.(0.0, -A'*y - c))
        @test dual_feasibility <= tol

        complementarity = abs(dot(c + A'*y, x))
        @test complementarity <= tol
    end

    n = 100 # primal dimension
    m = 80 # dual dimension (i.e. number of linear equalities)
    k = 50 # number of active dual constraints (must be 0 <= k <= n)
    x_star = vcat(rand(T, k), zeros(T, n-k)) # primal optimal point
    s_star = vcat(zeros(T, k), rand(T, n-k)) # dual optimal slack variable
    y_star = randn(T, m) # dual optimal point

    A = randn(T, m, n)
    b = A*x_star
    c = A'*y_star + s_star

    tol = 100*eps(T)
    maxit = 10_000

    @testset "AFBA" begin

        f = Linear(c)
        g = IndNonnegative()
        h = IndPoint(b)

        x0 = zeros(T, n)
        y0 = zeros(T, m)

        solver = ProximalAlgorithms.AFBA{T}(tol=tol, maxit=maxit)
        x, y, it = solver(x0, y0, f=f, g=g, h=h, L=A)

        @test eltype(x) == T
        @test eltype(y) == T

        @test it <= maxit

        assert_lp_solution(c, A, b, x, y, 1e2*tol)

    end

    @testset "VuCondat" begin

        f = Linear(c)
        g = IndNonnegative()
        h = IndPoint(b)

        x0 = zeros(T, n)
        y0 = zeros(T, m)

        solver = ProximalAlgorithms.VuCondat(T, tol=tol, maxit=maxit)
        x, y, it = solver(x0, y0, f=f, g=g, h=h, L=A)

        @test eltype(x) == T
        @test eltype(y) == T

        @test it <= maxit

        assert_lp_solution(c, A, b, x, y, 1e2*tol)

    end

    # TODO: add Chambolle-Pock (using separable sum of IndPoint and IndNonnegative)

    @testset "DavisYin" begin

        f = IndAffine(A, b)
        g = IndNonnegative()
        h = Linear(c)

        x0 = zeros(T, n)

        solver = ProximalAlgorithms.DavisYin{T}(gamma=T(1), tol=tol, maxit=maxit)
        xf, xg, it = solver(x0, f=f, g=g, h=h)

        @test eltype(xf) == T
        @test eltype(xg) == T

        @test it <= maxit

        @assert norm(xf - x_star) <= 1e2*tol
        @assert norm(xg - x_star) <= 1e2*tol

    end

end


text:
"""
This doesn't seem to be giving the right answer. Don't know why. The step function is either ill conditioned or doesn't kick the system like in the matlab simulation.
Solved. The heaviside function was ill conditioned. Made it a sigmoid and everything works
"""

function  NFKBModel(input)

@parameters t
D = Differential(t)

@parameters kprod kdeg k1 k2 k3 a1 a2 a3 t1 t2 c6a i1 kv c1 c2 c3 c4 c5 c4a c5a i1a e1a c1a c2a c3a e2a c1c c2c c3c
paramvars = [kprod, kdeg, k1, k2, k3, a1, a2, a3, t1, t2, c6a, i1, kv, c1, c2, c3, c4, c5, c4a, c5a, i1a, e1a, c1a, c2a, c3a, e2a, c1c, c2c, c3c]

@variables IKKN(t) IKKa(t) IKKi(t) IKKaIkBa(t) IKKaIkBaNfKb(t) NFkB(t) NFkBn(t) A20(t) A20t(t) IkBa(t) IkBan(t) IkBat(t) IkBaNfKb(t) IkBaNfKbn(t) Cgent(t)
statevars = [IKKN, IKKa, IKKi, IKKaIkBa, IKKaIkBaNfKb, NFkB, NFkBn, A20, A20t, IkBa, IkBan, IkBat, IkBaNfKb, IkBaNfKbn, Cgent]

eqs = [
D(IKKN) ~              kprod - kdeg*IKKN - k1*IKKN*input(t),
D(IKKa) ~                  -k3*IKKa-kdeg*IKKa-a2*IKKa*IkBa + t1*IKKaIkBa - a3*IKKa*IkBaNfKb + t2*IKKaIkBaNfKb + (k1*IKKN - k2*IKKa*A20)*input(t),
D(IKKi) ~                   k3*IKKa - kdeg*IKKi + k2*IKKa*A20*input(t),
D(IKKaIkBa) ~                   a2*IKKa*IkBa - t1*IKKaIkBa,
D(IKKaIkBaNfKb) ~                   a3*IKKa*IkBaNfKb - t2*IKKaIkBaNfKb,
D(NFkB) ~                   c6a*IkBaNfKb - a1*NFkB*IkBa + t2*IKKaIkBaNfKb - i1*NFkB,
D(NFkBn) ~                   i1*kv*NFkB - a1*IkBan*NFkBn,
D(A20) ~                   c4*A20t - c5*A20,
D(A20t) ~                   c2 + c1*NFkBn-c3*A20t,
D(IkBa) ~                  -a2*IKKa*IkBa - a1*IkBa*NFkB + c4a*IkBat - c5a*IkBa - i1a*IkBa + e1a*IkBan,
D(IkBan) ~                  -a1*IkBan*NFkBn + i1a*kv*IkBa-e1a*kv*IkBan,
D(IkBat) ~                  c2a + c1a*NFkBn - c3a*IkBat,
D(IkBaNfKb) ~                   a1*IkBa*NFkB - c6a*IkBaNfKb - a3*IKKa*IkBaNfKb + e2a*IkBaNfKbn,
D(IkBaNfKbn) ~                   a1*IkBan*NFkBn - e2a*kv*IkBaNfKbn,
D(Cgent) ~                   c2c + c1c*NFkBn - c3c*Cgent
]

ps = paramvars .=> [ 2.5e-5, 1.25e-4, 0.0025, 0.1, 0.0015, 0.5, 0.2, 1.0, 0.1, 0.1, 2.0e-5, 0.0025, 5.0, 5.0e-7, 0, 4.0e-4, 0.5, 3.0e-4, 0.5, 1.0e-4, 0.001, 5.0e-4, 5.0e-7, 0, 4.0e-4, 0.01, 5.0e-7, 0, 4.0e-4]

temp = zeros(15)
temp[13] = 0.06
tspan = (0., 50000.)
ic = statevars .=> temp

od = ODESystem(eqs, t, statevars, paramvars)
println("optional output map is [x[7], x[10] + x[13], x[9], x[1] + x[2] + x[3], x[2], x[12]]")
return od, ic, tspan, ps
end

function NFKB_output_map(x)
    return [x[7], x[10] + x[13], x[9], x[1] + x[2] + x[3], x[2], x[12]]
end
NFKB_output_map(x,t,integrator) = output_map(x)

# heaviside = soft_heaviside(0.01, 3600.)
# od0, ic0, tspan, ps0 = create(heaviside)
# prob = ODEProblem(od0, ic0, tspan, ps0)
# to_fix = ["c2c","c2","c2a","c3c", "c1c", "a2"]
# tstrct_fix = fix_params(last.(ps0), get_name_ids(ps0, to_fix))
# od2,ic2,ps2 = transform_problem(prob,tstrct_fix; unames = first.(ic0), pnames = first.(ps0))
# prob2 = ODEProblem(od2,ic2,tspan,ps2)
# tstrct_log = logabs_transform(last.(ps2))
# od3,ic3,ps3 = transform_problem(prob2,tstrct_log; unames = first.(ic2), pnames = first.(ps2))

# od,ic,ps = od3,ic3,ps3

# create() = (od,ic, tspan, ps)

# sol = solve(prob, Tsit5())
# plot(sol)

# saved_values = SavedValues(Float64, Array{Float64,1})
# output_callback = SavingCallback(output_map, saved_values)

# sol = solve(prob, AutoVern7(Rodas4()), callback = output_callback)

# c2,c2c,c2a,c3c,c1c,c5a


text:
abstract type AbstractSpatial end

include("line_and_plane.jl")
include("circle_and_sphere.jl")


text:
function (l::Conv1D)(li_1::Layer)
    l.prevLayer = li_1
    if ! in(l,li_1.nextLayers)
        push!(li_1.nextLayers, l)
    end
    padding = l.padding
    c = l.channels
    n_Hi, ci, m = l.inputS = li_1.outputS
    s_H = l.s
    f_H = l.f
    if l.padding == :same
        l.outputS = (n_Hi, c, m)
    elseif l.padding == :valid
        n_H = ((n_Hi - f_H) ÷ s_H) + 1
        l.outputS = (n_H, c, m)
    end
    return l
end

function (l::Conv2D)(li_1::Layer)
    l.prevLayer = li_1
    if ! in(l,li_1.nextLayers)
        push!(li_1.nextLayers, l)
    end
    padding = l.padding
    c = l.channels
    n_Hi, n_Wi, ci, m = l.inputS = li_1.outputS
    s_H, s_W = l.s
    f_H, f_W = l.f
    if l.padding == :same
        l.outputS = (n_Hi, n_Wi, c, m)
    elseif l.padding == :valid
        n_H = ((n_Hi - f_H) ÷ s_H) + 1
        n_W = ((n_Wi - f_W) ÷ s_W) + 1
        l.outputS = (n_H, n_W, c, m)
    end

    return l
end

function (l::Conv3D)(li_1::Layer)
    l.prevLayer = li_1
    if ! in(l,li_1.nextLayers)
        push!(li_1.nextLayers, l)
    end
    padding = l.padding
    c = l.channels
    n_Hi, n_Wi, n_Di, ci, m = l.inputS = li_1.outputS
    s_H, s_W, s_D = l.s
    f_H, f_W, f_D = l.f
    if l.padding == :same
        l.outputS = (n_Hi, n_Wi, n_Di, c, m)
    elseif l.padding == :valid
        n_H = ((n_Hi - f_H) ÷ s_H) + 1
        n_W = ((n_Wi - f_W) ÷ s_W) + 1
        n_D = ((n_Di - f_D) ÷ s_D) + 1
        l.outputS = (n_H, n_W, n_D, c, m)
    end
    return l
end

function (l::MaxPool1D)(li_1::Layer)
    l.prevLayer = li_1
    if ! in(l,li_1.nextLayers)
        push!(li_1.nextLayers, l)
    end
    padding = l.padding
    c = l.channels = li_1.channels
    n_Hi, ci, m = l.inputS = li_1.outputS
    s_H = l.s
    f_H = l.f
    if l.padding == :same
        l.outputS = (n_Hi, c, m)
    elseif l.padding == :valid
        n_H = ((n_Hi - f_H) ÷ s_H) + 1
        l.outputS = (n_H, c, m)
    end
    return l
end

function (l::MaxPool2D)(li_1::Layer)
    l.prevLayer = li_1
    if ! in(l,li_1.nextLayers)
        push!(li_1.nextLayers, l)
    end
    padding = l.padding
    c = l.channels = li_1.channels
    n_Hi, n_Wi, ci, m = l.inputS = li_1.outputS
    s_H, s_W = l.s
    f_H, f_W = l.f
    if l.padding == :same
        l.outputS = (n_Hi, n_Wi, c, m)
    elseif l.padding == :valid
        n_H = ((n_Hi - f_H) ÷ s_H) + 1
        n_W = ((n_Wi - f_W) ÷ s_W) + 1
        l.outputS = (n_H, n_W, c, m)
    end

    return l
end

function (l::MaxPool3D)(li_1::Layer)
    l.prevLayer = li_1
    if ! in(l,li_1.nextLayers)
        push!(li_1.nextLayers, l)
    end
    padding = l.padding
    c = l.channels = li_1.channels
    n_Hi, n_Wi, n_Di, ci, m = l.inputS = li_1.outputS
    s_H, s_W, s_D = l.s
    f_H, f_W, f_D = l.f
    if l.padding == :same
        l.outputS = (n_Hi, n_Wi, n_Di, c, m)
    elseif l.padding == :valid
        n_H = ((n_Hi - f_H) ÷ s_H) + 1
        n_W = ((n_Wi - f_W) ÷ s_W) + 1
        n_D = ((n_Di - f_D) ÷ s_D) + 1
        l.outputS = (n_H, n_W, n_D, c, m)
    end
    return l
end

function (l::AveragePool1D)(li_1::Layer)
    l.prevLayer = li_1
    if ! in(l,li_1.nextLayers)
        push!(li_1.nextLayers, l)
    end
    padding = l.padding
    c = l.channels = li_1.channels
    n_Hi, ci, m = l.inputS = li_1.outputS
    s_H = l.s
    f_H = l.f
    if l.padding == :same
        l.outputS = (n_Hi, c, m)
    elseif l.padding == :valid
        n_H = ((n_Hi - f_H) ÷ s_H) + 1
        l.outputS = (n_H, c, m)
    end
    return l
end

function (l::AveragePool2D)(li_1::Layer)
    l.prevLayer = li_1
    if ! in(l,li_1.nextLayers)
        push!(li_1.nextLayers, l)
    end
    padding = l.padding
    c = l.channels = li_1.channels
    n_Hi, n_Wi, ci, m = l.inputS = li_1.outputS
    s_H, s_W = l.s
    f_H, f_W = l.f
    if l.padding == :same
        l.outputS = (n_Hi, n_Wi, c, m)
    elseif l.padding == :valid
        n_H = ((n_Hi - f_H) ÷ s_H) + 1
        n_W = ((n_Wi - f_W) ÷ s_W) + 1
        l.outputS = (n_H, n_W, c, m)
    end

    return l
end

function (l::AveragePool3D)(li_1::Layer)
    l.prevLayer = li_1
    if ! in(l,li_1.nextLayers)
        push!(li_1.nextLayers, l)
    end
    padding = l.padding
    c = l.channels = li_1.channels
    n_Hi, n_Wi, n_Di, ci, m = l.inputS = li_1.outputS
    s_H, s_W, s_D = l.s
    f_H, f_W, f_D = l.f
    if l.padding == :same
        l.outputS = (n_Hi, n_Wi, n_Di, c, m)
    elseif l.padding == :valid
        n_H = ((n_Hi - f_H) ÷ s_H) + 1
        n_W = ((n_Wi - f_W) ÷ s_W) + 1
        n_D = ((n_Di - f_D) ÷ s_D) + 1
        l.outputS = (n_H, n_W, n_D, c, m)
    end
    return l
end


text:
import FFTW
using SparseArrays: sparse, spzeros
import ForwardDiff
import Zygote
using Statistics: mean
# import Plots
import CUDA: cu, allowscalar
allowscalar(false)
using Flux: gpu, cpu

using ProgressMeter
include("utils.jl")

γ = 0.5

function sumexp(x, γ)
    aP = exp.((x .- maximum(x)) / γ)
    aN = exp.(- (x .- minimum(x)) / γ)
    sum(x .* aP) / sum(aP) - sum(x .* aN) / sum(aN)
end

function We(net, xs, ys)
    # DEBUG using offx and offy
    x = [xs[n] - offx for (n, offx, offy) in net]
    y = [ys[n] - offy for (n, offx, offy) in net]
    wex = sumexp(x, γ)
    wey = sumexp(y, γ)
    return wex + wey
end

function W(nets, xs, ys)
    f = (e)->We(e, xs, ys)
    sum(f.(nets))
end

function We_grad_x_impl(x)
    aP = exp.((x .- maximum(x)) / γ)
    aN = exp.(- (x .- minimum(x)) / γ)
    bP = sum(aP)
    bN = sum(aN)
    cP = sum(x .* aP)
    cN = sum(x .* aN)
    res = (1 .+ x / γ) * bP .- 1/γ * cP
    res2 = (1 .- x / γ) * bN .+ 1 / γ * cN
    return res / (bP ^ 2) .* aP .- res2 / (bN .^ 2) .* aN
end

function We_grad_x(net, x, xory)
    # FIXME why sparse?
    # res = spzeros(length(x))
    res = zeros(length(x))
    # all the x values of net

    if xory == :x
        netx = [x[n] - offx for (n, offx, offy) in net]
    else
        netx = [x[n] - offy for (n, offx, offy) in net]
    end

    grad = We_grad_x_impl(netx)
    # map back to whole x
    for (i,(n,offx,offy)) in enumerate(net)
        res[n] = grad[i]
    end
    res
end

function W_grad_x(nets, x, xory)
    # batched by default
    # f = e->We_grad_x(e,x,xory)
    # res = spzeros(length(x))
    tmp = Array{Any}(nothing, length(nets))

    # export JULIA_NUM_THREADS=4
    #
    # Threads.@threads
    for i in 1:length(nets)
        # tmp[i] = f(nets[i])
        tmp[i] = We_grad_x(nets[i], x, xory)
    end
    sum(tmp)
end

function test()
    Threads.nthreads()
    @time W_grad_x(Es, xs, :x)
end

function rho_v_many(pxs, pys, pws, phs, bx1s, by1s, bx2s, by2s)
    x1 = max.(pxs .- pws / 2, bx1s)
    y1 = max.(pys .- phs / 2, by1s)
    x2 = min.(pxs .+ pws / 2, bx2s)
    y2 = min.(pys .+ phs / 2, by2s)
    # max.(x2 - x1, 0) .* max.(y2 - y1, 0)
    sum(max.(x2 - x1, 0) .* max.(y2 - y1, 0), dims=3)
end

function rho_cells(vs, R)
    xs, ys, ws, hs = vs

    # bx1s = reshape([R.xmin + R.bw * (i-1) for i in 1:R.M for j in 1:R.M], R.M, R.M)
    # bx2s = reshape([R.xmin + R.bw * i for i in 1:R.M for j in 1:R.M], R.M, R.M)
    # by1s = reshape([R.ymin + R.bh * (j-1) for i in 1:R.M for j in 1:R.M], R.M, R.M)
    # by2s = reshape([R.ymin + R.bh * j for i in 1:R.M for j in 1:R.M], R.M, R.M)

    pxs = reshape(xs, 1, 1, length(xs))
    pys = reshape(ys, 1, 1, length(ys))
    pws = reshape(ws, 1, 1, length(ws))
    phs = reshape(hs, 1, 1, length(hs))

    # FIXME moving data between cpu and gpu
    x1 = max.(pxs .- pws / 2, R.bx1s)
    y1 = max.(pys .- phs / 2, R.by1s)
    x2 = min.(pxs .+ pws / 2, R.bx2s)
    y2 = min.(pys .+ phs / 2, R.by2s)

    max.(x2 - x1, 0) .* max.(y2 - y1, 0)
end

function rho_fast(vs, R)
    xs, ys, ws, hs = vs
    # OK, I need to use batch
    pxs = reshape(xs, 1, 1, length(xs))
    pys = reshape(ys, 1, 1, length(ys))
    pws = reshape(ws, 1, 1, length(ws))
    phs = reshape(hs, 1, 1, length(hs))
    # FIXME whether to use GPU or not for small scale
    res = rho_v_many(pxs, pys, pws, phs,
                     R.bx1s, R.by1s, R.bx2s, R.by2s)
    dropdims(res, dims=3) |> cpu
end

function phi_b(rho, R)
    # auv = FFTW.dct(rho)
    auv = real.(FFTW.fft(rho))

    # R.M,R.M matrix
    phi = FFTW.idct(auv ./ R.wuv2)

    Ephix = imag.(FFTW.ifft(real.(FFTW.ifft(auv .* R.wv ./ R.wuv2, 1)), 2))
    Ephiy = real.(FFTW.ifft(imag.(FFTW.ifft(auv .* R.wu ./ R.wuv2, 1)), 2))

    phi, Ephix, Ephiy
end

function density(vs, R)
    # vs = xs, ys, ws, hs
    xs, ys, ws, hs = gpu.(vs)
    # 1. calculate rho
    @debug "calculating rho .."
    # rho = rho_all(vs, R)
    # rho = rho_fast(vs, R)
#     @time
    rho = rho_fast((xs, ys, ws, hs), R);
    # 2. calculate potential and field using FFT
    @debug "calculating potential and field .."
#     @time
    phi, Ephix, Ephiy = phi_b(rho, R);
    # calculate for each v

    # M,M,N
#     @time
    rrr = rho_cells((xs, ys, ws, hs), R);
    d = dropdims(sum(rrr .* gpu(phi), dims=(1,2)), dims=(1,2))
    dx = dropdims(sum(rrr .* gpu(Ephix), dims=(1,2)), dims=(1,2))
    dy = dropdims(sum(rrr .* gpu(Ephiy), dims=(1,2)), dims=(1,2))
    return d, dx, dy
end

struct Region
    # number of bins
    M
    # bin width and height
    bw
    bh
    # overall diearea
    xmin
    xmax
    ymin
    ymax
    # bounding box for bins (to compute only once)
    # TODO when printing out, omit these arrays
    bx1s
    by1s
    bx2s
    by2s
    # for FFT
    wu
    wv
    wuv2
end

function Region(xs, ys, ws, hs, diearea, M)
    xmin = 0
    ymin = 0
    xmax = diearea[1]
    ymax = diearea[2]

    # M = 10
    bw = (xmax-xmin) / M + 1e-8
    bh = (ymax-ymin) / M + 1e-8

    bx1s = reshape([xmin + bw * (i-1) for i in 1:M for j in 1:M], M, M)
    bx2s = reshape([xmin + bw * i for i in 1:M for j in 1:M], M, M)
    by1s = reshape([ymin + bh * (j-1) for i in 1:M for j in 1:M], M, M)
    by2s = reshape([ymin + bh * j for i in 1:M for j in 1:M], M, M)
    bx1s, by1s, bx2s, by2s = (bx1s, by1s, bx2s, by2s) .|> gpu

    # this is slow, but only compute once
    wuv_f(uv) = 2 * pi * uv / M
    wu = reshape([wuv_f(u) for u in 1:M for v in 1:M], M, M);
    wv = reshape([wuv_f(v) for u in 1:M for v in 1:M], M, M);
    wuv2 = reshape([wuv_f(u)^2 + wuv_f(v)^2 for u in 1:M for v in 1:M], M, M);

    Region(M, bw, bh, xmin, xmax, ymin, ymax,
           bx1s, by1s, bx2s, by2s,
           wu, wv, wuv2)
end

function hpwl(xs, ys, Es)
    res = map(Es) do e
        ns = [ei for (ei,offx,offy) in e]
        # max delta x
        abs(maximum((i->xs[i]).(ns)) - minimum((i->xs[i]).(ns))) +
            abs(maximum((i->ys[i]).(ns)) - minimum((i->ys[i]).(ns)))
    end
    sum(res)
end

function display_plot(p)
    path = tempname() * ".png"
    Plots.savefig(p, path)
    println("$(path)")
    println("#<Image: $(path)>")
end

function validate_region!(xs, ys, ws, hs, R)
    # FIXME consider rotation angles
    #
    # FIXME I belive this is pass-by-reference, so that I can modify the values
    # But strange enough, this doesn't work. I have to copy these inline in the
    # two places for now.
    #
    # UPDATE But the fixed locations should not change
    let idx = (xs .- ws ./ 2 .< R.xmin) .& (mask .== 1)
        xs[idx] .= (ws ./ 2 .+ R.xmin)[idx]
    end
    let idx = (xs .+ ws ./ 2 .> R.xmax) .& (mask .== 1)
        xs[idx] .= (ws ./ 2 .- R.xmax)[idx]
    end
    let idx = (ys .- hs ./ 2 .< R.ymin) .& (mask .== 1)
        ys[idx] .= (hs ./ 2 .+ R.ymin)[idx]
    end
    let idx = (ys .+ hs ./ 2 .> R.ymax) .& (mask .== 1)
        ys[idx] .= (hs ./ 2 .- R.ymax)[idx]
    end
end

# return a new pos
function place(xs, ys, ws, hs, Es, mask, diearea; vis=false, nsteps=50, nbins=300)
    xs = Float32.(xs)
    ys = Float32.(ys)
    ws = Float32.(ws)
    hs = Float32.(hs)

    # first, devide into bins
    # FIXME use a more formal way of deciding the bouding box
    # FIXME more bins
    R = Region(xs, ys, ws, hs, diearea, nbins)

    # loss: HPWL and density penalty
    hpwl(xs, ys, Es)
    mask

    # move all to the middle
    midx = (R.xmin + R.xmax) / 2
    midy = (R.ymin + R.ymax) / 2
    xs[mask.==1] .= midx
    ys[mask.==1] .= midy

    # iteratively solve the loss
    # 800 * 10 / 3600 = 2.2 hour
    # FIXME stop criteria: when the update is small enough for several epochs
    @showprogress 0.1 "Running for $nsteps steps .." for step in 1:nsteps
        @debug "step: $step"

        @debug "calculating W .."
        w = W(Es, xs, ys)
        @debug "calculating W gradient .."
        wgradx = W_grad_x(Es, xs, :x)
        wgrady = W_grad_x(Es, ys, :y)
        wgradx[mask .== 0] .= 0
        wgrady[mask .== 0] .= 0
        # remove the grad for fixed values

        # FIXME this density impl seems to be wrong
        @debug "Calculating fft density .."
        d, dx, dy = density((xs, ys, ws, hs), R) .|> cpu
        dx[mask .== 0] .= 0
        dy[mask .== 0] .= 0
        # weights HP here
        dx .*= 0.001
        dy .*= 0.001
        wgradx .*= 10
        wgrady .*= 10

        deltax = wgradx .- dx
        deltay = wgrady .- dy
        # do a cap
        # deltax[deltax .> 10] .= 10
        # deltax[deltax .< -10] .= -10
        # deltay[deltay .> 10] .= 10
        # deltay[deltay .< -10] .= -10
        loss = w + sum(d)

        # DEBUG use only wgrad
        # deltax = wgradx
        # deltay = wgrady
        # loss = w

        @debug("data", step, loss,
              mean(abs.(d)), mean(abs.(wgradx)),
              mean(abs.(dx)), hpwl(xs, ys, Es))

        # apply mask for fixed macros
        xs .-= deltax .* mask
        ys .-= deltay .* mask

        # map back to valid region
        # validate_region!(xs, ys, ws, hs, R)
        let idx = (xs .- ws ./ 2 .< R.xmin) .& (mask .== 1)
            xs[idx] .= (ws ./ 2 .+ R.xmin)[idx]
        end
        let idx = (xs .+ ws ./ 2 .> R.xmax) .& (mask .== 1)
            xs[idx] .= (ws ./ 2 .- R.xmax)[idx]
        end
        let idx = (ys .- hs ./ 2 .< R.ymin) .& (mask .== 1)
            ys[idx] .= (hs ./ 2 .+ R.ymin)[idx]
        end
        let idx = (ys .+ hs ./ 2 .> R.ymax) .& (mask .== 1)
            ys[idx] .= (hs ./ 2 .- R.ymax)[idx]
        end
        if vis display(visualize(xs, ys, ws, hs, R)) end
    end
    # if vis visualize_density((xs, ys, ws, hs), R) end
    xs, ys
end

function visualize_density(vs, R)
    xs, ys, ws, hs = vs
    visualize(xs, ys, ws, hs, R)
    rho = rho_fast(gpu.(vs), R)
    phi, Ephix, Ephiy = phi_b(rho, R)
    # display_plot(Plots.heatmap(rho))
    display_plot(Plots.heatmap(reverse(rho, dims=1)))
    # display_plot(Plots.heatmap(permutedims(rho)))
    display_plot(Plots.heatmap(reverse(Ephix, dims=1)))
    display_plot(Plots.heatmap(reverse(Ephiy, dims=1)))
    # display_plot(Plots.bar(sort(bxs)))
end

"""Calculate the number of overlap of the rectangle defined by (x,y,w,h) and all
the rectangles.

fasts but also deprecated because it is not easy to support angle.

"""
function cost_f_fast(xs, ys, ws, hs, x, y, w, h, R)
    x1 = x - w/2
    y1 = y - h/2
    x2 = x + w/2
    y2 = y + h/2

    x1s = xs .- ws ./ 2
    y1s = ys .- hs ./ 2

    x2s = xs .+ ws ./ 2
    y2s = ys .+ hs ./ 2

    return sum((x1s .< x2) .& (y1s .< y2) .& (x2s .> x1) .& (y2s .> y1))
end

function cost_f(xs, ys, as, ws, hs, x, y, a, w, h, R; except=[])
    # scale everything down
    #
    # UPDATE I don't need to scale it anymore due to use of Clipper
    l1, l2, l3, l4 = four_corners(x, y, a, w, h)

    # check for all the other rectangles
    ct = 0
    for i in 1:length(xs)
        if i in except continue end
        m1, m2, m3, m4 = four_corners(xs[i], ys[i], as[i], ws[i], hs[i])
        if isoverlap([l1, l2, l3, l4],
                     [m1, m2, m3, m4])
            ct += 1
        end
    end
    return ct
end

function test()
    isoverlap(four_corners(786, 617, 1.25, 122, 81),
              four_corners(707, 730, -0.314, 224, 81))
end

function four_corners(x, y, w, h)
    # compute the 4 corners of the rectangle
    l1 = (x - w/2, y - h/2)
    l2 = (x + w/2, y - h/2)
    l3 = (x + w/2, y + h/2)
    l4 = (x - w/2, y + h/2)
    return l1, l2, l3, l4
end

function four_corners(x, y, a, w, h)
    r = sqrt((h/2)^2 + (w/2)^2)
    return map([h/2, h/2, -h/2, -h/2],
               [w/2, -w/2, -w/2, w/2]) do Δh, Δw
                   # 1. compute theta
                   θ = atan(Δh, Δw)
                   x1 = x - cos(θ + a) * r
                   y1 = y + sin(θ + a) * r
                   (x1, y1)
               end
end

function test()
    four_corners(xs[1], ys[1], as[1], ws[1], hs[1])
    four_corners(xs[1], ys[1], ws[1], hs[1])
    four_corners(xs[1], ys[1], 1.2, ws[1], hs[1])
    map([1,2,3], [4,5,6]) do x, y
        (x,y)
    end
end

function accept(xs, ys, as, ws, hs, x, y, a, i, t, R)
    # calculate the xs[i], ys[i] conflicts with how many others
    c0 = cost_f(xs, ys, as, ws, hs,
                xs[i], ys[i], as[i], ws[i], hs[i],
                R, except=[i])
    # FIXME I should remove the conflict with itself
    c1 = cost_f(xs, ys, as, ws, hs,
                x, y, a, ws[i], hs[i],
                R, except=[i])
    if c1 < c0
        # @info "cost improves from $c0 to $c1"
        return true
    end
    p = exp(-(c1 - c0) / t)
    # @show p
    if rand() < p
        return true
    else
        return false
    end
end

function temperature(step)
    1 / log(step)
end

"""number of conflicts of all components
FIXME fixed position?
TODO report and debug this if there are still conflicts left?
"""
function compute_conflicts(xs, ys, as, ws, hs, mask, R)
    # conflicted items
    items = []
    for i in 1:length(xs)
        # TODO if the i-th item is fixed, then, skip it.
        if mask[i] == 0 continue end
        c = cost_f(xs, ys, as, ws, hs,
                   xs[i], ys[i], as[i], ws[i], hs[i],
                   R, except=[i])
        if c > 0
            append!(items, i)
        end
    end
    return items
end

function test()
    compute_conflicts(solxs2, solys2, solas, ws, hs, R)
end

function simulated_annealing_legalization(xs, ys, as, ws, hs, mask, diearea;
                                          vis=false, ncycles=20,
                                          nsteps=100, stepsize=50,
                                          theta_stepsize=0)
    # FIXME make a copy?
    xs = Float32.(xs)
    ys = Float32.(ys)
    ws = Float32.(ws)
    hs = Float32.(hs)
    as = Float32.(as)
    # FIXME this is for visulization and map back to valid region
    R = Region(xs, ys, ws, hs, diearea, 300)
    # actually I can probably parallelize the SA process, by trying to move
    # multiple components
    # t = 50
    # randomly choose a point
    @showprogress 0.1 "Running for $ncycles cycles .." for cycle in 2:ncycles+1
        t = temperature(cycle)
        for step in 1:nsteps
            # no movable components, not meaningful, for debugging purpose
            if isempty(findall(mask .== 1)) break end
            i = rand(findall(mask .== 1))
            if cost_f(xs, ys, as, ws, hs,
                      xs[i], ys[i], as[i], ws[i], hs[i],
                      R, except=[i]) == 0
                continue
            end
            # i = rand(1:length(xs))
            # FIXME the scale of the movement
            x = xs[i] + stepsize * randn()
            y = ys[i] + stepsize * randn()
            # FIXME maybe have a separate variable to control change angle or not?
            # FIXME what should be the step size? was 50 for x,y
            #
            # if theta_stepsize, this is effectively not changing angle
            a = as[i] + theta_stepsize * randn()

            # newxs = xs .+ randn()
            # newys = ys .+ randn()
            if accept(xs, ys, as, ws, hs,
                      x, y, a,
                      i, t, R)
                # @info "cycle $cycle step $step accepted"
                xs[i] = x
                ys[i] = y
                as[i] = a
            end
        end

        # map back to valid region
        # validate_region!(xs, ys, ws, hs, R)
        let idx = (xs .- ws ./ 2 .< R.xmin) .& (mask .== 1)
            xs[idx] .= (ws ./ 2 .+ R.xmin)[idx]
        end
        let idx = (xs .+ ws ./ 2 .> R.xmax) .& (mask .== 1)
            xs[idx] .= (ws ./ 2 .- R.xmax)[idx]
        end
        let idx = (ys .- hs ./ 2 .< R.ymin) .& (mask .== 1)
            ys[idx] .= (hs ./ 2 .+ R.ymin)[idx]
        end
        let idx = (ys .+ hs ./ 2 .> R.ymax) .& (mask .== 1)
            ys[idx] .= (hs ./ 2 .- R.ymax)[idx]
        end

        # print how many conflicts
        conflicts = compute_conflicts(xs, ys, as, ws, hs, mask, R)
        @info "cycle $cycle, remaining conflicts: $(length(conflicts))"

        # break if already no conflicts
        if length(conflicts) == 0 break end
        # FIXME visualize the angle
        if vis visualize(xs, ys, ws, hs, R) end
    end
    conflicts = compute_conflicts(xs, ys, as, ws, hs, mask, R)
    return xs, ys, as, conflicts
end

max_stars_repo_name: lihebi/bhdl

max_issues_repo_name: lihebi/bhdl

max_forks_repo_name: lihebi/bhdl


text:
using Flux: onehot, onehotbatch
using Random

const alphabet = [false, true]  # 0, 1

parity(x) = reduce(xor, x)

gendata(n::Int, k::Int) = gendata(n, k:k)
function gendata(n::Int, k::UnitRange{Int})
    X = bitrand.(rand(k, n))
    return [(onehotbatch(x, alphabet), onehot(y, alphabet)) for (x, y) in zip(X, parity.(X))]
end

max_stars_repo_name: c-p-murphy/model-zoo

max_issues_repo_name: c-p-murphy/model-zoo

max_forks_repo_name: c-p-murphy/model-zoo


text:
const TEST_NAME = "PCG_XSH_RS_128_64"

include("common.jl")

using RandomNumbers.PCG

r = PCGStateSetseq(UInt64, PCG_XSH_RS,
    (0x4e17a5abd5d47402ce332459f69eacfd, 0x96856028d0dc791c176537f21a77ab67))

test_all(r, 100_000_000)