text:
##==============================================================================
## AbstractDFG
##==============================================================================
##------------------------------------------------------------------------------
## Broadcasting
##------------------------------------------------------------------------------
# to allow stuff like `getFactorType.(dfg, [:x1x2f1;:x10l3f2])`
# https://docs.julialang.org/en/v1/manual/interfaces/#
Base.Broadcast.broadcastable(dfg::AbstractDFG) = Ref(dfg)

##==============================================================================
## Interface for an AbstractDFG
##==============================================================================
# Standard recommended fields to implement for AbstractDFG
# - `description::String`
# - `userId::String`
# - `robotId::String`
# - `sessionId::String`
# - `userData::Dict{Symbol, String}`
# - `robotData::Dict{Symbol, String}`
# - `sessionData::Dict{Symbol, String}`
# - `solverParams::T<:AbstractParams`
# - `addHistory::Vector{Symbol}`
# - `blobStores::Dict{Symbol, AbstractBlobStore}`
# AbstractDFG Accessors

##------------------------------------------------------------------------------
## Getters
##------------------------------------------------------------------------------
"""
    $(SIGNATURES)
Convenience function to get all the metadata of a DFG
"""
getDFGInfo(dfg::AbstractDFG) = (getDescription(dfg), getUserId(dfg), getRobotId(dfg), getSessionId(dfg), getUserData(dfg), getRobotData(dfg), getSessionData(dfg), getSolverParams(dfg))

"""
    $(SIGNATURES)
"""
getDescription(dfg::AbstractDFG) = dfg.description

"""
    $(SIGNATURES)
"""
getUserId(dfg::AbstractDFG) = dfg.userId

"""
    $(SIGNATURES)
"""
getRobotId(dfg::AbstractDFG) = dfg.robotId

"""
    $(SIGNATURES)
"""
getSessionId(dfg::AbstractDFG) = dfg.sessionId

"""
    $(SIGNATURES)
"""
getAddHistory(dfg::AbstractDFG) = dfg.addHistory

"""
    $(SIGNATURES)
"""
getSolverParams(dfg::AbstractDFG) = dfg.solverParams

"""
    $(SIGNATURES)

Method must be overloaded by the user for Serialization to work.  E.g. IncrementalInference uses `CommonConvWrapper <: FactorOperationalMemory`.
"""
getFactorOperationalMemoryType(dummy) = error("Please extend your workspace with function getFactorOperationalMemoryType(<:AbstractParams) for your usecase, e.g. IncrementalInference uses `CommonConvWrapper <: FactorOperationalMemory`")
getFactorOperationalMemoryType(dfg::AbstractDFG) = getFactorOperationalMemoryType(getSolverParams(dfg))

"""
    $(SIGNATURES)

Method must be overloaded by the user for Serialization to work.
"""
rebuildFactorMetadata!(dfg::AbstractDFG{<:AbstractParams}, factor::AbstractDFGFactor, neighbors=[]) = error("rebuildFactorMetadata! is not implemented for $(typeof(dfg))")

##------------------------------------------------------------------------------
## Setters
##------------------------------------------------------------------------------

"""
    $(SIGNATURES)
"""
setDescription!(dfg::AbstractDFG, description::String) = dfg.description = description

"""
    $(SIGNATURES)
"""
setUserId!(dfg::AbstractDFG, userId::String) = dfg.userId = userId

"""
    $(SIGNATURES)
"""
setRobotId!(dfg::AbstractDFG, robotId::String) = dfg.robotId = robotId

"""
    $(SIGNATURES)
"""
setSessionId!(dfg::AbstractDFG, sessionId::String) = dfg.sessionId = sessionId

"""
    $(SIGNATURES)
"""
#NOTE a MethodError will be thrown if solverParams type does not mach the one in dfg
# TODO Is it ok or do we want any abstract solver paramters
setSolverParams!(dfg::AbstractDFG, solverParams::AbstractParams) = dfg.solverParams = solverParams

# Accessors and CRUD for user/robot/session Data
"""
$SIGNATURES

Get the user data associated with the graph.
"""
getUserData(dfg::AbstractDFG) = return dfg.userData

"""
$SIGNATURES

Set the user data associated with the graph.
"""
function setUserData!(dfg::AbstractDFG, data::Dict{Symbol, String})
    dfg.userData = data #TODO keep memory? use clear and then add
    return dfg.userData
end

"""
$SIGNATURES

Get the robot data associated with the graph.
"""
getRobotData(dfg::AbstractDFG)::Union{Nothing, Dict{Symbol, String}} = return dfg.robotData

"""
$SIGNATURES

Set the robot data associated with the graph.
"""
function setRobotData!(dfg::AbstractDFG, data::Dict{Symbol, String})
    dfg.robotData = data
    return dfg.robotData
end

"""
$SIGNATURES

Get the session data associated with the graph.
"""
getSessionData(dfg::AbstractDFG)::Dict{Symbol, String} = return dfg.sessionData

"""
$SIGNATURES

Set the session data associated with the graph.
"""
function setSessionData!(dfg::AbstractDFG, data::Dict{Symbol, String})
    dfg.sessionData = data
    return dfg.sessionData
end

##==============================================================================
## User/Robot/Session Data CRUD
##==============================================================================

#NOTE with API standardization this should become something like:
getUserData(dfg::AbstractDFG, key::Symbol)::String = dfg.userData[key]
getRobotData(dfg::AbstractDFG, key::Symbol)::String = dfg.robotData[key]
getSessionData(dfg::AbstractDFG, key::Symbol)::String = dfg.sessionData[key]

updateUserData!(dfg::AbstractDFG, pair::Pair{Symbol,String}) = push!(dfg.userData, pair)
updateRobotData!(dfg::AbstractDFG, pair::Pair{Symbol,String}) = push!(dfg.robotData, pair)
updateSessionData!(dfg::AbstractDFG, pair::Pair{Symbol,String}) = push!(dfg.sessionData, pair)

deleteUserData!(dfg::AbstractDFG, key::Symbol) = pop!(dfg.userData, key)
deleteRobotData!(dfg::AbstractDFG, key::Symbol) = pop!(dfg.robotData, key)
deleteSessionData!(dfg::AbstractDFG, key::Symbol) = pop!(dfg.sessionData, key)

emptyUserData!(dfg::AbstractDFG) = empty!(dfg.userData)
emptyRobotData!(dfg::AbstractDFG) = empty!(dfg.robotData)
emptySessionData!(dfg::AbstractDFG) = empty!(dfg.sessionData)

#TODO add__Data!?

##==============================================================================
## AbstractBlobStore  CRUD
##==============================================================================
abstract type AbstractBlobStore{T} end
# AbstractBlobStore should have key or overwrite getKey
getKey(store::AbstractBlobStore) = store.key

getBlobStore(dfg::AbstractDFG, key::Symbol) = dfg.blobStores[key]
addBlobStore!(dfg::AbstractDFG, bs::AbstractBlobStore) = push!(dfg.blobStores, getKey(bs)=>bs)
updateBlobStore!(dfg::AbstractDFG, bs::AbstractBlobStore) = push!(dfg.blobStores, getKey(bs)=>bs)
deleteBlobStore!(dfg::AbstractDFG, key::Symbol) = pop!(dfg.blobStores, key)
emptyBlobStore!(dfg::AbstractDFG) = empty!(dfg.blobStores)
listBlobStores(dfg::AbstractDFG) = collect(keys(dfg.blobStores))

##==============================================================================
## CRUD Interfaces
##==============================================================================
##------------------------------------------------------------------------------
## Variable And Factor CRUD
##------------------------------------------------------------------------------
"""
    $(SIGNATURES)
True if the variable or factor exists in the graph.
"""
function exists(dfg::AbstractDFG, node::DFGNode)
    error("exists not implemented for $(typeof(dfg))")
end

function exists(dfg::AbstractDFG, label::Symbol)
    error("exists not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Add a DFGVariable to a DFG.
"""
function addVariable!(dfg::G, variable::V) where {G <: AbstractDFG, V <: AbstractDFGVariable}
    error("addVariable! not implemented for $(typeof(dfg))")
end

"""
Add a DFGFactor to a DFG.
    $(SIGNATURES)
"""
function addFactor!(dfg::AbstractDFG, factor::F) where F <: AbstractDFGFactor
    error("addFactor! not implemented for $(typeof(dfg))(dfg, factor)")
end

"""
    $(SIGNATURES)
Get a DFGVariable from a DFG using its label.
"""
function getVariable(dfg::G, label::Union{Symbol, String}) where G <: AbstractDFG
    error("getVariable not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Get a DFGFactor from a DFG using its label.
"""
function getFactor(dfg::G, label::Union{Symbol, String}) where G <: AbstractDFG
    error("getFactor not implemented for $(typeof(dfg))")
end

function Base.getindex(dfg::AbstractDFG, lbl::Union{Symbol, String})
    if isVariable(dfg, lbl)
        getVariable(dfg, lbl)
    elseif isFactor(dfg, lbl)
        getFactor(dfg, lbl)
    else
        error("Cannot find $lbl in this $(typeof(dfg))")
    end
end

"""
    $(SIGNATURES)
Update a complete DFGVariable in the DFG.
"""
function updateVariable!(dfg::G, variable::V) where {G <: AbstractDFG, V <: AbstractDFGVariable}
    error("updateVariable! not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Update a complete DFGFactor in the DFG.
"""
function updateFactor!(dfg::G, factor::F) where {G <: AbstractDFG, F <: AbstractDFGFactor}
    error("updateFactor! not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Delete a DFGVariable from the DFG using its label.
"""
function deleteVariable!(dfg::AbstractDFG, label::Symbol)
    error("deleteVariable! not implemented for $(typeof(dfg))")
end
"""
    $(SIGNATURES)
Delete a DFGFactor from the DFG using its label.
"""
function deleteFactor!(dfg::G, label::Symbol; suppressGetFactor::Bool=false) where G <: AbstractDFG
    error("deleteFactors not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
List the DFGVariables in the DFG.
Optionally specify a label regular expression to retrieves a subset of the variables.
Tags is a list of any tags that a node must have (at least one match).
"""
function getVariables(dfg::G, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[], solvable::Int=0) where G <: AbstractDFG
    error("getVariables not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
List the DFGFactors in the DFG.
Optionally specify a label regular expression to retrieves a subset of the factors.
"""
function getFactors(dfg::G, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[], solvable::Int=0) where G <: AbstractDFG
    error("getFactors not implemented for $(typeof(dfg))")
end

##------------------------------------------------------------------------------
## Checking Types
##------------------------------------------------------------------------------

"""
    $SIGNATURES

Return whether `sym::Symbol` represents a variable vertex in the graph DFG.
Checks whether it both exists in the graph and is a variable.
(If you rather want a quick for type, just do node isa DFGVariable)
"""
function isVariable(dfg::G, sym::Symbol) where G <: AbstractDFG
    error("isVariable not implemented for $(typeof(dfg))")
end

"""
    $SIGNATURES

Return whether `sym::Symbol` represents a factor vertex in the graph DFG.
Checks whether it both exists in the graph and is a factor.
(If you rather want a quicker for type, just do node isa DFGFactor)
"""
function isFactor(dfg::G, sym::Symbol) where G <: AbstractDFG
    error("isFactor not implemented for $(typeof(dfg))")
end

##------------------------------------------------------------------------------
## Neighbors
##------------------------------------------------------------------------------
"""
    $(SIGNATURES)
Checks if the graph is fully connected, returns true if so.
"""
function isConnected(dfg::AbstractDFG)
    error("isConnected not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor specified by its label.
"""
function getNeighbors(dfg::AbstractDFG, label::Symbol; solvable::Int=0)
    error("getNeighbors not implemented for $(typeof(dfg))")
end

##------------------------------------------------------------------------------
## copy and duplication
##------------------------------------------------------------------------------

#TODO use copy functions currently in attic
"""
    $(SIGNATURES)
Gets an empty and unique CloudGraphsDFG derived from an existing DFG.
"""
function _getDuplicatedEmptyDFG(dfg::G) where G <: AbstractDFG
    error("_getDuplicatedEmptyDFG not implemented for $(typeof(dfg))")
end

##------------------------------------------------------------------------------
## CRUD Aliases
##------------------------------------------------------------------------------

"""
    $(SIGNATURES)
Get a DFGVariable with a specific solver key.
In memory types still return a reference, other types returns a variable with only solveKey.
"""
function getVariable(dfg::AbstractDFG, label::Symbol, solveKey::Symbol)

    var = getVariable(dfg, label)

    if isa(var, DFGVariable) && !haskey(var.solverDataDict, solveKey)
        error("Solvekey '$solveKey' does not exists in the variable")
    elseif !isa(var, DFGVariable)
        @warn "getVariable(dfg, label, solveKey) only supported for type DFGVariable."
    end

    return var
end

"""
$(SIGNATURES)
"""
function addFactor!(dfg::AbstractDFG, variables::Vector{<:AbstractDFGVariable}, factor::F) where  F <: AbstractDFGFactor

    Base.depwarn("addFactor!(dfg, variables, factor) is deprecated, use addFactor!(dfg, factor)", :addFactor!)
    variableLabels = map(v->v.label, variables)

    if factor isa DFGFactor
        f = factor
        newfactor =  DFGFactor(f.label, f.timestamp, f.nstime, f.tags, f.solverData, f.solvable, Tuple(variableLabels))
        return addFactor!(dfg, newfactor)
    else
        resize!(factor._variableOrderSymbols, length(variableLabels))
        factor._variableOrderSymbols .= variableLabels
        return addFactor!(dfg, factor)
    end

end

"""
$(SIGNATURES)
"""
function addFactor!(dfg::AbstractDFG, variableLabels::Vector{Symbol}, factor::F) where F <: AbstractDFGFactor

    Base.depwarn("addFactor!(dfg, variables, factor) is deprecated, use addFactor!(dfg, factor)", :addFactor!)

    if factor isa DFGFactor
        f = factor
        newfactor =  DFGFactor(f.label, f.timestamp, f.nstime, f.tags, f.solverData, f.solvable, Tuple(variableLabels))
        return addFactor!(dfg, newfactor)
    else
        resize!(factor._variableOrderSymbols, length(variableLabels))
        factor._variableOrderSymbols .= variableLabels
        return addFactor!(dfg, factor)
    end

end

"""
    $(SIGNATURES)
Delete a referenced DFGVariable from the DFG.

Notes
- Returns `Tuple{AbstractDFGVariable, Vector{<:AbstractDFGFactor}}`
"""
function deleteVariable!(dfg::AbstractDFG, variable::AbstractDFGVariable)
    return deleteVariable!(dfg, variable.label)
end

"""
    $(SIGNATURES)
Delete the referened DFGFactor from the DFG.
"""
function deleteFactor!(dfg::G, factor::F; suppressGetFactor::Bool=false) where {G <: AbstractDFG, F <: AbstractDFGFactor}
    return deleteFactor!(dfg, factor.label, suppressGetFactor=suppressGetFactor)
end

# Alias - bit ridiculous but know it'll come up at some point. Does existential and type check.
function isVariable(dfg::G, node::N) where {G <: AbstractDFG, N <: DFGNode}
    return isVariable(dfg, node.label)
end
# Alias - bit ridiculous but know it'll come up at some point. Does existential and type check.
function isFactor(dfg::G, node::N) where {G <: AbstractDFG, N <: DFGNode}
    return isFactor(dfg, node.label)
end

##------------------------------------------------------------------------------
## Connectivity Alias
##------------------------------------------------------------------------------

function getNeighbors(dfg::AbstractDFG, node::DFGNode; solvable::Int=0)
    getNeighbors(dfg, node.label, solvable=solvable)
end

##==============================================================================
## Listing and listing aliases
##==============================================================================

##------------------------------------------------------------------------------
## Overwrite in driver for performance
##------------------------------------------------------------------------------
"""
    $(SIGNATURES)
Get a list of IDs of the DFGVariables in the DFG.
Optionally specify a label regular expression to retrieves a subset of the variables.
Tags is a list of any tags that a node must have (at least one match).

NOtes
- Returns `::Vector{Symbol}`

Example
```julia
listVariables(dfg, r"l", tags=[:APRILTAG;])
```

Related:
- ls
"""
function listVariables( dfg::AbstractDFG,
                        regexFilter::Union{Nothing, Regex}=nothing;
                        tags::Vector{Symbol}=Symbol[],
                        solvable::Int=0 )
    #
    vars = getVariables(dfg, regexFilter, tags=tags, solvable=solvable)
    return map(v -> v.label, vars)::Vector{Symbol}
end

# to be consolidated, see #612
function listVariables( dfg::AbstractDFG,
                        typeFilter::Type{<:InferenceVariable};
                        tags::Vector{Symbol}=Symbol[],
                        solvable::Int=0 )
    #
    retlist::Vector{Symbol} = ls(dfg, typeFilter)
    0 < length(tags) || solvable != 0 ? intersect(retlist, ls(dfg, tags=tags, solvable=solvable)) : retlist
end

"""
    $(SIGNATURES)
Get a list of the IDs (labels) of the DFGFactors in the DFG.
Optionally specify a label regular expression to retrieves a subset of the factors.
"""
function listFactors(dfg::G, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[], solvable::Int=0)::Vector{Symbol} where G <: AbstractDFG
    return map(f -> f.label, getFactors(dfg, regexFilter, tags=tags, solvable=solvable))
end

"""
    $TYPEDSIGNATURES
List all the solvekeys used amongst all variables in the distributed factor graph object.

Related

[`listSupersolves`](@ref), [`getSolverDataDict`](@ref), [`listVariables`](@ref)
"""
function listSolveKeys( variable::DFGVariable,
                        filterSolveKeys::Union{Regex,Nothing}=nothing,
                        skeys = Set{Symbol}() )
    #
    for ky in keys(getSolverDataDict(variable))
        push!(skeys, ky)
    end

    #filter the solveKey set with filterSolveKeys regex
    !isnothing(filterSolveKeys) && return filter!(k -> occursin(filterSolveKeys, string(k)), skeys)
    return skeys
end

listSolveKeys(  dfg::AbstractDFG, lbl::Symbol,
                filterSolveKeys::Union{Regex,Nothing}=nothing,
                skeys = Set{Symbol}() ) = listSolveKeys(getVariable(dfg, lbl), filterSolveKeys, skeys)
#

function listSolveKeys( dfg::AbstractDFG,
                        filterVariables::Union{Type{<:InferenceVariable},Regex, Nothing}=nothing;
                        filterSolveKeys::Union{Regex,Nothing}=nothing,
                        tags::Vector{Symbol}=Symbol[],
                        solvable::Int=0  )
    #
    skeys = Set{Symbol}()
    varList = listVariables(dfg, filterVariables, tags=tags, solvable=solvable)
    for vs in varList  #, ky in keys(getSolverDataDict(getVariable(dfg, vs)))
        listSolveKeys(dfg, vs, filterSolveKeys, skeys)
    end

    # done inside the loop
    # #filter the solveKey set with filterSolveKeys regex
    # !isnothing(filterSolveKeys) && return filter!(k -> occursin(filterSolveKeys, string(k)), skeys)

    return skeys
end
const listSupersolves = listSolveKeys

##------------------------------------------------------------------------------
## Aliases and Other filtered lists
##------------------------------------------------------------------------------

## Aliases
##--------
"""
    $(SIGNATURES)
List the DFGVariables in the DFG.
Optionally specify a label regular expression to retrieves a subset of the variables.
Tags is a list of any tags that a node must have (at least one match).

Notes:
- Returns `Vector{Symbol}`
"""
function ls(dfg::G, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[], solvable::Int=0) where G <: AbstractDFG
    return listVariables(dfg, regexFilter, tags=tags, solvable=solvable)
end

#TODO tags kwarg
"""
    $(SIGNATURES)
List the DFGFactors in the DFG.
Optionally specify a label regular expression to retrieves a subset of the factors.

Notes
- Return `Vector{Symbol}`
"""
function lsf(dfg::G, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[], solvable::Int=0) where G <: AbstractDFG
    return listFactors(dfg, regexFilter, tags=tags, solvable=solvable)
end

"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor.
"""
function ls(dfg::G, node::T; solvable::Int=0) where {G <: AbstractDFG, T <: DFGNode}
    return getNeighbors(dfg, node, solvable=solvable)
end
function ls(dfg::G, label::Symbol; solvable::Int=0) where G <: AbstractDFG
    return getNeighbors(dfg, label, solvable=solvable)
end

function lsf(dfg::G, label::Symbol; solvable::Int=0)::Vector{Symbol} where G <: AbstractDFG
  return getNeighbors(dfg, label, solvable=solvable)
end

## list by types
##--------------

function ls(dfg::G, ::Type{T}) where {G <: AbstractDFG, T <: InferenceVariable}
  xx = getVariables(dfg)
  mask = getVariableType.(xx) .|> typeof .== T
  vxx = view(xx, mask)
  map(x->x.label, vxx)
end

function ls(dfg::G, ::Type{T}) where {G <: AbstractDFG, T <: FunctorInferenceType}
  xx = getFactors(dfg)
  names = typeof.(getFactorType.(xx)) .|> nameof
  vxx = view(xx, names .== Symbol(T))
  map(x->x.label, vxx)
end

function lsf(dfg::G, ::Type{T}) where {G <: AbstractDFG, T <: FunctorInferenceType}
    ls(dfg, T)
end

"""
    $(SIGNATURES)
Helper to return neighbors at distance 2 around a given node label.
"""
function ls2(dfg::AbstractDFG, label::Symbol)
    l2 = getNeighborhood(dfg, label, 2)
    l1 = getNeighborhood(dfg, label, 1)
    return setdiff(l2, l1)
end

"""
    $SIGNATURES

Return vector of prior factor symbol labels in factor graph `dfg`.

Notes:
- Returns `Vector{Symbol}`
"""
function lsfPriors(dfg::G) where G <: AbstractDFG
  priors = Symbol[]
  fcts = lsf(dfg)
  for fc in fcts
    if isPrior(dfg, fc)
      push!(priors, fc)
    end
  end
  return priors
end

#TODO is this repeated functionality?

"""
    $(SIGNATURES)
Gives back all factor labels that fit the bill:
    lsWho(dfg, :Pose3)

Notes
- Returns `Vector{Symbol}`

Dev Notes
- Cloud versions will benefit from less data transfer
 - `ls(dfg::C, ::T) where {C <: CloudDFG, T <: ..}`

Related

ls, lsf, lsfPriors
"""
function lsWho(dfg::AbstractDFG, type::Symbol)
    vars = getVariables(dfg)
    labels = Symbol[]
    for v in vars
        varType = typeof(getVariableType(v)) |> nameof
        varType == type && push!(labels, v.label)
    end
    return labels
end

"""
    $(SIGNATURES)
Gives back all factor labels that fit the bill:
    lsfWho(dfg, :Point2Point2)

Notes
- Returns `Vector{Symbol}`

Dev Notes
- Cloud versions will benefit from less data transfer
 - `ls(dfg::C, ::T) where {C <: CloudDFG, T <: ..}`

Related

ls, lsf, lsfPriors
"""
function lsfWho(dfg::AbstractDFG, type::Symbol)
    facs = getFactors(dfg)
    labels = Symbol[]
    for f in facs
        facType = typeof(getFactorType(f)) |> nameof
        facType == type && push!(labels, f.label)
    end
    return labels
end

## list types
##-----------

"""
    $SIGNATURES

Return `Vector{Symbol}` of all unique variable types in factor graph.
"""
function lsTypes(dfg::AbstractDFG)
    vars = getVariables(dfg)
    alltypes = Set{Symbol}()
    for v in vars
        varType = typeof(getVariableType(v)) |> nameof
        push!(alltypes, varType)
    end
    return collect(alltypes)
end

"""
    $SIGNATURES

Return `::Dict{Symbol, Vector{Symbol}}` of all unique variable types with labels in a factor graph.
"""
function lsTypesDict(dfg::AbstractDFG)

    vars = getVariables(dfg)
    alltypes = Dict{Symbol,Vector{Symbol}}()
    for v in vars
        varType = typeof(getVariableType(v)) |> nameof
        d = get!(alltypes, varType, Symbol[])
        push!(d, v.label)
    end
  return alltypes
end

"""
    $SIGNATURES

Return `Vector{Symbol}` of all unique factor types in factor graph.
"""
function lsfTypes(dfg::AbstractDFG)
    facs = getFactors(dfg)
    alltypes = Set{Symbol}()
    for f in facs
        facType = typeof(getFactorType(f)) |> nameof
        push!(alltypes, facType)
    end
    return collect(alltypes)
end

"""
    $SIGNATURES

Return `::Dict{Symbol, Vector{Symbol}}` of all unique factors types with labels in a factor graph.
"""
function lsfTypesDict(dfg::AbstractDFG)
    facs = getFactors(dfg)
    alltypes = Dict{Symbol,Vector{Symbol}}()
    for f in facs
        facType = typeof(getFactorType(f)) |> nameof
        d = get!(alltypes, facType, Symbol[])
        push!(d, f.label)
    end
  return alltypes
end

##------------------------------------------------------------------------------
## tags
##------------------------------------------------------------------------------
"""
    $SIGNATURES

Determine if the variable or factor neighbors have the `tags:;Vector{Symbol}`, and `matchAll::Bool`.
"""
function hasTags(dfg::AbstractDFG,
                 sym::Symbol,
                 tags::Vector{Symbol};
                 matchAll::Bool=true  )
  #
  alltags = listTags(dfg, sym)
  length(filter(x->x in alltags, tags)) >= (matchAll ? length(tags) : 1)
end

"""
    $SIGNATURES

Determine if the variable or factor neighbors have the `tags:;Vector{Symbol}`, and `matchAll::Bool`.
"""
function hasTagsNeighbors(dfg::AbstractDFG,
                          sym::Symbol,
                          tags::Vector{Symbol};
                          matchAll::Bool=true  )
  #
  # assume only variables or factors are neighbors
  getNeiFnc = isVariable(dfg, sym) ? getFactor : getVariable
  alltags = union( (ls(dfg, sym) .|> x->getTags(getNeiFnc(dfg,x)))... )
  length(filter(x->x in alltags, tags)) >= (matchAll ? length(tags) : 1)
end

##==============================================================================
## Finding
##==============================================================================
# function findClosestTimestamp(setA::Vector{Tuple{DateTime,T}},
                              # setB::Vector{Tuple{DateTime,S}}) where {S,T}
"""
    $SIGNATURES

Find and return the closest timestamp from two sets of Tuples.  Also return the minimum delta-time (`::Millisecond`) and how many elements match from the two sets are separated by the minimum delta-time.
"""
function findClosestTimestamp(setA::Vector{Tuple{ZonedDateTime,T}},
                              setB::Vector{Tuple{ZonedDateTime,S}}) where {S,T}
  #
  # build matrix of delta times, ranges on rows x vars on columns
  DT = Array{Millisecond, 2}(undef, length(setA), length(setB))
  for i in 1:length(setA), j in 1:length(setB)
    DT[i,j] = setB[j][1] - setA[i][1]
  end

  DT .= abs.(DT)

  # absolute time differences
  # DTi = (x->x.value).(DT) .|> abs

  # find the smallest element
  mdt = minimum(DT)
  corrs = findall(x->x==mdt, DT)

  # return the closest timestamp, deltaT, number of correspondences
  return corrs[1].I, mdt, length(corrs)
end

"""
    $SIGNATURES

Find and return nearest variable labels per delta time.  Function will filter on `regexFilter`, `tags`, and `solvable`.

Notes
- Returns `Vector{Tuple{Vector{Symbol}, Millisecond}}`

DevNotes:
- TODO `number` should allow returning more than one for k-nearest matches.
- Future versions likely will require some optimization around the internal `getVariable` call.
  - Perhaps a dedicated/efficient `getVariableTimestamp` for all DFG flavors.

Related

ls, listVariables, findClosestTimestamp
"""
function findVariableNearTimestamp(dfg::AbstractDFG,
                                   timest::ZonedDateTime,
                                   regexFilter::Union{Nothing, Regex}=nothing;
                                   tags::Vector{Symbol}=Symbol[],
                                   solvable::Int=0,
                                   warnDuplicate::Bool=true,
                                   number::Int=1  )
  #
  # get the variable labels based on filters
  # syms = listVariables(dfg, regexFilter, tags=tags, solvable=solvable)
  syms = listVariables(dfg, regexFilter, tags=tags, solvable=solvable)
  # compile timestamps with label
  # vars = map( x->getVariable(dfg, x), syms )
  timeset = map(x->(getTimestamp(getVariable(dfg,x)), x), syms)
  mask = BitArray{1}(undef, length(syms))
  fill!(mask, true)

  RET = Vector{Tuple{Vector{Symbol},Millisecond}}()
  SYMS = Symbol[]
  CORRS = 1
  NUMBER = number
  while 0 < CORRS + NUMBER
    # get closest
    link, mdt, corrs = findClosestTimestamp([(timest,0)], timeset[mask])
    newsym = syms[link[2]]
    union!(SYMS, !isa(newsym, Vector) ? [newsym] : newsym)
    mask[link[2]] = false
    CORRS = corrs-1
    # last match, done with this delta time
    if corrs == 1
      NUMBER -=  1
      push!(RET, (deepcopy(SYMS),mdt))
      SYMS = Symbol[]
    end
  end
  # warn if duplicates found
  # warnDuplicate && 1 < corrs ? @warn("getVariableNearTimestamp found more than one variable at $timestamp") :   nothing

  return RET
end

findVariableNearTimestamp(dfg::AbstractDFG, timest::DateTime, regexFilter::Union{Nothing, Regex}=nothing;
                          timezone=localzone(), kwargs...) =
    findVariableNearTimestamp(dfg, ZonedDateTime(timest, timezone), regexFilter; kwargs...)

##==============================================================================
## Copy Functions
##==============================================================================

"""
    $(SIGNATURES)
Common function for copying nodes from one graph into another graph.
This is overridden in specialized implementations for performance.
Orphaned factors are not added, with a warning if verbose.
Set `overwriteDest` to overwrite existing variables and factors in the destination DFG.
NOTE: copyGraphMetadata not supported yet.
Related:
- [`deepcopyGraph`](@ref)
- [`deepcopyGraph!`](@ref)
- [`buildSubgraph`](@ref)
- [`getNeighborhood`](@ref)
- [`mergeGraph!`](@ref)
"""
function copyGraph!(destDFG::AbstractDFG,
                    sourceDFG::AbstractDFG,
                    variableLabels::Vector{Symbol},
                    factorLabels::Vector{Symbol};
                    copyGraphMetadata::Bool=false,
                    overwriteDest::Bool=false,
                    deepcopyNodes::Bool=false,
                    verbose::Bool = true)
    # Split into variables and factors
    sourceVariables = map(vId->getVariable(sourceDFG, vId), variableLabels)
    sourceFactors = map(fId->getFactor(sourceDFG, fId), factorLabels)

    # Now we have to add all variables first,
    for variable in sourceVariables
        variableCopy = deepcopyNodes ? deepcopy(variable) : variable
        if !exists(destDFG, variable)
            addVariable!(destDFG, variableCopy)
        elseif overwriteDest
            updateVariable!(destDFG, variableCopy)
        else
            error("Variable $(variable.label) already exists in destination graph!")
        end
    end
    # And then all factors to the destDFG.
    for factor in sourceFactors
        # Get the original factor variables (we need them to create it)
        sourceFactorVariableIds = getNeighbors(sourceDFG, factor)
        # Find the labels and associated variables in our new subgraph
        factVariableIds = Symbol[]
        for variable in sourceFactorVariableIds
            if exists(destDFG, variable)
                push!(factVariableIds, variable)
            end
        end
        # Only if we have all of them should we add it (otherwise strange things may happen on evaluation)
        if length(factVariableIds) == length(sourceFactorVariableIds)
            factorCopy = deepcopyNodes ? deepcopy(factor) : factor
            if !exists(destDFG, factor)
                addFactor!(destDFG, factorCopy)
            elseif overwriteDest
                updateFactor!(destDFG, factorCopy)
            else
                error("Factor $(factor.label) already exists in destination graph!")
            end
        elseif verbose
            @warn "Factor $(factor.label) will be an orphan in the destination graph, and therefore not added."
        end
    end

    if copyGraphMetadata
        setUserData(destDFG, getUserData(sourceDFG))
        setRobotData(destDFG, getRobotData(sourceDFG))
        setSessionData(destDFG, getSessionData(sourceDFG))
    end
    return nothing
end

"""
    $(SIGNATURES)
Copy nodes from one graph into another graph by making deepcopies.
see [`copyGraph!`](@ref) for more detail.
Related:
- [`deepcopyGraph`](@ref)
- [`buildSubgraph`](@ref)
- [`getNeighborhood`](@ref)
- [`mergeGraph!`](@ref)
"""
function deepcopyGraph!(destDFG::AbstractDFG,
                        sourceDFG::AbstractDFG,
                        variableLabels::Vector{Symbol} = ls(sourceDFG),
                        factorLabels::Vector{Symbol} = lsf(sourceDFG);
                        kwargs...)
    copyGraph!(destDFG, sourceDFG, variableLabels, factorLabels; deepcopyNodes=true, kwargs...)
end

"""
    $(SIGNATURES)
Copy nodes from one graph into a new graph by making deepcopies.
see [`copyGraph!`](@ref) for more detail.
Related:
- [`deepcopyGraph!`](@ref)
- [`buildSubgraph`](@ref)
- [`getNeighborhood`](@ref)
- [`mergeGraph!`](@ref)
"""
function deepcopyGraph( ::Type{T},
                        sourceDFG::AbstractDFG,
                        variableLabels::Vector{Symbol} = ls(sourceDFG),
                        factorLabels::Vector{Symbol} = lsf(sourceDFG);
                        sessionId::String = "",
                        kwargs...) where T <: AbstractDFG

    ginfo = [getDFGInfo(sourceDFG)...]
    if sessionId == ""
        ginfo[4] *= "_copy$(string(uuid4())[1:6])"
    else
        ginfo[4] = sessionId
    end
    destDFG = T(ginfo...)
    copyGraph!(destDFG, sourceDFG, variableLabels, factorLabels; deepcopyNodes=true, kwargs...)
    return destDFG
end

##==============================================================================
## Automated Graph Searching
##==============================================================================

function findShortestPathDijkstra end

"""
    $SIGNATURES

Relatively naive function counting linearly from-to

DevNotes
- Convert to using LightGraphs shortest path methods instead.
"""
function findFactorsBetweenNaive(   dfg::AbstractDFG,
                                    from::Symbol,
                                    to::Symbol,
                                    assertSingles::Bool=false)
  #
  @info "findFactorsBetweenNaive is naive linear number method -- improvements welcome"
  SRT = getVariableLabelNumber(from)
  STP = getVariableLabelNumber(to)
  prefix = string(from)[1]
  @assert prefix == string(to)[1] "from-to prefixes must match, one is $prefix, other $(string(to)[1])"
  prev = from
  fctlist = Symbol[]
  for num in (SRT+1):STP
    next = Symbol(prefix,num)
    fct = intersect(ls(dfg, prev),ls(dfg,next))
    !assertSingles ? nothing : @assert length(fct) == 1 "assertSingles=true, won't return multiple factors joining variables at this time"
    union!(fctlist, fct)
    prev = next
  end

  return fctlist
end

"""
    $SIGNATURES
Return (::Bool,::Vector{TypeName}) of types between two nodes in the factor graph

DevNotes
- Only works on LigthDFG at the moment.

Related

[`LightDFG.findShortestPathDijkstra`](@ref)
"""
function isPathFactorsHomogeneous(dfg::AbstractDFG, from::Symbol, to::Symbol)
    # FIXME, must consider all paths, not just shortest...
    pth = intersect(findShortestPathDijkstra(dfg, from, to), lsf(dfg))
    types = getFactorType.(dfg, pth) .|> typeof .|> x->(x).name #TODO this might not be correct in julia 1.6
    utyp = unique(types)
    (length(utyp) == 1), utyp
end

function existsPathOfFactorsType(dfg::AbstractDFG, from::Symbol, to::Symbol, ftype::FunctorInferenceType)
  error("WIP")
end

##==============================================================================
## Subgraphs and Neighborhoods
##==============================================================================

"""
    $(SIGNATURES)
Build a list of all unique neighbors inside 'distance'

Notes
- Returns `Vector{Symbol}`

Related:
- [`copyGraph!`](@ref)
- [`buildSubgraph`](@ref)
- [`deepcopyGraph`](@ref)
- [`mergeGraph!`](@ref)
"""
function getNeighborhood(dfg::AbstractDFG, label::Symbol, distance::Int)
    neighborList = Set{Symbol}([label])
    curList = Set{Symbol}([label])

    for dist in 1:distance
        newNeighbors = Set{Symbol}()
        for node in curList
            neighbors = getNeighbors(dfg, node)
            for neighbor in neighbors
                push!(neighborList, neighbor)
                push!(newNeighbors, neighbor)
            end
        end
        curList = newNeighbors
    end
    return collect(neighborList)
end

function getNeighborhood(dfg::AbstractDFG, variableFactorLabels::Vector{Symbol}, distance::Int; solvable::Int=0)
    # find neighbors at distance to add
    neighbors = Symbol[]
    if distance > 0
        for l in variableFactorLabels
            union!(neighbors, getNeighborhood(dfg, l, distance))
        end
    end

    allvarfacs = union(variableFactorLabels, neighbors)

    solvable != 0 && filter!(nlbl -> (getSolvable(dfg, nlbl) >= solvable), allvarfacs)

    return allvarfacs

end

"""
    $(SIGNATURES)
Build a deep subgraph copy from the DFG given a list of variables and factors and an optional distance.
Note: Orphaned factors (where the subgraph does not contain all the related variables) are not returned.
Related:
- [`copyGraph!`](@ref)
- [`getNeighborhood`](@ref)
- [`deepcopyGraph`](@ref)
- [`mergeGraph!`](@ref)
Dev Notes
- Bulk vs node for node: a list of labels are compiled and the sugraph is copied in bulk.
"""
function buildSubgraph(::Type{G},
                       dfg::AbstractDFG,
                       variableFactorLabels::Vector{Symbol},
                       distance::Int=0;
                       solvable::Int=0,
                       sessionId::String = "",
                       kwargs...) where G <: AbstractDFG

    if sessionId == ""
        sessionId = getSessionId(dfg) * "_sub_$(string(uuid4())[1:6])"
    end

    #build up the neighborhood from variableFactorLabels
    allvarfacs = getNeighborhood(dfg, variableFactorLabels, distance; solvable=solvable)

    variableLabels = intersect(listVariables(dfg), allvarfacs)
    factorLabels = intersect(listFactors(dfg), allvarfacs)
    # Copy the section of graph we want
    destDFG = deepcopyGraph(G, dfg, variableLabels, factorLabels; sessionId=sessionId, kwargs...)
    return destDFG
end

function buildSubgraph(dfg::AbstractDFG,
                       variableFactorLabels::Vector{Symbol},
                       distance::Int=0;
                       kwargs...)
    return buildSubgraph(DefaultDFG, dfg, variableFactorLabels, distance; kwargs...)
end

"""
    $(SIGNATURES)
Merger sourceDFG to destDFG given an optional list of variables and factors and distance.
Notes:
- Nodes already in the destination graph are updated from sourceDFG.
- Orphaned factors (where the subgraph does not contain all the related variables) are not included.
Related:
- [`copyGraph!`](@ref)
- [`buildSubgraph`](@ref)
- [`getNeighborhood`](@ref)
- [`deepcopyGraph`](@ref)
"""
function mergeGraph!(destDFG::AbstractDFG,
                     sourceDFG::AbstractDFG,
                     variableLabels::Vector{Symbol} = ls(sourceDFG),
                     factorLabels::Vector{Symbol} = lsf(sourceDFG),
                     distance::Int = 0;
                     solvable::Int = 0,
                     kwargs...)

    # find neighbors at distance to add
    allvarfacs = getNeighborhood(sourceDFG, union(variableLabels, factorLabels), distance; solvable=solvable)

    sourceVariables = intersect(listVariables(sourceDFG), allvarfacs)
    sourceFactors = intersect(listFactors(sourceDFG), allvarfacs)

    copyGraph!(destDFG, sourceDFG, sourceVariables, sourceFactors; deepcopyNodes=true, overwriteDest=true, kwargs...)

    return destDFG
end

##==============================================================================
## Variable Data: VND and PPE
##==============================================================================

#TODO API
"""
    $(SIGNATURES)
Merges and updates solver and estimate data for a variable (variable can be from another graph).
Note: Makes a copy of the estimates and solver data so that there is no coupling between graphs.
"""
function mergeVariableData!(dfg::AbstractDFG, sourceVariable::AbstractDFGVariable)

    var = getVariable(dfg, sourceVariable.label)

    mergePPEs!(var, sourceVariable)
    # If this variable has solverDataDict (summaries do not)
    :solverDataDict in fieldnames(typeof(var)) && mergeVariableSolverData!(var, sourceVariable)

    #update if its not a InMemoryDFGTypes, otherwise it was a reference
    # if satelite nodes are used it can be updated seprarately
    # !(isa(dfg, InMemoryDFGTypes)) && updateVariable!(dfg, var)

    return var
end

#TODO API
"""
    $(SIGNATURES)
Common function to update all solver data and estimates from one graph to another.
This should be used to push local solve data back into a cloud graph, for example.

Notes
- Returns `::Nothing`
"""
function mergeGraphVariableData!(destDFG::H, sourceDFG::G, varSyms::Vector{Symbol}) where {G <: AbstractDFG, H <: AbstractDFG}
    # Update all variables in the destination
    # (For now... we may change this soon)
    for variableId in varSyms
        mergeVariableData!(destDFG, getVariable(sourceDFG, variableId))
    end
end

##==============================================================================
## Graphs Structures (Abstract, overwrite for performance)
##==============================================================================
"""
    $(SIGNATURES)
Get a matrix indicating relationships between variables and factors. Rows are
all factors, columns are all variables, and each cell contains either nothing or
the symbol of the relating factor. The first row and first column are factor and
variable headings respectively.
Note:
- rather use getBiadjacencyMatrix
- Returns either of `::Matrix{Union{Nothing, Symbol}}`
"""
function getAdjacencyMatrixSymbols(dfg::AbstractDFG; solvable::Int=0)
    #
    varLabels = sort(map(v->v.label, getVariables(dfg, solvable=solvable)))
    factLabels = sort(map(f->f.label, getFactors(dfg, solvable=solvable)))
    vDict = Dict(varLabels .=> [1:length(varLabels)...].+1)

    adjMat = Matrix{Union{Nothing, Symbol}}(nothing, length(factLabels)+1, length(varLabels)+1)
    # Set row/col headings
    adjMat[2:end, 1] = factLabels
    adjMat[1, 2:end] = varLabels
    for (fIndex, factLabel) in enumerate(factLabels)
        factVars = getNeighbors(dfg, getFactor(dfg, factLabel), solvable=solvable)
        map(vLabel -> adjMat[fIndex+1,vDict[vLabel]] = factLabel, factVars)
    end
    return adjMat
end

# TODO API name get seems wrong maybe just biadjacencyMatrix
"""
    $(SIGNATURES)
Get a matrix indicating adjacency between variables and factors. Returned as
a named tuple: B::SparseMatrixCSC{Int}, varLabels::Vector{Symbol)
facLabels::Vector{Symbol). Rows are the factors, columns are the variables,
with the corresponding labels in varLabels,facLabels.

Notes
-  Returns `::NamedTuple{(:B, :varLabels, :facLabels), Tuple{SparseMatrixCSC, Vector{Symbol}, Vector{Symbol}}}`
"""
function getBiadjacencyMatrix(dfg::AbstractDFG; solvable::Int=0)
    varLabels = map(v->v.label, getVariables(dfg, solvable=solvable))
    factLabels = map(f->f.label, getFactors(dfg, solvable=solvable))

    vDict = Dict(varLabels .=> [1:length(varLabels)...])

    adjMat = spzeros(Int, length(factLabels), length(varLabels))

    for (fIndex, factLabel) in enumerate(factLabels)
        factVars = getNeighbors(dfg, getFactor(dfg, factLabel), solvable=solvable)
        map(vLabel -> adjMat[fIndex,vDict[vLabel]] = 1, factVars)
    end
    return (B=adjMat, varLabels=varLabels, facLabels=factLabels)
end

##==============================================================================
## DOT Files, falls back to LightDFG dot functions
##==============================================================================
"""
    $(SIGNATURES)
Produces a dot-format of the graph for visualization.

Notes
- Returns `::String`
"""
function toDot(dfg::AbstractDFG)
    #convert to LightDFG
    ldfg = LightDFG{NoSolverParams}()
    copyGraph!(ldfg, dfg, listVariables(dfg), listFactors(dfg))
    return toDot(ldfg)
end

"""
    $(SIGNATURES)
Produces a dot file of the graph for visualization.
Download XDot to see the data

Note
- Default location "/tmp/dfg.dot" -- MIGHT BE REMOVED
- Can be viewed with the `xdot` system application.
- Based on graphviz.org
"""
function toDotFile(dfg::AbstractDFG, fileName::String="/tmp/dfg.dot")

    #convert to LightDFG
    ldfg = LightDFG{NoSolverParams}()
    copyGraph!(ldfg, dfg, listVariables(dfg), listFactors(dfg))

    return toDotFile(ldfg, fileName)
end

##==============================================================================
## Summaries
##==============================================================================

"""
    $(SIGNATURES)
Get a summary of the graph (first-class citizens of variables and factors).
Returns a DFGSummary.

Notes
- Returns `::DFGSummary`
"""
function getSummary(dfg::G) where {G <: AbstractDFG}
    vars = map(v -> DFGVariableSummary(v), getVariables(dfg))
    facts = map(f -> DFGFactorSummary(f), getFactors(dfg))
    return DFGSummary(
        Dict(map(v->v.label, vars) .=> vars),
        Dict(map(f->f.label, facts) .=> facts),
        dfg.userId,
        dfg.robotId,
        dfg.sessionId)
end

"""
$(SIGNATURES)
Get a summary graph (first-class citizens of variables and factors) with the same structure as the original graph.

Notes
- this is a copy of the original.
- Returns `::LightDFG{NoSolverParams, DFGVariableSummary, DFGFactorSummary}`
"""
function getSummaryGraph(dfg::G) where {G <: AbstractDFG}
    summaryDfg = LightDFG{NoSolverParams, DFGVariableSummary, DFGFactorSummary}(
        description="Summary of $(getDescription(dfg))",
        userId=dfg.userId,
        robotId=dfg.robotId,
        sessionId=dfg.sessionId)
    deepcopyGraph!(summaryDfg, dfg)
    # for v in getVariables(dfg)
    #     newV = addVariable!(summaryDfg, DFGVariableSummary(v))
    # end
    # for f in getFactors(dfg)
    #     addFactor!(summaryDfg, getNeighbors(dfg, f), DFGFactorSummary(f))
    # end
    return summaryDfg
end


text:
"""
    $(SIGNATURES)

Check a scalar: should be real, finite, not nan, in optional bounds.
Real is implied by being `AbstractFloat`.
"""
function check_float(x :: T1; lb = nothing, ub = nothing) where T1 <: AbstractFloat
    isValid = true;
    if !isfinite(x)
        isValid = false;
        @warn "Not finite"
    end
    if !isnothing(lb)  &&  (x < lb - 1e-8)
        isValid = false;
        @warn "Below $lb"
    end
    if !isnothing(ub)  &&  (x > ub + 1e-8)
        isValid = false;
        @warn "Above $ub"
    end
    if !isValid
        println("Input value: $x");
    end
    return isValid
end

# Catches cases where the input is not `AbstractFloat`
function check_float(x; kwargs...)
    @warn "Not an AbstractFloat: $x"
    return false
end

"""
	$(SIGNATURES)

Check that a numeric array is finite, not NaN, and inside bounds.
"""
function check_float_array(m :: AbstractArray{T1},
    lb :: T2, ub :: T2) where {T1 <: AbstractFloat, T2 <: AbstractFloat}

    isValid = true;

    anyInf = any(x -> !isfinite(x), m);
    if anyInf
        isValid = false;
        nInf = sum(x -> !isfinite(x), m);
        @warn "$nInf elements not finite"
    end

    validLb = !isfinite(lb)  ||  all(x -> x >= lb, m);
    if !validLb
        isValid = false;
        mMin, idxMin = findmin(m);
        @warn """
            Array lower bound violated: $mMin < $lb
            at  $idxMin
            """
    end

    validUb = !isfinite(ub)  ||  all(x -> x <= ub, m);
    if !validUb
        isValid = false;
        mMax, idxMax = findmax(m);
        @warn """
            Array upper bound violated: $mMax < $ub
            at  $idxMax
            """
    end

    return isValid
end

"""
	$(SIGNATURES)

Check that an array is monotone along one dimension (weak or strict).

Based on `diff` in `Base`.
"""
function is_monotone(a :: AbstractArray{T,N}, d :: Integer;
    increasing, strict, atol = zero(T)) where {T,N}

    1 <= d <= N || throw(ArgumentError("dimension $d out of range (1:$N)"));

    r = axes(a);
    r0 = ntuple(i -> i == d ? UnitRange(1, last(r[i]) - 1) : UnitRange(r[i]), N);
    r1 = ntuple(i -> i == d ? UnitRange(2, last(r[i])) : UnitRange(r[i]), N);

    v0 = view(a, r0...);
    v1 = view(a, r1...);

    if increasing
        strict  ?  (op = >)  :  (op = >=);
    else
        strict  ?  (op = <)  :  (op = <=);
    end

    success = true;
    for idx in eachindex(v1)
        # E.g.: v1-v0 > -1e-8
        op(v1[idx] - v0[idx], -atol)  ||  (success = false; break);
    end
    return success
end

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

ext: jl


text:
# ---
# jupyter:
#   jupytext:
#     formats: jl:light,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.1
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # MintsHelper: Generating 1- and 2-electron Integrals with  Psi4
#
# In all of quantum chemistry, one process which is common to nearly every method is the evaluation of one-
# and two-electron integrals.  Fortunately, we can leverage infrastructure in
# Psi4  to perform this task for us.  This tutorial will discuss the [``psi4.core.MintsHelper``](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.MintsHelper "Go to API") class, which is an
# interface for the powerful Psi4 ``libmints`` library which wraps the `libint` library, where these integrals are actually computed.
#
# ## MintsHelper Overview
# In order to compute 1- and 2-electron integrals, we first need a molecule and basis set with which to work.  So,
# before diving into `MintsHelper`, we need to build these objects.  In the cell below, we have imported
#  Psi4  and NumPy, defined a water molecule, and set the basis to
# cc-pVDZ.  We've also set the memory available to  Psi4, as well as
# defined a variable `numpy_memory` which we will discuss later.

# +
# ==> Setup <==
# Import statements
using PyCall: pyimport
psi4 = pyimport("psi4")
np   = pyimport("numpy")
using Formatting: printfmt

# Memory & Output file
psi4.set_memory(Int(2e9))
numpy_memory = 2
psi4.core.set_output_file("output.dat", false)

# Molecule definition
h2o = psi4.geometry("""
O
H 1 0.96
H 1 0.96 2 104.5
""")

# Basis Set
psi4.set_options(Dict("basis" => "cc-pvdz"))
# -

# Now, we are ready to create an instance of the `MintsHelper` class.  To do this, we need to pass a `BasisSet`
# object to the `MintsHelper` initializer.  Fortunately, from the previous tutorial on the `Wavefunction` class, we know
# that we can obtain such an object from an existing wavefunction.  So, let's build a new wavefunction for our molecule,
# get the basis set object, and build an instance of `MintsHelper`:

# +
# ==> Build MintsHelper Instance <==
# Build new wavefunction
wfn = psi4.core.Wavefunction.build(h2o, psi4.core.get_global_option("basis"))

# Initialize MintsHelper with wavefunction's basis set
mints = psi4.core.MintsHelper(wfn.basisset())
# -

# Below are summarized several commonly computed quantities and how to obtain them using a `MintsHelper` class method:
#
# | Quantity | Function | Description |
# |----------|----------|-------------|
# | AO Overlap integrals | [mints.ao_overlap()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.MintsHelper.ao_overlap "Go to Documentation") | Returns AO overlap matrix as a `psi4.core.Matrix` object |
# | AO Kinetic Energy | [mints.ao_kinetic()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.MintsHelper.ao_kinetic "Go to Documentation") | Returns AO kinetic energy matrix as a `psi4.core.Matrix` object |
# | AO Potential Energy | [mints.ao_potential()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.MintsHelper.ao_potential "Go to Documentation") | Returns AO potential energy matrix as a `psi4.core.Matrix` object |
# | AO Electron Repulsion Integrals | [mints.ao_eri()](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.MintsHelper.ao_eri "Go to Documentation") | Returns AO electron repulsion integrals as a `psi4.core.Matrix` object

# As discussed previously, any of these `psi4.core.Matrix` objects can be accessed as Julia arrays, which is the preferred
# method in Psi4Julia.  For a Psi4 matrix `A`, we can access a Julia view using `psi4view(A)`, or we can make a
# copy of the matrix using `np.array(A)`.  This works as one would expect, converting square matrices into arrays of Array{Float64,2} type, for the overlap (S), kinetic energy (T), and potential energy (V) matrices.  In Psi4, the electron repulsion integrals
# (ERIs) are handled somewhat differently; `mints.ao_eri()` returns the rank-4 ERI tensor packed into a 2D matrix.  If the
# four indices of the ERI are p, q, r, s, then this element of the Psi4 Matrix can be accessed by first computing composite
# indices `pq = p * nbf + q` and `rs = r * nbf + s`, and then accessing element `A.get(pq,rs)`.  However, for convenience,
# the Julia view is a rank-4 tensor, and a particular ERI is more simply accessed like this:
# ~~~julia
# I = mints.ao_eri()
# I = psi4view(I)
# val = I[p,q,r,s]
# ~~~

function psi4view(psi4matrix)
   # Assumes Float64 type, C ordering
   if !hasproperty(psi4matrix,:__array_interface__)
      throw(ArgumentError("Input matrix cannot be accessed. Cannot be an rvalue"))
   end
   array_interface = psi4matrix.__array_interface__
   array_interface["data"][2] == false   || @warn "Not writable"
   array_interface["strides"] == nothing || @warn "Different ordering than C"
   array_interface["typestr"] == "<f8"   || @warn "Not little-endian Float64 eltype"
   ptr   = array_interface["data"][1]
   shape = reverse(array_interface["shape"])
   ndims = length(shape)
   unsafe_wrap(Array{Float64,ndims}, Ptr{Float64}(ptr), shape)
end

# In addition to these methods, another which is worth mentioning is the `MintsHelper.mo_eri()` ([Go to documentation](http://psicode.org/psi4manual/master/psi4api.html#psi4.core.MintsHelper.mo_eri)) function, which can transform
# the four-index, two-electron repulsion integrals from the atomic orbital (AO) to the molecular orbital (MO) basis,
# which will be important in MP2 theory.

# ## Memory Considerations

# Before moving forward to computing any 1- or 2-electron integrals, we must first discuss the memory requirements of
# these objects.  Whenever these quantities are computed, they are stored directly in memory (a.k.a. RAM,
# *not* on the hard drive) which, for a typical laptop or personal computer, usually tops out at around 16 GB of
# space.  The storage space required by the two-index AO overlap integrals and four-index ERIs scales as ${\cal O}(N^2)$
# and ${\cal O}(N^4)$, respectively, where $N$ is the number of AO basis functions.  This means that for a
# system with 500 AO basis functions, while the AO overlap integrals will only require 1 MB of memory to store,
# the ERIs will require a staggering **500 GB** of memory!! This can be reduced to **62.5 GB** of memory if integral permutational symmetry is used.
# However, this complicates the bookkeeping, and is not used in the `mints` functions discussed above.  For this reason, as well as the steep computational
# scaling of many of the methods demonstrated here, we limit ourselves to small systems ($\sim50$ basis functions)
# which should not require such egregious amounts of memory.  Additionally, we will employ a "memory check" to catch
# any case which could potentially try to use more memory than is available:
# ~~~julia
# # Memory check for ERI tensor
# I_size = nbf^4 * 8.e-9
# printfmt("Size of the ERI tensor will be {:4.2f} GB.\n", I_size)
# memory_footprint = I_size * 1.5
# if I_size > numpy_memory
#     psi4.core.clean()
#     throw(OutOfMemoryError("Estimated memory utilization ($memory_footprint GB) exceeds allotted memory " *
#                            "limit of $numpy_memory GB."))
# ~~~
# In this example, we have somewhat arbitrarily assumed that whatever other matrices we may need, in total their memory
# requirement will not exceed 50% of the size of the ERIs (hence, the total memory footprint of `I_size * 1.5`)
# Using the `numpy_memory` variable, we are able to control whether the ERIs will be computed, based on the amount of
# memory required to store them.
#
# **NOTE: DO NOT EXCEED YOUR SYSTEM'S MEMORY.  THIS MAY RESULT IN YOUR PROGRAM AND/OR COMPUTER CRASHING!**

# ## Examples: AO Overlap, AO ERIs, Core Hamiltonian
# The cell below demonstrates obtaining the AO overlap integrals, conducting the
# above memory check, and computing the ERIs and core Hamiltonian matrix for our water molecule.

# +
# ==> Integrals galore! <==
# AO Overlap
S = mints.ao_overlap()
S = psi4view(S) # psi4view(mints.ao_overlap()) would fail, first assign to a variable

# Number of basis functions
nbf = size(S)[1]

# Memory check
I_size = nbf^4 * 8.e-9
printfmt("Size of the ERI tensor will be {:4.2f} GB.\n", I_size)
memory_footprint = I_size * 1.5
if I_size > numpy_memory
    psi4.core.clean()
    throw(OutOfMemoryError("Estimated memory utilization ($memory_footprint GB) exceeds allotted memory " *
                           "limit of $numpy_memory GB."))
end

# Compute AO-basis ERIs
I = mints.ao_eri()

# Compute AO Core Hamiltonian
T = mints.ao_kinetic()
V = mints.ao_potential()
T = psi4view(T)
V = psi4view(V)
H = T + V


text:
using DocumentationGenerator

packages = DocumentationGenerator.installable_on_version(VERSION)

max_packages = if isempty(ARGS)
    length(packages)
else
    parse(Int, ARGS[1])
end

processes = parse(Int, get(ENV, "NUM_PKG_PROCESSES", "8"))
docspath = normpath(joinpath(@__DIR__, ".."))

DocumentationGenerator.build_documentation(collect(values(packages))[1:max_packages],
                                            processes = processes,
                                            basepath = docspath,
                                            update_only=true)


text:
type HeuristicSolver <: OPSolver end

function solve_op(solver::HeuristicSolver, op)
    return simple_heur(op)
end

type PathOptions
    path::Vector{Int}
    available::IntSet
    remaining::Float64
end

function simple_heur(op)
    d = op.distances
    start = op.start
    stop = op.stop
    # pick the furthest feasible point away
    feasible = collect(filter(j->within_range(op, j), 1:length(op)))
    if length(feasible) == 2 # only contains start and stop
        return [start, stop]
    end
    imax = indmax([d[start,i]+d[stop,i] for i in feasible])
    path = [start, feasible[imax], stop]
    available = setdiff(IntSet(filter(i->op.r[i]>0.0,1:length(op))),IntSet(path))
    remaining = op.distance_limit - distance(op, path)
    po = PathOptions(path, available, remaining)
    while insert_min!(po, d) end
    @assert abs(distance(op, path) + po.remaining - op.distance_limit) < 1e-5
    trade_better!(po, d, op.r, max_tries=div(length(available),2), op=op)
    @assert abs(distance(op, path) + po.remaining - op.distance_limit) < 1e-5
    return path
end

# start simple
# XXX rolled this into simple_heur
# function construct_cheap(op::SimpleOP)
#     d = op.distances
#     start = op.start
#     stop = op.stop
#     # pick the furthest feasible point away
#     feasible = collect(filter(j->within_range(op, j), 1:length(op)))
#     if length(feasible) == 2 # only contains start and stop
#         return [start, stop]
#     end
#     imax = indmax([d[start,i]+d[stop,i] for i in feasible])
#     path = [start, feasible[imax], stop]
#     available = setdiff(IntSet(1:length(op)),IntSet(path))
#     remaining = op.distance_limit - distance(op, path)
#     po = PathOptions(path, available, remaining)
#     while insert_min!(po, d) end
#     return path
# end

# insert the node with the minimum cost of insertion (regardless of price)
function insert_min!(po::PathOptions, d)
    p = po.path # for conciseness
    mindelta = po.remaining
    imin = 0
    jmin = 0
    for j in po.available
        for i in 2:length(p)
            delta = d[p[i-1],j] + d[j,p[i]] - d[p[i-1],p[i]]
            if delta <= mindelta
                imin = i
                jmin = j
                mindelta = delta
            end
        end
    end
    if jmin > 0
        delete!(po.available, jmin)
        insert!(p, imin, jmin)
        po.remaining -= mindelta
        return true
    end
    return false
end

# XXX get rid of op
function trade_better!(po::PathOptions, d, r; max_tries=typemax(Int), op=nothing)
    p = po.path # for conciseness
    potential = copy(po.available)
    trade = false
    tries = 0
    delta(i,b) = d[p[i-1],b] + d[b,p[i]] - d[p[i-1],p[i]] # change in distance if b is inserted
    deldelta(p,i) = d[p[i-1],p[i+1]] - (d[p[i-1],p[i]] + d[p[i],p[i+1]]) # change in distance if node at i is deleted
    while !isempty(potential) && tries <= max_tries
        b = reduce((i,j)->r[i]>=r[j]?i:j, potential)
        # find best place to insert
        i = reduce((i,j)->delta(i,b)>=delta(j,b)?i:j, 2:length(p))
        # insert
        bdelta = delta(i,b)
        insert!(p, i, b)
        # see if we can get back under budget by removing a cheaper node
        # first find all the removals that would put us back under budget
        under_budget = filter(i->po.remaining-bdelta-deldelta(p,i)>=0, 2:length(p)-1)
        # next find the cheapest of these removals
        j = reduce((i,j)->r[p[i]]<r[p[j]]?i:j, under_budget)
        # if this is not the one we just added, we have made a good trade
        if p[j] != b
            trade = true
            delete!(po.available, b)
            push!(po.available, p[j])
            po.remaining -= (bdelta+deldelta(p,j))
        end
        deleteat!(p,j)
        delete!(potential, b)
        tries += 1
    end
    return trade
end


text:
include("NumField/QQ.jl")
include("NumField/Elem.jl")
include("NumField/Hilbert.jl")
include("NumField/Field.jl")
include("NumField/Subfields.jl")
include("NumField/SimpleNumField.jl")
include("NumField/NonSimpleNumField.jl")
include("NumField/NfAbs.jl")
include("NumField/NfRel.jl")
include("NumField/CM.jl")
include("NumField/promote_rules.jl")


text:
# 1-D CNSOLT
using Random

function cnsoltValidConfigSet1D(seed = 610763893)
      rng = Random.seed!(seed)

      allcfgs = vec([((df,), (ord,), nch) for nch in 2:16, ord in 0:8, df in 1:4])
      cfgsTypeI  = filter(x ->  iseven(x[3]), allcfgs)
      cfgsTypeI  = randsubseq(cfgsTypeI,  15 / length( cfgsTypeI))
      cfgsTypeII = filter(x -> !iseven(x[3]), allcfgs)
      cfgsTypeII = randsubseq(cfgsTypeII, 20 / length(cfgsTypeII))

      cfgs = vcat(cfgsTypeI..., cfgsTypeII...)

      filter(cfgs) do cfg
            df, ord, nch = cfg
            flag = true
            flag = flag && prod(df) <= nch
            flag = flag && !(isodd(nch) && any(isodd.(ord)))
      end
end

# 2-D CNSOLT
function cnsoltValidConfigSet2D(seed = 625122358)
      rng = Random.seed!(seed)

      allcfgs = vec([(crdf.I, crord.I .- 1, nch) for nch in 2:20, crord in CartesianIndices((4,4) .+ 1), crdf in CartesianIndices((4,4))])

      cfgsTypeI  = filter(x ->  iseven(x[3]), allcfgs)
      cfgsTypeI  = randsubseq(cfgsTypeI,  30 / length( cfgsTypeI))
      cfgsTypeII = filter(x -> !iseven(x[3]), allcfgs)
      cfgsTypeII = randsubseq(cfgsTypeII, 40 / length(cfgsTypeII))

      cfgs = vcat(cfgsTypeI..., cfgsTypeII...)

      filter(cfgs) do cfg
            df, ord, nch = cfg
            flag = true
            flag = flag && prod(df) <= nch
            flag = flag && !(isodd(nch) && any(isodd.(ord)))
      end
end

# 3-D CNSOLT
function cnsoltValidConfigSet3D(seed = 229053284)
      rng = Random.seed!(seed)

      allcfgs = vec([(crdf.I, crord.I .- 1, nch) for nch in 2:10, crord in CartesianIndices((2,2,2) .+ 1), crdf in CartesianIndices((2,2,2))])

      cfgsTypeI  = filter(x ->  iseven(x[3]), allcfgs)
      cfgsTypeI  = randsubseq(cfgsTypeI,  20 / length( cfgsTypeI))
      cfgsTypeII = filter(x -> !iseven(x[3]), allcfgs)
      cfgsTypeII = randsubseq(cfgsTypeII, 40 / length(cfgsTypeII))

      cfgs = vcat(cfgsTypeI..., cfgsTypeII...)

      filter(cfgs) do cfg
            df, ord, nch = cfg
            flag = true
            flag = flag && prod(df) <= nch
            flag = flag && !(isodd(nch) && any(isodd.(ord)))
      end
end

# 1-D RNSOLT
function rnsoltValidConfigSet1D(seed = 75389923)
      rng = Random.seed!(seed)

      allcfgs = vec([((df,), (ord,), crnch.I) for crnch in CartesianIndices((10,10)), ord in 0:6, df in 1:4])

      cfgsTypeI  = filter(x -> x[3][1] == x[3][2], allcfgs)
      cfgsTypeI  = randsubseq(cfgsTypeI,  17 / length( cfgsTypeI))
      cfgsTypeII = filter(x -> x[3][1] != x[3][2], allcfgs)
      cfgsTypeII = randsubseq(cfgsTypeII, 20 / length(cfgsTypeII))

      cfgs = vcat(cfgsTypeI..., cfgsTypeII...)

      filter(cfgs) do cfg
            df, ord, nch = cfg
            flag = true
            flag = flag && prod(df) <= sum(nch)
            flag = flag && ((nch[1] == nch[2]) || all(iseven.(ord)))
            flag = flag && (cld(prod(df),2) <= nch[1] <= sum(nch) - fld(prod(df),2)) && (fld(prod(df),2) <= nch[2] <= sum(nch) - cld(prod(df),2))
      end
end

# 2-D RNSOLT
function rnsoltValidConfigSet2D(seed = 623122358)
      rng = Random.seed!(seed)

      allcfgs = vec([(crdf.I, crord.I .- 1, crnch.I) for crnch in CartesianIndices((10,10)), crord in CartesianIndices((4,4) .+ 1), crdf in CartesianIndices((4,4))])

      cfgsTypeI  = filter(x -> x[3][1] == x[3][2], allcfgs)
      cfgsTypeI  = randsubseq(cfgsTypeI,  20 / length( cfgsTypeI))
      cfgsTypeII = filter(x -> x[3][1] != x[3][2], allcfgs)
      cfgsTypeII = randsubseq(cfgsTypeII, 60 / length(cfgsTypeII))

      cfgs = vcat(cfgsTypeI..., cfgsTypeII...)

      filter(cfgs) do cfg
            df, ord, nch = cfg
            flag = true
            flag = flag && prod(df) <= sum(nch)
            flag = flag && ((nch[1] == nch[2]) || all(iseven.(ord)))
            flag = flag && (cld(prod(df),2) <= nch[1] <= sum(nch) - fld(prod(df),2)) && (fld(prod(df),2) <= nch[2] <= sum(nch) - cld(prod(df),2))
      end
end

# 3-D CNSOLT
function rnsoltValidConfigSet3D(seed = 385792992)
      rng = Random.seed!(seed)

      allcfgs = vec([(crdf.I, crord.I .- 1, crnch.I) for crnch in CartesianIndices((6,6)), crord in CartesianIndices((2,2,2) .+ 1), crdf in CartesianIndices((2,2,2))])

      cfgsTypeI  = filter(x -> x[3][1] == x[3][2], allcfgs)
      cfgsTypeI  = randsubseq(cfgsTypeI,  20 / length( cfgsTypeI))
      cfgsTypeII = filter(x -> x[3][1] != x[3][2], allcfgs)
      cfgsTypeII = randsubseq(cfgsTypeII, 70 / length(cfgsTypeII))

      cfgs = vcat(cfgsTypeI..., cfgsTypeII...)

      filter(cfgs) do cfg
            df, ord, nch = cfg
            flag = true
            flag = flag && prod(df) <= sum(nch)
            flag = flag && ((nch[1] == nch[2]) || all(iseven.(ord)))
            flag = flag && (cld(prod(df),2) <= nch[1] <= sum(nch) - fld(prod(df),2)) && (fld(prod(df),2) <= nch[2] <= sum(nch) - cld(prod(df),2))
      end
end

nothing

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
abstract type AbstractPoint{FT} end

abstract type Abstract1DPoint{FT} <: AbstractPoint{FT} end
abstract type Abstract2DPoint{FT} <: AbstractPoint{FT} end
abstract type Abstract3DPoint{FT} <: AbstractPoint{FT} end

"""
    @pointtype name fieldname1 ...

Define a subtype `name` of `AbstractPoint` with appropriate conversion functions.
"""
macro pointtype(name, fields...)
    if length(fields) == 1
        supertype = :(Abstract1DPoint{FT})
    elseif length(fields) == 2
        supertype = :(Abstract2DPoint{FT})
    elseif length(fields) == 3
        supertype = :(Abstract3DPoint{FT})
    end
    esc(
        quote
            struct $name{FT} <: $supertype
                $([:($field::FT) for field in fields]...)
            end
            Base.convert(::Type{$name{FT}}, pt2::$name) where {FT} =
                $name{FT}($((:(pt2.$field) for field in fields)...))
            $name{FT}(pt::$name) where {FT} = convert($name{FT}, pt)
            Base.eltype(::$name{FT}) where {FT} = FT
            Base.eltype(::Type{$name{FT}}) where {FT} = FT
            unionalltype(::$name{FT}) where {FT} = $name
            unionalltype(::Type{$name{FT}}) where {FT} = $name
            Base.promote_rule(
                ::Type{$name{FT1}},
                ::Type{$name{FT2}},
            ) where {FT1, FT2} = $name{promote_type(FT1, FT2)}
        end,
    )
end

@pointtype XPoint x
@pointtype YPoint y
@pointtype ZPoint z

@pointtype XYPoint x y
@pointtype XZPoint x z

@pointtype XYZPoint x y z

@pointtype Cartesian1Point x1
@pointtype Cartesian2Point x2
@pointtype Cartesian3Point x3

@pointtype Cartesian12Point x1 x2
@pointtype Cartesian123Point x1 x2 x3

const Cartesian2DPoint = Cartesian12Point
@deprecate Cartesian2DPoint Cartesian12Point

const Cartesian3DPoint = Cartesian123Point
@deprecate Cartesian3DPoint Cartesian123Point

product_coordinates(xp::XPoint, yp::YPoint) = XYPoint(promote(xp.x, yp.y)...)
product_coordinates(xp::XPoint, zp::ZPoint) = XZPoint(promote(xp.x, zp.z)...)

product_coordinates(xyp::XYPoint, zp::ZPoint) =
    XYZPoint(promote(xyp.x, xyp.y, zp.z)...)

# TODO: get rid of these and refactor to consistent point types
ZPoint(pt::Cartesian3Point{FT}) where {FT} = ZPoint{FT}(pt.x3)
XYPoint(pt::Cartesian12Point{FT}) where {FT} = XYPoint{FT}(pt.x1, pt.x2)
XYZPoint(pt::Cartesian123Point{FT}) where {FT} =
    XYZPoint{FT}(pt.x1, pt.x2, pt.x3)

Cartesian3Point(pt::ZPoint{FT}) where {FT} = Cartesian3Point{FT}(pt.z)
Cartesian12Point(pt::XYPoint{FT}) where {FT} = Cartesian12Point{FT}(pt.x, pt.y)
Cartesian123Point(pt::XYZPoint{FT}) where {FT} =
    Cartesian123Point{FT}(pt.x, pt.y, pt.z)

component(p::AbstractPoint{FT}, i::Symbol) where {FT} = getfield(p, i)::FT
component(p::AbstractPoint{FT}, i::Integer) where {FT} = getfield(p, i)::FT

ncomponents(p::AbstractPoint) = nfields(p)
components(p::AbstractPoint) =
    SVector(ntuple(i -> component(p, i), ncomponents(p)))

_coordinate_type(ptyp::Type{Abstract1DPoint}, ::Val{1}) = ptyp
_coordinate(p::Abstract1DPoint, ::Val{1}) = p

_coordinate_type(::Type{XYPoint{FT}}, ::Val{1}) where {FT} = XPoint{FT}
_coordinate_type(::Type{XYPoint{FT}}, ::Val{2}) where {FT} = YPoint{FT}
_coordinate(p::XYPoint, ::Val{1}) = XPoint(p.x)
_coordinate(p::XYPoint, ::Val{2}) = YPoint(p.y)

_coordinate_type(::Type{XZPoint{FT}}, ::Val{1}) where {FT} = XPoint{FT}
_coordinate_type(::Type{XZPoint{FT}}, ::Val{2}) where {FT} = ZPoint{FT}
_coordinate(p::XZPoint, ::Val{1}) = XPoint(p.x)
_coordinate(p::XZPoint, ::Val{2}) = ZPoint(p.z)

_coordinate_type(::Type{XYZPoint{FT}}, ::Val{1}) where {FT} = XPoint{FT}
_coordinate_type(::Type{XYZPoint{FT}}, ::Val{2}) where {FT} = YPoint{FT}
_coordinate_type(::Type{XYZPoint{FT}}, ::Val{3}) where {FT} = ZPoint{FT}
_coordinate(p::XYZPoint, ::Val{1}) = XPoint(p.x)
_coordinate(p::XYZPoint, ::Val{2}) = YPoint(p.y)
_coordinate(p::XYZPoint, ::Val{3}) = ZPoint(p.z)

coordinate_type(ptyp::Type{<:AbstractPoint}, ax::Int) =
    _coordinate_type(ptyp, Val(ax))
coordinate_type(ptyp::Type{<:AbstractPoint}, ax::Integer) =
    _coordinate_type(ptyp, Int(ax))

coordinate(pt::AbstractPoint, ax::Int) = _coordinate(pt, Val(ax))
coordinate(pt::AbstractPoint, ax::Integer) = _coordinate(pt, Int(ax))

# TODO: we need to rationalize point operations with vectors
#Base.:(-)(p1::T, p2::T) where {T <: AbstractPoint} =
#    unionalltype(T)((components(p1) - components(p2))...)

# the following are needed for linranges to work correctly with coordinate values
Base.:(+)(p1::T, p2::T) where {T <: AbstractPoint} =
    unionalltype(T)((components(p1) + components(p2))...)
Base.:(*)(p::T, x::Number) where {T <: AbstractPoint} =
    unionalltype(T)((components(p) * x)...)
Base.:(*)(x::Number, p::AbstractPoint) = p * x
Base.:(/)(p::T, x::Number) where {T <: AbstractPoint} =
    unionalltype(T)((components(p) / x)...)

function Base.:(==)(p1::T, p2::T) where {T <: AbstractPoint}
    return components(p1) == components(p2)
end

function Base.isapprox(p1::T, p2::T; kwargs...) where {T <: AbstractPoint}
    return isapprox(components(p1), components(p2); kwargs...)
end

"""
    interpolate(coords::NTuple{4}, ξ1, ξ2)

Interpolate between `coords` by parameters `ξ1, ξ2` in the interval `[-1,1]`.

The type of interpolation is determined by the element type of `coords`:
- `SVector`: use bilinear interpolation
"""
function interpolate(coords::NTuple{4, V}, ξ1, ξ2) where {V <: Abstract2DPoint}
    VV = unionalltype(V)
    VV(
        ((1 - ξ1) * (1 - ξ2)) / 4 * component(coords[1], 1) +
        ((1 + ξ1) * (1 - ξ2)) / 4 * component(coords[2], 1) +
        ((1 - ξ1) * (1 + ξ2)) / 4 * component(coords[3], 1) +
        ((1 + ξ1) * (1 + ξ2)) / 4 * component(coords[4], 1),
        ((1 - ξ1) * (1 - ξ2)) / 4 * component(coords[1], 2) +
        ((1 + ξ1) * (1 - ξ2)) / 4 * component(coords[2], 2) +
        ((1 - ξ1) * (1 + ξ2)) / 4 * component(coords[3], 2) +
        ((1 + ξ1) * (1 + ξ2)) / 4 * component(coords[4], 2),
    )
end

"""
    interpolate(coords::NTuple{2}, ξ1)

Interpolate between `coords` by parameters `ξ1` in the interval `[-1,1]`.
The type of interpolation is determined by the element type of `coords`
"""
function interpolate(coords::NTuple{2, V}, ξ1) where {V <: Abstract1DPoint}
    VV = unionalltype(V)
    VV(
        (1 - ξ1) / 2 * component(coords[1], 1) +
        (1 + ξ1) / 2 * component(coords[2], 1),
    )
end


text:
using Bio3DView
using Test
using BioStructures

testfile(path::AbstractString...) = normpath(
                    dirname(pathof(Bio3DView)), "..", "examples", path...)

style = Style("sphere", Dict("colorscheme"=> "greenCarbon"))
surface = Surface(Dict("colorscheme"=> "greenCarbon"))
isosurface = IsoSurface(testfile("benzene.cube"), 0.01)
box = Box([0.0, 0.0, 0.0], [6.0, 6.0, 6.0])
lines = [Line([0.0, 1.0, 2.0], [6.0, 7.0, 7.0]), Line([-3.0, -2.0, -1.0], [3.0, 4.0, 5.0])]
cylinders = [
    Cylinder([0.0, 1.0, 2.0], [6.0, 7.0, 7.0]),
    Cylinder([-3.0, -2.0, -1.0], [3.0, 4.0, 5.0]; startcap=RoundCap, stopcap=FlatCap),
]
axes = Axes(5, 0.3)
cameraangle = CameraAngle(0, 0, 0, 0, 0, 0, 0, 1)

viewfile(testfile("benzene.sdf"), "sdf", html=true)
viewpdb("1CRN", html=true)
viewfile(testfile("benzene.sdf"), style=style, html=true)

s = """
ATOM     72  N   ALA A  11      16.899  42.259  22.187  1.00 16.83           N
ATOM     73  CA  ALA A  11      15.960  42.201  23.284  1.00 18.31           C
ATOM     74  C   ALA A  11      15.625  43.630  23.738  1.00 17.96           C
ATOM     75  O   ALA A  11      14.821  43.804  24.675  1.00 22.53           O
ATOM     76  CB  ALA A  11      16.528  41.416  24.561  1.00 15.72           C
"""
viewstring(s, "pdb", style=Style("sphere"), html=true)

struc = read(testfile("1AKE.pdb"), PDB)
viewstruc(struc['A'], html=true)
viewstruc(struc, disorderselector, style=Style("sphere"), html=true)
viewstruc(struc['A'], surface=surface, html=true)

viewfile(testfile("benzene.sdf"), "sdf", isosurface=isosurface, html=true)
viewfile(testfile("benzene.sdf"), "sdf", box=box, html=true)
viewfile(testfile("benzene.sdf"), "sdf", lines=lines, html=true)
viewfile(testfile("benzene.sdf"), "sdf", lines=lines[1], html=true)
viewfile(testfile("benzene.sdf"), "sdf", cylinders=cylinders, html=true)
viewfile(testfile("benzene.sdf"), "sdf", cylinders=cylinders[1], html=true)

viewfile(testfile("IRMOF-1.xyz"), "xyz"; style=Style("stick"), html=true)
isosurface = IsoSurface(testfile("IRMOF-1.cube"), -5.0)
viewfile(testfile("IRMOF-1.xyz"), "xyz", isosurface=isosurface,
            vtkcell=testfile("IRMOF-1.vtk"), html=true)

isosurface = IsoSurface(testfile("SBMOF-1.cube"), 100.0, wireframe=true,
                color="green")
viewfile(testfile("SBMOF-1.xyz"), "xyz"; isosurface=isosurface,
            vtkcell=testfile("SBMOF-1.vtk"), axes=axes, cameraangle=cameraangle, html=true)


text:
sz = (4,5)
filename = tempname()
#filename = "/tmp/test-10.nc"
# The mode "c" stands for creating a new file (clobber)

NCDatasets.NCDataset(filename,"c") do ds
    # define the dimension "lon" and "lat"
    ds.dim["lon"] = sz[1]
    ds.dim["lat"] = sz[2]

    @test haskey(ds.dim,"lon")
    @test ds.dim["lon"] == sz[1]
    @test ds.dim["lat"] == sz[2]

    forecast = NCDatasets.defGroup(ds,"forecast")
    v = NCDatasets.defVar(forecast,"var",Float64,("lon","lat"))
    v[:,:] = fill(Float64(123),size(v))
end

NCDatasets.NCDataset(filename) do ds
    @test haskey(ds.group,"forecast")

    forecast = ds.group["forecast"]

    @test all(forecast["var"][:,:] .== 123)

    s = IOBuffer()
    show(s,ds)
    @test occursin("Groups",String(take!(s)))

end

#@show NCDatasets.NCDataset(filename)


text:
using Random
using Plots
using SparseArrays
using GeoStats
using SpectralGaussianSimulation
using LaTeXStrings
#reset past font scaling
Plots.scalefontsizes()
#scale fonts
Plots.scalefontsizes(1.75)
#Extracting stored varibles

lightblue = colorant"rgb(63%,74%,78%)"
orange = colorant"rgb(85%,55%,13%)"
silver = colorant"rgb(69%,67%,66%)"
rust = colorant"rgb(72%,26%,6%)"

include("./CovFuncs.jl")
include("./Utils.jl")

Random.seed!(123)

n = 1000
N = n^2
NTest = 20000
NTrain = N - NTest

l = 0.1
covfunc!(r) = matern32!(r, l)
vario = MaternVariogram(range=l, order=1.5)
x, Y = specSampleGrid2D(n, 1, 1, vario)

α = 0.02
β = 0.02

# random shuffle of indices
shuffleIDs = shuffle(1:N)
# Selecting groups of nodes as test data
innerInds = [i for i in 1:N if (norm(x[:, i] .- 0.3) < 0.05 || norm((x[:,i] .- 0.75) .* [1.0, 0.1]) < 0.01)]
shuffleIDs = vcat(shuffleIDs, innerInds)
shuffleIDs = unique(shuffleIDs[end:-1:1])[end:-1:1]

# Creating data to be used by each run of for loop
#Reordering data
xBase = x[:, shuffleIDs]

n_extrapolation = min(NTest, length(innerInds))
interpolation_indices_base = (NTrain + 1) : (NTrain + NTest - n_extrapolation)
extrapolation_indices_base = (NTrain + NTest - n_extrapolation + 1) : N

x_train = xBase[:, rand(1 : NTrain, round(Int, α * NTrain))]
x_scattered = xBase[:, interpolation_indices_base[rand(1 : length(interpolation_indices_base), round(Int, β * length(interpolation_indices_base)))]]
x_region = xBase[:, extrapolation_indices_base[rand(1 : length(extrapolation_indices_base), round(Int, β * length(extrapolation_indices_base)))]]

outplot = plot(size=(400,400), aspect_ratio=:equal, xticks=false, yticks=false, axis=false, legend=:bottomright)
scatter!(outplot, vec(x_train[1,:]), vec(x_train[2,:]), color=silver, label=L"\mathrm{training}", markerstrokewidth=0.0, markeralpha=0.5)
scatter!(outplot, vec(x_scattered[1,:]), vec(x_scattered[2,:]), color=rust, label=L"\mathrm{scattered}", markerstrokewidth=0.0)
scatter!(outplot, vec(x_region[1,:]), vec(x_region[2,:]), color=orange, label=L"\mathrm{region}", markerstrokewidth=0.0)

savefig(outplot, "./out/plots/FFTPoints.pdf")


text:
module Queens

export  solveNQueensProblemRBFS,
        solveNQueensProblemGSU,
        solveNQueensProblemGSBF,
        solveNQueensProblemGSAS,
        solveNQueensProblemHCS,
        solveNQueensProblemSAS,
        solveNQueensProblemGAS,
        solveNQueensProblemGAS9

using AIMACore

import AIMACore: search, goal_test, step_cost, actions,
                 result, heuristic, state_value, successor_states,
                 mutate, reproduce, random_state
import Base: ==, show, convert, hash

const E_INVALID_INPUT_LENGTH = "The input parameter is of a wrong length"
const E_INVALID_INPUT = "The input parameter is invalid"
const E_INVALID_STATE = "The state is invalid"
const E_INVALID_ACTION = "The action is invalid"

const SIZE = 8

mutable struct Grid{T} <: State
    qloc::Vector{Int}
    function Grid{T}(loc::Vector{Int}) where T
        @assert length(loc) == T E_INVALID_INPUT_LENGTH
        new(loc)
    end
end

# This heuristic is not relevant as this is the same for every step
h(g::Grid{T}) where T = T - count(x->x > 0, g.qloc)

Grid(T::Int) = Grid{T}(zeros(Int, T))

function convert(::Type{BitArray{2}}, g::Grid{T}) where T
    m = BitArray{2}(undef, T, T)
    l = g.qloc
    for i = 1:T
        l[i] == 0 && continue
        m[l[i],i] = true
    end
    return m
end

BitArray{2}(g::Grid{T}) where T = convert(BitArray{2}, g)

==(g1::Grid, g2::Grid) = hash(g1) == hash(g2)

hash(g::Grid{T}, h::UInt) where T =
    xor(sum([UInt(g.qloc[i])*UInt(T+1)^UInt(i) for i=1:T]), h)

function show(io::IO, g::Grid{T}) where T
    m = BitArray{2}(g)
    a = [ m[i,j] ? 'Q' : 'x' for i=1:T,j=1:T]
    for j = 1:T-1
        show(io, a[j,:])
        println(io, "")
    end
    show(io, a[T,:])
end

struct Place <: Action
    location::Tuple{Int,Int}
end

mutable struct NQueensProblem{SA <: SearchAlgorithm, T} <: Problem
    initial_state::Grid{T}
    search_algorithm::SA
    h::Function
end

function result(problem::NQueensProblem,
                state::Grid{T}, action::Place) where T
    qloc = copy(state.qloc)
    x, y = action.location
    qloc[x] = y
    return Grid{T}(qloc)
end

step_cost(problem::NQueensProblem, state::Grid, action::Place) = 1

heuristic(problem::NQueensProblem, state::Grid) = problem.h(state)

goal_test(problem::NQueensProblem, state::Grid{T}) where T =
    count(x->x > 0, state.qloc) == T && state_value(problem, state) == div(T*(T-1), 2)

cannot_place(x, y, i, j) = x == i || y == j || x + j == y + i || x + y == i + j

function cannot_place(state::Grid{T}, i, j) where T
    for k = 1:T
        state.qloc[k] == 0 && continue
        cannot_place(k, state.qloc[k], i, j) && return true
    end
    return false
end

actions(problem::NQueensProblem, state::Grid{T}) where T =
    Place.([(i, j) for i=1:T for j=1:T if !cannot_place(state, i, j)])

function solveNQueensProblem(obj::SA,
    h::Function, init_state::Grid=Grid(SIZE)) where {SA <: SearchAlgorithm}
    problem = NQueensProblem(init_state, obj, h)
    path = search(problem)
    path isa Symbol && return path
    path isa State && return path
    ret=[]
    count = 0
    for iter in path
        push!(ret, iter.state)
    end
    return ret
end

const GSD = GraphSearchDepth(Grid(SIZE))
const GSU = GraphSearchUniformCost(Grid(SIZE))
const GSBF = GraphSearchBestFirst(Grid(SIZE))
const GSAS = GraphSearchAStar(Grid(SIZE))
const RBFS = RecursiveBestFirstSearch(Grid(SIZE))

solveNQueensProblemGSU() = solveNQueensProblem(GSU, h)
solveNQueensProblemGSBF() = solveNQueensProblem(GSBF, h)
solveNQueensProblemGSAS() = solveNQueensProblem(GSAS, h)
solveNQueensProblemRBFS() = solveNQueensProblem(RBFS, h)

function queen_pair_score(g::Grid{T}) where T
    l = g.qloc
    return div(T*(T-1),2) -
        sum([l[j] > 0 && l[i] > 0 && cannot_place(j, l[j], i, l[i]) for i=1:T for j=i+1:T])
end

random_state(T::Int) = Grid{T}([rand(1:T) for i = 1:T])

state_value(problem::NQueensProblem, state::Grid) = queen_pair_score(state)

function successor_states(problem::NQueensProblem, state::Grid{T}) where T
    v = Vector{Grid{T}}()
    for i = 1:T
        s = [j for j = 1:T if j != state.qloc[i]]
        for n in s
            l = copy(state.qloc)
            l[i] = n
            push!(v, Grid{T}(l))
        end
    end
    return v
end

HCS = HillClimbingSearch(Grid(SIZE), 100)
SAS = SimulatedAnnealingSearch(Grid(SIZE))

solveNQueensProblemHCS() = solveNQueensProblem(HCS, h, random_state(SIZE))
solveNQueensProblemSAS() = solveNQueensProblem(SAS, h, random_state(SIZE))

random_state(s::Grid{T}) where T = random_state(T)

const BIT_OPERATION = true

function to_bits(s1::Grid{8}; mutate=false)
    b1 = BitArray(undef, 24)
    l = s1.qloc
    it = 0
    for i = 1:8
        n = l[i] - 1
        b1[it+=1] = (n & 1 != 0)
        b1[it+=1] = (n & 2 != 0)
        b1[it+=1] = (n & 4 != 0)
    end
    if (mutate)
        c = rand(1:24)
        b1[c] = !b1[c]
    end
    return b1
end

function from_bits(b::BitArray)
    l = zeros(Int, 8)
    it = 0
    for i = 1:8
        n = 1
        n += b[it+=1] ? 1 : 0
        n += b[it+=1] ? 2 : 0
        n += b[it+=1] ? 4 : 0
        l[i] = n
    end
    return Grid{8}(l)
end

bit_mutate(s::Grid{8}) = to_bits(s, mutate=true) |> from_bits

function bit_reproduce(s1::Grid{8}, s2::Grid{8})
    b1 = to_bits(s1)
    b2 = to_bits(s2)
    c = rand(1:24)
    copyto!(b1[c+1:24], b2[c+1:24])
    return from_bits(b1)
end

function mutate(s::Grid{T}) where T
    BIT_OPERATION && T == 8 && return bit_mutate(s)
    l = copy(s.qloc)
    l[rand(1:T)] = rand(1:T)
    return Grid{T}(l)
end

function reproduce(s1::Grid{T}, s2::Grid{T}) where T
    BIT_OPERATION && T == 8 && return bit_reproduce(s1, s2)
    c = rand(1:T)
    l = copy(s1.qloc)
    copyto!(l[c+1:T], s2.qloc[c+1:T])
    return Grid{T}(l)
end

GAS = GeneticAlgorithmSearch(Grid(SIZE), queen_pair_score)
solveNQueensProblemGAS() = solveNQueensProblem(GAS, h, random_state(SIZE))

GAS9 = GeneticAlgorithmSearch(Grid(9), queen_pair_score)
solveNQueensProblemGAS9() = solveNQueensProblem(GAS9, h, random_state(9))

end


text:
using LibFaceDetection
using Test
using FileIO
using GeometryBasics

@testset "LibFaceDetection.jl" begin
    img = load(joinpath(@__DIR__, "faces.jpg"))
    faces = detect_faces(img)
    @test length(faces) == 12
    @testset "faces: $(i)" for (i, face) in enumerate(faces)
        # mainly test against corruption!
        bb = face.boundingbox
        @test all(p-> p in bb, face.points)
        @test bb in Rect(Vec2f0(0, 0), Vec2f0(reverse(size(img))))
    end
end


text:
using AWS
using AWS.AWSExceptions: AWSException
using Dates
using HTTP: Response, StatusError
using JuliaOnLambda
using Mocking
using Test

@service ECR
Mocking.activate()

include("patches.jl")

function _now_formatted()
    return lowercase(Dates.format(now(Dates.UTC), dateformat"yyyymmdd\THHMMSSsss\Z"))
end

function docker_rmi(name::AbstractString)
    return read(`docker rmi $(name)`)
end

@testset "JuliaOnLambda.jl" begin
    @testset "_get_create_ecr_repo" begin
        repo_name = "juliaonlambda-" * _now_formatted()

        @testset "repo exists" begin
            apply(describe_repository_exists_patch) do
                response = JuliaOnLambda._get_create_ecr_repo(repo_name)

                @test response == REPO_URL
            end
        end

        @testset "repo dne" begin
            apply([describe_repository_dne_patch, create_repository_patch]) do
                response = JuliaOnLambda._get_create_ecr_repo(repo_name)

                @test response == REPO_URL
            end
        end

        @testset "alternative error" begin
            apply([describe_repository_dne_patch, create_repository_failure_patch]) do
                @test_throws AWSException JuliaOnLambda._get_create_ecr_repo(repo_name)
            end
        end
    end

    @testset "_build_docker_image" begin
        image_name = "julia-on-lambda-test-"
        image_count = 0

        function get_image_name()
            image_count += 1
            return string(image_name, image_count)
        end

        function image_exists(image_name::AbstractString)
            result = read(`docker images --format "{{.Repository}}"`, String)
            result = string.(split(result, "\n"))

            return image_name in result
        end

        @testset "default kwargs" begin
            try
                image_name = get_image_name()
                JuliaOnLambda._build_docker_image(image_name)
                @test image_exists(image_name)
            finally
                docker_rmi(image_name)
            end
        end

        @testset "dockerfile_path set" begin
            dockerfile_path = joinpath(@__DIR__, "resources", "Dockerfile")

            try
                image_name = get_image_name()
                JuliaOnLambda._build_docker_image(
                    image_name; dockerfile_path=dockerfile_path
                )
                @test image_exists(image_name)
            finally
                docker_rmi(image_name)
            end
        end
    end

    @testset "_upload_docker_image" begin
        repo_name = "juliaonlambda-" * _now_formatted()
        image_name = "bizbaz"
        tag = "latest"

        try
            repo_uri = JuliaOnLambda._get_create_ecr_repo(repo_name)
            JuliaOnLambda._create_docker_image(image_name, tag)
            JuliaOnLambda._tag_docker_image(image_name, tag, repo_uri)
            JuliaOnLambda._upload_docker_image(repo_uri)

            resp = ECR.describe_images(repo_name)["imageDetails"]

            @test !isempty(resp)
        finally
            ECR.delete_repository(repo_name, Dict("force" => true))
            docker_rmi(image_name)
        end
    end
end


text:
using Revise
using Random
using RDatasets
using ForwardDiff
using ForwardDiff: Dual
partiali(f, x, i, args...) = ForwardDiff.partials(f([Dual{}(x[j], 1.0*(i==j)) for j in eachindex(x)], args...))[]

using ZigZagBoomerang
g(y, μ, σ) = exp(-(y-μ)^2/(2σ^2))/sqrt(2*π*σ^2) # Gaussian
lo(x) = exp(x)/(1+exp(x))

galaxies = dataset("mass", "galaxies")
Yobs = galaxies.x1/1000
const N = 82
@assert N == length(Yobs)
Random.seed!(1)

const K = 3
const d = 3K - 1
const σ = 1.0
const R = 10 # subsampling

function probs(x)
    if K == 2
        (lo(x[end]), 1-lo(x[end]))
    elseif K == 3
        a, b = exp(x[end]), exp(x[end-1])
        c = 1.0 + a + b
        (a/c, b/c, 1 - a/c - b/c)
    else
        a = exp.(x[end-K+2:end])
        sa = sum(a)
        [a/(1 + sa); 1 - sa/(1 + sa)]
    end
end
lprior(x) = sum(log(g(x[k], 1.0, 10.0)) for k in 1:d)
f(y, x, ps=probs(x)) = sum(ps[k]*g(y, x[k], σ*x[K+k]) for k in 1:K)
ϕ(x, Y) = -lprior(x) - N/R*sum(log(f(y, x)) for y in rand(Y, R))
∇ϕ(x, i, Y) = partiali(ϕ, x, i, Y)

using LinearAlgebra
using SparseArrays
using Statistics

t0 = 0.0
x0 = collect(1:K)*10.0
x0 = [x0; ones(d-K)]
θ0 = rand([-1.0,1.0], d)

Γ = sparse(I, d, d)

c = [20.0 for i in 1:d]

μ = zeros(d)

Z = LocalZigZag(Γ, μ)
T = 2000.0

@time trace, (tT, xT, θT), (acc, num), c = pdmp(∇ϕ, t0, x0, θ0, T, c, Z, Yobs; adapt=true)
@show acc, acc/num, mean(c)

tr = collect(trace)
ts = first.(tr)
xs = last.(tr)

ms = [getindex.(xs, k) for k in 1:d]
ps_ = [probs(x) for x in xs]
ps = [getindex.(ps_, k) for k in 1:K]

using Makie
using Bridge
col = Bridge.viridis(1:K)
p1 = scatter(T*ones(K), xT[1:K])
for k in 1:K
    band!(p1, ts, ms[k] + 1.96ms[k+K], ms[k] - 1.96ms[k+K], color=(col[k],0.3) )
    scatter!(p1, ts, ms[k], markersize=1.8, color=col[k] )
    lines!(p1, ts, ms[k],  color=col[k] )

end

p2 = lines(ts, ps[1], color=col[1])
for k in 2:K
    lines!(p2, ts, ps[k], color=col[k])
end

m = median.(ms)
r = range(0, 40, length=200)
p3 = lines(r, [f(y, m) for y in r])
linesegments!(p3, [repeat(Yobs, inner=2) repeat([0,0.05], outer=N)])

p = hbox(title(p3, "est. density and obs."), title(p2, "Trace p[k]"), title(p1, "Trace μ[k] ± σ[k]"))
save(joinpath("figures", "galaxy.png"), title(p, "Galaxy dataset"))


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

# # Natural frequency analysis of 3d frame structure

# For general information about Euler-Bernoulli beam theory, see
# [this](https://en.wikipedia.org/wiki/Euler%E2%80%93Bernoulli_beam_theory)
# wikipedia page.

# The model is a 3d frame, shown in picture.

# ![](3d_frame/model.png)

using JuliaFEM, LinearAlgebra

# Reading mesh

datadir = Pkg.dir("JuliaFEM", "examples", "3d_frame")
mesh = aster_read_mesh(joinpath(datadir, "model.med"))
println("Number of nodes in a model: ", length(mesh.nodes))

# Create beam elements. For 3d model, we need to define at least
# [Young's modulus](https://en.wikipedia.org/wiki/Young%27s_modulus),
# [shear modulus](https://en.wikipedia.org/wiki/Shear_modulus),
# [density](https://en.wikipedia.org/wiki/Density)
# cross-section area, moment of inertia in local coordinate
# system and polar moment of inertia.

beam_elements = create_elements(mesh, "FRAME")
@info("Number of elements: ", length(beam_elements))
update!(beam_elements, "youngs modulus", 210.0e6)
update!(beam_elements, "shear modulus", 84.0e6)
update!(beam_elements, "density", 7850.0e-3)
update!(beam_elements, "cross-section area", 20.0e-2)
update!(beam_elements, "torsional moment of inertia 1", 10.0e-5)
update!(beam_elements, "torsional moment of inertia 2", 10.0e-5)
update!(beam_elements, "polar moment of inertia", 30.0e-5)

# The direction of beam is defined in same way than in ABAQUS.
# That is, we have a tangent direction and one normal direction.
# The third direction is then cross product of tangent and normal.
# Because the second area moment is same in both directions, we can
# choose normal direction freely.

for element in beam_elements
    X1, X2 = element("geometry", 0.0)
    t = (X2-X1)/norm(X2-X1)
    I = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
    k = indmax([norm(cross(t, I[:,k])) for k in 1:3])
    n = cross(t, I[:,k])/norm(cross(t, I[:,k]))
    update!(element, "normal", n)
end

# Create boundary conditions: fix all degrees of freedom for nodes in
# a set FIXED. Here we first create elements of type `Poi1` for each
# node j in set FIXED, update geometry field and then create new fields
# `fixed displacmeent 1`, `fixed displacement 2`, and so on, where the
# displacement / rotation is prescribed.

bc_elements = [Element(Poi1, [j]) for j in mesh.node_sets[:FIXED]]
update!(bc_elements, "geometry", mesh.nodes)
for i=1:3
    update!(bc_elements, "fixed displacement $i", 0.0)
    update!(bc_elements, "fixed rotation $i", 0.0)
end

# Create a problem, containing beam elements and boundary conditions:

frame = Problem(Beam, "3d frame", 6)
add_elements!(frame, beam_elements)
add_elements!(frame, bc_elements)

# Perform modal analysis

analysis = Analysis(Modal)
xdmf = Xdmf(joinpath(datadir, "3d_frame_results"); overwrite=true)
add_results_writer!(analysis, xdmf)
add_problems!(analysis, frame)
run!(analysis)
close(xdmf)

# Each `Analysis` can have properties, e.g. time, maximum number of iterations,
# convergence tolerance and so on. Eigenvalues of calculation are stored as a
# properties of analysis:

freqs = sqrt.(step.properties.eigvals) / (2*pi)
println("Natural frequencies [Hz]: $(round.(freqs, 2))")

# [![mode5](3d_frame/natfreq.png)](https://www.youtube.com/watch?v=GzktCqeASmo)


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Test
using Behavior.Gherkin
using Behavior.Gherkin: ScenarioStep
using Behavior
using Behavior: StepDefinitionContext, StepDefinition, StepDefinitionLocation
using Behavior: Executor, StepExecutionResult, QuietRealTimePresenter, executefeature
import Behavior: present

@testset "Feature Executor     " begin
    @testset "Execute a feature; Feature has one scenario; Feature result has one scenario result" begin
        presenter = QuietRealTimePresenter()
        given = Given("some precondition")
        matcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition))
        scenario = Scenario("some scenario", String[], ScenarioStep[given])
        featureheader = FeatureHeader("Some feature", [], [])
        feature = Feature(featureheader, [scenario])
        executor = Executor(matcher, presenter)

        featureresult = executefeature(executor, feature)

        @test length(featureresult.scenarioresults) == 1
    end

    @testset "Execute a feature; Feature has two scenarios; Result has two scenario results" begin
        presenter = QuietRealTimePresenter()
        given = Given("some precondition")
        matcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition))
        scenario1 = Scenario("some scenario", String[], ScenarioStep[given])
        scenario2 = Scenario("some other scenario", String[], ScenarioStep[given])
        featureheader = FeatureHeader("Some feature", [], [])
        feature = Feature(featureheader, [scenario1, scenario2])
        executor = Executor(matcher, presenter)

        featureresult = executefeature(executor, feature)

        @test length(featureresult.scenarioresults) == 2
    end

    @testset "Execute a feature; Feature has three scenarios; Scenarios are executed in order" begin
        presenter = QuietRealTimePresenter()
        given = Given("some precondition")
        matcher = FakeStepDefinitionMatcher(Dict(given => successful_step_definition))
        scenario1 = Scenario("some scenario", String[], ScenarioStep[given])
        scenario2 = Scenario("some other scenario", String[], ScenarioStep[given])
        scenario3 = Scenario("some third scenario", String[], ScenarioStep[given])
        featureheader = FeatureHeader("Some feature", [], [])
        feature = Feature(featureheader, [scenario1, scenario2, scenario3])
        executor = Executor(matcher, presenter)

        featureresult = executefeature(executor, feature)

        @test featureresult.scenarioresults[1].scenario == scenario1
        @test featureresult.scenarioresults[2].scenario == scenario2
        @test featureresult.scenarioresults[3].scenario == scenario3
    end

    @testset "Execute a feature; Feature has one failing scenario; Scenario result has a failing step" begin
        presenter = QuietRealTimePresenter()
        given = Given("some precondition")
        matcher = FakeStepDefinitionMatcher(Dict(given => failed_step_definition))
        scenario = Scenario("some scenario", String[], ScenarioStep[given])
        featureheader = FeatureHeader("Some feature", [], [])
        feature = Feature(featureheader, [scenario])
        executor = Executor(matcher, presenter)

        featureresult = executefeature(executor, feature)

        @test featureresult.scenarioresults[1].steps[1] isa Behavior.StepFailed
    end

    @testset "Execute a feature; One Scenario and an Outline with two examples; Three results" begin
        presenter = QuietRealTimePresenter()

        step1 = Given("step 1")
        step2 = Given("step 2")
        scenario = Scenario("some scenario", String[], ScenarioStep[step1])
        outline = ScenarioOutline("", String[], ScenarioStep[Given("step <stepnumber>")], ["stepnumber"], ["1" "2"])
        matcher = FakeStepDefinitionMatcher(Dict(step1 => successful_step_definition,
                                                 step2 => successful_step_definition))

        featureheader = FeatureHeader("Some feature", [], [])
        feature = Feature(featureheader, [scenario, outline])
        executor = Executor(matcher, presenter)

        featureresult = executefeature(executor, feature)

        @test length(featureresult.scenarioresults) == 3
    end
end