text:
"""
    $SIGNATURES

Based on a push model from child cliques that should have already completed their computation.
"""
getCliqInitUpMsgs(cliq::TreeClique)::Dict{Int, TempBeliefMsg} = getData(cliq).upInitMsgs

function setCliqUpInitMsgs!(cliq::TreeClique, childid::Int, msg::TempBeliefMsg)
  cd = getData(cliq)
  soco = getSolveCondition(cliq)
  lockUpStatus!(cd)
  cd.upInitMsgs[childid] = msg
  # notify cliq condition that there was a change
  notify(soco)
  unlockUpStatus!(cd)
  #hack for mitigating deadlocks, in case a user was not already waiting, but waiting on lock instead
  sleep(0.1)
  notify(soco)
  nothing
end

function isCliqInitialized(cliq::TreeClique)::Bool
  return getData(cliq).initialized in [:initialized; :upsolved]
end

function isCliqUpSolved(cliq::TreeClique)::Bool
  return getData(cliq).initialized == :upsolved
end

"""
    $SIGNATURES

Return the most likely  ordering for initializing factor (assuming up solve
sequence).

Notes:
- sorts id for increasing number of connected factors.
"""
function getCliqVarInitOrderUp(cliq::TreeClique)
  # rules to explore dimension from one to the other?

  # get all variable ids and number of associated factors
  allids = getCliqAllVarIds(cliq)
  nfcts = getCliqNumAssocFactorsPerVar(cliq)

  # get priors and singleton message variables (without partials)
  prids = getCliqVarIdsPriors(cliq, getCliqAllVarIds(cliq), false)

  # get current up msgs in the init process (now have all singletons)
  upmsgs = getCliqInitUpMsgs(cliq)
  upmsgids = collect(keys(upmsgs))

  # all singleton variables
  singids = union(prids, upmsgids)

  # add msg marginal prior (singletons) to number of factors
  for msid in upmsgids
    nfcts[msid .== allids] .+= 1
  end

  # sort permutation order for increasing number of factor association
  nfctsp = sortperm(nfcts)
  sortedids = allids[nfctsp]

  # organize the prior variables separately with asceding factor count
  initorder = Symbol[]
  for id in sortedids
    if id in singids
      push!(initorder, id)
    end
  end
  # in ascending order of number of factors
  for id in sortedids
    if !(id in initorder)
      push!(initorder, id)
    end
  end
  return initorder
end

lockUpStatus!(cdat::BayesTreeNodeData, idx::Int=1) = put!(cdat.lockUpStatus, idx)
lockUpStatus!(cliq::TreeClique, idx::Int=1) = lockUpStatus!(getData(cliq), idx)
unlockUpStatus!(cdat::BayesTreeNodeData) = take!(cdat.lockUpStatus)
unlockUpStatus!(cliq::TreeClique) = unlockUpStatus!(getData(cliq))

function lockDwnStatus!(cdat::BayesTreeNodeData, idx::Int=1; logger=ConsoleLogger())
  with_logger(logger) do
    @info "lockDwnStatus! isready=$(isready(cdat.lockDwnStatus)) with data $(cdat.lockDwnStatus.data)"
  end
    flush(logger.stream)
  # if isready(cdat.lockDwnStatus)
  #   with_logger(logger) do
  #     @info "lockDwnStatus! is locked with $(cdat.lockDwnStatus.data)"
  #   end
  #   flush(logger.stream)
  # end
  stf =  put!(cdat.lockDwnStatus, idx)
  with_logger(logger) do
    @info "lockDwnStatus! DONE: isready=$(isready(cdat.lockDwnStatus)) with data $(cdat.lockDwnStatus.data)"
  end
  flush(logger.stream)
  stf
end
unlockDwnStatus!(cdat::BayesTreeNodeData) = take!(cdat.lockDwnStatus)

"""
    $SIGNATURES

Update clique status and notify of the change

Notes
- Assumes users will lock the status state before getting status until after decision whether to update status.
- If so, only unlock after status and condition has been updated.

Dev Notes
- Should be made an atomic transaction
"""
function notifyCliqUpInitStatus!(cliq::TreeClique, status::Symbol; logger=ConsoleLogger())
  cd = getData(cliq)
  with_logger(logger) do
    tt = split(string(now()), 'T')[end]
    @info "$(tt) $(current_task()), cliq=$(cliq.index), notifyCliqUpInitStatus! -- pre-lock, $(cd.initialized)-->$(status)"
  end
  flush(logger.stream)

  ## TODO only notify if not data structure is not locked by other user (can then remove the hack)
  # Wait until lock can be aquired
  lockUpStatus!(cd)

  cd.initialized = status
  if isready(cd.initUpChannel)
    tkst = take!(cd.initUpChannel)
    # @info "dumping stale cliq=$(cliq.index) status message $(tkst), replacing with $(status)"
  end
  put!(cd.initUpChannel, status)
  cond = getSolveCondition(cliq)
  notify(cond)
    # hack to avoid a race condition  -- remove with atomic lock logic upgrade
    sleep(0.1)
    notify(cond) # getSolveCondition(cliq)

  # TODO unlock
  unlockUpStatus!(cd)
  with_logger(logger) do
    tt = split(string(now()), 'T')[end]
    @info "$(tt) $(current_task()), cliq=$(cliq.index), notifyCliqUpInitStatus! -- unlocked, $(cd.initialized)"
  end

  nothing
end

function notifyCliqDownInitStatus!(cliq::TreeClique, status::Symbol; logger=ConsoleLogger())
  cdat = getData(cliq)
  with_logger(logger) do
    @info "$(now()) $(current_task()), cliq=$(cliq.index), notifyCliqDownInitStatus! -- pre-lock, new $(cdat.initialized)-->$(status)"
  end

  # take lock for atomic transaction
  lockDwnStatus!(cdat, cliq.index, logger=logger)

  cdat.initialized = status

  if isready(cdat.initDownChannel)
    content = take!(cdat.initDownChannel)
    with_logger(logger) do
      @info "dumping stale cliq=$(cliq.index) status message $(content), replacing with $(status)"
    end
  end
  put!(cdat.initDownChannel, status)
  notify(getSolveCondition(cliq))
    # hack to avoid a race condition
    sleep(0.1)
    notify(getSolveCondition(cliq))

  # unlock for others to proceed
  unlockDwnStatus!(cdat)
  with_logger(logger) do
    @info "$(now()), cliq=$(cliq.index), notifyCliqDownInitStatus! -- unlocked, $(getCliqStatus(cliq))"
  end

  # flush(logger.stream)

  nothing
end

"""
    $SIGNATURES

Return true if clique has completed the local upward direction inference procedure.
"""
isUpInferenceComplete(cliq::TreeClique) = getData(cliq).upsolved

function areCliqVariablesAllInitialized(dfg::G, cliq::TreeClique) where {G <: AbstractDFG}
  allids = getCliqAllVarIds(cliq)
  isallinit = true
  for vid in allids
    var = DFG.getVariable(dfg, vid)
    isallinit &= isInitialized(var)
  end
  isallinit
end

"""
   $SIGNATURES

Determine if this `cliq` has been fully initialized and child cliques have completed their full upward inference.
"""
function isCliqReadyInferenceUp(fgl::FactorGraph, tree::BayesTree, cliq::TreeClique)
  isallinit = areCliqVariablesAllInitialized(fgl, cliq)

  # check that all child cliques have also completed full up inference.
  for chl in getChildren(tree, cliq)
    isallinit &= isUpInferenceComplete(chl)
  end
  return isallinit
end

"""
    $SIGNATURES

Blocking call until `cliq` upInit processes has arrived at a result.
"""
function getCliqInitUpResultFromChannel(cliq::TreeClique)
  status = take!(getData(cliq).initUpChannel)
  @info "$(current_task()) Clique $(cliq.index), dumping initUpChannel status, $status"
  return status
end

"""
    $SIGNATURES

Return `::Symbol` status a particular clique is in, with specific regard to solution
or numerical initialization status:
- :needdownmsg
- :upsolved
- :downsolved
- :initialized
- :marginalized
- :null

Notes:
- `:null` represents the first uninitialized state of a cliq.
"""
getCliqStatus(cliqdata::BayesTreeNodeData)::Symbol = cliqdata.initialized
getCliqStatus(cliq::TreeClique)::Symbol = getCliqStatus(getData(cliq))

getCliqStatusUp(cliq::TreeClique)::Symbol = getCliqStatus(cliq)

"""
    $SIGNATURES

Set up initialization or solve status of this `cliq`.
"""
function setCliqStatus!(cliq::TreeClique, status::Symbol)
  getData(cliq).initialized = status
end

"""
    $SIGNATURES

Return true if all variables in clique are considered marginalized (and initialized).
"""
function areCliqVariablesAllMarginalized(subfg::AbstractDFG,
                                         cliq::TreeClique)
  for vsym in getCliqAllVarIds(cliq)
    vert = getVariable(subfg, vsym)
    if !isMarginalized(vert) || !isInitialized(vert)
      return false
    end
  end
  return true
end

"""
    $SIGNATURES

Set all Bayes (Junction) tree cliques that have all marginalized and initialized variables.
"""
function setTreeCliquesMarginalized!(dfg::G,
                                     tree::BayesTree) where G <: AbstractDFG
  #
  for (cliid, cliq) in tree.cliques
    if areCliqVariablesAllMarginalized(dfg, cliq)
      # need to set the upward messages
      msgs = prepCliqInitMsgsUp(dfg, cliq)
      setUpMsg!(cliq, msgs)

      prnt = getParent(tree, cliq)
      if length(prnt) > 0
        # THIS IS FOR INIT PASSES ONLY
        setCliqUpInitMsgs!(prnt[1], cliq.index, msgs)
      end

      setCliqStatus!(cliq, :marginalized)
      setCliqDrawColor(cliq, "blue")
    end
  end
  nothing
end

function blockCliqUntilParentDownSolved(prnt::TreeClique; logger=ConsoleLogger())::Nothing
  #
  lbl = getLabel(prnt)

  with_logger(logger) do
    @info "blockCliqUntilParentDownSolved, prntcliq=$(prnt.index) | $lbl | going to fetch initdownchannel..."
  end
  flush(logger.stream)
  while fetch(getData(prnt).initDownChannel) != :downsolved
    # @sync begin
    #   @async begin
    #     sleep(1)
    #     notify(getSolveCondition(prnt))
    #   end
      with_logger(logger) do
        @info "blockCliqUntilParentDownSolved, prntcliq=$(prnt.index) | $lbl | waiting on solve condition..."
      end
      wait(getSolveCondition(prnt))
    # end
  end

  return nothing
end

"""
    $SIGNATURES

Block the thread until child cliques of `prnt::TreeClique` have finished
attempting upward initialization -- i.e. have status result.
Return `::Dict{Symbol}` indicating whether next action that should be taken
for each child clique.

Notes:
- See status options at `getCliqStatusUp(..)`.
- Can be called multiple times
"""
function blockCliqUntilChildrenHaveUpStatus(tree::BayesTree,
                                            prnt::TreeClique,
                                            logger=ConsoleLogger() )::Dict{Int, Symbol}
  #
  ret = Dict{Int, Symbol}()
  chlr = getChildren(tree, prnt)
  for ch in chlr
    # either wait to fetch new result, or report or result
    chst = getCliqStatusUp(ch)
    with_logger(logger) do
      @info "cliq $(prnt.index), child $(ch.index) status is $(chst), isready(initUpCh)=$(isready(getData(ch).initUpChannel))."
    end
    flush(logger.stream)
    ret[ch.index] = fetch(getData(ch).initUpChannel)
  end
  with_logger(logger) do
      @info "cliq $(prnt.index), fetched all."
  end
  return ret
end

"""
    $SIGNATURES

Wait here if all siblings and the parent status are `:needdownmsg`.
Return true when parent is `:initialized` after all were `:needdownmsg`

Notes
- used for regulating long need down message chains.
- exit strategy is parent becomes status `:initialized`.
"""
function blockCliqSiblingsParentNeedDown(tree::BayesTree,
                                         cliq::TreeClique; logger=ConsoleLogger())
  #
  with_logger(logger) do
    @info "cliq $(cliq.index), blockCliqSiblingsParentNeedDown -- start of function"
  end
  # ret = Dict{Int, Symbol}()
  # flush(logger.stream)
  prnt = getParent(tree, cliq)
  allneeddwn = true
  if length(prnt) > 0
    prstat = getCliqStatus(prnt[1])
    if prstat == :needdownmsg
      for ch in getChildren(tree, prnt[1])
        chst = getCliqStatusUp(ch)
        if chst != :needdownmsg
          allneeddwn = false
          break;
        end
      end

      if allneeddwn
        with_logger(logger) do
          tt = split(string(now()), 'T')[end]
          @warn "$(tt) | $(current_task()), cliq $(cliq.index), block since all siblings/parent($(prnt[1].index)) :needdownmsg."
        end
        flush(logger.stream)
        # do actual fetch
        prtmsg = fetch(getData(prnt[1]).initDownChannel)
        with_logger(logger) do
            tt = split(string(now()), 'T')[end]
          @info "$tt | $(current_task()) clique $(prnt[1].index), blockCliqSiblingsParentNeedDown -- after fetch $prstat, $prtmsg"
        end
        if prtmsg == :initialized
          return true
        else
            with_logger(logger) do
                tt = split(string(now()), 'T')[end]
                @warn "$tt | $(current_task()) Clique $(prnt[1].index), maybe clear down init message $prtmsg"
            end
          # take!(getData(prnt[1]).initDownChannel)
        end
      end
    end
  end
  return false
end

"""
    $SIGNATURES

Cycle through var order and initialize variables as possible in `subfg::FactorGraph`.
Return true if something was updated.

Notes:
- assumed `subfg` is a subgraph containing only the factors that can be used.
  - including the required up or down messages
- intended for both up and down initialization operations.

Dev Notes
- Should monitor updates based on the number of inferred & solvable dimensions
"""
function cycleInitByVarOrder!(subfg::G,
                              varorder::Vector{Symbol};
                              logger=ConsoleLogger()  )::Bool where G <: AbstractDFG
  #
  with_logger(logger) do
    @info "cycleInitByVarOrder! -- varorder=$(varorder)"
  end
  retval = false
  count = 1
  while count > 0
    count = 0
    for vsym in varorder
      var = DFG.getVariable(subfg, vsym)
      isinit = isInitialized(var)
      with_logger(logger) do
        @info "var.label=$(var.label) is initialized=$(isinit)"
      end
      doautoinit!(subfg, [var;], logger=logger)
      if isinit != isInitialized(var)
        count += 1
        retval = true
      end
    end
  end
  with_logger(logger) do
    @info "cycleInitByVarOrder!, retval=$(retval)"
  end
  flush(logger.stream)
  return retval
end

"""
    $SIGNATURES

Update `subfg<:AbstractDFG` according to internal computations for a full upsolve.
"""
function doCliqUpSolve!(subfg::G,
                        tree::BayesTree,
                        cliq::TreeClique;
                        multiproc::Bool=true,
                        logger=ConsoleLogger()  )::Symbol where G <: AbstractDFG
  #
  csym = getCliqFrontalVarIds(cliq)[1]
  # csym = DFG.getVariable(subfg, getCliqFrontalVarIds(cliq)[1]).label # ??
  approxCliqMarginalUp!(subfg, tree, csym, false, N=getSolverParams(subfg).N, logger=logger, multiproc=multiproc)
  getData(cliq).upsolved = true
  return :upsolved
end

"""
    $SIGNATURES

Prepare the upward inference messages from clique to parent and return as `Dict{Symbol}`.
"""
function prepCliqInitMsgsUp(subfg::G,
                            cliq::TreeClique,
                            logger=ConsoleLogger() )::TempBeliefMsg  where G <: AbstractDFG
  #
  # construct init's up msg to place in parent from initialized separator variables
  msg = TempBeliefMsg()
  seps = getCliqSeparatorVarIds(cliq)
  with_logger(logger) do
    @info "prepCliqInitMsgsUp, seps=$seps"
  end
  for vid in seps
    var = DFG.getVariable(subfg, vid)
    if isInitialized(var)
      msg[Symbol(var.label)] = (getKDE(var), solverData(var).inferdim)
    end
  end
  return msg
end

function prepCliqInitMsgsUp(subfg::G, tree::BayesTree, cliq::TreeClique)::TempBeliefMsg  where G <: AbstractDFG
  @warn "deprecated, use prepCliqInitMsgsUp(subfg::FactorGraph, cliq::TreeClique) instead"
  prepCliqInitMsgsUp(subfg, cliq)
end

"""
    $SIGNATURES

Perform cliq initalization calculation based on current state of the tree and factor graph,
using upward message passing logic.

Notes
- adds msg priors added to clique subgraph
- Return either of (:initialized, :upsolved, :needdownmsg, :badinit)
- must use factors in cliq only, ensured by using subgraph -- TODO general case.
"""
function doCliqAutoInitUpPart1!(subfg::G,
                                tree::BayesTree,
                                cliq::TreeClique;
                                up_solve_if_able::Bool=true,
                                multiproc::Bool=true,
                                logger=ConsoleLogger() ) where {G <: AbstractDFG}
  #

  # init up msg has special procedure for incomplete messages
  varorder = Int[]

  # attempt initialize if necessary
  if !areCliqVariablesAllInitialized(subfg, cliq)
    # structure for all up message densities computed during this initialization procedure.
    varorder = getCliqVarInitOrderUp(cliq)
    # do physical inits, ignore cycle return value
    with_logger(logger) do
      @info "cliq $(cliq.index), doCliqAutoInitUpPart1! -- going for up cycle order"
    end

    cycleInitByVarOrder!(subfg, varorder, logger=logger)
    with_logger(logger) do
      @info "cliq $(cliq.index), doCliqAutoInitUpPart1! -- finished with up cycle order"
    end
  end
  flush(logger.stream)

  return nothing
end

"""
    $SIGNATURES

Follows cliq initalization calculation and attempts full upsolve
based on current state of the tree and factor graph,
using upward message passing logic.

Notes
- Removes msg priors added to clique subgraph
- Return either of (:initialized, :upsolved, :needdownmsg, :badinit)
- must use factors in cliq only, ensured by using subgraph -- TODO general case.
"""
function doCliqAutoInitUpPart2!(subfg::G,
                                tree::BayesTree,
                                cliq::TreeClique;
                                # msgfcts;
                                up_solve_if_able::Bool=true,
                                multiproc::Bool=true,
                                logger=ConsoleLogger()  )::Symbol where {G <: AbstractDFG}
  #
  cliqst = getCliqStatus(cliq)
  status = (cliqst == :initialized || length(getParent(tree, cliq)) == 0) ? cliqst : :needdownmsg

  # print out the partial init status of all vars in clique
  varids = getCliqAllVarIds(cliq)
  initstatus = Vector{Bool}(undef, length(varids))
  initpartial = Vector{Float64}(undef, length(varids))
  for i in 1:length(varids)
    initstatus[i] = solverData(getVariable(subfg, varids[i])).initialized
    initpartial[i] = solverData(getVariable(subfg, varids[i])).inferdim
  end
  with_logger(logger) do
    tt = split(string(now()),'T')[end]
    @info "$tt, cliq $(cliq.index), PARINIT: $varids | $initstatus | $initpartial"
  end

  # check if all cliq vars have been initialized so that full inference can occur on clique
  if areCliqVariablesAllInitialized(subfg, cliq)
    with_logger(logger) do
      tt = split(string(now()),'T')[end]
      @info "$(tt), cliq $(cliq.index), doCliqUpSolvePart2!, clique status = $(status)"
    end
    status = doCliqUpSolve!(subfg, tree, cliq, multiproc=multiproc, logger=logger)
  else
    with_logger(logger) do
      @info "cliq $(cliq.index), all variables not initialized, status = $(status)"
    end
  end

  # construct init's up msg to place in parent from initialized separator variables
  with_logger(logger) do
    @info "cliq $(cliq.index), going to prepCliqInitMsgsUp"
  end
  msg = prepCliqInitMsgsUp(subfg, cliq) # , tree

  # put the init result in the parent cliq.
  prnt = getParent(tree, cliq)
  with_logger(logger) do
    @info "cliq $(cliq.index), prnt = getParent(tree, cliq) = $(prnt)"
  end
  if length(prnt) > 0
    # not a root clique
    with_logger(logger) do
      tt = split(string(now()),'T')[end]
      @info "$tt, cliq $(cliq.index), doCliqAutoInitUpPart2! -- umsg in prnt=$(prnt[1].index), with $(collect(keys(msg)))"
    end
    # does internal notify on parent
    setCliqUpInitMsgs!(prnt[1], cliq.index, msg)
  end

  return status
end

# Helper function for prepCliqInitMsgsDown!
# populate products with products of upward messages
function condenseDownMsgsProductOnly!(fgl::G,
                                      products,
                                      msgspervar  ) where G <: AbstractDFG
  #
  # multiply multiple messages together
  for (msgsym, msgsBo) in msgspervar
    # check if this particular down message requires msgsym
    if exists(fgl, msgsym) # DFG.hasVariable(fgl, msgsym)
      if length(msgspervar[msgsym]) > 1
        msgs = getindex.(msgsBo, 1)
        haspars = 0.0
        for mb in msgsBo, val in mb[2]
            haspars += val
        end
        products[msgsym] = (manifoldProduct(msgs, getManifolds(fgl, msgsym)), haspars)
      else
        # transfer if only have a single belief
        products[msgsym] = (msgsBo[1][1], msgsBo[1][2])
      end
    else
      # not required, therefore remove from message to avoid confusion
      if haskey(products, msgsym)
        delete!(products, msgsym)
      end
    end
  end
  nothing
end

# currently for internal use only
# initialize variables based on best current achievable ordering
# OBVIOUSLY a lot of refactoring and consolidation needed with cliqGibbs / approxCliqMarginalUp
function initSolveSubFg!(subfg::G,
                         logger=ConsoleLogger() ) where G <: AbstractDFG
  #
  varorder = getSubFgPriorityInitOrder(subfg, logger)
  with_logger(logger) do
    @info "initSolveSubFg! -- varorder=$varorder"
  end
  cycleInitByVarOrder!(subfg, varorder, logger=logger)
  nothing
end

# Helper function for prepCliqInitMsgsDown!
# future, be used in a cached system with parent in one location only for all siblings
function condenseDownMsgsProductPrntFactors!(fgl::G,
                                             products,
                                             msgspervar,
                                             prnt::TreeClique,
                                             cliq::TreeClique,
                                             logger=ConsoleLogger()  ) where G <: AbstractDFG
  #

  # determine the variables of interest
  reqMsgIds = collect(keys(msgspervar))
  # unique frontals per cliq
  prntvars = intersect(getCliqSeparatorVarIds(cliq), getCliqAllVarIds(prnt))
  lvarids = union(prntvars, reqMsgIds)
  # determine allowable factors, if any (only from parent cliq)
  awfcts = getCliqAllFactIds(prnt)
  with_logger(logger) do
      @info "condenseDownMsgsProductPrntFactors! -- reqMsgIds=$(reqMsgIds),"
      @info "condenseDownMsgsProductPrntFactors! -- vars=$(lvarids),"
      @info "condenseDownMsgsProductPrntFactors! -- allow factors $awfcts"
  end

  # build required subgraph for parent/sibling down msgs
  lsfg = buildSubgraphFromLabels(fgl, lvarids)
  tempfcts = lsf(lsfg)
  dellist = setdiff(awfcts, tempfcts)
  for delf in dellist
    # TODO -- double check this deletefactor method is leaving the right parent sharing factor graph behind
    if exists(lsfg, delf) # hasFactor
      deleteFactor!(lsfg,delf)
    end
  end
  with_logger(logger) do
      @info "condenseDownMsgsProductPrntFactors! -- lsfg fcts=$(lsf(lsfg)),"
      @info "condenseDownMsgsProductPrntFactors! -- excess factors $dellist"
  end

  # add message priors
  addMsgFactors!(lsfg, msgspervar)

  # perform initialization/inference
  # ....uhhh TODO
  initSolveSubFg!(lsfg, logger)

      # QUICK DBG CODE
      vars = ls(lsfg)
      len = length(vars)
      tdims = Vector{Float64}(undef, len)
      isinit = Vector{Bool}(undef, len)
      for idx in 1:len
        tdims[idx] = getVariableSolvableDim(lsfg, vars[idx])
        isinit[idx] = isInitialized(lsfg, vars[idx])
      end
      with_logger(logger) do
          @info "condenseDownMsgsProductPrntFactors! -- after cycle init: vars=$vars"
          @info "condenseDownMsgsProductPrntFactors! -- after cycle init: tdims=$tdims"
          @info "condenseDownMsgsProductPrntFactors! -- after cycle init: isinit=$isinit"
      end
      # QUICK DBG CODE

  # extract complete downward marginal msg priors
  for id in intersect(getCliqSeparatorVarIds(cliq), lvarids)
    vari = getVariable(lsfg, id)
    products[id] = (getKDE(vari), getVariableInferredDim(vari))
  end

  # if recycling lsfg
  # deleteMsgFactors!(...)

  nothing
end

"""
    $SIGNATURES

Initialization downward message passing is different from regular inference since
it is possible that none of the child cliq variables have been initialized.

Notes
- init msgs from child upward passes are individually stored in this `cliq`.
- fresh product of overlapping beliefs are calculated on each function call.
- Assumed that `prnt` of siblings

Dev Notes
- This should be the initialization cycle of parent, build up bit by bit...
"""
function prepCliqInitMsgsDown!(fgl::G,
                               tree::BayesTree,
                               prnt::TreeClique,
                               cliq::TreeClique;
                               logger=ConsoleLogger(),
                               dbgnew::Bool=true  ) where G <: AbstractDFG
  #
  tt = split(string(now()), 'T')[end]
  with_logger(logger) do
    @info "$(tt) prnt $(prnt.index), prepCliqInitMsgsDown! --"
  end
  # get the current messages stored in the parent
  currmsgs = getCliqInitUpMsgs(prnt)
  with_logger(logger) do
    @info "prnt $(prnt.index), prepCliqInitMsgsDown! -- prnt ids::Int=$(collect(keys(currmsgs)))"
  end

  # check if any msgs should be multiplied together for the same variable
  msgspervar = Dict{Symbol, Vector{Tuple{BallTreeDensity,Float64}}}()
  for (prntid, msgs) in currmsgs
    for (msgsym, msg) in msgs
      if !haskey(msgspervar, msgsym)
        msgspervar[msgsym] = Vector{Tuple{BallTreeDensity,Float64}}()
      end
      with_logger(logger) do
        @info "prepCliqInitMsgsDown! -- prntid=$(prntid), msgsym $(msgsym), inferdim=$(msg[2])"
      end
      push!(msgspervar[msgsym], msg)
    end
  end

  with_logger(logger) do
    @info "cliq $(prnt.index), prepCliqInitMsgsDown! -- vars fw/ down msgs=$(collect(keys(msgspervar)))"
  end

  # reference to default allocated dict location
  products = getData(prnt).downInitMsg

  ## TODO use parent factors too
  # intersect with the asking clique's seperator variables

    # products only method
    if dbgnew
        condenseDownMsgsProductPrntFactors!(fgl, products, msgspervar, prnt, cliq, logger) # WIP -- not ready yet
    else
        condenseDownMsgsProductOnly!(fgl, products, msgspervar) # BASELINE deprecated
    end

  # remove msgs that have no data
  rmlist = Symbol[]
  for (prsym,belmsg) in products
    if belmsg[2] < 1e-10
      # no information so remove
      push!(rmlist, prsym)
    end
  end
  with_logger(logger) do
    @info "cliq $(prnt.index), prepCliqInitMsgsDown! -- rmlist, no inferdim, keys=$(rmlist)"
  end
  for pr in rmlist
    delete!(products, pr)
  end

  with_logger(logger) do
    @info "cliq $(prnt.index), prepCliqInitMsgsDown! -- product keys=$(collect(keys(products)))"
  end
  return products
end

"""
    $SIGNATURES

Special function to do initialization in downward direction, assuming that not all
variables can be initialized.  Relies on outside down messages.

Notes:
- assumed this `cliq` is being initialized from a previous `:needdownmsg` status.
- will use all possible local factors of cliquq in initilization process
- similar to upward initialization, but uses different message structure
  - first draft assumes upward messages will not be used,
  - full up solve still required which explicitly depends on upward messages.

Dev Notes
- Streamline add/delete msg priors from calling function and csm.
- TODO replace with nested 'minimum degree' type variable ordering.
"""
function getCliqInitVarOrderDown(dfg::G,
                                 cliq::TreeClique,
                                 downmsgs::TempBeliefMsg )::Vector{Symbol} where G <: AbstractDFG
  #
  allsyms = getCliqAllVarIds(cliq)
  # convert input downmsg var symbols to integers (also assumed as prior beliefs)
  # make sure ids are in the clique set, since parent may have more variables.
  dwnmsgsym = intersect(collect(keys(downmsgs)), DFG.getVariableIds(dfg)) #dfg.IDs
  # dwnmsgids =  map(x -> dfg.IDs[x], dwnmsgsym )
  dwnvarids = intersect(allsyms, dwnmsgsym)

  # find any other prior factors (might have partials)
  prvarids = getCliqVarIdsPriors(cliq, allsyms, true)
  hassinglids = union(dwnvarids, prvarids)

  # Get all other variable factor counts
  nfcts = getCliqNumAssocFactorsPerVar(cliq)
  # add msg marginal prior (singletons) to number of factors
  for msid in dwnmsgsym
    nfcts[msid .== allsyms] .+= 1
  end

  # sort permutation order for increasing number of factor association
  nfctsp = sortperm(nfcts)
  sortedids = allsyms[nfctsp]

  # all singleton variables
  singids = union(prvarids, dwnvarids)

  # organize the prior variables separately with asceding factor count
  initorder = Symbol[] #zeros(Int, 0)
  for id in sortedids
    if id in singids
      push!(initorder, id)
    end
  end
  # sort remaining variables for increasing associated factors
  for id in sortedids
    if !(id in initorder)
      push!(initorder, id)
    end
  end

  # return variable order
  return initorder
end

"""
    $SIGNATURES

Modify the `subfg::FactorGraph` to include `msgs` as priors that are used
during clique inference.

Notes
- May be used initialization or inference, in both upward and downward directions.
- `msgs` are identified by variable label `::Symbol`, and my consist of multiple beliefs.
- Message sets from different cliques are identified by clique id `::Int`.

Related

`deleteMsgFactors!`
"""
function addMsgFactors!(subfg::G,
                        msgs::TempBeliefMsg)::Vector{DFGFactor} where G <: AbstractDFG
  # add messages as priors to this sub factor graph
  msgfcts = DFGFactor[]
  svars = DFG.getVariableIds(subfg)
  for (msym, dm) in msgs
    if msym in svars
      # TODO prior missing manifold information
      fc = addFactor!(subfg, [msym], MsgPrior(dm[1], dm[2]), autoinit=false)
      push!(msgfcts, fc)
    end
  end
  return msgfcts
end
function addMsgFactors!(subfg::G, msgs::Dict{Symbol, BallTreeDensity}) where G <: AbstractDFG
  @warn "addMsgFactors! use TempBeliefMsg format instead"

  # assemble new temporary list
  tmpmsgs = TempBeliefMsg()
  for (id, val) in msgs
    tmpmsgs[id] = (val, 0.0)
  end

  # call intended function with temporary message work around
  return addMsgFactors!(subfg, tmpmsgs)
end
function addMsgFactors!(subfg::G,
                        msgs::Dict{Symbol, Vector{Tuple{BallTreeDensity, Float64}}})::Vector{DFGFactor} where G <: AbstractDFG
  # add messages as priors to this sub factor graph
  msgfcts = DFGFactor[]
  svars = ls(subfg)
  for (msym, dms) in msgs
    for dm in dms
      if msym in svars
        # TODO should be on manifold prior, not just generic euclidean prior -- okay since variable on manifold, but not for long term
        fc = addFactor!(subfg, [msym], MsgPrior(dm[1], dm[2]), autoinit=false)
        push!(msgfcts, fc)
      end
    end
  end
  return msgfcts
end
function addMsgFactors!(subfg::G,
                        allmsgs::Dict{Int,TempBeliefMsg} )::Vector{DFGFactor} where G <: AbstractDFG
  #
  allfcts = DFGFactor[]
  for (cliqid, msgs) in allmsgs
    # do each dict in array separately
    newfcts = addMsgFactors!(subfg, msgs)
    union!( allfcts, newfcts )
  end
  return allfcts
end

"""
    $SIGNATURES

Delete from the subgraph`::FactorGraph` prior belief `msgs` that could/would be used
during clique inference.

Related

`addMsgFactors!`
"""
function deleteMsgFactors!(subfg::G,
                           fcts::Vector{DFGFactor}) where G <: AbstractDFG
  #
  for fc in fcts
    deleteFactor!(subfg, fc.label)
  end
end

"""
    $SIGNATURES

Return true or false depending on whether child cliques are all up solved.
"""
function areCliqChildrenAllUpSolved(treel::BayesTree,
                                    prnt::TreeClique)::Bool
  #
  for ch in getChildren(treel, prnt)
    if !isCliqUpSolved(ch)
      return false
    end
  end
  return true
end

"""
    $SIGNATURES

Initialization requires down message passing of more specialized down init msgs.
This function performs any possible initialization of variables and retriggers
children cliques that have not yet initialized.

Notes:
- Assumed this function is only called after status from child clique up inits completed.
- Assumes cliq has parent.
  - will fetch message from parent
- Will perform down initialization if status == `:needdownmsg`.
- might be necessary to pass furhter down messges to child cliques that also `:needdownmsg`.
- Will not complete cliq solve unless all children are `:upsolved` (upward is priority).
- `dwinmsgs` assumed to come from parent initialization process.
- assume `subfg` as a subgraph that can be modified by this function (add message factors)
  - should remove message prior factors from subgraph before returning.
- May modify `cliq` values.
  - `setCliqUpInitMsgs!(cliq, cliq.index, msg)`
  - `setCliqStatus!(cliq, status)`
  - `setCliqDrawColor(cliq, "sienna")`
  - `notifyCliqDownInitStatus!(cliq, status)`

Algorithm:
- determine which downward messages influence initialization order
- initialize from singletons to most connected non-singletons
- revert back to needdownmsg if cycleInit does nothing
- can only ever return :initialized or :needdownmsg status
"""
function doCliqInitDown!(subfg::G,
                         cliq::TreeClique,
                         dwinmsgs::TempBeliefMsg;
                         dbg::Bool=false,
                         logpath::String="/tmp/caesar/",
                         logger=ConsoleLogger() ) where G <: AbstractDFG
  #
  with_logger(logger) do
    @info "cliq $(cliq.index), doCliqInitDown! -- 1, dwinmsgs=$(collect(keys(dwinmsgs)))"
  end
  status = :needdownmsg #:badinit

  # get down variable initialization order
  initorder = getCliqInitVarOrderDown(subfg, cliq, dwinmsgs)
  with_logger(logger) do
    @info "cliq $(cliq.index), doCliqInitDown! -- 4, initorder=$(initorder))"
  end

  # add messages as priors to this sub factor graph
  msgfcts = addMsgFactors!(subfg, dwinmsgs)

  # store the cliqSubFg for later debugging
  if dbg
    DFG.saveDFG(subfg, joinpath(logpath,"logs/cliq$(cliq.index)/fg_beforedowninit"))
  end

  # cycle through vars and attempt init
  with_logger(logger) do
    @info "cliq $(cliq.index), doCliqInitDown! -- 5, cycle through vars and attempt init"
  end
  if cycleInitByVarOrder!(subfg, initorder)
    status = :initialized
  end

  with_logger(logger) do
    @info "cliq $(cliq.index), doCliqInitDown! -- 6, current status: $status"
  end
  # remove msg factors previously added
  deleteMsgFactors!(subfg, msgfcts)

  # store the cliqSubFg for later debugging
  if dbg
      DFG.saveDFG(subfg, joinpath(logpath,"logs/cliq$(cliq.index)/fg_afterdowninit"))
  end

  return status
end

function doCliqInitDown!(subfg::G,
                         tree::BayesTree,
                         cliq::TreeClique;
                         dbg::Bool=false ) where G <: AbstractDFG
  #
  @error "deprecated doCliqInitDown!(subfg, tree, cliq) use doCliqInitDown!(subfg, cliq, dwinmsgs) instead."
  prnt = getParent(tree, cliq)[1]
  dwinmsgs = prepCliqInitMsgsDown!(subfg, tree, prnt)
  status = doCliqInitDown!(subfg, cliq, dwinmsgs, dbg=dbg)

  return status
end

"""
    $SIGNATURES

Return `true` if any of the children cliques have status `:needdownmsg`.
"""
function areCliqChildrenNeedDownMsg(children::Vector{TreeClique})::Bool
  for ch in children
    if getCliqStatus(ch) == :needdownmsg
      return true
    end
  end
  return false
end

function areCliqChildrenNeedDownMsg(tree::BayesTree, cliq::TreeClique)::Bool
  areCliqChildrenNeedDownMsg( getChildren(tree, cliq) )
end

"""
    $SIGNATURES

Return true if has parent with status `:needdownmsg`.
"""
function isCliqParentNeedDownMsg(tree::BayesTree, cliq::TreeClique, logger=ConsoleLogger())
  prnt = getParent(tree, cliq)
  if length(prnt) == 0
    return false
  end
  prstat = getCliqStatus(prnt[1])
  with_logger(logger) do
    @info "$(current_task()) Clique $(cliq.index), isCliqParentNeedDownMsg -- parent status: $(prstat)"
  end
  return prstat == :needdownmsg
end


text:
module PlutoUI

import Base: show, get
import Markdown: htmlesc, withtag

using Reexport

const PKG_ROOT_DIR = normpath(joinpath(@__DIR__, ".."))

@reexport module BuiltinsNotebook
    include("./Builtins.jl")
end
include("./Resource.jl")
include("./DisplayTricks.jl")

@reexport module RangeSliderNotebook
    include("./RangeSlider.jl")
end
@reexport module TerminalNotebook
    include("./TerminalNotebook.jl")
end
@reexport module MultiCheckBoxNotebook
    include("./MultiCheckBox.jl")
end
@reexport module TableOfContentsNotebook
    include("./TableOfContents.jl")
end
@reexport module ClockNotebook
    include("./Clock.jl")
end
@reexport module ScrubbableNotebook
    include("./Scrubbable.jl")
end
@reexport module ConfirmNotebook
    include("./Confirm.jl")
end
module CombineNotebook
    include("./Combine.jl")
end
# not exporting to avoid clash with DataFrames.combine
const combine = CombineNotebook.combine

# this is a submodule
module ExperimentalLayout
    include("./Layout.jl")
end

"""
Components in this module are still experimental, and after a trial period (with breaking changes!), they will be added to PlutoUI as public API. Use with caution.

Currently included: [`PlutoUI.Experimental.transformed_value`](@ref).
"""
module Experimental
module TransformedValueNotebook
    include("./TransformedValue.jl")
end

const transformed_value = TransformedValueNotebook.transformed_value

module WrappedNotebook
    import ...PlutoUI: combine
    import ..Experimental: transformed_value
    include("./Wrapped.jl")
end
const wrapped = WrappedNotebook.wrapped

end
end

ext: jl


text:
# Practical Machine learning
# Reinforcement learning - Q learning example
# Chapter 12

using DeepQLearning
using Base.Test

# write your own tests here
@test 1 == 1


text:
module Counterfactuals

using ..Models, ..Generators, ..DataPreprocessing

include("functions.jl")
export CounterfactualExplanation, initialize!, update!,
    total_steps, converged, terminated, path, target_probs

end


text:
using MCMC
using Test

using BangBang
using Bijectors
using Distributions
using StatsFuns
import Random
import LinearAlgebra

@testset "MCMC" begin
  include("models/simple_model.jl")
  include("models/gibbs_models.jl")
  include("metropolis.jl")
  include("mcmc.jl")
  include("misc.jl")
  include("models/gmm.jl")
  include("metrics.jl")
end


text:
#
# VUMPS code
#

struct Hᶜ
  ∑h::InfiniteSum{MPO}
  Hᴸ::InfiniteMPS
  Hᴿ::InfiniteMPS
  ψ::InfiniteCanonicalMPS
  n::Int
end

struct Hᴬᶜ
  ∑h::InfiniteSum{MPO}
  Hᴸ::InfiniteMPS
  Hᴿ::InfiniteMPS
  ψ::InfiniteCanonicalMPS
  n::Int
end

#Assume all local Hamiltonians have the same nrange
function Base.:*(H::Hᶜ, v::ITensor)
  ∑h = H.∑h
  Hᴸ = H.Hᴸ
  Hᴿ = H.Hᴿ
  ψ = H.ψ
  ψ′ = dag(ψ)'
  Nsites = nsites(ψ)
  range_∑h = nrange(∑h, 1)
  n = H.n
  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)
  s = siteinds(only, ψ)
  s′ = siteinds(only, ψ′)
  δˢ(n) = δ(dag(s[n]), prime(s[n]))
  δʳ(n) = δ(dag(r[n]), prime(r[n]))
  δˡ(n) = δ(l[n], l′[n])

  #Build the contribution of the left environment
  Hᶜᴸv = v * Hᴸ[n] * δʳ(n)
  #Build the contribution of the right environment
  Hᶜᴿv = v * δˡ(n) * Hᴿ[n]
  #We now start building terms where C overlap with the local Hamiltonian
  # We start with the tensor AL[n] - v - AR[n+1] ... AR[n + range_∑h - 1]
  Hᶜʰv = v * ψ.AL[n] * δˡ(n - 1) * ψ′.AL[n] * ∑h[n][1] #left extremity
  common_sites = findsites(ψ, ∑h[n])
  idx = 2 #list the sites Σh, we start at 2 because n is already taken into account
  for k in 1:(range_∑h - 2)
    if n + k == common_sites[idx]
      Hᶜʰv = Hᶜʰv * ψ.AR[n + k] * ψ′.AR[n + k] * ∑h[n][idx]
      idx += 1
    else
      Hᶜʰv = Hᶜʰv * ψ.AR[n + k] * ψ′.AR[n + k] * δˢ(n + k)
    end
  end
  Hᶜʰv =
    Hᶜʰv *
    ψ.AR[n + range_∑h - 1] *
    δʳ(n + range_∑h - 1) *
    ψ′.AR[n + range_∑h - 1] *
    ∑h[n][end]    #right most extremity
  #Now we are building contributions of the form AL[n - j] ... AL[n] - v - AR[n + 1] ... AR[n + range_∑h - 1 - j]
  for j in 1:(range_∑h - 2)
    temp_Hᶜʰv = ψ.AL[n - j] * δˡ(n - 1 - j) * ψ′.AL[n - j] * ∑h[n - j][1]
    common_sites = findsites(ψ, ∑h[n - j])
    idx = 2
    for k in 1:j
      if n - j + k == common_sites[idx]
        temp_Hᶜʰv = temp_Hᶜʰv * ψ.AL[n - j + k] * ψ′.AL[n - j + k] * ∑h[n - j][idx]
        idx += 1
      else
        temp_Hᶜʰv = temp_Hᶜʰv * ψ.AL[n - j + k] * ψ′.AL[n - j + k] * δˢ(n - j + k)
      end
    end
    # Finished the AL part
    temp_Hᶜʰv = temp_Hᶜʰv * v
    for k in (j + 1):(range_∑h - 2)
      if n - j + k == common_sites[idx]
        temp_Hᶜʰv = temp_Hᶜʰv * ψ.AR[n - j + k] * ψ′.AR[n - j + k] * ∑h[n - j][idx]
        idx += 1
      else
        temp_Hᶜʰv = temp_Hᶜʰv * ψ.AR[n - j + k] * ψ′.AR[n - j + k] * δˢ(n - j + k)
      end
    end
    temp_Hᶜʰv =
      temp_Hᶜʰv *
      ∑h[n - j][end] *
      ψ.AR[n - j + range_∑h - 1] *
      δʳ(n - j + range_∑h - 1) *
      ψ′.AR[n - j + range_∑h - 1]
    Hᶜʰv = Hᶜʰv + temp_Hᶜʰv
  end
  Hᶜv = Hᶜᴸv + Hᶜʰv + Hᶜᴿv
  return Hᶜv * dag(δˡ(n)) * dag(δʳ(n))
end

function Base.:*(H::Hᴬᶜ, v::ITensor)
  ∑h = H.∑h
  Hᴸ = H.Hᴸ
  Hᴿ = H.Hᴿ
  ψ = H.ψ
  ψ′ = dag(ψ)'
  Nsites = nsites(ψ)
  range_∑h = nrange(∑h, 1)
  n = H.n
  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)
  s = siteinds(only, ψ)
  s′ = siteinds(only, ψ′)

  δˢ(n) = δ(dag(s[n]), prime(s[n]))
  δʳ(n) = δ(dag(r[n]), prime(r[n]))
  δˡ(n) = δ(l[n], l′[n])

  #Contribution of the left environment
  Hᴬᶜᴸv = v * Hᴸ[n - 1] * δˢ(n) * δʳ(n)
  #Contribution of the right environment
  Hᴬᶜᴿv = v * δˡ(n - 1) * δˢ(n) * Hᴿ[n]
  #We now start building terms where AC overlap with the local Hamiltonian
  # We start with the tensor v - AR[n+1] ... AR[n + range_∑h - 1]
  Hᴬᶜʰv = v * δˡ(n - 1) * ∑h[n][1]
  common_sites = findsites(ψ, ∑h[n])
  idx = 2 #list the sites Σh, we start at 2 because n is already taken into account
  for k in 1:(range_∑h - 2)
    if n + k == common_sites[idx]
      Hᴬᶜʰv = Hᴬᶜʰv * ψ.AR[n + k] * ψ′.AR[n + k] * ∑h[n][idx]
      idx += 1
    else
      Hᴬᶜʰv = Hᴬᶜʰv * ψ.AR[n + k] * ψ′.AR[n + k] * δˢ(n + k)
    end
  end
  Hᴬᶜʰv =
    Hᴬᶜʰv *
    ∑h[n][end] *
    ψ.AR[n + range_∑h - 1] *
    ψ′.AR[n + range_∑h - 1] *
    δʳ(n + range_∑h - 1) #rightmost extremity
  #Now we are building contributions of the form AL[n - j] ... AL[n-1] - v - AR[n + 1] ... AR[n + range_∑h - 1 - j]
  for j in 1:(range_∑h - 1)
    temp_Hᴬᶜʰv = ψ.AL[n - j] * δˡ(n - j - 1) * ψ′.AL[n - j] * ∑h[n - j][1]
    common_sites = findsites(ψ, ∑h[n - j])
    idx = 2
    for k in 1:(j - 1)
      if n - j + k == common_sites[idx]
        temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * ψ.AL[n - j + k] * ψ′.AL[n - j + k] * ∑h[n - j][idx]
        idx += 1
      else
        temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * ψ.AL[n - j + k] * ψ′.AL[n - j + k] * δˢ(n - j + k)
      end
    end
    #Finished with AL, treating the center AC = v
    if j == range_∑h - 1
      temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * v * δʳ(n - j + range_∑h - 1) * ∑h[n - j][end]
    else
      if n == common_sites[idx] #need to check whether we need to branch v
        temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * v * ∑h[n - j][idx]
        idx += 1
      else
        temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * v * δˢ(n)
      end
      for k in (j + 1):(range_∑h - 2)
        if n + k - j == common_sites[idx]
          temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * ψ.AR[n + k - j] * ψ′.AR[n + k - j] * ∑h[n - j][idx]
          idx += 1
        else
          temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * ψ.AR[n + k - j] * ψ′.AR[n + k - j] * δˢ(n + k - j)
        end
      end
      temp_Hᴬᶜʰv =
        temp_Hᴬᶜʰv *
        ∑h[n - j][end] *
        ψ.AR[n + range_∑h - 1 - j] *
        ψ′.AR[n + range_∑h - 1 - j] *
        δʳ(n - j + range_∑h - 1)
    end
    Hᴬᶜʰv = Hᴬᶜʰv + temp_Hᴬᶜʰv
  end
  Hᴬᶜv = Hᴬᶜᴸv + Hᴬᶜʰv + Hᴬᶜᴿv
  return Hᴬᶜv * dag(δˡ(n - 1)) * dag(δˢ(n)) * dag(δʳ(n))
end

function (H::Hᶜ)(v)
  Hᶜv = H * v
  ## return Hᶜv * δˡ * δʳ
  return noprime(Hᶜv)
end

function (H::Hᴬᶜ)(v)
  Hᴬᶜv = H * v
  ## return Hᶜv * δˡ⁻¹ * δˢ * δʳ
  return noprime(Hᴬᶜv)
end

# Struct for use in linear system solver
struct Aᴸ
  ψ::InfiniteCanonicalMPS
  n::Int
end

function (A::Aᴸ)(x)
  ψ = A.ψ
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  ψ̃ = prime(linkinds, ψᴴ)
  n = A.n

  N = length(ψ)
  #@assert n == N

  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)

  xT = translatecell(x, -1)
  for k in (n - N + 1):n
    xT = xT * ψ.AL[k] * ψ̃.AL[k]
  end
  δˡ = δ(l[n], l′[n])
  δʳ = δ(r[n], r′[n])
  xR = x * ψ.C[n] * ψ′.C[n] * dag(δʳ) * denseblocks(δˡ)
  return xT - xR
end

# TODO Generate all environments, why? Only one is needed in the sequential version
function left_environment(hᴸ, 𝕙ᴸ, ψ; tol=1e-15)
  ψ̃ = prime(linkinds, dag(ψ))
  N = nsites(ψ)

  Aᴺ = Aᴸ(ψ, N)
  Hᴸᴺ¹, info = linsolve(Aᴺ, 𝕙ᴸ[N], 1, -1; tol=tol)
  # Get the rest of the environments in the unit cell
  Hᴸ = InfiniteMPS(Vector{ITensor}(undef, N))
  Hᴸ[N] = Hᴸᴺ¹
  for n in 1:(N - 1)
    Hᴸ[n] = Hᴸ[n - 1] * ψ.AL[n] * ψ̃.AL[n] + hᴸ[n]
  end
  return Hᴸ
end

function left_environment(∑h::InfiniteSum{MPO}, ψ::InfiniteCanonicalMPS; tol=1e-15)
  Nsites = nsites(ψ)
  range_∑h = nrange(∑h, 1)
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  ψ̃ = prime(linkinds, ψᴴ)

  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)
  s = siteinds(only, ψ)
  δʳ(n) = δ(dag(r[n]), prime(r[n]))
  δˡ(n) = δ(l[n], l′[n])
  δˢ(n) = δ(dag(s[n]), prime(s[n]))
  hᴸ = Vector{ITensor}(undef, Nsites)
  for k in 1:Nsites
    hᴸ[k] =
      δˡ(k - range_∑h) *
      ψ.AL[k - range_∑h + 1] *
      ∑h[(k - range_∑h + 1, k - range_∑h + 2)][1] *
      ψ′.AL[k - range_∑h + 1]
    common_sites = findsites(ψ, ∑h[(k - range_∑h + 1, k - range_∑h + 2)])
    idx = 2
    for j in 2:range_∑h
      if k - range_∑h + j == common_sites[idx]
        hᴸ[k] =
          hᴸ[k] *
          ψ.AL[k - range_∑h + j] *
          ψ′.AL[k - range_∑h + j] *
          ∑h[(k - range_∑h + 1, k - range_∑h + 2)][idx]
        idx += 1
      else
        hᴸ[k] =
          hᴸ[k] * ψ.AL[k - range_∑h + j] * ψ′.AL[k - range_∑h + j] * δˢ(k - range_∑h + j)
      end
    end
  end
  hᴸ = InfiniteMPS(hᴸ)
  eᴸ = [(hᴸ[k] * ψ.C[k] * δʳ(k) * ψ′.C[k])[] for k in 1:Nsites]
  for k in 1:Nsites
    # TODO: remove `denseblocks` once BlockSparse + DiagBlockSparse is supported
    hᴸ[k] -= eᴸ[k] * denseblocks(δ(inds(hᴸ[k])))
  end

  𝕙ᴸ = copy(hᴸ)
  # TODO restrict to the useful ones only?
  for n in 1:Nsites
    for k in 1:(Nsites - 1)
      temp = copy(hᴸ[n - k])
      for kp in reverse(0:(k - 1))
        temp = temp * ψ.AL[n - kp] * ψ̃.AL[n - kp]
      end
      𝕙ᴸ[n] = temp + 𝕙ᴸ[n]
    end
  end
  Hᴸ = left_environment(hᴸ, 𝕙ᴸ, ψ; tol=tol)

  return Hᴸ, eᴸ
end
# Struct for use in linear system solver
struct Aᴿ
  hᴿ::InfiniteMPS
  ψ::InfiniteCanonicalMPS
  n::Int
end

function (A::Aᴿ)(x)
  hᴿ = A.hᴿ
  ψ = A.ψ
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  ψ̃ = prime(linkinds, ψᴴ)
  n = A.n

  N = length(ψ)
  @assert n == N

  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)

  xT = x
  for k in reverse(1:N)
    xT = xT * ψ.AR[k] * ψ̃.AR[k]
  end
  xT = translatecell(xT, 1)
  δˡ = δ(l[n], l′[n])
  δʳ = δ(r[n], r′[n])
  xR = x * ψ.C[n] * ψ′.C[n] * δˡ * denseblocks(dag(δʳ))
  return xT - xR
end

# TODO Generate all environments, why? Only one is needed in the sequential version
function right_environment(hᴿ, 𝕙ᴿ, ψ; tol=1e-15)
  ψ̃ = prime(linkinds, dag(ψ))
  N = nsites(ψ)

  A = Aᴿ(hᴿ, ψ, N)
  Hᴿᴺ¹, info = linsolve(A, 𝕙ᴿ[N], 1, -1; tol=tol)
  # Get the rest of the environments in the unit cell
  Hᴿ = InfiniteMPS(Vector{ITensor}(undef, N))
  Hᴿ[N] = Hᴿᴺ¹
  for n in reverse(1:(N - 1))
    Hᴿ[n] = Hᴿ[n + 1] * ψ.AR[n + 1] * ψ̃.AR[n + 1] + hᴿ[n]
  end
  return Hᴿ
end

function right_environment(∑h::InfiniteSum{MPO}, ψ::InfiniteCanonicalMPS; tol=1e-15)
  Nsites = nsites(ψ)
  range_∑h = nrange(∑h, 1)
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  ψ̃ = prime(linkinds, ψᴴ)

  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)
  s = siteinds(only, ψ)
  δʳ(n) = δ(dag(r[n]), prime(r[n]))
  δˡ(n) = δ(l[n], l′[n])
  δˢ(n) = δ(dag(s[n]), prime(s[n]))

  hᴿ = Vector{ITensor}(undef, Nsites)
  for k in 1:Nsites
    hᴿ[k] = ψ.AR[k + range_∑h] * ∑h[k + 1][end] * ψ′.AR[k + range_∑h] * δʳ(k + range_∑h)
    common_sites = findsites(ψ, ∑h[k + 1])
    idx = length(common_sites) - 1
    for j in (range_∑h - 1):-1:1
      if k + j == common_sites[idx]
        hᴿ[k] = hᴿ[k] * ψ.AR[k + j] * ψ′.AR[k + j] * ∑h[k + 1][idx]
        idx -= 1
      else
        hᴿ[k] = hᴿ[k] * ψ.AR[k + j] * ψ′.AR[k + j] * δˢ(k + j)
      end
    end
  end
  hᴿ = InfiniteMPS(hᴿ)
  eᴿ = [(hᴿ[k] * ψ.C[k] * δˡ(k) * ψ′.C[k])[] for k in 1:Nsites]
  for k in 1:Nsites
    hᴿ[k] -= eᴿ[k] * denseblocks(δ(inds(hᴿ[k])))
  end

  𝕙ᴿ = copy(hᴿ)
  #TODO restrict to the useful ones only
  for n in 1:Nsites
    for k in 1:(Nsites - 1)
      temp = copy(hᴿ[n + k])
      for kp in reverse(1:k)
        temp = temp * ψ.AR[n + kp] * ψ̃.AR[n + kp]
      end
      𝕙ᴿ[n] = temp + 𝕙ᴿ[n]
    end
  end
  Hᴿ = right_environment(hᴿ, 𝕙ᴿ, ψ; tol=tol)
  return Hᴿ, eᴿ
end

function tdvp_iteration(args...; multisite_update_alg="sequential", kwargs...)
  if multisite_update_alg == "sequential"
    return tdvp_iteration_sequential(args...; kwargs...)
  elseif multisite_update_alg == "parallel"
    return tdvp_iteration_parallel(args...; kwargs...)
  else
    error(
      "Multisite update algorithm multisite_update_alg = $multisite_update_alg not supported, use \"parallel\" or \"sequential\"",
    )
  end
end

#TODO put these functions somewhere else
function ortho_overlap(AC, C)
  AL, _ = polar(AC * dag(C), uniqueinds(AC, C))
  return noprime(AL)
end

function ortho_polar(AC, C)
  UAC, _ = polar(AC, uniqueinds(AC, C))
  UC, _ = polar(C, commoninds(C, AC))
  return noprime(UAC) * noprime(dag(UC))
end

#In principle, could share even more code with vumps_mpo or with parallel
function tdvp_iteration_sequential(
  solver::Function,
  ∑h::InfiniteSum,
  ψ::InfiniteCanonicalMPS;
  (ϵᴸ!)=fill(1e-15, nsites(ψ)),
  (ϵᴿ!)=fill(1e-15, nsites(ψ)),
  time_step,
  solver_tol=(x -> x / 100),
)
  Nsites = nsites(ψ)
  range_∑h = nrange(∑h, 1)
  ϵᵖʳᵉˢ = max(maximum(ϵᴸ!), maximum(ϵᴿ!))
  _solver_tol = solver_tol(ϵᵖʳᵉˢ)

  ψ = copy(ψ)
  C̃ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  Ãᶜ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  Ãᴸ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  Ãᴿ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  eᴸ = Vector{Float64}(undef, Nsites)
  eᴿ = Vector{Float64}(undef, Nsites)

  for n in 1:Nsites
    ψᴴ = dag(ψ)
    ψ′ = ψᴴ'
    # XXX: make this prime the center sites
    ψ̃ = prime(linkinds, ψᴴ)

    l = linkinds(only, ψ.AL)
    l′ = linkinds(only, ψ′.AL)
    r = linkinds(only, ψ.AR)
    r′ = linkinds(only, ψ′.AR)
    s = siteinds(only, ψ)
    δʳ(n) = δ(dag(r[n]), prime(r[n]))
    δˡ(n) = δ(l[n], l′[n])
    δˢ(n) = δ(dag(s[n]), prime(s[n]))

    Hᴸ, eᴸ = left_environment(∑h, ψ; tol=_solver_tol)
    Hᴿ, eᴿ = right_environment(∑h, ψ; tol=_solver_tol)

    Cvalsₙ₋₁, Cvecsₙ₋₁, Cinfoₙ₋₁ = solver(
      Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n - 1), time_step, ψ.C[n - 1], _solver_tol
    )
    Cvalsₙ, Cvecsₙ, Cinfoₙ = solver(Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.C[n], _solver_tol)
    Avalsₙ, Avecsₙ, Ainfoₙ = solver(
      Hᴬᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.AL[n] * ψ.C[n], _solver_tol
    )

    C̃[n - 1] = Cvecsₙ₋₁
    C̃[n] = Cvecsₙ
    Ãᶜ[n] = Avecsₙ

    Ãᴸ[n] = ortho_polar(Ãᶜ[n], C̃[n])
    Ãᴿ[n] = ortho_polar(Ãᶜ[n], C̃[n - 1])
    # Update state for next iteration
    #ψ = InfiniteCanonicalMPS(Ãᴸ, C̃, Ãᴿ)
    ψ.AL[n] = Ãᴸ[n]
    ψ.AR[n] = Ãᴿ[n]
    ψ.C[n - 1] = C̃[n - 1]
    ψ.C[n] = C̃[n]
  end

  for n in 1:Nsites
    ϵᴸ![n] = norm(Ãᶜ[n] - Ãᴸ[n] * C̃[n])
    ϵᴿ![n] = norm(Ãᶜ[n] - C̃[n - 1] * Ãᴿ[n])
  end
  return ψ, (eᴸ, eᴿ)
end

function tdvp_iteration_parallel(
  solver::Function,
  ∑h::InfiniteSum,
  ψ::InfiniteCanonicalMPS;
  (ϵᴸ!)=fill(1e-15, nsites(ψ)),
  (ϵᴿ!)=fill(1e-15, nsites(ψ)),
  time_step,
  solver_tol=(x -> x / 100),
)
  Nsites = nsites(ψ)
  range_∑h = nrange(∑h, 1)
  ϵᵖʳᵉˢ = max(maximum(ϵᴸ!), maximum(ϵᴿ!))
  _solver_tol = solver_tol(ϵᵖʳᵉˢ)
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  # XXX: make this prime the center sites
  ψ̃ = prime(linkinds, ψᴴ)

  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)
  s = siteinds(only, ψ)
  δʳ(n) = δ(dag(r[n]), prime(r[n]))
  δˡ(n) = δ(l[n], l′[n])
  δˢ(n) = δ(dag(s[n]), prime(s[n]))
  Hᴸ, eᴸ = left_environment(∑h, ψ; tol=_solver_tol)
  Hᴿ, eᴿ = right_environment(∑h, ψ; tol=_solver_tol)

  C̃ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  Ãᶜ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  for n in 1:Nsites
    Cvalsₙ, Cvecsₙ, Cinfoₙ = solver(Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.C[n], _solver_tol)
    Avalsₙ, Avecsₙ, Ainfoₙ = solver(
      Hᴬᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.AL[n] * ψ.C[n], _solver_tol
    )

    C̃[n] = Cvecsₙ
    Ãᶜ[n] = Avecsₙ
  end

  Ãᴸ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  Ãᴿ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  for n in 1:Nsites
    Ãᴸ[n] = ortho_polar(Ãᶜ[n], C̃[n])
    Ãᴿ[n] = ortho_polar(Ãᶜ[n], C̃[n - 1])
  end

  for n in 1:Nsites
    ϵᴸ![n] = norm(Ãᶜ[n] - Ãᴸ[n] * C̃[n])
    ϵᴿ![n] = norm(Ãᶜ[n] - C̃[n - 1] * Ãᴿ[n])
  end
  return InfiniteCanonicalMPS(Ãᴸ, C̃, Ãᴿ), (eᴸ, eᴿ)
end

function tdvp(
  solver::Function,
  ∑h,
  ψ;
  maxiter=10,
  tol=1e-8,
  outputlevel=1,
  multisite_update_alg="sequential",
  time_step,
  solver_tol=(x -> x / 100),
)
  N = nsites(ψ)
  (ϵᴸ!) = fill(tol, nsites(ψ))
  (ϵᴿ!) = fill(tol, nsites(ψ))
  outputlevel > 0 &&
    println("Running VUMPS with multisite_update_alg = $multisite_update_alg")
  for iter in 1:maxiter
    ψ, (eᴸ, eᴿ) = tdvp_iteration(
      solver,
      ∑h,
      ψ;
      (ϵᴸ!)=(ϵᴸ!),
      (ϵᴿ!)=(ϵᴿ!),
      multisite_update_alg=multisite_update_alg,
      solver_tol=solver_tol,
      time_step=time_step,
    )
    ϵᵖʳᵉˢ = max(maximum(ϵᴸ!), maximum(ϵᴿ!))
    maxdimψ = maxlinkdim(ψ[0:(N + 1)])
    outputlevel > 0 && println(
      "VUMPS iteration $iter (out of maximum $maxiter). Bond dimension = $maxdimψ, energy = $((eᴸ, eᴿ)), ϵᵖʳᵉˢ = $ϵᵖʳᵉˢ, tol = $tol",
    )
    if ϵᵖʳᵉˢ < tol
      println(
        "Precision error $ϵᵖʳᵉˢ reached tolerance $tol, stopping VUMPS after $iter iterations (of a maximum $maxiter).",
      )
      break
    end
  end
  return ψ
end

function vumps_solver(M, time_step, v₀, solver_tol)
  λ⃗, v⃗, info = eigsolve(M, v₀, 1, :SR; ishermitian=true, tol=solver_tol)
  return λ⃗[1], v⃗[1], info
end

return function tdvp_solver(M, time_step, v₀, solver_tol)
  v, info = exponentiate(M, time_step, v₀; ishermitian=true, tol=solver_tol)
  v = v / norm(v)
  return nothing, v, info
end

function vumps(
  args...; time_step=-Inf, eigsolve_tol=(x -> x / 100), solver_tol=eigsolve_tol, kwargs...
)
  @assert isinf(time_step) && time_step < 0
  println("Using VUMPS solver with time step $time_step")
  return tdvp(vumps_solver, args...; time_step=time_step, solver_tol=solver_tol, kwargs...)
end

function tdvp(args...; time_step, solver_tol=(x -> x / 100), kwargs...)
  solver = if !isinf(time_step)
    println("Using TDVP solver with time step $time_step")
    tdvp_solver
  elseif time_step < 0
    # Call VUMPS instead
    println("Using VUMPS solver with time step $time_step")
    vumps_solver
  else
    error("Time step $time_step not supported.")
  end
  return tdvp(solver, args...; time_step=time_step, solver_tol=solver_tol, kwargs...)
end

##################################################################
# Old functionality, only used for testing

## function left_environment_recursive(hᴸ, ψ; niter=10)
##   ψ̃ = prime(linkinds, dag(ψ))
##   # XXX: replace with `nsites`
##   #N = nsites(ψ)
##   N = length(ψ)
##   Hᴸᴺ¹ = hᴸ[N]
##   for _ in 1:niter
##     Hᴸᴺ¹ = translatecell(Hᴸᴺ¹, -1)
##     for n in 1:N
##       Hᴸᴺ¹ = Hᴸᴺ¹ * ψ.AL[n] * ψ̃.AL[n]
##     end
##     # Loop over the Hamiltonian terms in the unit cell
##     for n in 1:N
##       hᴸⁿ = hᴸ[n]
##       for k in (n + 1):N
##         hᴸⁿ = hᴸⁿ * ψ.AL[k] * ψ̃.AL[k]
##       end
##       Hᴸᴺ¹ += hᴸⁿ
##     end
##   end
##   # Get the rest of the environments in the unit cell
##   Hᴸ = InfiniteMPS(Vector{ITensor}(undef, N))
##   Hᴸ[N] = Hᴸᴺ¹
##   Hᴸᴺ¹ = translatecell(Hᴸᴺ¹, -1)
##   for n in 1:(N - 1)
##     Hᴸ[n] = Hᴸ[n - 1] * ψ.AL[n] * ψ̃.AL[n] + hᴸ[n]
##   end
##   return Hᴸ
## end
##
## function right_environment_recursive(hᴿ, ψ; niter=10)
##   ψ̃ = prime(linkinds, dag(ψ))
##   # XXX: replace with `nsites`
##   #N = nsites(ψ)
##   N = length(ψ)
##   Hᴿᴺ¹ = hᴿ[0]
##   for _ in 1:niter
##     Hᴿᴺ¹ = translatecell(Hᴿᴺ¹, 1)
##     for n in reverse(1:N)
##       Hᴿᴺ¹ = Hᴿᴺ¹ * ψ.AR[n] * ψ̃.AR[n]
##     end
##     # Loop over the Hamiltonian terms in the unit cell
##     for n in reverse(0:(N - 1))
##       hᴿⁿ = hᴿ[n]
##       for k in reverse(1:n)
##         hᴿⁿ = hᴿⁿ * ψ.AR[k] * ψ̃.AR[k]
##       end
##       Hᴿᴺ¹ += hᴿⁿ
##     end
##   end
##   Hᴿᴺ¹ = translatecell(Hᴿᴺ¹, 1)
##   # Get the rest of the environments in the unit cell
##   Hᴿ = InfiniteMPS(Vector{ITensor}(undef, N))
##   Hᴿ[N] = Hᴿᴺ¹
##   for n in reverse(1:(N - 1))
##     Hᴿ[n] = Hᴿ[n + 1] * ψ.AR[n + 1] * ψ̃.AR[n + 1] + hᴿ[n]
##   end
##   return Hᴿ
## end


text:
using NativeDebugger
"""
Attempt to emulate the first instruction in insts. Return value indicates whether
the instruction was known.
"""
function instemulate!(insts, vm, RC)
    if insts[1] == 0x55 # pushq %rbp
        rsp = get_dwarf(RC, inverse_dwarf[:rsp])
        set_dwarf!(RC, inverse_dwarf[:rsp], rsp-8)
        NativeDebugger.store!(vm, Gallium.RemotePtr{UInt64}(rsp-8),
            get_dwarf(RC, inverse_dwarf[:rbp]))
        set_ip!(RC, UInt(Gallium.ip(RC))+1)
        return true
    end
    @show insts
    return false
end


text:
module Findpeaks

export findpeaks

"""
`findpeaks(y::Array{T},
x::Array{S}=collect(1:length(y))
;minHeight::T=minimum(y), minProm::T=minimum(y),
minDist::S=0, threshold::T=0 ) where {T<:Real,S}`\n
Returns indices of local maxima (sorted from highest peaks to lowest)
in 1D array of real numbers. Similar to MATLAB's findpeaks().\n
*Arguments*:\n
`y` -- data\n
*Optional*:\n
`x` -- x-data\n
*Keyword*:\n
`minHeight` -- minimal peak height\n
`minProm` -- minimal peak prominence\n
`minDist` -- minimal peak distance (keeping highest peaks)\n
`threshold` -- minimal difference (absolute value) between
 peak and neighboring points\n
"""
function findpeaks(y::Array{T},
    x::Array{S}=collect(1:length(y))
    ;minHeight::T=minimum(y), minProm::T=minimum(y),
    minDist::S=zero(x[1]), threshold::T=zero(y[1]) ) where {T<:Real,S}

    peaks = Array{Int64,1}()
    dy = diff(y)
    #differences inside threshold
    for i = 2:length(dy)
        if dy[i] <= -threshold && dy[i-1] >= threshold
            push!(peaks,i)
        end
    end

    #calculate prominences
    yP = y[peaks]
    prominence = zeros(Float64,length(yP))
    for (i,p) in enumerate(peaks)
        lP, rP = 1, length(y)
        for j = (i-1):-1:1
            if yP[j] > yP[i]
                lP = peaks[j]
                break
            end
        end
        ml = minimum(y[lP:p])
        for j = (i+1):length(yP)
            if yP[j] > yP[i]
                rP = peaks[j]
                break
            end
        end
        mr = minimum(y[p:rP])
        ref = max(mr,ml)
        prominence[i] = yP[i] - ref
    end

    #minimal prominence refinement
    peaks = peaks[prominence .> minProm]

    #minimal height refinement
    peaks = peaks[y[peaks] .> minHeight]
    yP = y[peaks]

    #minimal distance refinement
    peaks2del = zeros(Bool,length(peaks))
    inds = sortperm(y[peaks],rev=true)
    permute!(peaks,inds)
    for i = 1:length(peaks)
        for j = 1:(i-1)
            if abs(x[peaks[i]] - x[peaks[j]]) <= minDist
                if !peaks2del[j]
                    peaks2del[i] = true
                end
            end
        end
    end

    peaks[.!peaks2del]
end

end # module


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# test/callback.jl
# Testing callbacks
# Must be run as part of runtests.jl, as it needs a list of solvers.
#############################################################################
using JuMP, MathProgBase, FactCheck

facts("[callback] Test lazy constraints") do
for lazysolver in lazy_solvers
context("With solver $(typeof(lazysolver))") do
    entered = [false,false]

    mod = Model(solver=lazysolver)
    @variable(mod, 0 <= x <= 2, Int)
    @variable(mod, 0 <= y <= 2, Int)
    @objective(mod, Max, y + 0.5x)
    function corners(cb)
        x_val = getvalue(x)
        y_val = getvalue(y)
        TOL = 1e-6
        # Check top right
        if y_val + x_val > 3 + TOL
            @lazyconstraint(cb, y + 0.5x + 0.5x <= 3)
        end
        entered[1] = true
        @fact_throws ErrorException @variable(cb, z)
        @fact_throws ErrorException @lazyconstraint(cb, x^2 <= 1)
    end
    addlazycallback(mod, corners)
    addlazycallback(mod, cb -> (entered[2] = true))
    @fact solve(mod) --> :Optimal
    @fact entered --> [true,true]
    @fact getvalue(x) --> roughly(1.0, 1e-6)
    @fact getvalue(y) --> roughly(2.0, 1e-6)
end; end; end

facts("[callback] Test vectorized lazy constraints") do
for lazysolver in lazy_solvers
context("With solver $(typeof(lazysolver))") do
    entered = [false,false]

    mod = Model(solver=lazysolver)
    @variable(mod, 0 <= x <= 2, Int)
    @variable(mod, 0 <= y <= 2, Int)
    @objective(mod, Max, y + 0.5x)
    function corners(cb)
        x_val = getvalue(x)
        y_val = getvalue(y)
        TOL = 1e-6
        # Check top right
        if y_val + x_val > 3 + TOL
            A = [1.0 1.0
                 1.0 0.5]
            b = [3, 4]
            @lazyconstraint(cb, A*[x,y] .<= b)
        end
        entered[1] = true
        @fact_throws ErrorException @variable(cb, z)
        @fact_throws ErrorException @lazyconstraint(cb, x^2 <= 1)
    end
    addlazycallback(mod, corners)
    addlazycallback(mod, cb -> (entered[2] = true))
    @fact solve(mod) --> :Optimal
    @fact entered --> [true,true]
    @fact getvalue(x) --> roughly(1.0, 1e-6)
    @fact getvalue(y) --> roughly(2.0, 1e-6)
end; end; end

facts("[callback] Test local lazy constraints") do
for lazylocalsolver in lazylocal_solvers
context("With solver $(typeof(lazylocalsolver))") do
    entered = [false,false, false]

    weights = [24714888; 21118272; 5063487; 23450813; 5598179; 4049178; 13516450; 8385365; 9684076; 31317634; 14084148; 21750211; 29261668; 17996589; 12115244]
    values = weights
    W = floor(sum(weights)/2) # knsapsack instance generated according to criteria in 'Hard knapsack Instances' (Chvatal, Op. Res., 1980), with weights composed of random integers in [1:10^(nbItems / 2)]
    mod = Model(solver=lazylocalsolver)
    @variable(mod, 0 <= x[1:length(weights)] <= 1, Int)
    @objective(mod, Max, dot(x, values))
    @constraint(mod, dot(x, weights) <=  W)

    function mycb_localzero(cb)
        nodesexpl = MathProgBase.cbgetexplorednodes(cb)
        if  entered[1] == false && nodesexpl >= 1
            # the following lazy cut  constrains all x[i] to be zero, but applies only locally at the node of the first feasible solution found: it doesn't preclude the existence of "optimal" non-trival solutions
            @lazyconstraint(cb, sum{x[i], i=1:length(weights)} <= 0, localcut=true)
            # @lazyconstraint(cb, sum{x[i], i=1:length(weights)} <= 0) # applying the cut globally would lead the solver to x=0 as the optimal solution
            @fact macroexpand(:(@lazyconstraint(cb, sum{x[i], i=1:length(weights)} <= 0, badkwarg=true))).head --> :error
            entered[1] = true
        end
        entered[2] = true
    end
    addlazycallback(mod, mycb_localzero)
    addlazycallback(mod, cb -> (entered[3] = true))
    @fact solve(mod) --> :Optimal
    @fact entered --> [true,true,true]
    @fact sum(getvalue(x)) --> greater_than(0)
end; end; end

facts("[callback] Test user cuts") do
for cutsolver in cut_solvers
context("With solver $(typeof(cutsolver))") do
    entered = [false,false]

    N = 1000
    # Include explicit data from srand(234) so that we can reproduce across platforms
    include(joinpath("data","usercut.jl"))
    mod = Model(solver=cutsolver)
    @variable(mod, x[1:N], Bin)
    @objective(mod, Max, dot(r1,x))
    @constraint(mod, c[i=1:10], dot(r2[i],x) <= rhs[i]*N/10)
    function mycutgenerator(cb)
        # add a trivially valid cut
        @usercut(cb, sum{x[i], i=1:N} <= N)
        entered[1] = true
    end
    addcutcallback(mod, mycutgenerator)
    addcutcallback(mod, cb -> (entered[2] = true))
    @fact solve(mod) --> :Optimal
    @fact entered --> [true,true]
    @fact find(getvalue(x)[:]) --> [35,38,283,305,359,397,419,426,442,453,526,553,659,751,840,865,878,978]
end; end; end

facts("[callback] Test local user cuts") do
for cutlocalsolver in cutlocal_solvers
context("With solver $(typeof(cutlocalsolver))") do
    entered = [false,false, false]

    weights = [24714888; 21118272; 5063487; 23450813; 5598179; 4049178; 13516450; 8385365; 9684076; 31317634; 14084148; 21750211; 29261668; 17996589; 12115244]
    values = weights
    W = floor(sum(weights)/2) # knsapsack instance generated according to criteria in 'Hard knapsack Instances' (Chvatal, Op. Res., 1980), with weights composed of random integers in [1:10^(nbItems / 2)]
    mod = Model(solver=cutlocalsolver)
    @variable(mod, 0 <= x[1:length(weights)] <= 1, Int)
    @objective(mod, Max, dot(x, values))
    @constraint(mod, dot(x, weights) <=  W)

    function mycb_localzero(cb)
        nodesexpl = MathProgBase.cbgetexplorednodes(cb)
        if  entered[1] == false && nodesexpl >= 1
            # the following user cut  constrains all x[i] to be zero, but applies only locally at the first node after the root node, and doesn't preclude the existence non-trival "optimal" solutions
            @usercut(cb, sum{x[i], i=1:length(weights)} <= 0, localcut=true)
            # @usercut(cb, sum{x[i], i=1:length(weights)} <= 0) # applying the cut globally would lead the solver to x=0 as the optimal solution
            @fact macroexpand(:(@usercut(cb, sum{x[i], i=1:length(weights)} <= 0, badkwarg=true))).head --> :error
            entered[1] = true
        end
        entered[2] = true
    end
    addcutcallback(mod, mycb_localzero)
    addcutcallback(mod, cb -> (entered[3] = true))
    @fact solve(mod) --> :Optimal
    @fact entered --> [true,true,true]
    @fact sum(getvalue(x)) --> greater_than(0)
end; end; end

facts("[callback] Test heuristics") do
for heursolver in heur_solvers
context("With solver $(typeof(heursolver))") do
    entered = [false,false]

    N = 100
    # Include explicit data from srand(250) so that we can reproduce across platforms
    include(joinpath("data","heuristic.jl"))
    mod = Model(solver=heursolver)
    @variable(mod, x[1:N], Bin)
    @objective(mod, Max, dot(r1,x))
    @constraint(mod, dot(ones(N),x) <= rhs*N)
    function myheuristic1(cb)
        entered[1] == true && return
        entered[1] = true
        for i in 1:100
            if i in [9,10,11,14,15,16,25,30,32,41,44,49,50,53,54,98,100]
                setsolutionvalue(cb, x[i], 0)
            else
                setsolutionvalue(cb, x[i], 1)
            end
        end
        addsolution(cb)
    end
    addheuristiccallback(mod, myheuristic1)
    addheuristiccallback(mod, cb -> (entered[2] = true))
    @fact solve(mod) --> :Optimal
    @fact entered --> [true,true]
    @fact find(getvalue(x)[:]) --> setdiff(1:N,[9,10,11,14,15,16,25,30,32,41,44,49,50,53,54,98,100])

    empty!(mod.callbacks)
    entered[1] = false
    # Test that solver rejects infeasible partial solutions...
    # ...the second solution has higher objective value, but is infeasible
    function myheuristic2(cb)
        entered[1] == true && return
        entered[1] = true
        for i in 1:90 # not every component, but close
            setsolutionvalue(cb, x[i], 1)
        end
        addsolution(cb)
    end
    addheuristiccallback(mod, myheuristic2)
    addheuristiccallback(mod, cb -> (entered[2] = true))
    @fact solve(mod) --> :Optimal
    @fact entered --> [true,true]
    @fact find(getvalue(x)[:]) --> setdiff(1:N,[9,10,11,14,15,16,25,30,32,41,44,49,50,53,54,98,100])
end; end; end

facts("[callback] Test informational callback") do
for infosolver in info_solvers
context("With solver $(typeof(infosolver))") do
    nodes      = Int[]
    objs       = Float64[]
    objvals    = Float64[]
    bestbounds = Float64[]
    entered = fill(false, 2)

    N = 10000
    include(joinpath("data","informational.jl"))
    mod = Model(solver=infosolver)
    @variable(mod, x[1:N], Bin)
    @objective(mod, Max, dot(r1,x))
    @constraint(mod, c[i=1:10], dot(r2[i],x) <= rhs[i]*N/10)
    # Test that solver fills solution correctly
    function myinfo1(cb)
        @assert MathProgBase.cbgetstate(cb) == :Intermediate
        entered[1] = true
        push!(nodes,      MathProgBase.cbgetexplorednodes(cb))
        push!(objs,       MathProgBase.cbgetobj(cb))
        push!(bestbounds, MathProgBase.cbgetbestbound(cb))
    end
    function myinfo2(cb)
        @assert MathProgBase.cbgetstate(cb) == :MIPSol
        push!(objvals, dot(r1,JuMP.getvalue(x[:])))
    end
    addinfocallback(mod, myinfo1, when = :Intermediate)
    addinfocallback(mod, myinfo2, when = :MIPSol)
    addinfocallback(mod, cb -> (entered[2] = true), when = :Intermediate)
    @fact solve(mod) --> :Optimal
    @fact entered --> fill(true, 2)
    mono_node, mono_obj, mono_bestbound = true, true, true
    for n in 2:length(nodes)
        mono_node &= (nodes[n-1] <= nodes[n] + 1e-8)
        if nodes[n] > 0 # all bets are off at monotonicity at root node
            mono_obj &= (objs[n-1] <= objs[n] + 1e-8)
            mono_bestbound &= (bestbounds[n-1] >= bestbounds[n] - 1e-8)
        end
    end
    @fact mono_node      --> true
    @fact mono_obj       --> true
    @fact mono_bestbound --> true
    for (s1,s2) in zip(objvals[1:end-1], objvals[2:end])
        @fact  s1 <= s2 --> true
    end
end; end; end

facts("[callback] Callback exit on CallbackAbort") do
for solver in lazy_solvers
context("With solver $(typeof(solver))") do
    mod = Model(solver=solver)
    @variable(mod, 0 <= x <= 2, Int)
    @variable(mod, 0 <= y <= 2, Int)
    @objective(mod, Max, x + 2y)
    @constraint(mod, y + x <= 3.5)

    mycallback = _ -> throw(CallbackAbort())
    addlazycallback(mod, mycallback)
    @fact solve(mod) --> :UserLimit
end; end; end

cbc && facts("[callback] Solver doesn't support callbacks") do
    mycb(cb) = nothing
    mod = Model(solver=Cbc.CbcSolver())
    addlazycallback(mod, mycb)
    @fact_throws ErrorException solve(mod)
    mod = Model(solver=Cbc.CbcSolver())
    addcutcallback(mod, mycb)
    @fact_throws ErrorException solve(mod)
    mod = Model(solver=Cbc.CbcSolver())
    addheuristiccallback(mod, mycb)
    @fact_throws ErrorException solve(mod)
    mod = Model(solver=Cbc.CbcSolver())
    addinfocallback(mod, mycb, when = :Intermediate)
    @fact_throws ErrorException solve(mod)
end


text:
# MIT License

# Copyright (c) 2021 Fraunhofer-Chalmers Centre

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

########################################
# Runs the LR algorithm and Gurobi on all categories of instances in the "instance" folder and writes the result data to the "result" folder
# Results from our run is already in the results.7z folder
########################################
using DelimitedFiles
cd(@__DIR__ )
folders = [name  for  name in readdir("instances") if isnothing(findfirst('.',name))]
instance_names = Dict()
n_tot = 0
T_limit = 300
for folder in folders
    global  n_tot
    instance_name = []
    for file in readdir("instances/"*folder)
        if file[end-3:end] != ".csv"
            continue
        end
        inst_path = folder*"/"*file
        push!(instance_name,file[1:end-4])
        n_tot += 1
    end
    instance_names[folder]=instance_name
end
cp("../binary/RCmax.exe","RCmax.exe")
##
n_solved = 0
start = time()
for folder in folders
    global  n_solved
    out_path = "results/"*folder*"/res.txt"
    if isfile(out_path)
        rm(out_path)
    end

    for name in instance_names[folder]
        inst_path = "instances/"*folder*"/"*name*".csv"
        A = readdlm(inst_path,',')
        (m,n) = size(A)
        println(inst_path)
        if n <= 200
            run(`RCmax.exe $inst_path $T_limit -o $out_path`)
        else
            run(`RCmax.exe $inst_path $T_limit -o $out_path -nL -nS`)
        end
        n_solved += 1
        println(time()-start,"s ", n_solved,"/",n_tot)
    end
    name_v = vcat("name", instance_names[folder])
    A = readdlm(out_path, ',')
    writedlm(out_path, [name_v A],',')
end

## with Gurobi
using JuMP, Gurobi
function solveRCmax(inst_path)
    p = readdlm(inst_path,',')
    m,n = size(p)
    model = Model(Gurobi.Optimizer)
    @variable(model,z)
    @objective(model,Min,z)
    @variable(model, x[1:m,1:n], Bin)
    @constraint(model, [i in 1:m], p[i,:]'*x[i,:] <= z)
    @constraint(model, [j in 1:n], sum(x[:,j]) == 1)
    set_time_limit_sec(model,T_limit)
    set_silent(model)
    set_optimizer_attribute(model,"Threads", 1)
    set_optimizer_attribute(model, "Cuts", 2)
    set_optimizer_attribute(model,"MIPGap", 1e-8)
    set_optimizer_attribute(model,"MIPGapAbs", 0.999) # integer costs, ensure that it is not a round off error
    optimize!(model)

    zv = objective_value(model)
    zb = objective_bound(model)
    st = solve_time(model)
    si = simplex_iterations(model)
    nc = node_count(model)
    solved = st < T_limit

    unset_binary.(x)
    set_lower_bound.(x,0.0)
    optimize!(model)
    z_lp = objective_value(model)

    gap = (zv-zb)/zv
    lb_root_gap = (zv-z_lp)/zv

    return [m n nc si st solved gap lb_root_gap zv]
end

start = time()
n_solved = 0
for folder in folders
    global  n_solved
    out_path = "results/"*folder*"/res_gur_cut.txt"
    A = ["name" "m" "n" "nodes" "iter" "time" "solved" "gap" "lbrootgap" "z"]
    for name in instance_names[folder]
        inst_path = "instances/"*folder*"/"*name*".csv"
        println(inst_path)
        res = solveRCmax(inst_path)
        n_solved += 1
        println(time()-start,"s ", n_solved,"/",n_tot)
        A = vcat(A, hcat(name, res))
        writedlm(out_path, A, ",")
    end
end

max_stars_repo_name: Fraunhofer-Chalmers-Centre/RCmax

max_issues_repo_name: Fraunhofer-Chalmers-Centre/RCmax

max_forks_repo_name: Fraunhofer-Chalmers-Centre/RCmax


text:
# AUTO GENERATED FILE - DO NOT EDIT

export vtk_algorithm

"""
    vtk_algorithm(;kwargs...)
    vtk_algorithm(children::Any;kwargs...)
    vtk_algorithm(children_maker::Function;kwargs...)

An Algorithm component.
Algorithm is exposing a source or filter to a downstream filter
It takes the following set of properties:
  - vtkClass: vtkClassName
  - state: {}
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional)
- `id` (String; optional): The ID used to identify this component.
- `port` (Real; optional): downstream connection port
- `vtkClass` (String; optional): vtkClass name
- `state` (Dict; optional): set of property values for vtkClass
"""
function vtk_algorithm(; kwargs...)
        available_props = Symbol[:children, :id, :port, :vtkClass, :state]
        wild_props = Symbol[]
        return Component("vtk_algorithm", "Algorithm", "dash_vtk", available_props, wild_props; kwargs...)
end

vtk_algorithm(children::Any; kwargs...) = vtk_algorithm(;kwargs..., children = children)
vtk_algorithm(children_maker::Function; kwargs...) = vtk_algorithm(children_maker(); kwargs...)

max_stars_repo_name: ericbdaniels/dash-vtk

max_issues_repo_name: ericbdaniels/dash-vtk

max_forks_repo_name: ericbdaniels/dash-vtk


text:
numbers = [0,14,1,3,7,9]
# numbers = [1,3,2]
# numbers = [0,3,6]
# numbers = [2,1,3]

function get_next_number!(spoken, last_number, current_idx)
    if !(last_number in keys(spoken))
        spoken[last_number] = current_idx - 1
        return 0
    else
        result = current_idx - spoken[last_number] - 1
        spoken[last_number] = current_idx - 1
        return result
    end
end

spoken_history = Dict(n[2] => n[1] for n in enumerate(numbers[1:end-1]))
last_spoken = numbers[end]
for i in length(numbers)+1:30000000
    global spoken_history
    global last_spoken = get_next_number!(spoken_history, last_spoken, i)
end

println(last_spoken)


text:
struct LeastSquaresMonteCarlo{M<:MonteCarloModel}
    m::M
    degree::Int
end

function value(m::LeastSquaresMonteCarlo, c::AnytimeBefore{C}) where C
    uc = c.c # underlying contract
    uc1 = uc.c1
    uc2 = uc.c2
    N = date2index(m.m, maturitydate(c))
    R = discount(m.m.core.yieldcurve, m.m.dates[2]) # TODO: could be better...

    V = [valueat(ms, uc, N)*R for ms in scenarios(m.m)]
    for n = N-1:-1:2
        I = V .> 0 # in the money options
        A = [x^d for d = 0:m.degree, x in m.m.paths[:,n]] # design matrix
        β = A[:,I]' \ V[I] # least squares regression
        cV = A'*β          # estimated continuation values
        for (p,ms) in enumerate(scenarios(m.m))
            ev = valueat(ms,uc,n)  # exercise value
            if I[p] && cV[p] < ev
                # in-the-money and exercise is greater than continuation
                V[p] = ev*R
            else
                V[p] *= R
            end
        end
    end
    return max(mean(V), mean(valueat(ms, uc, N) for ms in scenarios(m.m)))
end

function value(m::GeomBMModel, c::Contract, ::Type{LeastSquaresMonteCarlo}, npaths::Integer, degree::Integer)
    mcm = montecarlo(m, startdate(m):Day(1):maturitydate(c), npaths)
    value(LeastSquaresMonteCarlo(mcm, degree), c)
end


text:
"""
    test_constvel(; kwargs...)

Advects a gaussian concentration c0(x, y, t) with a constant velocity flow
u(x, y) = uvel and v(x, y) = vvel and compares the final state with
cfinal = c0(x-uvel*tfinal, y-vvel*tfinal)
"""
function test_constvel(stepper, dt, nsteps, dev::Device=CPU())

  nx, Lx = 128, 2π
  uvel, vvel = 0.2, 0.1
  u(x, y) = uvel
  v(x, y) = vvel

  prob = TracerAdvectionDiffusion.Problem(dev; nx=nx, Lx=Lx, κ=0.0, u=u, v=v, dt=dt, stepper=stepper, steadyflow=true)
  sol, cl, vs, pr, gr = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

  x, y = gridpoints(gr)

  σ = 0.1
  c0ampl = 0.1
  c0func(x, y) = c0ampl * exp(-(x^2+y^2)/(2σ^2))

  c0 = c0func.(x, y)
  tfinal = nsteps*dt
  cfinal = @. c0func(x - uvel*tfinal, y - vvel*tfinal)

  TracerAdvectionDiffusion.set_c!(prob, c0)

  stepforward!(prob, nsteps)
  TracerAdvectionDiffusion.updatevars!(prob)

  isapprox(cfinal, vs.c, rtol=gr.nx*gr.ny*nsteps*1e-12)
end

"""
    test_timedependenttvel(; kwargs...)

Advects a gaussian concentration c0(x, y, t) with a time-varying velocity flow
u(x, y, t) = uvel and v(x, y, t) = vvel*sign(-t+tfinal/2) and compares the final
state with cfinal = c0(x-uvel*tfinal, y)
"""
function test_timedependentvel(stepper, dt, tfinal, dev::Device=CPU(); uvel=0.5, αv=0.5)

  nx, Lx = 128, 2π
  nsteps = round(Int, tfinal/dt)

  if !isapprox(tfinal, nsteps*dt, rtol=rtol_traceradvectiondiffusion)
    error("tfinal is not multiple of dt")
  end

  u(x, y, t) = uvel
  v(x, y, t) = αv*t + αv*dt/2

  prob = TracerAdvectionDiffusion.Problem(dev; nx=nx, Lx=Lx, κ=0.0, u=u, v=v, dt=dt, stepper=stepper)
  sol, cl, vs, pr, gr = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
  x, y = gridpoints(gr)

  σ = 0.2
  c0func(x, y) = 0.1*exp(-(x^2+y^2)/(2σ^2))

  c0 = @. c0func(x, y)
  tfinal = nsteps*dt
  cfinal = @. c0func(x - uvel*tfinal, y - 0.5αv*tfinal^2)

  TracerAdvectionDiffusion.set_c!(prob, c0)

  stepforward!(prob, nsteps)
  TracerAdvectionDiffusion.updatevars!(prob)
  isapprox(cfinal, vs.c, rtol=gr.nx*gr.ny*nsteps*1e-12)
end

"""
    test_diffusion(; kwargs...)

Diffuses a gaussian concentration c0(x, y, t) and compares the final state with
the analytic solution of the heat equation, cfinal
"""
function test_diffusion(stepper, dt, tfinal, dev::Device=CPU(); steadyflow = true)

  nx = 128
  Lx = 2π
   κ = 0.01
  nsteps = round(Int, tfinal/dt)

  if !isapprox(tfinal, nsteps*dt, rtol=rtol_traceradvectiondiffusion)
    error("tfinal is not multiple of dt")
  end

  prob = TracerAdvectionDiffusion.Problem(dev; steadyflow=steadyflow, nx=nx,
    Lx=Lx, κ=κ, dt=dt, stepper=stepper)
  sol, cl, vs, pr, gr = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
  x, y = gridpoints(gr)

  c0ampl, σ = 0.1, 0.1
  c0func(x, y) = c0ampl*exp(-(x^2+y^2)/(2σ^2))

  c0 = @. c0func(x, y)
  tfinal = nsteps*dt
  σt = sqrt(2*κ*tfinal + σ^2)
  cfinal = @. c0ampl*(σ^2/σt^2)*exp(-(x^2+y^2)/(2*σt^2))

  TracerAdvectionDiffusion.set_c!(prob, c0)

  stepforward!(prob, nsteps)
  TracerAdvectionDiffusion.updatevars!(prob)

  isapprox(cfinal, vs.c, rtol=gr.nx*gr.ny*nsteps*1e-12)
end

"""
    test_hyperdiffusion(; kwargs...)

Diffuses a gaussian concentration c0(x, y, t) using hyperdiffusivity and
compares the final state with the analytic solution of the heat equation, cfinal
"""
function test_hyperdiffusion(stepper, dt, tfinal, dev::Device=CPU(); steadyflow = true)

   nx = 128
   Lx = 2π
    κ = 0.0   # no diffusivity
    η = κ     # no diffusivity
   κh = 0.01  # hyperdiffusivity coeff
  nκh = 1     # nκh=1 converts hyperdiffusivity to plain diffusivity
              # so we can compare with the analytic solution of heat equation

  nsteps = round(Int, tfinal/dt)

  if !isapprox(tfinal, nsteps*dt, rtol=rtol_traceradvectiondiffusion)
    error("tfinal is not multiple of dt")
  end

  gr = TwoDGrid(dev, nx, Lx)
  x, y = gridpoints(gr)

  u, v = zero(x), zero(x) #0*x, 0*x

  vs = TracerAdvectionDiffusion.Vars(dev, gr)
  pr = TracerAdvectionDiffusion.ConstDiffSteadyFlowParams(η, κ, κh, nκh, u, v)
  eq = TracerAdvectionDiffusion.Equation(pr, gr)
  prob = FourierFlows.Problem(eq, stepper, dt, gr, vs, pr, dev)

  c0ampl, σ = 0.1, 0.1
  c0func(x, y) = c0ampl*exp(-(x^2+y^2)/(2σ^2))

  c0 = @. c0func(x, y)
  tfinal = nsteps*dt
  σt = sqrt(2*κh*tfinal + σ^2)
  cfinal = @. c0ampl*σ^2/σt^2 * exp(-(x^2+y^2)/(2*σt^2))

  TracerAdvectionDiffusion.set_c!(prob, c0)

  stepforward!(prob, nsteps)
  TracerAdvectionDiffusion.updatevars!(prob)

  isapprox(cfinal, vs.c, rtol=gr.nx*gr.ny*nsteps*1e-12)
end

max_stars_repo_name: FourierFlows/PassiveTracerFlows

max_issues_repo_name: FourierFlows/PassiveTracerFlows

max_forks_repo_name: FourierFlows/PassiveTracerFlows


text:
# lanczos.jl

mutable struct LanczosFactorization{T, S<:Real} <: KrylovFactorization{T,S}
    k::Int # current Krylov dimension
    V::OrthonormalBasis{T} # basis of length k
    αs::Vector{S}
    βs::Vector{S}
    r::T
end

Base.length(F::LanczosFactorization) = F.k
Base.sizehint!(F::LanczosFactorization, n) = begin
    sizehint!(F.V, n)
    sizehint!(F.αs, n)
    sizehint!(F.βs, n)
    return F
end
Base.eltype(F::LanczosFactorization) = eltype(typeof(F))
Base.eltype(::Type{<:LanczosFactorization{<:Any,S}}) where {S} = S

basis(F::LanczosFactorization) = length(F.V) == F.k ? F.V :
    error("Not keeping vectors during Lanczos factorization")
rayleighquotient(F::LanczosFactorization) = SymTridiagonal(F.αs, F.βs)
residual(F::LanczosFactorization) = F.r
@inbounds normres(F::LanczosFactorization) = F.βs[F.k]
rayleighextension(F::LanczosFactorization) = SimpleBasisVector(F.k, F.k)

# Lanczos iteration for constructing the orthonormal basis of a Krylov subspace.
struct LanczosIterator{F,T,O<:Orthogonalizer} <: KrylovIterator{F,T}
    operator::F
    x₀::T
    orth::O
    keepvecs::Bool
    function LanczosIterator{F,T,O}(operator::F, x₀::T, orth::O, keepvecs::Bool) where
        {F,T,O<:Orthogonalizer}

        if !keepvecs && isa(orth, Reorthogonalizer)
            error("Cannot use reorthogonalization without keeping all Krylov vectors")
        end
        new{F,T,O}(operator, x₀, orth, keepvecs)
    end
end
LanczosIterator(operator::F, x₀::T, orth::O = KrylovDefaults.orth,
                keepvecs::Bool = true) where {F,T,O<:Orthogonalizer} =
    LanczosIterator{F,T,O}(operator, x₀, orth, keepvecs)
LanczosIterator(A::AbstractMatrix, x₀::AbstractVector, orth::O = KrylovDefaults.orth,
                keepvecs::Bool = true) where {O<:Orthogonalizer} =
    LanczosIterator(x->A*x, x₀, orth, keepvecs)

Base.IteratorSize(::Type{<:LanczosIterator}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:LanczosIterator}) = Base.EltypeUnknown()

function Base.iterate(iter::LanczosIterator)
    state = initialize(iter)
    return state, state
end
function Base.iterate(iter::LanczosIterator, state::LanczosFactorization)
    nr = normres(state)
    if nr < eps(typeof(nr))
        return nothing
    else
        state = expand!(iter, deepcopy(state))
        return state, state
    end
end

function initialize(iter::LanczosIterator; verbosity::Int = 0)
    β₀ = norm(iter.x₀)
    β₀ == 0 && throw(ArgumentError("initial vector should not have norm zero"))
    invβ₀ = one(eltype(iter.x₀))/β₀
    T = typeof(invβ₀) # division might change eltype
    x₀ = mul!(similar(iter.x₀, T), iter.x₀, invβ₀)
    w = iter.operator(x₀) # applying the operator might change eltype
    v = eltype(x₀) == eltype(w) ? x₀ : copyto!(similar(w), x₀)
    r, α = orthogonalize!(w, v, iter.orth)
    β = norm(r)
    n = hypot(α,2*β)
    imag(α) <= sqrt(max(eps(n),eps(one(n)))) ||
        error("operator does not appear to be hermitian: $(imag(α)) vs $n")

    V = OrthonormalBasis([v])
    αs = [real(α)]
    βs = [β]
    if verbosity > 0
        @info "Lanczos iteration step 1: normres = $β"
    end

    return LanczosFactorization(1, V, αs, βs, r)
end
function initialize!(iter::LanczosIterator, state::LanczosFactorization; verbosity::Int = 0)
    x₀ = iter.x₀
    V = state.V
    while length(V) > 1
        pop!(V)
    end
    αs = empty!(state.αs)
    βs = empty!(state.βs)

    v = mul!(V[1], x₀, 1/norm(x₀))
    w = iter.operator(v)
    r, α = orthogonalize!(w, v, iter.orth)
    β = norm(r)
    n = hypot(α,β)
    imag(α) <= sqrt(max(eps(n),eps(one(n)))) ||
        error("operator does not appear to be hermitian: $(imag(α)) vs $n")

    state.k = 1
    push!(αs, real(α))
    push!(βs, β)
    state.r = r
    if verbosity > 0
        @info "Lanczos iteration step 1: normres = $β"
    end
    return state
end
function expand!(iter::LanczosIterator, state::LanczosFactorization; verbosity::Int = 0)
    βold = normres(state)
    V = state.V
    r = state.r
    V = push!(V, rmul!(r, 1/βold))
    r, α, β = lanczosrecurrence(iter.operator, V, βold, iter.orth)
    n = hypot(α, β, βold)
    imag(α) <= sqrt(max(eps(n),eps(one(n)))) ||
        error("operator does not appear to be hermitian: $(imag(α)) vs $n")

    αs = push!(state.αs, real(α))
    βs = push!(state.βs, β)

    !iter.keepvecs && popfirst!(state.V) # remove oldest V if not keepvecs

    state.k += 1
    state.r = r
    if verbosity > 0
        @info "Lanczos iteration step $(state.k): normres = $β"
    end
    return state
end
function shrink!(state::LanczosFactorization, k)
    length(state) == length(state.V) ||
        error("we cannot shrink LanczosFactorization without keeping Lanczos vectors")
    length(state) <= k && return state
    V = state.V
    while length(V) > k+1
        pop!(V)
    end
    r = pop!(V)
    resize!(state.αs, k)
    resize!(state.βs, k)
    state.k = k
    state.r = rmul!(r, normres(state))
    return state
end

# Exploit hermiticity to "simplify" orthonormalization process:
# Lanczos three-term recurrence relation
function lanczosrecurrence(operator, V::OrthonormalBasis, β, orth::ClassicalGramSchmidt)
    v = V[end]
    w = apply(operator, v)
    α = dot(v, w)
    w = axpy!(-β, V[end-1], w)
    w = axpy!(-α, v, w)
    β = norm(w)
    return w, α, β
end
function lanczosrecurrence(operator, V::OrthonormalBasis, β, orth::ModifiedGramSchmidt)
    v = V[end]
    w = apply(operator, v)
    w = axpy!( -β, V[end-1], w)
    α = dot(v, w)
    w = axpy!(-α, v, w)
    β = norm(w)
    return w, α, β
end
function lanczosrecurrence(operator, V::OrthonormalBasis, β, orth::ClassicalGramSchmidt2)
    v = V[end]
    w = operator(v)
    α = dot(v, w)
    w = axpy!(-β, V[end-1], w)
    w = axpy!(-α, v, w)

    w, s = orthogonalize!(w, V, ClassicalGramSchmidt())
    α += s[end]
    β = norm(w)
    return w, α, β
end
function lanczosrecurrence(operator, V::OrthonormalBasis, β, orth::ModifiedGramSchmidt2)
    v = V[end]
    w = operator(v)
    w = axpy!( -β, V[end-1], w)
    w, α = orthogonalize!(w, v, ModifiedGramSchmidt())

    s = α
    for q in V
        w, s = orthogonalize!(w, q, ModifiedGramSchmidt())
    end
    α += s
    β = norm(w)
    return w, α, β
end
function lanczosrecurrence(operator, V::OrthonormalBasis, β, orth::ClassicalGramSchmidtIR)
    v = V[end]
    w = apply(operator, v)
    α = dot(v, w)
    w = axpy!(-β, V[end-1], w)
    w = axpy!(-α, v, w)

    ab2 = abs2(α) + abs2(β)
    β = norm(w)
    nold = sqrt(abs2(β)+ab2)
    while eps(one(β)) < β < orth.η * nold
        nold = β
        w, s = orthogonalize!(w, V, ClassicalGramSchmidt())
        α += s[end]
        β = norm(w)
    end
    return w, α, β
end
function lanczosrecurrence(operator, V::OrthonormalBasis, β, orth::ModifiedGramSchmidtIR)
    v = V[end]
    w = apply(operator, v)
    w = axpy!( -β, V[end-1], w)

    w, α = orthogonalize!(w, v, ModifiedGramSchmidt())
    ab2 = abs2(α) + abs2(β)
    β = norm(w)
    nold = sqrt(abs2(β)+ab2)
    while eps(one(β)) < β < orth.η * nold
        nold = β
        s = zero(α)
        for q in V
            w, s = orthogonalize!(w, q, ModifiedGramSchmidt())
        end
        α += s
        β = norm(w)
    end
    return w, α, β
end


text:
"""
    degrees(G::SimpleGraph)

Compute the sequence of degrees
"""
function degrees(G::SimpleGraph)
    degrees(CombinatorialAdjacency(adjacency_matrix(G)))
end

"""
    neighbors_edges(G::SimpleGraph)

Returns deg the sequence of degrees and neig a list of list of direcetd edges numbered from 1 to 2*ne(G)
neig[i] is the list of edges j->i
"""
function neighbors_edges(G::SimpleGraph)
    deg = degrees(G)
    neig = [Vector{Int64}() for _ in vertices(G)]
    c = 1
    for (i,d) in enumerate(deg)
        neig[i] = c:(c+d-1)
        c +=d
    end
    deg, neig
end

"""
    Network

`deg` is a `2*ne(G)` vectors with the degree of the j for each edge `j->i`
`edges_arr` is a `2*ne(G)x2` array containing all directed edges `j->i` where `i` is increasing.
`neig` is a list of vectors where `neig[e]` are all indices of edges `u->j` for `u!=i` where edge `e`
 is `j->i`
"""
struct Network
    G::SimpleGraph
    deg::Vector{Int64}
    edges_arr::Array{Int64,2}
    neig::Vector{Vector{Int64}}
end

"""
    create_edges(G::SimpleGraph)

returns a 2*ne(G)x2 array containing all directed edges j->i where i is increasing.
"""
function create_edges(G::SimpleGraph)
    edges_array = Array{Int64,2}(undef,2*ne(G),2)
    c = 1
    for i in vertices(G)
        for j in G.fadjlist[i]
            edges_array[c,:] = [j i]
            c +=1
        end
    end
    edges_array
end

"""
    neig_remove(G::SimpleGraph,i::Int64,j::Int64)

returns the neighbors of i in G whitout j
"""
function neig_remove(G::SimpleGraph,i::Int64,j::Int64)
    neig = neighbors(G,i)
    if j ∉ neig
        error(j, "is not a neighbor of ", i, "in ", neig)
    end
    filter(s -> s!=j, neig)
end

"""
    create_dic(edges_arr::Array{Int64,2})

returns a dic such that dic[(i,j)] = index of edge i->j in edges_arr
"""
function create_dic(edges_arr::Array{Int64,2})
    ne, _ = size(edges_arr)
    local_dic = Dict{Tuple,Int64}()
    for i in 1:ne
        local_dic[(edges_arr[i,1],edges_arr[i,2])] = i
    end
    local_dic
end

"""
    create_neig_r(G::SimpleGraph,edges_arr::Array{Int64,2})

returns a list of vectors where neig_r[i] are all indices of edges u->v for u!=w where edge i is v->w
"""
function create_neig_r(G::SimpleGraph,edges_arr::Array{Int64,2})
    ne, _ = size(edges_arr)
    neig_r = [Vector{Int64}() for _ in 1:ne]
    local_dic = create_dic(edges_arr)
    for i in 1:ne
        neig_r[i] = map(s -> local_dic[(s,edges_arr[i,1])], neig_remove(G,edges_arr[i,1],edges_arr[i,2]))
    end
    neig_r
end

"""
    create_degrees(G::SimpleGraph,edges_arr::Array{Int64,2})

retunrs a list where degrees_G[i] is the degree of u where edge i is u->v
"""
function create_degrees(G::SimpleGraph,edges_arr::Array{Int64,2})
    ne, _ = size(edges_arr)
    degrees_G = zeros(ne)
    for i in 1:ne
        degrees_G[i] = degree(G,edges_arr[i,1])
    end
    degrees_G
end

function Network(G::SimpleGraph)
    edges_arr = create_edges(G)
    neig = create_neig_r(G,edges_arr)
    deg = create_degrees(G,edges_arr)
    return Network(G, deg, edges_arr, neig)
end

"""
    Pair_ER

for `i` in `G1`, `perm_G1G2[i]` is the corresponding vertex in `G2`. It is `0` if no corresponding vertex in `G2`.
`edge_G1G2` used for init of bias message passing (TBC).
"""
struct Pair_ER
    N1::Network
    N2::Network
    Ginter::SimpleGraph
    n::Int64
    λ::Float64
    s::Float64
    perm_G1G2::Array{Int64,1}
    perm_G2G1::Array{Int64,1}
    edge_G1G2::Array{Int64,2}
    dic_comb::Dict{Tuple{Int64, Int64}, Vector{Vector{Int64}}}
    dic_perm::Dict{Int64, Vector{Vector{Int64}}}
end

"""
    compute_ab(n,λ,s)

returns the solution of the system of equations:
a+(1-a)*b^2 = λ*s/n and (1-a)*b*(1-b) = λ*(1-s)/n
"""
function compute_ab(n,λ,s)
    λ*(s-λ/n)/(λ*(s-2)+n), λ*(1-s)/(n-λ)
end

"""
    create_cor_graphs(n,λ,s)

returns correlated ER graphs G1, G2 and their intersection. All graphs are aligned.
"""
function create_cor_graphs(n,λ,s)
    a,b = compute_ab(n,λ,s)
    G_inter = erdos_renyi(n,a)
    G_add1 = erdos_renyi(n,b)
    G_add2 = erdos_renyi(n,b)
    G_1 = union(G_inter,G_add1)
    G_2 = union(G_inter,G_add2)
    return G_1, G_2, intersect(G_1,G_2)
end

function create_perm(n, vmap1, vmap2)
    revert_vmap2 = zeros(Int64,n)
    for (i,v) in enumerate(vmap2)
        revert_vmap2[v] = i
    end
    [revert_vmap2[i] for i in vmap1]
end

function Pair_ER(n,λ,s)
    G1,G2,Ginter = create_cor_graphs(n,λ,s)
    GG1, vmap1 = largest_comp(G1)
    GG2, vmap2 = largest_comp(G2)
    GGinter, vmapinter = largest_comp(Ginter)
    perm_GinterG1 = create_perm(n, vmapinter, vmap1)
    perm_GinterG2 = create_perm(n, vmapinter, vmap2)
    perm_G1G2 = create_perm(n, vmap1, vmap2)
    perm_G2G1 = create_perm(n, vmap2, vmap1)
    N1 = Network(GG1)
    N2 = Network(GG2)
    dic_G1 = create_dic(N1.edges_arr)
    dic_G2 = create_dic(N2.edges_arr)
    Ninter= Network(GGinter)
    ne , _ = size(Ninter.edges_arr)
    edge_G1G2 = Array{Int64,2}(undef,ne,2)
    for (i,e) in enumerate(eachrow(Ninter.edges_arr))
        #println(i,e)
        edge_G1G2[i,:] = [dic_G1[(perm_GinterG1[e[1]],perm_GinterG1[e[2]])], dic_G2[(perm_GinterG2[e[1]],perm_GinterG2[e[2]])]]
    end
    deg1 = unique(sort(N1.deg))
    deg2 = unique(sort(N2.deg))
    min_deg = min(deg1[end],deg2[end])
    max_deg = max(deg1[end],deg2[end])

    dic_comb = Dict([((d,k),collect(combinations(range(1,d),k))) for d in 0:max_deg for k in 0:d])
    dic_perm = Dict([(k,collect(permutations(range(1,k)))) for k in 0:min_deg])
    return Pair_ER(N1,N2,GGinter,n,λ,s,perm_G1G2,perm_G2G1,edge_G1G2,dic_comb,dic_perm)
end

"""
    largest_comp(G::SimpleGraph)

returns the largest connected component of `G` and
for i a vertex of this component `vmap[i]` is the corresponding index in the original `G`.
"""
function largest_comp(G::SimpleGraph)
    cc = connected_components(G)
    i = argmax([length(c) for c in cc])
    sg, vmap = induced_subgraph(G,cc[i])
    return sg, vmap
end


text:
# -----------------------------------------------
# Total Welfare
# -----------------------------------------------

@defcomp welfare begin

    regions  = Index()                         # 12 regions in RICE2010.

    ρ       = Parameter()                      # Pure rate of time preference.
    η       = Parameter()                      # Elasticity of marginal utility of consumption.
    pop     = Parameter(index=[time, regions]) # Regional population levels (millions of people).
    cpc     = Parameter(index=[time, regions]) # Per capita consumption levels ($1000s).

    UTILITY = Variable()                       # Total economic welfare over model time horizon (using the same name as the default RICE welfare component).

    function run_timestep(p, v, d, t)

        # Calculate total welfare for each period. Note, if η = 1, the social welfare function becomes log(x).
        if is_first(t)
            if p.η == 1.0
                v.UTILITY = sum(log.(p.cpc[t,:]) .* p.pop[t,:]) / (1.0 + p.ρ)^(10*(t.t - 1))
            else
                v.UTILITY = sum((p.cpc[t,:] .^ (1.0 - p.η)) ./ (1.0 - p.η) .* p.pop[t,:]) / (1.0 + p.ρ)^(10*(t.t - 1))
            end
        else
            if p.η == 1.0
                v.UTILITY = v.UTILITY + sum(log.(p.cpc[t,:]) .* p.pop[t,:]) / (1.0 + p.ρ)^(10*(t.t - 1))
            else
                v.UTILITY = v.UTILITY + sum((p.cpc[t,:] .^ (1.0 - p.η)) ./ (1.0 - p.η) .* p.pop[t,:]) / (1.0 + p.ρ)^(10*(t.t - 1))
            end
        end
    end
end

max_stars_repo_name: Environment-Research/Utilitarianism

max_issues_repo_name: Environment-Research/Utilitarianism

max_forks_repo_name: Environment-Research/Utilitarianism