text:
using AccretionFormulae
include("../../arbitrary_f/alt-emissivity-deviations.jl")

a = 0.97
M = 1.0

# epsilon
# creating arbitrary function
fs, r_range = init()
fs_epsneg, r_range_epsneg = init(ϵ_3=-2)
fs_epspos, r_range_epspos = init(ϵ_3=2)

r_vals = collect(r_range)
pop!(r_vals)

f_vals_analytic = AccretionFormulae.flux.(r_range, a, M)
f_vals_analytic ./= maximum(f_vals_analytic)
f_vals_approx = []
f_vals_approx_epsneg = []
f_vals_approx_epspos = []
for r in r_range
    push!(f_vals_approx, f_approx_func(r, fs, r_range, a=a, M=M))
    push!(f_vals_approx_epsneg, f_approx_func(r, fs_epsneg, r_range_epsneg, a=a, M=M))
    push!(f_vals_approx_epspos, f_approx_func(r, fs_epspos, r_range_epspos, a=a, M=M))
end

plt = plot(r_range, f_vals_analytic, label="Analytic Solution", grid=false, framestyle=:box, ylims=(0,1.05), xlabel="Radius (M)", ylabel="f (Arbitrary Units)")
plot!(r_range, f_vals_approx, label="ϵ₃=0", linestyle=:dash)
plot!(r_range, f_vals_approx_epsneg, label="ϵ₃=-2")
plot!(r_range, f_vals_approx_epspos, label="ϵ₃=2")

png(plt, "emission_comparison_ana_arb_eps.png")

# alpha 22
# creating arbitrary function
fs, r_range = init()
fs_epsneg, r_range_epsneg = init(α_22 =-5)
fs_epspos, r_range_epspos = init(α_22 =5)

r_vals = collect(r_range)
pop!(r_vals)

f_vals_approx = []
f_vals_approx_epsneg = []
f_vals_approx_epspos = []
for r in r_range
    push!(f_vals_approx, f_approx_func(r, fs, r_range, a=a, M=M))
    push!(f_vals_approx_epsneg, f_approx_func(r, fs_epsneg, r_range_epsneg, a=a, M=M))
    push!(f_vals_approx_epspos, f_approx_func(r, fs_epspos, r_range_epspos, a=a, M=M))
end

plt = plot(r_range, f_vals_analytic, label="Analytic Solution", grid=false, framestyle=:box, ylims=(0,1.05), xlabel="Radius (M)", ylabel="f (Arbitrary Units)")
plot!(r_range, f_vals_approx, label="α₂₂=0", linestyle=:dash)
plot!(r_range, f_vals_approx_epsneg, label="α₂₂=-5")
plot!(r_range, f_vals_approx_epspos, label="α₂₂=5")

png(plt, "emission_comparison_ana_arb_alpha22.png")

# alpha 13
# creating arbitrary function
fs, r_range = init()
fs_epsneg, r_range_epsneg = init(α_13 =-2)
fs_epspos, r_range_epspos = init(α_13 =2)

r_vals = collect(r_range)
pop!(r_vals)

f_vals_approx = []
f_vals_approx_epsneg = []
f_vals_approx_epspos = []
for r in r_range
    push!(f_vals_approx, f_approx_func(r, fs, r_range, a=a, M=M))
    push!(f_vals_approx_epsneg, f_approx_func(r, fs_epsneg, r_range_epsneg, a=a, M=M))
    push!(f_vals_approx_epspos, f_approx_func(r, fs_epspos, r_range_epspos, a=a, M=M))
end

plt = plot(r_range, f_vals_analytic, label="Analytic Solution", grid=false, framestyle=:box, ylims=(0,1.05), xlabel="Radius (M)", ylabel="f (Arbitrary Units)")
plot!(r_range, f_vals_approx, label="α₁₃ =0", linestyle=:dash)
plot!(r_range, f_vals_approx_epsneg, label="α₁₃=-2")
plot!(r_range, f_vals_approx_epspos, label="α₁₃=2")

png(plt, "emission_comparison_ana_arb_alpha13.png")

# alpha 52
# creating arbitrary function
fs, r_range = init()
fs_epsneg, r_range_epsneg = init(α_52 =-2)
fs_epspos, r_range_epspos = init(α_52 =2)

r_vals = collect(r_range)
pop!(r_vals)

f_vals_approx = []
f_vals_approx_epsneg = []
f_vals_approx_epspos = []
for r in r_range
    push!(f_vals_approx, f_approx_func(r, fs, r_range, a=a, M=M))
    push!(f_vals_approx_epsneg, f_approx_func(r, fs_epsneg, r_range_epsneg, a=a, M=M))
    push!(f_vals_approx_epspos, f_approx_func(r, fs_epspos, r_range_epspos, a=a, M=M))
end

plt = plot(r_range, f_vals_analytic, label="Analytic Solution", grid=false, framestyle=:box, ylims=(0,1.05), xlabel="Radius (M)", ylabel="f (Arbitrary Units)")
plot!(r_range, f_vals_approx, label="α₅₂ =0", linestyle=:dash)
plot!(r_range, f_vals_approx_epsneg, label="α₅₂ =-2")
plot!(r_range, f_vals_approx_epspos, label="α₅₂ =2")

png(plt, "emission_comparison_ana_arb_alpha52.png")


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.

using Base.Meta
# multivarate "dictionary" used for collections of variables/constraints

abstract JuMPContainer{T,N}

include("JuMPArray.jl")

type IndexPair
    idxvar
    idxset
end

type JuMPDict{T,N} <: JuMPContainer{T,N}
    tupledict::Dict{NTuple{N,Any},T}
    meta::Dict{Symbol,Any}

    JuMPDict() = new(Dict{NTuple{N,Any},T}(), Dict{Symbol,Any}())
end

function JuMPDict{T,N}(d::Dict{NTuple{N,Any},T})
    tmp = JuMPDict{T,N}()
    tmp.tupledict = d
    tmp
end

function JuMPDict{T,N}(d::Dict{NTuple{N,Any},T}, meta::Dict{Symbol,Any})
    tmp = JuMPDict{T,N}()
    tmp.tupledict = d
    tmp.meta = meta
    tmp
end

type JuMPContainerData
    name
    indexsets
    indexexprs::Vector{IndexPair}
    condition
end

#JuMPDict{T,N}(name::AbstractString) =
#    JuMPDict{T,N}(Dict{NTuple{N},T}(), name)

Base.getindex(d::JuMPDict, t...) = d.tupledict[t]
Base.setindex!(d::JuMPDict, value, t...) = (d.tupledict[t] = value)

Base.map(f::Function, d::JuMPArray) =
    JuMPArray(map(f, d.innerArray), d.indexsets, d.lookup, d.meta)

function Base.map{T,N}(f::Function, d::JuMPDict{T,N})
    ret = Base.return_types(f, Tuple{T})
    R = (length(ret) == 1 ? ret[1] : Any)
    x = JuMPDict{R,N}()
    for (k,v) in d.tupledict
        x.tupledict[k] = f(v)
    end
    return x
end

Base.isempty(d::JuMPContainer) = isempty(_innercontainer(d))

coloncheck(args::Number...) = nothing

function coloncheck(args...)
    if any(t -> t == Colon(), args)
        error("Colons not allowed as keys in JuMP containers.")
    end
end

# generate and instantiate a type which is indexed by the given index sets
# the following types of index sets are allowed:
# 0:K -- range with compile-time starting index
# S -- general iterable set
function gendict(instancename,T,idxsets...)
    N = length(idxsets)
    truearray = true
    for idxset in idxsets
        s = isexpr(idxset,:escape) ? idxset.args[1] : idxset
        if !(isexpr(s,:(:)) && length(s.args) == 2 && s.args[1] == 1)
            truearray = false
            break
        end
    end
    sizes = Expr(:tuple, [:(length($rng)) for rng in idxsets]...)
    if truearray
        :($instancename = Array($T, $sizes))
    else
        indexsets = Expr(:tuple, idxsets...)
        :($instancename = JuMPArray(Array($T, $sizes), $indexsets))
    end
end

pushmeta!(x::JuMPContainer, sym::Symbol, val) = (x.meta[sym] = val)
getmeta(x::JuMPContainer, sym::Symbol) = x.meta[sym]

# duck typing approach -- if eltype(innerArray) doesn't support accessor, will fail
for accessor in (:getdual, :getlowerbound, :getupperbound)
    @eval $accessor(x::Union{JuMPContainer,Array}) = map($accessor,x)
end

_similar(x::Array) = Array(Float64,size(x))
_similar{T}(x::Dict{T}) = Dict{T,Float64}()

_innercontainer(x::JuMPArray) = x.innerArray
_innercontainer(x::JuMPDict)  = x.tupledict

function _getValueInner(x)
    vars = _innercontainer(x)
    vals = _similar(vars)
    data = printdata(x)
    warnedyet = false
    for I in eachindex(vars)
        tmp = _getValue(vars[I])
        if isnan(tmp) && !warnedyet
            warn("Variable value not defined for entry of $(data.name). Check that the model was properly solved.")
            warnedyet = true
        end
        vals[I] = tmp
    end
    vals
end

JuMPContainer_from(x::JuMPDict,inner) = JuMPDict(inner)
JuMPContainer_from(x::JuMPArray,inner) = JuMPArray(inner, x.indexsets)

function getvalue(x::JuMPContainer)
    getvalue_warn(x)
    ret = JuMPContainer_from(x,_getValueInner(x))
    # I guess copy!(::Dict, ::Dict) isn't defined, so...
    for (key,val) in x.meta
        ret.meta[key] = val
    end
    m = getmeta(x, :model)
    # cache indexing info for new container for printing purposes
    m.varData[ret] = printdata(x)
    ret
end

# delegate zero-argument functions
for f in (:(Base.ndims), :(Base.length), :(Base.abs))
    @eval $f(x::JuMPArray) = $f(x.innerArray)
end

Base.length(x::JuMPDict) = length(x.tupledict)

Base.ndims{T,N}(x::JuMPDict{T,N}) = N
Base.abs(x::JuMPDict) = map(abs, x)
# avoid dangerous behavior with "end" (#730)
Base.endof(x::JuMPArray) = error("endof() (and \"end\" syntax) not implemented for JuMPArray objects.")
Base.size(x::JuMPArray) = error("size (and \"end\" syntax) not implemented for JuMPArray objects. Use JuMP.size if you want to access the dimensions.")
Base.size(x::JuMPArray,k) = error("size (and \"end\" syntax) not implemented for JuMPArray objects. Use JuMP.size if you want to access the dimensions.")
size(x::JuMPArray) = size(x.innerArray)
size(x::JuMPArray,k) = size(x.innerArray,k)
# for uses of size() within JuMP
size(x) = Base.size(x)
size(x,k) = Base.size(x,k)
# delegate one-argument functions
Compat.issymmetric(x::JuMPArray) = Compat.issymmetric(x.innerArray)

Base.eltype{T}(x::JuMPContainer{T}) = T

Base.full(x::JuMPContainer) = x

# keys/vals iterations for JuMPContainers
Base.keys(d::JuMPDict)    = keys(d.tupledict)
Base.values(d::JuMPDict)  = values(d.tupledict)

Base.keys(d::JuMPArray)   = KeyIterator(d)
Base.values(d::JuMPArray) = ValueIterator(d.innerArray)

# Wrapper type so that you can't access the values directly
type ValueIterator{T,N}
    x::Array{T,N}
end
Base.start(it::ValueIterator)   =  start(it.x)
Base.next(it::ValueIterator, k) =   next(it.x, k)
Base.done(it::ValueIterator, k) =   done(it.x, k)
Base.length(it::ValueIterator)  = length(it.x)

type KeyIterator{JA<:JuMPArray}
    x::JA
end
Base.start(it::KeyIterator)   =  start(it.x.innerArray)
@generated __next{T,N,NT}(x::JuMPArray{T,N,NT}, k) =
    quote
        subidx = ind2sub(size(x),k)
        $(Expr(:tuple, [:(x.indexsets[$i][subidx[$i]]) for i in 1:N]...)), next(x.innerArray,k)[2]
    end
Base.next(it::KeyIterator, k) = __next(it.x,k)
Base.done(it::KeyIterator, k) =   done(it.x.innerArray, k)
Base.length(it::KeyIterator)  = length(it.x.innerArray)

ext: jl


text:
# @testset "transform" begin

#     # Tests for Circularised Kernels.
#     @test lb(Circularised(EQ(), 0.0, 1.0)) == 0.0
#     @test ub(Circularised(EQ(), 0.0, 1.0)) == 1.0
#     @test Circularised(EQ(), 0.0, 1.0) == Circularised(EQ(), 0.0, 1.0)
#     @test Circularised(RQ(1.0), 0.0, 1.0) != Circularised(EQ(), 0.0, 1.0)
#     @test Circularised(EQ(), 0.5, 1.0) != Circularised(EQ(), 0.0, 1.0)
#     @test Circularised(EQ(), 0.0, 0.5) != Circularised(EQ(), 0.0, 1.0)
#     @test Circularised(EQ(), 0.0, 1.0)(0.0, 0.0) == 1.0
#     @test Circularised(EQ(), 0.0, 1.0)(0.0, 0.0) == EQ()(0.0, 0.0)

#     # Test some invariances.
#     let k = Circularised(EQ(), 0.0, 1.0)
#         @test k(0.0, 0.2) ≈ k(0.8, 1.0)
#         @test k(0.2, 0.0) ≈ k(0.0, 0.2)
#         @test k(0.0, 0.2) ≈ k(1.0, 0.8)
#         @test k(0.0, 1e-5) ≈ k(1.0 - 1e-5, 1.0)
#     end
# end


text:
# DESCRIPTION
# This file contains the aircraft parameters. Always re-run DI_design.jl after
# modifying the aircraft parameters.

# ------------------------- MASS AND AERO PARAMETERS ---------------------------

# CG location
XCG=0.35
# reference CG location
XCGR=0.35
# wing area [ft²]
S=300.0
# wing span [ft]
B=30.0
# wing chord [ft]
CBAR=11.32
# weight [lb]
WEIGHT=20500.0
# momentum of engine turbine [lb-ft]
HE=160.0
# inertia constants
C1=-0.770
C2=0.02755
C3=1.055E-4
C4= 1.642E-6
C5=0.9604
C6=1.759E-2
C7=1.792E-5
C8=-0.7336
C9=1.587E-5
# longitudinal location of vertical accelerometer [ft]
XA=15.0

# ------------------------ MISCELLANEOUS CONSTANTS------------------------------

# radians to degrees
R2D=180.0/pi
# degrees to radians
D2R=pi/180.0
# gravitational acceleration [ft²]
G=32.17

# ------------------- TRIM AND LINEARIZATION PARAMETERS ------------------------

# state vector perturbations
DELXLIN=[0.1;0.001;0.001;0.001;0.001;0.001;0.001;0.001;0.001;0.1;0.1;0.1;0.1]
# control vector perturbations
DELCLIN=[0.1;0.1;0.1;0.1]
# trim target tolerance
TOL=1e-5*[1.0;D2R;D2R;D2R;D2R;D2R;D2R;D2R;D2R;1.0;1.0;1.0;1.0;0.1]
# trim variables: [VT α β ϕ θ ψ P Q R pow controls]
TRIMVARS=[1:9; 13; 14; 15; 16; 17]
# trim targets
TRIMTARG=[collect(1:13); 15]
# number of states
NSTATES=13
# number of controls
NCTRLS=4
# number of outputs
NOUT=6


text:
function RangeRate(arg0::GroundStation, arg1::AbsoluteDate, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jboolean, arg6::ObservableSatellite)
    return RangeRate((GroundStation, AbsoluteDate, jdouble, jdouble, jdouble, jboolean, ObservableSatellite), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function get_station(obj::RangeRate)
    return jcall(obj, "getStation", GroundStation, ())
end

function is_two_way(obj::RangeRate)
    return jcall(obj, "isTwoWay", jboolean, ())
end


text:
using SafeTestsets

@safetestset GraphTopo = "Graph topology" begin include("graph.jl") end
@safetestset PacketDrop = "Packet drop"    begin include("packet_drop.jl") end
@safetestset FlowLimiter = "Flow limiter"   begin include("flow_limiter.jl") end
@safetestset BindStreams = "Bind streams" begin include("bind_streams.jl") end

# TODO: test that close() is called
# TODO: test error system once implemented


text:
#
println("Dcp) Test of the Photoionization and PlasmaShift module with ASF from an internally generated initial- and final-state multiplet.")
#
setDefaults("print summary: open", "zzz-PhotoIonization-Plasma.sum")
setDefaults("method: continuum, Galerkin")
setDefaults("method: normalization, pure sine")   ## setDefaults("method: normalization, pure Coulomb")    setDefaults("method: normalization, pure sine")

wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=JAC.Radial.Grid(false), nuclearModel=Nuclear.Model(36.),
                        initialConfigs=[Configuration("1s^2 2s^2 2p^6")],
                        finalConfigs  =[Configuration("1s^2 2s^2 2p^5"), Configuration("1s^2 2s 2p^6") ],
                        process = PhotoInPlasma(),
                        processSettings=PlasmaShift.PhotoSettings(PlasmaShift.DebyeHueckel(), 0.25, 0., 0, [E1], [JAC.UseCoulomb],
                                                                  [3000., 3200.], true, false, Tuple{Int64,Int64}[]) )

wb = perform(wa)
setDefaults("print summary: close", "")


text:
import MittagLefflerFunctions.MLF.mlf2
import MittagLefflerFunctions.MLF.radius

α = 0.5
β = 1.0
z = 8.7 * exp(Complex(0.0, 3π/4))
p = 10
@assert (   (  abs(angle(z)) > π*α )
         && ( abs(abs(angle(z)) - (π*α)) > 10.0^(-p) ) )

r0 = radius(α, β, z, p)
val = mlf2(α, β, z, p, r0)
ref = erfcx(-z)

@test abs(val-ref) ≤ abs(ref)


text:
"""
```
impulse_responses(m, system)

impulse_responses(system, horizon)
```

Compute impulse responses for a single draw.

### Inputs

- `m::AbstractModel`: model object
- `system::System{S}`: state-space system matrices
- `horizon::Int`: number of periods ahead to forecast

where `S<:AbstractFloat`

### Outputs

- `states::Array{S, 3}`: matrix of size `nstates` x `horizon` x `nshocks` of
  state impulse response functions
- `obs::Array{S, 3}`: matrix of size `nobs` x `horizon` x `nshocks` of
  observable impulse response functions
- `pseudo::Array{S, 3}`: matrix of size `npseudo` x `horizon` x `nshocks` of
  pseudo-observable impulse response functions
"""
function impulse_responses{S<:AbstractFloat}(m::AbstractModel, system::System{S})
    horizon = impulse_response_horizons(m)
    impulse_responses(system, horizon)
end

function impulse_responses{S<:AbstractFloat}(system::System{S}, horizon::Int;
                                             flip_shocks::Bool = false)
    # Setup
    nshocks      = size(system[:RRR], 2)
    nstates      = size(system[:TTT], 1)
    nobs         = size(system[:ZZ], 1)
    npseudo      = size(system[:ZZ_pseudo], 1)

    states = zeros(S, nstates, horizon, nshocks)
    obs    = zeros(S, nobs,    horizon, nshocks)
    pseudo = zeros(S, npseudo, horizon, nshocks)

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    z0 = zeros(S, nstates)

    for i = 1:nshocks
        # Isolate single shock
        shocks = zeros(S, nshocks, horizon)
        if flip_shocks
            shocks[i, 1] = sqrt(system[:QQ][i, i])  # a positive 1 s.d. shock
        else
            shocks[i, 1] = -sqrt(system[:QQ][i, i]) # a negative 1 s.d. shock
        end

        # Iterate state space forward
        states[:, :, i], obs[:, :, i], pseudo[:, :, i], _ = forecast(system, z0, shocks)
    end

    return states, obs, pseudo
end

# Method for specifying the subset of shocks, and the size of each shock
function impulse_responses{S<:AbstractFloat}(m::AbstractModel, system::System{S},
                                             horizon::Int, shock_names::Vector{Symbol},
                                             shock_values::Vector{Float64})

    # Must provide a name and value for each shock
    @assert length(shock_names) == length(shock_values)

    # Setup
    exo          = m.exogenous_shocks
    nshocks      = size(system[:RRR], 2)
    nstates      = size(system[:TTT], 1)
    nobs         = size(system[:ZZ], 1)
    npseudo      = size(system[:ZZ_pseudo], 1)

    states = zeros(S, nstates, horizon, nshocks)
    obs    = zeros(S, nobs,    horizon, nshocks)
    pseudo = zeros(S, npseudo, horizon, nshocks)

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    s_0 = zeros(S, nstates)

    for (i, shock) in enumerate(shock_names)
        # Isolate single shock
        shocks = zeros(S, nshocks, horizon)
        shocks[exo[shock], 1] = shock_values[i]

        # Iterate state space forward
        states[:, :, exo[shock]], obs[:, :, exo[shock]], pseudo[:, :, exo[shock]], _ = forecast(system, s_0, shocks)
    end

    return states, obs, pseudo
end

# Method for picking a specific shock and the size of the desired shift in the state or
# observed variable using that shock. (Omitting the feature to do pseudo-observables
# now, for simplicity)
function impulse_responses{S<:AbstractFloat}(m::AbstractModel, system::System{S},
                                             horizon::Int, shock_name::Symbol,
                                             var_name::Symbol, var_value::Float64)

    # Setup
    var_names, var_class =
    if var_name in keys(m.endogenous_states)
        m.endogenous_states, :states
    elseif var_name in keys(m.observables)
        m.observables, :obs
    end
    exo          = m.exogenous_shocks
    nshocks      = size(system[:RRR], 2)
    nstates      = size(system[:TTT], 1)
    nobs         = size(system[:ZZ], 1)
    npseudo      = size(system[:ZZ_pseudo], 1)

    states = zeros(S, nstates, horizon, nshocks)
    obs    = zeros(S, nobs,    horizon, nshocks)
    pseudo = zeros(S, npseudo, horizon, nshocks)

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    s_0 = zeros(S, nstates)

    # Isolate single shock
    shocks = zeros(S, nshocks, horizon)
    if var_class == :states
        shocks[exo[shock_name], 1] = obtain_shock_from_desired_state_value(var_value,
                                                                      var_names[var_name],
                                                                      exo[shock_name],
                                                                      system[:RRR])
    else # == :obs
        shocks[exo[shock_name], 1] = obtain_shock_from_desired_obs_value(var_value,
                                                                    var_names[var_name],
                                                                    exo[shock_name],
                                                                    system[:ZZ],
                                                                    system[:RRR])
    end

    # Iterate state space forward
    states[:, :, exo[shock_name]], obs[:, :, exo[shock_name]], pseudo[:, :, exo[shock_name]], _ = forecast(system, s_0, shocks)

    return states, obs, pseudo
end

# Given a desired value, x, for some state s^i_1, back out the necessary value of
# ϵ^j_1 needed s.t. s^i_1 = x.
# E.g. if I want the impulse response of output (s^i_1) to be 50 basis points
# based on a shock to productivity (ϵ^j_1), how large does that shock
# to productivity need to be?
# The answer is ϵ^j_1 = x/R_{i,j}
function obtain_shock_from_desired_state_value(state_value::Float64, state_ind::Int,
                                               shock_ind::Int, RRR::Matrix{Float64})
    return state_value/RRR[state_ind, shock_ind]
end

# Similar principle to the function for states
# The answer is ϵ^j_1 = x/{Σ_{k=1}^n Z_{i, k} R_{k, j}},
# where n is the total number of states, i is the index of the observable of interest
# y^i_1, and j is the index of the shock of interest, ϵ^j_1
function obtain_shock_from_desired_obs_value(obs_value::Float64, obs_ind::Int,
                                             shock_ind::Int, ZZ::Matrix{Float64},
                                             RRR::Matrix{Float64})
    return obs_value/dot(ZZ[obs_ind, :], RRR[:, shock_ind])
end


text:
# TODO docstring
struct PairDictionary{I, T, D <: AbstractDictionary{I, T}} <: AbstractDictionary{I, Pair{I, T}}
    d::D
end

"""
    pairs(dict::AbstractDictionary)

Return a new dictionary, wrapping `dict`, that shares the same `keys` but containing
key-value pairs.

# Example

```julia
julia> dict = HashDictionary(["a", "b", "c"], [1, 2, 3])
3-element HashDictionary{String,Int64}
 "c" │ 3
 "b" │ 2
 "a" │ 1

julia> pairs(dict)
3-element Dictionaries.PairDictionary{String,Int64,HashDictionary{String,Int64}}
 "c" │ "c" => 3
 "b" │ "b" => 2
 "a" │ "a" => 1
```
"""
Base.pairs(d::AbstractDictionary) = PairDictionary{keytype(d), eltype(d), typeof(d)}(d)

Base.parent(d::PairDictionary) = d.d
Base.keys(d::PairDictionary) = keys(parent(d))

# Length
Base.IteratorSize(d::PairDictionary) = Base.IteratorSize(parent(d))
Base.length(d::PairDictionary) = Base.length(parent(d))

# Standard interface
@propagate_inbounds Base.getindex(d::PairDictionary{I}, i::I) where {I} = i => parent(d)[i]
Base.isassigned(d::PairDictionary{I}, i::I) where {I} = isassigned(parent(d), i)

# Token interface
@propagate_inbounds gettoken(d::PairDictionary{I}, i::I) where {I} = gettoken(parent(d), i)

@propagate_inbounds function gettokenvalue(pd::PairDictionary, t)
    d = parent(pd)
    gettokenvalue(keys(d), t) => gettokenvalue(d, t)
end

istokenassigned(pd::PairDictionary, t) = istokenassigned(parent(pd), t)

iteratetoken(pd::PairDictionary, s...) = iteratetoken(parent(pd), s...)

# Factories
Base.similar(dict::PairDictionary, ::Type{T}, indices) where {T} = similar(parent(dict), T, indices)
Base.empty(dict::PairDictionary, ::Type{I}, ::Type{T}) where {I, T} = similar(parent(dict), I, T)


text:
using Test
using Documenter

using BioGenerics
using FormatSpecimens
using GenomicFeatures
using XAM

import BioAlignments: Alignment, AlignmentAnchor, OP_START, OP_MATCH, OP_DELETE
import BGZFStreams: BGZFStream
import BioGenerics.Exceptions: MissingFieldException
import BioSequences: @dna_str, @aa_str

# Generate a random range within `range`.
function randrange(range)
    x = rand(range)
    y = rand(range)
    if x < y
        return x:y
    else
        return y:x
    end
end

@testset "XAM" begin
    include("test_sam.jl")
    include("test_bam.jl")
    include("test_crosscheck.jl")

    # Include doctests.
    DocMeta.setdocmeta!(XAM, :DocTestSetup, :(using XAM); recursive=true)
    doctest(XAM; manual = false)
end


text:
using Documenter
using EnergySystemModeling

makedocs(
    sitename = "EnergySystemModeling.jl",
    format = Documenter.HTML(
        assets = ["assets/favicon.ico"]
    ),
    modules = [EnergySystemModeling],
    authors = "Lucas Condeixa, Fabricio Oliveira, Jaan Tollander de Balsch",
    pages = [
        "Home" => "index.md",
        "aggreg.md",
        "plotting.md",
        "datageneration.md",
        "api.md",
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/gamma-opt/EnergySystemModeling.jl.git"
)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "QuantReg"
version = v"0.1.2"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/fogarty-ben/QuantReg.jl-Builder.git", "7800e5d4e4c2dfd89e8d9076d764b11a51f55081")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd QuantReg.jl-Builder/
mkdir -p ${libdir}
gfortran -fPIC -shared -std=legacy rqbr.f -o ${libdir}/librqbr.${dlext}
gfortran -fPIC -shared -std=legacy rqfnb.f -o ${libdir}/librqfnb.${dlext}
install_license LICENSE
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_gfortran_versions(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("librqfnb", :librqfnb),
    LibraryProduct("librqbr", :librqbr)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
## Section definitions
#==========================================================================================#

"""
HalfWingSection(; span, dihedral, sweep, taper, root_chord,
                  root_twist, tip_twist, root_foil, tip_foil,
                  position, angle, axis)

Define a `HalfWing` consisting of a single trapezoidal section.

# Arguments
- `span       :: Real         = 1.`: Span length
- `dihedral   :: Real         = 1.`: Dihedral angle (degrees)
- `sweep      :: Real         = 0.`: Leading-edge sweep angle (degrees)
- `taper      :: Real         = 1.`: Taper ratio of tip to root chord
- `root_chord :: Real         = 1.`: Root chord length
- `root_twist :: Real         = 0.`: Twist angle at root (degrees)
- `tip_twist  :: Real         = 0.`: Twist angle at tip (degrees)
- `root_foil  :: Foil         = naca4((0,0,1,2))`: Foil coordinates at root
- `tip_foil   :: Foil         = root_foil`: Foil coordinates at tip
- `position   :: Vector{Real} = zeros(3)`: Position
- `angle      :: Real         = 0.`: Angle of rotation (degrees)
- `axis       :: Vector{Real} = [0.,1.,0.]`: Axis of rotation
"""
function HalfWingSection(;
        span         = 1.,
        dihedral     = 0.,
        sweep        = 0.,
        w_sweep      = 0.,
        taper        = 1.,
        root_chord   = 1.,
        root_twist   = 0.,
        tip_twist    = 0.,
        root_control = (0., 0.75),
        tip_control  = root_control,
        root_foil    = naca4(0,0,1,2),
        tip_foil     = root_foil,
        position     = zeros(3),
        angle        = 0.,
        axis         = SVector(0., 1., 0.),
    )
    # Control surface deflections at root and tip
    root_foil = control_surface(root_foil;
                                angle = root_control[1],
                                hinge = root_control[2])

    tip_foil  = control_surface(tip_foil;
                                angle = tip_control[1],
                                hinge = tip_control[2])

    section = HalfWing(
        foils     = [root_foil, tip_foil],
        chords    = [root_chord, taper * root_chord],
        twists    = [root_twist, tip_twist],
        spans     = [span],
        dihedrals = [dihedral],
        sweeps    = [sweep],
        w_sweep   = w_sweep,
        position  = position,
        angle     = angle,
        axis      = axis
    )

    return section
end

function HalfWingSection(S, AR, λ, Λ, δ, τ_r, τ_t, w = 0.25)
    b    = S^2 / AR
    c_r  = 2 * S / (b * (1 + λ))

    HalfWingSection(span = AR, dihedral = δ, sweep = Λ, w_sweep = w, taper = λ, root_chord = c_r, root_twist = τ_r, tip_twist = τ_t)
end

"""
    WingSection(; span, dihedral, sweep, taper, root_chord,
                  root_twist, tip_twist, root_foil, tip_foil,
                  position, angle, axis)

Define a `Wing` consisting of two trapezoidal sections with reflection symmetry in the ``x``-``z`` plane.

# Arguments
- `span       :: Real         = 1.`: Span length
- `dihedral   :: Real         = 1.`: Dihedral angle (degrees)
- `LE_sweep   :: Real         = 0.`: Leading-edge sweep angle (degrees)
- `taper      :: Real         = 1.`: Taper ratio of tip to root chord
- `root_chord :: Real         = 1.`: Root chord length
- `root_twist :: Real         = 0.`: Twist angle at root (degrees)
- `tip_twist  :: Real         = 0.`: Twist angle at tip (degrees)
- `root_foil  :: Foil         = naca4((0,0,1,2))`: Foil coordinates at root
- `tip_foil   :: Foil         = root_foil`: Foil coordinates at tip
- `position   :: Vector{Real} = zeros(3)`: Position
- `angle      :: Real         = 0.`: Angle of rotation (degrees)
- `axis       :: Vector{Real} = [0.,1.,0.]`: Axis of rotation
"""
function WingSection(;
        span         = 1.,
        dihedral     = 0.,
        sweep        = 0.,
        w_sweep      = 0.,
        taper        = 1.,
        root_chord   = 1.,
        root_twist   = 0.,
        tip_twist    = 0.,
        root_control = (0., 0.75),
        tip_control  = root_control,
        root_foil    = naca4((0,0,1,2)),
        tip_foil     = root_foil,
        position     = zeros(3),
        angle        = 0.,
        axis         = SVector(1., 0., 0.),
    )

    right = HalfWingSection(
                span         = span / 2,
                dihedral     = dihedral,
                sweep        = sweep,
                w_sweep      = w_sweep,
                taper        = taper,
                root_chord   = root_chord,
                root_twist   = root_twist,
                tip_twist    = tip_twist,
                root_control = root_control,
                tip_control  = tip_control,
                root_foil    = root_foil,
                tip_foil     = tip_foil,
                position     = position,
                angle        = angle,
                axis         = axis
            )

    Wing(right)
end

## Standard stabilizers (non-exhaustive)
#==========================================================================================#

HorizontalTail(; root_chord = 1.0, taper = 1.0, span = 1.0, sweep = 0., w_sweep = 0., root_twist = 0., tip_twist = 0., root_foil = naca4(0,0,1,2), tip_foil = root_foil, angle = 0.) = WingSection(root_chord = root_chord, taper = taper, span = span, sweep = sweep, w_sweep = w_sweep, root_twist = root_twist, tip_twist = tip_twist, root_foil = naca4(0,0,1,2), tip_foil = tip_foil, angle = angle, axis = [0., 1., 0.])

VerticalTail(; root_chord = 1.0, taper = 1.0, span = 1.0, sweep = 0., w_sweep = 0., root_foil = naca4(0,0,1,2), tip_foil = root_foil) = HalfWingSection(root_chord = root_chord, taper = taper, span = span, sweep = sweep, w_sweep = w_sweep, tip_foil = tip_foil, angle = 90., axis = [1., 0., 0.])


text:
using ITensors.HDF5

function save_structs(struc, path::String)
    function Name(arg)
        string(arg)
    end
    fnames = fieldnames(typeof(struc))
    for fn in fnames
        n = Name(fn)
        d = getfield(struc, fn)

        ## DEBUGGING ##
        if n=="optimized_basis"#typeof(d) == Vector{ITensor}
            try
                h5open(path, "r+") do file
                    for i in 1:length(d)
                        write(file, n*"_"*string(i), d[i])
                    end
                end
            catch e
                if !isfile(path)
                    h5open(path, "w") do file
                        for i in 1:length(d)
                            write(file, n*"_"*string(i), d[i])
                        end
                    end
                else
                    h5open(path, "r+") do file
                        for i in 1:length(d)
                            if haskey(file, n*"_"*string(i))
                                delete_object(file, n*"_"*string(i))
                                file[n*"_"*string(i)] = d[i]
                            end
                        end
                    end
                end
            end
        else
            try
                h5open(path, "r+") do file
                    write(file, n, d)
                end
            catch
                if !isfile(path)
                    h5open(path, "w") do file
                        write(file, n, d)
                    end
                else
                    h5open(path, "r+") do file
                        if haskey(file, n)
                            delete_object(file, n)
                            file[n] = d
                        end
                    end
                end
            end
        end
    end
end

parse_numbers(s, delimiter) = parse(Float64, split(s, delimiter, keepempty=false)[end])

function read_ITensor_vector(d, f, prefix, delimiter)
    # Read in the gates
    Ts = []
    idx = []
    for key in keys(d)
        if startswith(key, prefix)
            T = read(f, key, ITensor)
            push!(Ts, T)
            push!(idx, parse_numbers(key, delimiter))
        end
    end
    return Ts[sortperm(idx)]
end

function load_structs(loadpath::String)
    f = h5open(loadpath,"r")
    d = read(f)

    # Read in the states
    d["ground_state"] = read(f, "ground_state", MPS)
    d["phi_t"] = read(f, "phi_t", MPS)
    d["psi_t"] = read(f, "psi_t", MPS)
    d["mpo"] = read(f, "mpo", MPO)

    # Read in the ITensor vectors
    if d["DMRG_LBO"]
        optimized_basis = read_ITensor_vector(d, f, "optimized_basis_", "_")
    else
        optimized_basis = nothing
    end

    close(f)

    # Make the structs
    sites = siteinds(d["ground_state"])
    params = Parameters(d["N"], d["t"], d["U"], d["ω"], d["g0"], d["g1"],d["λ"],
                        d["doping"],d["max_phonons"], d["init_phonons"],
                        d["DMRG_numsweeps"], d["DMRG_noise"],d["DMRG_maxdim"],
                        d["DMRG_cutoff"], d["DMRG_LBO"], d["max_LBO_dim"],d["min_LBO_dim"],
                        d["mid"], d["T"], d["τ"],
                        d["TEBD_cutoff"], d["TEBD_maxdim"], d["TEBD_LBO"])
    hubbholst = HubbardHolsteinModel(params, sites)
    dmrg_results = DMRGResults(d["ground_state"], d["ground_state_energy"],
                            d["ground_state_entropy"], optimized_basis)
    equilibrium_corr = EquilibriumCorrelations(d["spin"], d["charge"],
                                                d["sSC"], d["pSC"], d["dSC"])
    tebd_results = TEBDResults(d["entropy"], d["self_overlap"],
                                d["corrs"], d["phi_t"], d["psi_t"])

    return params, hubbholst, dmrg_results, equilibrium_corr, tebd_results
end


text:
include("/Users/kurama/.julia/dev/SeisNetwork/src/SeisNetwork.jl")
using .SeisNetwork
#---Input parameters---#
finame = "./src/input_stationlist.jld2"
fogmtname = "./output/station_triangluration.jld2"
fodvvname = "./output/station_dvvtriangluration.jld2"
maxdistance = 60e3 #[m] thresholding triangle mesh with this distance
mindistance = 1e3 #[m] thresholding triangle mesh with this distance
maxangle = 160.0 #[deg] thresholding triangle mesh with this angle
minangle = 20.0 #[deg] thresholding triangle mesh with this angle
studyarea = [37, -119.5, 35, -122]
plotfig = true
figname = "stationloc.png"
#-----------------------#

InputDict = Dict("finame" => finame,
                 "fogmtname" => fogmtname,
                 "fodvvname" => fodvvname,
                 "maxdistance" => maxdistance,
                 "mindistance" => mindistance,
                 "maxangle" => maxangle,
                 "minangle" => minangle,
                 "studyarea" => studyarea,
                 "plotfig"  =>plotfig,
                 "figname"  =>figname
                 )
triangulation(InputDict)


text:
using PseudospectraView
using QML
using Pseudospectra
const PSA=Pseudospectra
using Test

@testset "Basic" begin
    A = PSA.landau_fox_li(200,16)

    # set the time bomb
    PseudospectraView._doomed[] = true

    psagui(A; backend=:gr)

    app_data = PseudospectraView.PSApp.ps_data
    @test size(app_data.ps_dict[:schur_mtx]) == size(A)
end


text:
module FunctionalStateMachine

using Dates
using Requires
using DocStringExtensions
using ProgressMeter

export
  StateMachine,
  emptyState,
  exitStateMachine,
  getIterCount,
  sandboxStateMachineStep,
  getStateLabel,
  histStateMachineTransitions

include("StateMachine.jl")

function __init__()
  @require Graphs="86223c79-3864-5bf0-83f7-82e725a168b6" include("StateMachineAnimation.jl")
end

end


text:
using Packmol
using Test

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


text:
# Autogenerated wrapper script for SwiftShader_jll for x86_64-apple-darwin
export libvulkan

JLLWrappers.@generate_wrapper_header("SwiftShader")
JLLWrappers.@declare_library_product(libvulkan, "@rpath/libvk_swiftshader.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libvulkan,
        "lib/libvulkan.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text: # manova.jl


text:
#!/usr/bin/env julia

using Luxor, Test, Random

Random.seed!(42)

w, h = 600, 600

fname = "polycross-test.png"
Drawing(w, h, fname)
background("white")
origin()
setline(0.5)
for i in 100:-1:1
    rotate(0.01)
    polycross(O, 250, 4, i/100, 0, :stroke)
end

@test finish() == true

println("...finished test: output in $(fname)")


text:
"""
### function plot_nyquist(real_fr, imag_fr, num_u, num_y)

Plot the Nyquist diagram.

##### Args

* real_fr: Real part of the frequency response.
* imag_fr: Imaginary part of the frequency response.
* num_u: Number of inputs.
* num_y: Number of outputs.

"""

function plot_nyquist(real_fr, imag_fr, num_u, num_y)
    if (plot_system == :pyplot)
        # Get the number of graphics.
        num_plots = length(real_fr[1])

        fig, ax = plt.subplots(num_y, num_u, sharex=true, sharey=true)

        ( num_plots == 1 ) && (ax = [ax])

        for i = 1:num_plots
            real_fr_i = map(x->x[i], real_fr)
            imag_fr_i = map(x->x[i], imag_fr)

            # Get the default color for the plots.
            color = pybuiltin("next")(ax[i][:_get_lines][:prop_cycler])["color"]

            ax[i][:plot](real_fr_i, +imag_fr_i, color=color, linewidth=plw)
            ax[i][:plot](real_fr_i, -imag_fr_i, color=color, linewidth=plw)

            ax[i][:grid]("on", which="major",
                         color="#AAAAAA", linewidth=glw, zorder=1)

            # Mark the -1 point but do not resize the plot.
            ax[i][:autoscale](false)
            ax[i][:scatter]([-1.0],[0.0],
                            s=mksz, marker="+", linewidths=plw,
                            edgecolor="red", color="red", zorder=3)

            # Create arrows at:
            #     1) 1/4 of the range of the real part;
            #     2) 1/2 of the range of the real part;
            #     3) 3/4 of the range of the real part.
            min_real_fr = minimum(real_fr_i[isfinite.(real_fr_i)])
            max_real_fr = maximum(real_fr_i[isfinite.(real_fr_i)])
            Δp_real     = max_real_fr - min_real_fr

            fq_x  = min_real_fr + 1*Δp_real/4
            mid_x = min_real_fr + 2*Δp_real/4
            tq_x  = min_real_fr + 3*Δp_real/4

            for p in [fq_x, mid_x, tq_x]
                # Search for the freq. that provides the closest real part to p.
                ind = indmin(abs.(real_fr_i-p))

                # Compute the tangent.
                Δ = [real_fr_i[ind+1] - real_fr_i[ind];
                     imag_fr_i[ind+1] - imag_fr_i[ind];]

                ax[i][:annotate]("",
                                 xytext = [real_fr_i[ind];
                                           imag_fr_i[ind]],
                                 xy     = [real_fr_i[ind]+Δ[1]*0.01;
                                           imag_fr_i[ind]+Δ[2]*0.01],
                                 arrowprops = Dict("arrowstyle"=>"simple"),
                                 size       = arrsz)

                ax[i][:annotate]("",
                                 xytext = [ real_fr_i[ind+1];
                                           -imag_fr_i[ind+1]],
                                 xy     = [ real_fr_i[ind+1]-Δ[1]*0.01;
                                           -imag_fr_i[ind+1]+Δ[2]*0.01],
                                 arrowprops = Dict("arrowstyle"=>"simple"),
                                 size       = arrsz)
            end

            # Change figure border width.
            ax[i][:spines]["top"][:set_linewidth](blw)
            ax[i][:spines]["bottom"][:set_linewidth](blw)
            ax[i][:spines]["left"][:set_linewidth](blw)
            ax[i][:spines]["right"][:set_linewidth](blw)
        end

        for i = 1:num_u
            (num_u > 1) && ax[num_y*(i-1)+1][:set_title]("From Input $i")
            ax[num_y*i][:set_xlabel]("Real")
        end

        for i = 1:num_y
            ax[i][:set_ylabel]("Imag.")

            ax_aux = ax[(num_u-1)*num_y+(i-1)+1][:twinx]()
            ax_aux[:yaxis][:set_label_position]("right")
            ax_aux[:set_yticks]([])
            ax_aux[:set_ylabel]("To Output $i")
        end

        fig[:suptitle]("Nyquist Diagram", fontsize=tifsz)
    end
end


text:
function economicdispatch(system::MarketClearingProblem{Z, T},
        solver::JuMP.MOI.OptimizerWithAttributes,
        resultfile::String="") where {Z, T}

        lines = [line.name for line in system.lines] # Lines
        projects = [project.name for project in system.projects] # Projects
        invperiods = 1:length(system.inv_periods)     # Investment periods
        opperiods = 1:T
        rep_hour_weight = system.rep_hour_weight

        end_of_day = collect(24:24:T)

        project_type,                             # is the project of storage type
        marginal_energy_cost,                      # $/MW
        marginal_reserveup_cost,                # $/MW
        marginal_reservedown_cost,                # $/MW
        min_gen, max_gen,                           # MW/unit
        min_input, max_input,                       # MW/unit
        efficiency_in, efficiency_out,            # percentage
        min_storage, max_storage,                   # MWh/unit
        init_storage,                              # MWh/unit
        availability,                             # hourly availability factor
        ramp_limits,                               # MW/unit/hour
        max_reserveup, max_reservedown,             # MW/unit/hour
        location =                                # Project zone location
        make_parameter_vector.(
        Ref(system.projects), :name,
        [:project_type, :marginal_energy_cost, :marginal_reserveup_cost, :marginal_reservedown_cost,
        :min_gen, :max_gen, :min_input, :max_input,
        :efficiency_in, :efficiency_out, :min_storage, :max_storage, :init_storage,
        :availability,
        :ramp_limits, :max_reserveup, :max_reservedown, :zone]
        )

        zones = system.zones
        price_cap_e = AxisArrays.AxisArray(zeros(length(zones), length(invperiods)), zones, invperiods)
        demand_e = AxisArrays.AxisArray(zeros(length(zones), length(invperiods), T), zones, invperiods, opperiods)

        demand_ru = AxisArrays.AxisArray(zeros(length(zones), length(invperiods), T), zones, invperiods, opperiods)

        price_cap_rd = AxisArrays.AxisArray(zeros(length(zones), length(invperiods)), zones, invperiods)
        demand_rd = AxisArrays.AxisArray(zeros(length(zones), length(invperiods), T), zones, invperiods, opperiods)

        for p in invperiods
                for z in zones
                        price_cap_e[z, p] = getproperty(getproperty(system.inv_periods[p], :energy_market), :price_cap)[z]
                        demand_e[z, p, :] = getproperty(getproperty(system.inv_periods[p], :energy_market), :demand)[z, :]

                        demand_ru[z, p, :] = getproperty(getproperty(system.inv_periods[p], :reserveup_market), :demand)[z, :]

                        price_cap_rd[z, p] = getproperty(getproperty(system.inv_periods[p], :reservedown_market), :price_cap)[z]
                        demand_rd[z, p, :] = getproperty(getproperty(system.inv_periods[p], :reservedown_market), :demand)[z, :]
                end
        end

        ru_segmentsize, ru_segmentgrad, ru_price_points, ru_numsegments = make_ORDC_vectors(getproperty.(system.inv_periods, :reserveup_market))

        zone_projects = Dict(z => Vector{String}() for z in zones)

        generator_projects = Vector{String}()
        storage_projects = Vector{String}()

        zone_projects = Dict(z => Vector{String}() for z in zones)
        zone_storage = Dict(z => Vector{String}() for z in zones)

        ramp_lim_projects = Vector{String}()

        for project in projects

                if ramp_limits[project] !== nothing
                        push!(ramp_lim_projects, project)
                end

            for z in zones
                if location[project] == z
                    push!(zone_projects[z], project)

                    if project_type[project] == "generator"
                        push!(generator_projects, project)
                    elseif project_type[project] == "storage"
                        push!(storage_projects, project)
                        push!(zone_storage[z], project)
                    end

                end
            end
        end

        # Populate line data
        from_zone,
        to_zone,
        linepowerlimit =
        make_parameter_vector.(
        Ref(system.lines), :name,
        [:from_zone, :to_zone, :active_power_limit]
        )

        lines_from_zone = Dict(z => Vector{String}() for z in zones)
        lines_to_zone = Dict(z => Vector{String}() for z in zones)

        for line in lines
                for zone in zones
                        if from_zone[line] == zone
                                push!(lines_from_zone[zone], line)
                        elseif to_zone[line] == zone
                                push!(lines_to_zone[zone], line)
                        end
                end
        end

        m = JuMP.Model(solver)

        # OPERATION DECISION VARIABLES

        JuMP.@variable(m, p_e[g in projects, p in invperiods, t in opperiods] >= 0) # Unit energy production [MW]
        JuMP.@variable(m, p_in[g in storage_projects, p in invperiods, t in opperiods] >= 0) # Storage charging [MW]
        JuMP.@variable(m, p_ru[g in projects, p in invperiods, t in opperiods] >= 0) # reserve up provided [MW]
        JuMP.@variable(m, p_rd[g in projects, p in invperiods, t in opperiods] >= 0) # reserve down provided [MW]

        JuMP.@variable(m, p_in_ru[g in storage_projects, p in invperiods, t in opperiods] >= 0) # reserve up provided by storage charging [MW]
        JuMP.@variable(m, p_in_rd[g in storage_projects, p in invperiods, t in opperiods] >= 0) # reserve down provided by storage charging [MW]

        JuMP.@variable(m, p_out_ru[g in storage_projects, p in invperiods, t in opperiods] >= 0) # reserve up provided by storage discharging [MW]
        JuMP.@variable(m, p_out_rd[g in storage_projects, p in invperiods, t in opperiods] >= 0) # reserve down provided by storage discharging [MW]

        JuMP.@variable(m, flow[l in lines, p in invperiods, t in opperiods])     # Line flow

        # ORDC DEMAND VARIABLE
        JuMP.@variable(m, d_ru[z in zones, p in invperiods, t in opperiods, s in 1:ru_numsegments[z, p]] >= 0)

        JuMP.@variable(m, v_e[z in zones, p in invperiods, t in opperiods] >= 0) # Load shortfall [MW]
        JuMP.@variable(m, v_rd[z in zones, p in invperiods, t in opperiods] >= 0) # Reserve Down shortfall [MW]

        JuMP.@expression(m, annual_project_costs[g in projects, p in invperiods], 0
                              + sum((p_e[g, p, t] * marginal_energy_cost[g]
                                  + p_ru[g, p, t] * marginal_reserveup_cost[g]
                                  + p_rd[g, p, t] * marginal_reservedown_cost[g]) * rep_hour_weight[t]  for t in opperiods))    # project operating costs

        # Storage level evolution
        JuMP.@expression(m, storage_level[g in storage_projects, p in invperiods, t in opperiods],
                    p_in[g, p, t] * efficiency_in[g] - p_e[g, p, t] / efficiency_out[g])

        for g in storage_projects
                for p in invperiods
                        for t in opperiods
                                if t == 1 || in(t, end_of_day + ones(length(end_of_day)))
                                        JuMP.add_to_expression!(storage_level[g, p, t], init_storage[g] * max_storage[g])
                                else
                                        JuMP.add_to_expression!(storage_level[g, p, t], storage_level[g, p, t - 1])
                                end
                        end
                end

        end

        # OBJECTIVE FUNCTION

        JuMP.@objective(m, Max,

        # ORDC welfare
        sum(sum(d_ru[z, p, t, s] * demand_ru[z, p, t] * ru_price_points[z, p][s] * rep_hour_weight[t] for z in zones, t in opperiods, s in 1:ru_numsegments[z, p])
                for p in invperiods)
        + 1/2 * sum(sum(ru_segmentgrad[z, p][s] * (d_ru[z, p, t, s] * demand_ru[z, p, t])^2 * rep_hour_weight[t] for z in zones, t in opperiods, s in 1:ru_numsegments[z, p])
                for p in invperiods)

        # Annual costs incurred by the projects
        - sum(annual_project_costs[g, p]
            for g in projects, p in invperiods)

        # weight operations according to actual duration
        # operations market shortfall penalties
        + ( - sum((v_e[z, p, t] * price_cap_e[z, p] + v_rd[z, p, t] * price_cap_rd[z, p]) * rep_hour_weight[t]
                for z in zones, p in invperiods, t in opperiods)
         )

            )

    # CONSTRAINTS

    # Operations:

        # Energy market
        JuMP.@constraint(m, energy_market[z in zones, p in invperiods, t in opperiods], # Power balance
                sum(p_e[g ,p, t] for g in zone_projects[z])
                + sum(flow[l, p, t] for l in lines_to_zone[z])
                - sum(flow[l, p, t] for l in lines_from_zone[z])
                + v_e[z, p, t] == demand_e[z, p, t] + sum(p_in[g, p, t] for g in zone_storage[z]))

        # reserve up market
        JuMP.@constraint(m, reserve_up_market[z in zones, p in invperiods, t in opperiods], # reserve up balance
                sum(p_ru[g ,p, t] for g in zone_projects[z]) == sum(d_ru[z, p, t, s] for s in 1:ru_numsegments[z, p]) * demand_ru[z, p, t] )

        JuMP.@constraint(m, ru_segment_limit[z in zones, p in invperiods, t in opperiods, s in 1:ru_numsegments[z, p]],
                d_ru[z, p, t, s] <= ru_segmentsize[z, p][s])

        # reserve down market
        JuMP.@constraint(m, reserve_down_market[z in zones, p in invperiods, t in opperiods], # reserve down balance
                sum(p_rd[g ,p, t] for g in zone_projects[z]) + v_rd[z, p, t]  == demand_rd[z, p, t])

        # Line flow constraints
        JuMP.@constraint(m, negative_flow[l in lines, p in invperiods, t in opperiods],
                        flow[l, p, t] >= -linepowerlimit[l])

        JuMP.@constraint(m, positive_flow[l in lines, p in invperiods, t in opperiods],
                        flow[l, p, t] <= linepowerlimit[l])

        # Unit operating constraints (dispatch and ramping)

        # Generator technical constraints
        JuMP.@constraint(m, min_generation[g in generator_projects, p in invperiods, t in opperiods],
                   p_e[g, p, t] - p_rd[g, p, t] >= 0.0) # Minimum dispatch

        JuMP.@constraint(m, max_generation[g in generator_projects, p in invperiods, t in opperiods],
                   p_e[g, p, t] + p_ru[g, p, t] <= max_gen[g] * availability[g][p, t]) # Maximum dispatch

        JuMP.@constraint(m, reserve_down_max[g in generator_projects, p in invperiods, t in opperiods],
                   p_rd[g, p, t] <= max_reservedown[g]) # Maximum Reserve Down

        JuMP.@constraint(m, reserve_up_max[g in generator_projects, p in invperiods, t in opperiods],
                   p_ru[g, p, t] <= max_reserveup[g]) # Maximum Reserve Up

        JuMP.@constraint(m, rampdown[g in ramp_lim_projects, p in invperiods, t in opperiods[2:end]], # Down-ramp limit
                      p_e[g, p ,t] - p_rd[g, p ,t] >= p_e[g, p, t-1] - ramp_limits[g][:down])

        JuMP.@constraint(m, rampup[g in ramp_lim_projects, p in invperiods, t in opperiods[2:end]], # Up-ramp limit
                          p_e[g, p, t] + p_ru[g, p, t] <= p_e[g, p, t-1] + ramp_limits[g][:up])

        # Storage technical constraints
        JuMP.@constraint(m, minoutput[g in storage_projects, p in invperiods, t in opperiods],
                   p_e[g, p, t] - p_out_rd[g, p, t] >= 0.0) # Minimum output dispatch

        JuMP.@constraint(m, maxoutput[g in storage_projects, p in invperiods, t in opperiods],
                   p_e[g, p, t] + p_out_ru[g, p, t] <= max_gen[g]) # Maximum output dispatch

        JuMP.@constraint(m, min_input[g in storage_projects, p in invperiods, t in opperiods],
                   p_in[g, p, t] - p_in_ru[g, p, t] >= 0.0) # Minimum input dispatch

        JuMP.@constraint(m, max_input[g in storage_projects, p in invperiods, t in opperiods],
                   p_in[g, p, t] + p_in_rd[g, p, t] <= max_input[g]) # Maximum input dispatch

        JuMP.@constraint(m, storage_reserve_up[g in storage_projects, p in invperiods, t in opperiods],
                   p_ru[g, p, t] <= p_out_ru[g, p, t] + p_in_ru[g, p, t]) # Maximum storage reserve up

        JuMP.@constraint(m, storage_reserve_down[g in storage_projects, p in invperiods, t in opperiods],
                   p_rd[g, p, t] <= p_out_rd[g, p, t] + p_in_rd[g, p, t]) # Maximum storage reserve down

        JuMP.@constraint(m, min_storage_level[g in storage_projects, p in invperiods, t in opperiods],
                   storage_level[g, p, t] >= min_storage[g]) # Minimum storage level

        JuMP.@constraint(m, max_storage_level[g in storage_projects, p in invperiods, t in opperiods],
                   storage_level[g, p, t] <= max_storage[g]) # Maximum storage level

        JuMP.@constraint(m, end_day_storage_level[g in storage_projects, p in invperiods, t in end_of_day],
                   storage_level[g, p, t] == init_storage[g] * max_storage[g]) # Maximum storage level after reserve provision

        println("Actual Energy and A/S Market Clearing:")
        JuMP.optimize!(m)
        println(JuMP.termination_status(m))
        println(JuMP.objective_value(m))

        energy_dual = JuMP.dual.(energy_market)
        reserve_up_dual = JuMP.dual.(reserve_up_market)
        reserve_down_dual = JuMP.dual.(reserve_down_market)

        energy_price = AxisArrays.AxisArray(zeros(size(energy_dual)), zones, invperiods, opperiods)
        reserve_up_price = AxisArrays.AxisArray(zeros(size(reserve_up_dual)), zones, invperiods, opperiods)
        reserve_down_price = AxisArrays.AxisArray(zeros(size(reserve_down_dual)), zones, invperiods, opperiods)

        for p in invperiods
                for t in opperiods
                        energy_price[:, p, t] = round.((energy_dual[:, p, t] / rep_hour_weight[t]), digits = 5)
                        reserve_up_price[:, p, t] = round.((reserve_up_dual[:, p, t] / rep_hour_weight[t]), digits = 5)
                        reserve_down_price[:, p, t] = round.((reserve_down_dual[:, p, t] / rep_hour_weight[t]), digits = 5)
                end
        end

        capacity_factor = Dict([g => zeros(length(invperiods), length(opperiods)) for g in projects])
        reserve_up_perc = Dict([g => zeros(length(invperiods), length(opperiods)) for g in projects])
        reserve_down_perc = Dict([g => zeros(length(invperiods), length(opperiods)) for g in projects])

        for g in projects
                for p in invperiods
                        for t in opperiods
                                capacity_factor[g][p, t] = value.(p_e[g, p, t]) / max_gen[g]
                                reserve_up_perc[g][p, t] = value.(p_ru[g, p, t]) / max_gen[g]
                                reserve_down_perc[g][p, t] = value.(p_rd[g, p, t]) / max_gen[g]
                        end

                end
        end

        return energy_price, reserve_up_price, reserve_down_price, capacity_factor, reserve_up_perc, reserve_down_perc;

    end


text:
Base.:*(u::AbstractVector, v::SVector{1,T}) where T = u * v[1]

Base.vcat(v::Adjoint{T,SVector{N,T}}, E::UniformScaling) where {T, N} = vcat(v, SMatrix{N,N,T,N*N}(E))
Base.vcat(A::SMatrix{N1,N2,T,N1N2}, E::UniformScaling) where {T, N1, N2, N1N2} = vcat(A, SMatrix{N2,N2,T,N2*N2}(E))

svcat(a::T, b::T) where T = SA[a; b]
svcat(a::StaticArray, b::StaticArray) = vcat(a, b)
svcat(a::StaticArray{Tuple{N},T,1}, b::T) where {T,N} = vcat(a, SA[b])
svcat(a::T, b::StaticArray{Tuple{N},T,1}) where {T,N} = vcat(SA[a], b)

szeros(::Type{T}, N) where T = @SVector zeros(T, N)
szeros(N)= @SVector zeros(N)
szeros(::Type{T}, N1, N2) where T = @SMatrix zeros(T, N1, N2)
szeros(N1, N2)= @SMatrix zeros(N1, N2)

sones(::Type{T}, N) where T = @SVector ones(T, N)
sones(N)= @SVector ones(N)

srand(N)= @SVector rand(N)

# TODO: check StaticArray bug fix, then remove
zerodimstaticadjoint(A) = A'
zerodimstaticadjoint(::SMatrix{0,N,T,0}) where {T,N} = SMatrix{N,0,T,0}()


text:
"""
    $(SIGNATURES)
Gets an empty and unique CloudGraphsDFG derived from an existing DFG.
"""
function _getDuplicatedEmptyDFG(dfg::CloudGraphsDFG)::CloudGraphsDFG
    count = 0
    sessionId = dfg.sessionId*"_$count"
    while true #do..while loop
        count += 1
        sessionId = dfg.sessionId*"_$count"
        _getNodeCount(dfg.neo4jInstance, [dfg.userId, dfg.robotId, sessionId]) == 0 && break
    end
    @debug "Unique+empty copy session name: $sessionId"
    return CloudGraphsDFG{typeof(dfg.solverParams)}(
        dfg.neo4jInstance.connection,
        dfg.userId,
        dfg.robotId,
        sessionId,
        dfg.description,
        solverParams=deepcopy(dfg.solverParams))
end

# Accessors
function getDescription(dfg::CloudGraphsDFG)
    return _getNodeProperty(
        dfg.neo4jInstance,
        [dfg.sessionId, dfg.robotId, dfg.userId, "SESSION"],
        "description")
end
function setDescription!(dfg::CloudGraphsDFG, description::String)
    count = _setNodeProperty(dfg.neo4jInstance,
                            [dfg.sessionId, dfg.robotId, dfg.userId, "SESSION"],
                            "description",
                            description)
    @assert(count == 1)
    dfg.description = description
    return getDescription(dfg)
end

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

function getUserData(dfg::CloudGraphsDFG)::Dict{Symbol, String}
    propVal = _getNodeProperty(dfg.neo4jInstance, _getLabelsForType(dfg, User), "data")
    return JSON2.read(String(base64decode(propVal)), Dict{Symbol, String})
end
function setUserData!(dfg::CloudGraphsDFG, data::Dict{Symbol, String})::Dict{Symbol, String}
    count = _setNodeProperty(dfg.neo4jInstance, [dfg.userId, "USER"], "data", base64encode(JSON2.write(data)))
    @assert(count == 1)
    return  getUserData(dfg)
end

function getRobotData(dfg::CloudGraphsDFG)::Dict{Symbol, String}
    propVal = _getNodeProperty(dfg.neo4jInstance, [dfg.userId, dfg.robotId, "ROBOT"], "data")
    return JSON2.read(String(base64decode(propVal)), Dict{Symbol, String})
end
function setRobotData!(dfg::CloudGraphsDFG, data::Dict{Symbol, String})::Dict{Symbol, String}
    count = _setNodeProperty(dfg.neo4jInstance, [dfg.userId, dfg.robotId, "ROBOT"], "data", base64encode(JSON2.write(data)))
    @assert(count == 1)
    return  getRobotData(dfg)
end

function getSessionData(dfg::CloudGraphsDFG)::Dict{Symbol, String}
    propVal = _getNodeProperty(dfg.neo4jInstance, [dfg.userId, dfg.robotId, dfg.sessionId, "SESSION"], "data")
    return JSON2.read(String(base64decode(propVal)), Dict{Symbol, String})
end
function setSessionData!(dfg::CloudGraphsDFG, data::Dict{Symbol, String})::Dict{Symbol, String}
    count = _setNodeProperty(dfg.neo4jInstance, [dfg.userId, dfg.robotId, dfg.sessionId, "SESSION"], "data", base64encode(JSON2.write(data)))
    @assert(count == 1)
    return  getSessionData(dfg)
end

# New API
getUserData(dfg::CloudGraphsDFG, key::Symbol) = getUserData(dfg::CloudGraphsDFG)[key]
getRobotData(dfg::CloudGraphsDFG, key::Symbol)::String = getRobotData(dfg::CloudGraphsDFG)[key]
getSessionData(dfg::CloudGraphsDFG, key::Symbol)::String = getSessionData(dfg::CloudGraphsDFG)[key]

function updateUserData!(dfg::CloudGraphsDFG, pair::Pair{Symbol,String})
    data = getUserData(dfg::CloudGraphsDFG)
    push!(data, pair)
    setUserData!(dfg, data)
end

function updateRobotData!(dfg::CloudGraphsDFG, pair::Pair{Symbol,String})
    data = getRobotData(dfg::CloudGraphsDFG)
    push!(data, pair)
    setRobotData!(dfg, data)
end

function updateSessionData!(dfg::CloudGraphsDFG, pair::Pair{Symbol,String})
    data = getSessionData(dfg::CloudGraphsDFG)
    push!(data, pair)
    setSessionData!(dfg, data)
end

function deleteUserData!(dfg::CloudGraphsDFG, key::Symbol)
    data = getUserData(dfg::CloudGraphsDFG)
    delval = pop!(data, key)
    setUserData!(dfg, data)
    return delval
end

function deleteRobotData!(dfg::CloudGraphsDFG, key::Symbol)
    data = getRobotData(dfg::CloudGraphsDFG)
    delval = pop!(data, key)
    setRobotData!(dfg, data)
    return delval
end

function deleteSessionData!(dfg::CloudGraphsDFG, key::Symbol)
    data = getSessionData(dfg::CloudGraphsDFG)
    delval = pop!(data, key)
    setSessionData!(dfg, data)
    return delval
end

function emptyUserData!(dfg::CloudGraphsDFG)
    return setUserData!(dfg, Dict{Symbol, String}())
end

function emptyRobotData!(dfg::CloudGraphsDFG)
    return setRobotData!(dfg, Dict{Symbol, String}())
end

function emptySessionData!(dfg::CloudGraphsDFG)
    return setSessionData!(dfg, Dict{Symbol, String}())
end

##==============================================================================
## CRUD Interfaces
##==============================================================================
##------------------------------------------------------------------------------
## Variable And Factor CRUD
##------------------------------------------------------------------------------

function exists(dfg::CloudGraphsDFG, label::Symbol)
    # Otherwise try get it
    nodeId = _tryGetNeoNodeIdFromNodeLabel(dfg.neo4jInstance, dfg.userId, dfg.robotId, dfg.sessionId, label)
    nodeId != nothing && return true
    return false
end
function exists(dfg::CloudGraphsDFG, node::N) where N <: DFGNode
    return exists(dfg, node.label)
end

isVariable(dfg::CloudGraphsDFG, sym::Symbol)::Bool =
    _getNodeCount(dfg.neo4jInstance, ["VARIABLE", dfg.userId, dfg.robotId, dfg.sessionId, String(sym)]) == 1

isFactor(dfg::CloudGraphsDFG, sym::Symbol)::Bool =
    _getNodeCount(dfg.neo4jInstance, ["FACTOR", dfg.userId, dfg.robotId, dfg.sessionId, String(sym)]) == 1

function addVariable!(dfg::CloudGraphsDFG, variable::DFGVariable)::DFGVariable
    if exists(dfg, variable)
        error("Variable '$(variable.label)' already exists in the factor graph")
    end
    props = packVariable(dfg, variable)

    neo4jNode = Neo4j.createnode(dfg.neo4jInstance.graph, props);
    Neo4j.updatenodelabels(neo4jNode, union([string(variable.label), "VARIABLE", dfg.userId, dfg.robotId, dfg.sessionId], variable.tags))

    # Make sure that if there exists a SESSION sentinel that it is attached.
    _bindSessionNodeToInitialVariable(dfg.neo4jInstance, dfg.userId, dfg.robotId, dfg.sessionId, string(variable.label))

    # Track insertion
    push!(dfg.addHistory, variable.label)

    return variable
end

function addFactor!(dfg::CloudGraphsDFG, factor::DFGFactor)::DFGFactor
    if exists(dfg, factor)
        error("Factor '$(factor.label)' already exists in the factor graph")
    end

    variableIds = factor._variableOrderSymbols
    variables = map(vId -> getVariable(dfg, vId), variableIds)

    # Construct the properties to save
    props = packFactor(dfg, factor)

    # TODO - Pack this into a transaction.
    neo4jNode = Neo4j.createnode(dfg.neo4jInstance.graph, props);
    Neo4j.updatenodelabels(neo4jNode, union([string(factor.label), "FACTOR", dfg.userId, dfg.robotId, dfg.sessionId], factor.tags))

    # Add all the relationships - get them to cache them + make sure the links are correct
    for variable in variables
        v = getVariable(dfg, variable.label)
        vNode = Neo4j.getnode(
            dfg.neo4jInstance.graph,
            _tryGetNeoNodeIdFromNodeLabel(dfg.neo4jInstance, dfg.userId, dfg.robotId, dfg.sessionId, variable.label))
        Neo4j.createrel(neo4jNode, vNode, "FACTORGRAPH")
    end

    return factor
end

function getVariable(dfg::CloudGraphsDFG, label::Union{Symbol, String})::DFGVariable
    if typeof(label) == String
        label = Symbol(label)
    end
    nodeId = _tryGetNeoNodeIdFromNodeLabel(dfg.neo4jInstance, dfg.userId, dfg.robotId, dfg.sessionId, label, nodeType="VARIABLE")
    if nodeId == nothing
        error("Unable to retrieve the ID for variable '$label'. Please check your connection to the database and that the variable exists.")
    end

    props = getnodeproperties(dfg.neo4jInstance.graph, nodeId)
    variable = unpackVariable(dfg, props)

    # TODO - make this get PPE's in batch
    for ppe in listPPEs(dfg, label)
        variable.ppeDict[ppe] = getPPE(dfg, label, ppe)
    end
    return variable
end

function getFactor(dfg::CloudGraphsDFG, label::Union{Symbol, String})::DFGFactor
    if typeof(label) == String
        label = Symbol(label)
    end
    nodeId = _tryGetNeoNodeIdFromNodeLabel(dfg.neo4jInstance, dfg.userId, dfg.robotId, dfg.sessionId, label, nodeType="FACTOR")
    if nodeId == nothing
        error("Unable to retrieve the ID for factor '$label'. Please check your connection to the database and that the factor exists.")
    end

    props = getnodeproperties(dfg.neo4jInstance.graph, nodeId)
    factor = unpackFactor(dfg, props)

    # Lastly, rebuild the metadata
    factor = rebuildFactorMetadata!(dfg, factor)

    return factor
end

function updateVariable!(dfg::CloudGraphsDFG, variable::DFGVariable)::DFGVariable
    if !exists(dfg, variable)
        @warn "Variable label '$(variable.label)' does not exist in the factor graph, adding"
        return addVariable!(dfg, variable)
    end
    nodeId = _tryGetNeoNodeIdFromNodeLabel(dfg.neo4jInstance, dfg.userId, dfg.robotId, dfg.sessionId, variable.label)

    # TODO: Do this with a single query and/or a single transaction.
    neo4jNode = Neo4j.getnode(dfg.neo4jInstance.graph, nodeId)
    props = packVariable(dfg, variable)
    Neo4j.updatenodeproperties(neo4jNode, props)
    Neo4j.updatenodelabels(neo4jNode, union([string(variable.label), "VARIABLE", dfg.userId, dfg.robotId, dfg.sessionId], variable.tags))
    return variable
end

function mergeVariableData!(dfg::CloudGraphsDFG, sourceVariable::DFGVariable)::DFGVariable
    if !exists(dfg, sourceVariable)
        error("Source variable '$(sourceVariable.label)' doesn't exist in the graph.")
    end
    for (k,v) in sourceVariable.ppeDict
        # TODO what is happening inside, is this truely an update or is it a merge? (API consistency hounds are apon you)
        updatePPE!(dfg, getLabel(sourceVariable), v, k)
    end
    for (k,v) in sourceVariable.solverDataDict
        updateVariableSolverData!(dfg, getLabel(sourceVariable), v, k)
    end
    return sourceVariable
end

function updateFactor!(dfg::CloudGraphsDFG, factor::DFGFactor)::DFGFactor
    if !exists(dfg, factor)
        @warn "Factor label '$(factor.label)' does not exist in the factor graph, adding"
        return addFactor!(dfg, factor)
    end

    neighborsList = JSON2.read(
        _getNodeProperty(dfg.neo4jInstance, _getLabelsForInst(dfg, factor), "_variableOrderSymbols"),
        Vector{Symbol})

    # Confirm that we're not updating the neighbors
    neighborsList != factor._variableOrderSymbols && error("Cannot update the factor, the neighbors are not the same.")

    # TODO: Optimize this as a single query with a single transaction.
    nodeId = _tryGetNeoNodeIdFromNodeLabel(dfg.neo4jInstance, dfg.userId, dfg.robotId, dfg.sessionId, factor.label)
    neo4jNode = Neo4j.getnode(dfg.neo4jInstance.graph, nodeId)
    props = packFactor(dfg, factor)
    Neo4j.updatenodeproperties(neo4jNode, props)
    Neo4j.updatenodelabels(neo4jNode, union([string(factor.label), "FACTOR", dfg.userId, dfg.robotId, dfg.sessionId], factor.tags))

    return factor
end

function deleteVariable!(dfg::CloudGraphsDFG, label::Symbol)#::Tuple{AbstractDFGVariable, Vector{<:AbstractDFGFactor}}
    variable = getVariable(dfg, label)
    if variable == nothing
        error("Unable to retrieve the ID for variable '$label'. Please check your connection to the database and that the variable exists.")
    end

    deleteNeighbors = true # reserved, orphaned factors are not supported at this time
    if deleteNeighbors
        neigfacs = map(l->deleteFactor!(dfg, l), getNeighbors(dfg, label))
    end

    # Perform detach+deletion
    _getNeoNodesFromCyphonQuery(dfg.neo4jInstance, "(node:$label:$(join(_getLabelsForType(dfg, DFGVariable),':'))) detach delete node ")

    # Clearing history
    dfg.addHistory = symdiff(dfg.addHistory, [label])
    return variable, neigfacs
end

#Alias
deleteVariable!(dfg::CloudGraphsDFG, variable::DFGVariable) = deleteVariable!(dfg, variable.label)

function deleteFactor!(dfg::CloudGraphsDFG, label::Symbol)::DFGFactor
    factor = getFactor(dfg, label)
    if factor == nothing
        error("Unable to retrieve the ID for factor '$label'. Please check your connection to the database and that the factor exists.")
    end

    # Perform detach+deletion
    _getNeoNodesFromCyphonQuery(dfg.neo4jInstance, "(node:$label:$(join(_getLabelsForType(dfg, DFGFactor),':'))) detach delete node ")

    # Clearing history
    dfg.addHistory = symdiff(dfg.addHistory, [label])
    return factor
end

# Alias
deleteFactor!(dfg::CloudGraphsDFG, factor::DFGFactor)::DFGFactor = deleteFactor!(dfg, factor.label)

function getVariables(dfg::CloudGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[], solvable::Int=0)::Vector{DFGVariable}
    variableIds = listVariables(dfg, regexFilter, tags=tags, solvable=solvable)
    # TODO: Optimize to use tags in query here!
    variables = map(vId->getVariable(dfg, vId), variableIds)
    if length(tags) > 0
        mask = map(v -> length(intersect(v.tags, tags)) > 0, variables )
        return variables[mask]
    end
    return variables
end

function listVariables(dfg::CloudGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[], solvable::Int=0)::Vector{Symbol}
    # Optimized for DB call
    tagsFilter = length(tags) > 0 ? " and "*join("node:".*String.(tags), " or ") : ""
    if regexFilter == nothing
        return _getLabelsFromCyphonQuery(dfg.neo4jInstance, "(node:$(join(_getLabelsForType(dfg, DFGVariable),':'))) where node.solvable >= $solvable $tagsFilter")
    else
        return _getLabelsFromCyphonQuery(dfg.neo4jInstance, "(node:$(join(_getLabelsForType(dfg, DFGVariable),':'))) where node.label =~ '$(regexFilter.pattern)' and node.solvable >= $solvable $tagsFilter")
    end
end

function getFactors(dfg::CloudGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing; solvable::Int=0)::Vector{DFGFactor}
    # TODO: Use optimized cypher query.
    factorIds = listFactors(dfg, regexFilter, solvable=solvable)
    return map(vId->getFactor(dfg, vId), factorIds)
end

function listFactors(dfg::CloudGraphsDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[], solvable::Int=0)::Vector{Symbol}
    # Optimized for DB call
    length(tags) > 0 && (@error "Filter on tags not implemented for CloudGraphsDFG")
    if regexFilter == nothing
        return _getLabelsFromCyphonQuery(dfg.neo4jInstance, "(node:$(join(_getLabelsForType(dfg, DFGFactor),':'))) where node.solvable >= $solvable")
    else
        return _getLabelsFromCyphonQuery(dfg.neo4jInstance, "(node:$(join(_getLabelsForType(dfg, DFGFactor),':'))) where node.label =~ '$(regexFilter.pattern)' and node.solvable >= $solvable")
    end
end

function isConnected(dfg::CloudGraphsDFG)::Bool
    # If the total number of nodes == total number of distinct connected nodes, then it is fully connected
    # Total nodes
    varIds = listVariables(dfg)
    factIds = listFactors(dfg)
    length(varIds) + length(factIds) == 0 && return false

    # Total distinct connected nodes - thank you Neo4j for 0..* awesomeness!!
    # TODO: Deprecated matching technique and it's technically an expensive call - optimize.
    query = """
        match (n:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId):$(varIds[1]))-[FACTORGRAPH*]-(node:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId))
        WHERE (n:VARIABLE OR n:FACTOR OR node:VARIABLE OR node:FACTOR) and not (node:SESSION or n:SESSION) and not (n:PPE) and not (node:PPE)
        WITH collect(n)+collect(node) as nodelist
        unwind nodelist as nodes
        return count(distinct nodes)"""
    @debug "[Query] $query"
    result = _queryNeo4j(dfg.neo4jInstance, query)
    # Neo4j.jl data structure sometimes feels brittle... like below
    return result.results[1]["data"][1]["row"][1] == length(varIds) + length(factIds)
end

function getNeighbors(dfg::CloudGraphsDFG, node::T; solvable::Int=0)::Vector{Symbol}  where T <: DFGNode
    query = "(n:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId):$(node.label))-[r:FACTORGRAPH]-(node) where (node:VARIABLE or node:FACTOR) and node.solvable >= $solvable"
    @debug "[Query] $query"
    neighbors = _getLabelsFromCyphonQuery(dfg.neo4jInstance, query)
    # If factor, need to do variable ordering TODO, Do we? does it matter if we always use _variableOrderSymbols in calculations?
    if T <: DFGFactor
        neighbors = intersect(node._variableOrderSymbols, neighbors)
    end
    return neighbors
end

function getNeighbors(dfg::CloudGraphsDFG, label::Symbol; solvable::Int=0)::Vector{Symbol}
    query = "(n:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId):$(label))-[r:FACTORGRAPH]-(node) where (node:VARIABLE or node:FACTOR) and node.solvable >= $solvable"
    @debug "[Query] $query"
    neighbors = _getLabelsFromCyphonQuery(dfg.neo4jInstance, query)
    # If factor, need to do variable ordering TODO, Do we? does it matter if we always use _variableOrderSymbols in calculations?
    if isFactor(dfg, label)
        # Server is authority
        serverOrder = Symbol.(JSON2.read(_getNodeProperty(dfg.neo4jInstance, [dfg.userId, dfg.robotId, dfg.sessionId, String(label)], "_variableOrderSymbols")))
        neighbors = intersect(serverOrder, neighbors)
    end
    return neighbors
end

function getSubgraphAroundNode(dfg::CloudGraphsDFG, node::DFGNode, distance::Int64=1, includeOrphanFactors::Bool=false, addToDFG::AbstractDFG=_getDuplicatedEmptyDFG(dfg); solvable::Int=0)::AbstractDFG
    distance < 1 && error("getSubgraphAroundNode() only works for distance > 0")

    # Thank you Neo4j for 0..* awesomeness!!
    neighborList = _getLabelsFromCyphonQuery(dfg.neo4jInstance,
        """
        (n:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId):$(node.label))-[FACTORGRAPH*0..$distance]-(node:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId))
        WHERE (n:VARIABLE OR n:FACTOR OR node:VARIABLE OR node:FACTOR)
        and not (node:SESSION)
        and (node.solvable >= $solvable or node:$(dfg.userId):$(dfg.robotId):$(dfg.sessionId):$(node.label))""" # Always return the root node
        )

    # Copy the section of graph we want
    _copyIntoGraph!(dfg, addToDFG, neighborList, includeOrphanFactors)
    return addToDFG
end

function getSubgraph(dfg::CloudGraphsDFG,
                     variableFactorLabels::Vector{Symbol},
                     includeOrphanFactors::Bool=false,
                     addToDFG::G=_getDuplicatedEmptyDFG(dfg) )::G where {G <: AbstractDFG}
    # Making a copy session if not specified

    _copyIntoGraph!(dfg, addToDFG, variableFactorLabels, includeOrphanFactors)

    return addToDFG
end

function getBiadjacencyMatrix(dfg::CloudGraphsDFG; solvable::Int=0)::NamedTuple{(:B, :varLabels, :facLabels),Tuple{SparseMatrixCSC,Vector{Symbol}, Vector{Symbol}}}
    varLabels = listVariables(dfg, solvable=solvable)
    factLabels = listFactors(dfg, solvable=solvable)
    vDict = Dict(varLabels .=> [1:length(varLabels)...])
    fDict = Dict(factLabels .=> [1:length(factLabels)...])

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

    # Now ask for all relationships for this session graph
    loadtx = transaction(dfg.neo4jInstance.connection)
    query = """
    START n=node(*)
    MATCH (n:$(join(_getLabelsForType(dfg, DFGVariable),':')))-[r:FACTORGRAPH]-(m:$(join(_getLabelsForType(dfg, DFGFactor),':')))
    where n.solvable >= $solvable and m.solvable >= $solvable
    RETURN n.label as variable, m.label as factor;"""
    @debug "[Query] $query"
    nodes = loadtx(query; submit=true)
    # Have to finish the transaction
    commit(loadtx)
    if length(nodes.errors) > 0
        error(string(nodes.errors))
    end
    # Add in the relationships
    varRel = Symbol.(map(node -> node["row"][1], nodes.results[1]["data"]))
    factRel = Symbol.(map(node -> node["row"][2], nodes.results[1]["data"]))
    for i = 1:length(varRel)
        adjMat[fDict[factRel[i]], vDict[varRel[i]]] = 1
    end

    return (B=adjMat, varLabels=varLabels, facLabels=factLabels)
end

### PPEs with DB calls

function listPPEs(dfg::CloudGraphsDFG, variablekey::Symbol; currentTransaction::Union{Nothing, Neo4j.Transaction}=nothing)::Vector{Symbol}
    query = "match (ppe:$(join(_getLabelsForType(dfg, MeanMaxPPE, parentKey=variablekey),':'))) return ppe.solverKey"
    @debug "[Query] PPE read query:\r\n$query"
    result = nothing
    if currentTransaction != nothing
        result = currentTransaction(query; submit=true)
    else
        tx = transaction(dfg.neo4jInstance.connection)
        tx(query)
        result = commit(tx)
    end
    length(result.errors) > 0 && error(string(result.errors))
    vals = map(d -> d["row"][1], result.results[1]["data"])
    return Symbol.(vals)
end

function getPPE(dfg::CloudGraphsDFG, variablekey::Symbol, ppekey::Symbol=:default; currentTransaction::Union{Nothing, Neo4j.Transaction}=nothing)::AbstractPointParametricEst
    query = "match (ppe:$(join(_getLabelsForType(dfg, MeanMaxPPE, parentKey=variablekey),':')):$(ppekey)) return properties(ppe)"
    @debug "[Query] PPE read query:\r\n$query"
    result = nothing
    if currentTransaction != nothing
        result = currentTransaction(query; submit=true)
    else
        tx = transaction(dfg.neo4jInstance.connection)
        tx(query)
        result = commit(tx)
    end
    length(result.errors) > 0 && error(string(result.errors))
    length(result.results[1]["data"]) != 1 && error("Cannot find PPE '$ppekey' for variable '$variablekey'")
    length(result.results[1]["data"][1]["row"]) != 1 && error("Cannot find PPE '$ppekey' for variable '$variablekey'")
    return unpackPPE(dfg, result.results[1]["data"][1]["row"][1])
end

function addPPE!(dfg::CloudGraphsDFG, variablekey::Symbol, ppe::P, ppekey::Symbol=:default; currentTransaction::Union{Nothing, Neo4j.Transaction}=nothing)::AbstractPointParametricEst where P <: AbstractPointParametricEst
    if ppekey in listPPEs(dfg, variablekey, currentTransaction=currentTransaction)
        error("PPE '$(ppekey)' already exists")
    end
    return _matchmergePPE!(dfg, variablekey, ppe, ppekey, currentTransaction=currentTransaction)
end

#TODO clean this, just to match api for update
function _matchmergePPE!(dfg::CloudGraphsDFG, variablekey::Symbol, ppe::P, ppekey::Symbol=:default; currentTransaction::Union{Nothing, Neo4j.Transaction}=nothing)::P where P <: AbstractPointParametricEst
    packed = packPPE(dfg, ppe)
    query = """
                MATCH (var:$variablekey:$(join(_getLabelsForType(dfg, DFGVariable, parentKey=variablekey),':')))
                MERGE (ppe:$(join(_getLabelsForInst(dfg, ppe, parentKey=variablekey),':')))
                SET ppe = $(_dictToNeo4jProps(packed))
                CREATE UNIQUE (var)-[:PPE]->(ppe)
                RETURN properties(ppe)"""
    @debug "[Query] PPE update query:\r\n$query"
    result = nothing
    if currentTransaction != nothing
        result = currentTransaction(query; submit=true) # TODO: Maybe we should submit (; submit = true) for the results to fail early?
    else
        tx = transaction(dfg.neo4jInstance.connection)
        tx(query)
        result = commit(tx)
    end
    length(result.errors) > 0 && error(string(result.errors))
    length(result.results[1]["data"]) != 1 && error("Cannot find PPE '$(ppe.solverKey)' for variable '$variablekey'")
    length(result.results[1]["data"][1]["row"]) != 1 && error("Cannot find PPE '$(ppe.solverKey)' for variable '$variablekey'")
    return unpackPPE(dfg, result.results[1]["data"][1]["row"][1])
end

function updatePPE!(dfg::CloudGraphsDFG, variablekey::Symbol, ppe::P, ppekey::Symbol=:default; currentTransaction::Union{Nothing, Neo4j.Transaction}=nothing)::P where P <: AbstractPointParametricEst
    if !(ppekey in listPPEs(dfg, variablekey, currentTransaction=currentTransaction))
        @warn "PPE '$(ppekey)' does not exist, adding"
    end
    return _matchmergePPE!(dfg, variablekey, ppe, ppekey, currentTransaction=currentTransaction)
end

function updatePPE!(dfg::CloudGraphsDFG, sourceVariables::Vector{<:DFGVariable}, ppekey::Symbol=:default; currentTransaction::Union{Nothing, Neo4j.Transaction}=nothing)
    tx = currentTransaction == nothing ? transaction(dfg.neo4jInstance.connection) : currentTransaction
    for var in sourceVariables
        updatePPE!(dfg, var.label, getPPE(var, ppekey), ppekey, currentTransaction=tx)
    end
    if currentTransaction == nothing
        result = commit(tx)
    end
    return nothing
end

function deletePPE!(dfg::CloudGraphsDFG, variablekey::Symbol, ppekey::Symbol=:default; currentTransaction::Union{Nothing, Neo4j.Transaction}=nothing)::AbstractPointParametricEst
    query = """
    match (ppe:$ppekey:$(join(_getLabelsForType(dfg, MeanMaxPPE, parentKey=variablekey),':')))
    with ppe, properties(ppe) as props
    detach delete ppe
    return props
    """
    @debug "[Query] PPE delete query:\r\n$query"
    result = nothing
    if currentTransaction != nothing
        result = currentTransaction(query; submit=true) # TODO: Maybe we should submit (; submit = true) for the results to fail early?
    else
        tx = transaction(dfg.neo4jInstance.connection)
        tx(query)
        result = commit(tx)
    end
    length(result.errors) > 0 && error(string(result.errors))
    length(result.results[1]["data"]) != 1 && error("Cannot find PPE '$ppekey' for variable '$variablekey'")
    length(result.results[1]["data"][1]["row"]) != 1 && error("Cannot find PPE '$ppekey' for variable '$variablekey'")
    return unpackPPE(dfg, @show result.results[1]["data"][1]["row"][1])
end

### Updated functions from AbstractDFG
### These functions write back as you add the data.

function addVariableSolverData!(dfg::CloudGraphsDFG,
                                variablekey::Symbol,
                                vnd::VariableNodeData,
                                solvekey::Symbol=:default)::VariableNodeData
    # TODO: Switch out to their own nodes, don't get the whole variable
    var = getVariable(dfg, variablekey)
    if haskey(var.solverDataDict, solvekey)
        error("VariableNodeData '$(solvekey)' already exists")
    end
    var.solverDataDict[solvekey] = vnd

    # TODO: Cleanup and consolidate as one function.
    solverDataDict = JSON2.write(Dict(keys(var.solverDataDict) .=> map(vnd -> packVariableNodeData(dfg, vnd), values(var.solverDataDict))))
    _setNodeProperty(
        dfg.neo4jInstance,
        _getLabelsForInst(dfg, var),
        "solverDataDict",
        solverDataDict)
    return vnd
end

function updateVariableSolverData!(dfg::CloudGraphsDFG,
                                   variablekey::Symbol,
                                   vnd::VariableNodeData,
                                   solvekey::Symbol=:default,
                                   useCopy::Bool=true,
                                   fields::Vector{Symbol}=Symbol[])::VariableNodeData
    # TODO: Switch out to their own nodes, don't get the whole variable
    var = getVariable(dfg, variablekey)
    if !haskey(var.solverDataDict, solvekey)
        @warn "VariableNodeData '$(solvekey)' does not exist, adding"
    end

    # Unnecessary for cloud, but (probably) being (too) conservative (just because).
    usevnd = useCopy ? deepcopy(vnd) : vnd
    # should just one, or many pointers be updated?
    if haskey(var.solverDataDict, solvekey) && isa(var.solverDataDict[solvekey], VariableNodeData) && length(fields) != 0
      # change multiple pointers inside the VND var.solverDataDict[solvekey]
      for field in fields
        destField = getfield(var.solverDataDict[solvekey], field)
        srcField = getfield(usevnd, field)
        if isa(destField, Array) && size(destField) == size(srcField)
          # use broadcast (in-place operation)
          destField .= srcField
        else
          # change pointer of destination VND object member
          setfield!(var.solverDataDict[solvekey], field, srcField)
        end
      end
    else
      # change a single pointer in var.solverDataDict
      var.solverDataDict[solvekey] = usevnd
    end

    # TODO: Cleanup and consolidate
    solverDataDict = JSON2.write(Dict(keys(var.solverDataDict) .=> map(vnd -> packVariableNodeData(dfg, vnd), values(var.solverDataDict))))
    _setNodeProperty(
        dfg.neo4jInstance,
        _getLabelsForInst(dfg, var),
        "solverDataDict",
        solverDataDict)
    return var.solverDataDict[solvekey]
end

function updateVariableSolverData!(dfg::CloudGraphsDFG,
                                   sourceVariables::Vector{<:DFGVariable},
                                   solvekey::Symbol=:default,
                                   useCopy::Bool=true,
                                   fields::Vector{Symbol}=Symbol[] )
    # TODO: Switch out to their own nodes, don't get the whole variable
    #TODO: Do in bulk for speed.
    for var in sourceVariables
        updateVariableSolverData!(dfg, var.label, getSolverData(var, solvekey), solvekey, useCopy, fields)
    end
end

function deleteVariableSolverData!(dfg::CloudGraphsDFG, variablekey::Symbol, solvekey::Symbol=:default)::VariableNodeData
    # TODO: Switch out to their own nodes, don't get the whole variable
    var = getVariable(dfg, variablekey)

    if !haskey(var.solverDataDict, solvekey)
        error("VariableNodeData '$(solvekey)' does not exist")
    end
    vnd = pop!(var.solverDataDict, solvekey)
    # TODO: Cleanup and consolidate... yadda yadda just do it already :)
    solverDataDict = JSON2.write(Dict(keys(var.solverDataDict) .=> map(vnd -> packVariableNodeData(dfg, vnd), values(var.solverDataDict))))
    _setNodeProperty(
        dfg.neo4jInstance,
        _getLabelsForInst(dfg, var),
        "solverDataDict",
        solverDataDict)
    return vnd
end

function mergeVariableSolverData!(dfg::CloudGraphsDFG, sourceVariable::DFGVariable)
    # TODO: Switch out to their own nodes, don't get the whole variable
    var = getVariable(dfg, sourceVariable.label)

    merge!(var.solverDataDict, deepcopy(sourceVariable.solverDataDict))

    # TODO: Cleanup and consolidate
    solverDataDict = JSON2.write(Dict(keys(var.solverDataDict) .=> map(vnd -> packVariableNodeData(dfg, vnd), values(var.solverDataDict))))
    _setNodeProperty(
        dfg.neo4jInstance,
        _getLabelsForInst(dfg, var),
        "solverDataDict",
        solverDataDict)
    return var
end

function getSolvable(dfg::CloudGraphsDFG, sym::Symbol)
    prop = _getNodeProperty(
        dfg.neo4jInstance,
        union(_getLabelsForType(dfg, isVariable(dfg, sym) ? DFGVariable : DFGFactor),[String(sym)]),
        "solvable")
    return prop
end

function setSolvable!(dfg::CloudGraphsDFG, sym::Symbol, solvable::Int)::Int
    prop = _setNodeProperty(
        dfg.neo4jInstance,
        union(_getLabelsForType(dfg, isVariable(dfg, sym) ? DFGVariable : DFGFactor),[String(sym)]),
        "solvable",
        solvable)
    return solvable
end

##==============================================================================
## TAGS as a set, list, merge, remove, empty
## CloudGraphsDFG functions
##==============================================================================
function mergeTags!(dfg::CloudGraphsDFG, sym::Symbol, tags::Vector{Symbol})

    if isVariable(dfg,sym)
        getNode = getVariable
        updateNode! = updateVariable!
    else
        getNode = getFactor
        updateNode! = updateFactor!
    end

    node = getNode(dfg, sym)
    union!(getTags(node), tags)

    updateNode!(dfg, node)

    return getTags(getNode(dfg, sym))

end

function removeTags!(dfg::CloudGraphsDFG, sym::Symbol, tags::Vector{Symbol})

    if isVariable(dfg,sym)
        getNode = getVariable
        updateNode! = updateVariable!
    else
        getNode = getFactor
        updateNode! = updateFactor!
    end

    node = getNode(dfg, sym)
    setdiff!(getTags(node), tags)

    updateNode!(dfg, node)

    return getTags(getNode(dfg, sym))

end

function emptyTags!(dfg::CloudGraphsDFG, sym::Symbol)

    if isVariable(dfg,sym)
        getNode = getVariable
        updateNode! = updateVariable!
    else
        getNode = getFactor
        updateNode! = updateFactor!
    end

    node = getNode(dfg, sym)
    empty!(getTags(node))

    updateNode!(dfg, node)

    return getTags(getNode(dfg, sym))

end

function RESERVED_mergeTags!(dfg::CloudGraphsDFG, sym::Symbol, tags::Vector{Symbol})

    nodeId = _tryGetNeoNodeIdFromNodeLabel(dfg.neo4jInstance,
                                            dfg.userId,
                                            dfg.robotId,
                                            dfg.sessionId,
                                            sym)

    neo4jNode = Neo4j.getnode(dfg.neo4jInstance.graph, nodeId)

    addnodelabels(neo4jNode, string.(tags))

    return Set(setdiff(Symbol.(getnodelabels(neo4jNode)), Symbol.([dfg.userId, dfg.robotId, dfg.sessionId]), [sym]))

end

function RESERVED_removeTags!(dfg::CloudGraphsDFG, sym::Symbol, tags::Vector{Symbol})

    nodeId = _tryGetNeoNodeIdFromNodeLabel(dfg.neo4jInstance,
                                            dfg.userId,
                                            dfg.robotId,
                                            dfg.sessionId,
                                            sym)

    neo4jNode = Neo4j.getnode(dfg.neo4jInstance.graph, nodeId)

    shouldStay = Symbol.([dfg.userId, dfg.robotId, dfg.sessionId]) ∪ [sym, :VARIABLE, :FACTOR]

    for tag in tags
        if tag in shouldStay
            @warn("Label:$tag is not allowed to be removed from tags and will be ignored.")
        else
            deletenodelabel(neo4jNode, string(tag))
        end
    end

    return Set(setdiff(Symbol.(getnodelabels(neo4jNode)), shouldStay))

end

function RESERVED_emptyTags!(dfg::CloudGraphsDFG, sym::Symbol)

    nodeId = _tryGetNeoNodeIdFromNodeLabel(dfg.neo4jInstance,
                                            dfg.userId,
                                            dfg.robotId,
                                            dfg.sessionId,
                                            sym)

    neo4jNode = Neo4j.getnode(dfg.neo4jInstance.graph, nodeId)

    shouldStay = Symbol.([dfg.userId, dfg.robotId, dfg.sessionId]) ∪ [sym, :VARIABLE, :FACTOR]
    tags = setdiff(Symbol.(getnodelabels(neo4jNode)), shouldStay)
    # tags = Symbol.(getnodelabels(neo4jNode))

    return removeTags!(dfg, sym, tags)
end