text:
export add, add!, remove, remove!, posteriorParameters, logpostpred

"""
	add(d::ConjugatePostDistribution, X)

Add data to Posterior Distribution.
"""
function add(d::ConjugatePostDistribution, X)
	dist = deepcopy(d)
	add!(dist, X)
	return dist
end

"""
	add!(d::ConjugatePostDistribution, X)

Add data to Posterior Distribution (inplace).
"""
function add!(d::MultivariateConjugatePostDistribution, X::AbstractMatrix)
    @simd for i in 1:size(X, 1)
			@inbounds add!(d, @view X[i,:])
		end
end

function add!(d::UnivariateConjugatePostDistribution, X::AbstractMatrix)
	(N,D) = size(X)
	@assert D == 1
	add!(d, vec(X))
end

function add!(d::UnivariateConjugatePostDistribution, X::AbstractVector)
	@simd for i in 1:size(X, 1)
		@inbounds add!(d, X[i])
	end
end

function add!(d::WishartGaussian, X::AbstractVector)
	@assert length(X) == d.D

	d.n += 1
  d.sums += X
  d.ssums += X * X'
end

function add!(d::GaussianDiagonal, X::AbstractVector)
	for (dim, dist) in enumerate(d.dists)
		add!(dist, X[dim])
	end
end

function add!(d::ContinuousUnivariateConjugatePostDistribution, X::AbstractFloat)
	d.n += 1
	d.sums += X
	d.ssums += X^2
end

function add!(d::DirichletMultinomial, X)

    # process samples
    Dmin = minimum(X)
    Dmax = maximum(X)

    if Dmax > d.D
        throw(ArgumentError("Value of X and is larger than Multinomial Distribution!"))
    end

    if Dmin < 1
        throw(ArgumentError("Value of X and is smaller than 1!"))
    end

    d.n += length(X)

    for x in X
        d.counts[x] += 1
    end

    d.dirty = true
end

function add!(d::DirichletMultinomial, X::SparseVector{Int, Int})
		@assert length(X) == d.D

    d.n += 1
    d.counts += X
    d.dirty = true
end

function add!(d::BetaBernoulli, X::Bool)
	d.successes += X
	d.n += 1
end

"""
	remove(d::ConjugatePostDistribution, X::AbstractArray)

Remove data from Posterior Distribution.
"""
function remove(d::ConjugatePostDistribution, X)
		dist = deepcopy(d)
		remove!(dist, X)
		return dist
end

"""
	remove!(d::ConjugatePostDistribution, X)

Remove data from Posterior Distribution (inplace).
"""
function remove!(d::MultivariateConjugatePostDistribution, X::AbstractMatrix)
    @simd for i in 1:size(X, 1)
			@inbounds remove!(d, @view X[i,:])
		end
end

function remove!(d::UnivariateConjugatePostDistribution, X::AbstractMatrix)
	(N,D) = size(X)
	@assert D == 1
	remove!(d, vec(X))
end

function remove!(d::UnivariateConjugatePostDistribution, X::AbstractVector)
	@simd for i in 1:size(X, 1)
		@inbounds remove!(d, X[i])
	end
end

function remove!(d::WishartGaussian, X::AbstractVector)
		@assert length(X) == d.D

    d.n -= 1
    d.sums -= X
    d.ssums -= X * X'
end

function remove!(d::ContinuousUnivariateConjugatePostDistribution, X::AbstractFloat)
    d.n -= 1
    d.sums -= X
    d.ssums -= X^2
end

function remove!(d::GaussianDiagonal, X::AbstractVector)
	for (dim, dist) in enumerate(d.dists)
		remove!(dist, X[dim])
	end
end

function remove!(d::DirichletMultinomial, X)

    # process samples
    Dmin = minimum(X)
    Dmax = maximum(X)

    if Dmax > d.D
        throw(ArgumentError("Value of X and is larger than Multinomial Distribution!"))
    end

    if Dmin < 1
        throw(ArgumentError("Value of X and is smaller than 1!"))
    end

    d.n -= length(X)

    for x in X
        d.counts[x] -= 1
    end

    d.dirty = true
end

function remove!(d::DirichletMultinomial, X::SparseVector{Int, Int})

   # process samples
   @assert length(X) == d.D

   d.n -= 1
   d.counts -= X
   d.dirty = true
end

function remove!(d::BetaBernoulli, X::Bool)
		d.successes -= X
		d.n -= 1
end

"""
	isdistempty(d::ConjugatePostDistribution)

Check if distribution is empty (contains no samples).
"""
function isdistempty(d::ConjugatePostDistribution)
    return d.n <= 0
end

"""
	posteriorParameters(d::ConjugatePostDistribution)

Compute posterior distribution parameters.
"""
function posteriorParameters(d::WishartGaussian)
	if isdistempty(d)
	  return (d.μ0, d.κ0, d.ν0, d.Σ0)
	else

		# statistics
		sample_mu = d.sums / d.n

		# make sure values are not NaN
		sample_mu[sample_mu .!= sample_mu] .= 0

		# compute posterior parameters
		κ = d.κ0 + d.n
		ν = d.ν0 + d.n
		μ = (d.κ0 * d.μ0 + d.n * sample_mu) / κ
		Σ = d.Σ0 + d.ssums - κ * (μ * μ') + ( d.κ0 * (d.μ0 * d.μ0') )

		return (μ, κ, ν, Σ)
	end
end

function posteriorParameters(d::GammaNormal)
	if isdistempty(d)
	  return (d.μ0, d.λ0, d.α0, d.β0)
	else

		# statistics
    sample_mu = d.sums / d.n

    # make sure values are not NaN
    sample_mu = sample_mu != sample_mu ? 0 : sample_mu

    # compute posterior parameters
    μ = (d.λ0 * d.μ0 + d.sums) / (d.λ0 + d.n)
    λ = d.λ0 + d.n
    α = d.α0 + (d.n / 2)
    s = (d.ssums / d.n) - (sample_mu * sample_mu)
    β = d.β0 + 1/2 * (d.n * s + (d.λ0 * d.n * (sample_mu - d.μ0)^2 ) / (d.λ0 + d.n) )

		return (μ, λ, α, β)
	end
end

function posteriorParameters(d::NormalNormal)
	if isdistempty(d)
	  return (d.μ0, d.σ0)
	else

		# statistics
	  sample_mu = d.sums / d.n
		sample_var = (d.n * d.ssums - d.sums^2) / (d.n*(d.n-1))
		sample_var = (d.ssums - (d.sums^2)/d.n) / (d.n - 1)
		sample_var = 10.0

		# make sure values are not NaN
		sample_mu = sample_mu != sample_mu ? 0 : sample_mu

		σ = (sample_var * d.σ0) / (d.n * d.σ0 + sample_var)
	  μ = σ * (d.n*sample_mu/sample_var + d.μ0/d.σ0)

		return (μ, σ)
	end
end

function posteriorParameters(d::BetaBernoulli)
	α = d.α0 + d.successes
	β = d.β0 + d.n - d.successes
	return (α, β)
end

"""
	logpostpred(d::ConjugatePostDistribution, X)

Compute log posterior predictive.
"""
function logpostpred(d::WishartGaussian, x)

	(μ, κ, ν, Σ) = posteriorParameters(d)

  # posterior predictive of Normal Inverse Wishart is student-t Distribution, see:
  # K. Murphy, Conjugate Bayesian analysis of the Gaussian distribution. Eq. 258
	C = Σ * ((κ + 1) / (κ * (ν - d.D + 1)))
	@assert isposdef(C)

  return Distributions.logpdf(Distributions.MvTDist(ν - d.D + 1, μ, C), x)
end

function logpostpred(d::GammaNormal, x)
	return Float64[logpostpred(d, xi) for xi in x]
end

function logpostpred(d::GammaNormal, x::Number)

   (μ, λ, α, β) = posteriorParameters(d)

    # posterior predictive of Normal Gamma is student-t Distribution
    df = 2 * α
    mean = μ
    sigma = ( β * (λ + 1) ) / (λ * α)

		return tlogpdf(x, df, mean, sigma)
end

function logpostpred(d::GaussianDiagonal, x)
	return vec(sum([logpostpred(di, x[dim, :]) for (dim, di) in enumerate(d.dists)]))
end

function logpostpred(d::NormalNormal, x)
	return Float64[logpostpred(d, xi) for xi in x]
end

function logpostpred(d::NormalNormal, x::Number)

	(μ, σ) = posteriorParameters(d)

	return normlogpdf(μ, sqrt(σ), x)
end

"Log PMF for MultinomialDirichlet."
function logpostpred(d::DirichletMultinomial, x)

    N = length(x)

	  # TODO: This is bad and slow code, improve!!!
    # construct sparse vector
    xx = spzeros(d.D, N)
    for i in 1:N
      xx[x[i]] += 1
    end

		m = sum(xx)
		mi = nnz(xx)

		if d.dirty
			d.Z2 = lgamma(d.α0 + d.n)
			d.Z3 = lgamma( d.α0/d.D + d.counts )
			d.dirty = false
		end

	 l1 = lgamma(m + 1) - sum(lgamma(mi + 1))
	 l2 = d.Z2 - lgamma(d.α0 + d.n + m)
	 l3 = sum( lgamma( d.α0 / d.D + d.counts + xx ) - d.Z3 )

	 return [l1 + l2 + l3]
end

"Log PMF for MultinomialDirichlet."
function logpostpred(d::DirichletMultinomial, x::SparseVector{Int, Int})

   D = size(x, 1)
   N = size(x, 2)

   if N > 1
      throw(ErrorException("Multiple samples are not supported yet!"))
   end

   m = sum(x)
   mi = nnz(x)

   if d.dirty
      d.Z2 = lgamma(d.α0 + d.n)
      d.Z3 = lgamma( d.α0/d.D + d.counts )

      d.dirty = false
   end

	l1 = lgamma(m + 1) - sum(lgamma(mi + 1))
	l2 = d.Z2 - lgamma(d.α0 + d.n + m)
	l3 = sum( lgamma( d.α0 / d.D + d.counts + x ) - d.Z3 )

   return [l1 + l2 + l3]
end

"Log PMF for BernoulliBeta."
function logpostpred(d::BetaBernoulli, X::AbstractArray)
	return Float64[logpostpred(d, x) for x in X]
end

function logpostpred(d::BetaBernoulli, X::Bool)

  	# posterior
		(α, β) = posteriorParameters(d)

   return log(α) - log(α + β)
end


text:
module operators_lazysum

export LazySum

import Base: ==, *, /, +, -
import ..operators

using ..bases, ..states, ..operators, ..operators_dense

"""
    LazySum([factors,] operators)

Lazy evaluation of sums of operators.

All operators have to be given in respect to the same bases. The field
`factors` accounts for an additional multiplicative factor for each operator
stored in the `operators` field.
"""
mutable struct LazySum <: Operator
    basis_l::Basis
    basis_r::Basis
    factors::Vector{Complex128}
    operators::Vector{Operator}

    function LazySum(factors::Vector{Complex128}, operators::Vector{Operator})
        @assert length(operators)>0
        @assert length(operators)==length(factors)
        for i = 2:length(operators)
            @assert operators[1].basis_l == operators[i].basis_l
            @assert operators[1].basis_r == operators[i].basis_r
        end
        new(operators[1].basis_l, operators[1].basis_r, factors, operators)
    end
end
LazySum(factors::Vector{T}, operators::Vector) where {T<:Number} = LazySum(complex(factors), Operator[op for op in operators])
LazySum(operators::Operator...) = LazySum(ones(Complex128, length(operators)), Operator[operators...])

Base.copy(x::LazySum) = LazySum(copy(x.factors), [copy(op) for op in x.operators])

Base.full(op::LazySum) = sum(op.factors .* full.(op.operators))
Base.sparse(op::LazySum) = sum(op.factors .* sparse.(op.operators))

==(x::LazySum, y::LazySum) = (x.basis_l == y.basis_l) && (x.basis_r == y.basis_r) && x.operators==y.operators && x.factors==y.factors

# Arithmetic operations
+(a::LazySum, b::LazySum) = (check_samebases(a,b); LazySum([a.factors; b.factors], [a.operators; b.operators]))

-(a::LazySum) = LazySum(-a.factors, a.operators)
-(a::LazySum, b::LazySum) = (check_samebases(a,b); LazySum([a.factors; -b.factors], [a.operators; b.operators]))

*(a::LazySum, b::Number) = LazySum(b*a.factors, a.operators)
*(a::Number, b::LazySum) = LazySum(a*b.factors, b.operators)

/(a::LazySum, b::Number) = LazySum(a.factors/b, a.operators)

operators.dagger(op::LazySum) = LazySum(conj.(op.factors), dagger.(op.operators))

operators.trace(op::LazySum) = sum(op.factors .* trace.(op.operators))

function operators.ptrace(op::LazySum, indices::Vector{Int})
    operators.check_ptrace_arguments(op, indices)
    rank = length(op.basis_l.shape) - length(indices)
    D = Operator[ptrace(op_i, indices) for op_i in op.operators]
    LazySum(op.factors, D)
end

operators.normalize!(op::LazySum) = (op.factors /= trace(op))

operators.permutesystems(op::LazySum, perm::Vector{Int}) = LazySum(op.factors, Operator[permutesystems(op_i, perm) for op_i in op.operators])

operators.identityoperator(::Type{LazySum}, b1::Basis, b2::Basis) = LazySum(identityoperator(b1, b2))

# Fast in-place multiplication
function operators.gemv!(alpha, a::LazySum, b::Ket, beta, result::Ket)
    operators.gemv!(alpha*a.factors[1], a.operators[1], b, beta, result)
    for i=2:length(a.operators)
        operators.gemv!(alpha*a.factors[i], a.operators[i], b, 1, result)
    end
end

function operators.gemv!(alpha, a::Bra, b::LazySum, beta, result::Bra)
    operators.gemv!(alpha*b.factors[1], a, b.operators[1], beta, result)
    for i=2:length(b.operators)
        operators.gemv!(alpha*b.factors[i], a, b.operators[i], 1, result)
    end
end

end # module

ext: jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

#####################
# structs/constants #
#####################

# N.B.: Const/PartialStruct/InterConditional are defined in Core, to allow them to be used
# inside the global code cache.
#
# # The type of a value might be constant
# struct Const
#     val
# end
#
# struct PartialStruct
#     typ
#     fields::Vector{Any} # elements are other type lattice members
# end
import Core: Const, PartialStruct

# The type of this value might be Bool.
# However, to enable a limited amount of back-propagation,
# we also keep some information about how this Bool value was created.
# In particular, if you branch on this value, then may assume that in
# the true branch, the type of `var` will be limited by `vtype` and in
# the false branch, it will be limited by `elsetype`. Example:
# ```
# cond = isa(x::Union{Int, Float}, Int)::Conditional(x, Int, Float)
# if cond
#    # May assume x is `Int` now
# else
#    # May assume x is `Float` now
# end
# ```
struct Conditional
    var::SlotNumber
    vtype
    elsetype
    function Conditional(
                var,
                @nospecialize(vtype),
                @nospecialize(nottype))
        return new(var, vtype, nottype)
    end
end

# # Similar to `Conditional`, but conveys inter-procedural constraints imposed on call arguments.
# # This is separate from `Conditional` to catch logic errors: the lattice element name is InterConditional
# # while processing a call, then Conditional everywhere else. Thus InterConditional does not appear in
# # CompilerTypes—these type's usages are disjoint—though we define the lattice for InterConditional.
# struct InterConditional
#     slot::Int
#     vtype
#     elsetype
# end
import Core: InterConditional
const AnyConditional = Union{Conditional,InterConditional}

struct PartialTypeVar
    tv::TypeVar
    # N.B.: Currently unused, but would allow turning something back
    # into Const, if the bounds are pulled out of this TypeVar
    lb_certain::Bool
    ub_certain::Bool
    PartialTypeVar(tv::TypeVar, lb_certain::Bool, ub_certain::Bool) = new(tv, lb_certain, ub_certain)
end

# Wraps a type and represents that the value may also be undef at this point.
# (only used in optimize, not abstractinterpret)
# N.B. in the lattice, this is epsilon bigger than `typ` (even Any)
struct MaybeUndef
    typ
    MaybeUndef(@nospecialize(typ)) = new(typ)
end

struct StateUpdate
    var::SlotNumber
    vtype::VarState
    state::VarTable
    conditional::Bool
end

# Represent that the type estimate has been approximated, due to "causes"
# (only used in abstract interpretion, doesn't appear in optimization)
# N.B. in the lattice, this is epsilon smaller than `typ` (except Union{})
struct LimitedAccuracy
    typ
    causes::IdSet{InferenceState}
    function LimitedAccuracy(@nospecialize(typ), causes::IdSet{InferenceState})
        @assert !isa(typ, LimitedAccuracy) "malformed LimitedAccuracy"
        return new(typ, causes)
    end
end

"""
    struct NotFound end
    const NOT_FOUND = NotFound()

A special sigleton that represents a variable has not been analyzed yet.
Particularly, all SSA value types are initialized as `NOT_FOUND` when creating a new `InferenceState`.
Note that this is only used for `smerge`, which updates abstract state `VarTable`,
and thus we don't define the lattice for this.
"""
struct NotFound end

const NOT_FOUND = NotFound()

const CompilerTypes = Union{MaybeUndef, Const, Conditional, NotFound, PartialStruct}
==(x::CompilerTypes, y::CompilerTypes) = x === y
==(x::Type, y::CompilerTypes) = false
==(x::CompilerTypes, y::Type) = false

#################
# lattice logic #
#################

# `Conditional` and `InterConditional` are valid in opposite contexts
# (i.e. local inference and inter-procedural call), as such they will never be compared
function issubconditional(a::C, b::C) where {C<:AnyConditional}
    if is_same_conditionals(a, b)
        if a.vtype ⊑ b.vtype
            if a.elsetype ⊑ b.elsetype
                return true
            end
        end
    end
    return false
end

is_same_conditionals(a::Conditional,      b::Conditional)      = slot_id(a.var) === slot_id(b.var)
is_same_conditionals(a::InterConditional, b::InterConditional) = a.slot === b.slot

is_lattice_bool(@nospecialize(typ)) = typ !== Bottom && typ ⊑ Bool

maybe_extract_const_bool(c::Const) = (val = c.val; isa(val, Bool)) ? val : nothing
function maybe_extract_const_bool(c::AnyConditional)
    (c.vtype === Bottom && !(c.elsetype === Bottom)) && return false
    (c.elsetype === Bottom && !(c.vtype === Bottom)) && return true
    nothing
end
maybe_extract_const_bool(@nospecialize c) = nothing

"""
    a ⊑ b -> Bool

The non-strict partial order over the type inference lattice.
"""
@nospecialize(a) ⊑ @nospecialize(b) = begin
    if isa(b, LimitedAccuracy)
        if !isa(a, LimitedAccuracy)
            return false
        end
        if b.causes ⊈ a.causes
            return false
        end
        b = b.typ
    end
    isa(a, LimitedAccuracy) && (a = a.typ)
    if isa(a, MaybeUndef) && !isa(b, MaybeUndef)
        return false
    end
    isa(a, MaybeUndef) && (a = a.typ)
    isa(b, MaybeUndef) && (b = b.typ)
    b === Any && return true
    a === Any && return false
    a === Union{} && return true
    b === Union{} && return false
    @assert !isa(a, TypeVar) "invalid lattice item"
    @assert !isa(b, TypeVar) "invalid lattice item"
    if isa(a, AnyConditional)
        if isa(b, AnyConditional)
            return issubconditional(a, b)
        elseif isa(b, Const) && isa(b.val, Bool)
            return maybe_extract_const_bool(a) === b.val
        end
        a = Bool
    elseif isa(b, AnyConditional)
        return false
    end
    if isa(a, PartialStruct)
        if isa(b, PartialStruct)
            if !(length(a.fields) == length(b.fields) && a.typ <: b.typ)
                return false
            end
            for i in 1:length(b.fields)
                # XXX: let's handle varargs later
                ⊑(a.fields[i], b.fields[i]) || return false
            end
            return true
        end
        return isa(b, Type) && a.typ <: b
    elseif isa(b, PartialStruct)
        if isa(a, Const)
            nfields(a.val) == length(b.fields) || return false
            widenconst(b).name === widenconst(a).name || return false
            # We can skip the subtype check if b is a Tuple, since in that
            # case, the ⊑ of the elements is sufficient.
            if b.typ.name !== Tuple.name && !(widenconst(a) <: widenconst(b))
                return false
            end
            for i in 1:nfields(a.val)
                # XXX: let's handle varargs later
                isdefined(a.val, i) || continue # since ∀ T Union{} ⊑ T
                ⊑(Const(getfield(a.val, i)), b.fields[i]) || return false
            end
            return true
        end
        return false
    end
    if isa(a, PartialOpaque)
        if isa(b, PartialOpaque)
            (a.parent === b.parent && a.source === b.source) || return false
            return (widenconst(a) <: widenconst(b)) &&
                ⊑(a.env, b.env)
        end
        return widenconst(a) ⊑ b
    end
    if isa(a, Const)
        if isa(b, Const)
            return a.val === b.val
        end
        # TODO: `b` could potentially be a `PartialTypeVar` here, in which case we might be
        # able to return `true` in more cases; in the meantime, just returning this is the
        # most conservative option.
        return isa(b, Type) && isa(a.val, b)
    elseif isa(b, Const)
        if isa(a, DataType) && isdefined(a, :instance)
            return a.instance === b.val
        end
        return false
    elseif isa(a, PartialTypeVar) && b === TypeVar
        return true
    elseif isa(a, Type) && isa(b, Type)
        return a <: b
    else # handle this conservatively in the remaining cases
        return a === b
    end
end

"""
    a ⊏ b -> Bool

The strict partial order over the type inference lattice.
This is defined as the irreflexive kernel of `⊑`.
"""
@nospecialize(a) ⊏ @nospecialize(b) = a ⊑ b && !⊑(b, a)

"""
    a ⋤ b -> Bool

This order could be used as a slightly more efficient version of the strict order `⊏`,
where we can safely assume `a ⊑ b` holds.
"""
@nospecialize(a) ⋤ @nospecialize(b) = !⊑(b, a)

# Check if two lattice elements are partial order equivalent. This is basically
# `a ⊑ b && b ⊑ a` but with extra performance optimizations.
function is_lattice_equal(@nospecialize(a), @nospecialize(b))
    a === b && return true
    if isa(a, PartialStruct)
        isa(b, PartialStruct) || return false
        length(a.fields) == length(b.fields) || return false
        widenconst(a) == widenconst(b) || return false
        for i in 1:length(a.fields)
            is_lattice_equal(a.fields[i], b.fields[i]) || return false
        end
        return true
    end
    isa(b, PartialStruct) && return false
    if a isa Const
        if issingletontype(b)
            return a.val === b.instance
        end
        return false
    end
    if b isa Const
        if issingletontype(a)
            return a.instance === b.val
        end
        return false
    end
    if isa(a, PartialOpaque)
        isa(b, PartialOpaque) || return false
        widenconst(a) == widenconst(b) || return false
        a.source === b.source || return false
        a.parent === b.parent || return false
        return is_lattice_equal(a.env, b.env)
    end
    return a ⊑ b && b ⊑ a
end

# compute typeintersect over the extended inference lattice,
# as precisely as we can,
# where v is in the extended lattice, and t is a Type.
function tmeet(@nospecialize(v), @nospecialize(t))
    if isa(v, Const)
        if !has_free_typevars(t) && !isa(v.val, t)
            return Bottom
        end
        return v
    elseif isa(v, PartialStruct)
        has_free_typevars(t) && return v
        widev = widenconst(v)
        if widev <: t
            return v
        end
        ti = typeintersect(widev, t)
        valid_as_lattice(ti) || return Bottom
        @assert widev <: Tuple
        new_fields = Vector{Any}(undef, length(v.fields))
        for i = 1:length(new_fields)
            vfi = v.fields[i]
            if isvarargtype(vfi)
                new_fields[i] = vfi
            else
                new_fields[i] = tmeet(vfi, widenconst(getfield_tfunc(t, Const(i))))
                if new_fields[i] === Bottom
                    return Bottom
                end
            end
        end
        return tuple_tfunc(new_fields)
    elseif isa(v, Conditional)
        if !(Bool <: t)
            return Bottom
        end
        return v
    end
    ti = typeintersect(widenconst(v), t)
    valid_as_lattice(ti) || return Bottom
    return ti
end

widenconst(c::AnyConditional) = Bool
widenconst((; val)::Const) = isa(val, Type) ? Type{val} : typeof(val)
widenconst(m::MaybeUndef) = widenconst(m.typ)
widenconst(c::PartialTypeVar) = TypeVar
widenconst(t::PartialStruct) = t.typ
widenconst(t::PartialOpaque) = t.typ
widenconst(t::Type) = t
widenconst(t::TypeVar) = error("unhandled TypeVar")
widenconst(t::TypeofVararg) = error("unhandled Vararg")
widenconst(t::LimitedAccuracy) = error("unhandled LimitedAccuracy")

issubstate(a::VarState, b::VarState) = (a.typ ⊑ b.typ && a.undef <= b.undef)

function smerge(sa::Union{NotFound,VarState}, sb::Union{NotFound,VarState})
    sa === sb && return sa
    sa === NOT_FOUND && return sb
    sb === NOT_FOUND && return sa
    issubstate(sa, sb) && return sb
    issubstate(sb, sa) && return sa
    return VarState(tmerge(sa.typ, sb.typ), sa.undef | sb.undef)
end

@inline tchanged(@nospecialize(n), @nospecialize(o)) = o === NOT_FOUND || (n !== NOT_FOUND && !(n ⊑ o))
@inline schanged(@nospecialize(n), @nospecialize(o)) = (n !== o) && (o === NOT_FOUND || (n !== NOT_FOUND && !issubstate(n::VarState, o::VarState)))

function widenconditional(@nospecialize typ)
    if isa(typ, AnyConditional)
        if typ.vtype === Union{}
            return Const(false)
        elseif typ.elsetype === Union{}
            return Const(true)
        else
            return Bool
        end
    end
    return typ
end
widenconditional(t::LimitedAccuracy) = error("unhandled LimitedAccuracy")

widenwrappedconditional(@nospecialize(typ))   = widenconditional(typ)
widenwrappedconditional(typ::LimitedAccuracy) = LimitedAccuracy(widenconditional(typ.typ), typ.causes)

ignorelimited(@nospecialize typ) = typ
ignorelimited(typ::LimitedAccuracy) = typ.typ

function stupdate!(state::Nothing, changes::StateUpdate)
    newst = copy(changes.state)
    changeid = slot_id(changes.var)
    newst[changeid] = changes.vtype
    # remove any Conditional for this slot from the vtable
    # (unless this change is came from the conditional)
    if !changes.conditional
        for i = 1:length(newst)
            newtype = newst[i]
            if isa(newtype, VarState)
                newtypetyp = ignorelimited(newtype.typ)
                if isa(newtypetyp, Conditional) && slot_id(newtypetyp.var) == changeid
                    newtypetyp = widenwrappedconditional(newtype.typ)
                    newst[i] = VarState(newtypetyp, newtype.undef)
                end
            end
        end
    end
    return newst
end

function stupdate!(state::VarTable, changes::StateUpdate)
    newstate = nothing
    changeid = slot_id(changes.var)
    for i = 1:length(state)
        if i == changeid
            newtype = changes.vtype
        else
            newtype = changes.state[i]
        end
        oldtype = state[i]
        # remove any Conditional for this slot from the vtable
        # (unless this change is came from the conditional)
        if !changes.conditional && isa(newtype, VarState)
            newtypetyp = ignorelimited(newtype.typ)
            if isa(newtypetyp, Conditional) && slot_id(newtypetyp.var) == changeid
                newtypetyp = widenwrappedconditional(newtype.typ)
                newtype = VarState(newtypetyp, newtype.undef)
            end
        end
        if schanged(newtype, oldtype)
            newstate = state
            state[i] = smerge(oldtype, newtype)
        end
    end
    return newstate
end

function stupdate!(state::VarTable, changes::VarTable)
    newstate = nothing
    for i = 1:length(state)
        newtype = changes[i]
        oldtype = state[i]
        if schanged(newtype, oldtype)
            newstate = state
            state[i] = smerge(oldtype, newtype)
        end
    end
    return newstate
end

stupdate!(state::Nothing, changes::VarTable) = copy(changes)

stupdate!(state::Nothing, changes::Nothing) = nothing

function stupdate1!(state::VarTable, change::StateUpdate)
    changeid = slot_id(change.var)
    # remove any Conditional for this slot from the catch block vtable
    # (unless this change is came from the conditional)
    if !change.conditional
        for i = 1:length(state)
            oldtype = state[i]
            if isa(oldtype, VarState)
                oldtypetyp = ignorelimited(oldtype.typ)
                if isa(oldtypetyp, Conditional) && slot_id(oldtypetyp.var) == changeid
                    oldtypetyp = widenconditional(oldtypetyp)
                    if oldtype.typ isa LimitedAccuracy
                        oldtypetyp = LimitedAccuracy(oldtypetyp, (oldtype.typ::LimitedAccuracy).causes)
                    end
                    state[i] = VarState(oldtypetyp, oldtype.undef)
                end
            end
        end
    end
    # and update the type of it
    newtype = change.vtype
    oldtype = state[changeid]
    if schanged(newtype, oldtype)
        state[changeid] = smerge(oldtype, newtype)
        return true
    end
    return false
end

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
using SimplePaint
using Test

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


text:
# Basically taken from Flux. TODO: just use Tuple?
struct Composition{T<:Tuple}
  funcs::T
  Composition(xs...) = new{typeof(xs)}(xs)
end
Composition(x::Nothing) = nothing
Composition(c::Composition) = c
Base.getindex(x::Composition) = getindex(x.funcs)
Base.getindex(c::Composition, i::AbstractArray) = Composition(c.funcs[i]...)
Base.length(x::Composition) = length(x.funcs)
Base.first(x::Composition) = first(x.funcs)
Base.last(x::Composition) = last(x.funcs)
Base.tail(x::Composition) = Base.tail(x.funcs)
Base.iterate(x::Composition) = iterate(x.funcs)
Base.lastindex(x::Composition) = lastindex(x.funcs)

extend(c::Composition, el) = Composition(c.funcs..., el)
extend(el, c::Composition) = Composition(el, c.funcs...)
extend(a::Nothing, b::Nothing) = nothing
extend(c1::Composition, c2::Composition) = Composition(c1.funcs..., c2.funcs...)
extend(c::Composition, el::Nothing) = c
extend(el::Nothing, c::Composition) = c
extend(f, el::Nothing) = Composition(f)
extend(el::Nothing, f) = Composition(f)
extend(f1, f2) = Composition(f1, f2)

function Base.show(io::IO, c::Composition)
    print(io, "Composition(")
    join(io, c.funcs, ", ")
    print(io, ")")
end


text:
# MadNLP.jl.
# Created by Sungho Shin (sungho.shin@wisc.edu)

mutable struct MonolevelPartition
    g::Graph
    nparts::Int
    part::Vector{Int}
end

mutable struct MonolevelStruc
    V::Vector{Int}
    new_nbr::Vector{Int}
end

mutable struct BilevelPartition
    g_lower::Graph
    nparts_lower::Int
    part_lower::Vector{Int}
    V_lower::Vector{Vector{Int}}

    g_upper::Graph
    nparts_upper::Int
    part_upper::Vector{Int}
end

mutable struct BilevelStruc
    V_lower::Vector{Int}
    V_upper::Vector{Int}
    new_nbr_upper::Vector{Int}
end

mutable struct TwoStagePartition
    nparts::Int
    part::Vector{Int}
end

get_current_V(mls::MonolevelStruc) = mls.V
get_current_size(mls::MonolevelStruc) = length(mls.V)
get_full_size(mlp::MonolevelPartition) = nv(mlp.g)

function MonolevelPartition(csc::SparseMatrixCSC,part,nparts;max_size=0.)
    g   = Graph(csc)
    isempty(part) && (part=partition(g,nparts,alg=:KWAY))
    return MonolevelPartition(g,nparts,part)
end

function MonolevelStruc(mlp::MonolevelPartition,k;max_size=0.)
    V = findall(mlp.part.==k)
    new_nbr = expand!(V,mlp.g,max_size)
    return MonolevelStruc(V,new_nbr)
end

function expand!(mls::MonolevelStruc,mlp::MonolevelPartition,max_size)
    mls.new_nbr = expand!(mls.V,mlp.g,max_size,new_nbr=mls.new_nbr)
    return
end

get_current_V(bls::BilevelStruc) = bls.V_lower
get_current_size(bls::BilevelStruc) = length(bls.V_upper)
get_full_size(blp::BilevelPartition) = blp.nparts_lower

function BilevelPartition(csc,part_lower,nparts_lower,part_upper,nparts_upper;max_size=0.)
    g_lower   = Graph(csc)
    isempty(part_lower) && (part_lower= partition(g_lower,nparts_lower,alg=:KWAY))
    V_lower = Vector{Vector{Int}}(undef,nparts_lower)
    @blas_safe_threads for k=1:nparts_lower
        V_lower[k] = findall(part_lower.==k)
    end

    g_upper = Graph(nparts_lower)
    for e in edges(g_lower)
        add_edge!(g_upper,part_lower[src(e)],part_lower[dst(e)])
    end
    isempty(part_upper) && (part_upper = partition(g_upper,nparts_upper,alg=:KWAY))

    return BilevelPartition(g_lower,nparts_lower,part_lower,V_lower,g_upper,nparts_upper,part_upper)
end

function BilevelStruc(blp::BilevelPartition,k;max_size=0.)
    V_upper = findall(blp.part_upper.==k)
    new_nbr_upper = expand!(V_upper,blp.g_upper,max_size)
    V_lower = vcat(blp.V_lower[V_upper]...)

    return BilevelStruc(V_lower,V_upper,new_nbr_upper)
end

function TwoStagePartition(csc::SparseMatrixCSC,part,nparts)
    if isempty(part) || findfirst(x->x==0.,part) == nothing
        g = Graph(csc)
        isempty(part) && (part = partition(g,nparts,alg=:KWAY))
        mark_boundary!(g,part)
    end
    return TwoStagePartition(nparts,part)
end

Graph(csc::SparseMatrixCSC) = Graph(getelistcsc(csc.colptr,csc.rowval))
getelistcsc(colptr,rowval) = [Edge(i,Int(j)) for i=1:length(colptr)-1 for j in @view rowval[colptr[i]:colptr[i+1]-1]]

function expand!(bls::BilevelStruc,blp::BilevelPartition,max_size)
    orig_size = length(bls.V_upper)
    bls.new_nbr_upper = expand!(bls.V_upper,blp.g_upper,max_size,new_nbr=bls.new_nbr_upper)
    bls.V_lower = vcat(blp.V_lower[bls.V_upper]...)
    return
end

function expand!(V_om,g::Graph,max_size;
                 new_nbr=[])
    if isempty(new_nbr)
        new_nbr = Int[]
        for v in V_om
            append!(new_nbr,neighbors(g,v))
        end
        unique!(new_nbr)
        setdiff!(new_nbr,V_om)
    end

    old_nbr = V_om

    while (length(V_om) + length(new_nbr) < max_size) && length(V_om) < nv(g) && !isempty(new_nbr)
        append!(V_om,new_nbr)
        old_old_nbr = old_nbr
        old_nbr=new_nbr
        new_nbr = Int[]
        for v in old_nbr
            append!(new_nbr,neighbors(g,v))
        end
        unique!(new_nbr)
        setdiff!(new_nbr,old_old_nbr)
        setdiff!(new_nbr,old_nbr)
    end

    return new_nbr
end

function mark_boundary!(g,part)
    for e in edges(g)
        (part[src(e)]!=part[dst(e)] && part[src(e)]!= 0 && part[dst(e)] != 0) &&
            (part[src(e)] = 0; part[dst(e)] = 0)
    end
end


text:
struct MTable
    symels
    table
end

struct Group
    idxs
    elements
    mtable
    classes
    order
    ctable
end

function multifly(symels, A::Symel, B::Symel)
    Crrep = A.rrep * B.rrep
    for (i,g) in enumerate(symels)
        s = sum(abs.(g.rrep - Crrep))
        if abs(s) < tol
            return i,g
        end
    end
    println(i,g)
    throw(ArgumentError("No match found!"))
end

function build_mult_table(symels)
    h = length(symels)
    t = zeros(Int64,h,h)
    for (i,a) in enumerate(symels)
        for (j,b) in enumerate(symels)
            t[i,j] = multifly(symels, a, b)[1]
        end
    end
    return MTable(symels, t)
end

function find_inv(mtable, g)
    for i = 1:length(mtable.symels)
        if mtable.table[i,g] == 1
            return i
        end
    end
end

function find_classes(mtable)
    h = length(mtable.symels)
    assigned_elements = []
    classes = []
    for a = 1:h
        if a ∉ assigned_elements
            results = []
            for b = 1:h
                if a != b
                    binv = find_inv(mtable, b)
                    r1 = mtable.table[binv,a]
                    r = mtable.table[r1,b]
                    push!(results, r)
                end
            end
            push!(classes,unique(results))
            for i in unique(results)
                push!(assigned_elements, i)
            end
        end
    end
    return classes
end

function crst(mtable, R, S, T)
    t = T[1] # Choice of t doesn't matter
    counter = 0
    for r in R
        for s in S
            if mtable.table[r,s] == t
                counter += 1
            end
        end
    end
    return counter
end

function build_Mk(mtable, classes)
    N = length(classes)
    M = []
    for k = 1:N
        Mk = zeros(Float64,N,N)
        for i = 1:N
            for j = 1:N
                Mk[i,j] = crst(mtable, classes[k], classes[i], classes[j])
            end
        end
        push!(M, Hermitian(Mk))
    end
    return ℍVector(M)
end

function dixons_alg(classes, M, h)
    class_orders = []
    for c in classes
        push!(class_orders, length(c))
    end
    g = (2(h)^0.5)
    e = lcm(class_orders...)
    if e == 1
        throw(ArgumentError("I can't handle this yet"))
    end
    primes = get_primes(100) # This will probably need to be adjusted in case the lcm is large
    p = 0
    for prm in primes
        if (prm-1) % e == 0
            if prm > g
                p = prm
                break
            end
        end
    end
    λ = 1:p-1
    println("prime:$p")
    c = 1
    jimbo = nothing
    flag = false
    V = Array{Int64}(undef, length(M), 0)
    println("Doin stuff")
    for k = 1:length(M)
        for l in λ
            pizza = (M[k]-l*I) .% p
            nl = nullspace(pizza)
            sz = size(nl)[2]
            if sz == 0
                continue
            elseif sz == 1
                V = hcat(V, mynorm(nl))
            else
                #jimbo = mapslices(mynorm, nl, dims=(1))
                for v in 1:sz
                    mart = mynorm(nl[:,v])
                    V = hcat(V, mart)
                end
            end
        end
    end
    for v = 1:3
        A = (M[3]-3*I) .% p
        println(A*V[:,v])
    end
end

function mynorm(v)
    b = 0
    for i in v
        if i != 0
            b = i
            break
        end
    end
    return round.(v / b)
end

function get_primes(n)
    primes = []
    for i = 2:n
        if is_prime(i)
            push!(primes, i)
        end
    end
    return primes
end

function is_prime(n)
    if n == 2
        return true
    elseif n == 3
        return true
    end
    sn = floor((n)^0.5)
    for i = 2:sn
        if n % i == 0
            return false
        end
    end
    return true
end


text:
using DoubleDouble

n = 1_000_000
X = 64*rand(n)
Y = 64*rand(n)

Zbig = Array(Float64,n)
Zdouble = Array(Float64,n)

function bigrem!(X,Y,Z)
    for i = 1:length(Z)
        Z[i] = float64(rem(big(X[i])*big(Y[i]),1.0))
    end
end

function doublerem!(X,Y,Z)
    for i = 1:length(Z)
        Z[i] = float64(rem(Single(X[i])*Single(Y[i]),1.0))
    end
end

bigrem!(X,Y,Zbig)
doublerem!(X,Y,Zdouble)

X = 64*rand(n)
Y = 64*rand(n)

@time bigrem!(X,Y,Zbig)
@time doublerem!(X,Y,Zdouble)
@assert Zbig == Zdouble


text:
"""
    Abstract quantization type. Its subtypes are two singletons
corresponding to additive/residual and orthogonal quantization.
"""
abstract type AbstractQuantization end

struct OrthogonalQuantization <: AbstractQuantization end

struct AdditiveQuantization <: AbstractQuantization end

"""
    ArrayQuantizer{Q,U,D,T,N}

The array quantizer object. It transforms an array with elements
of type `T` into a 'quantized' version with elemetns of type `U`.

# Fields
  * `quantization::Q` type of quantization employed i.e. additive, orthogonal
  * `dims::NTuple{N,Int}` the original array dimensionality
  * `codebooks::Vector{CodeBook{U,T}}` the codebooks
  * `k::Int` the number of vector prototypes in each codebooks
  * `distance::D` the distance employed
  * `rot::Matrix{T}` rotation matrix
"""
struct ArrayQuantizer{Q,U,D,T,N}
    quantization::Q                   # type of quantization
    dims::NTuple{N, Int}              # original array size
    codebooks::Vector{CodeBook{U,T}}  # codebooks
    k::Int                            # number of codes/quantizer
    distance::D
    rot::Matrix{T}
end

# Aliases
const OrthogonalQuantizer{U,D,T,N} = ArrayQuantizer{OrthogonalQuantization,U,D,T,N}
const AdditiveQuantizer{U,D,T,N} = ArrayQuantizer{AdditiveQuantization,U,D,T,N}

# show methods
Base.show(io::IO, aq::ArrayQuantizer{Q,U,D,T,N}) where {Q,U,D,T,N} = begin
    m = length(codebooks(aq))
    quantstr = ifelse(Q <: OrthogonalQuantization, "OrthogonalQuantizer",
                      "AdditiveQuantizer")
    qstr = ifelse(m==1, "quantizer", "quantizers")
    cstr = ifelse(aq.k==1, "code", "codes")
    print(io, "$quantstr{$U,$D,$T,$N}, $m $qstr, $(aq.k) $cstr")
end

# Constructors
ArrayQuantizer(aa::AbstractMatrix{T};
               k::Int=DEFAULT_K,
               m::Int=DEFAULT_M,
               method::Symbol=DEFAULT_METHOD,
               distance::Distances.PreMetric=DEFAULT_DISTANCE,
               kwargs...) where {T} = begin
    U = quantized_eltype(k)
    quant = quantization_type(method)
    cbooks, rot = build_codebooks(aa, k, m, U, method=method,
                                  distance=distance; kwargs...)
    return ArrayQuantizer(quant, size(aa), cbooks, k, distance, rot)
end

ArrayQuantizer(aa::AbstractVector{T};
               k::Int=DEFAULT_K,
               m::Int=DEFAULT_M,
               method::Symbol=DEFAULT_METHOD,
               distance::Distances.PreMetric=DEFAULT_DISTANCE,
               kwargs...) where {T} = begin
    aq = ArrayQuantizer(aa', k=k, m=m, method=method, distance=distance; kwargs...)
    return ArrayQuantizer(aq.quantization, size(aa), codebooks(aq), k, distance, aq.rot)
end

"""
    build_quantizer(aa [;kwargs])

Builds an array quantizer using the input array `aa`.

# Keyword arguments
  * `k::Int` the number of vector prototypes in each codebook
  * `m::Int` the number of codebooks
  * `method::Symbol` the algorithm to be employed for codebook
generation; possible values are `:sample` (default), `:pq` for
classical k-means clustering codebooks and `:opq` for 'cartesian'
k-means clustering codebooks
  * `distance::PreMetric` the distance to be used in the
codebook generation methods and data encoding

Other codebook generation algorithm specific keyword arguments
such as `maxiter::Int` can be specified as well.
"""
build_quantizer(aa::AbstractArray; kwargs...) where{T,N} = ArrayQuantizer(aa; kwargs...)

"""
    codebooks(aq)

Access the codebooks field of the array quantizer `aq`.
"""
codebooks(aq::ArrayQuantizer) = aq.codebooks

"""
    quantize_data(aq, aa)

Returns a quantized version of the array `aa` using the array quantizer `aq`.
"""
function quantize_data(aq::OrthogonalQuantizer{U,D,T,2}, aa::AbstractMatrix{T}) where {U,D,T}
    nrows, ncols = size(aa)
    @assert nrows == aq.dims[1] "Quantized matrix needs to have $nrows rows"
    cbooks = codebooks(aq)
    m = length(cbooks)
    qa = Matrix{U}(undef, m, ncols)
    _aa = aq.rot * aa
    @inbounds @simd for i in 1:m
        rr = rowrange(nrows, m, i)
        qa[i,:] = encode(cbooks[i], _aa[rr,:], distance=aq.distance)
    end
    return qa
end

function quantize_data(aq::AdditiveQuantizer{U,D,T,2}, aa::AbstractMatrix{T}) where {U,D,T}
    nrows, ncols = size(aa)
    @assert nrows == aq.dims[1] "Quantized matrix needs to have $nrows rows"
    cbooks = codebooks(aq)
    m = length(cbooks)
    qa = Matrix{U}(undef, m, ncols)
    aac = copy(aa)
    @inbounds @simd for i in 1:m
        qa[i,:] = encode(cbooks[i], aac, distance=aq.distance)
        aac .-= decode(cbooks[i], qa[i,:])
    end
    return qa
end

function quantize_data(aq::ArrayQuantizer{Q,U,D,T,1}, aa::AbstractVector{T}) where {Q,U,D,T}
    nrows = aq.dims[1]
    @assert nrows == length(aa) "Quantized vector needs to have $nrows elements"
    aat = aa'  # use transpose as a single row matrix and quantize that
    aqt = ArrayQuantizer(aq.quantization, size(aat), codebooks(aq), aq.k, aq.distance, aq.rot)
    qat = quantize_data(aqt, aat)
    return qat
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/functions/nonlinear/reverse_pass.jl
# Functions used to compute reverse pass of nonlinear functions.
#############################################################################

# maximum number to perform reverse operation on associative term by summing
# and evaluating pairs remaining terms not reversed
const MAX_ASSOCIATIVE_REVERSE = 6
const REVERSE_DEBUG = false

"""
    reverse_plus_binary!

Updates storage tapes with reverse evalution of node representing `n = x + y` which updates x & y.
"""
function reverse_plus_binary!(k::Int64, children_arr::Vector{Int64}, children_idx::UnitRange{Int64},
                              numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N,T}},
                              x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, sparsity::Vector{Int},
                              subgrad_tol::Float64, is_post::Bool) where {N, T<:RelaxTag}

    REVERSE_DEBUG && println("--- start reverse plus binary ---")

    # extract values for k
    argk_is_number =  numvalued[k]
    if !argk_is_number
        setk =  setstorage[k]
    end

    # get row indices
    idx1 = first(children_idx)
    idx2 = last(children_idx)

    # extract values for argument 1
    arg1_index = children_arr[idx1]
    arg1_is_number = numvalued[arg1_index]
    if arg1_is_number
        set1 = zero(MC{N,T})
        num1 = numberstorage[arg1_index]
    else
        num1 = 0.0
        set1 = setstorage[arg1_index]
    end

    # extract values for argument 2
    arg2_index =  children_arr[idx2]
    arg2_is_number =  numvalued[arg2_index]
    if arg2_is_number
        num2 = numberstorage[arg2_index]
        set2 = zero(MC{N,T})
    else
        set2 = setstorage[arg2_index]
        num2 = 0.0
    end

    if !arg1_is_number && arg2_is_number
        c, a, b = IntervalContractors.plus_rev(setk.Intv, set1.Intv, num2)

    elseif arg1_is_number && !arg2_is_number
        c, a, b = IntervalContractors.plus_rev(setk.Intv, num1, set2.Intv)

    else
        c, a, b = IntervalContractors.plus_rev(setk.Intv, set1.Intv, set2.Intv)
    end

    #=
    if isnan(a) || isnan(b)
        if !arg1_is_number && arg2_is_number
            c, a, b = plus_rev(setk.Intv, MC{N,T}(set1.Intv), num2)

        elseif arg1_is_number && !arg2_is_number
            c, a, b = plus_rev(MC{N,T}(setk.Intv), num1, MC{N,T}(set2.Intv))

        else
            c, a, b = plus_rev(MC{N,T}(setk.Intv), MC{N,T}(set1.Intv), MC{N,T}(set2.Intv))
        end
    end
    =#

    if REVERSE_DEBUG
        println("val out = $(c)")
        println("arg1 out = $(a)")
        println("arg2 out = $(b)")
    end

    if !arg1_is_number
        if isempty(a)
            return false
        else
            setstorage[arg1_index] = MC{N,T}(a) #set1 ∩ a
        end
        REVERSE_DEBUG && println("setstorage[arg1_index] = $(setstorage[arg1_index])")
    end

    if !arg2_is_number
        if isempty(b)
            return false
        else
            setstorage[arg2_index] = MC{N,T}(b) # set2 ∩ b
        end
        REVERSE_DEBUG && println("setstorage[arg2_index] = $(setstorage[arg2_index])")
    end

    return true
end

"""
    reverse_plus_narity!

Updates storage tapes with reverse evalution of node representing `n = +(x,y,z...)` which updates x, y, z and so on.
"""
function reverse_plus_narity!(k::Int64, children_arr::Vector{Int64}, children_idx::UnitRange{Int64},
                              numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N,T}},
                              x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, sparsity::Vector{Int},
                              subgrad_tol::Float64, is_post::Bool) where {N, T<:RelaxTag}

    REVERSE_DEBUG && println("--- start reverse plus narity ---")

    continue_flag = true
    argk_is_number = numvalued[k]
    if !argk_is_number
        setk =  setstorage[k]
    end

    # out loops makes a temporary sum (minus one argument)
    # a reverse is then compute with respect to this argument
    active_count_number = 0
    for idx in children_idx
        active_idx = children_arr[idx]

        # don't contract a number valued argument
        active_arg_is_number = numvalued[active_idx]
        active_arg_is_number && continue

        if active_count_number >= MAX_ASSOCIATIVE_REVERSE
            break
        end

        tmp_sum = zero(MC{N,T})
        active_count_number += 1
        for nidx in children_idx
            inactive_idx =  children_arr[nidx]
            if inactive_idx != active_idx
                if  numvalued[inactive_idx]
                    tmp_sum += numberstorage[inactive_idx]
                else
                    tmp_sum += setstorage[inactive_idx]
                end
            end
        end

        active_set = setstorage[active_idx]
        c, a, b = IntervalContractors.plus_rev(setk.Intv, active_set.Intv, tmp_sum.Intv)

        if REVERSE_DEBUG
            println("val out = $(c)")
            println("arg1 out = $(a)")
            println("arg2 out = $(b)")
        end

        if isempty(a)
            return false
        else
            setstorage[active_idx] = MC{N,T}(a) # setk ∩ a
        end
    end

    return true
end

"""
    reverse_multiply_binary!

Updates storage tapes with reverse evalution of node representing `n = x * y` which updates x & y.
"""
function reverse_multiply_binary!(k::Int64, children_arr::Vector{Int64}, children_idx::UnitRange{Int64},
                                  numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N,T}},
                                  x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, sparsity::Vector{Int},
                                  subgrad_tol::Float64, is_post::Bool) where {N, T<:RelaxTag}

    REVERSE_DEBUG && println(" is_post = $is_post")

    REVERSE_DEBUG && println("--- start reverse mult binary ---")

    # extract values for k
    argk_is_number = numvalued[k]
    if !argk_is_number
        setk = setstorage[k]
        REVERSE_DEBUG && println("setk = $setk")
    end

    # get row indices
    idx1 = first(children_idx)
    idx2 = last(children_idx)

    # extract values for argument 1
    arg1_index = children_arr[idx1]
    arg1_is_number = numvalued[arg1_index]
    if arg1_is_number
        set1 = zero(MC{N,T})
        num1 =  numberstorage[arg1_index]
    else
        num1 = 0.0
        set1 = setstorage[arg1_index]
    end

    # extract values for argument 2
    arg2_index =  children_arr[idx2]
    arg2_is_number = numvalued[arg2_index]
    if arg2_is_number
        num2 = numberstorage[arg2_index]
        set2 = zero(MC{N,T})
    else
        set2 = setstorage[arg2_index]
        num2 = 0.0
    end

    if !arg1_is_number && arg2_is_number
        c, a, b = IntervalContractors.mul_rev(setk.Intv, set1.Intv, num2)

    elseif arg1_is_number && !arg2_is_number
        c, a, b = IntervalContractors.mul_rev(setk.Intv, num1, set2.Intv)

    else
        c, a, b = IntervalContractors.mul_rev(setk.Intv, set1.Intv, set2.Intv)
    end

    #=
    if isnan(a) || isnan(b)
        if !arg1_is_number && arg2_is_number
            c, a, b = mult_rev(MC{N,T}(setk.Intv), MC{N,T}(set1.Intv), num2)

        elseif arg1_is_number && !arg2_is_number
            c, a, b = mult_rev(MC{N,T}(setk.Intv), num1, MC{N,T}(set2.Intv))

        else
            c, a, b = mult_rev(MC{N,T}(setk.Intv), MC{N,T}(set1.Intv), MC{N,T}(set2.Intv))
        end
    end
    =#

    if REVERSE_DEBUG
        println("val out = $(c)")
        println("arg1 out = $(a)")
        println("arg2 out = $(b)")
    end

    if !arg1_is_number
        if isempty(a)
            return false
        else
            setstorage[arg1_index] = MC{N,T}(a) #set1 ∩ a
        end
        REVERSE_DEBUG && println("setstorage[arg1_index] = $(setstorage[arg1_index])")
    end

    if !arg2_is_number
        if isempty(b)
            return false
        else
            setstorage[arg2_index] = MC{N,T}(b) #set2 ∩ b
        end
        REVERSE_DEBUG && println("setstorage[arg2_index] = $(setstorage[arg2_index])")
    end

    return true
end

"""
    reverse_multiply_narity!

Updates storage tapes with reverse evalution of node representing `n = *(x,y,z...)` which updates x, y, z and so on.
"""
function reverse_multiply_narity!(k::Int64, children_arr::Vector{Int64}, children_idx::UnitRange{Int64},
                              numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N,T}},
                              x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, sparsity::Vector{Int},
                              subgrad_tol::Float64, is_post::Bool) where {N, T<:RelaxTag}

    REVERSE_DEBUG && println("--- start reverse mult narity ---")
    continue_flag = true
    argk_is_number = numvalued[k]
    if !argk_is_number
        setk = setstorage[k]
    end

    # out loops makes a temporary sum (minus one argument)
    # a reverse is then compute with respect to this argument
    active_count_number = 0
    for idx in children_idx
        active_idx = children_arr[idx]

        # don't contract a number valued argument
        active_arg_is_number = numvalued[active_idx]
        active_arg_is_number && continue

        if active_count_number >= MAX_ASSOCIATIVE_REVERSE
            break
        end

        tmp_mul = one(MC{N,T})
        active_count_number += 1
        for nidx in children_idx
            inactive_idx = children_arr[nidx]
            if inactive_idx != active_idx
                if  numvalued[inactive_idx]
                    tmp_mul *= numberstorage[inactive_idx]
                else
                    tmp_mul *= setstorage[inactive_idx]
                end
            end
        end

        active_set = setstorage[active_idx]
        c, a, b = IntervalContractors.mul_rev(setk.Intv, active_set.Intv, tmp_mul.Intv)

        if REVERSE_DEBUG
            println("val out = $(c)")
            println("arg1 out = $(a)")
            println("arg2 out = $(b)")
        end

        if isempty(a)
            return false
        else
            setstorage[active_idx] = MC{N,T}(a) #setstorage[active_idx] ∩ a
        end
    end

    return true
end

function reverse_minus!(k::Int64, children_arr::Vector{Int64}, children_idx::UnitRange{Int64},
                        numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N,T}},
                        x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, sparsity::Vector{Int},
                        subgrad_tol::Float64, is_post::Bool) where {N, T<:RelaxTag}

    REVERSE_DEBUG && println("--- start reverse minus ---")

    #argk_index =  children_arr[k]
    argk_is_number = numvalued[k]
    if !argk_is_number
        setk = setstorage[k]
    end

    # don't perform a reverse pass if the output was a number
    if argk_is_number
        return true
    end

    # get row indices
    idx1 = first(children_idx)
    idx2 = last(children_idx)

    # extract values for argument 1
    arg1_index = children_arr[idx1]
    arg1_is_number = numvalued[arg1_index]
    if arg1_is_number
        set1 = zero(MC{N,T})
        num1 = numberstorage[arg1_index]
    else
        num1 = 0.0
        set1 = setstorage[arg1_index]
    end

    # extract values for argument 2
    arg2_index = children_arr[idx2]
    arg2_is_number = numvalued[arg2_index]
    if arg2_is_number
        num2 = numberstorage[arg2_index]
        set2 = zero(MC{N,T})
    else
        set2 = setstorage[arg2_index]
        num2 = 0.0
    end

    if !arg1_is_number && arg2_is_number
        #c, a, b = minus_rev(setk, set1, num2)
        c, a, b = IntervalContractors.minus_rev(setk.Intv, set1.Intv, num2)

    elseif arg1_is_number && !arg2_is_number
        #c, a, b = minus_rev(setk, num1, set2)
        c, a, b = IntervalContractors.minus_rev(setk.Intv, num1, set2.Intv)

    else
        #c, a, b = minus_rev(setk, set1, set2)
        c, a, b = IntervalContractors.minus_rev(setk.Intv, set1.Intv, set2.Intv)
    end

    #=
    if isnan(a) || isnan(b)
        if !arg1_is_number && arg2_is_number
            c, a, b = minus_rev(MC{N,T}(setk.Intv), MC{N,T}(set1.Intv), num2)

        elseif arg1_is_number && !arg2_is_number
            c, a, b = minus_rev(MC{N,T}(setk.Intv), num1, MC{N,T}(set2.Intv))

        else
            c, a, b = minus_rev(MC{N,T}(setk.Intv), MC{N,T}(set1.Intv), MC{N,T}(set2.Intv))
        end
    end
    =#

    if REVERSE_DEBUG
        println("val out = $(c)")
        println("arg1 out = $(a)")
        println("arg2 out = $(b)")
    end

    if !arg1_is_number
        if isempty(a)
            return false
        else
            setstorage[arg1_index] = MC{N,T}(a) #set1 ∩ a
        end
        REVERSE_DEBUG && println("setstorage[arg1_index] = $(setstorage[arg1_index])")
    end

    if !arg2_is_number
        if isempty(b)
            return false
        else
            setstorage[arg2_index] = MC{N,T}(b) #set2 ∩ b
        end
        REVERSE_DEBUG && println("setstorage[arg2_index] = $(setstorage[arg2_index])")
    end

    return true
end

function reverse_power!(k::Int64, children_arr::Vector{Int64}, children_idx::UnitRange{Int64},
                        numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N,T}},
                        x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, sparsity::Vector{Int},
                        subgrad_tol::Float64, is_post::Bool) where {N, T<:RelaxTag}

    REVERSE_DEBUG && println("--- start reverse power ---")
    # extract values for k
    argk_is_number =  numvalued[k]
    if !argk_is_number
        setk = setstorage[k]
    end

    # don't perform a reverse pass if the output was a number
    if argk_is_number
        return true
    end

    # get row indices
    idx1 = first(children_idx)
    idx2 = last(children_idx)

    # extract values for argument 1
    arg1_index = children_arr[idx1]
    arg1_is_number = numvalued[arg1_index]
    if arg1_is_number
        set1 = zero(MC{N,T})
        num1 = numberstorage[arg1_index]
    else
        num1 = 0.0
        set1 = setstorage[arg1_index]
    end

    # extract values for argument 2
    arg2_index = children_arr[idx2]
    arg2_is_number = numvalued[arg2_index]
    if arg2_is_number
        num2 = numberstorage[arg2_index]
        set2 = zero(MC{N,T})
    else
        set2 = setstorage[arg2_index]
        num2 = 0.0
    end

    if !arg1_is_number && arg2_is_number
        #c, a, b = power_rev(setk, set1, num2)
        c, a, b = IntervalContractors.power_rev(setk.Intv, set1.Intv, num2)

    elseif arg1_is_number && !arg2_is_number
        #c, a, b = power_rev(setk, num1, set2)
        c, a, b = IntervalContractors.power_rev(setk.Intv, num1, set2.Intv)

    else
        #c, a, b = power_rev(setk, set1, set2)
        c, a, b = IntervalContractors.power_rev(setk.Intv, set1.Intv, set2.Intv)
    end

    #=
    if isnan(a) || isnan(b)
        if !arg1_is_number && arg2_is_number
            c, a, b = power_rev(MC{N,T}(setk.Intv), MC{N,T}(set1.Intv), num2)

        elseif arg1_is_number && !arg2_is_number
            c, a, b = power_rev(MC{N,T}(setk.Intv), num1, MC{N,T}(set2.Intv))

        else
            c, a, b = power_rev(MC{N,T}(setk.Intv), MC{N,T}(set1.Intv), MC{N,T}(set2.Intv))
        end
    end
    =#

    if REVERSE_DEBUG
        println("val out = $(c)")
        println("arg1 out = $(a)")
        println("arg2 out = $(b)")
    end

    if !arg1_is_number
        if isempty(a)
            return false
        else
            setstorage[arg1_index] = MC{N,T}(a) #setstorage[arg1_index] ∩ a
        end
    end

    if !arg2_is_number
        if isempty(b)
            return false
        else
            setstorage[arg2_index] = MC{N,T}(b) #setstorage[arg2_index] ∩  b
        end
    end

    return true
end

function reverse_divide!(k::Int64, children_arr::Vector{Int64}, children_idx::UnitRange{Int64},
                        numvalued::Vector{Bool}, numberstorage::Vector{Float64}, setstorage::Vector{MC{N,T}},
                        x::Vector{Float64}, lbd::Vector{Float64}, ubd::Vector{Float64}, sparsity::Vector{Int},
                        subgrad_tol::Float64, is_post::Bool) where {N, T<:RelaxTag}

    REVERSE_DEBUG && println("--- start reverse divide ---")

    # extract values for k
    argk_is_number =  numvalued[k]
    if !argk_is_number
        setk = setstorage[k]
    end

    # don't perform a reverse pass if the output was a number
    if argk_is_number
        return true
    end

    # get row indices
    idx1 = first(children_idx)
    idx2 = last(children_idx)

    # extract values for argument 1
    arg1_index = children_arr[idx1]
    arg1_is_number = numvalued[arg1_index]
    if arg1_is_number
        set1 = zero(MC{N,T})
        num1 = numberstorage[arg1_index]
    else
        num1 = 0.0
        set1 = setstorage[arg1_index]
    end

    # extract values for argument 2
    arg2_index = children_arr[idx2]
    arg2_is_number = numvalued[arg2_index]
    if arg2_is_number
        num2 = numberstorage[arg2_index]
        set2 = zero(MC{N,T})
    else
        set2 =  setstorage[arg2_index]
        num2 = 0.0
    end

    if !arg1_is_number && arg2_is_number
        c, a, b = IntervalContractors.div_rev(setk.Intv, set1.Intv, num2)

    elseif arg1_is_number && !arg2_is_number
        c, a, b = IntervalContractors.div_rev(setk.Intv, num1, set2.Intv)

    else
        c, a, b = IntervalContractors.div_rev(setk.Intv, set1.Intv, set2.Intv)
    end

    #=
    if isnan(a) || isnan(b)
        if !arg1_is_number && arg2_is_number
            c, a, b = div_rev(MC{N,T}(setk.Intv), MC{N,T}(set1.Intv), num2)

        elseif arg1_is_number && !arg2_is_number
            c, a, b = div_rev(MC{N,T}(setk.Intv), num1, MC{N,T}(set2.Intv))

        else
            c, a, b = div_rev(MC{N,T}(setk.Intv), MC{N,T}(set1.Intv), MC{N,T}(set2.Intv))
        end
    end
    =#

    if REVERSE_DEBUG
        println("val out = $(c)")
        println("arg1 out = $(a)")
        println("arg2 out = $(b)")
    end

    if !arg1_is_number
        if isempty(a)
            return false
        else
            setstorage[arg1_index] = MC{N,T}(a) #setstorage[arg1_index] ∩ a
        end
    end

    if !arg2_is_number
        if isempty(b)
            return false
        else
            setstorage[arg2_index] = MC{N,T}(b) #setstorage[arg2_index] ∩ b
        end
    end

    return true
end

function reverse_univariate!(k::Int64, op::Int64, arg_indx::Int64, setstorage::Vector{MC{N,T}}, x::Vector{Float64},
                             lbd::Vector{Float64}, ubd::Vector{Float64}, sparsity::Vector{Int},
                             subgrad_tol::Float64, is_post::Bool) where {N, T<:RelaxTag}

    REVERSE_DEBUG && println("--- start reverse reverse_univariate ---")
    valset = setstorage[k]
    argset = setstorage[arg_indx]

    REVERSE_DEBUG && println("val set = $(valset)")
    REVERSE_DEBUG && println("arg set = $(argset)")

    #a, b = eval_univariate_set_reverse(op, valset, argset)

    a, b = eval_univariate_set_reverse(op, valset.Intv, argset.Intv)

    REVERSE_DEBUG && println("val out = $(a)")
    REVERSE_DEBUG && println("arg out = $(b)")

    if isempty(b)
        return false
    else
        setstorage[arg_indx] = MC{N,T}(b) #setstorage[arg_indx] ∩ b
    end

    return true
end

function reverse_set_subexpression!(k::Int64, op::Int64, subexpressions::Vector{NonlinearExpression},
                                    numvalued::Vector{Bool}, numberstorage::Vector{Float64},
                                    setstorage::Vector{MC{N,T}}, cv_buffer::Vector{Float64},
                                    cc_buffer::Vector{Float64}, func_sparsity::Vector{Int64}) where {N, T<:RelaxTag}

    subexpression = subexpressions[op]
    isa_number = subexpression.is_number[1]
    numvalued[k] = isa_number

    return nothing
end

"""
    reverse_pass_kernel!

Performs a reverse McCormick/interval pass. If a NaN value is computed for the McCormick relaxation then the
routine defaults to the interval value instead.

There is a tacit assumption here that an abstract tree structure is used for propagation. If a cse structure is
eventually used then each reverse_xxx function will need to keep track of each variables state.
"""
function reverse_pass_kernel!(nd::Vector{JuMP.NodeData}, adj::SparseMatrixCSC{Bool,Int64}, x::Vector{Float64},
                              lbd::Vector{Float64}, ubd::Vector{Float64}, sparsity::Vector{Int},
                              setstorage::Vector{MC{N,T}}, subgrad_tol::Float64,
                              numberstorage::Vector{Float64}, numvalued::Vector{Bool},
                              is_post_input::Bool) where {N, T<:RelaxTag}

    children_arr = rowvals(adj)
    continue_flag = true
    is_post = is_post_input

    for k = 1:length(nd)

        REVERSE_DEBUG && println(" ")
        nod = nd[k]
        ntype = nod.nodetype
        nvalued = numvalued[k]

        if ntype == JuMP._Derivatives.VALUE      || ntype == JuMP._Derivatives.LOGIC     ||
           ntype == JuMP._Derivatives.COMPARISON || ntype == JuMP._Derivatives.PARAMETER ||
           ntype == JuMP._Derivatives.EXTRA      || ntype == JuMP._Derivatives.SUBEXPRESSION
           continue

        elseif nod.nodetype == JuMP._Derivatives.VARIABLE
            op = nod.index
            REVERSE_DEBUG && println("--- start reverse reverse variable ---")
            variable_interval = setstorage[k].Intv
            lower_interval = variable_interval.lo
            upper_interval = variable_interval.hi

            # update
            prior_x = x[op]
            if lower_interval > prior_x
                x[op] = 0.5*(lower_interval + upper_interval)
                is_post = false
            end
            if upper_interval < prior_x
                x[op] = 0.5*(lower_interval + upper_interval)
                is_post = false
            end
            lbd[op] = lower_interval
            ubd[op] = upper_interval

            REVERSE_DEBUG && println("variable_rev[$op][$k] at k = $k -> $(setstorage[k])")

        elseif nvalued
            continue

        elseif nod.nodetype == JuMP._Derivatives.CALL
            op = nod.index
            parent_index = nod.parent
            children_idx = nzrange(adj, k)
            parent_value = setstorage[k]
            n_children = length(children_idx)

            # SKIPS USER DEFINE OPERATORS NOT BRIDGED INTO JuMP Tree Representation
            if op >= JuMP._Derivatives.USER_OPERATOR_ID_START
                continue

            # :+
            elseif op === 1
                REVERSE_DEBUG && println("plus_rev[$n_children][$k]   at k = $k -> $(setstorage[k])")
                if n_children === 2
                    continue_flag &= reverse_plus_binary!(k, children_arr, children_idx, numvalued, numberstorage,
                                                          setstorage, x, lbd, ubd, sparsity, subgrad_tol, is_post)
                else
                    continue_flag &= reverse_plus_narity!(k, children_arr, children_idx, numvalued, numberstorage,
                                                          setstorage, x, lbd, ubd, sparsity, subgrad_tol, is_post)
                end
                REVERSE_DEBUG && !continue_flag && println("Infeasible node encountered.")

            # :-
            elseif op === 2
                REVERSE_DEBUG && println("minus_rev[$k]        at k = $k -> $(setstorage[k])")
                continue_flag &= reverse_minus!(k, children_arr, children_idx, numvalued, numberstorage,
                                                setstorage, x, lbd, ubd, sparsity, subgrad_tol, is_post)
                REVERSE_DEBUG && !continue_flag && println("Infeasible node encountered.")

            elseif op === 3 # :*
                REVERSE_DEBUG && println("mult_rev[$n_children][$k]     at k = $k -> $(setstorage[k])")
                if n_children === 2
                    continue_flag &= reverse_multiply_binary!(k, children_arr, children_idx, numvalued, numberstorage,
                                                              setstorage, x, lbd, ubd, sparsity, subgrad_tol, is_post)
                else
                    continue_flag &= reverse_multiply_narity!(k, children_arr, children_idx, numvalued, numberstorage,
                                                              setstorage, x, lbd, ubd, sparsity, subgrad_tol, is_post)
                end
                REVERSE_DEBUG && !continue_flag && println("Infeasible node encountered.")

             # :^
            elseif op === 4
                REVERSE_DEBUG && println("power_rev[$k]      at k = $k -> $(setstorage[k])")
                continue_flag &= reverse_power!(k, children_arr, children_idx, numvalued, numberstorage,
                                                setstorage, x, lbd, ubd, sparsity, subgrad_tol, is_post)
                REVERSE_DEBUG && !continue_flag && println("Infeasible node encountered.")

            # :/
            elseif op === 5
                REVERSE_DEBUG && println("power_div[$k]      at k = $k -> $(setstorage[k])")
                continue_flag &= reverse_divide!(k, children_arr, children_idx, numvalued, numberstorage,
                                                 setstorage, x, lbd, ubd, sparsity, subgrad_tol, is_post)
                REVERSE_DEBUG && !continue_flag && println("Infeasible node encountered.")

            # ifelse
            elseif op === 6
                continue
            end

        # assumes that child is set-valued and thus parent is set-valued (since isnumber already checked)
        elseif nod.nodetype == JuMP._Derivatives.CALLUNIVAR
            op = nod.index
            REVERSE_DEBUG && println("fop_rev[$op][$k]       at k = $k -> $(setstorage[k])")
            if op <= JuMP._Derivatives.USER_UNIVAR_OPERATOR_ID_START
                arg_indx =  children_arr[adj.colptr[k]]
                continue_flag &= reverse_univariate!(k, op, arg_indx, setstorage, x, lbd, ubd, sparsity,
                                                     subgrad_tol, is_post)
            end
            REVERSE_DEBUG && !continue_flag && println("Infeasible node encountered.")
        end

        !continue_flag && break
    end

    return continue_flag
end

"""
    reverse_pass!

A reverse_pass! on a `BufferedNonlinear` structure `d` intersects the existing value of the `d` with
constraint bounds then reverse propagates a set-valued operator (by default McCormick operator) along the
computational tape. The tapes are updated in place and boolean value is returned indicating whether the
reverse propagation yeilded a infeasible point (true = still feasible, false is proved infeasible).
"""
function reverse_pass!(evaluator::Evaluator, d::NonlinearExpression{V}) where V

    return reverse_pass_kernel!(d.nd, d.adj, evaluator.x, evaluator.lower_variable_bound,
                                evaluator.upper_variable_bound, d.grad_sparsity,
                                d.setstorage, evaluator.subgrad_tol, d.numberstorage,
                                d.isnumber, evaluator.reverse_subgrad_tighten)
end

function reverse_pass!(evaluator::Evaluator, d::BufferedNonlinearFunction{V}) where V
    d.last_past_reverse = true
    set_intersect_value!(d.expr, Interval(d.lower_bound, d.upper_bound))
    return reverse_pass!(evaluator, d.expr)
end


text:
"""
    Waveguide

A waveguide propagation path with a background `BField`, ionosphere `Species`, and `Ground`.
"""
abstract type Waveguide end

"""
    HomogeneousWaveguide{S} <: Waveguide

Defines a homogeneous segment of waveguide.

# Fields

- `bfield::BField`: background magnetic field.
- `species::S`: ionosphere constituents.
- `ground::Ground`: waveguide ground.
- `distance::Float64`: distance from the `Emitter` at the start of the segment in meters.
"""
struct HomogeneousWaveguide{S} <: Waveguide
    bfield::BField
    species::S  # `Species` or `Vector{Species}`
    ground::Ground
    distance::Float64
end

"""
    HomogeneousWaveguide(bfield, species, ground)

By default, `distance` is `0.0`.
"""
HomogeneousWaveguide(bfield, species::S, ground) where S =
    HomogeneousWaveguide{S}(bfield, species, ground, 0.0)

"""
    adjoint(w::HomogeneousWaveguide)

Return `w` with an adjoint `BField` having an `x` component of opposite sign.
"""
function adjoint(w::HomogeneousWaveguide)
    @unpack bfield, species, ground, distance = w
    adjoint_bfield = BField(bfield.B, -bfield.dcl, bfield.dcm, bfield.dcn)
    return HomogeneousWaveguide(adjoint_bfield, species, ground, distance)
end

"""
    SegmentedWaveguide{T<:Vector{<:Waveguide}} <: Waveguide

A collection of `Waveguide`s make up an inhomogeneous segmented waveguide.
"""
struct SegmentedWaveguide{T<:Vector{<:Waveguide}} <: Waveguide
    v::T
end

Base.length(w::SegmentedWaveguide) = length(w.v)
Base.sort!(w::SegmentedWaveguide) = sort!(w.v; by=x->getfield(x, :distance))
Base.sort(w::SegmentedWaveguide) = SegmentedWaveguide(sort!(copy(w.v)))
Base.issorted(w::SegmentedWaveguide) = issorted(w.v; by=x->getfield(x, :distance))
Base.size(w::SegmentedWaveguide) = size(w.v)
Base.getindex(w::SegmentedWaveguide, i::Int) = w.v[i]
Base.setindex(w::SegmentedWaveguide, x, i::Int) = (w.v[i] = x)
Base.push!(w::SegmentedWaveguide, x) = push!(w.v, x)
Base.iterate(w::SegmentedWaveguide) = iterate(w.v)
Base.iterate(w::SegmentedWaveguide, state) = iterate(w.v, state)

# TODO: WKBWaveguide (that's why we can't call SegmentedWaveguide -> InhomogeneousWaveguide)


text:
function _tostring(io; strip_summary=false)
    contents = map(readlines(seek(io,0))) do line
        replace(strip(line), "$Int" => "Int64")
    end

    # ignore summary
    strip_summary ? contents[2:end] : contents
end

@testset "enable/disable encoding" begin
    old_colormode = ImageInTerminal.colormode[1]
    old_should_render_image = ImageInTerminal.should_render_image[1]

    ImageInTerminal.enable_encoding()
    @test ImageInTerminal.colormode[1] == old_colormode
    @test ImageInTerminal.should_render_image[1] == true

    ImageInTerminal.enable_encoding()
    ImageInTerminal.disable_encoding()
    @test ImageInTerminal.colormode[1] == old_colormode
    @test ImageInTerminal.should_render_image[1] == false

    ImageInTerminal.disable_encoding()
    ImageInTerminal.use_256()
    @test ImageInTerminal.colormode[1] == ImageInTerminal.TermColor256()
    @test ImageInTerminal.should_render_image[1] == true

    ImageInTerminal.disable_encoding()
    ImageInTerminal.use_24bit()
    @test ImageInTerminal.colormode[1] == ImageInTerminal.TermColor24bit()
    @test ImageInTerminal.should_render_image[1] == true

    ImageInTerminal.colormode[1] = old_colormode
    ImageInTerminal.should_render_image[1] = old_should_render_image
end

@testset "no encoding" begin
    ImageInTerminal.disable_encoding()
    io = IOBuffer()
    img = fill(RGB(1.0, 1.0, 1.0), 4, 4)
    show(io, MIME"text/plain"(), img)
    @test_reference "reference/2d_show_raw.txt" _tostring(io; strip_summary=true)
    io = IOBuffer()
    show(io, MIME"text/plain"(), collect(rgb_line))
    @test_reference "reference/rgbline_show_raw.txt" _tostring(io; strip_summary=true)
    io = IOBuffer()
    show(io, MIME"text/plain"(), RGB(0.5,0.1,0.9))
    @test_reference "reference/colorant_show_raw.txt" _tostring(io)
end

@testset "256 colors" begin
    ImageInTerminal.use_256()
    io = IOBuffer()
    ensurecolor(show, io, MIME"text/plain"(), lena)
    @test_reference "reference/lena_show_256.txt" _tostring(io; strip_summary=true)
    io = IOBuffer()
    ensurecolor(show, io, MIME"text/plain"(), rgb_line)
    @test_reference "reference/rgbline_show_256.txt" _tostring(io; strip_summary=true)
    io = IOBuffer()
    ensurecolor(show, io, MIME"text/plain"(), RGB(0.5,0.1,0.9))
    @test_reference "reference/colorant_show_256.txt" _tostring(io)
end

@testset "24 bit" begin
    ImageInTerminal.use_24bit()
    io = IOBuffer()
    ensurecolor(show, io, MIME"text/plain"(), lena)
    @test_reference "reference/lena_show_24bit.txt" _tostring(io; strip_summary=true)
    io = IOBuffer()
    ensurecolor(show, io, MIME"text/plain"(), rgb_line)
    @test_reference "reference/rgbline_show_24bit.txt" _tostring(io; strip_summary=true)
    io = IOBuffer()
    ensurecolor(show, io, MIME"text/plain"(), RGB(0.5,0.1,0.9))
    @test_reference "reference/colorant_show_24bit.txt" _tostring(io)
end

ImageInTerminal.use_256()


text:
# import Pkg
# Pkg.add("YAML")
import YAML

const CONVENTIONAL_WEIGHT = Dict{Int,Float64}(
    1 => 1.008,
    3 => 6.94,
    5 => 10.81,
    6 => 12.011,
    7 => 14.007,
    8 => 15.999,
    12 => 24.305,
    14 => 28.085,
    16 => 32.06,
    17 => 35.45,
    35 => 79.904,
    81 => 204.38
)

const path = joinpath(
    dirname(@__FILE__), "../assets/raw/nist_atomic_weights20200410.txt")
const dest = joinpath(
    dirname(@__FILE__), "../assets/const/atomicweights.yaml")
const revdest = joinpath(
    dirname(@__FILE__), "../assets/const/symboltonumber.yaml")

function run()
    data = []
    revdict = Dict{String,Int}()
    lines = eachline(path)
    state = nothing
    next = iterate(lines, state)
    while true
        rcd = Dict()
        iso = Dict()
        while true
            (line, state) = next
            line == "" && break
            k, v = split(line, " = ")
            if k == "Atomic Number"
                rcd["Number"] = parse(Int, v)
            elseif k == "Atomic Symbol"
                rcd["Symbol"] = v
            elseif k == "Mass Number"
                iso["Number"] = parse(Int, v)
            elseif k == "Relative Atomic Mass"
                m = match(r"(.*?)\((.*?)#?\)", v)
                value = parse(Float64, m[1])
                unc = parse(Float64, string(
                    replace(m[1], r"[0-9]" => "0")[1:end-length(m[2])], m[2]))
                iso["Mass"] = value
                iso["MassUncertainty"] = unc
            elseif k == "Isotopic Composition"
                m = match(r"(.*?)\((.*?)#?\)", v)
                if m !== nothing
                    value = parse(Float64, m[1])
                    unc = parse(Float64, string(
                        replace(m[1], r"[0-9]" => "0")[1:end-length(m[2])], m[2]))
                elseif tryparse(Float64, v) !== nothing
                    value = parse(Float64, v)
                    unc = 0.0
                else
                    value = NaN
                    unc = NaN
                end
                iso["Composition"] = value
                iso["CompositionUncertainty"] = unc
            elseif k == "Standard Atomic Weight"
                m1 = match(r"\[(.*?),(.*?)\]", v)
                m2 = match(r"(.*?)\((.*?)\)", v)
                m3 = match(r"\[(.*?)\]", v)
                if m1 !== nothing
                    rcd["Weight"] = CONVENTIONAL_WEIGHT[rcd["Number"]]
                    rcd["WeightType"] = "Interval"
                    rcd["WeightLower"] = parse(Float64, m1[1])
                    rcd["WeightHigher"] = parse(Float64, m1[2])
                elseif m2 !== nothing
                    rcd["Weight"] = parse(Float64, m2[1])
                    unc = parse(Float64, string(
                        replace(m2[1], r"[0-9]" => "0")[1:end-length(m2[2])], m2[2]))
                    rcd["WeightType"] = "Uncertainty"
                    rcd["WeightUncertainty"] = unc
                elseif m3 !== nothing
                    rcd["Weight"] = parse(Int, m3[1])
                    rcd["WeightType"] = "MostStable"
                else
                    rcd["Weight"] = NaN
                    rcd["WeightType"] = "NotAvailable"
                end
            elseif k == "Notes"
                rcd["Notes"] = []
                v != " " && append!(rcd["Notes"], split(v, ","))
            end
            next = iterate(lines, state)
            next === nothing && break
        end
        revdict[rcd["Symbol"]] = rcd["Number"]
        if length(data) + 1 == rcd["Number"]
            if !isempty(data)
                e = data[end]
                # Monoisotopic mass
                comp = [
                    r["Composition"] for r in e["Isotopes"]
                    if r["Composition"] !== NaN]
                if !isempty(comp)
                    val, i = findmax(comp)
                    e["Monoisotopic"] = e["Isotopes"][i]["Mass"]
                    e["MonoisotopicUncertainty"
                        ] = e["Isotopes"][i]["MassUncertainty"]
                else
                    e["Monoisotopic"] = NaN
                    e["MonoisotopicUncertainty"] = NaN
                end
            end
            rcd["Isotopes"] = [iso]
            push!(data, rcd)
            # New entry
            if next === nothing
                e = data[end]
                # Monoisotopic mass
                comp = [
                    r["Composition"] for r in e["Isotopes"]
                    if r["Composition"] !== NaN]
                if !isempty(comp)
                    val, i = findmax(comp)
                    e["Monoisotopic"] = e["Isotopes"][i]["Mass"]
                    e["MonoisotopicUncertainty"
                        ] = e["Isotopes"][i]["MassUncertainty"]
                else
                    e["Monoisotopic"] = NaN
                    e["MonoisotopicUncertainty"] = NaN
                end
                break
            end
        else
            if rcd["Symbol"] != data[end]["Symbol"]
                iso["Symbol"] = rcd["Symbol"]
            end
            push!(data[end]["Isotopes"], iso)
        end
        next = iterate(lines, state)
    end
    println(length(data))
    YAML.write_file(dest, data)
    YAML.write_file(revdest, revdict)
end

run()


text:
import HTTP

dirs = eachline("dirs_wordlist.txt")

function scan_dirs(url)
    print("\nRunning common files and directory scan in $url\n")
    try
        for dir in dirs
            status = 0
            try
                request = HTTP.request("GET", "$url/$dir")
                status = request.status
            catch e
                if isa(e, HTTP.ExceptionRequest.StatusError)
                    status = e.status
                else
                    throw(e)
                end
            end
            if status == 200 || status == 400 || status == 403
                println(dir)
            end
        end
        return true
    catch e
        return false
    end
end

max_stars_repo_name: thecrypters/nis

max_issues_repo_name: thecrypters/nis

max_forks_repo_name: thecrypters/julia-host-scanner


text:
export rseis, wseis

# SeisIO file format version changes
# 0.54  2020-07-14    added Types: NodalChannel, NodalData, NodalLoc
# 0.53  2019-09-11    removed :i, :o from CoeffResp
#                     added :i, :o to MultiStageResp
# 0.52  2019-09-03    added Types: CoeffResp, MultiStageResp
# 0.51  2019-08-01    added :f0 to PZResp, PZResp64
# 0.50  all custom types can use write(); rseis, wseis no longer required
#       String arrays and :misc are written in a completely different way
#       Type codes for :misc changed
#       deprecated BigFloat/BigInt support in :misc
#       :n is no longer stored as a UInt32
#       :x compression no longer automatic and changed from Blosc to lz4
# ======== File versions below 0.5 are no longer supported ===================
# 0.41  :loc, :resp are now custom types with their own io subroutines
#       Julia version is no longer written to file
#       (likely misidentified in file header as 0.50)
# 0.40  SeisData.id[i] no longer needs to be a length ≤ 15 ASCII string
#       files have a new file TOC field: number of channels in each record
# 0.30  SeisData.loc[i] no longer is assumed to be length 5
# 0.20  First reliable file format

const TNames = Type[  EventChannel,
                      SeisChannel,
                      NodalChannel,
                      EventTraceData,
                      SeisData,
                      NodalData,
                      GenLoc,
                      GeoLoc,
                      UTMLoc,
                      XYLoc,
                      NodalLoc,
                      GenResp,
                      PZResp64,
                      PZResp,
                      CoeffResp,
                      MultiStageResp,
                      PhaseCat,
                      SeisEvent,
                      SeisHdr,
                      SeisPha,
                      SeisSrc,
                      SourceTime,
                      EQLoc,
                      EQMag ]

const TCodes = UInt32[ 0x20474330, # " GC0"  EventChannel
                       0x20474331, # " GC1"  SeisChannel
                       0x20474332, # " GC2"  NodalChannel
                       0x20474430, # " GD0"  EventTraceData
                       0x20474431, # " GD1"  SeisData
                       0x20474432, # " GD2"  NodalData
                       0x20495030, # " IP0"  GenLoc
                       0x20495031, # " IP1"  GeoLoc
                       0x20495032, # " IP2"  UTMLoc
                       0x20495033, # " IP3"  XYLoc
                       0x20495034, # " IP4"  NodalLoc
                       0x20495230, # " IR0"  GenResp
                       0x20495231, # " IR1"  PZResp64
                       0x20495232, # " IR2"  PZResp
                       0x20495233, # " IR3"  CoeffResp
                       0x20495234, # " IR4"  MultiStageResp
                       0x20504330, # " PC0"  PhaseCat = Dict{String, SeisPha}
                       0x20534530, # " SE0"  SeisEvent
                       0x20534830, # " SH0"  SeisHdr
                       0x20535030, # " SP0"  SeisPha
                       0x20535330, # " SS0"  SeisSrc
                       0x20535430, # " ST0"  SourceTime
                       0x45514c30, # "EQL0"  EQLoc
                       0x45514d30  # "EQM0"  EQMag
                    ]

#= check:
L = length(TCodes)
u = reinterpret(UInt8, TCodes)
S = Array{String,1}(undef, L)
for i = 1:L
  S[i] = reverse(join([Char(u[4*i-3]), Char(u[4*i-2]), Char(u[4*i-1]), Char(u[4*i])]))
end

should yield the character codes in the comments above
=#

# ===========================================================================
# Auxiliary file read functions
function read_rec(io::IO, ver::Float32, c::UInt32, b::UInt64)
  i = 0
  while i < length(TCodes)
    i = i + 1
    if c == getindex(TCodes, i)
      if (ver < vSeisIO) && (c == 0x20474431)
        return read_legacy(io, ver)
      else
        return read(io, getindex(TNames, i))
      end
    end
  end
  @warn("Non-SeisIO data at byte ", fastpos(io), "; skipped.")
  fastseek(io, b)
  return nothing
end

function build_file_list(patts::Union{String,Array{String,1}})
  plist = String[]
  if isa(patts, String)
    if safe_isfile(patts)
      return [patts]
    else
      plist = [patts]
    end
  else
    plist = patts
  end

  file_list = String[]
  for pat in plist
    if safe_isfile(pat)
      push!(file_list, pat)
    else
      append!(file_list, ls(pat))
    end
  end
  return file_list
end

"""
    rseis(fstr::String[, c::Array{Int64,1}=C, v::Integer=0, memmap::Bool=false])

Read SeisIO files matching file pattern ``fstr`` into memory. If an array of
record indices is passed to keyword c, only those record indices are read from
each file.

* Set v>0 to control verbosity.
* Set memmap=true to use memory mapping. Faster but potentially unsafe.
"""
function rseis(patts::Union{String,Array{String,1}};
  c::Union{Int64,Array{Int64,1}}  = Int64[],
  memmap::Bool                    = false,
  v::Integer                      = KW.v)

  A     = []
  files = build_file_list(patts)
  sbuf  = Array{UInt8, 1}(undef, 65535)
  chans = isa(c, Int64) ? [c] : c
  ver   = vSeisIO
  L     = zero(Int64)

  nf = 0
  for f in files
    nf  = nf + 1
    io  = memmap ? IOBuffer(Mmap.mmap(f)) : open(f, "r")

    # All SeisIO files begin with "SEISIO"
    if fastread(io, 6) != UInt8[0x53, 0x45, 0x49, 0x53, 0x49, 0x4f]
      @warn string("Skipped ", f, ": not a SeisIO file!")
      close(io)
      continue
    end

    ver = fastread(io, Float32)
    L   = fastread(io, Int64)
    C   = read!(io, Array{UInt32,1}(undef, L))
    B   = read!(io, Array{UInt64,1}(undef, L))

    # DND -- faster to avoid seek
    @inbounds if isempty(chans)
      (v > 0) && println("Reading ", L, " objects from ", f)
      (v > 1) && println("C = ", C)
      for n = 1:L
        (v > 1) && println("Reading object with code ", repr(getindex(C, n)), " (", n, "/", L, ")")
        R = read_rec(io, ver, getindex(C, n), getindex(B, n == L ? n : n+1))
        push!(A, R)
        (v > 1) && println("Read ", typeof(getindex(A, n)), " object (", n, "/", L, ")")
      end

    else
      if minimum(chans) > L
        (v > 0) && println("Skipped ", f, ": no matching record numbers.")
        close(io)
        continue
      end
      for n in chans
        if n in 1:L
          fastseek(io, getindex(B, n))
          R = read_rec(io, ver, getindex(C, n), getindex(B, n == L ? n : n+1))
          push!(A, R)
          (v > 1) && println("Read ", typeof(last(A)), " object from ", f, ", bytes ", getindex(B, n), ":", ((n == L) ? fastpos(io) : getindex(B, n+1)))
        else
          (v > 0) && println(n > L ? "No" : "Skipped", " record ", n, " in ", f)
        end
      end
    end

    close(io)
  end
  (v > 0) && println("Processed ", nf, " files.")
  return A
end

"""
    wseis(fname, S)

Write SeisIO objects S to file. S can be a single object, multiple comma-delineated objects, or an array of objects.
"""
function wseis(fname::String, S...)
    L = Int64(lastindex(S))
    (L == zero(Int64)) && return nothing

    C = zeros(UInt32, L)                  # Codes
    B = zeros(UInt64, L)                  # Byte indices
    ID = Array{UInt64,1}(undef, 0)        # IDs
    TS = Array{Int64,1}(undef, 0)         # Start times
    TE = Array{Int64,1}(undef, 0)         # End times
    P  = Array{Int64,1}(undef, 0)         # Parent object indices in C, B

    # Buffer checks
    checkbuf!(BUF.int64_buf, 8)

    # open file for writing
    io = open(fname, "w")
    write(io, "SEISIO")
    write(io, vSeisIO)
    write(io, L)
    p = fastpos(io)
    skip(io, sizeof(C) + sizeof(B))

    # Write all objects
    i = 0
    @inbounds while i < L
      i = i + 1
      setindex!(B, UInt64(fastpos(io)), i)
      seis = getindex(S, i)
      write(io, seis)
      T = typeof(seis)

      # store type code to C
      j = 0
      while j < length(TNames)
        j = j + 1
        if T == getindex(TNames, j)
          setindex!(C, getindex(TCodes, j), i)
          break
        end
      end

      # Add to id, time index; log to notes
      if T <: GphysChannel
        push!(ID, hash(getfield(seis, :id)))
        push!(P, i)
        fs = getfield(seis, :fs)
        if !isempty(seis.t)
          t = getfield(seis, :t)
          push!(TS, starttime(t, seis.fs))
          push!(TE, fs == zero(Float64) ? t[end,2] : sum(t, dims=1)[2] +
            round(Int64, sμ*lastindex(seis.x)/fs))
        end
        fwrite_note!(seis, "wseis", fname, "")
      elseif T <: GphysData
        append!(ID, hash.(getfield(seis, :id)))
        append!(P, ones(Int64, seis.n).*i)
        k = 0
        ts = Array{Int64,1}(undef,seis.n)
        te = similar(ts)
        @inbounds while k < seis.n
          k = k + 1
          t = seis.t[k]
          fs = seis.fs[k]
          if !isempty(seis.t)
            setindex!(ts, starttime(t, fs), i)
            setindex!(te, fs == zero(Float64) ? t[end,2] : sum(t, dims=1)[2] +
              round(Int64, sμ*lastindex(seis.x)/fs), i)
          end
          fwrite_note!(seis, k, "wseis", fname, "")
        end
        append!(TS, ts)
        append!(TE, te)
      elseif T in (SeisHdr, SeisSrc)
        fwrite_note_quake!(seis, "wseis", fname, "")
      end
    end

    # Write TOC.
    # format: array of object types, array of byte indices
    fastseek(io, p)
    write(io, C)
    write(io, B)

    # File index added 2017-02-23; changed 2019-05-27
    # index format: ID hash, TS, TE, P, positions
    fastseekend(io)
    b = zeros(Int64, 4)
    if isempty(ID)
      fill!(b, fastpos(io))
    else
      setindex!(b, Int64(fastpos(io)), 1)
      write(io, ID)
      setindex!(b, Int64(fastpos(io)), 2)
      write(io, TS)
      setindex!(b, Int64(fastpos(io)), 3)
      write(io, TE)
      setindex!(b, Int64(fastpos(io)), 4)
      write(io, P)
    end
    write(io, b)
    close(io)
    return nothing
end


text:
# using Distributed
using JLD
using DelimitedFiles
using Plots
using LinearAlgebra
# gr()
# using Arpack
ENV["GKSwstype"]="100"
gr()

# @everywhere begin
#     # using ProgressMeter
#     using Distributed
#     using SharedArrays

    include("kicked_Ising_ed.jl")

    function otoc_one_shot(Jz::Float64, Jx::Float64, hz::Float64, Lt::Int64, Lr::Int64, phi_arr::Array{Float64, 1})
        Jz_arr = Jz*ones(Lr)
        hz_arr = hz*ones(Lr)
        Jx_arr = Jx*ones(Lr)

        hamJz = Ham_Jz_ops(Jz_arr, Lr)
        hamhz = Ham_hz_ops(hz_arr, Lr)
        hamJx = Ham_Jx_ops(Jx_arr, Lr)

        UF = exp(1im*(hamJz+hamhz)) * exp(1im*(hamJx))

        # observable
        # at site 1

        op1 = zeros(ComplexF64, 2^Lr, 2^Lr)
        a = 3
        term_ops = fill(Idmatrix, Lr)
        term_ops[1] = PauliMatrices[a]
        op1 = foldl(⊗, term_ops)

        # boundary operator

        # phi = philist[ii]
        # Z_op_sum = Ham_hz_ops(ones(Lr), Lr)
        # boundary_op = exp(i*phi*Z_op_sum)

        # Calculate different version of OTOC
        Utemp = Matrix{ComplexF64}(I, 2^Lr, 2^Lr)

        otoc_type_list = zeros(Float64, Lt+1, length(phi_arr), 4)
        saturate_val_list = zeros(Float64, length(phi_arr))
        # Type 1 Tr(e^{iHt} a e^{iHt} e^{iϕW}e^{iHt}a e^{iHt}e^{-iϕW})
        # Type 2 Tr(e^{iH(t+1)} a e^{iHt} e^{iϕW}e^{iHt}a e^{iHt}e^{-iϕW})
        # Type 3 Tr(e^{iHt} a e^{iHt} e^{iϕW}e^{iHt}a e^{iHt}e^{-iϕW})

        Z_op_sum = Ham_hz_ops(ones(Lr), Lr)

        for t in 1:Lt+1
            print("t=$(t), ")

            # Type 1
            u1 = Utemp*UF*op1*adjoint(Utemp)
            u2 = adjoint(UF)*Utemp*op1*adjoint(Utemp)
            for (phi_index, phi) in enumerate(phi_arr)
                # boundary operator
                boundary_op = exp(1im*phi*Z_op_sum)
                otoc_type_list[t, phi_index, 1] += real(tr(u1*boundary_op*u2*adjoint(boundary_op) ) ) / (2^Lr)
            end

            # Type 2
            u1 = Utemp*UF*op1*adjoint(UF)*op1*adjoint(Utemp)*UF
            u2 = adjoint(UF)
            for (phi_index, phi) in enumerate(phi_arr)
                # boundary operator
                boundary_op = exp(1im*phi*Z_op_sum)
                otoc_type_list[t, phi_index, 2] += real(tr(u1*boundary_op*u2*adjoint(boundary_op) ) ) / (2^Lr)
            end

            # Type 3 all four e^{iHt}
            u1 = Utemp*op1*Utemp
            u2 = Utemp*op1*Utemp
            for (phi_index, phi) in enumerate(phi_arr)
                # boundary operator
                boundary_op = exp(1im*phi*Z_op_sum)
                otoc_type_list[t, phi_index, 3] += real(tr(u1*boundary_op*u2*adjoint(boundary_op) ) ) / (2^Lr)
            end

            # Type 4 normal OTOC
            u1 = Utemp*op1*adjoint(Utemp)
            u2 = Utemp*op1*adjoint(Utemp)
            for (phi_index, phi) in enumerate(phi_arr)
                # boundary operator
                boundary_op = exp(1im*phi*Z_op_sum)
                otoc_type_list[t, phi_index, 4] += real(tr(u1*boundary_op*u2*adjoint(boundary_op) ) ) / (2^Lr)
            end

            Utemp = Utemp*UF

            for (phi_index, phi) in enumerate(phi_arr)
                # boundary operator
                boundary_op = exp(1im*phi*Z_op_sum)
                saturate_val_list[phi_index] = real(tr(UF*boundary_op*adjoint(UF)*adjoint(boundary_op))/(2^Lr) * tr(UF*op1*adjoint(UF)*op1)/(2^Lr)  )
            end
        end

        print("")

        return otoc_type_list, saturate_val_list
    end

# end

Lr = 10
Lt = 15
phi_arr = [0.0, 0.2, 0.4, 0.6]

data_dirname = "dual_circuits_otoc_20211125"
if !isdir(data_dirname)
    mkdir(data_dirname)
end

# Jz_list = [2.0]
# hz_list = [0.9]
# Jx_list = [0.5]

Jz_list = []
hz_list = []
Jx_list = []

for i in 1:3
    for j in 1:3
        for z in 1:3
            append!(Jz_list, 1.5+0.5*i)
            append!(hz_list, 0.5+0.2*j)
            append!(Jx_list, 0.3+0.1*z)
        end
    end
end

otoc_type_arr_data = zeros(Lt+1, length(phi_arr), 4, length(Jz_list))
saturate_val_list_data = zeros(length(phi_arr), length(Jz_list))

#
#
# otoc_aver_listall = JLD.load("otoc_aver.jld", "otoc_aver")

t_list = collect(0:Lt)

label_phi = Array{String, 2}(undef, 1, length(phi_arr))
for ll in 1:length(phi_arr)
    label_phi[1,ll] = "\$ \\phi=$(phi_arr[ll]) \$"
end
using LaTeXStrings

for i in 1:length(Jz_list)
    Jz = Jz_list[i]
    hz = hz_list[i]
    Jx = Jx_list[i]
    println("index=$(i), J_z=$(Jz), h_z=$(hz), J_x=$(Jx)")

    otoc_type_arr_data[:,:,:,i], saturate_val_list_data[:,i] =  otoc_one_shot(Jz, Jx, hz, Lt, Lr, phi_arr)
    JLD.save("$(data_dirname)/otoc_Jz$(Jz)_hz$(hz)_Jx$(Jx)_Lr$(Lr)_Lt$(Lt).jld", "otoc_type_arr_data",otoc_type_arr_data)

    p1 = plot(t_list, otoc_type_arr_data[:,:,1,i], title = "\$ O(t+\\Delta t) W(\\Delta t) O(t) W(0) \$", xlabel="\$ L_t \$",ylabel=L"F(\phi, L_t)",    markershape = :auto, markersize = 3, label=label_phi, titlefont=9, legendfontsize=6, xguidefontsize=6,yguidefontsize=6)
    plot!(saturate_val_list_data[:,i], seriestype="hline")
    p2 = plot(t_list, otoc_type_arr_data[:,:,2,i], title = "\$ O(t+\\Delta t) O(t) W(\\Delta t) W(0) \$", xlabel="\$ L_t \$",ylabel=L"F(\phi, L_t)",    markershape = :auto, markersize = 3, label=label_phi, titlefont=9, legendfontsize=6, xguidefontsize=6,yguidefontsize=6)
    plot!(saturate_val_list_data[:,i], seriestype="hline")
    p3 = plot(t_list, otoc_type_arr_data[:,:,3,i], title = "\$ Tr[e^{iHt}O e^{iHt}W e^{iHt} O e^{iHt}W] \$", xlabel="\$ L_t \$",ylabel=L"F(\phi, L_t)", markershape = :auto, markersize = 3, label=label_phi, titlefont=9, legendfontsize=6, xguidefontsize=6,yguidefontsize=6)
    p4 = plot(t_list, otoc_type_arr_data[:,:,4,i], title = "\$ O(t) W(0) O(t) W(0) \$", xlabel="\$ L_t \$",ylabel=L"F(\phi, L_t)",                      markershape = :auto, markersize = 3, label=label_phi, titlefont=9, legendfontsize=6, xguidefontsize=6,yguidefontsize=6)
    # plot(p1,p2,p3,p4, layout=(2,2), title = "\$ J_z=$(Jz), h_z=$(hz), J_x=$(Jx) \$")
    plot(p1,p2,p3,p4, layout=(2,2))

    # JLD.save("$(data_dirname)/otoc_aver.jld", "otoc_aver", otoc_aver_listall)
    savefig("$(data_dirname)/otoc_Jz$(Jz)_hz$(hz)_Jx$(Jx)_Lr$(Lr)_Lt$(Lt).pdf")
    # savefig("thermalotoc.pdf")
end