text:
"""
$SIGNATURES

Initialize an empty in-memory DistributedFactorGraph `::DistributedFactorGraph` object.
"""
function initfg(dfg::T=LocalDFG(solverParams=SolverParams());
                                sessionname="NA",
                                robotname="",
                                username="",
                                cloudgraph=nothing) where T <: AbstractDFG
#
  #
  return dfg
end

#init an empty fg with a provided type and SolverParams
function initfg(::Type{T};solverParams=SolverParams(),
                      sessionname="NA",
                      robotname="",
                      username="",
                      cloudgraph=nothing) where T <: AbstractDFG
  #
  return T(solverParams=solverParams)
end

function initfg(::Type{T},solverParams::S;
                      sessionname="NA",
                      robotname="",
                      username="",
                      cloudgraph=nothing) where {T <: AbstractDFG, S <: SolverParams}
  #
  return T{S}(solverParams=solverParams)
end

# Should deprecate in favor of TensorCast.jl
reshapeVec2Mat(vec::Vector, rows::Int) = reshape(vec, rows, round(Int,length(vec)/rows))

## ==============================================================================================
## MOVE TO / CONSOLIDATE WITH DFG
## ==============================================================================================

"""
    $(SIGNATURES)

Fetch the variable marginal joint sampled points.  Use [`getBelief`](@ref) to retrieve the full Belief object.
"""
getVal(v::DFGVariable; solveKey::Symbol=:default) = v.solverDataDict[solveKey].val
getVal(v::DFGVariable, idx::Int; solveKey::Symbol=:default) = v.solverDataDict[solveKey].val[:,idx]
getVal(vnd::VariableNodeData) = vnd.val
getVal(vnd::VariableNodeData, idx::Int) = vnd.val[:, idx]
getVal(dfg::AbstractDFG, lbl::Symbol; solveKey::Symbol=:default) = getVariable(dfg, lbl).solverDataDict[solveKey].val

"""
    $(SIGNATURES)

Get the number of points used for the current marginal belief estimate represtation for a particular variable in the factor graph.
"""
function getNumPts(v::DFGVariable; solveKey::Symbol=:default)::Int
  return length(getVal(getSolverData(v, solveKey)))
end

function AMP.getBW(vnd::VariableNodeData)
  return vnd.bw
end

# setVal! assumes you will update values to database separate, this used for local graph mods only
function getBWVal(v::DFGVariable; solveKey::Symbol=:default)
  return getSolverData(v, solveKey).bw
end
function setBW!(vd::VariableNodeData, bw::Array{Float64,2}; solveKey::Symbol=:default)
  vd.bw = bw
  nothing
end
function setBW!(v::DFGVariable, bw::Array{Float64,2}; solveKey::Symbol=:default)
  setBW!(getSolverData(v, solveKey), bw)
  nothing
end

function setVal!(vd::VariableNodeData, val::AbstractVector{P}) where P
  vd.val = val
  nothing
end
function setVal!(v::DFGVariable, val::AbstractVector{P}; solveKey::Symbol=:default) where P
    setVal!(getSolverData(v, solveKey), val)
    nothing
end
function setVal!(vd::VariableNodeData, val::AbstractVector{P}, bw::Array{Float64,2}) where P
    setVal!(vd, val)
    setBW!(vd, bw)
    nothing
end
function setVal!(v::DFGVariable, val::AbstractVector{P}, bw::Array{Float64,2}; solveKey::Symbol=:default) where P
  setVal!(v, val, solveKey=solveKey)
  setBW!(v, bw, solveKey=solveKey)
  nothing
end
function setVal!(vd::VariableNodeData, val::AbstractVector{P}, bw::Vector{Float64}) where P
  setVal!(vd, val, reshape(bw,length(bw),1))
  nothing
end
function setVal!(v::DFGVariable, val::AbstractVector{P}, bw::Vector{Float64}; solveKey::Symbol=:default) where P
  setVal!(getSolverData(v, solveKey), val, bw)
  nothing
end
function setVal!(dfg::AbstractDFG, sym::Symbol, val::AbstractVector{P}; solveKey::Symbol=:default) where P
  setVal!(getVariable(dfg, sym), val, solveKey=solveKey)
end

"""
    $SIGNATURES

Set the point centers and bandwidth parameters of a variable node, also set `isInitialized=true` if `setinit::Bool=true` (as per default).

Notes
- `initialized` is used for initial solve of factor graph where variables are not yet initialized.
- `inferdim` is used to identify if the initialized was only partial.
"""
function setValKDE!(vd::VariableNodeData,
                    pts::AbstractVector{P},
                    bws::Vector{Float64},
                    setinit::Bool=true,
                    ipc::AbstractVector{<:Real}=[0.0;] ) where P
  #

  setVal!(vd, pts, bws) # BUG ...al!(., val, . ) ## TODO -- this can be a little faster
  setinit ? (vd.initialized = true) : nothing
  vd.infoPerCoord = ipc
  nothing
end

function setValKDE!(vd::VariableNodeData,
                    val::AbstractVector{P},
                    setinit::Bool=true,
                    ipc::AbstractVector{<:Real}=[0.0;]  ) where P
  # recover variableType information
  varType = getVariableType(vd)
  p = AMP.manikde!(varType, val)
  setValKDE!(vd, p, setinit, ipc)
  nothing
end

function setValKDE!(v::DFGVariable,
                    val::AbstractVector{P},
                    bws::Array{<:Real,2},
                    setinit::Bool=true,
                    ipc::AbstractVector{<:Real}=[0.0;];
                    solveKey::Symbol=:default) where P
  # recover variableType information
  setValKDE!(getSolverData(v, solveKey), val, bws[:,1], setinit, ipc )

  nothing
end

function setValKDE!(v::DFGVariable,
                    val::AbstractVector{P},
                    setinit::Bool=true,
                    ipc::AbstractVector{<:Real}=[0.0;];
                    solveKey::Symbol=:default) where P
  # recover variableType information
  setValKDE!(getSolverData(v, solveKey),val, setinit, ipc )
  nothing
end
function setValKDE!(v::DFGVariable,
                    em::TreeBelief,
                    setinit::Bool=true;
                    # inferdim::Union{Float32, Float64, Int32, Int64}=0;
                    solveKey::Symbol=:default  )
  #
  setValKDE!(v, em.val, em.bw, setinit, em.infoPerCoord, solveKey=solveKey)
  nothing
end
function setValKDE!(v::DFGVariable,
                    mkd::ManifoldKernelDensity,
                    setinit::Bool=true,
                    ipc::AbstractVector{<:Real}=[0.0;];
                    solveKey::Symbol=:default  )
  #
  # @error("TESTING setValKDE! ", solveKey, string(listSolveKeys(v)))
  setValKDE!(getSolverData(v,solveKey),mkd,setinit, Float64.(ipc))
  nothing
end
function setValKDE!(dfg::AbstractDFG,
                    sym::Symbol,
                    mkd::ManifoldKernelDensity,
                    setinit::Bool=true,
                    ipc::AbstractVector{<:Real}=[0.0;];
                    solveKey::Symbol=:default  )
  #
  setValKDE!(getVariable(dfg, sym), mkd, setinit, ipc, solveKey=solveKey)
  nothing
end

function setValKDE!(vnd::VariableNodeData,
                    mkd::ManifoldKernelDensity{M,B,Nothing},
                    setinit::Bool=true,
                    ipc::AbstractVector{<:Real}=[0.0;] ) where {M,B}
  #
  # L==Nothing means no partials
  ptsArr = AMP.getPoints(mkd) # , false) # for not partial
  # also set the bandwidth
  bws = getBW(mkd)[:,1]
  setValKDE!(vnd,ptsArr,bws,setinit,ipc )
  nothing
end

function setValKDE!(vnd::VariableNodeData,
                    mkd::ManifoldKernelDensity{M,B,L},
                    setinit::Bool=true,
                    ipc::AbstractVector{<:Real}=[0.0;] ) where {M,B,L<:AbstractVector}
  #
  oldBel = getBelief(vnd)

  # New infomation might be partial
  newBel = replace(oldBel, mkd)

  # Set partial dims as Manifold points
  ptsArr = AMP.getPoints(newBel, false)

  # also get the bandwidth
  bws = getBandwidth(newBel, false)

  # update values in graph
  setValKDE!(vnd,ptsArr,bws,setinit,ipc )
  nothing
end

"""
    $SIGNATURES

Set variable initialized status.
"""
function setVariableInitialized!( varid::VariableNodeData,
                                  status::Bool)
  #
  varid.initialized = status
end

setVariableInitialized!(vari::DFGVariable, status::Bool) = setVariableInitialized!(getSolverData(vari), status)

"""
    $SIGNATURES

Set method for the inferred dimension value in a variable.
"""
setIPC!(varid::VariableNodeData, val::AbstractVector{<:Real}) = varid.infoPerCoord = val
setIPC!(vari::DFGVariable, val::AbstractVector{<:Real}, solveKey::Symbol=:default) = setVariableIPC!(getSolverData(vari, solveKey), val)

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

"""
    $(SIGNATURES)

Get a ManifoldKernelDensity estimate from variable node data.
"""
getBelief(vnd::VariableNodeData) = manikde!(getManifold(getVariableType(vnd)), getVal(vnd), bw=getBW(vnd)[:,1] )

getBelief(v::DFGVariable, solvekey::Symbol=:default) = getBelief(getSolverData(v, solvekey))
getBelief(dfg::AbstractDFG, lbl::Symbol, solvekey::Symbol=:default) = getBelief(getVariable(dfg, lbl), solvekey)

"""
    $SIGNATURES

Reset the solve state of a variable to uninitialized/unsolved state.
"""
function resetVariable!(varid::VariableNodeData;
                        solveKey::Symbol=:default  )::Nothing
  #
  val = getBelief(varid)
  pts = AMP.getPoints(val)
  # TODO not all manifolds will initialize to zero
  for pt in pts
    fill!(pt, 0.0)
  end
  pn = manikde!(getManifold(varid), pts, bw=zeros(Ndim(val)))
  setValKDE!(varid, pn, false, [0.0;])
  # setVariableInferDim!(varid, 0)
  # setVariableInitialized!(vari, false)
  nothing
end

resetVariable!(vari::DFGVariable; solveKey::Symbol=:default  ) = resetVariable!(getSolverData(vari), solveKey=solveKey)

function resetVariable!(dfg::G,
                        sym::Symbol;
                        solveKey::Symbol=:default  )::Nothing where G <: AbstractDFG
  #
  resetVariable!(getVariable(dfg, sym), solveKey=solveKey)
end

function DefaultNodeDataParametric( dodims::Int,
                                    dims::Int,
                                    variableType::InferenceVariable;
                                    initialized::Bool=true,
                                    dontmargin::Bool=false)::VariableNodeData

  # this should be the only function allocating memory for the node points
  if false && initialized
    error("not implemented yet")
    # pN = AMP.manikde!(variableType.manifold, randn(dims, N));
    #
    # sp = Int[0;] #round.(Int,range(dodims,stop=dodims+dims-1,length=dims))
    # gbw = getBW(pN)[:,1]
    # gbw2 = Array{Float64}(undef, length(gbw),1)
    # gbw2[:,1] = gbw[:]
    # pNpts = getPoints(pN)
    # #initval, stdev
    # return VariableNodeData(pNpts,
    #                         gbw2, Symbol[], sp,
    #                         dims, false, :_null, Symbol[], variableType, true, 0.0, false, dontmargin)
  else
    sp = round.(Int,range(dodims,stop=dodims+dims-1,length=dims))
    ϵ = getPointIdentity(variableType)
    return VariableNodeData([ϵ],
                            zeros(dims,dims), Symbol[], sp,
                            dims, false, :_null, Symbol[], variableType, false, zeros(dims), false, dontmargin, 0, 0, :parametric)
  end

end

function setDefaultNodeDataParametric!(v::DFGVariable, variableType::InferenceVariable; kwargs...)
  vnd = DefaultNodeDataParametric(0, variableType |> getDimension, variableType; kwargs...)
  setSolverData!(v, vnd, :parametric)
  return nothing
end

"""
    $SIGNATURES

Create new solverData.

Notes
- Used during creation of new variable, as well as in CSM unique `solveKey`.
"""
function setDefaultNodeData!( v::DFGVariable,
                              dodims::Int,
                              N::Int,
                              dims::Int;
                              solveKey::Symbol=:default,
                              gt=Dict(),
                              initialized::Bool=true,
                              dontmargin::Bool=false,
                              varType=nothing )
  #
  # TODO review and refactor this function, exists as legacy from pre-v0.3.0
  # this should be the only function allocating memory for the node points (unless number of points are changed)
  data = nothing
  isinit = false
  sp = Int[0;]
  (valpts, bws) = if initialized
    pN = resample(getBelief(v))
    bws = getBW(pN)[:,1:1]
    pNpts = getPoints(pN)
    isinit = true
    (pNpts, bws)
  else
    sp = round.(Int,range(dodims,stop=dodims+dims-1,length=dims))
    @assert getPointType(varType) != DataType "cannot add manifold point type $(getPointType(varType)), make sure the identity element argument in @defVariable $varType arguments is correct"
    valpts = Vector{getPointType(varType)}(undef,N)
    for i in 1:length(valpts)
      valpts[i] = getPointIdentity(varType)
    end
    bws = zeros(dims,1)
    #
    (valpts, bws)
  end
  # make and set the new solverData
  setSolverData!(v, VariableNodeData( valpts, bws,
                                      Symbol[], sp,
                                      dims, false, :_null, Symbol[],
                                      varType, isinit, zeros(getDimension(v)), false, dontmargin,0,0,solveKey), solveKey)
  return nothing
end
# if size(initval,2) < N && size(initval, 1) == dims
#   @warn "setDefaultNodeData! -- deprecated use of stdev."
#   p = manikde!(varType.manifold, initval,diag(stdev));
#   pN = resample(p,N)
# if size(initval,2) < N && size(initval, 1) != dims
  # @info "Node value memory allocated but not initialized"
# else
#   pN = manikde!(varType.manifold, initval)
# end
# dims = size(initval,1) # rows indicate dimensions

"""
    $SIGNATURES

Reference data can be stored in the factor graph as a super-solve.

Notes
- Intended as a mechanism to store reference data alongside the numerical computations.
"""
function setVariableRefence!( dfg::AbstractDFG,
                              sym::Symbol,
                              val::AbstractVector;
                              refKey::Symbol=:reference)
  #
  # which variable to update
  var = getVariable(dfg, sym)

  # Construct an empty VND object
  vnd = VariableNodeData( val,
                          zeros(getDimension(var),1),
                          Symbol[],
                          Int[0;],
                          getDimension(var),
                          false,
                          :_null,
                          Symbol[],
                          getVariableType(var),
                          true,
                          zeros(getDimension(var)),
                          false,
                          true  )
  #
  # set the value in the DFGVariable
  setSolverData!(var, vnd, refKey)
end

# get instance from variableType
_variableType(varType::InferenceVariable) = varType
_variableType(varType::Type{<:InferenceVariable}) = varType()

## ==================================================================================================
## DFG Overloads on addVariable! and addFactor!
## ==================================================================================================

"""
$(SIGNATURES)

Add a variable node `label::Symbol` to `dfg::AbstractDFG`, as `varType<:InferenceVariable`.

Notes
-----
- keyword `nanosecondtime` is experimental and intended as the whole subsection portion -- i.e. accurateTime = (timestamp MOD second) + Nanosecond

Example
-------

```julia
fg = initfg()
addVariable!(fg, :x0, Pose2)
```
"""
function addVariable!(dfg::AbstractDFG,
                      label::Symbol,
                      varTypeU::Union{T, Type{T}};
                      N::Int=getSolverParams(dfg).N,
                      solvable::Int=1,
                      timestamp::Union{DateTime,ZonedDateTime}=now(localzone()),
                      nanosecondtime::Union{Nanosecond,Int64,Nothing}=nothing,
                      dontmargin::Bool=false,
                      labels::Union{Vector{Symbol},Nothing}=nothing,
                      tags::Vector{Symbol}=Symbol[],
                      smalldata=Dict{Symbol, DFG.SmallDataTypes}(),
                      checkduplicates::Bool=true,
                      initsolvekeys::Vector{Symbol}=getSolverParams(dfg).algorithms ) where T<:InferenceVariable
  #
  varType = _variableType(varTypeU)
  # TODO Remove deprecation in v0.16
  if :ut in fieldnames(T)
    Base.depwarn("Field `ut` (microseconds) for variable type ($T) has been deprecated please use DFGVariable.nstime, kwarg: nanosecondtime", :addVariable!)
    if isnothing(nanosecondtime)
      varType.ut == -(2^(Sys.WORD_SIZE-1)-1) && error("please define a time for type $(T), use FGVariable.nstime, kwarg: nanosecondtime")
      nanosecondtime = Nanosecond(varType.ut*1000)
    else
      @warn "Nanosecond time has been specified as $nanosecondtime, ignoring `ut` field value: $(varType.ut)."
    end
  elseif isnothing(nanosecondtime)
    nanosecondtime = Nanosecond(0)
  end

  # deprecate in v0.16
  labels isa Vector ? (union!(tags, labels); @warn("labels is deprecated, use tags instead")) : nothing

  union!(tags, [:VARIABLE])
  v = DFGVariable(label, varType; tags=Set(tags), smallData=smalldata, solvable=solvable, timestamp=timestamp, nstime=Nanosecond(nanosecondtime))

  (:default in initsolvekeys) &&
    setDefaultNodeData!(v, 0, N, getDimension(varType), initialized=false, varType=varType, dontmargin=dontmargin) # dodims

  (:parametric in initsolvekeys) &&
    setDefaultNodeDataParametric!(v, varType, initialized=false, dontmargin=dontmargin)

  DFG.addVariable!(dfg, v)

  return v
end

function parseusermultihypo(multihypo::Nothing, nullhypo::Float64)
  verts = Symbol[]
  mh = nothing
  return mh, nullhypo
end
function parseusermultihypo(multihypo::Vector{Float64}, nullhypo::Float64)
  mh = nothing
  if 0 < length(multihypo)
    multihypo2 = multihypo
    multihypo2[1-1e-10 .< multihypo] .= 0.0
    # check that terms sum to full probability
    @assert abs(sum(multihypo2) % 1) < 1e-10  || 1-1e-10 < sum(multihypo2) % 1 "ensure multihypo sums to a (or nearly, 1e-10) interger, see #1086"
    # check that only one variable broken into fractions
    @assert sum(multihypo2[1e-10 .< multihypo2]) ≈ 1
    # force normalize something that is now known to be close
    multihypo2 ./= sum(multihypo2)

    mh = Categorical(Float64[multihypo2...] )
  end
  return mh, nullhypo
end

# return a BitVector masking the fractional portion, assuming converted 0's on 100% confident variables
_getFractionalVars(varList::Union{<:Tuple, <:AbstractVector}, mh::Nothing) = zeros(length(varList)) .== 1
_getFractionalVars(varList::Union{<:Tuple, <:AbstractVector}, mh::Categorical) = 0 .< mh.p

function _selectHypoVariables(allVars::Union{<:Tuple, <:AbstractVector},
                              mh::Categorical,
                              sel::Integer = rand(mh) )
  #
  mask = mh.p .≈ 0.0
  mask[sel] = true
  (1:length(allVars))[mask]
end

_selectHypoVariables(allVars::Union{<:Tuple, <:AbstractVector},mh::Nothing,sel::Integer=0 ) = collect(1:length(allVars))

"""
    $SIGNATURES

Overload for specific factor preamble usage.

Notes:
- See https://github.com/JuliaRobotics/IncrementalInference.jl/issues/1462

DevNotes
- Integrate into CalcFactor
  - Add threading

Example:

```julia
import IncrementalInference: preableCache

preableCache(dfg::AbstractDFG, vars::AbstractVector{<:DFGVariable}, usrfnc::MyFactor) = MyFactorCache(randn(10))

# continue regular use, e.g.
mfc = MyFactor(...)
addFactor!(fg, [:a;:b], mfc)
# ...
```
"""
preambleCache(dfg::AbstractDFG, vars::AbstractVector{<:DFGVariable}, usrfnc::AbstractFactor) = nothing

# TODO perhaps consolidate with constructor?
"""
$SIGNATURES

Generate the default factor data for a new DFGFactor.
"""
function getDefaultFactorData(dfg::AbstractDFG,
                              Xi::Vector{<:DFGVariable},
                              usrfnc::T;
                              multihypo::Vector{<:Real}=Float64[],
                              nullhypo::Float64=0.0,
                              threadmodel=SingleThreaded,
                              eliminated::Bool = false,
                              potentialused::Bool = false,
                              edgeIDs = Int[],
                              solveInProgress = 0,
                              inflation::Real=getSolverParams(dfg).inflation,
                              _blockRecursion::Bool=false ) where T <: AbstractFactor
  #

  # prepare multihypo particulars
  # storeMH::Vector{Float64} = multihypo == nothing ? Float64[] : [multihypo...]
  mhcat, nh = parseusermultihypo(multihypo, nullhypo)

  # allocate temporary state for convolutional operations (not stored)
  userCache = preambleCache(dfg, Xi, usrfnc)
  ccw = _prepCCW(Xi, usrfnc; multihypo=mhcat, nullhypo=nh, threadmodel, inflation, _blockRecursion, userCache)

  # and the factor data itself
  return FunctionNodeData{typeof(ccw)}(eliminated, potentialused, edgeIDs, ccw, multihypo, ccw.certainhypo, nullhypo, solveInProgress, inflation)
end

"""
    $SIGNATURES

Return `::Bool` on whether at least one hypothesis is available for intended computations (assuming direction `sfidx`).
"""
function isLeastOneHypoAvailable( sfidx::Int,
                                  certainidx::Vector{Int},
                                  uncertnidx::Vector{Int},
                                  isinit::Vector{Bool})
  #
  # @show isinit
  # @show sfidx in certainidx, sum(isinit[uncertnidx])
  # @show sfidx in uncertnidx, sum(isinit[certainidx])
  return  sfidx in certainidx && 0 < sum(isinit[uncertnidx]) ||
          sfidx in uncertnidx && sum(isinit[certainidx]) == length(certainidx)
end

function assembleFactorName(dfg::AbstractDFG,
                            Xi::Vector{<:DFGVariable} )
  #

  existingFactorLabels = listFactors(dfg)
  existingFactorLabelDict = Dict(existingFactorLabels .=> existingFactorLabels)
  namestring = ""
  for vert in Xi #f.Xi
    namestring = string(namestring,vert.label)
  end
  opt = getSolverParams(dfg)
  for i in 1:opt.maxincidence
    tempnm = string(namestring, "f$i")
    if !haskey(existingFactorLabelDict, Symbol(tempnm))
      namestring = tempnm
      break
    end
    i != opt.maxincidence ? nothing : error("Artificial restriction to not connect more than $(opt.maxincidence) factors to a variable (bad for sparsity), try setting getSolverParams(fg).maxincidence=1000 to adjust this restriction.")
  end
  return Symbol(namestring)
end

"""
    $(SIGNATURES)

Add factor with user defined type `<:AbstractFactor`` to the factor graph
object. Define whether the automatic initialization of variables should be
performed.  Use order sensitive `multihypo` keyword argument to define if any
variables are related to data association uncertainty.

Experimental
- `inflation`, to better disperse kernels before convolution solve, see IIF #1051.
"""
function DFG.addFactor!(dfg::AbstractDFG,
                        Xi::AbstractVector{<:DFGVariable},
                        usrfnc::AbstractFactor;
                        multihypo::Vector{Float64}=Float64[],
                        nullhypo::Float64=0.0,
                        solvable::Int=1,
                        tags::Vector{Symbol}=Symbol[],
                        timestamp::Union{DateTime,ZonedDateTime}=now(localzone()),
                        graphinit::Bool=getSolverParams(dfg).graphinit,
                        threadmodel=SingleThreaded,
                        suppressChecks::Bool=false,
                        inflation::Real=getSolverParams(dfg).inflation,
                        namestring::Symbol = assembleFactorName(dfg, Xi),
                        _blockRecursion::Bool=!getSolverParams(dfg).attemptGradients  )
  #

  varOrderLabels = Symbol[v.label for v=Xi]
  solverData = getDefaultFactorData(dfg,
                                    Xi,
                                    deepcopy(usrfnc);
                                    multihypo,
                                    nullhypo,
                                    threadmodel,
                                    inflation,
                                    _blockRecursion)
  #
  newFactor = DFGFactor(Symbol(namestring),
                        varOrderLabels,
                        solverData;
                        tags=Set(union(tags, [:FACTOR])),
                        solvable,
                        timestamp)
  #

  success = addFactor!(dfg, newFactor)

  # TODO: change this operation to update a conditioning variable
  graphinit && doautoinit!(dfg, Xi, singles=false)

  return newFactor
end

function _checkFactorAdd(usrfnc, xisyms)
  if length(xisyms) == 1 && !(usrfnc isa AbstractPrior) && !(usrfnc isa Mixture)
    @warn("Listing only one variable $xisyms for non-unary factor type $(typeof(usrfnc))")
  end
  nothing
end

function DFG.addFactor!(dfg::AbstractDFG,
                        xisyms::AbstractVector{Symbol},
                        usrfnc::AbstractFactor;
                        suppressChecks::Bool=false,
                        kw...  )
  #

  # basic sanity check for unary vs n-ary
  if !suppressChecks
    _checkFactorAdd(usrfnc, xisyms)
  end

  # variables = getVariable.(dfg, xisyms)
  variables = map(vid -> getVariable(dfg, vid), xisyms)
  addFactor!(dfg, variables, usrfnc; suppressChecks=suppressChecks, kw... )
end

#


text:
Pkg.activate("./docs")

using Documenter, DeviationsLH, FilesLH

makedocs(
    modules = [DeviationsLH],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    checkdocs = :exports,
    sitename = "DeviationsLH",
    pages = Any["index.md"]
)

pkgDir = rstrip(normpath(@__DIR__, ".."), '/');
@assert endswith(pkgDir, "DeviationsLH")
deploy_docs(pkgDir);

Pkg.activate(".")

# ------------

ext: jl


text:
function edges(e::Ellipse{T,T,Nothing}; n = 4) where {T}
    φs = range(0, stop = 2π, length = n + 1)
    pts = [CartesianPoint(CylindricalPoint{T}(e.r, φ, zero(T))) for φ in φs]
    pts = map(p -> _transform_into_global_coordinate_system(p, e), pts)
    edges = [Edge(pts[i], pts[i+1]) for i in 1:n]
end
function edges(e::Ellipse{T,T,T}; n = 4) where {T}
    φs = range(T(0), stop = e.φ, length = n + 1)
    pts = [CartesianPoint(CylindricalPoint{T}(e.r, φ, zero(T))) for φ in φs]
    pts = map(p -> _transform_into_global_coordinate_system(p, e), pts)
    edges = [Edge(pts[i], pts[i+1]) for i in 1:n]
end
function edges(e::Ellipse{T,NTuple{2,Tuple{T}},Nothing}; n = 4) where {T}
    φs = range(0, stop = 2π, length = n + 1)
    pts = [CartesianPoint{T}(e.r[1][1] * cos(φ), e.r[2][1] * sin(φ), zero(T)) for φ in φs]
    pts = map(p -> _transform_into_global_coordinate_system(p, e), pts)
    edges = [Edge(pts[i], pts[i+1]) for i in 1:n]
end

@recipe function f(e::Ellipse; n = 40)

    @series begin
        label --> "Ellipse"
        edges(e, n = n)
    end
end


text:
# module MultivariateCalculus

export vech, trilind, triuind,
  commutation, spcommutation,
  duplication, spduplication,
  chol_gradient, chol_gradient!,
  kron_gradient, kron_gradient!, kronaxpy!,
  bump_diagonal!, clamp_diagonal!

"""
Vectorize the lower triangular part of a matrix.
"""
function vech(a::Union{Number, AbstractVecOrMat})
  m, n = size(a)
  out = similar(a, convert(Int, (2m - n + 1) * n / 2))
  ooffset, aoffset = 1, 1
  for j = 1:n
    len = m - j + 1 # no. elements to copy in column j
    copy!(out, ooffset, a, aoffset, len)
    ooffset += m - j + 1
    aoffset += m + 1
  end
  out
end
# a is a scalar or (column) vector
vech(a::Union{Number, AbstractVector}) = copy(a)

"""
    commutation(type, m[, n])

Create the `mn x mn` commutation matrix `K`, defined by
`K * vec(A) = vec(A')` for any `m x n` matrix A.
"""
function commutation(t::Type, m::Integer, n::Integer)
  ((m < 0) || (n < 0)) && throw(ArgumentError("invalid Array dimensions"))
  mn = m * n
  reshape(kron(vec(eye(t, m)), eye(t, n)), mn, mn)
end
commutation(m::Integer, n::Integer) = commutation(Float64, m, n)
commutation(t::Type, m::Integer) = commutation(t, m, m)
commutation(m::Integer) = commutation(m, m)
commutation(M::AbstractMatrix) = commutation(eltype(M), size(M, 1), size(M, 2))

"""
    spcommutation(type, m[, n])

Create the sparse `mn`-by-`mn` commutation matrix `K`, defined by
`K * vec(A) = vec(A')` for any `m x n` matrix A.
"""
function spcommutation(t::Type, m::Integer, n::Integer)
  ((m < 0) || (n < 0)) && throw(ArgumentError("invalid Array dimensions"))
  mn = m * n
  reshape(kron(vec(speye(t, m)), speye(t, n)), mn, mn)
end
spcommutation(m::Integer, n::Integer) = spcommutation(Float64, m, n)
spcommutation(t::Type, m::Integer) = spcommutation(t, m, m)
spcommutation(m::Integer) = spcommutation(m, m)
spcommutation(M::AbstractMatrix) = spcommutation(eltype(M), size(M, 1), size(M, 2))

"""
    trilind(m, n,[ k])

Linear indices of the lower triangular part of an `m x n` array.
"""
function trilind(m::Integer, n::Integer, k::Integer)
  find(tril(trues(m, n), k))
end
function trilind(m::Integer, n::Integer)
  find(tril(trues(m, n)))
end
function trilind(m::Integer)
  find(tril(trues(m, m)))
end
trilind(M::AbstractArray) = trilind(size(M, 1), size(M, 2))
trilind(M::AbstractArray, k::Integer) = trilind(size(M, 1), size(M, 2), k)

"""
    triuind(m, n,[ k])

Linear indices of the upper triangular part of an `m x n` array.
"""
function triuind(m::Integer, n::Integer, k::Integer)
  find(triu(trues(m, n), k))
end
function triuind(m::Integer, n::Integer)
  find(triu(trues(m, n)))
end
function triuind(m::Integer)
  find(triu(trues(m, m)))
end
triuind(M::AbstractArray) = triuind(size(M, 1), size(M, 2))
triuind(M::AbstractArray, k::Integer) = triuind(size(M, 1), size(M, 2), k)

"""
    spduplication(type, n)

Create the sparse `n^2 x n(n+1)/2` duplication matrix, defined by
`D * vech(A) = vec(A)` for any symmetric matrix.
"""
function spduplication(t::Type, n::Integer)
  imatrix = zeros(typeof(n), n, n)
  imatrix[trilind(n, n)] = 1:binomial(n + 1, 2)
  imatrix = imatrix + tril(imatrix, -1).'
  sparse(1:n^2, vec(imatrix), one(t))
end
spduplication(n::Integer) = spduplication(Float64, n)
spduplication(M::AbstractMatrix) = spduplication(eltype(M), size(M, 1))
duplication(t::Type, n::Integer) = full(spduplication(t, n))
duplication(n::Integer) = duplication(Float64, n)
duplication(M::AbstractMatrix) = duplication(eltype(M), size(M, 1))

"""
    kron_gradient!(g, dM, Y)

Compute the gradient `d / d vec(X)` from a vector of derivatives `dM` where
`M=X⊗Y`, `n, q = size(X)`, and `p, r = size(Y)`.
"""
function kron_gradient!(g::VecOrMat{T}, dM::VecOrMat{T},
  Y::Matrix{T}, n::Integer, q::Integer) where {T <: Real}
  p, r = size(Y)
  A_mul_B!(g, kron(speye(n * q), vec(Y)'),
    (kron(speye(q), spcommutation(n, r), speye(p)) * dM))
end
function kron_gradient(dM::VecOrMat{T}, Y::Matrix{T},
  n::Integer, q::Integer) where {T <: Real}
  if ndims(dM) == 1
    g = zeros(T, n * q)
  else
    g = zeros(T, n * q, size(dM, 2))
  end
  kron_gradient!(g, dM, Y, n, q)
end

"""
    chol_gradient!(g, dM, L)

Compute the gradient `d / d vech(L)` from a vector of derivatives `dM` where
`M=L*L'`.
# TODO make it more memory efficient
"""
function chol_gradient!(g::AbstractVecOrMat{T},
  dM::AbstractVecOrMat{T}, L::AbstractMatrix{T}) where {T <: Real}
  n = size(L, 1)
  At_mul_B!(g, spduplication(n),
    kron(L', speye(n)) * (dM + spcommutation(n) * dM))
end
function chol_gradient(dM::AbstractVecOrMat{T}, L::AbstractMatrix{T}) where {T <: Real}
  n = size(L, 1)
  if ndims(dM) == 1 # vector
    g = zeros(T, binomial(n + 1, 2))
  else # matrix
    g = zeros(T, binomial(n + 1, 2), size(dM, 2))
  end
  chol_gradient!(g, dM, L)
end

"""
    kronaxpy!(A, X, Y)

Overwrites `Y` with `A ⊗ X + Y`. Same as `Y += kron(A, X)` but more efficient.
"""
function kronaxpy!(A::AbstractVecOrMat{T},
  X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) where {T <: Real}

  # retrieve matrix sizes
  m, n = size(A, 1), size(A, 2)
  p, q = size(X, 1), size(X, 2)
  # loop over (i,j) blocks of Y
  irange, jrange = 1:p, 1:q
  @inbounds for j in 1:n, i in 1:m
    a = A[i, j]
    irange = ((i - 1) * p + 1):(i * p)
    jrange = ((j - 1) * q + 1):(j * q)
    Yij = view(Y, irange, jrange)  # view of (i, j)-block
    @simd for k in eachindex(Yij)
      Yij[k] += a * X[k]
    end
  end
  Y
end

"""
Add `ϵ` to the diagonal entries of matrix `A`.
"""
function bump_diagonal!(A::Matrix{T}, ϵ::T) where {T}
  @inbounds @simd for i in 1:minimum(size(A))
    A[i, i] += ϵ
  end
  A
end

"""
Clamp the diagonal entries of matrix `A` to `[lo, hi]`.
"""
function clamp_diagonal!(A::Matrix{T}, lo::T, hi::T) where {T}
  @inbounds @simd for i in 1:minimum(size(A))
    A[i, i] = clamp(A[i, i], lo, hi)
  end
  A
end

#end # module MultivariateCalculus


text:
module PlotsTests

include("imgcomp.jl")

# don't actually show the plots
srand(1234)
default(show=false, reuse=true)
img_eps = 5e-2

facts("Gadfly") do
    @fact gadfly() --> Plots.GadflyBackend()
    @fact backend() --> Plots.GadflyBackend()

    @fact typeof(plot(1:10)) --> Plots.Plot{Plots.GadflyBackend}
    @fact plot(Int[1,2,3], rand(3)) --> not(nothing)
    @fact plot(sort(rand(10)), rand(Int, 10, 3)) --> not(nothing)
    @fact plot!(rand(10,3), rand(10,3)) --> not(nothing)

    image_comparison_facts(:gadfly, skip=[4,6,19,23,24,27], eps=img_eps)
end

facts("PyPlot") do
    @fact pyplot() --> Plots.PyPlotBackend()
    @fact backend() --> Plots.PyPlotBackend()

    image_comparison_facts(:pyplot, skip=[19], eps=img_eps)
end

facts("GR") do
    @fact gr() --> Plots.GRBackend()
    @fact backend() --> Plots.GRBackend()

    @linux_only image_comparison_facts(:gr, skip=[24], eps=img_eps)
end

facts("Plotly") do
    @fact plotly() --> Plots.PlotlyBackend()
    @fact backend() --> Plots.PlotlyBackend()

    # # until png generation is reliable on OSX, just test on linux
    # @linux_only image_comparison_facts(:plotly, only=[1,3,4,7,8,9,10,11,12,14,15,20,22,23,27], eps=img_eps)
end

# facts("Immerse") do
#     @fact immerse() --> Plots.ImmerseBackend()
#     @fact backend() --> Plots.ImmerseBackend()
#
#     # as long as we can plot anything without error, it should be the same as Gadfly
#     image_comparison_facts(:immerse, only=[1], eps=img_eps)
# end

# facts("PlotlyJS") do
#     @fact plotlyjs() --> Plots.PlotlyJSBackend()
#     @fact backend() --> Plots.PlotlyJSBackend()
#
#     # as long as we can plot anything without error, it should be the same as Plotly
#     image_comparison_facts(:plotlyjs, only=[1], eps=img_eps)
# end

facts("UnicodePlots") do
    @fact unicodeplots() --> Plots.UnicodePlotsBackend()
    @fact backend() --> Plots.UnicodePlotsBackend()

    # lets just make sure it runs without error
    @fact isa(plot(rand(10)), Plot) --> true
end

FactCheck.exitstatus()
end # module


text:
using PyCall

# Change that to whatever packages you need.
const PACKAGES = ["numpy", "cma"]

# Use eventual proxy info
proxy_arg=String[]
if haskey(ENV, "http_proxy")
    push!(proxy_arg, "--proxy")
    push!(proxy_arg, ENV["http_proxy"])
end

# Import pip
if !PyCall.pyexists("pip")
    # If it is not found, install it
    println("Pip not found on your system. Downloading it.")
    get_pip = joinpath(dirname(@__FILE__), "get-pip.py")
    download("https://bootstrap.pypa.io/get-pip.py", get_pip)
    run(`$(PyCall.python) $(proxy_arg) $get_pip --user`)
end

println("Installing required python packages using pip")
run(`$(PyCall.python) $(proxy_arg) -m pip install --user $(PACKAGES)`)


text:
include("prop_logic_theory.jl")
include("prover.jl")

using Test

# Metatheory.options.printiter = true
# Metatheory.options.verbose = true

ex = rewrite(:(((p => q) ∧ (r => s) ∧ (p ∨ r)) => (q ∨ s)), impl)
@test prove(t, ex, 3, 10, 5000)

@test @areequal t true ((¬p == p) == false)
@test @areequal t true ((¬p == ¬p) == true)
@test @areequal t true ((¬p ∨ ¬p) == ¬p) (¬p ∨ p) ¬(¬p ∧ p)
@test @areequal t true ((p => (p ∨ p)))
@test @areequal t true ((p => (p ∨ p)) == ((¬(p) ∧ q) => q)) == true

# Frege's theorem
@test @areequal t true (p => (q => r)) => ((p => q) => (p => r))

# Demorgan's
@test @areequal t true (¬(p ∨ q) == (¬p ∧ ¬q))

# Consensus theorem
@test @areequal t ((x ∧ y) ∨ (¬x ∧ z) ∨ (y ∧ z))   ((x ∧ y) ∨ (¬x ∧ z))


text:
module Seqeval

export NETag, get_tags, f1score, SeqEval

include("./tags.jl")

function curr_isstart(prev::NETag, curr::NETag)
    (curr.loc == B_Tag || curr.loc == S_Tag) && return true
    (curr.loc == E_Tag || curr.loc == I_Tag) &&
        (prev.loc == E_Tag || prev.loc == S_Tag || prev.loc == O_Tag) && return true
    curr.loc != O_Tag && prev.type != curr.type && return true
    return false
end

function prev_isend(prev::NETag, curr::NETag)
    (prev.loc == E_Tag || prev.loc == S_Tag) && return true
    (prev.loc == B_Tag || prev.loc == I_Tag) &&
        (curr.loc == B_Tag || curr.loc == S_Tag || curr.loc == O_Tag) && return true
    prev.loc != O_Tag && prev.type != curr.type && return true
    return false
end

function entity_ranges(tags::AbstractVector{NETag})
    result = Array{Tuple{NEType, UnitRange}}(undef, 0)
    sizehint!(result, length(tags))
    entity_ranges!(result, tags)
    sizehint!(result, length(result))
    return result
end

function entity_ranges!(result, tags::AbstractVector{NETag})
    len = length(tags)

    @assert len > 0
    prev = first(tags)
    base, i = 1, 2
    @inbounds while i <= len
        tag = tags[i]

        prev_isend(prev, tag) && push!(result, (prev.type, base:i-1))
        curr_isstart(prev, tag) && (base = i)

        prev = tag
        i += 1
    end
    prev_isend(prev, nulltag) && push!(result, (prev.type, base:i-1))
    return result
end

function plus1!(d::Dict, key)
    d[key] = get(d, key, 0) + 1
end

function compute_correct!(result, p::AbstractVector{<:AbstractVector{NETag}}, gt::AbstractVector{<:AbstractVector{NETag}})
    for (pi, gti) in zip(p, gt)
        compute_correct!(result, pi, gti)
    end
    return result
end

compute_correct(p::AbstractVector{<:AbstractVector{NETag}}, gt::AbstractVector{<:AbstractVector{NETag}}) = compute_correct!(SeqEval(), p, gt)

compute_correct(p::AbstractVector{NETag}, gt::AbstractVector{NETag}) = compute_correct!(SeqEval(), p, gt)
function compute_correct!(result, p::AbstractVector{NETag}, gt::AbstractVector{NETag})
    plen = length(p)
    gtlen = length(gt)
    maxlen = max(plen, gtlen)

    @assert plen > 0 && gtlen > 0
    p_prev = first(p)
    gt_prev = first(gt)
    pbase = gtbase = 1
    i = 2
    c = p_prev == gt_prev
    @inbounds while i <= maxlen
        ptag = i > plen ? nulltag : p[i]
        gttag = i > gtlen ? nulltag : gt[i]

        ptag == gttag && (c += 1)

        p_isend = prev_isend(p_prev, ptag)
        gt_isend = prev_isend(gt_prev, gttag)

        p_isstart = curr_isstart(p_prev, ptag)
        gt_isstart = curr_isstart(gt_prev, gttag)

        p_isend && plus1!(result.p, p_prev.type)
        gt_isend && plus1!(result.gt, gt_prev.type)

        if p_isend && gt_isend && pbase == gtbase && p_prev.type == gt_prev.type
            plus1!(result.correct, gt_prev.type)
        end

        p_isstart && (pbase = i)
        gt_isstart && (gtbase =  i)

        p_prev = ptag
        gt_prev = gttag
        i += 1
    end

    p_isend = prev_isend(p_prev, nulltag)
    gt_isend = prev_isend(gt_prev, nulltag)

    p_isend && plus1!(result.p, p_prev.type)
    gt_isend && plus1!(result.gt, gt_prev.type)

    if p_isend && gt_isend && pbase == gtbase && p_prev.type == gt_prev.type
        plus1!(result.correct, gt_prev.type)
    end

    result.count[] += c
    result.total[] += maxlen

    return result
end

function micro_average(result, β)
    p = sum(values(result.p))
    gt = sum(values(result.gt))
    correct = sum(values(result.correct))

    precision = iszero(p) ? 0 : correct / p
    recall = iszero(gt) ? 0 : correct / gt
    acc = result.count[] / result.total[]

    if isinf(β) && β > 0
        f1 = recall
    else
        denominator = iszero(precision) && iszero(recall) ? one(recall) : recall + precision * β^2
        f1 = (1 + β^2) * recall * precision / denominator
    end
    return (accuracy = acc, precision = precision, recall = recall, f1score = f1)
end

function precision_recall_f1score_support(p, gt; β = true)
    result = compute_correct(p, gt)
    micro_average(result, β)
end

accuracy(p, gt) = precision_recall_f1score_support(p, gt).accuracy
precision(p, gt) = precision_recall_f1score_support(p, gt).precision
recall(p, gt) = precision_recall_f1score_support(p, gt).recall
f1score(p, gt; β = true) = precision_recall_f1score_support(p, gt; β).f1score

struct SeqEval
    correct::Dict{NEType,Int}
    p::Dict{NEType, Int}
    gt::Dict{NEType, Int}
    count::Ref{Int}
    total::Ref{Int}
end

SeqEval() = SeqEval(Dict{NEType,Int}(), Dict{NEType, Int}(), Dict{NEType, Int}(), Ref(0), Ref(0))

Base.push!(s::SeqEval, p, gt) = compute_correct!(s, p, gt)

precision_recall_f1score_support(s::SeqEval; β = true) = micro_average(s, β)
accuracy(s::SeqEval) = precision_recall_f1score_support(s).accuracy
precision(s::SeqEval) = precision_recall_f1score_support(s).precision
recall(s::SeqEval) = precision_recall_f1score_support(s).recall
f1score(s::SeqEval; β = true) = precision_recall_f1score_support(s; β).f1score

# y_true = [["O", "O", "O", "B-MISC", "I-MISC", "I-MISC", "O"], ["B-PER", "I-PER", "O"]]
# y_pred = [["O", "O", "B-MISC", "I-MISC", "I-MISC", "I-MISC", "O"], ["B-PER", "I-PER", "O"]]
# y_true = [["O", "O", "B-MISC", "I-MISC", "B-MISC", "O", "O"], ["B-PER", "I-PER", "O"]]
# y_pred = [["O", "O", "B-MISC", "I-MISC", "B-MISC", "I-MISC", "O"], ["B-PER", "I-PER", "O"]]

end


text:
include("Yfind.jl")
include("is_it.jl")

ODE2(z,w)=conj(  im * z.*( 1 ./(w.^2-z.^2)+1 ./(1+z.^2) ))
ODE1(z,w)=conj(  im * w.*( 1 ./(z.^2-w.^2)+1 ./(1+w.^2) ))
function Eq_of_M(du,u,p,t)
    du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[5]=0
return
end
condition(u,t,integrator)= u[5]>max_hit
function condition2(u,t,integrator) # Event when event_f(u,t) == 0
u[1]
end
affect!(integrator) = terminate!(integrator)
function affect2!(integrator)
integrator.u[5]=integrator.u[5]+1
end
cb2 = ContinuousCallback(condition2,affect2!,nothing)
cb1 = DiscreteCallback(condition,affect!)

cb=CallbackSet(cb2,cb1)
H_test(u)=( (u[3]-u[1])^2+(u[2]-u[4])^2 )*( (u[3]+u[1])^2+(u[2]+u[4])^2 )/ ((u[2]^4+2*u[2]^2*(u[1]^2-1)+(1+u[1]^2)^2 )*(u[3]^4+2*u[3]^2*(u[4]^2+1)+(u[4]^2-1)^2 ))

function escape_exit_function(Q,P, Energy, t_end, max_hit)

    tol_dist=1e-5
    H=(2*Energy)^2

    Q_0=0
    P_0=0

    Q_1=0
    P_1=0

    Q_2=0
    P_2=0

    Q_3=0
    P_3=0

    Y=Yfind(Q,P,H)
    if ~isempty(Y)
        u0=zeros(5)
        u0[1]=0 #X
        u0[2]=P #P
        u0[3]=Q #Q
        u0[4]=Y #Y
        u0[5]=0 #Y
        if true
            prob = ODEProblem(Eq_of_M,u0,(0., t_end))
            # sol=solve(prob,RK4(),reltol=1e-6,abstol=1e-8,callback=cb)
            sol=solve(prob, Vern9(),maxiters=1e10, reltol=1e-13,abstol=1e-15,callback=cb)
            uf=zeros(5)
            uf[1]=sol[1,end] #X
            uf[2]=sol[2,end] #P
            uf[3]=sol[3,end] #Q
            uf[4]=sol[4,end] #Y
            uf[5]=0
            dH=abs(H_test(uf)-H)
            #
            if dH<1e-5
                if is_it(uf,H,tol_dist)
                    return sol.u[end][5]
                else
                    return -1
                end
            end
        end
    end
    return -2
end


text:
# This file contains the types and definitions for the
# downsampled lippmann-Schwinger equation
using Interpolations

type FastDualDownSampled
    # type to encapsulate the fast application of M = I + omega^2G*spadiagm(nu)
    # Fast:: FastMslow
    Fast:: FastMslowDual
    n  :: Int64
    m  :: Int64
    # downsampling factor in each dimension
    downsampleX :: Int64
    downsampleY :: Int64
    # Sampling Matrix
    S :: SparseMatrixCSC{Float64,Int64}

    function FastDualDownSampled(Fast, downsampleX, downsampleY)
        # Building the Sampling matrix
        index1 = 1:downsampleX:Fast.n
        index2 = 1:downsampleY:Fast.m
        Sindex = spzeros(n,m)
        for ii = 0:round(Integer,(n-1)/downsampleX)
            for jj = 0:round(Integer,(m-1)/downsampleY)
                Sindex[1+ii*downsampleX,1+jj*downsampleY ] = 1
            end
        end

        index = find(Sindex[:])
        S = speye(n*m, n*m)
        Sampling = S[index,:];
        return new( Fast, Fast.n, Fast.m, downsampleX, downsampleY,Sampling)
    end
end

# we need to overload these functions for gmres
import Base.*
import Base.A_mul_B!
import Base.size

function size(M::FastDualDownSampled, dim::Int64)
  # function to returns the size of the underliying matrix i.e. the rank of M.S
  return size(M.S,1)
end

function A_mul_B!(Y,
                  M::FastDualDownSampled,
                  V)
    # in place matrix matrix multiplication
    @assert(size(Y) == size(V))
    # print(size(V))
    for ii = 1:size(V,2)
        Y[:,ii] = M*V[:,ii]
    end
end

function *(M::FastDualDownSampled, b::Array{Complex128,1})
    knots = (collect(1:M.downsampleX:M.n), collect(1:M.downsampleY:M.m))
    nDown = round(Integer, (M.n-1)/M.downsampleX +1)
    mDown = round(Integer, (M.m-1)/M.downsampleY +1)
    B = reshape(b, nDown,mDown);

    itp_real = interpolate(real(B), BSpline(Quadratic(Reflect())), OnGrid())
    itp_imag = interpolate(imag(B), BSpline(Quadratic(Reflect())), OnGrid())

    interpU =     itp_real[collect(M.downsampleX:M.n+M.downsampleX-1)/M.downsampleX,
                           collect(M.downsampleY:M.m+M.downsampleY-1)/M.downsampleY] +
              1im*itp_imag[collect(M.downsampleX:M.n+M.downsampleX-1)/M.downsampleX,
                           collect(M.downsampleY:M.m+M.downsampleY-1)/M.downsampleY];

    u = M.Fast*interpU[:]
    size(u)
    # Subsampling the wavefield
    return (M.S*u)[:]
end

max_stars_repo_name: Forgotten/Lippmann-Schingwer

max_issues_repo_name: Forgotten/Lippmann-Schingwer

max_forks_repo_name: Forgotten/Lippmann-Schingwer


text:
using Plots

x = collect(1:100)
y1 = rand(100)
y1p = rand(100)
y1pp = rand(100)
y2 = rand(100)
y3 = rand(100)
y4 = rand(100)

p1 = plot(x, y1)
plot!(x, y1p)
plot!(x, y1pp)
p2 = plot(x, y2)
p3 = plot(x, y3)
p4 = plot(x, y4)

plot(p1, p2, p3, p4, layout=4)

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
# After first running snoop_on(), open a new julia session and run this

if !isinteractive()
    error("Because Gadfly uses the display, you must run this interactively")
end

using SnoopCompile, Gadfly

SnoopCompile.@snoop1 "/tmp/gadfly_compiles.csv" begin
    include(Pkg.dir("Gadfly", "test","runtests.jl"))
end

snoop_off()

data = SnoopCompile.read("/tmp/gadfly_compiles.csv")
pc, discards = SnoopCompile.parcel(data[end:-1:1,2])
SnoopCompile.write("/tmp/precompile", pc)


text:
import Distributed

"""
    addprocs_existing(x; kwargs...)

Equivalent to `addprocs(ExistingProcessManager(x); kwargs...)`.

## Examples

### Example 1

```julia
julia> hosts_and_ports = [
       ("127.0.0.1", 9601), # the host is "127.0.0.1", and the port number is 9601
       ("127.0.0.1", 9602), # the host is "127.0.0.1", and the port number is 9602
       ]
2-element Vector{Tuple{String, Int64}}:
 ("127.0.0.1", 9601)
 ("127.0.0.1", 9602)

julia> addprocs_existing(hosts_and_ports; kwargs...)
```

### Example 2

```julia
julia> w1 = WorkerConfig();

julia> w1.host = "127.0.0.1";

julia> w1.port = 9601;

julia> w2 = WorkerConfig();

julia> w2.host = "127.0.0.1";

julia> w2.port = 9602;

julia> wconfigs = WorkerConfig[w1, w2];

julia> addprocs_existing(wconfigs; kwargs...)
```

### Example 3

```julia
julia> worker_output = \"""
       julia_worker:9960#192.168.1.151
       julia_worker:9960#192.168.1.151
       julia_worker:9960#192.168.1.151
       julia_worker:9966#192.168.1.157
       julia_worker:9968#192.168.1.159
       \""";

julia> addprocs_existing(worker_output; kwargs...)
```
"""
@inline function addprocs_existing(x; kwargs...)
    return Distributed.addprocs(ExistingProcessManager(x); kwargs...)
end


text:
# Core structure and functions for decision diagrams
# This structure stores both node and arc layers in arrays as fields of DD
# It is fast to compute longest path, but (possibly?!) slow in removing nodes and arcs!
# Commented out sections are not adjusted...

# All higher level functions should use these functions instead of accessing
# the internal structure of a decision diagram directly

# Functions in this file must be careful to keep graph and layers synchronized

using LightGraphs
#using ResumableFunctions
#import LightGraphs.SimpleGraphs.SimpleEdge
import LightGraphs.inneighbors

"""Node structure"""
mutable struct Node{S<:AbstractFloat}
	layer::Int				# node layer
	state::S				# node state
#	out_neighbor::Array{Int,1}		# outgoing arcs of the node, stores arc id in the arc layer
#	in_neighbor::Array{Int,1}		# outgoing arcs of the node, stores arc id in the arc layer
	out_neighbor::Dict{Int, Array{Int, 1}}		# key: the id of node neighbor in the next layer, value: set of arcs (id) connected to that node
	in_neighbor::Dict{Int, Array{Int, 1}}		# key: the id of node neighbor in the previous layer, value: set of arcs (id) connected to that node
	aux_val::Float64		# auxiliary value at the node e.g., the longest path value
end

"""Decision diagram structure"""
struct DecisionDiagram{S<:AbstractFloat}
	graph::SimpleDiGraph{Int}            # represents the directed graph of the DD
    node_layers::Array{Array{Node{S},1},1}        # stores instances of the Node class in each node layer
	arc_layers::Array{Array{Tuple{Int, Int, Float64},1},1}		# stores arc info as tuples (tail, head, label) in arc layers
end

"""
A constrcter for type `DecisionDiagram`

# Input
- `t::Type`: Data type of the state values of nodes
- `nvars::Int`: Number of variables (arc layers)
"""

DecisionDiagram(nvars::Int; t::DataType = Float64)::DecisionDiagram{t} = DecisionDiagram{t}(SimpleDiGraph{Int}(), [Array{Node{t},1}(  0) for i=1:nvars+1], [Array{Tuple{Int, Int, Float64},1}(  0) for i=1:nvars])

"""Return the number of variables (or number of arc layers) in the decision diagram"""
nvars(dd::DecisionDiagram)::Int = length(dd.node_layers) - 1

#nnodes(dd::DecisionDiagram) = nv(dd.graph)
nnodes(dd::DecisionDiagram)::Int = sum(length(dd.node_layers[i]) for i = 1:length(dd.node_layers))

#narcs(dd::DecisionDiagram) = ne(dd.graph)
narcs(dd::DecisionDiagram)::Int = sum(length(dd.arc_layers[i]) for i = 1:length(dd.arc_layers))

"""Set state of an existing node"""
function set_node_state!(dd::DecisionDiagram, layer::Int, id::Int, state::S) where S<:AbstractFloat
	dd.node_layers[layer][id].state = state
end

"""Get state of a node"""
get_node_state(dd::DecisionDiagram, layer::Int, id::Int)::S where S<:AbstractFloat = dd.node_layers[layer][id].state

"""Return arc label"""
get_arc_label(dd::DecisionDiagram, layer::Int, id::Int)::Float64 = dd.arc_layers[layer][id][3]
#get_arc_labels(dd::DecisionDiagram, arc::LightGraphs.SimpleGraphs.SimpleEdge{Int64}) = collect(dd.nodes[src(e)].out_neighbor[dst(e)])

"""Return arc tail and head"""
get_arc_tail(dd::DecisionDiagram, layer::Int, id::Int)::Int = dd.arc_layers[layer][id][1]
get_arc_head(dd::DecisionDiagram, layer::Int, id::Int)::Int = dd.arc_layers[layer][id][2]

get_node_layer_size(dd::DecisionDiagram, layer::Int)::Int = length(dd.node_layers[layer])
get_arc_layer_size(dd::DecisionDiagram, layer::Int)::Int = length(dd.arc_layers[layer])

get_width(dd::DecisionDiagram)::Int = maximum(length(dd.node_layers[i]) for i=1:length(dd.node_layers))

#=
function get_arc_labels(dd::DecisionDiagram, node1::Int, node2::Int)
	i = findfirst(dd.nodes[node1].out_neighbor_nodes, node2)
	#if i == 0
	#	error("No arc exists between node ", node1, " and node ", node2)
	#else
		return dd.nodes[node1].out_neighbor_nodes[i]
	#end
end
=#

"""Return the array of incoming arcs of a node"""
inneighbors(dd::DecisionDiagram, layer::Int, id::Int)::Dict{Int, Array{Int, 1}} = dd.node_layers[layer][id].in_neighbor
#inneighbors(dd::DecisionDiagram, layer::Int, id::Int)::Array{Int, 1} = dd.node_layers[layer][id].in_neighbor

"""Return the array of outgoing nodes of a node"""
outneighbors(dd::DecisionDiagram, layer::Int, id::Int)::Dict{Int, Array{Int, 1}} = dd.node_layers[layer][id].out_neighbor
#outneighbors(dd::DecisionDiagram, layer::Int, id::Int)::Array{Int, 1} = dd.node_layers[layer][id].out_neighbor

"""Add a new node to the decision diagram with a given state and return the node's id at its layer"""
function add_node!(dd::DecisionDiagram, layer::Int, state::S)::Int where S<:AbstractFloat
	#g = dd.graph
	#add_vertex!(g)
	#aux_node::Node{S} = Node(layer, state, Array{Int, 1}(0), Array{Int,1}(0), 0.0)		#This is for the case where neighbor field stores array of arcs only
	aux_node::Node{S} = Node(layer, state, Dict{Int, Array{Int, 1}}(), Dict{Int, Array{Int, 1}}(), 0.0)
	push!(dd.node_layers[layer], aux_node)
	id::Int = length(dd.node_layers[layer])
	return id
end

#=
# TODO: rem_node! is untested
"""Remove a node from the decision diagram"""
function rem_node!(dd::DecisionDiagram, node::Int)

	# ensure the node belongs to the vertex set
	# node in vertices(g) || return false
	node <= length(dd.nodes) || return false

	# Remove incoming and outgoing arcs of the node
	for outnode in keys(dd.nodes[node].out_neighbor)
		delete!(dd.nodes[outnode].in_neighbor, node)
	end
	for innode in keys(dd.nodes[node].in_neighbor)
		delete!(dd.nodes[innode].out_neighbor, node)
	end

	# Move end node in the layer to node's position within the layer
	layer = dd.nodes[node].layer
	layerid = dd.nodes[node].layerid
	dd.layers[layer][layerid] = pop!(dd.layers[layer])
	dd.nodes[dd.layers[layer][layerid]].layerid = layerid

	# Renumber the last node of the graph to the removing node according to LightGraphs implementation
	lastid = length(dd.nodes)
	dd.nodes[node] = pop!(dd.nodes)			# swap the last node with the removed node
	dd.layers[dd.nodes[node].layer][dd.nodes[node].layerid] = node			# update the last node number in dd.layers
	# Update the head/tail number in the incoming/outgoing arcs of the last node
	for outnode in keys(dd.nodes[node].out_neighbor)
		dd.nodes[outnode].in_neighbor[node] = dd.nodes[outnode].in_neighbor[lastid]
		delete!(dd.nodes[outnode].in_neighbor, lastid)
	end
	for innode in keys(dd.nodes[node].in_neighbor)
		dd.nodes[innode].out_neighbor[node] = dd.nodes[innode].out_neighbor[lastid]
		delete!(dd.nodes[innode].out_neighbor, lastid)
	end

	# Remove vertex from graph structure
	#removed = rem_vertex!(dd.graph, node)
	#@assert(removed)
end
=#

#=
"""Remove a node from the decision diagram given its layer and layerid"""
function rem_node!(dd::DecisionDiagram, layer::Int, layerid::Int)
	rem_node!(dd, dd.layers[layer][layerid])
end
=#

"""Add a new arc to the decision diagram between two nodes and retunrs arc id in its layer"""
function add_arc!(dd::DecisionDiagram, layer1::Int, node1_id::Int, layer2::Int, node2_id::Int, label::Float64)::Int
	#@assert((dd.nodes[node1].layer < dd.nodes[node2].layer) && (node1 <= length(dd.nodes)) && (node2 <= length(dd.nodes)))
	#g = dd.graph
	#e = Edge(node1, node2)
	@assert(layer1 == layer2 - 1)
#	if isa(label, Int)					# converts Int label type to floating point
#		lbl = float(label)
#	else
#		lbl = label
#	end
	push!(dd.arc_layers[layer1], (node1_id, node2_id, label))
	id::Int = length(dd.arc_layers[layer1])
#	push!(dd.node_layers[layer1][node1_id].out_neighbor, id)	#for case where neighbors are array of arcs only
#	push!(dd.node_layers[layer2][node2_id].in_neighbor, id)	#for case where neighbors are array of arcs only
	if haskey(dd.node_layers[layer1][node1_id].out_neighbor, node2_id)
		push!(dd.node_layers[layer1][node1_id].out_neighbor[node2_id], id)
		push!(dd.node_layers[layer2][node2_id].in_neighbor[node1_id], id)
	else
		dd.node_layers[layer1][node1_id].out_neighbor[node2_id] = [id]
		dd.node_layers[layer2][node2_id].in_neighbor[node1_id] = [id]
	end
	return id
end

"""
	Add a new arc to the decision diagram between two nodes, where the arc_merging is on. Returns the arc id in the layer.
	If the arc is parallel, it does not add the arc if its label is not lower or upper bound among other arcs.
"""
#NOTE: this code considers the case where inneighbor/outneighbor field is a Dict (NOT array)
function add_reduced_arc!(dd::DecisionDiagram, layer1::Int, node1_id::Int, layer2::Int, node2_id::Int, label::Float64)::Int
	#g = dd.graph
	#e = Edge(node1, node2)
	id::Int = 0			# if the arc is not added, the returning id will be zero
	# check if an arc between two nodes already exists
	if !haskey(dd.node_layers[layer1][node1_id].out_neighbor, node2_id)
		#add_edge!(g, e)
		push!(dd.arc_layers[layer1], (node1_id, node2_id, label))
		id = length(dd.arc_layers[layer1])
		dd.node_layers[layer1][node1_id].out_neighbor[node2_id] = [id]
		dd.node_layers[layer2][node2_id].in_neighbor[node1_id] = [id]
	# if there is only a single arc, the arc is added to in/outneighbor in an increasing order wrt label value
	elseif length(dd.node_layers[layer1][node1_id].out_neighbor[node2_id]) == 1
		push!(dd.arc_layers[layer1], (node1_id, node2_id, label))
		id = length(dd.arc_layers[layer1])
		if dd.arc_layers[layer1][dd.node_layers[layer1][node1_id].out_neighbor[node2_id][1]][3] <= label
			 push!(dd.node_layers[layer1][node1_id].out_neighbor[node2_id], id)	# adds the new label to the parallel arcs
			 push!(dd.node_layers[layer2][node2_id].in_neighbor[node1_id], id)
		 else
			 pushfirst!(dd.node_layers[layer1][node1_id].out_neighbor[node2_id], id)	# adds the new label to the parallel arcs
 			 pushfirst!(dd.node_layers[layer2][node2_id].in_neighbor[node1_id], id)
 		 end
	# if there are two arcs between the nodes, the new arc is added only if it does not belong to the label interval of the current arcs
	# To add the new arc, we only need to change the label of the arc that is being replaced
	else
		@assert(length(dd.node_layers[layer1][node1_id].out_neighbor[node2_id]) == 2)
		id1 = dd.node_layers[layer1][node1_id].out_neighbor[node2_id][1]
		id2 = dd.node_layers[layer1][node1_id].out_neighbor[node2_id][2]
		if label < dd.arc_layers[layer1][id1][3]
			id = id1
			dd.arc_layers[layer1][id] = (node1_id, node2_id, label)
		elseif label > dd.arc_layers[layer1][id2][3]
			id = id2
			dd.arc_layers[layer1][id] = (node1_id, node2_id, label)
		end
	end
	return id
end

#=
# TODO rem_arc! is untested
"""Remove a labeled arc from a decision diagram"""
function rem_arc!(dd::DecisionDiagram, node1::Int, node2::Int, label::Int)
	# checks if the arc exists
	if haskey(dd.nodes[node1].out_neighbor, node2)
		setdiff!(dd.nodes[node1].out_neighbor[node2], label)
		setdiff!(dd.nodes[node2].in_neighbor[node1], label)
		if isempty(dd.nodes[node1].out_neighbor[node2])
			#e = Edge(node1, node2)
			#rem_edge!(dd.graph, e)
			delete!(dd.nodes[node1].out_neighbor, node2)
			delete!(dd.nodes[node2].in_neighbor, node1)
		end
	end
end
=#

#=
# TODO rem_arcs! is untested
"""Remove all (parallel) arcs between two nodes from a decision diagram"""
function rem_arcs!(dd::DecisionDiagram, node1::Int, node2::Int)
	# checks if the arc exists
	if haskey(dd.nodes[node1].out_neighbor, node2)
		#e = Edge(node1, node2)
		#rem_edge!(dd.graph, e)
		delete!(dd.nodes[node1].out_neighbor, node2)
		delete!(dd.nodes[node2].in_neighbor, node1)
	end
end
=#

max_stars_repo_name: ctjandra/ddopt-ecp

max_issues_repo_name: ctjandra/ddopt-ecp

max_forks_repo_name: ctjandra/ddopt-ecp


text:
# default branch fill strategy
# returns the tuple:
# the number of 1st dim slices and the number of child nodes in each slice
function omt_branch_fill(tree::RTree; fill_factor::Number=1.0) where T<:Node
    @assert 0.0 <= fill_factor <= 1.0
    maxlen = ceil(Int, capacity(Branch, tree) * fill_factor)
    slicelen = ceil(Int, 0.75*sqrt(maxlen))
    nslices = fld1(maxlen, slicelen)
    return (nslices, slicelen)
end

"""
    load!(tree::RTree{T,N,V}, data::Any;
          convertel = identity, method = :OMT,
          leaf_fill = capacity(Leaf, tree),
          branch_fill::Tuple{Integer, Integer} = omt_branch_fill(tree)) where {T,N,V}

Bulk-load `data` into `tree`.
  * `tree`: an *empty* R-tree for storing elements of type `V`
  * `data`: iterable container with the elements to put into `tree`
  * `convertel`: function to convert elements of `data` to type `V`
  * `method`: bulk-loading method
  * `leaf_fill`: the average number of elements to store in R-tree leaves (1-level nodes)
  * `branch_fill`: the tuple of the number of slices and the number of subtrees per slice
    in the R-tree nodes (level ≥ 1).

The supported bulk-loading methods are:
  * `:OMT`: *Overlap Minimizing Top-down method* by Taewon Lee and Sukho Lee
"""
function load!(tree::RTree{T,N,V}, data::Any;
               convertel = identity,
               method::Symbol = :OMT,
               leaf_fill::Integer = capacity(Leaf, tree),
               branch_fill::Tuple{Integer, Integer} = omt_branch_fill(tree)) where {T, N, V}
    # load from arbitrary data into Elem-backed R-tree
    isempty(tree) || throw(ArgumentError("Cannot bulk-load into non-empty tree"))
    if isempty(data)
        @warn "No bulk-load data provided"
        return tree
    end
    0 < leaf_fill <= capacity(Leaf, tree) ||
        throw(ArgumentError("Leaf fill should be positive and not exceed leaf capacity"))
    prod(branch_fill) > 1 || throw(ArgumentError("Branch fill should be > 1"))
    if method == :OMT
        return load_omt!(tree, V[convertel(x) for x in data],
                         leaf_fill=leaf_fill, branch_fill=branch_fill)
    else
        throw(ArgumentError("Unsupported bulk-loading method $method"))
    end
end

# See "OMT: Overlap Minimizing Top-down Bulk Loading Algorithm for R-tree"
# by Taewon Lee and Sukho Lee, http://ceur-ws.org/Vol-74/files/FORUM_18.pdf
# Modified to take into account that leaf and branch capacity is different
# and allow branch slices at any tree level (specified by `branch_fill`)
function load_omt!(tree::RTree, elems::AbstractVector;
                   leaf_fill::Integer = capacity(Leaf, tree),
                   branch_fill::Tuple{Integer, Integer} = omt_branch_fill(tree))
    # calculate the tree properties
    nbranch_subtrees = min(branch_fill[1] * branch_fill[2], capacity(Branch, tree))
    height = max(ceil(Int, (log(length(elems)) - log(leaf_fill)) / log(nbranch_subtrees)), 0) + 1 # 1 for leaves
    if height > 1
        maxelems_subtree = leaf_fill*nbranch_subtrees^(height-2) # max elements in each root subtree
        nroot_subtrees = fld1(length(elems), maxelems_subtree)
        nroot_slices = floor(Int, sqrt(nroot_subtrees))
    else # root === single leaf
        maxelems_subtree = nroot_subtrees = length(elems)
        nroot_slices = 1
    end
    #@show length(elems) leaf_fill branch_fill height maxelems_subtree nsubtrees nslices
    # sort by the center of the first dim
    dim = mod1(height, ndims(tree)) # root level dimension
    elems_sorted = sort(elems, by = elem -> mbr(elem).low[dim] + mbr(elem).high[dim])
    resize!(tree.nnodes_perlevel, height)
    fill!(tree.nnodes_perlevel, 0)
    tree.root = omt_subtree(elems_sorted, tree, height,
                            nroot_slices, nroot_subtrees, leaf_fill, branch_fill)
    tree.nnodes_perlevel[end] = 1
    return tree
end

# put the elems into the subtree of height lev
# the root node of the subtree should have nsubtree children
# organized into nslices slices along the current elems order
# recursively calls itself
function omt_subtree(elems::AbstractVector, tree::RTree,
                     lev::Integer, nslices::Integer, nsubtrees::Integer,
                     leaf_fill::Integer, branch_fill::Tuple{Integer, Integer})
    @debug "omt_subtree(): lev=$lev nslices=$nslices nsubtrees=$nsubtrees nelems=$(length(elems))"
    @assert lev > 0
    if length(elems) <= nsubtrees
        # if fewer elements than the number of subtrees,
        # then all elements should be put into single leaf
        if lev == 1 # create a Leaf and attach all elements
            node = acquire(tree, Leaf, lev)
            for elem in elems
                _attach!(node, elem, tree)
            end
            tree.nelems += length(elems)
            return node
        else # not leaf level yet, create a branch with a single child
            child = omt_subtree(elems, tree, lev - 1, nslices, nsubtrees,
                                leaf_fill, branch_fill)
            node = _attach!(acquire(tree, Branch, lev), child, tree)
            tree.nnodes_perlevel[lev - 1] += 1
            return node
        end
    end

    # subtree root
    @assert lev > 1
    node = acquire(tree, Branch, lev)

    # create subtrees
    nelems_subtree = fld1(length(elems), nsubtrees) # actual fill of the subtree
    nelems_slice = nelems_subtree * fld1(nsubtrees, nslices)
    dim = mod1(lev, ndims(tree)) # cycle sorting dimensions through tree levels
    for i in 1:nelems_slice:length(elems) # slice using the external elements order (sort-dim of the level above)
        # sort the elements of the slice by the current dim
        elems_slice = sort(view(elems, i:min(i+nelems_slice-1, length(elems))),
                           by = elem -> mbr(elem).low[dim] + mbr(elem).high[dim])
        # create slice subtrees and attach to the node
        for j in 1:nelems_subtree:length(elems_slice)
            subtree_elems_range = j:min(j+nelems_subtree-1, length(elems_slice))
            child = omt_subtree(view(elems_slice, subtree_elems_range), tree,
                                lev - 1, branch_fill[1],
                                min(branch_fill[1] * branch_fill[2], capacity(Branch, tree)),
                                leaf_fill, branch_fill)
            _attach!(node, child, tree)
            tree.nnodes_perlevel[lev - 1] += 1
        end
    end
    return node
end


text:
using Test
using Diderot

@testset "internals" begin
    include("internals.jl")
end

@testset "Knapsack" begin
    include("knapsack.jl")
end

@testset "Set Cover" begin
    include("setcover.jl")
end

@testset "Index Fund" begin
    include("indexfund.jl")
end


text:
using Base: Indices, tail

abstract type DataSeries{T,N} <: AbstractArray{T,N} end

for (TDataSeries, TArray) in
    ((:SDataSeries, :Array),
     (:PDataSeries, :SharedArray))
    @eval begin
        struct $TDataSeries{T,N} <: DataSeries{T,N}
            nd::Int
            nt::Int
            ni::Int
            d::$TArray{T,N}

            function $TDataSeries{T,N}(nd, nt, ni) where {T,N}
                @assert T <: Number
                @assert nd > 0
                @assert nt ≥ 0
                @assert ni > 0

                @assert N ∈ (1,2,3)

                if N == 1
                    d = $TArray{T}(undef, nt+1)
                elseif N == 2
                    d = $TArray{T}(undef, nd, nt+1)
                elseif N == 3
                    d = $TArray{T}(undef, nd, nt+1, ni)
                end

                fill!(d,zero(T))

                new(nd, nt, ni, d)
            end
        end

        function $TDataSeries(T, nt::Int)
            return $TDataSeries{T,1}(1, nt, 1)
        end

        function $TDataSeries(T, nd::Int, nt::Int)
            return $TDataSeries{T,2}(nd, nt, 1)
        end

        function $TDataSeries(T, nd::Int, nt::Int, ni::Int)
            ni == 1 ? N = 2 : N = 3
            return $TDataSeries{T,N}(nd, nt, ni)
        end

        function $TDataSeries(d::AbstractArray{T,1}) where {T}
            nd = 1
            nt = size(d,1)-1
            ni = 1
            ds = $TDataSeries{T,1}(nd, nt, ni)
            copy!(ds.d, d[1,:,1])
            return ds
        end

        function $TDataSeries(d::AbstractArray{T,2}) where {T}
            nd = size(d,1)
            nt = size(d,2)-1
            ni = 1
            ds = $TDataSeries{T,2}(nd, nt, ni)
            copy!(ds.d, d[:,:,1])
            return ds
        end

        function $TDataSeries(d::AbstractArray{T,3}) where {T}
            nd = size(d,1)
            nt = size(d,2)-1
            ni = size(d,3)
            ds = $TDataSeries{T,3}(nd, nt, ni)
            copy!(ds.d, d)
            return ds
        end
    end
end

Base.:(==)(ds1::DataSeries, ds2::DataSeries) = (
                                ds1.nd == ds2.nd
                             && ds1.nt == ds2.nt
                             && ds1.ni == ds2.ni
                             && ds1.d  == ds2.d)

function Base.show(io::IO, ds::DataSeries{T,N}) where {T,N}
    print(io, "DataSeries with data type ", T, " and ", N, " dimensions:\n")
    print(io, ds.d)
end

function similar(ds::DataSeries{T,N}) where {T,N}
    typeof(ds){T,N}(ds.nd, ds.nt, ds.ni)
end

Base.parent(ds::DataSeries) = ds.d
Base.eltype(ds::DataSeries{T,N}) where {T,N} = T
Base.ndims(ds::DataSeries{T,N}) where {T,N} = N

# errmsg_size(ds::DataSeries) = error("size not supported for data series with axes $(axes(ds))")
# Base.size(ds::DataSeries) = errmsg_size(ds)
# Base.size(ds::DataSeries, d) = errmsg_size(ds)
Base.size(ds::DataSeries) = size(ds.d)
Base.size(ds::DataSeries, d) = size(ds.d, d)

Base.eachindex(::IndexCartesian, ds::DataSeries) = CartesianIndices(axes(ds))
Base.eachindex(::IndexLinear, ds::DataSeries) = axes(ds, 1)

Base.firstindex(ds::DataSeries{T,1}) where {T}   = 0
Base.firstindex(ds::DataSeries{T,N}) where {T,N} = 1

Base.firstindex(ds::DataSeries{T,1}, d) where {T} = d ≥ 1 && d ≤ 1 ? 0 : 1
Base.firstindex(ds::DataSeries{T,2}, d) where {T} = d ≥ 1 && d ≤ 2 ? (1,0)[d] : 1
Base.firstindex(ds::DataSeries{T,3}, d) where {T} = d ≥ 1 && d ≤ 3 ? (1,0,1)[d] : 1

Base.lastindex(ds::DataSeries{T,1}) where {T} = ds.nt
Base.lastindex(ds::DataSeries{T,2}) where {T} = ds.nd*(ds.nt+1)
Base.lastindex(ds::DataSeries{T,3}) where {T} = ds.nd*(ds.nt+1)*ds.ni

Base.lastindex(ds::DataSeries{T,1}, d) where {T} = d ≥ 1 && d ≤ 1 ? ds.nt : 1
Base.lastindex(ds::DataSeries{T,2}, d) where {T} = d ≥ 1 && d ≤ 2 ? (ds.nd, ds.nt)[d] : 1
Base.lastindex(ds::DataSeries{T,3}, d) where {T} = d ≥ 1 && d ≤ 3 ? (ds.nd, ds.nt, ds.ni)[d] : 1

@inline Base.axes(ds::DataSeries{T,1}) where {T} = (0:ds.nt,)
@inline Base.axes(ds::DataSeries{T,2}) where {T} = (1:ds.nd, 0:ds.nt)
@inline Base.axes(ds::DataSeries{T,3}) where {T} = (1:ds.nd, 0:ds.nt, 1:ds.ni)
@inline Base.axes(ds::DataSeries{T,N}, d) where {T,N} = d ≥ 1 && d ≤ N ? axes(ds)[d] : (1:1)

Base.strides(ds::DataSeries) = strides(ds.d)

function get_data!(ds::DataSeries{T,1}, n, k=1) where {T}
    @assert n ≥ 0 && n ≤ ds.nt
    @assert k == 1
    @inbounds return ds.d[n+1]
end

function get_data!(ds::DataSeries{T,2}, x::Union{Array{T,1}, Array{TwicePrecision{T},1}}, n, k=1) where {T}
    @assert size(x,1) == ds.nd
    @assert n ≥ 0 && n ≤ ds.nt
    @assert k == 1
    j = n+1
    @inbounds for i in eachindex(x)
        x[i] = ds.d[i,j]
    end
end

function get_data!(ds::DataSeries{T,3}, x::Union{Array{T,2}, Array{TwicePrecision{T},2}}, n) where {T}
    @assert size(x,1) == ds.nd
    @assert size(x,2) == ds.ni
    @assert n ≥ 0 && n ≤ ds.nt
    j = n+1
    @inbounds for k in axes(x, 2)
        for i in axes(x, 1)
            x[i,k] = ds.d[i,j,k]
        end
    end
end

function get_data!(ds::DataSeries{T,3}, x::Union{Array{T,1}, Array{TwicePrecision{T},1}}, n, k) where {T}
    @assert size(x,1) == ds.nd
    @assert n ≥ 0 && n ≤ ds.nt
    @assert k ≥ 1 && k ≤ ds.ni
    j = n+1
    @inbounds for i in eachindex(x)
        x[i] = ds.d[i,j,k]
    end
end

function set_data!(ds::DataSeries{T,1}, x::Union{T, TwicePrecision{T}}, n, k=1) where {T}
    @assert n ≥ 0 && n ≤ ds.nt
    @assert k == 1
    j = n+1
    @inbounds ds.d[j] = x[i]
end

function set_data!(ds::DataSeries{T,2}, x::Union{Array{T,1}, Array{TwicePrecision{T},1}}, n, k=1) where {T}
    @assert size(x,1) == ds.nd
    @assert n ≥ 0 && n ≤ ds.nt
    @assert k == 1
    j = n+1
    @inbounds for i in axes(ds.d, 1)
        ds.d[i,j] = x[i]
    end
end

function set_data!(ds::DataSeries{T,3}, x::Union{Array{T,2}, Array{TwicePrecision{T},2}}, n) where {T}
    @assert size(x,1) == ds.nd
    @assert size(x,2) == ds.ni
    @assert n ≥ 0 && n ≤ ds.nt
    j = n+1
    @inbounds for k in axes(ds.d, 3)
        for i in axes(ds.d, 1)
            ds.d[i,j,k] = x[i,k]
        end
    end
end

function set_data!(ds::DataSeries{T,3}, x::Union{Array{T,1}, Array{TwicePrecision{T},1}}, n, k) where {T}
    @assert size(x,1) == ds.nd
    @assert n ≥ 0 && n ≤ ds.nt
    @assert k ≥ 1 && k ≤ ds.ni
    j = n+1
    @inbounds for i in axes(ds.d, 1)
        ds.d[i,j,k] = x[i]
    end
end

function CommonFunctions.reset!(ds::DataSeries{T,1}) where {T}
    @inbounds ds[0] = ds[end]
end

function CommonFunctions.reset!(ds::DataSeries{T,2}) where {T}
    @inbounds for i in axes(ds, 1)
        ds[i,0] = ds[i,end]
    end
end

function CommonFunctions.reset!(ds::DataSeries{T,3}) where {T}
    @inbounds for k in axes(ds.d, 3)
        for i in axes(ds.d, 1)
            ds[i,0,k] = ds[i,end,k]
        end
    end
end

@inline function Base.getindex(ds::DataSeries{T,1}, j::Union{UnitRange,Int}) where {T}
    @boundscheck checkbounds(ds.d, j.+1)
    @inbounds r = getindex(ds.d, j.+1)
    return r
end

@inline function Base.getindex(ds::DataSeries{T,2}, i::Union{UnitRange,Int}, j::Union{UnitRange,Int}) where {T}
    @boundscheck checkbounds(ds.d, i, j.+1)
    @inbounds r = getindex(ds.d, i, j.+1)
    return r
end

@inline function Base.getindex(ds::DataSeries{T,2}, j::Union{UnitRange,Int}) where {T}
    @boundscheck checkbounds(ds.d, :, j.+1)
    @inbounds r = getindex(ds.d, :, j+1)
    return r
end

@inline function Base.getindex(ds::DataSeries{T,3}, i::Union{UnitRange,Int}, j::Union{UnitRange,Int}, k::Union{UnitRange,Int}) where {T}
    @boundscheck checkbounds(ds.d, i, j.+1, k)
    @inbounds r = getindex(ds.d, i, j.+1, k)
    return r
end

@inline function Base.getindex(ds::DataSeries{T,3}, j::Union{UnitRange,Int}, k::Union{UnitRange,Int}) where {T}
    @boundscheck checkbounds(ds.d, :, j.+1, k)
    @inbounds r = getindex(ds.d, :, j.+1, k)
    return r
end

@inline function Base.getindex(ds::DataSeries{T,3}, k::Union{UnitRange,Int}) where {T}
    @boundscheck checkbounds(ds.d, :, :, k)
    @inbounds r = getindex(ds.d, :, :, k)
    return r
end

@inline function Base.setindex!(ds::DataSeries{T,1}, x, j::Union{UnitRange,Int}) where {T}
    @assert length(x) == 1
    @boundscheck checkbounds(ds.d, j.+1)
    @inbounds setindex!(ds.d, x, j.+1)
end

@inline function Base.setindex!(ds::DataSeries{T,2}, x, i::Union{UnitRange,Int}, j::Union{UnitRange,Int}) where {T}
    @assert length(x) == 1
    @boundscheck checkbounds(ds.d, i, j.+1)
    @inbounds setindex!(ds.d, x, i, j.+1)
end

@inline function Base.setindex!(ds::DataSeries{T,2}, x, j::Union{UnitRange,Int}) where {T}
    @assert ndims(x) == 1
    @assert length(x) == size(ds.d, 1)
    @boundscheck checkbounds(ds.d, :, j.+1)
    @inbounds setindex!(ds.d, x, :, j.+1)
end

@inline function Base.setindex!(ds::DataSeries{T,3}, x, i::Union{UnitRange,Int}, j::Union{UnitRange,Int}, k::Union{UnitRange,Int}) where {T}
    @assert length(x) == 1
    @boundscheck checkbounds(ds.d, i, j.+1, k)
    @inbounds setindex!(ds.d, x, i, j.+1, k)
end

@inline function Base.setindex!(ds::DataSeries{T,3}, x, j::Union{UnitRange,Int}, k::Union{UnitRange,Int}) where {T}
    @assert ndims(x) == 1
    @assert length(x) == size(ds.d, 1)
    @boundscheck checkbounds(ds.d, :, j.+1, k)
    @inbounds setindex!(ds.d, x, :, j.+1, k)
end

@inline function Base.setindex!(ds::DataSeries{T,3}, x, k::Union{UnitRange,Int}) where {T}
    @assert ndims(x) == 2
    @assert size(x,1) == size(ds.d, 1)
    @assert size(x,2) == size(ds.d, 2)
    @boundscheck checkbounds(ds.d, :, :, k)
    @inbounds setindex!(ds.d, x, :, :, k)
end

# TODO Implement convert() to/from array.


text:
function testSMC(model::SMCModel, N::Int64, numTrials::Int64)
  smcio = SMCIO{model.particle, model.pScratch}(N, model.maxn)
  println("Running SMC. N = $N")
  for i=1:numTrials
    @time smc!(model, smcio)
    left = max(length(smcio.logZhats)-10,1)
    right = length(smcio.logZhats)
    println(smcio.logZhats[left:right])
  end
end

function testCSMC(model::SMCModel, N::Int64, numTrials::Int64)
  smcio = SMCIO{model.particle, model.pScratch}(N, model.maxn)
  println("Running CSMC. N = $N")
  smc!(model, smcio)
  v::Vector{model.particle} = Vector{model.particle}(undef, smcio.n)
  for i in 1:smcio.n
    v[i] = model.particle()
  end
  SimpleSMC.pickParticle!(v, smcio)
  for i=1:numTrials
    @time csmc!(model, smcio, v, v)
    left = max(length(smcio.logZhats)-10,1)
    right = length(smcio.logZhats)
    println(smcio.logZhats[left:right])
  end
end

function testCSMCBS(model::SMCModel, lM::F, N::Int64, numTrials::Int64) where F<:Function
  smcio = SMCIO{model.particle, model.pScratch}(N, model.maxn)
  println("Running CSMC-BS. N = $N")
  smc!(model, smcio)
  v::Vector{model.particle} = Vector{model.particle}(undef, smcio.n)
  for i in 1:smcio.n
    v[i] = model.particle()
  end
  SimpleSMC.pickParticle!(v, smcio)
  for i=1:numTrials
    @time csmc!(model, lM, smcio, v, v)
    left = max(length(smcio.logZhats)-10,1)
    right = length(smcio.logZhats)
    println(smcio.logZhats[left:right])
  end
end


text:
module DemoModule
    using ExportPublic

    function add(a, b)
        # `add` is exported automatically
        return _private_helper(a) + _private_helper(b)
    end

    function _private_helper(a)
        # `_private_helper` should not be exported, so it's marked with a `_`
        return a
    end

    @exportPublic()
end


text:
isdefined(Base, :__precompile__) && __precompile__()

"""
Main module for QuantizedStateSystems.jl
"""
module QuantizedStateSystems

  using MacroTools
  using TaylorSeries

  export @model

  include("model.jl")
  include("quantizer.jl")
  include("integrator.jl")

end


text:
using DiffEqJump, DiffEqBase, OrdinaryDiffEq
using Test
# using StableRNGs
# rng = StableRNG(12345)

# k10, G -> G*
# k01, G* -> G
# ρ,  G -> G + P
# d, P -> 0

u0 = [0,1,0]
reactstoch = [[1=>1],[2=>1],[1=>1],[3=>1]]
netstoch = [[1=>-1,2=>1],[1=>1,2=>-1],[3=>1],[3=>-1]]
rates = [.1,0.25,10.,1.2]
prob = DiscreteProblem(u0,(0.0,200.0))
majumps = MassActionJump(rates,reactstoch,netstoch)

jump_prob_MNRM = JumpProblem(prob,MNRM(),majumps,save_positions=(true,false))
jump_prob_NRM = JumpProblem(prob,NRM(),majumps,save_positions=(false,false))
jump_prob_Direct = JumpProblem(prob,Direct(),majumps,save_positions=(true,false))

integrator = DiffEqBase.__init(jump_prob_Direct,SSAStepper(),seed=1234)
aggregator = integrator.cb.initialize
# aggregator(integrator)
aggregator
# integrator.sol.t
integrator.tstop

sol =@time solve(jump_prob_Direct,SSAStepper())

sol_NRM =@time solve(jump_prob_NRM,SSAStepper(),saveat=1.)
sol.u
sol_NRM.u

ens_sol_MNRM =@time solve(EnsembleProblem(jump_prob_MNRM), SSAStepper(), EnsembleThreads(); trajectories=5e4)
ens_sol_NRM =@time solve(EnsembleProblem(jump_prob_NRM), SSAStepper(), EnsembleThreads(); trajectories=5e4)
ens_sol_Direct =@time solve(EnsembleProblem(jump_prob_Direct), SSAStepper(), EnsembleThreads(); trajectories=5e4)

using DifferentialEquations.EnsembleAnalysis
end_P = componentwise_vectors_timepoint(ens_sol_MNRM,200.)[3]
end_P_D = componentwise_vectors_timepoint(ens_sol_Direct,200.)[3]

using Plots
histogram(end_P, bins = 0:1:20, normalize = :pdf)
histogram!(end_P_D, bins = 0:1:20, normalize = :pdf)


text:
export FiniteStrainContinuum

"""
    FiniteStrainContinuum

Implements the finite strain element.

- `eledim`: spatial dimension of the element (default = 2).
- `mat`: constitutive law, a length `#elem` vector of materials such as [`PlaneStress`](@ref)
- `elnodes`: the node indices in this finite element, an integer array
- `props`: property dictionary
- `coords`: coordinates of the vertices of the element
- `dhdx`, `weights`, `hs`: data for integral
- `stress`: stress at each quadrature points
"""
mutable struct FiniteStrainContinuum
    eledim::Int64
    mat  # constitutive law
    elnodes::Array{Int64}   # the node indices in this finite element
    props::Dict{String, Any}
    coords::Array{Float64}
    dhdx::Array{Array{Float64}}  # 4nPointsx2 matrix
    weights::Array{Float64}
    hs::Array{Array{Float64}}
    stress::Array{Array{Float64}}
end

function FiniteStrainContinuum(coords::Array{Float64}, elnodes::Array{Int64}, props::Dict{String, Any}, ngp::Int64=2)
    eledim = 2
    dhdx, weights, hs = get2DElemShapeData( coords, ngp )
    nGauss = length(weights)
    name = props["name"]
    if name=="PlaneStrain"
        mat = [PlaneStrain(props) for i = 1:nGauss]
    elseif name=="PlaneStress"
        mat = [PlaneStress(props) for i = 1:nGauss]
    elseif name=="PlaneStressPlasticity"
        mat = [PlaneStressPlasticity(props) for i = 1:nGauss]
    elseif name=="PlaneStressIncompressibleRivlinSaunders"
        mat = [PlaneStressIncompressibleRivlinSaunders(props) for i = 1:nGauss]
    elseif name=="NeuralNetwork2D"
        mat = [NeuralNetwork2D(props) for i = 1:nGauss]
    else
        error("Not implemented yet: $name")
    end
    strain = Array{Array{Float64}}(undef, length(weights))
    FiniteStrainContinuum(eledim, mat, elnodes, props, coords, dhdx, weights, hs, strain)
end

function getStiffAndForce(self::FiniteStrainContinuum, state::Array{Float64}, Dstate::Array{Float64}, Δt::Float64)
    ndofs = dofCount(self);
    nnodes = length(self.elnodes)
    fint = zeros(Float64, ndofs)
    stiff = zeros(Float64, ndofs,ndofs)
    out = Array{Float64}[]
    u = state[1:nnodes]; v = state[nnodes+1:2*nnodes]
    Du = Dstate[1:nnodes]; Dv = Dstate[nnodes+1:2*nnodes]
    for k = 1:length(self.weights)
        g1 = self.dhdx[k][:,1]; g2 = self.dhdx[k][:,2]

        ux = u'*g1; uy = u'*g2; vx = v'*g1; vy = v'*g2
        Dux = Du'*g1; Duy = Du'*g2; Dvx = Dv'*g1; Dvy = Dv'*g2
        # compute  ∂E∂u.T, 8 by 3 array
        ∂E∂u = [g1+ux*g1 uy*g2    g2 + g2*ux+g1*uy;
                vx*g1    g2+vy*g2 g1 + g1*vy+g2*vx;]

        E = [ux+0.5*(ux*ux+vx*vx); vy+0.5*(uy*uy+vy*vy); uy+vx+ux*uy+vx*vy]
        DE = [Dux+0.5*(Dux*Dux+Dvx*Dvx); Dvy+0.5*(Duy*Duy+Dvy*Dvy); Duy+Dvx+Dux*Duy+Dvx*Dvy]

        # #@show "+++",E, DE
        S, dS_dE = getStress(self.mat[k], E, DE, Δt)
        # error()
        self.stress[k] = S

        fint += ∂E∂u * S * self.weights[k] # 1x8

        S∂∂E∂∂u = [g1*g1'*S[1]+g2*g2'*S[2]+g1*g2'*S[3]+g2*g1'*S[3] zeros(nnodes, nnodes);
                  zeros(nnodes, nnodes)  g1*g1'*S[1]+g2*g2'*S[2]+g1*g2'*S[3]+g2*g1'*S[3]]

        stiff += (∂E∂u * dS_dE * ∂E∂u' + S∂∂E∂∂u)*self.weights[k] # 8x8
    end
    return fint, stiff
end

function getInternalForce(self::FiniteStrainContinuum, state::Array{Float64}, Dstate::Array{Float64}, Δt::Float64)
    n = dofCount(self)
    fint = zeros(Float64,n)
    out = Array{Float64}[]
    u = state[1:4]; v = state[5:8]
    Du = Dstate[1:4]; Dv = Dstate[5:8]
    for k = 1:length(self.weights)
        g1 = self.dhdx[k][:,1]; g2 = self.dhdx[k][:,2]
        ux = u'*g1; uy = u'*g2; vx = v'*g1; vy = v'*g2
        Dux = Du'*g1; Duy = Du'*g2; Dvx = Dv'*g1; Dvy = Dv'*g2
        # compute  ∂E∂u.T, 8 by 3 array
        ∂E∂u = [g1+ux*g1 uy*g2    g2 + g2*ux+g1*uy;
                vx*g1    g2+vy*g2 g1 + g1*vy+g2*vx;]

        E = [ux+0.5*(ux*ux+vx*vx); vy+0.5*(uy*uy+vy*vy); uy+vx+ux*uy+vx*vy]
        DE = [Dux+0.5*(Dux*Dux+Dvx*Dvx); Dvy+0.5*(Duy*Duy+Dvy*Dvy); Duy+Dvx+Dux*Duy+Dvx*Dvy]
        S, _ = getStress(self.mat[k],E, DE, Δt)

        self.stress[k] = S

        fint += ∂E∂u * S * self.weights[k] # 1x8
    end
    return fint
end

function getMassMatrix(self::FiniteStrainContinuum)
    ndofs = dofCount(self)
    nnodes = length(self.elnodes)
    mass = zeros(ndofs,ndofs)
    for k = 1:length(self.weights)
        rho = self.mat[k].ρ
        mass += [self.hs[k]*self.hs[k]' zeros(nnodes, nnodes)
                 zeros(nnodes, nnodes)  self.hs[k]*self.hs[k]']  * rho * self.weights[k]
    end
    lumped = sum(mass, dims=2)
    mass, lumped
end

function getStrain(self::FiniteStrainContinuum, state::Array{Float64})
    ndofs = dofCount(self);
    nnodes = length(self.elnodes)
    fint = zeros(Float64, ndofs)
    stiff = zeros(Float64, ndofs,ndofs)
    out = Array{Float64}[]
    u = state[1:nnodes]; v = state[nnodes+1:2*nnodes]

    nGauss = length(self.weights)
    E = zeros(nGauss, 3)
    w∂E∂u = zeros(nGauss, 2nnodes, 3)
    for k = 1:length(self.weights)
        g1 = self.dhdx[k][:,1]; g2 = self.dhdx[k][:,2]

        ux = u'*g1; uy = u'*g2; vx = v'*g1; vy = v'*g2
        # compute  ∂E∂u.T, 8 by 3 array
        ∂E∂u = [g1+ux*g1 uy*g2    g2 + g2*ux+g1*uy;
                vx*g1    g2+vy*g2 g1 + g1*vy+g2*vx;]  # 8x3
        E[k,:] = [ux+0.5*(ux*ux+vx*vx); vy+0.5*(uy*uy+vy*vy); uy+vx+ux*uy+vx*vy] # 3x1
        w∂E∂u[k,:,:] = ∂E∂u*self.weights[k] # 8x3
    end

    return E, w∂E∂u
end

function getNodes(self::FiniteStrainContinuum)
    return self.elnodes
end

function dofCount(self::FiniteStrainContinuum)
    return 2length(self.elnodes)
end

function commitHistory(self::FiniteStrainContinuum)
    for m in self.mat
        commitHistory(m)
    end
end

function  getForceAndDforceDstress(self::FiniteStrainContinuum, state::Array{Float64},
    stress::Array{Float64,2})
    ndofs = dofCount(self);
    nnodes = length(self.elnodes)
    nStrain = 3
    nGauss = length(self.weights)
    fint = zeros(Float64, ndofs)
    dfint_dstress = zeros(Float64,  ndofs, nGauss * nStrain)

    u = state[1:nnodes]; v = state[nnodes+1:2*nnodes]

    for k = 1:length(self.weights)
        g1 = self.dhdx[k][:,1]; g2 = self.dhdx[k][:,2]

        ux = u'*g1; uy = u'*g2; vx = v'*g1; vy = v'*g2
        # compute  ∂E∂u, 3 by 2nnodes array
        ∂E∂u = [g1+ux*g1 uy*g2    g2 + g2*ux+g1*uy;
                vx*g1    g2+vy*g2 g1 + g1*vy+g2*vx;]

        S = stress[k, :]

        self.stress[k] = S

        fint += ∂E∂u * S * self.weights[k] # 1x8

        dfint_dstress[:, (k-1)*nStrain+1:k*nStrain] = ∂E∂u * self.weights[k]

    end
    return fint , dfint_dstress
end


text:
using Test
using LinearAlgebra, SparseArrays, Random
using LuxurySparse

Random.seed!(2)

p1 = IMatrix{4}()
sp = sprand(ComplexF64, 4, 4, 0.5)
ds = rand(ComplexF64, 4, 4)
pm = PermMatrix([2, 3, 4, 1], randn(4))
v = [0.5, 0.3im, 0.2, 1.0]
dv = Diagonal(v)

@testset "invdet" begin
    ####### linear algebra  ######
    @test inv(p1) == inv(Matrix(p1))
    @test det(p1) == det(Matrix(p1))
    @test diag(p1) == diag(Matrix(p1))
    @test logdet(p1) == 0
    @test inv(pm) == inv(Matrix(pm))

    for m in [pm, sp, p1, dv]
        @test !(m |> isdense)
        @test !(m' |> isdense)
        @test !(transpose(m) |> isdense)
    end
    for m in [ds, v]
        @test m |> isdense
        @test m' |> isdense
        @test transpose(m) |> isdense
    end
end

@testset "multiply" begin
    for source_ in [p1, sp, ds, dv, pm]
        for target in [p1, sp, ds, dv, pm]
            for source in [source_, source_', transpose(source_)]
                lres = source * target
                rres = target * source
                flres = Matrix(source) * Matrix(target)
                frres = Matrix(target) * Matrix(source)
                @test lres ≈ flres
                @test rres ≈ frres
                if !(target === p1 || parent(source) === p1)
                    @test eltype(lres) == eltype(flres)
                    @test eltype(rres) == eltype(frres)
                end
                if !(target |> isdense) && !(parent(source) |> isdense)
                    @test lres |> isdense == false
                    @test lres |> isdense == false
                end
            end
        end
    end
end

@testset "mul-vector" begin
    # permutation multiply
    lres = transpose(conj(v)) * pm  #! v' not realized!
    rres = pm * v
    flres = v' * Matrix(pm)
    frres = Matrix(pm) * v
    @test lres == flres
    @test rres == frres
    @test eltype(lres) == eltype(flres)
    @test eltype(rres) == eltype(frres)

    # IMatrix
    @test v' * p1 == v'
    @test p1 * v == v
end

@testset "sparse-diag" begin
    N = 100
    dg = Diagonal(randn(ComplexF64, 1000))
    sp = SparseMatrixCSC(dg)
    @test sp * dg == sp * sp
    @test dg * sp == sp * sp
end

@testset "randn" begin
    Random.seed!(2)
    T = ComplexF64
    for m in [sprand(T, 5, 5, 0.5), pmrand(T, 5), Diagonal(randn(T, 5))]
        zm = zero(m)
        @test zm ≈ zeros(T, 5, 5)
        rand!(zm)
        @test !(zm ≈ zeros(T, 5, 5))
        zm = zero(m)
        randn!(zm)
        @test !(zm ≈ zeros(T, 5, 5))
    end
end

@testset "multiply rectangular matrix" begin
    pm = pmrand(4)
    sp = sparse(reshape([1.0 2 3 4], 4, 1))
    res = pm * sp
    @test res ≈ Matrix(pm) * Matrix(sp)
end

@testset "findnz" begin
    for m in [p1, sp, ds, dv, pm]
        for _m in [m, staticize(m)]
            out = zeros(eltype(m), size(m)...)
            for (i, j, v) in zip(LuxurySparse.findnz(_m)...)
                out[i, j] = v
            end
            @test out ≈ m
        end
    end
end


text:
struct Lineage <: AbstractVector{AbstractTaxon}
    line::Vector{AbstractTaxon}
    index::Dict{Symbol,Int}
end

function Lineage(taxon::Taxon)
    line = Taxon[]
    current_taxon = taxon
    push!(line,current_taxon)
    while parent(current_taxon) !== nothing
        current_taxon = parent(current_taxon)
        push!(line, current_taxon)
    end
    reverse!(line)
    rankline = map(rank, line)
    index = Dict{Symbol,Int}()
    for crank in CanonicalRank
        position = findfirst(x -> x == crank, rankline)
        position === nothing ? continue : index[crank] = position
    end
    return Lineage(line,index)
end

Base.size(l::Lineage) = size(l.line)
Base.getindex(l::Lineage, i::Int) = getindex(l.line, i)
Base.lastindex(l::Lineage) = lastindex(l.line)

Base.getindex(l::Lineage, s::Symbol) = l.line[l.index[s]]

function Base.getindex(l::Lineage, range::UnitRange{Int})
    line = l.line[range]
    idx = Dict(key => value - range.start + 1 for (key, value) in l.index if value in range)
    return Lineage(line,idx)
end

function Base.getindex(l::Lineage, idx::All)
    if isempty(idx.cols)
        return l
    else
        return getindex(l, Cols(idx.cols))
    end
end

function Base.getindex(l::Lineage, idx::Cols)
    line = AbstractTaxon[]
    index = Dict{Symbol,Int}()
    count = 0
    for rank in idx.cols
        count += 1
        taxon = getindex(l, rank)
        push!(line, taxon)
        index[rank]=count
    end
    return Lineage(line,index)
end

Base.getindex(l::Lineage, idx::Between{Int,Int}) = l[idx.first:idx.last]
Base.getindex(l::Lineage, idx::Between{Symbol,Int}) = getindex(l, Between(l.index[idx.first], idx.last))
Base.getindex(l::Lineage, idx::Between{Int,Symbol}) = getindex(l, Between(idx.first, l.index[idx.last]))
Base.getindex(l::Lineage, idx::Between{Symbol,Symbol}) = getindex(l, Between(l.index[idx.first], l.index[idx.last]))
Base.getindex(l::Lineage, idx::From{Int}) = l[idx.first:end]
Base.getindex(l::Lineage, idx::From{Symbol}) = getindex(l, From(l.index[idx.first]))
Base.getindex(l::Lineage, idx::Until{Int}) = l[1:idx.last]
Base.getindex(l::Lineage, idx::Until{Symbol}) = getindex(l, Until(l.index[idx.last]))

function Base.get(l::Lineage, idx::Union{Int,Symbol}, default::Any)
    try
        return getindex(l,idx)
    catch
        return default
    end
end

"""
    reformat(l::Lineage, ranks::Vector{Symbol})

Return the `Lineage` object reformatted according to the given ranks.
"""

function reformat(l::Lineage, ranks::Vector{Symbol})
    line = AbstractTaxon[]
    idx = Dict{Symbol,Int}()
    count = 0
    for rank in ranks
        count += 1
        taxon = try
            getindex(l, rank)
        catch
            filtered_line = filter(x -> typeof(x) == Taxon, line)
            edge = try
                filtered_line[end]
            catch #if there is no taxon corresponding to the ranks
                l[end]
            end
            UnclassifiedTaxon(rank, edge)
        end
        push!(line, taxon)
        idx[rank]=count
    end
    return Lineage(line, idx)
end

"""
    print_lineage(lineage::Lineage; kwargs...)
    print_lineage(io::IO, lineage::Lineage; kwargs...)

Print a formatted representation of the lineage to the given `IO` object.

# Arguments

* `delim::AbstractString = ";"` - The delimiter between taxon fields.
* `fill::Bool = false` - If true, prints UnclassifiedTaxon. only availavle when skip is false
* `skip::Bool`= false` - If true, skip printing `UnclassifiedTaxon` and delimiter.
"""

function print_lineage(io::IO, lineage::Lineage; delim::AbstractString=";", fill::Bool=false, skip::Bool=false)
    name_line = String[]
    for taxon in lineage
        if typeof(taxon) == UnclassifiedTaxon
            if skip
                continue
            end

            if fill
                push!(name_line, name(taxon))
            else
                push!(name_line, "")
            end
        else
            push!(name_line, name(taxon))
        end
    end

    if isempty(name_line)
        return nothing
    end

    l = foldl((x,y) -> x*delim*y, name_line)
    print(io, l)
    return nothing
end

print_lineage(lineage::Lineage; kwargs...) = print_lineage(stdout::IO, lineage; kwargs...)

"""
    isdescendant(descendant::Taxon, ancestor::Taxon)

Return true if the former taxon is a descendant of the latter taxon.
"""

isdescendant(descendant::Taxon, ancestor::Taxon) = ancestor in Lineage(descendant)

"""
    isancestor(ancestor::Taxon, descendant::Taxon)

Return true if the former taxon is an ancestor of the latter taxon.
"""

isancestor(ancestor::Taxon, descendant::Taxon) = isdescendant(descendant, ancestor)


text:
using POMDPs
using ARDESPOT
using POMDPToolbox
using POMDPModels
using ProgressMeter

T = 50
N = 100

pomdp = BabyPOMDP()

bounds = IndependentBounds(reward(pomdp, true, false)/(1-discount(pomdp)), 0.0)

solver = DESPOTSolver(epsilon_0=0.1,
                      bounds=bounds,
                      T_max=Inf,
                      rng=MersenneTwister(4)
                     )

rsum = 0.0
fwc_rsum = 0.0
@showprogress for i in 1:N
    planner = solve(solver, pomdp)
    sim = RolloutSimulator(max_steps=T, rng=MersenneTwister(i))
    fwc_sim = deepcopy(sim)
    rsum += simulate(sim, pomdp, planner)
    fwc_rsum += simulate(sim, pomdp, FeedWhenCrying())
end

@show rsum/N
@show fwc_rsum/N


text:
function GRASP(
    cost, #Une array représentant les couts, de taille n
    matrix, #une matrice de taille m*n representant les contrantes
    n, #nombre de variables
    m, #nombbre de contraintes
    alpha # le alpha que nous allons utiliser
     )

     #intitialisation des listes utilisé
     (desactive_condition, stop1, variables_actives, stop2, util, SOL) = initaliser(m,n)

    #Création de la solution
    while desactive_condition!=stop1 && variables_actives!=stop2

        util = Utilite(cost, matrix, desactive_condition, n, m , variables_actives, stop1, stop2) #Foction d'utilité

        PosCandidat::Int64 = choixcandidat(util, alpha) #Choix parmis les candidats de util

        SOL[PosCandidat] = true #On ajoute le candidat à la solution

        Desactive!(PosCandidat, matrix, desactive_condition, m, variables_actives,n) #Desactive! le candidat selectionné et les conditions où il apparait
    end

    Z = calculz(SOL,cost,n)
    return(SOL, Z, desactive_condition)

end

# =========================================================================== #
function initaliser(m,n)
    #Creating a set of lines that will be avaluated
    desactive_condition= zeros(Bool, m) #si desactive_condition[j]=0 la ligne sera évalué
    stop1= ones(Bool,m) #la une condition d'arret, s'active quand toutes les lignes ont été traité

    #Creation the an array of activated variables
    variables_actives = ones(Bool, n) #Si variables_actives[i]=1, la variable sera traité
    stop2= zeros(Bool , n) #condition d'arret, s'active quand toutes les variables sont traité

    #On crée une liste contenant nos utilités
    util = zeros(Float64, n)

    #Initialisation de la Solution
    SOL = zeros(Bool, n)

    return(desactive_condition, stop1, variables_actives, stop2, util, SOL)
end

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

#Détermine l'utilité des éléments basé sur leur nombres d'apparition dans les condition et leur valeur dans les couts
function Utilite(cost, matrix, desactive_condition, n, m , variables_actives, stop1, stop2)

    util = zeros(Float64, n) #réinitialisation du vecteur

    for j=1:n #Pour chaque variable
            if variables_actives[j] == 1 #on verifie si la variable est active
                    for i=1:m #Pour chaque condition
                        if desactive_condition[i]==0 #On verifie si la condition doit être évalué
                            K=matrix[i,j]
                            util[j]=util[j]+K
                        end
                    end

                #On divise le le cout de chqaue variable par son nombre d'appartition dans les conditions
                if util[j]!=false #On évide de divisé par 0...
                    util[j] =  cost[j] / util[j]
                else #Si la variable n'a pas été décompté, on la désactive
                    variables_actives[j] = 0
                end
            end
    end
    return util
end
#============================================================#

function Desactive!(PosCandidat, matrix, desactive_condition, m , variables_actives , n) #On désactive les lignes où est le candidat
    for i=1:m
        if matrix[i,PosCandidat] == 1
            desactive_condition[i] = 1
            for j=1:n
                if matrix[i,j] ==1
                    variables_actives[j] = 0
                end
            end
        end
    end
end

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

function calculz(x,costs,m)
    z = 0
    for i in 1:m
        z+= x[i]*costs[i]
    end
    return z
end

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

function choixcandidat(util, alpha) #Permet de construire et de choisir un candidat dans la liste de candidat restreint
    max = maximum(util)
    min = minimum(util)
    Candidat_Restreint = Float64[]
    cpt::Int=1
    for i in 1:length(util)
        if util[i]>= min+(alpha*(max-min))
            push!(Candidat_Restreint,i)
        end
    end
    prob = rand(1:length(Candidat_Restreint))
    return(Candidat_Restreint[prob])
end

max_stars_repo_path: Archives/Triste vie 10-08/DM2_1.jl

max_issues_repo_path: Archives/Triste vie 10-08/DM2_1.jl

max_forks_repo_path: Archives/Triste vie 10-08/DM2_1.jl