text:
module OptimizationAlgorithms
using LinearAlgebra

# principle of this package:
# reduce optimization problems to fixed point iterations
# to this end, design fixed point operator f! and stopping rule isfixed

######################## General Fixed Point Algorithm #########################
# IDEA: allocating (for immutable types), via 0D array?
# attempts to find a fixed point of f! by iterating f! on x
# in light of ode's, this could be a gradient flow
# general enough: iterate! - ?
# f! should be a fixed point operator, so ∃x s.t. x = f!(x)
function fixedpoint!(f!, x, isfixed)
    t = 1
    while true
        x = f!(x, t)
        t += 1
        isfixed(x, t) && break
    end
    x, t
end

# using default stopping criterion
fixedpoint!(f!, x) = fixedpoint!(f!, x, StoppingCriterion(x))

# # TODO:
# # f!(x, t) should return both x and y for generality with scalar case
# this is not a problem for Directions, but for Updates
# either rewrite update! to do this, thereby breaking it's use in fixedpoint!,
# or write new function update_evaluate!()
# function optimize!(f!, x, isfixed)
#     t = 0
#     while true
#         x, y = f!(x, t)
#         t += 1
#         isfixed(x, y, t) && break
#     end
#     x, y, t
# end
# optimize!(f!, x) = optimize!(f!, x, StoppingCriterion(x))

########################### Stopping Criterion ##################################
# IDEA: have input output be single argument as tuple?
# abstract type StoppingCriterion{T} end
mutable struct StoppingCriterion{T, S}
    x::T # holds last value of parameters
    dx::S # minimum absolute change in parameters x for termination
    rx::S # minimum relative change in parameters x for termination

    y::S # holds last function value
    dy::S
    ry::S

    maxiter::Int # maximum iterations
    verbose::Bool
    function StoppingCriterion(x, y::Real = Inf; dx::Real = 1e-6, rx::Real = 0.,
        dy::Real = 1e-6, ry::Real = 0., maxiter::Int = 1024, verbose::Bool = false)
        new{typeof(x), typeof(y)}(copy(x), dx, rx, y, dy, ry, maxiter, verbose)
    end
end

(T::StoppingCriterion)(x) = input_criterion!(T, x)
function (T::StoppingCriterion)(x, t::Integer)
    T.verbose && println("in iteration $t:")
    1 < t && (input_criterion!(T, x) || T.maxiter < t)
end
function (T::StoppingCriterion)(x, y::Real, t::Integer)
    T.verbose && println("in iteration $t:")
    1 < t && (input_criterion!(T, x) || output_criterion!(T, y) || T.maxiter < t)
end

function input_criterion!(T::StoppingCriterion, x::Union{Real, AbstractVector})
    any(isnan, x) && error("x contains NaN: x = $x")
    dx = norm(difference(x, T.x))
    rx = dx / norm(x)
    if T.verbose
        println("absolute change in x = $dx")
        println("relative change in x = $rx")
    end
    isfixed = dx < T.dx || rx < T.rx
    if x isa Real
        T.x = x
    else # AbstractVector
        T.x .= x
    end
    return isfixed
end
function output_criterion!(T::StoppingCriterion, y::Real)
    isnan(y) && error("y = NaN")
    dy = abs(y - T.y)
    ry = dy / abs(y)
    if T.verbose
        println("absolute change in y = $dy")
        println("relative change in y = $ry")
    end
    isfixed = dy < T.dy || ry < T.ry
    T.y = y
    return isfixed
end

######################## Canonical Descent Directions ##########################
# abstract type FixedPointIterator{T} end # makes sense as long as we stick to optimization
abstract type Update{T} end # <: FixedPointIterator{T} end
abstract type Direction{T} <: Update{T} end
# NOTE: all directions are decreasing unless a negative step size is chosen

# update functor
(D::Update)(x, t::Int) = update!(D, x, t)
(D::Update)(x) = update!(D, x)

# direction functor
(D::Direction)(x, t::Int) = direction(D, x) # since directions do not depend on time
(D::Direction)(x) = direction(D, x)

update!(D::Update, x, t::Int) = update!(D, x)
update!(D::Direction, x) = (x .+= D(x))
update!(D::Direction, x, t::Int) = (x .+= D(x, t))

# converting a direction to an update
update!(D::Direction) = (x, t::Int) -> (x .+= D(x, t))
fixedpoint!(D::Direction, x, isfixed) = fixedpoint!(update!(D), x, isfixed)

function objective end # returns objective function
objective(D::Direction, x) = D.f(x)
objective(D::Direction) = x->objective(D, x)

const value = objective

function direction end # returns direction
direction(D::Direction, x) = valdir(D, x)[2]
# direction(D::Direction, x, t::Int) = valdir(D, x)[2]

# calculates direction and value at the same time
# value_direction(D::Direction, x, t::Int) = valdir(D::Direction, x)
value_direction(D::Direction) = x -> valdir(D, x)
const valdir = value_direction

################################################################################
include("util.jl")
include("descentdirections.jl")
include("gaussnewton.jl")
include("descentupdates.jl")
include("stochastic.jl")
include("ode.jl")
include("stepsize.jl")
include("submodular.jl")

end # Optimization


text:
#==========================================================================================#

# AUXILIARY FUNCTIONS

macro component!(x)
    SM.terms!(SM.sort_terms!(SM.parse!(x)))
end

function component(x::Expr)

    (x.args[1] == :(=>)) || throw(ArgumentError("invalid expression encountered: $x"))

    y = @eval @component! $(x.args[3])

    if (typeof(y) <: Tuple) && (typeof(y[1]) <: SM.ConstantTerm)
        y = (y[2:end]..., y[1])
    end

    return y
end

macro micromodel(x...)
    Dict{Symbol, SM.TermOrTerms}(xi.args[2] => component(xi) for xi in x)
end

terms(y::Dict{Symbol, SM.TermOrTerms})    = mapreduce(terms, union, values(y))
termvars(y::Dict{Symbol, SM.TermOrTerms}) = mapreduce(termvars, union, values(y))
eterms(y::Dict{Symbol, SM.TermOrTerms})   = Tuple(mapreduce(SM.vectorize, union, values(y)))

function schema(y::Dict{Symbol, SM.TermOrTerms}, df::DataFrame, hints::Dict)
    schema(filter(SM.needs_schema, terms(y)), df, hints)
end

function apply_schema(y::Dict{Symbol, SM.TermOrTerms}, schema::SM.Schema)
    Dict{Symbol, SM.TermOrTerms}(i => apply_schema(v, schema) for (i, v) in y)
end

has_schema(y::Dict{Symbol, SM.TermOrTerms}) = all(has_schema, values(y))

function modelcols(y::Dict{Symbol, SM.TermOrTerms}, df::DataFrame)
    modelcols(SM.collect_matrix_terms(eterms(y)), df)
end

function asgn(y::Dict{Symbol, SM.TermOrTerms})

    Y = eterms(y)
    Y = Y[SM.asgn(Y)]

    return Dict{Symbol, Vector{Int}}(i => findall((in)(SM.vectorize(v)), Y) for (i, v) in y)
end

function asgn(y::Dict{Symbol, SM.TermOrTerms}, t::SM.TermOrTerms)

    Y = eterms(y)
    Y = Y[SM.asgn(Y)]

    return [findfirst(isequal(ti), Y) for ti in t]
end

#==========================================================================================#

# TYPE

mutable struct Microdata
    nonmissing::BitVector
    model::Dict{Symbol, SM.TermOrTerms}
    mapping::Dict{Symbol, Vector{Int}}
    matrix::Matrix{Float64}
    weights::AbstractWeights{Int, Float64, Vector{Float64}}
    corr::CorrStructure
end

# CONSTRUCTORS

function Microdata(
        df::DataFrame, model::Dict{Symbol, SM.TermOrTerms};
        hints::Dict{Symbol} = Dict{Symbol, Any}(),
        subset::AbstractVector{Bool} = trues(size(df, 1)),
        vcov::CorrStructure = Heteroscedastic(),
        weights::AbstractWeights = UnitWeights(Float64, size(df, 1))
    )

    nonmissing  = completecases(df, termvars(model))
    nonmissing .= nonmissing .* .!iszero.(weights) .* subset
    corr        = parse_corr!(nonmissing, vcov)
    weights     = parse_weights(nonmissing, weights)
    columns     = termvars(model)
    data        = DataFrame(map(x -> disallowmissing(df[nonmissing, x]), columns), columns)
    model       = apply_schema(model, schema(model, data, hints))
    mapping     = asgn(model)
    matrix      = modelcols(model, data)
    mapping     = asgn(model)

    return Microdata(nonmissing, model, mapping, matrix, weights, corr)
end

function Microdata(MD::Microdata)
    Microdata(
        copy(MD.nonmissing),
        copy(MD.model),
        copy(MD.mapping),
        MD.matrix,
        MD.weights,
        MD.corr
    )
end

#==========================================================================================#

# COPY

copy(m::Microdata) = Microdata([copy(getfield(m, k))     for k in fieldnames(m)]...)

ext: jl


text:
#=
solve a simple LMI problem

min y
st. sum(x) = 1
    Qₖ + y I + sumᵢ xᵢ Pₖᵢ ⪰ 0, ∀k = 1..K
where Qₖ ≺ 0, Pₖᵢ ∈ 𝕊ˢ

formulations: PSD, sparse PSD (sparse Pₖᵢ, with sparse or dense chol),
LMI (sparse or dense Pₖᵢ)
=#

using SparseArrays

struct SparseLMIJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    num_lmis::Int
    side_Ps::Int
    num_Ps::Int
    sparse_Ps::Bool
    use_psd::Bool # use PSD cone
    use_linmatrixineq::Bool # use linmatrixineq cone
    use_sparsepsd::Bool # use sparsePSD cone
    use_cholmod_impl::Bool # use sparsePSD cone with CHOLMOD, else dense impl
end

function build(inst::SparseLMIJuMP{T}) where {T <: Float64}
    @assert +(inst.use_psd, inst.use_sparsepsd, inst.use_linmatrixineq) == 1
    (num_lmis, side_Ps, num_Ps) = (inst.num_lmis, inst.side_Ps, inst.num_Ps)

    function rand_symm()
        if inst.sparse_Ps
            sparsity = min(3.0 / side_Ps, 1.0) # sparsity factor
            M = sprandn(side_Ps, side_Ps, sparsity)
            for idx in rand(1:side_Ps, div(side_Ps, 3))
                M[idx, idx] = rand()
            end
            return Symmetric(M)
        else
            return Symmetric(randn(side_Ps, side_Ps))
        end
    end
    rand_psd() = (M = rand_symm(); Symmetric(M * M'))

    Ps = [rand_symm() for k in 1:num_lmis, i in 1:num_Ps]
    Qs = [-rand_psd() for k in 1:num_lmis]
    matI = Symmetric(one(Qs[1]))

    model = JuMP.Model()
    JuMP.@variable(model, y)
    JuMP.@objective(model, Min, y)
    JuMP.@variable(model, x[1:num_Ps])
    JuMP.@constraint(model, sum(x) == 1)

    if inst.use_psd || inst.use_sparsepsd
        for k in 1:num_lmis
            Sk = Symmetric(Qs[k] + y * matI +
                sum(x[i] * Ps[k, i] for i in 1:num_Ps))
            if inst.use_psd
                JuMP.@constraint(model, Sk in JuMP.PSDCone())
            else
                impl = (inst.use_cholmod_impl ? Cones.PSDSparseCholmod :
                    Cones.PSDSparseDense)
                cone = Hypatia.PosSemidefTriSparseCone{impl, T, T}
                (row_idxs, col_idxs, vals) = findnz(sparse(LowerTriangular(Sk)))
                JuMP.@constraint(model, vals in cone(
                    side_Ps, row_idxs, col_idxs, false))
            end
        end
    elseif inst.use_linmatrixineq
        JuMP.@constraint(model, [k in 1:num_lmis], vcat(y, x, 1) in
            Hypatia.LinMatrixIneqCone{T}([matI, Ps[k, :]..., Qs[k]]))
    else
        error()
    end

    return model
end


text:
const default_sound = :boing
const default_error = :blip

macro beep(ex)
    return quote
        try
            local res = $(esc(ex))
            beep(default_sound)
            return res
        catch
            beep(default_error)
            rethrow()
        end
    end
end

macro beep(ex, sound)
    return quote
        try
            local res = $(esc(ex))
            beep($sound)
            return res
        catch
            beep(default_error)
            rethrow()
        end
    end
end

macro beep(ex, sound, error_sound)
    return quote
        try
            local res = $(esc(ex))
            beep($sound)
            return res
        catch
            beep($error_sound)
            rethrow()
        end
    end
end


text:
export DiffEqSystem, ODEFunction

using Base: RefValue

isintermediate(eq::Equation) = !(isa(eq.lhs, Operation) && isa(eq.lhs.op, Differential))

function flatten_differential(O::Operation)
    @assert is_derivative(O) "invalid differential: $O"
    is_derivative(O.args[1]) || return (O.args[1], O.op.x, 1)
    (x, t, order) = flatten_differential(O.args[1])
    isequal(t, O.op.x) || throw(ArgumentError("non-matching differentials on lhs: $t, $(O.op.x)"))
    return (x, t, order + 1)
end

struct DiffEq  # dⁿx/dtⁿ = rhs
    x::Expression
    t::Variable
    n::Int
    rhs::Expression
end
function Base.convert(::Type{DiffEq}, eq::Equation)
    isintermediate(eq) && throw(ArgumentError("intermediate equation received"))
    (x, t, n) = flatten_differential(eq.lhs)
    return DiffEq(x, t, n, eq.rhs)
end
Base.:(==)(a::DiffEq, b::DiffEq) = isequal((a.x, a.t, a.n, a.rhs), (b.x, b.t, b.n, b.rhs))
get_args(eq::DiffEq) = Expression[eq.x, eq.t, eq.rhs]

struct DiffEqSystem <: AbstractSystem
    eqs::Vector{DiffEq}
    iv::Variable
    dvs::Vector{Variable}
    ps::Vector{Variable}
    jac::RefValue{Matrix{Expression}}
    function DiffEqSystem(eqs, iv, dvs, ps)
        jac = RefValue(Matrix{Expression}(undef, 0, 0))
        new(eqs, iv, dvs, ps, jac)
    end
end

function DiffEqSystem(eqs)
    dvs, = extract_elements(eqs, [_is_dependent])
    ivs = unique(vcat((dv.dependents for dv ∈ dvs)...))
    length(ivs) == 1 || throw(ArgumentError("one independent variable currently supported"))
    iv = first(ivs)
    ps, = extract_elements(eqs, [_is_parameter(iv)])
    DiffEqSystem(eqs, iv, dvs, ps)
end

function DiffEqSystem(eqs, iv)
    dvs, ps = extract_elements(eqs, [_is_dependent, _is_parameter(iv)])
    DiffEqSystem(eqs, iv, dvs, ps)
end

function calculate_jacobian(sys::DiffEqSystem)
    isempty(sys.jac[]) || return sys.jac[]  # use cached Jacobian, if possible
    rhs = [eq.rhs for eq in sys.eqs]

    jac = expand_derivatives.(calculate_jacobian(rhs, sys.dvs))
    sys.jac[] = jac  # cache Jacobian
    return jac
end

function generate_jacobian(sys::DiffEqSystem; version::FunctionVersion = ArrayFunction)
    jac = calculate_jacobian(sys)
    return build_function(jac, sys.dvs, sys.ps, (sys.iv.name,); version = version)
end

function generate_function(sys::DiffEqSystem; version::FunctionVersion = ArrayFunction)
    rhss = [eq.rhs for eq ∈ sys.eqs]
    return build_function(rhss, sys.dvs, sys.ps, (sys.iv.name,); version = version)
end

function generate_ode_iW(sys::DiffEqSystem, simplify=true; version::FunctionVersion = ArrayFunction)
    jac = calculate_jacobian(sys)

    gam = Variable(:gam; known = true)

    W = LinearAlgebra.I - gam*jac
    W = SMatrix{size(W,1),size(W,2)}(W)
    iW = inv(W)

    if simplify
        iW = simplify_constants.(iW)
    end

    W = inv(LinearAlgebra.I/gam - jac)
    W = SMatrix{size(W,1),size(W,2)}(W)
    iW_t = inv(W)
    if simplify
        iW_t = simplify_constants.(iW_t)
    end

    vs, ps = sys.dvs, sys.ps
    iW_func   = build_function(iW  , vs, ps, (:gam,:t); version = version)
    iW_t_func = build_function(iW_t, vs, ps, (:gam,:t); version = version)

    return (iW_func, iW_t_func)
end

function DiffEqBase.ODEFunction(sys::DiffEqSystem; version::FunctionVersion = ArrayFunction)
    expr = generate_function(sys; version = version)
    if version === ArrayFunction
        ODEFunction{true}(eval(expr))
    elseif version === SArrayFunction
        ODEFunction{false}(eval(expr))
    end
end


text:
# ----------------------------------------------------
# --- minimalTransceiver.jl
# ----------------------------------------------------
# This file is intented to give an example on how we can use tree based network architecture with SDR
# This file has to be run on a SDR based SoC that communicated with  a remote PC
# ---------
# |
# For instance, a USRP e310 with a julia session runs minimalTransceiver.jl
# On the PC side, the backend SDROverNetworks can be used to recover data from the E310
@info "Julia-based minimal transceiver";

using Distributed
module E310

# ----------------------------------------------------
# --- Module loading
# ----------------------------------------------------
# --- Module dependency
using UHDBindings
using Distributed
using ZMQ
@everywhere using Sockets

mutable struct Configuration
	carrierFreq::Float64;
	samplingRate::Float64;
	gain::Union{Float64,Int};
	Antenna::String;
	nbSamples::Int;
end
mutable struct MD
	error::Cint;
	timeStamp::Timestamp;
end

# Setting max piority to avoid CPU congestion
function setMaxPriority();
	pid = getpid();
	run(`renice -n -20 -p $pid`);
	run(`chrt -p 99 $pid`)
end

# ----------------------------------------------------
# --- Main call
# ----------------------------------------------------
function main(mode,carrierFreq, samplingRate, gain, nbSamples)
	# --- Setting a very first configuration
	radio = openUHD(carrierFreq, samplingRate, gain);
	# --- Configuration socket
	rtcSocket   = ZMQ.Socket(REP);
	bind(rtcSocket, "tcp://*:5555");
	# --- RTT socket for Tx
	rttSocket   = ZMQ.Socket(REQ);
	bind(rttSocket, "tcp://*:9999");
	# --- Socket for broadcast Rx
	brSocket = ZMQ.Socket(PUB);
	bind(brSocket, "tcp://*:1111");
	# --- Get samples
	sig		  = zeros(Complex{Cfloat}, radio.rx.packetSize);
	cnt		  = 0;
	# --- Mode used
	# mode    = :rx;
	# --- Processing
	try
		print(radio);
		# --- Second order loop setup
		flag      = false;
		# --- Interruption to update radio config
		@async begin
			while (true)
				# --- We wait in this @async for a reception
				receiver = ZMQ.recv(rtcSocket);
				# --- Here, we have receive something
				# Raise a flag because something happens
				flag = true;
				# we create an evaluation here
				res = Meta.parse(String(receiver));
				# and we update the radio and get back the desired feeback level
                (requestConfig, requestMD, mode, buffer,updateBuffer) = updateUHD!(radio, res,mode);
				if requestConfig
					# --- Sending the  configuration to Host
					sendConfig(rtcSocket, radio.rx, nbSamples);
				end
				if requestMD
					sendMD(rtcSocket, radio.rx);
				end
                if updateBuffer
                    # --- Replace sig by obtained buffer
                    sig = buffer;
                end
			end
			print(radio)
		end
        # if  mode == :txbuffer
            # --- Send data to radio
            # UHDBindings.send(radio, sig,true);
        # else
        while (true)
            if mode == :rx
                # --- Direct call to avoid allocation
                recv!(sig, radio);
                # --- To UDP socket
                ZMQ.send(brSocket, sig)
                yield();
            elseif mode == :tx
                # --- We now transmit data !
                # Wait for RTT from host
                ZMQ.send(rttSocket,0x01);
                # --- Get the data
                sig = convert.(Complex{Cfloat},ZMQ.recv(rttSocket));
                # --- Send data to radio
                UHDBindings.send(radio, sig,false);
                yield();
            elseif mode == :txbuffer
                # --- Send data to radio using cyclic mode
                nbE = UHDBindings.send(radio, sig,false);
                (nbE == 0) && (break);
                yield();
            end
        end
	catch exception;
		# --- Release USRP
		show(exception);
    end
    # --- Close UHD
    close(radio);
    # --- Close sockets
    close(rtcSocket);
    close(rttSocket);
    close(brSocket);
    return sig;
end

# --- To effectively update the radio config
function updateUHD!(radio, res,mode)
	# --- Default output
	requestConfig   = true;
	requestMD       = false;
    updateBuffer = false;
	# --- We create the dictionnary entry to update the radio config
	D = eval(res);
	# --- Apply the changes
    buffer = zeros(Complex{Cfloat},1024);
	for key in keys(D)
		elem = D[key];
		if key == :requestMD
			# --- We only ask for MD and not config
			requestMD = true;
			requestConfig  = false;
		elseif key == :requestConfig
			# --- Nothing to do
		elseif key == :updateCarrierFreq
			# --- Update the carrier freq
			updateCarrierFreq!(radio, elem);
		elseif key == :updateSamplingRate
			# --- Update sampling frequency
			updateSamplingRate!(radio, elem);
		elseif key == :updateGain
			# --- Update Gain
            println("$elem");
			updateGain!(radio, elem);
        elseif key == :buffer
            buffer = elem;
            updateBuffer = true;
		elseif key == :mode
			# --- We change mode
			mode = elem;
			requestConfig = false;
			requestMD = true;
			@info "Change mode to $mode";
		else
			@warn "Unknown Host order. Ask to update $key field with value $elem which is unknwown"
		end
	end
    print(radio)
	return (requestConfig, requestMD, mode, buffer, updateBuffer);
end

function sendConfig(rtcSocket, radio, nbSamples)
	# --- get Configuration from radio
	config = (radio.carrierFreq, radio.samplingRate, radio.gain, radio.antenna, nbSamples);
	# --- Send config
	strF = "$(config)";
	ZMQ.send(rtcSocket, strF);
end
function sendMD(rtcSocket, radio)
	md = (getTimestamp(radio)..., Cint(getError(radio)));
	# --- Send config
	strF = "$(md)";
	ZMQ.send(rtcSocket, strF);
end

end

tx()       = E310.main(:tx,868e6,4e6,10,(512 + 36) * 2 * 32);
txBuffer() = E310.main(:txbuffer,868e6,4e6,10,(512 + 36) * 2 * 32);
rx()       = E310.main(:rx,868e6,4e6,10,(512 + 36) * 2 * 32);

# call main function
# E310.main(:rx,868e6,4e6,10,(512 + 36) * 2 * 32);
# E310.main(868e6,4e6,10,32768);


text:
using Distributed

using CSV
using DataFrames
@everywhere using FileIO
using JLD2
using ProgressMeter
using Statistics

@assert length(ARGS) >= 1 "no experiment path given"

prefix_path = ARGS[1]
num_trajectories = length(ARGS)==2 ? parse(Int, ARGS[2]) : 1000

df = CSV.read(joinpath(prefix_path, "parameters_map.csv"), DataFrame)
output_dirs = map( x-> joinpath(prefix_path, x, "output"), df.path)
nans = fill(NaN, num_trajectories)

results = @showprogress pmap(output_dirs) do path
  try
    load(joinpath(path, "summary.jld2"), "last_infections", "num_infections", "peak_daily_detections", "peak_daily_infections")
  catch e
    @warn "could not load summary" path
    nans, nans, nans, nans
  end
end

last_infections = getindex.(results, 1)
num_infections = getindex.(results, 2)
peak_daily_detections = getindex.(results, 3)
peak_daily_infections = getindex.(results, 4)

df.mean_last_infection_time = last_infections .|> mean
df.mean_num_infected = num_infections .|> mean
df.peak_daily_detections = peak_daily_detections .|> mean
df.peak_daily_infections = peak_daily_infections .|> mean

save(joinpath(prefix_path,"summary.jld2"),
  "summary", df,
  "last_infections", last_infections,
  "num_infections", num_infections,
  "peak_daily_detections", peak_daily_detections,
  "peak_daily_infections", peak_daily_infections,
  compress=true
)

save(joinpath(prefix_path,"summary.csv"), df)


text:
struct Dynamics{T} <: Constraint
    val
    jac
    hess
    ny::Int
    nx::Int
    nu::Int
    nw::Int
    nj::Int
    nh::Int
    sp_jac::Vector{Vector{Int}}
    sp_hess::Vector{Vector{Int}}
    val_cache::Vector{T}
    jac_cache::Vector{T}
    hess_cache::Vector{T}
end

function Dynamics(f::Function, ny::Int, nx::Int, nu::Int; nw::Int=0, eval_hess=false)
    #TODO: option to load/save methods
    @variables y[1:ny], x[1:nx], u[1:nu], w[1:nw]
    val = f(y, x, u, w)
    jac = Symbolics.sparsejacobian(val, [x; u; y]);
    val_func = eval(Symbolics.build_function(val, y, x, u, w)[2]);
    jac_func = eval(Symbolics.build_function(jac.nzval, y, x, u, w)[2]);
    nj = length(jac.nzval)
    sp_jac = [findnz(jac)[1:2]...]
    if eval_hess
        @variables λ[1:ny]
        lag_con = dot(λ, val)
        hess = Symbolics.sparsehessian(lag_con, [x; u; y])
        hess_func = eval(Symbolics.build_function(hess.nzval, y, x, u, w, λ)[2])
        sp_hess = [findnz(hess)[1:2]...]
        nh = length(hess.nzval)
    else
        hess_func = Expr(:null)
        sp_hess = [Int[]]
        nh = 0
    end

    return Dynamics(val_func, jac_func, hess_func, ny, nx, nu, nw, nj, nh,
        sp_jac, sp_hess, zeros(ny), zeros(nj), zeros(nh))
end

function eval_con!(c, idx, cons::Vector{Dynamics{T}}, x, u, w) where T
    for (t, con) in enumerate(cons)
        con.val(con.val_cache, x[t+1], x[t], u[t], w[t])
        @views c[idx[t]] .= con.val_cache
        fill!(con.val_cache, 0.0) # TODO: confirm this is necessary
    end
end

function eval_jac!(j, idx, cons::Vector{Dynamics{T}}, x, u, w) where T
    for (t, con) in enumerate(cons)
        con.jac(con.jac_cache, x[t+1], x[t], u[t], w[t])
        @views j[idx[t]] .= con.jac_cache
        fill!(con.jac_cache, 0.0) # TODO: confirm this is necessary
    end
end

function eval_hess_lag!(h, idx, cons::Vector{Dynamics{T}}, x, u, w, λ) where T
    for (t, con) in enumerate(cons)
        con.hess(con.hess_cache, x[t+1], x[t], u[t], w[t], λ[t])
        @views h[idx[t]] .+= con.hess_cache
        fill!(con.hess_cache, 0.0) # TODO: confirm this is necessary
    end
end

function sparsity_jacobian(cons::Vector{Dynamics{T}}, nx::Vector{Int}, nu::Vector{Int}; row_shift=0) where T
    row = Int[]
    col = Int[]
    for (t, con) in enumerate(cons)
        col_shift = (t > 1 ? (sum(nx[1:t-1]) + sum(nu[1:t-1])) : 0)
        push!(row, (con.sp_jac[1] .+ row_shift)...)
        push!(col, (con.sp_jac[2] .+ col_shift)...)
        row_shift += con.ny
    end
    return collect(zip(row, col))
end

function sparsity_hessian(cons::Vector{Dynamics{T}}, nx::Vector{Int}, nu::Vector{Int}) where T
    row = Int[]
    col = Int[]
    for (t, con) in enumerate(cons)
        if !isempty(con.sp_hess[1])
            shift = (t > 1 ? (sum(nx[1:t-1]) + sum(nu[1:t-1])) : 0)
            push!(row, (con.sp_hess[1] .+ shift)...)
            push!(col, (con.sp_hess[2] .+ shift)...)
        end
    end
    return collect(zip(row, col))
end

num_con(cons::Vector{Dynamics{T}}) where T = sum([con.ny for con in cons])
num_xuy(cons::Vector{Dynamics{T}}) where T = sum([con.nx + con.nu for con in cons]) + cons[end].ny
num_jac(cons::Vector{Dynamics{T}}) where T = sum([con.nj for con in cons])

function constraint_indices(cons::Vector{Dynamics{T}}; shift=0) where T
    [collect(shift + (t > 1 ? sum([cons[s].ny for s = 1:(t-1)]) : 0) .+ (1:cons[t].ny)) for t = 1:length(cons)]
end

function jacobian_indices(cons::Vector{Dynamics{T}}; shift=0) where T
    [collect(shift + (t > 1 ? sum([cons[s].nj for s = 1:(t-1)]) : 0) .+ (1:cons[t].nj)) for t = 1:length(cons)]
end

function hessian_indices(cons::Vector{Dynamics{T}}, key::Vector{Tuple{Int,Int}}, nx::Vector{Int}, nu::Vector{Int}) where T
    idx = Vector{Int}[]
    for (t, con) in enumerate(cons)
        if !isempty(con.sp_hess[1])
            shift = (t > 1 ? (sum(nx[1:t-1]) + sum(nu[1:t-1])) : 0)
            row = collect(con.sp_hess[1] .+ shift)
            col = collect(con.sp_hess[2] .+ shift)
            rc = collect(zip(row, col))
            push!(idx, [findfirst(x -> x == i, key) for i in rc])
        end
    end
    return idx
end

function x_indices(cons::Vector{Dynamics{T}}) where T
    [[collect((t > 1 ? sum([cons[s].nx + cons[s].nu for s = 1:(t-1)]) : 0) .+ (1:cons[t].nx)) for t = 1:length(cons)]...,
        collect(sum([cons[s].nx + cons[s].nu for s = 1:length(cons)]) .+ (1:cons[end].ny))]
end

function u_indices(cons::Vector{Dynamics{T}}) where T
    [collect((t > 1 ? sum([cons[s].nx + cons[s].nu for s = 1:(t-1)]) : 0) + cons[t].nx .+ (1:cons[t].nu)) for t = 1:length(cons)]
end

function xu_indices(cons::Vector{Dynamics{T}}) where T
    [[collect((t > 1 ? sum([cons[s].nx + cons[s].nu for s = 1:(t-1)]) : 0) .+ (1:(+ cons[t].nx + cons[t].nu))) for t = 1:length(cons)]...,
        collect(sum([cons[s].nx + cons[s].nu for s = 1:length(cons)]) .+ (1:cons[end].ny))]
end

function xuy_indices(cons::Vector{Dynamics{T}}) where T
    [collect((t > 1 ? sum([cons[s].nx + cons[s].nu for s = 1:(t-1)]) : 0) .+ (1:(+ cons[t].nx + cons[t].nu + cons[t].ny))) for t = 1:length(cons)]
end

struct DynamicsIndices
    x::Vector{Vector{Int}}
    u::Vector{Vector{Int}}
    xu::Vector{Vector{Int}}
    xuy::Vector{Vector{Int}}
end

struct DynamicsDimensions
    x::Vector{Int}
    u::Vector{Int}
    w::Vector{Int}
end

function DynamicsDimensions(cons::Vector{Dynamics{T}};
    w_dim=[0 for t = 1:(length(cons) + 1)]) where T
    x_dim = [[con.nx for con in cons]..., cons[end].ny]
    u_dim = [con.nu for con in cons]
    return DynamicsDimensions(x_dim, u_dim, w_dim)
end

struct DynamicsModel{T}
    dyn::Vector{Dynamics{T}}
    idx::DynamicsIndices
    dim::DynamicsDimensions
end

function DynamicsModel(cons::Vector{Dynamics{T}};
    w_dim=[0 for t = 1:(length(cons) + 1)]) where T
    idx = DynamicsIndices(
            x_indices(cons),
            u_indices(cons),
            xu_indices(cons),
            xuy_indices(cons))
    dim = DynamicsDimensions(cons, w_dim=w_dim)
    DynamicsModel(cons, idx, dim)
end


text:
using Espresso: simplify, matchex, rewrite, set_default_placeholders
using Match: @match

set_default_placeholders(Set([:x, :y, :a, :b, :n]))

# Get the integration variable
integrationvar(diff::Symbol) = diff
function integrationvar(diff::Expr)
    var = nothing
    f(x) = nothing
    f(x::Symbol) = (var = x)
    walk(f, diff)

    @assert var ≠ nothing
    var
end

function evalIntegral(int::Integral)
    antiderivative = integrate(int.expr, int.diff)

    if int.bounds == nothing
        return antiderivative
    end

    var = integrationvar(int.diff)
    result = @eval Context begin
        F($var) = $antiderivative

        a, b = eval($int.bounds.a), eval($int.bounds.b)
        ε = eps(Float64)
        if !isfinite(F(a))
            a = (a < b) ? (a + ε) : (a - ε)
        end
        if !isfinite(F(b))
            b = (b > a) ? (b - ε) : (b + ε)
        end
        F(b) - F(a)
    end

    if !isfinite(result)
        error("function $(int.expr) cannot be integrated at ($(int.bounds.a), $(int.bounds.b)).")
    end

    result
end

function integrate(n, x)
    @match n begin
        0 => 0
        1 => x
        n::Number => Expr(:call, :*, n, x)
        _ => tableIntegral(n, x)
    end
end

function integrate(ex::Expr, x)
    if ex.head ≠ :call
        ArgumentError("ex must be a call expression") |> throw
    end

    op = ex.args[1]
    @match op begin
        :+ || :- => Expr(:call, op, [integrate(a, x) for a ∈ ex.args[2:end]]...)
        :* || :/ =>
            @match consts(ex) begin
                (0, _)  => 0
                (1, ex) => tableIntegral(ex, x)
                (n, ex) => Expr(:call, :*, n, tableIntegral(ex, x))
            end
        _ => tableIntegral(ex, x)
    end
end

function tableIntegral(ex::Value, x)
    int = :($(ex) * d($(x)))

    for rule ∈ integralTable
        pat = Expr(:call, :*, rule.left, :(d(x)))
        match = matchex(pat, int; allow_ex = true, exact = false)

        if match ≠ nothing && rule.condition(match)
            return rewrite(ex, rule.left, rule.right)
        end
    end

    nothing
end


text:
using CSV
using DataFrames

export personality2scores

ans_mapping = Dict("HO" => 1, "O" => 2, "N" => 3, "E" => 4, "HE" => 5)
reverse(ans::Number)::Number = 6 - ans
reversed_questions = [
    # Page 113
    1, 61, 121, 181,
    36, 96, 156,
    11, 71,
    46, 106, 166,
    21, 81, 141, 231,
    56, 116, 176, 206, 236,
    # Page 114
    32, 92,
    7, 67, 187,
    42, 102, 162, 222,
    17, 77, 137,
    52, 112,
    # Page 115
    33, 93, 153, 183,
    8, 68, 128,
    43, 103, 163,
    18, 78, 138, 198, 228,
    53, 113, 173,
    # Page 116
    4, 64, 124,
    39, 99, 159, 189, 219,
    14, 74, 134,
    49, 109, 169, 199, 229,
    24, 84, 144, 234,
    59, 119,
    # Page 117
    95, 155,
    10, 70, 130, 190, 220,
    45, 105,
    20, 80, 140,
    55, 115, 205,
    30, 90, 150,
]
is_reversed(question::Number)::Bool = question in reversed_questions

function ans2num(ans::AbstractString, question::Number)::Number
    ans = string(ans)
    ans = ans_mapping[ans]
    return is_reversed(question) ? reverse(ans) : ans
end

function ans2num(ans::AbstractString, question::AbstractString)
    return ans2num(string(ans), parse(Int64, question[2:end]))
end

ans2num(ans::Missing, question::AbstractString) = ans2num("N", question)
ans2num(ans::AbstractString, question::Symbol) = ans2num(string(ans), string(question))

function personality2digits(lima::DataFrame)::DataFrame
    lima_digits = select(lima, [:id, :completed_at])

    for question in names(lima)[3:end]
        digits = map(ans -> ans2num(ans, question), lima[!, question])
        lima_digits = hcat(lima_digits, DataFrame(question => digits))
    end
    return lima_digits
end

facets = Dict(
    "N1" => [1, 31, 61, 91, 121, 151, 181, 211],
    "N2" => [6, 36, 66, 96, 126, 156, 186, 216],
    "N3" => [11, 41, 71, 101, 131, 161, 191, 221],
    "N4" => [16, 46, 76, 106, 136, 166, 196, 226],
    "N5" => [21, 51, 81, 111, 141, 171, 201, 231],
    "N6" => [26, 56, 86, 116, 146, 176, 206, 236],
    "E1" => [2, 32, 62, 92, 122, 152, 182, 212],
    "E2" => [7,37, 67, 97, 127, 157, 187, 217],
    "E3" => [12, 42, 72, 102, 132, 162, 192, 222],
    "E4" => [17, 47, 77, 107, 137, 167, 197, 227],
    "E5" => [22, 52, 82, 112, 142, 172, 202, 232],
    "E6" => [27, 57, 87, 117, 147, 177, 207, 237],
    "O1" => [3, 33, 63, 93, 123, 153, 183, 213],
    "O2" => [8, 38, 68, 98, 128, 158, 188, 218],
    "O3" => [13, 43, 73, 103, 133, 163, 193, 223],
    "O4" => [18, 48, 78, 108, 138, 168, 198, 228],
    "O5" => [23, 53, 83, 113, 143, 173, 203, 233],
    "O6" => [28, 58, 88, 118, 148, 178, 208, 238],
    "A1" => [4, 34, 64, 94, 124, 154, 184, 214],
    "A2" => [9, 39, 69, 99, 129, 159, 189, 219],
    "A3" => [14, 44, 74, 104, 134, 164, 194, 224],
    "A4" => [19, 49, 79, 109, 139, 169, 199, 229],
    "A5" => [24, 54, 84, 114, 144, 174, 204, 234],
    "A6" => [29, 59, 89, 119, 149, 179, 209, 239],
    "C1" => [5, 35, 65, 95, 125, 155, 185, 215],
    "C2" => [10, 40, 70, 100, 130, 160, 190, 220],
    "C3" => [15, 45, 75, 105, 135, 165, 195, 225],
    "C4" => [20, 50, 80, 110, 140, 170, 200, 230],
    "C5" => [25, 55, 85, 115, 145, 175, 205, 235],
    "C6" => [30, 60, 90, 120, 150, 180, 210, 240]
)

function facet(question::Number)::String
    for key in keys(facets)
        if question in facets[key]
            return string(key)
        end
    end
    return "question not listed in facets"
end

domain(question::Number)::Char = facet(question)[1]
answer(lima::DataFrame, row::Number, question::Symbol) = lima[row, question]
answer(lima::DataFrame, row::Number, question::Number) = answer(lima, row, Symbol("v" * string(question)))

function score_everyone(lima_digits::DataFrame, facet::String)
    answers(row::Number) = map(question -> answer(lima_digits, row, question), facets[facet])
    map(row -> sum(answers(row)), 1:nrow(lima_digits))
end

function personality2facets(lima_digits::DataFrame)::DataFrame
    lima_facets = select(lima_digits, [:id, :completed_at])

    for facet in sort(collect(keys(facets)))
        df = DataFrame(Symbol(facet) => score_everyone(lima_digits, facet))
        lima_facets = hcat(lima_facets, df)
    end
    return lima_facets
end

function personality2scores(lima::DataFrame)::DataFrame
    lima_digits = personality2digits(lima)
    lima_facets = personality2facets(lima_digits)
    lima_scores = copy(lima_facets)

    answers(row::Number, facet::String) = map(question -> answer(lima_digits, row, question), facets[facet])
    score_everyone(facet::String) = map(row -> sum(answers(row, facet)), 1:nrow(lima_digits))

    domains = Set(map(key -> string(key)[1], collect(keys(facets))))
    get_facets(domain::Char)::Array{String,1} = filter(key -> string(key)[1] == domain, collect(keys(facets)))

    facet_scores(domain::Char) = map(score_everyone, get_facets(domain))
    score_everyone(domain::Char) = sum(facet_scores(domain))

    for domain in sort(collect(domains))
        df = DataFrame(Symbol(domain) => score_everyone(domain))
        lima_scores = hcat(lima_scores, df)
    end

    return lima_scores
end


text:
using Profile
using BenchmarkTools

"""
    build_sigma!(sigma_vector::CompressedTuckerState, ci_vector::CompressedTuckerState, cluster_ops, clustered_ham)
"""
function build_sigma_serial!(sigma_vector::CompressedTuckerState, ci_vector::CompressedTuckerState, cluster_ops, clustered_ham; nbody=4, cache=false)
    #={{{=#

    for (fock_bra, configs_bra) in sigma_vector
        for (fock_ket, configs_ket) in ci_vector
            fock_trans = fock_bra - fock_ket

            # check if transition is connected by H
            haskey(clustered_ham, fock_trans) == true || continue

            for (config_bra, coeff_bra) in configs_bra
                for (config_ket, coeff_ket) in configs_ket

                    for term in clustered_ham[fock_trans]

                        length(term.clusters) <= nbody || continue

                        check_term(term, fock_bra, config_bra, fock_ket, config_ket) || continue

                        coeff_bra.core .= form_sigma_block!(term, cluster_ops, fock_bra, config_bra,
                                                            fock_ket, config_ket,
                                                            coeff_bra, coeff_ket,
                                                            cache=cache)

                    end
                end
            end
        end
    end
    return
    #=}}}=#
end

"""
    build_sigma!(sigma_vector::CompressedTuckerState, ci_vector::CompressedTuckerState, cluster_ops, clustered_ham)
"""
function cache_hamiltonian_old(sigma_vector::CompressedTuckerState, ci_vector::CompressedTuckerState, cluster_ops, clustered_ham; nbody=4)
    #={{{=#

    println(" Cache hamiltonian terms")

    for (fock_bra, configs_bra) in sigma_vector
        for (fock_ket, configs_ket) in ci_vector
            fock_trans = fock_bra - fock_ket

            # check if transition is connected by H
            haskey(clustered_ham, fock_trans) == true || continue

            for (config_bra, coeff_bra) in configs_bra
                for (config_ket, coeff_ket) in configs_ket

                    for term in clustered_ham[fock_trans]

                        length(term.clusters) <= nbody || continue

                        check_term(term, fock_bra, config_bra, fock_ket, config_ket) || continue

                        cache_key = OperatorConfig((fock_bra, fock_ket, config_bra, config_ket))
                        term.cache[cache_key] = build_dense_H_term(term, cluster_ops,
                                                                   fock_bra, config_bra, coeff_bra,
                                                                   fock_ket, config_ket, coeff_ket)

                    end
                end
            end
        end
    end
    return
    #=}}}=#
end

function cache_hamiltonian(bra::CompressedTuckerState, ket::CompressedTuckerState, cluster_ops, clustered_ham; nbody=4)
    println(" Cache hamiltonian terms")

    keys_to_loop = [keys(clustered_ham.trans)...]
    println(" Number of threaded jobs:", length(keys_to_loop))
    Threads.@threads for ftrans in keys_to_loop
        terms = clustered_ham[ftrans]
        for term in terms

            length(term.clusters) <= nbody || continue

            for (fock_ket, configs_ket) in ket
                #fock_bra = [fock_ket.config...]
                #for (cii,ci) in enumerate(term.clusters)
                #    fock_bra[ci.idx] = (fock_ket[ci.idx][1] + ftrans[cii][1], fock_ket[ci.idx][2] + ftrans[cii][2])
                #end

                fock_bra = ftrans + fock_ket
                #fock_bra = FockConfig(fock_bra)

                haskey(bra.data, fock_bra) == true || continue

                for (config_ket, tuck_ket) in configs_ket
                    for (config_bra, tuck_bra) in bra[fock_bra]

                        check_term(term, fock_bra, config_bra, fock_ket, config_ket) || continue

                        cache_key = OperatorConfig((fock_bra, fock_ket, config_bra, config_ket))

                        term.cache[cache_key] = build_dense_H_term(term, cluster_ops,
                                                                   fock_bra, config_bra, tuck_bra,
                                                                   fock_ket, config_ket, tuck_ket)
                    end
                end
            end
        end
    end
end

"""
    build_sigma_parallel!(sigma_vector::CompressedTuckerState, ci_vector::CompressedTuckerState, cluster_ops, clustered_ham)
"""
function build_sigma!(sigma_vector::CompressedTuckerState, ci_vector::CompressedTuckerState, cluster_ops, clustered_ham; nbody=4, cache=false)
    #={{{=#

    jobs = []
    output = [[] for i in 1:Threads.nthreads()]
    for (fock_bra, configs_bra) in sigma_vector
        for (config_bra, tuck_bra) in configs_bra
            push!(jobs, [fock_bra, config_bra])
        end
    end

    function do_job(job)

        fock_bra = job[1]
        config_bra = job[2]
        coeff_bra = sigma_vector[fock_bra][config_bra]

        for (fock_ket, configs_ket) in ci_vector
            fock_trans = fock_bra - fock_ket

            # check if transition is connected by H
            haskey(clustered_ham, fock_trans) == true || continue

            for (config_ket, coeff_ket) in configs_ket

                for term in clustered_ham[fock_trans]

                    length(term.clusters) <= nbody || continue

                    check_term(term, fock_bra, config_bra, fock_ket, config_ket) || continue

                    # these methods dispatched on type of term
                    #coeff_bra.core .= form_sigma_block!(term, cluster_ops, fock_bra, config_bra,
                    #                              fock_ket, config_ket,
                    #                              coeff_bra, coeff_ket,
                    #                              cache=cache)
                    out = form_sigma_block!(term, cluster_ops, fock_bra, config_bra,
                                                  fock_ket, config_ket,
                                                  coeff_bra, coeff_ket,
                                                  cache=cache)

                    push!(output[Threads.threadid()], (fock_bra, config_bra, out))

                end
            end
        end
    end

    Threads.@threads for job in jobs
        do_job(job)
    end

    flush(stdout)

    for tid in output
        for out in tid
            fock_bra = out[1]
            config_bra = out[2]
            core = out[3]
            sigma_vector[fock_bra][config_bra].core .+= core
        end
    end
    return
    #=}}}=#
end

#
# form_sigma_block computes the action of the term on a Tucker compressed state,
# projected into the space defined by bra. This is used to work with H within a subspace defined by a compression
#
#
function form_sigma_block!(term::C,
                            cluster_ops::Vector{ClusterOps},
                            fock_bra::FockConfig, bra::TuckerConfig,
                            fock_ket::FockConfig, ket::TuckerConfig,
                            coeffs_bra::Tucker{T,N}, coeffs_ket::Tucker{T,N};
                            cache=false ) where {T,N, C<:ClusteredTerm}
    #={{{=#
    check_term(term, fock_bra, bra, fock_ket, ket) || throw(Exception)
    #
    # determine sign from rearranging clusters if odd number of operators
    state_sign = compute_terms_state_sign(term, fock_ket)

    # todo: add in 2e integral tucker decomposition and compress gamma along 1st index first

    op = Array{Float64}[]
    cache_key = OperatorConfig((fock_bra, fock_ket, bra, ket))
    #if cache && haskey(term.cache, cache_key)
    if cache

        #
        # read the dense H term
        op = term.cache[cache_key]

    else

        #cache == false || println(" couldn't find:", cache_key)

        #
        # build the dense H term
        op = build_dense_H_term(term, cluster_ops, fock_bra, bra, coeffs_bra, fock_ket, ket, coeffs_ket)

        #if cache
        #    term.cache[cache_key] = op
        #end
    end

    return contract_dense_H_with_state(term, op, state_sign, coeffs_bra, coeffs_ket)
end
#=}}}=#

"""
    Contract integrals and ClusterOps to form dense 4-body Hamiltonian matrix (tensor) in Tucker basis
"""
function build_dense_H_term(term::ClusteredTerm1B, cluster_ops, fock_bra, bra, coeffs_bra::Tucker, fock_ket, ket, coeffs_ket::Tucker)
#={{{=#
    c1 = term.clusters[1]
    op = Array{Float64}[]

    op1 = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])]

    #
    # Get 1body operator and compress it using the cluster's Tucker factors
    op = coeffs_bra.factors[c1.idx]' * (op1[bra[c1.idx],ket[c1.idx]] * coeffs_ket.factors[c1.idx])

    return op
end
#=}}}=#
function build_dense_H_term(term::ClusteredTerm2B, cluster_ops, fock_bra, bra, coeffs_bra::Tucker, fock_ket, ket, coeffs_ket::Tucker)
#={{{=#
    c1 = term.clusters[1]
    c2 = term.clusters[2]
    op = Array{Float64}[]

    #
    # Compress Gammas using the cluster's Tucker factors
    # e.g.,
    #   Gamma(pqr, I, J) Ul(I,k) Ur(J,l) = Gamma(pqr, k, l) where k and l are compressed indices
    @views gamma1 = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])][:,bra[c1.idx],ket[c1.idx]]
    Ul = coeffs_bra.factors[c1.idx]
    Ur = coeffs_ket.factors[c1.idx]
    @tensor begin
        tmp[p,k,J] := Ul[I,k] * gamma1[p,I,J]
        g1[p,k,l] := Ur[J,l] * tmp[p,k,J]
    end
    #g1 = _compress_local_operator(gamma1, Ul, Ur)
    #g1 = @ncon([gamma1, U1, U2], [[-1,2,3], [2,-2], [3,-3]])

    @views gamma2 = cluster_ops[c2.idx][term.ops[2]][(fock_bra[c2.idx],fock_ket[c2.idx])][:,bra[c2.idx],ket[c2.idx]]
    Ul = coeffs_bra.factors[c2.idx]
    Ur = coeffs_ket.factors[c2.idx]
    @tensor begin
        tmp[p,k,J] := Ul[I,k] * gamma2[p,I,J]
        g2[p,k,l] := Ur[J,l] * tmp[p,k,J]
    end
    #g2 = @ncon([gamma2, U1, U2], [[-1,2,3], [2,-2], [3,-3]])
    #display(("g1/2", size(g1), size(g2)))

    @tensor begin
        op[q,J,I] := term.ints[p,q] * g1[p,I,J]
        op[J,L,I,K] := op[q,J,I] * g2[q,K,L]
    end

    return op
end
#=}}}=#
function build_dense_H_term(term::ClusteredTerm3B, cluster_ops, fock_bra, bra, coeffs_bra::Tucker, fock_ket, ket, coeffs_ket::Tucker)
#={{{=#
    c1 = term.clusters[1]
    c2 = term.clusters[2]
    c3 = term.clusters[3]
    op = Array{Float64}[]

    #
    # Compress Gammas using the cluster's Tucker factors
    # e.g.,
    #   Gamma(pqr, I, J) Ul(I,k) Ur(J,l) = Gamma(pqr, k, l) where k and l are compressed indices
    @views gamma1 = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])][:,bra[c1.idx],ket[c1.idx]]
    Ul = coeffs_bra.factors[c1.idx]
    Ur = coeffs_ket.factors[c1.idx]
    @tensor begin
        tmp[p,k,J] := Ul[I,k] * gamma1[p,I,J]
        g1[p,k,l] := Ur[J,l] * tmp[p,k,J]
    end

    @views gamma2 = cluster_ops[c2.idx][term.ops[2]][(fock_bra[c2.idx],fock_ket[c2.idx])][:,bra[c2.idx],ket[c2.idx]]
    Ul = coeffs_bra.factors[c2.idx]
    Ur = coeffs_ket.factors[c2.idx]
    @tensor begin
        tmp[p,k,J] := Ul[I,k] * gamma2[p,I,J]
        g2[p,k,l] := Ur[J,l] * tmp[p,k,J]
    end
    #display(("g1/2", size(g1), size(g2)))

    @views gamma3 = cluster_ops[c3.idx][term.ops[3]][(fock_bra[c3.idx],fock_ket[c3.idx])][:,bra[c3.idx],ket[c3.idx]]
    Ul = coeffs_bra.factors[c3.idx]
    Ur = coeffs_ket.factors[c3.idx]
    @tensor begin
        tmp[p,k,J] := Ul[I,k] * gamma3[p,I,J]
        g3[p,k,l] := Ur[J,l] * tmp[p,k,J]
    end

    #
    # Now contract into 3body term
    #
    # h(p,q) * g1(p,I,J) * g2(q,K,L) = op(J,L,I,K)
    @tensor begin
        op[q,r,I,J] := term.ints[p,q,r] * g1[p,I,J]
        op[r,I,J,K,L] := op[q,r,I,J] * g2[q,K,L]
        op[J,L,N,I,K,M] := op[r,I,J,K,L] * g3[r,M,N]
    end

    return op
end
#=}}}=#
function build_dense_H_term(term::ClusteredTerm4B, cluster_ops, fock_bra, bra, coeffs_bra::Tucker, fock_ket, ket, coeffs_ket::Tucker)
#={{{=#
    c1 = term.clusters[1]
    c2 = term.clusters[2]
    c3 = term.clusters[3]
    c4 = term.clusters[4]
    op = Array{Float64}[]

    #
    # Compress Gammas using the cluster's Tucker factors
    # e.g.,
    #   Gamma(pqr, I, J) Ul(I,k) Ur(J,l) = Gamma(pqr, k, l) where k and l are compressed indices
    @views gamma1 = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])][:,bra[c1.idx],ket[c1.idx]]
    Ul = coeffs_bra.factors[c1.idx]
    Ur = coeffs_ket.factors[c1.idx]
    @tensor begin
        tmp[p,k,J] := Ul[I,k] * gamma1[p,I,J]
        g1[p,k,l] := Ur[J,l] * tmp[p,k,J]
    end

    @views gamma2 = cluster_ops[c2.idx][term.ops[2]][(fock_bra[c2.idx],fock_ket[c2.idx])][:,bra[c2.idx],ket[c2.idx]]
    Ul = coeffs_bra.factors[c2.idx]
    Ur = coeffs_ket.factors[c2.idx]
    @tensor begin
        tmp[p,k,J] := Ul[I,k] * gamma2[p,I,J]
        g2[p,k,l] := Ur[J,l] * tmp[p,k,J]
    end
    #display(("g1/2", size(g1), size(g2)))

    @views gamma3 = cluster_ops[c3.idx][term.ops[3]][(fock_bra[c3.idx],fock_ket[c3.idx])][:,bra[c3.idx],ket[c3.idx]]
    Ul = coeffs_bra.factors[c3.idx]
    Ur = coeffs_ket.factors[c3.idx]
    @tensor begin
        tmp[p,k,J] := Ul[I,k] * gamma3[p,I,J]
        g3[p,k,l] := Ur[J,l] * tmp[p,k,J]
    end

    @views gamma4 = cluster_ops[c4.idx][term.ops[4]][(fock_bra[c4.idx],fock_ket[c4.idx])][:,bra[c4.idx],ket[c4.idx]]
    Ul = coeffs_bra.factors[c4.idx]
    Ur = coeffs_ket.factors[c4.idx]
    @tensor begin
        tmp[p,k,J] := Ul[I,k] * gamma4[p,I,J]
        g4[p,k,l] := Ur[J,l] * tmp[p,k,J]
    end

    #
    # Now contract into 4body term
    #
    # h(p,q) * g1(p,I,J) * g2(q,K,L) = op(J,L,I,K)
    @tensor begin
        op[q,r,s,J,I] := term.ints[p,q,r,s] * g1[p,I,J]
        op[r,s,J,L,I,K] := op[q,r,s,J,I] * g2[q,K,L]
        op[s,J,L,N,I,K,M] := op[r,s,J,L,I,K] * g3[r,M,N]
        op[J,L,N,P,I,K,M,O] := op[s,J,L,N,I,K,M] * g4[s,O,P]
    end
    return op
end
#=}}}=#

"""
    Contract  Hamiltonian matrix (tensor) in Tucker basis with trial vector (Tucker)
"""
function contract_dense_H_with_state(term::ClusteredTerm1B, op, state_sign, coeffs_bra::Tucker{T,N}, coeffs_ket::Tucker{T,N}) where {T,N}
#={{{=#
    c1 = term.clusters[1]

    n_clusters = N
    #
    # form overlaps - needed when TuckerConfigs aren't the same because each does their own compression and has
    # distinct Tucker factors
    overlaps = Dict{Int,Matrix{T}}()
    s = 1.0 # this is the product of scalar overlaps that don't need tensor contractions
    for ci in 1:n_clusters
        ci != c1.idx || continue

        S = coeffs_bra.factors[ci]' * coeffs_ket.factors[ci]

        # if overlap not just scalar, form and prepare for contraction
        # if it is a scalar, then just update the sign with the value
        if length(S) == 1
            s *= S[1]
        else
            overlaps[ci] = S
        end
    end

    indices = collect(1:n_clusters)
    indices[c1.idx] = 0
    perm,_ = bubble_sort(indices)

    coeffs_bra2 = copy(coeffs_bra.core)
    #coeffs_ket2 = copy(coeffs_ket.core)

    #
    # multiply by overlaps first if the bra side is smaller,
    # otherwise multiply by Hamiltonian term first
    coeffs_ket2 = transform_basis(coeffs_ket.core, overlaps, trans=true)

    #
    # Transpose to get contig data for blas (explicit copy?)
    coeffs_ket2 = permutedims(coeffs_ket2, perm)
    coeffs_bra2 = permutedims(coeffs_bra2, perm)

    #
    # Reshape for matrix multiply, shouldn't do any copies, right?
    dim1 = size(coeffs_ket2)
    coeffs_ket2 = reshape(coeffs_ket2, dim1[1], prod(dim1[2:end]))
    dim2 = size(coeffs_bra2)
    coeffs_bra2 = reshape(coeffs_bra2, dim2[1], prod(dim2[2:end]))

    #
    # Reshape Hamiltonian term operator
    # ... not needed for 1b term

    #
    # Multiply
    coeffs_bra2 .+= s .* (op * coeffs_ket2)

    # now untranspose
    perm,_ = bubble_sort(perm)
    coeffs_bra2 = reshape(coeffs_bra2, dim2)
    coeffs_bra2 = permutedims(coeffs_bra2,perm)

    #
    # multiply by overlaps now if the bra side is larger,
    if length(coeffs_bra2) >= length(coeffs_ket2)
        #coeffs_ket2 = transform_basis(coeffs_ket2, overlaps, trans=true)
    end

    return coeffs_bra2
end
#=}}}=#
function contract_dense_H_with_state(term::ClusteredTerm2B, op, state_sign, coeffs_bra::Tucker{T,N}, coeffs_ket::Tucker{T,N}) where {T,N}
#={{{=#
    c1 = term.clusters[1]
    c2 = term.clusters[2]

    n_clusters = N
    #
    # form overlaps - needed when TuckerConfigs aren't the same because each does their own compression and has
    # distinct Tucker factors
    overlaps = Dict{Int,Matrix{T}}()
    s = state_sign # this is the product of scalar overlaps that don't need tensor contractions
    for ci in 1:n_clusters
        ci != c1.idx || continue
        ci != c2.idx || continue

        S = coeffs_bra.factors[ci]' * coeffs_ket.factors[ci]

        # if overlap not just scalar, form and prepare for contraction
        if length(S) == 1
            s *= S[1]
        else
            overlaps[ci] = S
        end
    end

    indices = collect(1:n_clusters)
    indices[c1.idx] = 0
    indices[c2.idx] = 0
    perm,_ = bubble_sort(indices)

    coeffs_bra2 = copy(coeffs_bra.core)
    coeffs_ket2 = copy(coeffs_ket.core)

    #
    # multiply by overlaps first if the bra side is smaller,
    # otherwise multiply by Hamiltonian term first
    if length(coeffs_bra2) < length(coeffs_ket2)
        #coeffs_ket2 = transform_basis(coeffs_ket2, overlaps, trans=true)
    end
    coeffs_ket2 = transform_basis(coeffs_ket2, overlaps, trans=true)

    #
    # Transpose to get contig data for blas (explicit copy?)
    coeffs_ket2 = permutedims(coeffs_ket2, perm)
    coeffs_bra2 = permutedims(coeffs_bra2, perm)

    #
    # Reshape for matrix multiply, shouldn't do any copies, right?
    dim1 = size(coeffs_ket2)
    coeffs_ket2 = reshape(coeffs_ket2, dim1[1]*dim1[2], prod(dim1[3:end]))
    dim2 = size(coeffs_bra2)
    coeffs_bra2 = reshape(coeffs_bra2, dim2[1]*dim2[2], prod(dim2[3:end]))

    #
    # Reshape Hamiltonian term operator
    op = reshape(op, prod(size(op)[1:2]), prod(size(op)[3:4]))

    #display(term)
    ##display(overlaps)
    #display((size(coeffs_bra2), size(coeffs_bra.core)))
    #display(size(op'))
    #display((size(coeffs_ket2), size(coeffs_ket.core)))
    coeffs_bra2 .+= s .* (op' * coeffs_ket2)

    coeffs_bra2 = reshape(coeffs_bra2, dim2)

    # now untranspose
    perm,_ = bubble_sort(perm)
    coeffs_bra2 = permutedims(coeffs_bra2,perm)

    #
    # multiply by overlaps now if the bra side is larger,
    if length(coeffs_bra2) >= length(coeffs_ket2)
        #coeffs_ket2 = transform_basis(coeffs_ket2, overlaps, trans=true)
    end

    return coeffs_bra2
end
#=}}}=#
function contract_dense_H_with_state(term::ClusteredTerm3B, op, state_sign, coeffs_bra::Tucker{T,N}, coeffs_ket::Tucker{T,N}) where {T,N}
#={{{=#
    c1 = term.clusters[1]
    c2 = term.clusters[2]
    c3 = term.clusters[3]

    n_clusters = N
    #
    # form overlaps - needed when TuckerConfigs aren't the same because each does their own compression and has
    # distinct Tucker factors
    overlaps = Dict{Int,Matrix{T}}()
    s = state_sign # this is the product of scalar overlaps that don't need tensor contractions
    for ci in 1:n_clusters
        ci != c1.idx || continue
        ci != c2.idx || continue
        ci != c3.idx || continue

        S = coeffs_bra.factors[ci]' * coeffs_ket.factors[ci]

        # if overlap not just scalar, form and prepare for contraction
        if length(S) == 1
            s *= S[1]
        else
            overlaps[ci] = S
        end
    end

    indices = collect(1:n_clusters)
    indices[c1.idx] = 0
    indices[c2.idx] = 0
    indices[c3.idx] = 0
    perm,_ = bubble_sort(indices)

    coeffs_bra2 = copy(coeffs_bra.core)
    coeffs_ket2 = copy(coeffs_ket.core)

    #
    # multiply by overlaps first if the bra side is smaller,
    # otherwise multiply by Hamiltonian term first
    if length(coeffs_bra2) < length(coeffs_ket2)
        #coeffs_ket2 = transform_basis(coeffs_ket2, overlaps, trans=true)
    end
    coeffs_ket2 = transform_basis(coeffs_ket2, overlaps, trans=true)

    #
    # Transpose to get contig data for blas (explicit copy?)
    coeffs_ket2 = permutedims(coeffs_ket2, perm)
    coeffs_bra2 = permutedims(coeffs_bra2, perm)

    #
    # Reshape for matrix multiply, shouldn't do any copies, right?
    dim1 = size(coeffs_ket2)
    coeffs_ket2 = reshape(coeffs_ket2, dim1[1]*dim1[2]*dim1[3], prod(dim1[4:end]))
    dim2 = size(coeffs_bra2)
    coeffs_bra2 = reshape(coeffs_bra2, dim2[1]*dim2[2]*dim2[3], prod(dim2[4:end]))

    #
    # Reshape Hamiltonian term operator
    op = reshape(op, prod(size(op)[1:3]), prod(size(op)[4:6]))

    #
    # Multiply
    coeffs_bra2 .+= s .* (op' * coeffs_ket2)

    # now untranspose
    perm,_ = bubble_sort(perm)
    coeffs_bra2 = reshape(coeffs_bra2, dim2)
    coeffs_bra2 = permutedims(coeffs_bra2,perm)

    #
    # multiply by overlaps now if the bra side is larger,
    if length(coeffs_bra2) >= length(coeffs_ket2)
        #coeffs_ket2 = transform_basis(coeffs_ket2, overlaps, trans=true)
    end

    return coeffs_bra2
end
#=}}}=#
function contract_dense_H_with_state(term::ClusteredTerm4B, op, state_sign, coeffs_bra::Tucker{T,N}, coeffs_ket::Tucker{T,N}) where {T,N}
#={{{=#
    c1 = term.clusters[1]
    c2 = term.clusters[2]
    c3 = term.clusters[3]
    c4 = term.clusters[4]

    n_clusters = N

    #
    # form overlaps - needed when TuckerConfigs aren't the same because each does their own compression and has
    # distinct Tucker factors
    overlaps = Dict{Int,Matrix{T}}()
    s = state_sign # this is the product of scalar overlaps that don't need tensor contractions
    for ci in 1:n_clusters
        ci != c1.idx || continue
        ci != c2.idx || continue
        ci != c3.idx || continue
        ci != c4.idx || continue

        S = coeffs_bra.factors[ci]' * coeffs_ket.factors[ci]

        # if overlap not just scalar, form and prepare for contraction
        if length(S) == 1
            s *= S[1]
        else
            overlaps[ci] = S
        end
    end

    indices = collect(1:n_clusters)
    indices[c1.idx] = 0
    indices[c2.idx] = 0
    indices[c3.idx] = 0
    indices[c4.idx] = 0
    perm,_ = bubble_sort(indices)

    coeffs_bra2 = copy(coeffs_bra.core)
    coeffs_ket2 = copy(coeffs_ket.core)

    #
    # multiply by overlaps first if the bra side is smaller,
    # otherwise multiply by Hamiltonian term first
    if length(coeffs_bra2) < length(coeffs_ket2)
        #coeffs_ket2 = transform_basis(coeffs_ket2, overlaps, trans=true)
    end
    coeffs_ket2 = transform_basis(coeffs_ket2, overlaps, trans=true)

    #
    # Transpose to get contig data for blas (explicit copy?)
    coeffs_ket2 = permutedims(coeffs_ket2, perm)
    coeffs_bra2 = permutedims(coeffs_bra2, perm)

    #
    # Reshape for matrix multiply, shouldn't do any copies, right?
    dim1 = size(coeffs_ket2)
    coeffs_ket2 = reshape(coeffs_ket2, dim1[1]*dim1[2]*dim1[3]*dim1[4], prod(dim1[5:end]))
    dim2 = size(coeffs_bra2)
    coeffs_bra2 = reshape(coeffs_bra2, dim2[1]*dim2[2]*dim2[3]*dim2[4], prod(dim2[5:end]))

    #
    # Reshape Hamiltonian term operator
    op = reshape(op, prod(size(op)[1:4]), prod(size(op)[5:8]))

    #
    # Multiply
    coeffs_bra2 .+= s .* (op' * coeffs_ket2)

    # now untranspose
    perm,_ = bubble_sort(perm)
    coeffs_bra2 = reshape(coeffs_bra2, dim2)
    coeffs_bra2 = permutedims(coeffs_bra2,perm)

    #
    # multiply by overlaps now if the bra side is larger,
    if length(coeffs_bra2) >= length(coeffs_ket2)
        #coeffs_ket2 = transform_basis(coeffs_ket2, overlaps, trans=true)
    end
    return coeffs_bra2
end
#=}}}=#

function _compress_local_operator(gamma, Ul::Matrix{T}, Ur::Matrix{T}) where T
# this is way slower than @tensor
#={{{=#
    # gamma has 3 indices (orbital indices, cluster indices (left), cluster indices (right)

    #
    # out(i,jp) = gamma(p,I,J) Ul(I,i)
    out = Ul' * reshape(permutedims(gamma, [2,3,1]), size(gamma,2), size(gamma,3)*size(gamma,1))

    #
    # out(j,pi) = out(J,pi) Ur(J,j)
    out = Ur' * reshape(out', size(gamma,3), size(gamma,1)*size(Ul,2))

    # out(j,pi) -> out(p,i,j)
    return reshape(out', size(gamma,1), size(Ul,2), size(Ur,2))

#    # out(i,pJ) = gamma(I,pJ) U(I,i)
#    out = Ul' * unfold(gamma, 2)
#    # out(ip,J)
#    out = reshape(out, size(out,1) * size(gamma,1), size(gamma,3))
#
#    # out(ip,j) = gamma(ip,J) U(J,j)
#    out = out * Ur
#    # out(i,p,j)
#    out = reshape(out, size(Ul,2), size(gamma,1), size(Ur,2))
#
#    # out(p,i,j)
#    return permutedims(out, [2,1,3])
end
#=}}}=#

#
# form_sigma_block_expand computes the action of the term on a Tucker compressed state,
# NOT projected into the space defined by bra, but rather prescreen each operator defined by prescreen.
# This is used to find the first order interactiong space (FOIS) from a Tucker compressed state.
#

function form_sigma_block_expand(term::ClusteredTerm1B,
                            cluster_ops::Vector{ClusterOps},
                            fock_bra::FockConfig, bra::TuckerConfig,
                            fock_ket::FockConfig, ket::TuckerConfig,
                            coeffs_ket::Tucker{T,N};
                            max_number=nothing, prescreen=1e-4) where {T,N}
#={{{=#
    #display(term)
    c1 = term.clusters[1]
    n_clusters = length(bra)

    #
    # make sure active clusters are correct transitions
    fock_bra[c1.idx] == fock_ket[c1.idx] .+ term.delta[1] || throw(Exception)

    #
    # determine sign from rearranging clusters if odd number of operators
    state_sign = compute_terms_state_sign(term, fock_ket)

    #
    # Get 1body operator and compress it using the cluster's Tucker factors,
    # but since we are expanding the compression space
    # only compress the right hand side
    op1 = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])]
    op =  (op1[bra[c1.idx],ket[c1.idx]] * coeffs_ket.factors[c1.idx])

    tensors = Vector{Array{T}}()
    indices = Vector{Vector{Int16}}()
    state_indices = -collect(1:n_clusters)
    s = 1.0 # this is the product of scalar overlaps that don't need tensor contractions

    # if the compressed operator becomes a scalar, treat it as such
    if length(op) == 1
        s *= op[1]
    else
        op_indices = [-c1.idx, c1.idx]
        state_indices[c1.idx] = c1.idx
        push!(tensors, op)
        push!(indices, op_indices)
    end

    push!(tensors, coeffs_ket.core)
    push!(indices, state_indices)

    length(tensors) == length(indices) || error(" mismatch between operators and indices")

    bra_core = zeros(1,1)
    if length(tensors) == 1
        bra_core = coeffs_ket.core .* s
    else
        bra_core = @ncon(tensors, indices)
        bra_core .= bra_core .* s
    end

    # first decompose the already partially decomposed core tensor
    #
    # Vket ~ Aket x U1 x U2 x ...
    #
    # then we modified the compressed coefficients in the ket, Aket
    #
    # to get Abra, which we then compress again.
    #
    # The active cluster tucker factors become identity matrices
    #
    # Abra ~ Bbra x UU1 x UU2 x ....
    #
    #
    # e.g, V(IJK) = C(ijk) * U1(Ii) * U2(Jj) * U3(Kk)
    #
    # then C get's modified and furhter compressed
    #
    # V(IJK) = C(abc) * U1(ia) * U2(jb) * U3(kc) * U1(Ii) * U2(Jj) * U3(Kk)
    # V(IJK) = C(abc) * (U1(ia) * U1(Ii)) * (U2(jb) * U2(Jj)) * (U3(kc) * U3(Kk))
    # V(IJK) = C(abc) * U1(Ia)  * U2(Jb) * U3(Kc)
    #

    new_factors = [coeffs_ket.factors[i] for i in 1:N]
    new_factors[c1.idx] = Matrix(1.0I, size(bra_core,c1.idx), size(bra_core,c1.idx))
    return Tucker(bra_core, NTuple{N}(new_factors))

end
#=}}}=#

function form_sigma_block_expand(term::ClusteredTerm2B,
                            cluster_ops::Vector{ClusterOps},
                            fock_bra::FockConfig, bra::TuckerConfig,
                            fock_ket::FockConfig, ket::TuckerConfig,
                            coeffs_ket::Tucker{T,N};
                            max_number=nothing, prescreen=1e-4) where {T,N}
    #={{{=#
    #display(term)
    #display.((fock_bra, fock_ket))
    c1 = term.clusters[1]
    c2 = term.clusters[2]
    n_clusters = length(bra)
    #bound1 = norm(term.ints)*norm(coeffs_ket.core)
    ##bound1 *= norm(gamma1)*norm(gamma2)
    ##bound1 *= norm(coeffs_ket.factors[c1.idx])*norm(coeffs_ket.factors[c2.idx])
    #if bound1 < sqrt(prescreen)
    #    return FermiCG.Tucker(Array{Float64}(undef,[0 for i in 1:N]...), ntuple(i->Array{Float64}(undef,0,0), N))
    #end

    #
    # make sure active clusters are correct transitions
    fock_bra[c1.idx] == fock_ket[c1.idx] .+ term.delta[1] || throw(Exception)
    fock_bra[c2.idx] == fock_ket[c2.idx] .+ term.delta[2] || throw(Exception)

    #
    # determine sign from rearranging clusters if odd number of operators
    state_sign = compute_terms_state_sign(term, fock_ket)

    #
    # op[IK,JL] = <I|p'|J> h(pq) <K|q|L>

    # todo: add in 2e integral tucker decomposition and compress gamma along 1st index first

    #
    # Compress Gammas using the cluster's Tucker factors, but since we are expanding the compression space
    # only compress the right hand side
    # e.g.,
    #   Gamma(pqr, I, J) Ur(J,l) = Gamma(pqr, I, l) where k and l are compressed indices
    @views gamma1 = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])][:,bra[c1.idx],ket[c1.idx]]
    @views gamma2 = cluster_ops[c2.idx][term.ops[2]][(fock_bra[c2.idx],fock_ket[c2.idx])][:,bra[c2.idx],ket[c2.idx]]

    Ur = coeffs_ket.factors[c1.idx]
    @tensor begin
        g1[p,I,l] := Ur[J,l] * gamma1[p,I,J]
    end
    #@printf("Norm: %12.8f %12.8f %12.8f\n",norm(gamma1),norm(coeffs_ket.core),norm(gamma1)*norm(coeffs_ket.core))

    Ur = coeffs_ket.factors[c2.idx]
    @tensor begin
        g2[p,I,l] := Ur[J,l] * gamma2[p,I,J]
    end

    #
    # Decompose the local operators. Since gamma[p,I,l] has indices (small, large, small),
    # we only need at most p*l number of new vectors for the index we are searching over

    new_factor1 = Matrix(1.0I, size(g1,2), size(g1,2))
    new_factor2 = Matrix(1.0I, size(g2,2), size(g2,2))

    D = permutedims(g1, [2,1,3])
    F = svd(reshape(D, size(D,1), size(D,2)*size(D,3)))
    nkeep = 0
    for si in F.S
        if si > 1e-12
            nkeep += 1
        end
    end
    new_factor1 = F.U[:,1:nkeep]
    g1 = Diagonal(F.S[1:nkeep]) * F.Vt[1:nkeep,:]
    g1 = reshape(g1, size(g1,1), size(D,2), size(D,3))
    g1 = permutedims(g1, [2,1,3])

    D = permutedims(g2, [2,1,3])
    F = svd(reshape(D, size(D,1), size(D,2)*size(D,3)))
    nkeep = 0
    for si in F.S
        if si > 1e-12
            nkeep += 1
        end
    end
    new_factor2 = F.U[:,1:nkeep]
    g2 = Diagonal(F.S[1:nkeep]) * F.Vt[1:nkeep,:]
    g2 = reshape(g2, size(g2,1), size(D,2), size(D,3))
    g2 = permutedims(g2, [2,1,3])

    #
    # Now contract into 2body term
    #
    # h(p,q) * g1(p,I,J) * g2(q,K,L) = op(J,L,I,K)
    op = Array{Float64}[]
    @tensor begin
        op[q,J,I] := term.ints[p,q] * g1[p,I,J]
        op[J,L,I,K] := op[q,J,I] * g2[q,K,L]
    end

    #
    # form overlaps - needed when TuckerConfigs aren't the same because each does their own compression and has
    # distinct Tucker factors
    tensors = Vector{Array{T}}()
    indices = Vector{Vector{Int16}}()
    state_indices = -collect(1:n_clusters)
    s = state_sign # this is the product of scalar overlaps that don't need tensor contractions

    # if the compressed operator becomes a scalar, treat it as such
    if length(op) == 1
        s *= op[1]
    else
        op_indices = [c1.idx, c2.idx, -c1.idx, -c2.idx]
        state_indices[c1.idx] = c1.idx
        state_indices[c2.idx] = c2.idx
        push!(tensors, op)
        push!(indices, op_indices)
    end

    push!(tensors, coeffs_ket.core)
    push!(indices, state_indices)

    length(tensors) == length(indices) || error(" mismatch between operators and indices")

    bra_core = zeros(1,1)
    if length(tensors) == 1
        # this means that all the overlaps and the operator is a scalar
        bra_core = coeffs_ket.core .* s
    else
        #display.(("a", size(coeffs_bra), size(coeffs_ket), "sizes: ", size.(overlaps), indices))
        #display.(("a", size(coeffs_bra), size(coeffs_ket), "sizes: ", overlaps, indices))
        bra_core = @ncon(tensors, indices)
        bra_core .= bra_core .* s
    end

    # first decompose the already partially decomposed core tensor
    #
    # Vket ~ Aket x U1 x U2 x ...
    #
    # then we modified the compressed coefficients in the ket, Aket
    #
    # to get Abra, which we then compress again.
    #
    # The active cluster tucker factors become identity matrices
    #
    # Abra ~ Bbra x UU1 x UU2 x ....
    #
    #
    # e.g, V(IJK) = C(ijk) * U1(Ii) * U2(Jj) * U3(Kk)
    #
    # then C get's modified and furhter compressed
    #
    # V(IJK) = C(abc) * U1(ia) * U2(jb) * U3(kc) * U1(Ii) * U2(Jj) * U3(Kk)
    # V(IJK) = C(abc) * (U1(ia) * U1(Ii)) * (U2(jb) * U2(Jj)) * (U3(kc) * U3(Kk))
    # V(IJK) = C(abc) * U1(Ia)  * U2(Jb) * U3(Kc)
    #

    new_factors = [coeffs_ket.factors[i] for i in 1:N]
    #new_factors[c1.idx] = Matrix(1.0I, size(bra_core,c1.idx), size(bra_core,c1.idx))
    #new_factors[c2.idx] = Matrix(1.0I, size(bra_core,c2.idx), size(bra_core,c2.idx))
    new_factors[c1.idx] = new_factor1
    new_factors[c2.idx] = new_factor2
    return Tucker(bra_core, NTuple{N}(new_factors))
end
#=}}}=#

function form_sigma_block_expand(term::ClusteredTerm3B,
                            cluster_ops::Vector{ClusterOps},
                            fock_bra::FockConfig, bra::TuckerConfig,
                            fock_ket::FockConfig, ket::TuckerConfig,
                            coeffs_ket::Tucker{T,N};
                            max_number=nothing, prescreen=1e-4) where {T,N}
    #={{{=#
    #display(term)
    #display.((fock_bra, fock_ket))
    c1 = term.clusters[1]
    c2 = term.clusters[2]
    c3 = term.clusters[3]
    n_clusters = length(bra)

    #
    # make sure active clusters are correct transitions
    fock_bra[c1.idx] == fock_ket[c1.idx] .+ term.delta[1] || throw(Exception)
    fock_bra[c2.idx] == fock_ket[c2.idx] .+ term.delta[2] || throw(Exception)
    fock_bra[c3.idx] == fock_ket[c3.idx] .+ term.delta[3] || throw(Exception)

    #
    # determine sign from rearranging clusters if odd number of operators
    state_sign = compute_terms_state_sign(term, fock_ket)

    #
    # op[IK,JL] = <I|p'|J> h(pq) <K|q|L>
    # todo: add in 2e integral tucker decomposition and compress gamma along 1st index first

    #
    # Compress Gammas using the cluster's Tucker factors, but since we are expanding the compression space
    # only compress the right hand side
    # e.g.,
    #   Gamma(pqr, I, J) Ur(J,l) = Gamma(pqr, I, l) where k and l are compressed indices
    @views gamma1 = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])][:,bra[c1.idx],ket[c1.idx]]
    Ur = coeffs_ket.factors[c1.idx]
    @tensor begin
        g1[p,I,l] := Ur[J,l] * gamma1[p,I,J]
    end

    @views gamma2 = cluster_ops[c2.idx][term.ops[2]][(fock_bra[c2.idx],fock_ket[c2.idx])][:,bra[c2.idx],ket[c2.idx]]
    Ur = coeffs_ket.factors[c2.idx]
    @tensor begin
        g2[p,I,l] := Ur[J,l] * gamma2[p,I,J]
    end

    @views gamma3 = cluster_ops[c3.idx][term.ops[3]][(fock_bra[c3.idx],fock_ket[c3.idx])][:,bra[c3.idx],ket[c3.idx]]
    Ur = coeffs_ket.factors[c3.idx]
    @tensor begin
        g3[p,I,l] := Ur[J,l] * gamma3[p,I,J]
    end

    #
    # Decompose the local operators. Since gamma[p,I,l] has indices (small, large, small),
    # we only need at most p*l number of new vectors for the index we are searching over

    new_factor1 = Matrix(1.0I, size(g1,2), size(g1,2))
    new_factor2 = Matrix(1.0I, size(g2,2), size(g2,2))
    new_factor3 = Matrix(1.0I, size(g3,2), size(g3,2))

    D = permutedims(g1, [2,1,3])
    F = svd(reshape(D, size(D,1), size(D,2)*size(D,3)))
    nkeep = 0
    for si in F.S
        if si > prescreen
            nkeep += 1
        end
    end
    new_factor1 = F.U[:,1:nkeep]
    g1 = Diagonal(F.S[1:nkeep]) * F.Vt[1:nkeep,:]
    g1 = reshape(g1, size(g1,1), size(D,2), size(D,3))
    g1 = permutedims(g1, [2,1,3])

    D = permutedims(g2, [2,1,3])
    F = svd(reshape(D, size(D,1), size(D,2)*size(D,3)))
    nkeep = 0
    for si in F.S
        if si > prescreen
            nkeep += 1
        end
    end
    new_factor2 = F.U[:,1:nkeep]
    g2 = Diagonal(F.S[1:nkeep]) * F.Vt[1:nkeep,:]
    g2 = reshape(g2, size(g2,1), size(D,2), size(D,3))
    g2 = permutedims(g2, [2,1,3])

    D = permutedims(g3, [2,1,3])
    F = svd(reshape(D, size(D,1), size(D,2)*size(D,3)))
    nkeep = 0
    for si in F.S
        if si > prescreen
            nkeep += 1
        end
    end
    #
    # for now, let's just keep the full space, then maybe later start threshing
    new_factor3 = F.U[:,1:nkeep]
    g3 = Diagonal(F.S[1:nkeep]) * F.Vt[1:nkeep,:]
    g3 = reshape(g3, size(g3,1), size(D,2), size(D,3))
    g3 = permutedims(g3, [2,1,3])

    #
    # Now contract into 3body term
    #
    # h(p,q) * g1(p,I,J) * g2(q,K,L) = op(J,L,I,K)
    op = Array{Float64}[]
    @tensor begin
        op[q,r,I,J] := term.ints[p,q,r] * g1[p,I,J]
    end
    @tensor begin
        op[r,I,J,K,L] := op[q,r,I,J] * g2[q,K,L]
    end
    @tensor begin
        op[J,L,N,I,K,M] := op[r,I,J,K,L] * g3[r,M,N]
    end

    #
    # form overlaps - needed when TuckerConfigs aren't the same because each does their own compression and has
    # distinct Tucker factors
    tensors = Vector{Array{T}}()
    indices = Vector{Vector{Int16}}()
    state_indices = -collect(1:n_clusters)
    s = state_sign # this is the product of scalar overlaps that don't need tensor contractions

    # if the compressed operator becomes a scalar, treat it as such
    if length(op) == 1
        s *= op[1]
    else
        op_indices = [c1.idx, c2.idx, c3.idx, -c1.idx, -c2.idx, -c3.idx]
        state_indices[c1.idx] = c1.idx
        state_indices[c2.idx] = c2.idx
        state_indices[c3.idx] = c3.idx
        push!(tensors, op)
        push!(indices, op_indices)
    end

    push!(tensors, coeffs_ket.core)
    push!(indices, state_indices)

    length(tensors) == length(indices) || error(" mismatch between operators and indices")

    bra_core = zeros(1,1)
    if length(tensors) == 1
        # this means that all the overlaps and the operator is a scalar
        bra_core = coeffs_ket.core .* s
    else
        #display.(("a", size(coeffs_bra), size(coeffs_ket), "sizes: ", size.(overlaps), indices))
        #display.(("a", size(coeffs_bra), size(coeffs_ket), "sizes: ", overlaps, indices))
        bra_core = @ncon(tensors, indices)
        bra_core .= bra_core .* s
    end

    new_factors = [coeffs_ket.factors[i] for i in 1:N]
    new_factors[c1.idx] = new_factor1
    new_factors[c2.idx] = new_factor2
    new_factors[c3.idx] = new_factor3
    return Tucker(bra_core, NTuple{N}(new_factors))

end
#=}}}=#

function form_sigma_block_expand(term::ClusteredTerm4B,
                            cluster_ops::Vector{ClusterOps},
                            fock_bra::FockConfig, bra::TuckerConfig,
                            fock_ket::FockConfig, ket::TuckerConfig,
                            coeffs_ket::Tucker{T,N};
                            max_number=nothing, prescreen=1e-4) where {T,N}
#={{{=#
    #display(term)
    #display.((fock_bra, fock_ket))
    c1 = term.clusters[1]
    c2 = term.clusters[2]
    c3 = term.clusters[3]
    c4 = term.clusters[4]
    n_clusters = length(bra)

    #
    # make sure active clusters are correct transitions
    fock_bra[c1.idx] == fock_ket[c1.idx] .+ term.delta[1] || throw(Exception)
    fock_bra[c2.idx] == fock_ket[c2.idx] .+ term.delta[2] || throw(Exception)
    fock_bra[c3.idx] == fock_ket[c3.idx] .+ term.delta[3] || throw(Exception)
    fock_bra[c4.idx] == fock_ket[c4.idx] .+ term.delta[4] || throw(Exception)

    #
    # determine sign from rearranging clusters if odd number of operators
    state_sign = compute_terms_state_sign(term, fock_ket)

    #
    # op[IK,JL] = <I|p'|J> h(pq) <K|q|L>
    # todo: add in 2e integral tucker decomposition and compress gamma along 1st index first

    #
    # Compress Gammas using the cluster's Tucker factors, but since we are expanding the compression space
    # only compress the right hand side
    # e.g.,
    #   Gamma(pqr, I, J) Ur(J,l) = Gamma(pqr, I, l) where k and l are compressed indices
    @views gamma1 = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])][:,bra[c1.idx],ket[c1.idx]]
    Ur = coeffs_ket.factors[c1.idx]
    @tensor begin
        g1[p,I,l] := Ur[J,l] * gamma1[p,I,J]
    end

    @views gamma2 = cluster_ops[c2.idx][term.ops[2]][(fock_bra[c2.idx],fock_ket[c2.idx])][:,bra[c2.idx],ket[c2.idx]]
    Ur = coeffs_ket.factors[c2.idx]
    @tensor begin
        g2[p,I,l] := Ur[J,l] * gamma2[p,I,J]
    end

    @views gamma3 = cluster_ops[c3.idx][term.ops[3]][(fock_bra[c3.idx],fock_ket[c3.idx])][:,bra[c3.idx],ket[c3.idx]]
    Ur = coeffs_ket.factors[c3.idx]
    @tensor begin
        g3[p,I,l] := Ur[J,l] * gamma3[p,I,J]
    end

    @views gamma4 = cluster_ops[c4.idx][term.ops[4]][(fock_bra[c4.idx],fock_ket[c4.idx])][:,bra[c4.idx],ket[c4.idx]]
    Ur = coeffs_ket.factors[c4.idx]
    @tensor begin
        g4[p,I,l] := Ur[J,l] * gamma4[p,I,J]
    end

    #
    # Decompose the local operators. Since gamma[p,I,l] has indices (small, large, small),
    # we only need at most p*l number of new vectors for the index we are searching over

    new_factor1 = Matrix(1.0I, size(g1,2), size(g1,2))
    new_factor2 = Matrix(1.0I, size(g2,2), size(g2,2))
    new_factor3 = Matrix(1.0I, size(g3,2), size(g3,2))
    new_factor4 = Matrix(1.0I, size(g4,2), size(g4,2))

    D = permutedims(g1, [2,1,3])
    F = svd(reshape(D, size(D,1), size(D,2)*size(D,3)))
    nkeep = 0
    for si in F.S
        if si > prescreen
            nkeep += 1
        end
    end
    new_factor1 = F.U[:,1:nkeep]
    g1 = Diagonal(F.S[1:nkeep]) * F.Vt[1:nkeep,:]
    g1 = reshape(g1, size(g1,1), size(D,2), size(D,3))
    g1 = permutedims(g1, [2,1,3])

    D = permutedims(g2, [2,1,3])
    F = svd(reshape(D, size(D,1), size(D,2)*size(D,3)))
    nkeep = 0
    for si in F.S
        if si > prescreen
            nkeep += 1
        end
    end
    new_factor2 = F.U[:,1:nkeep]
    g2 = Diagonal(F.S[1:nkeep]) * F.Vt[1:nkeep,:]
    g2 = reshape(g2, size(g2,1), size(D,2), size(D,3))
    g2 = permutedims(g2, [2,1,3])

    D = permutedims(g3, [2,1,3])
    F = svd(reshape(D, size(D,1), size(D,2)*size(D,3)))
    nkeep = 0
    for si in F.S
        if si > prescreen
            nkeep += 1
        end
    end
    new_factor3 = F.U[:,1:nkeep]
    g3 = Diagonal(F.S[1:nkeep]) * F.Vt[1:nkeep,:]
    g3 = reshape(g3, size(g3,1), size(D,2), size(D,3))
    g3 = permutedims(g3, [2,1,3])

    D = permutedims(g4, [2,1,3])
    F = svd(reshape(D, size(D,1), size(D,2)*size(D,3)))
    nkeep = 0
    for si in F.S
        if si > prescreen
            nkeep += 1
        end
    end
    new_factor4 = F.U[:,1:nkeep]
    g4 = Diagonal(F.S[1:nkeep]) * F.Vt[1:nkeep,:]
    g4 = reshape(g4, size(g4,1), size(D,2), size(D,3))
    g4 = permutedims(g4, [2,1,3])

    #
    # Now contract into 4body term
    #
    # h(p,q) * g1(p,I,J) * g2(q,K,L) = op(J,L,I,K)
    op = Array{Float64}[]
    @tensor begin
        op[q,r,s,J,I] := term.ints[p,q,r,s] * g1[p,I,J]
        op[r,s,J,L,I,K] := op[q,r,s,J,I] * g2[q,K,L]
        op[s,J,L,N,I,K,M] := op[r,s,J,L,I,K] * g3[r,M,N]
        op[J,L,N,P,I,K,M,O] := op[s,J,L,N,I,K,M] * g4[s,O,P]
    end

    #
    # form overlaps - needed when TuckerConfigs aren't the same because each does their own compression and has
    # distinct Tucker factors
    tensors = Vector{Array{T}}()
    indices = Vector{Vector{Int16}}()
    state_indices = -collect(1:n_clusters)
    s = state_sign # this is the product of scalar overlaps that don't need tensor contractions

    # if the compressed operator becomes a scalar, treat it as such
    if length(op) == 1
        s *= op[1]
    else
        op_indices = [c1.idx, c2.idx, c3.idx, c4.idx, -c1.idx, -c2.idx, -c3.idx, -c4.idx]
        state_indices[c1.idx] = c1.idx
        state_indices[c2.idx] = c2.idx
        state_indices[c3.idx] = c3.idx
        state_indices[c4.idx] = c4.idx
        push!(tensors, op)
        push!(indices, op_indices)
    end

    push!(tensors, coeffs_ket.core)
    push!(indices, state_indices)

    length(tensors) == length(indices) || error(" mismatch between operators and indices")

    bra_core = zeros(1,1)
    if length(tensors) == 1
        # this means that all the overlaps and the operator is a scalar
        bra_core = coeffs_ket.core .* s
    else
        #display.(("a", size(coeffs_bra), size(coeffs_ket), "sizes: ", size.(overlaps), indices))
        #display.(("a", size(coeffs_bra), size(coeffs_ket), "sizes: ", overlaps, indices))
        bra_core = @ncon(tensors, indices)
        bra_core .= bra_core .* s
    end

    new_factors = [coeffs_ket.factors[i] for i in 1:N]
    new_factors[c1.idx] = new_factor1
    new_factors[c2.idx] = new_factor2
    new_factors[c3.idx] = new_factor3
    new_factors[c4.idx] = new_factor4
    return Tucker(bra_core, NTuple{N}(new_factors))

end
#=}}}=#

function form_sigma_block_expand2(term::ClusteredTerm2B,
                            cluster_ops::Vector{ClusterOps},
                            fock_bra::FockConfig, bra::TuckerConfig,
                            fock_ket::FockConfig, ket::TuckerConfig,
                            coeffs_ket::Tucker{T,N},
                            scr::Vector{Vector{Float64}};
                            max_number=nothing, prescreen=1e-4) where {T,N}
    #={{{=#
    #display(term)
    #display.((fock_bra, fock_ket))
    c1 = term.clusters[1]
    c2 = term.clusters[2]
    n_clusters = length(bra)

    #
    # make sure active clusters are correct transitions
    fock_bra[c1.idx] == fock_ket[c1.idx] .+ term.delta[1] || throw(Exception)
    fock_bra[c2.idx] == fock_ket[c2.idx] .+ term.delta[2] || throw(Exception)

    #
    # determine sign from rearranging clusters if odd number of operators
    state_sign = compute_terms_state_sign(term, fock_ket)

    scr1 = scr[1]
    scr2 = scr[2]
    scr3 = scr[3]
    scr4 = scr[4]

    np = size(term.ints,1)
    nq = size(term.ints,2)

    #
    # op[IK,JL] = <I|p'|J> h(pq) <K|q|L>

    # todo: add in 2e integral tucker decomposition and compress gamma along 1st index first

    #
    # Compress Gammas using the cluster's Tucker factors, but since we are expanding the compression space
    # only compress the right hand side
    # e.g.,
    #   Gamma(pqr, I, J) Ur(J,l) = Gamma(pqr, I, l) where k and l are compressed indices
    g1a::Array{Float64,3} = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])]
    #@btime @views gamma1 = $g1a[:,$bra[$c1.idx],$ket[$c1.idx]]
    #@views gamma1 = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])][:,bra[c1.idx],ket[c1.idx]]
    @views gamma1 = g1a[:,bra[c1.idx],ket[c1.idx]]
    Ur = coeffs_ket.factors[c1.idx]

    # gamma1(p,I,J) U(J,j)  = g1(p,I,j)
    nI = size(gamma1, 2)
    nJ = size(gamma1, 3)
    nj = size(Ur, 2)

    #resize!(scr1, np*nI*nj)
    #scr1 = reshape2(scr1, (np*nI,nj))
    #mul!(scr1, reshape2(gamma1,(np*nI,nJ)), Ur)
    #g1 = reshape2(scr1, (np,nI,nj))

    g1 = scr[1]
    resize!(g1, np*nI*nj)
    g1 = reshape2(g1, (np,nI,nj))
    @tensor begin
        g1[p,I,l] = Ur[J,l] * gamma1[p,I,J]
        #g1[p,I,l] := Ur[J,l] * gamma1[p,I,J]
    end

    @views gamma2 = cluster_ops[c2.idx][term.ops[2]][(fock_bra[c2.idx],fock_ket[c2.idx])][:,bra[c2.idx],ket[c2.idx]]
    Ur = coeffs_ket.factors[c2.idx]
    @tensor begin
        g2[p,I,l] := Ur[J,l] * gamma2[p,I,J]
    end

    #
    # Decompose the local operators. Since gamma[p,I,l] has indices (small, large, small),
    # we only need at most p*l number of new vectors for the index we are searching over

    new_factor1 = Matrix(1.0I, size(g1,2), size(g1,2))
    new_factor2 = Matrix(1.0I, size(g2,2), size(g2,2))

    D = permutedims(g1, [2,1,3])
    F = svd(reshape(D, size(D,1), size(D,2)*size(D,3)))
    nkeep = 0
    for si in F.S
        if si > prescreen
            nkeep += 1
        end
    end
    new_factor1 = F.U[:,1:nkeep]
    g1 = Diagonal(F.S[1:nkeep]) * F.Vt[1:nkeep,:]
    g1 = reshape(g1, size(g1,1), size(D,2), size(D,3))
    g1 = permutedims(g1, [2,1,3])

    D = permutedims(g2, [2,1,3])
    F = svd(reshape(D, size(D,1), size(D,2)*size(D,3)))
    nkeep = 0
    for si in F.S
        if si > prescreen
            nkeep += 1
        end
    end
    new_factor2 = F.U[:,1:nkeep]
    g2 = Diagonal(F.S[1:nkeep]) * F.Vt[1:nkeep,:]
    g2 = reshape(g2, size(g2,1), size(D,2), size(D,3))
    g2 = permutedims(g2, [2,1,3])

    #
    # Now contract into 2body term
    #
    # h(p,q) * g1(p,I,J) * g2(q,K,L) = op(J,L,I,K)
    op = Array{Float64}[]
    @tensor begin
        op[q,J,I] := term.ints[p,q] * g1[p,I,J]
        op[J,L,I,K] := op[q,J,I] * g2[q,K,L]
    end

    #
    # form overlaps - needed when TuckerConfigs aren't the same because each does their own compression and has
    # distinct Tucker factors
    tensors = Vector{Array{T}}()
    indices = Vector{Vector{Int16}}()
    state_indices = -collect(1:n_clusters)
    s = state_sign # this is the product of scalar overlaps that don't need tensor contractions

    # if the compressed operator becomes a scalar, treat it as such
    if length(op) == 1
        s *= op[1]
    else
        op_indices = [c1.idx, c2.idx, -c1.idx, -c2.idx]
        state_indices[c1.idx] = c1.idx
        state_indices[c2.idx] = c2.idx
        push!(tensors, op)
        push!(indices, op_indices)
    end

    push!(tensors, coeffs_ket.core)
    push!(indices, state_indices)

    length(tensors) == length(indices) || error(" mismatch between operators and indices")

    bra_core = zeros(1,1)
    if length(tensors) == 1
        # this means that all the overlaps and the operator is a scalar
        bra_core = coeffs_ket.core .* s
    else
        #display.(("a", size(coeffs_bra), size(coeffs_ket), "sizes: ", size.(overlaps), indices))
        #display.(("a", size(coeffs_bra), size(coeffs_ket), "sizes: ", overlaps, indices))
        bra_core = @ncon(tensors, indices)
        bra_core .= bra_core .* s
    end

    # first decompose the already partially decomposed core tensor
    #
    # Vket ~ Aket x U1 x U2 x ...
    #
    # then we modified the compressed coefficients in the ket, Aket
    #
    # to get Abra, which we then compress again.
    #
    # The active cluster tucker factors become identity matrices
    #
    # Abra ~ Bbra x UU1 x UU2 x ....
    #
    #
    # e.g, V(IJK) = C(ijk) * U1(Ii) * U2(Jj) * U3(Kk)
    #
    # then C get's modified and furhter compressed
    #
    # V(IJK) = C(abc) * U1(ia) * U2(jb) * U3(kc) * U1(Ii) * U2(Jj) * U3(Kk)
    # V(IJK) = C(abc) * (U1(ia) * U1(Ii)) * (U2(jb) * U2(Jj)) * (U3(kc) * U3(Kk))
    # V(IJK) = C(abc) * U1(Ia)  * U2(Jb) * U3(Kc)
    #

    new_factors = [coeffs_ket.factors[i] for i in 1:N]
    #new_factors[c1.idx] = Matrix(1.0I, size(bra_core,c1.idx), size(bra_core,c1.idx))
    #new_factors[c2.idx] = Matrix(1.0I, size(bra_core,c2.idx), size(bra_core,c2.idx))
    new_factors[c1.idx] = new_factor1
    new_factors[c2.idx] = new_factor2
    return Tucker(bra_core, NTuple{N}(new_factors))
end
#=}}}=#

function calc_bound(term::ClusteredTerm1B,
                            cluster_ops::Vector{ClusterOps},
                            fock_bra::FockConfig, bra::TuckerConfig,
                            fock_ket::FockConfig, ket::TuckerConfig,
                            coeffs_ket::Tucker{T,N};
                            prescreen=1e-4) where {T,N}
    c1 = term.clusters[1]

    bound1 = norm(term.ints)*norm(coeffs_ket.core)
    if bound1 < sqrt(prescreen)
        return false
    end
    return true
end

function calc_bound(term::ClusteredTerm2B,
                            cluster_ops::Vector{ClusterOps},
                            fock_bra::FockConfig, bra::TuckerConfig,
                            fock_ket::FockConfig, ket::TuckerConfig,
                            coeffs_ket::Tucker{T,N};
                            prescreen=1e-4) where {T,N}
    c1 = term.clusters[1]
    c2 = term.clusters[2]

    #@views gamma1 = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])][:,bra[c1.idx],ket[c1.idx]]
    #@views gamma2 = cluster_ops[c2.idx][term.ops[2]][(fock_bra[c2.idx],fock_ket[c2.idx])][:,bra[c2.idx],ket[c2.idx]]
    bound1 = norm(term.ints)*norm(coeffs_ket.core)
    #bound1 *= norm(gamma1)*norm(gamma2)
    bound1 *= norm(coeffs_ket.factors[c1.idx])*norm(coeffs_ket.factors[c2.idx])
    if bound1 < sqrt(prescreen)
        return false
    end
    return true
end

function calc_bound(term::ClusteredTerm3B,
                            cluster_ops::Vector{ClusterOps},
                            fock_bra::FockConfig, bra::TuckerConfig,
                            fock_ket::FockConfig, ket::TuckerConfig,
                            coeffs_ket::Tucker{T,N};
                            prescreen=1e-4) where {T,N}
    c1 = term.clusters[1]
    c2 = term.clusters[2]
    c3 = term.clusters[3]

    @views gamma1 = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])][:,bra[c1.idx],ket[c1.idx]]
    @views gamma2 = cluster_ops[c2.idx][term.ops[2]][(fock_bra[c2.idx],fock_ket[c2.idx])][:,bra[c2.idx],ket[c2.idx]]
    @views gamma3 = cluster_ops[c3.idx][term.ops[3]][(fock_bra[c3.idx],fock_ket[c3.idx])][:,bra[c3.idx],ket[c3.idx]]
    bound1 = norm(term.ints)*norm(coeffs_ket.core)
    #bound1 *= norm(gamma1)*norm(gamma2)*norm(gamma3)
    bound1 *= norm(coeffs_ket.factors[c1.idx])*norm(coeffs_ket.factors[c2.idx])*norm(coeffs_ket.factors[c3.idx])
    if bound1 < sqrt(prescreen)
        return false
    end
    return true
end

function calc_bound(term::ClusteredTerm4B,
                            cluster_ops::Vector{ClusterOps},
                            fock_bra::FockConfig, bra::TuckerConfig,
                            fock_ket::FockConfig, ket::TuckerConfig,
                            coeffs_ket::Tucker{T,N};
                            prescreen=1e-4) where {T,N}
    c1 = term.clusters[1]
    c2 = term.clusters[2]
    c3 = term.clusters[3]
    c4 = term.clusters[4]

    @views gamma1 = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])][:,bra[c1.idx],ket[c1.idx]]
    @views gamma2 = cluster_ops[c2.idx][term.ops[2]][(fock_bra[c2.idx],fock_ket[c2.idx])][:,bra[c2.idx],ket[c2.idx]]
    @views gamma3 = cluster_ops[c3.idx][term.ops[3]][(fock_bra[c3.idx],fock_ket[c3.idx])][:,bra[c3.idx],ket[c3.idx]]
    @views gamma4 = cluster_ops[c4.idx][term.ops[4]][(fock_bra[c4.idx],fock_ket[c4.idx])][:,bra[c4.idx],ket[c4.idx]]
    bound1 = norm(term.ints)*norm(coeffs_ket.core)
    #bound1 *= norm(gamma1)*norm(gamma2)*norm(gamma3)*norm(gamma4)
    bound1 *= norm(coeffs_ket.factors[c1.idx])*norm(coeffs_ket.factors[c2.idx])*norm(coeffs_ket.factors[c3.idx])*norm(coeffs_ket.factors[c4.idx])
    if bound1 < sqrt(prescreen)
        return false
    end
    return true
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
""" Cálculo de IMC
    Leia o peso e a altura do usuário
    imprima seu Índice de Massa Corporal.
        * IMC = peso / (altura^2)
    OBS: use ponto final para informar o peso
"""

println("informe seu peso")
weight = parse(Float32, readline())

println("informe sua altura")
height = parse(Float32, readline())

imc = (weight / (height^2))
println("IMC = ", imc)

max_stars_repo_name: codinginbrazil/functional-programming

max_issues_repo_name: codinginbrazil/functional-programming

max_forks_repo_name: codinginbrazil/functional-programming