text:
# OccamNet
abstract type AbstractProbabilityLayer end
abstract type AbstractOccam end

## Overload softmax and logsoftmax with Temperature
NNlib.softmax(x, T; dims = 1) = T > eps() ? softmax(x ./ T, dims = dims) : softmax(x, dims = dims)
NNlib.logsoftmax(x, T; dims =1) = T > eps() ? x ./ T .- log.(sum(exp, x ./ T, dims = dims)) : logsoftmax(x, dims)

"""
$(TYPEDEF)

Defines a basic `ProbabilityLayer` in which the parameters act as probabilities via the `softmax` function for an
array of functions.

The layer is callable either via `layer(x)`, using all weights to form the output or by `layer(x, route)` where route is
the result of `rand(layer)` which samples the function arguments from the underlying distribution.

# Fields

$(FIELDS)

"""
mutable struct ProbabilityLayer{F,W,T,A} <: AbstractProbabilityLayer
    "Nonlinear functions forming the basis of the layer"
    op::F
    "Weights"
    weight::W
    "Temperature controlling the shape of the distribution"
    t::T
    "Arities of the functions"
    arieties::A
    "Skip connection"
    skip::Bool
end

function ProbabilityLayer(inp::Int, f::Function, t::Real = 1.0; skip = false, init_w = ones)
    arieties = ariety(f, typeof(t))
    w = init_w(eltype(t), sum(arieties), inp)
    return ProbabilityLayer([f], w, t, [arieties], skip)
end

function ProbabilityLayer(inp::Int, out::Int, t::Real = 1.0; skip = false, init_w = ones)
    arieties = [out]
    w = init_w(eltype(t), out, inp)
    return ProbabilityLayer(identity, w, t, arieties, skip)
end

function ProbabilityLayer(inp::Int, f::AbstractVector{Function}, t::Real = 1.0; skip = false, init_w = ones)
    arieties = [ariety(fi,typeof(t)) for fi in f]
    w = init_w(eltype(t), sum(arieties), inp)
    return ProbabilityLayer(f, w, t, arieties, skip)
end

function (p::AbstractProbabilityLayer)(x::AbstractVector)
     h̃ = softmax(p.weight, p.t, dims = 2)*x
     idx = 1
     res = map(1:length(p.arieties)) do i
         # Index of input
         args_ = h̃[idx:idx+p.arieties[i]-1]
         idx += p.arieties[i]
         p.op[i](args_...)
     end
     !p.skip && return res
     [res;x]
end

function (p::AbstractProbabilityLayer)(x::AbstractVector, sample::Vector{Int64})
    h̃ = [x[i] for i in sample]
    idx = 1
    res = map(1:length(p.arieties)) do i
        # Index of input
        args_ = h̃[idx:idx+p.arieties[i]-1]
        idx += p.arieties[i]
        p.op[i](args_...)
    end
    !p.skip && return res
    [res;x]
end

function (p::ProbabilityLayer{typeof(identity), W, T, A})(x::AbstractVector) where {W, T, A}
    res = softmax(p.weight, p.t, dims = 2)*x
    !p.skip && return res
    [res;x]
end

function (p::ProbabilityLayer{typeof(identity), W, T, A})(x::AbstractVector, sample::Vector{Int64}) where {W, T, A}
    res = [x[i] for i in sample]
    !p.skip && return res
    [res;x]
end

(p::AbstractProbabilityLayer)(x::AbstractMatrix) = hcat(map(p, eachcol(x))...)

"""
$(SIGNATURES)

Return the probability associated with the `ProbabilityLayer` or `OccamNet` by applying `softmax` on the weights.
"""
probabilities(p::AbstractProbabilityLayer) = softmax(p.weight, p.t, dims = 2)

function probabilities(p::AbstractProbabilityLayer, sample::Vector{Int})
    probs = probabilities(p)
    map(i->getindex(probs, i, sample[i]), 1:length(sample))
end

"""
$(SIGNATURES)

Return the logprobability associated with the `ProbabilityLayer` or `OccamNet` by applying `logsoftmax` on the weights.
"""
logprobabilities(p::AbstractProbabilityLayer) = logsoftmax(p.weight, p.t, dims = 2)

function logprobabilities(p::AbstractProbabilityLayer, sample::Vector{Int})
    logprobs = logprobabilities(p)
    map(i->getindex(logprobs, i, sample[i]), 1:length(sample))
end

# The resulting logprob
function logprobability(p::AbstractProbabilityLayer, sample, x)
    h̃ = logprobabilities(p, sample) + [x[i] for i in sample] # Add the input logprob to the output
    idx = 1
    res = map(1:length(p.arieties)) do i
        # Index of input
        logp = sum(h̃[idx:idx+p.arieties[i]-1])
        idx += p.arieties[i]
        logp
    end
    !p.skip && return res
    [res;x]
end

function logprobability(p::ProbabilityLayer{typeof(identity), W, T, A}, sample, x) where {W,T,A}
    h̃ = logprobabilities(p, sample) + [x[i] for i in sample]
    !p.skip && return h̃
    [h̃;x]
end

"""
$(SIGNATURES)

Return the probability of a pathway defined by sample associated with the object by applying `softmax` on the weights.
"""
resprobability(p::AbstractProbabilityLayer, sample) = exp.(reslogprobability(p, sample))

"""
$(SIGNATURES)

Set the temperature of the `ProbabilityLayer` or `OccamNet`.
"""
set_temp!(p::AbstractProbabilityLayer, t::Real) = p.t = t

Base.rand(p::AbstractProbabilityLayer) = [rand(Categorical(Vector(w))) for w in eachrow(probabilities(p))]

function Base.show(io::IO, l::ProbabilityLayer)
  print(io, "ProbabilityLayer(", size(l.weight, 2), ", ", size(l.weight, 1))
  print(io, ", $(l.op)")
  l.skip ? print(io, ", Skip") : nothing
  print(io, ")")
end

Flux.@functor ProbabilityLayer
Flux.trainable(u::ProbabilityLayer) = (u.weight,)

## NETWORK

"""
$(TYPEDEF)

Defines an `OccamNet` which learns symbolic expressions from data using a probabilistic approach.
See [Interpretable Neuroevolutionary Models for Learning Non-Differentiable Functions and Programs
](https://arxiv.org/abs/2007.10784) for more details.

It get constructed via:

```julia
net = OccamNet(inp::Int, outp::Int, layers::Int, f::Vector{Function}, t::Real = 1.0; constants = typeof(t)[], parameters::Int = 0, skip::Bool = false, init_w = ones, init_p = Flux.glorot_uniform)
```

`inp` describes the size of the input domain, `outp` the size of the output domain, `layers` the number of layers (including the input layer and excluding the linear output layer) and
`f` the functions to be used. Optional is the temperature `t` which is set to `1.0` at the beginning.

Keyworded arguments are `constants`, a vector of constants like π, ℯ which can concatenated to the input, the number of trainable `parameters` and if `skip` connections should be used.
The constructors to the weights and parameters can be passed in via `init_w` and `init_p`.

`OccamNet` is callable with and without a specific route, which can be sampled from the networks weights via `rand(net)`.

# Fields
$(FIELDS)

"""
mutable struct OccamNet{F, C, P} <: AbstractOccam
    "The `Chain` representing the network"
    c::F # The chain
    "Additional constants added to the input which are not trainable."
    constants::C # Additional constants which are fixed
    "Additional parameters added to the input which are trainable."
    parameters::P # Additional learnable parameters
end

function OccamNet(inp::Int, outp::Int, layers::Int, f::Vector{Function}, t::Real = 1.0; constants = typeof(t)[], parameters::Int = 0, skip::Bool = false, init_w = ones, init_p = Flux.glorot_uniform)
    inp += parameters
    inp += length(constants)
    x0 = rand(typeof(t), inp)
    ls = []
    for l in 1:layers
        l_ = ProbabilityLayer(inp, f, t, skip = skip, init_w = init_w)
        push!(ls, l_)
        x0 = l_(x0)
        inp = length(x0)
    end
    # Output layer
    push!(ls, ProbabilityLayer(inp, outp, t, skip = false, init_w = init_w))
    b = parameters > 0 ? convert.(typeof(t), init_p(parameters)) : typeof(t)[]
    return OccamNet(Chain(ls...),constants, b)
end

function Base.show(io::IO, l::OccamNet)
  print(io, "OccamNet(", length(l.c))
  !isempty(l.constants) ? print(io, ", Constants ", length(l.constants)) : nothing
  !isempty(l.constants) ? print(io, ", Parameters ", length(l.parameters)) : nothing
  print(io, ")")
end

(o::OccamNet)(x) = o.c([x;o.constants;o.parameters])
(o::OccamNet)(x::AbstractMatrix) = hcat(map(o, eachcol(x))...)
set_temp!(o::OccamNet, t::Real) = map(x->set_temp!(x, t), o.c)
probabilities(o::OccamNet) = map(probabilities, o.c)
probabilities(o::OccamNet, route) = map(probabilities, o.c, route)

logprobabilities(o::OccamNet) = map(logprobabilities, o.c)
logprobabilities(o::OccamNet, route) = map(logprobabilities, o.c, route)

Base.rand(o::OccamNet) = map(rand, o.c)

function (o::OccamNet)(x::AbstractVector, route::Vector{Vector{Int64}})
    res = [x;o.constants;o.parameters]
    for i in 1:length(route)
        res = o.c[i](res,route[i])
    end
    return res
end

(o::OccamNet)(x::AbstractMatrix, route::Vector{Vector{Int64}}) = hcat(map(xi->o(xi, route), eachcol(x))...)

"""
$(SIGNATURES)

Returns the logprobability of the result of the `OccamNet` using the specific route.
"""
function logprobability(o::OccamNet, route::Vector{Vector{Int64}})
    # Input prob = 1 -> logprob = 0
    res = zeros(eltype(o.c[1].t), size(o.c[1].weight, 2) + length(o.constants))
    for i in 1:length(route)
        res = logprobability(o.c[i], route[i], res)
    end
    return res
end

"""
$(SIGNATURES)

Returns the probability of the result of the `OccamNet` using the specific route.
"""
probability(o::OccamNet, route::Vector{Vector{Int64}}) = exp.(logprobability(o, route))

Flux.@functor OccamNet
Flux.trainable(u::OccamNet) = (u.parameters, Flux.trainable.(u.c)...,)

gkernel(x,y,σ = one(eltype(x))) = 1/sqrt(2π*σ^2) .* exp.(-(x-y).^2 ./(2*σ))

function gkernel(x::AbstractMatrix, y::AbstractMatrix, σ::AbstractVector = ones(eltype(x), size(x, 1)))
    vcat(map(i->gkernel(x[i,:], y[i,:], σ[i])', 1:size(x, 1))...)
end

"""
$(SIGNATURES)

Overloads `Flux.train!` method to be used with an `OccamNet`.
"""
function Flux.train!(net::OccamNet, X, Y, opt, maxiters = 10; routes = 10, nbest = 1, cb = ()->(), progress = false)
    ny = size(Y,1)
    vary = [var(Y, dims = 2)...]
    ps_prob = Flux.params(net)

    if progress
        prog = Progress(
            maxiters, "Training $(net)"
        )
    end

    for k in 1:maxiters
        ls = map(1:routes) do i
            route = rand(net)
            res = net(X, route)
            route, sum(gkernel(res, Y, vary))
        end
        ls = sort(ls, by = last, rev = true)

        candidates = first.(ls[1:nbest])

        # Update probabilities
        gs = Flux.gradient(ps_prob) do
            l_ = zero(eltype(X))
            for c in candidates
                res = net(X, c)
                p = logprobability(net, c)
                l_ += dot(p, sum(gkernel(res, Y, vary), dims =2))
            end
            -l_
        end
        Flux.Optimise.update!(opt, ps_prob, gs)

        if progress
            rp = round.(exp.(logprobability(net, first(first(ls)))), digits = 5)
            loss = sum(abs2, net(X, first(first(ls))) - Y) / size(Y, 2)

            ProgressMeter.next!(
            prog;
            showvalues = [
                (:Probabilities, rp), (Symbol("Equivalent L2-Loss"), loss)
                ]
                )
        end
    end
    return
end

## SR WRAPPER

"""
$(TYPEDEF)

Options for using OccamNet within the `solve` function. Automatically creates a network with the given specification.

# Fields
$(FIELDS)

"""
struct OccamSR{F, C, T} <: AbstractSymbolicRegression "Functions used within the network"
    functions::F
    "Constants added to the input"
    constants::C
    "Number of layers"
    layers::Int
    "Number of parameters"
    parameters::Int
    "Activate skip connections"
    skip::Bool
end

function Base.show(io::IO, l::OccamSR{F, C, T}) where {F,C,T}
  T ? print(io, "Implicit ") : nothing
  print(io, "OccamSR(", length(l.functions))
  !isempty(l.constants) ? print(io, ", Constants ", length(l.constants)) : nothing
  l.parameters > 0 ? print(io, ", Parameters ", l.parameters) : nothing
  print(io, ")")
end

function Base.summary(io::IO, l::OccamSR)
  print(io, "OccamSR\n")
  print(io, "Functions ", l.functions, "\n")
  !isempty(l.constants) ? print(io, "Constants ", l.constants, "\n") : nothing
  l.parameters > 0 ? print(io, "Parameters ", l.parameters, "\n") : nothing
end

function Base.print(io::IO, l::OccamSR)
    summary(io, l)
end

function OccamSR(;functions = [+,*,sin,exp], constants = [π, ℯ],
    layers = 2, parameters = 0, skip = true, implicit = false,
    kwargs...)
    implicit && throw(error("Implicit OccamSR is not supported at the moment."))
    return OccamSR{typeof(functions), typeof(constants), implicit}(functions, constants, layers, parameters, skip)
end

## SOLVE

function DiffEqBase.solve(
    p, o::OccamSR{F,C, false}, opt;
    max_iter = 1000, cb = ()->(), progress = false, routes = 10, nbest = 1, temperature = 1.0
    ) where {F,C}

    # Target variables
    Y = get_target(p)
    # Inputs
    X̂, _, t, U = get_oop_args(p)

    t = iszero(t) ? [] : t
    # Cat the inputs
    X = vcat([x for x in (X̂, U, permutedims(t)) if !isempty(x)]...)

    inp = size(X,1)
    outp = size(Y,1)

    net = OccamNet(inp, outp, o.layers, o.functions, temperature, constants = o.constants, parameters = o.parameters, skip = o.skip)

    Flux.train!(net, X, Y, opt, max_iter, routes = routes, nbest = nbest, cb = cb, progress = progress)
    build_solution(p, net, o, opt)
end

## SOLUTION
function build_solution(prob::DataDrivenProblem, net::OccamNet, o::OccamSR, opt;
    eval_expression = false)

    @variables x[1:size(prob.X, 1)] u[1:size(prob.U,1)] t
    x = collect(x)
    u = collect(u)
    x_ = [x;u;t]

    inp = size(net.c[1].weight, 2)- length(net.constants)
    temp_ = net.c[1].t

    # Draw routers
    set_temp!(net, 0.1)
    route = rand(net)
    eqs = simplify.(net(x_[1:inp], route))
    set_temp!(net, temp_)

    lhs, dt = assert_lhs(prob)

    # Build the lhs
    if (length(eqs) == length(x)) && (lhs != :direct)
            if lhs == :continuous
                d = Differential(t)
            elseif lhs == :discrete
                d = Difference(t, dt = dt)
            end
            eqs = [d(xs[i]) ~ eq for (i,eq) in enumerate(eqs)]
    end

    # Build a basis
    res_ = Basis(
        eqs, x, iv = t,
        controls = u,
        eval_expression = eval_expression
    )

    X = get_target(prob)
    Y = res_(get_oop_args(prob)...)

    # Build the metrics
    pb = exp(sum(logprobability(net, route)))
    #pbs = probability(net, route)
    retcode = Symbol("$(pb)")

    error = sum(abs2, X-Y, dims = 2)[:,1]
    aic = 2*(-size(error, 2) .* log.(error ./ size(error, 2)) .+ sum(length, route))

    return DataDrivenSolution(
        false, res_, [], retcode, o, net, prob, error, aic
    )
end


text:
function initialize_filter!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::DynamicWrapper{PSY.DynamicInverter{C, O, IC, DC, P, PSY.LCLFilter}},
    inner_vars::AbstractVector,
) where {
    C <: PSY.Converter,
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
}
    #PowerFlow Data
    P0 = PSY.get_active_power(static)
    Q0 = PSY.get_reactive_power(static)
    Vm = PSY.get_magnitude(PSY.get_bus(static))
    θ = PSY.get_angle(PSY.get_bus(static))
    S0 = P0 + Q0 * 1im
    V_R = Vm * cos(θ)
    V_I = Vm * sin(θ)
    V = V_R + V_I * 1im
    I = conj(S0 / V)
    I_R = real(I)
    I_I = imag(I)

    #Get Parameters
    filter = PSY.get_filter(dynamic_device)
    lf = PSY.get_lf(filter)
    rf = PSY.get_rf(filter)
    cf = PSY.get_cf(filter)
    lg = PSY.get_lg(filter)
    rg = PSY.get_rg(filter)

    #Set parameters
    Ir_filter = I_R
    Ii_filter = I_I
    ω_sys = get_ω_ref(dynamic_device)

    #To solve Vr_cnv, Vi_cnv, Ir_cnv, Ii_cnv, Vr_filter, Vi_filter
    function f!(out, x)
        Vr_cnv = x[1]
        Vi_cnv = x[2]
        Ir_cnv = x[3]
        Ii_cnv = x[4]
        Vr_filter = x[5]
        Vi_filter = x[6]

        #𝜕Ir_cnv/𝜕t
        out[1] = Vr_cnv - Vr_filter - rf * Ir_cnv + ω_sys * lf * Ii_cnv
        #𝜕Ii_cnv/𝜕t
        out[2] = Vi_cnv - Vi_filter - rf * Ii_cnv - ω_sys * lf * Ir_cnv
        #𝜕Vr_filter/𝜕t
        out[3] = Ir_cnv - Ir_filter + ω_sys * cf * Vi_filter
        #𝜕Vi_filter/𝜕t
        out[4] = Ii_cnv - Ii_filter - ω_sys * cf * Vr_filter
        #𝜕Ir_filter/𝜕t
        out[5] = Vr_filter - V_R - rg * Ir_filter + ω_sys * lg * Ii_filter
        #𝜕Ii_filter/𝜕t
        out[6] = Vi_filter - V_I - rg * Ii_filter - ω_sys * lg * Ir_filter
    end
    x0 = [V_R, V_I, I_R, I_I, V_R, V_I]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Filter failed")
    else
        sol_x0 = sol.zero
        #Update terminal voltages
        inner_vars[Vr_inv_var] = V_R
        inner_vars[Vi_inv_var] = V_I
        #Update Converter voltages
        inner_vars[Vr_cnv_var] = sol_x0[1]
        inner_vars[Vi_cnv_var] = sol_x0[2]
        #Update filter voltages
        inner_vars[Vr_filter_var] = sol_x0[5]
        inner_vars[Vi_filter_var] = sol_x0[6]
        #Update states
        filter_ix = get_local_state_ix(dynamic_device, PSY.LCLFilter)
        filter_states = @view device_states[filter_ix]
        filter_states[1] = sol_x0[3] #Ir_cnv
        filter_states[2] = sol_x0[4] #Ii_cnv
        filter_states[3] = sol_x0[5] #Vr_filter
        filter_states[4] = sol_x0[6] #Vi_filter
        filter_states[5] = I_R #Ir_filter
        filter_states[6] = I_I #Ii_filter
    end
end

function initialize_filter!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::DynamicWrapper{PSY.DynamicInverter{C, O, IC, DC, P, PSY.RLFilter}},
    inner_vars::AbstractVector,
) where {
    C <: PSY.Converter,
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
}
    #PowerFlow Data
    P0 = PSY.get_active_power(static)
    Q0 = PSY.get_reactive_power(static)
    Vm = PSY.get_magnitude(PSY.get_bus(static))
    θ = PSY.get_angle(PSY.get_bus(static))
    S0 = P0 + Q0 * 1im
    V_R = Vm * cos(θ)
    V_I = Vm * sin(θ)
    V = V_R + V_I * 1im
    # Doesn't match PSS/e initialization because PTI doesn't use the conjugate.
    I = conj(S0 / V)
    # PSS/e names I_R as Ip. But is calculated as P/Vt
    I_R = real(I)
    # PSS/e names I_I as Iq. But is calculated as Q/Vt
    I_I = imag(I)

    # Update Control References
    PSY.set_Q_ref!(PSY.get_converter(dynamic_device), Q0)
    set_Q_ref(dynamic_device, Q0)
    PSY.set_Q_ref!(PSY.get_reactive_power(PSY.get_outer_control(dynamic_device)), Q0)
    PSY.set_P_ref!(PSY.get_active_power(PSY.get_outer_control(dynamic_device)), P0)
    set_P_ref(dynamic_device, P0)
    PSY.set_V_ref!(PSY.get_reactive_power(PSY.get_outer_control(dynamic_device)), Vm)
    set_V_ref(dynamic_device, Vm)

    #Get Parameters
    filt = PSY.get_filter(dynamic_device)
    rf = PSY.get_rf(filt)
    lf = PSY.get_lf(filt)
    converter = PSY.get_converter(dynamic_device)
    R_source = PSY.get_R_source(converter)
    X_source = PSY.get_X_source(converter)

    #Update terminal voltages
    inner_vars[Vr_inv_var] = V_R
    inner_vars[Vi_inv_var] = V_I
    #Update filter currents (output of converter)
    inner_vars[Ir_inv_var] = I_R
    inner_vars[Ii_inv_var] = I_I
    #Update converter currents
    V_cnv = V + (rf + lf * 1im) * I
    I_aux = V_cnv / (R_source + X_source * 1im)
    I_cnv = I + I_aux

    inner_vars[Vr_cnv_var] = real(V_cnv)
    inner_vars[Vi_cnv_var] = imag(V_cnv)
    inner_vars[Ir_cnv_var] = real(I_cnv)
    inner_vars[Ii_cnv_var] = imag(I_cnv)
    return
end

ext: jl


text:
import Base.show

"""

# DAG

Directed acyclic graph struct

### Struct
```julia
DAG(
* `name::AbstractString`                    : Name for the DAG object
* `d::OrderedDictOrNothing`                 : DAG definition as an OrderedDict
* `a::NamedArrayOrNothing`                  : Adjacency matrix
* `e::NamedArrayOrNothing`                  : Edge matrix
* `s::NamedArrayOrNothing`                  : Covariance matrix
* `df::DataFrameOrNothing`                  : Variable observations
* `vars::Vector{Symbol}`                    : Names of variables in DAG
)
```

Part of API, exported.
"""
mutable struct DAG
  name::AbstractString
  d::OrderedDictOrNothing                    # Name of DAG
  a::NamedArrayOrNothing                     # Adjacency matrix
  e::NamedArrayOrNothing                     # Edge matrix
  s::NamedArrayOrNothing                     # Covariance matrix as NamedArray
  df::DataFrameOrNothing                     # DataFrame with variables
  vars::Vector{Symbol}                       # Names of variables in DAG
end

# Constructor

"""

# DAG

Directed acyclic graph constructor

$(SIGNATURES)

### Required arguments
```julia
* `name::AbstractString`               : Name for the DAG object
* `d::ModelDefinition`                 : DAG definition
```

where
```
ModelDefinition = Union{OrderedDict, AbstractString, NamedArray}
```

See the extended help for a usage example.

### Keyword arguments
```julia
* `df::DataFrame`                      : DataFrame with observations
```

### Returns
```julia
* `dag::DAG`                           : Boolean result of test
```

# Extended help

In the definition of the OrderedDict, read `=>` as `~` in regression models
or `<-` in causal models, e.g.:

```julia
d = OrderedDict(
  :u => [:x, :v],
  :s1 => [:u],
  :w => [:v, :y],
  :s2 => [:w]
);
dag = DAG("my_name", d)
```

Coming from R's dagitty:

amat <- dagitty("dag { {X V} -> U; S1 <- U; {Y V} -> W; S2 <- W}”)
```julia
dag = DAG("my_name", "dag { {X V} -> U; S1 <- U; {Y V} -> W; S2 <- W}”)
display(dag) # Show the DAG
```

Coming from R's ggm:

amat <- DAG(U~X+V, S1~U, W~V+Y, S2~W, order=FALSE)
```julia
dag = DAG("my_name", "DAG(U~X+V, S1~U, W~V+Y, S2~W”)
display(dag) # Show the DAG
```

### Acknowledgements

Original author:                       Giovanni M. Marchetti

Translated to Julia:                   Rob J Goedman

### License

The R package ggm is licensed under License: GPL-2.

The Julia translation is licenced under: MIT.

Part of API, exported.
"""
function DAG(name::AbstractString, model::ModelDefinition; df::DataFrameOrNothing=nothing)

  local d
  if typeof(model) <: OrderedDict
    d = model
  elseif typeof(model) <: AbstractString
    ds = strip(model)
    if ds[1:3] == "DAG"
      d = from_ggm(ds)
    elseif ds[1:3] == "dag"
      d = from_dagitty(ds)
    else
      @error "Unrecognized model string: $(ds))"
    end
  elseif typeof(model) <: NamedArray
    d = adjacency_matrix_to_dict(model)
  end

  vars = dag_vars(d)
  a = adjacency_matrix(d)
  e = edge_matrix(d)

  if isnothing(df)
    s = nothing
  else
    # Compute covariance matrix and store as NamedArray
    @assert length(names(df)) == length(vars) "DataFrame has different number of columns"
    s = NamedArray(cov(Array(df)), (names(df), names(df)), ("Rows", "Cols"))
  end

  # Create object

  DAG(name, d, a, e, s, df, vars)

end

"""

# `set_dag_df!`

Set or update Dataframe associated to DAG

$(SIGNATURES)

### Required arguments
```julia
* `d::DAG`                                  : Previously defined DAG object
* `df::DataFrameOrNothing`                  : DataFrame associated with DAG
)
```

### Optional arguments
```julia
* `force=false`                             : Force assignment of df
)
```

The `force = true` option can be used if the DAG involves unobserved nodes.

Part of API, exported.
"""
function set_dag_df!(d::DAG, df::DataFrameOrNothing; force=false)
  # Compute covariance matrix and store as NamedArray

  if !(force || df == nothing)
    @assert length(names(df)) == length(d.vars) "DataFrame has different number of columns"
    @assert names(df, 1) !== d.vars "DataFrame names differ from DAG variable names"
  end

  d.df = df
  if df == nothing || nrow(df) == 0
    d.s = nothing
  else
    if nrow(df) > 1 && length(names(df)) > 0
      d.s = NamedArray(cov(Array(df)), (names(df), names(df)), ("Rows", "Cols"))
    else
      d.s = nothing
    end
  end

end

"""

# `set_dag_cov_matrix!`

Set or update the covariance matrix associated to DAG

$(SIGNATURES)

### Required arguments
```julia
* `d::DAG`                                  : Previously defined DAG object
* `cm::NamedArrayOrNothing`                 : Covariance matrix in NamedArray format
)
```

### Optional arguments
```julia
* `force=false`                             : Force assignment of df
)
```

The `force = true` option can be used if the DAG involves unobserved nodes.

Part of API, exported.
"""
function set_dag_cov_matrix!(d::DAG, cm::NamedArrayOrNothing; force=false)
  # Compute covariance matrix and store as NamedArray

  if !(force || cm == nothing)
    @assert length(names(cm)) == length(d.vars) "Covariance matrix has different number of columns"
    @assert names(cm) !== d.vars "Covariance matrix names differ from DAG variable names"
  end

  if cm == nothing
    d.s = nothing
  else
    d.s = cm
  end

end

"""

# `dag_show`

$(SIGNATURES)

Internal
"""
function dag_show(io::IO, d::DAG)
  println("\nDAG object:\n")
  println(io, "name = \"$(d.name)\"")
  println(io, "vars = $(d.vars)")
  println()
  !isnothing(d.d) && display(d.d)
  println()
end

show(io::IO, d::DAG) = dag_show(io, d)

export
  DAG,
  set_dag_df!,
  set_dag_cov_matrix!


text:
# functions to write an MP3 file, uses LAME under the hood

include("lame.jl")

function initialize_writers()
    # nop
end

# number of frames in 1 MP3 block
const MP3_BLOCKFRAMES = 1152
# max size needed to hold the encoded mp3 data for a full block
const MP3_BUFBYTES = ceil(Int, 1.25 * MP3_BLOCKFRAMES + 7200)

mutable struct MP3FileSink <: SampleSink
    lame::LAME
    samplerate::Int
    nchannels::Int
    output::IO
    nframes::Int
    mp3buf::Vector{UInt8}
end

function MP3FileSink(lame, samplerate, nchannels, output)
    MP3FileSink(lame, samplerate, nchannels, output, 0, Array{UInt8}(undef, MP3_BUFBYTES))
end

@inline nchannels(sink::MP3FileSink) = sink.nchannels
@inline samplerate(sink::MP3FileSink) = sink.samplerate
@inline nframes(sink::MP3FileSink) = sink.nframes
@inline Base.eltype(sink::MP3FileSink) = PCM16Sample

"""
save an MP3 file, using parameters as specified

# Arguments
* `file::File`: the MP3 file to save
* `buf::SampleBuf`: the audio buffer to save
* `nchannels::Int`: whether the output should be mono (1) or stereo (2)
* `samplerate::Int`: the sample rate, default is to retain the original;
                     will raise an error if the value is not one of
                     8000, 11025, 12000, 16000, 22050, 24000, 32000, 44100, 48000
* `bitrate::Int`: the bitrate in kbps, used for CBR
* `VBR::Bool`: whether to use VBR (false by default)
* `quality::Int`: 0 for highest quality, 9 for lowest quality (default 4)
                  used as algorithm selector for CBR and quality level for VBR.
* others: ID3v2 tag items to be added
"""
function save(file::File, buf::SampleBuf;
        samplerate=SampledSignals.samplerate(buf), nchannels=SampledSignals.nchannels(buf),
        kwargs...)
    savestream(file; samplerate=samplerate, nchannels=nchannels, kwargs...) do stream
        frameswritten = write(stream, buf)
        if frameswritten != nframes(buf)
            error("Only wrote $frameswritten frames, expected $(nframes(buf))")
        end
    end

    nothing
end

@inline savestream(path::AbstractString, args...; kwargs...) =
    savestream(query(path), args...; kwargs...)

function savestream(f::Function, args...; kwargs...)
    stream = savestream(args...; kwargs...)
    try
        f(stream)
    finally
        close(stream)
    end
end

function savestream(path::File;
                    nchannels = 2,
                    samplerate = 44100,
                    bitrate = 320,
                    VBR = false,
                    quality = 4,
                    title = "",
                    artist = "",
                    album = "",
                    year = "",
                    comment = "")

    lame = lame_init()
    # lame_set_num_samples(lame, info.nframes)

    if nchannels == 1
        lame_set_num_channels(lame, 1)
        lame_set_mode(lame, LAME_MONO)
    elseif nchannels == 2
        lame_set_num_channels(lame, 2)
        lame_set_mode(lame, LAME_JOINT_STEREO)
    else
        error("the output channels should be either mono (1) or stereo (2)")
    end

    if !(samplerate in [8000, 11025, 12000, 16000, 22050, 24000, 32000, 44100, 48000])
        error("sample rate $samplerate Hz is not supported")
    end
    # resampling is handled by SampledSignals
    lame_set_in_samplerate(lame, Int(samplerate))
    lame_set_out_samplerate(lame, Int(samplerate))

    if VBR == false
        # CBR mode
        lame_set_brate(lame, Int(bitrate))
        lame_set_quality(lame, Int(quality))
    else
        # VBR mode
        lame_set_VBR(lame)
        lame_set_VBR_quality(lame, Int(quality))
    end

    id3tag_init(lame)
    id3tag_add_v2(lame)
    id3tag_v2_only(lame)
    isempty(title)   || id3tag_set_title(lame, title)
    isempty(artist)  || id3tag_set_artist(lame, artist)
    isempty(album)   || id3tag_set_album(lame, album)
    isempty(year)    || id3tag_set_year(lame, year)
    isempty(comment) || id3tag_set_comment(lame, comment)
    lame_set_write_id3tag_automatic(lame, 0)

    lame_init_params(lame)

    id3size = lame_get_id3v2_tag(lame, UInt8[], 0)
    id3buffer = Array{UInt8}(undef, id3size)
    lame_get_id3v2_tag(lame, id3buffer, id3size) == id3size || error("failed to get ID3v2 buffer")

    output = open(filename(path), "w")
    write(output, id3buffer)

    MP3FileSink(lame, samplerate, nchannels, output)
end

function unsafe_write(sink::MP3FileSink, buf::Array, frameoffset, framecount)
    lame = sink.lame

    # the data in the buffer is not interleaved; we pass them separately
    encsize = sizeof(PCM16Sample)
    left = channelptr(buf, 1, frameoffset)
    # the right-channel pointer only gets accessed for a stereo file, in which
    # case we should have a stereo buffer (because of the SampledSignals
    # conversion)
    right = channelptr(buf, 2, frameoffset)

    mp3buf = pointer(sink.mp3buf)

    written = 0
    while written < framecount
        nframes = min(MP3_BLOCKFRAMES, framecount - written)
        l = left + written * encsize
        r = right + written * encsize
        bytes = lame_encode_buffer!(lame, l, r, nframes, mp3buf, MP3_BUFBYTES)
        Base.unsafe_write(sink.output, mp3buf, bytes)

        written += nframes
    end

    sink.nframes += written
    written
end

function Base.close(sink::MP3FileSink)
    if sink.lame != C_NULL
        mp3buf = pointer(sink.mp3buf)
        bytes = lame_encode_flush_nogap(sink.lame, mp3buf, MP3_BUFBYTES)
        Base.unsafe_write(sink.output, mp3buf, bytes)

        err = lame_close(sink.lame)
        close(sink.output)
        if err != 0
            error("Could not close LAME handle: ", err)
        end
        sink.lame = C_NULL
    end
end


text:
# For all types that pack their type into their own structure (e.g. PPE)
const TYPEKEY = "_type"

## Custom serialization
using JSON
import JSON.show_json
import JSON.Writer: StructuralContext, JSONContext, show_json
import JSON.Serializations: CommonSerialization, StandardSerialization
JSON.show_json(io::JSONContext, serialization::CommonSerialization, uuid::UUID) = print(io.io, "\"$uuid\"")

## Version checking
function _getDFGVersion()
    if haskey(Pkg.dependencies(), Base.UUID("b5cc3c7e-6572-11e9-2517-99fb8daf2f04"))
        return string(Pkg.dependencies()[Base.UUID("b5cc3c7e-6572-11e9-2517-99fb8daf2f04")].version)
    else
        # This is arguably slower, but needed for Travis.
        return Pkg.TOML.parse(read(joinpath(dirname(pathof(@__MODULE__)), "..", "Project.toml"), String))["version"]
    end
end

function _versionCheck(props::Dict{String, Any})
    if haskey(props, "_version")
        if props["_version"] != _getDFGVersion()
            @warn "This data was serialized using DFG $(props["_version"]) but you have $(_getDFGVersion()) installed, there may be deserialization issues."
        end
    else
        @warn "There isn't a version tag in this data so it's older than v0.10, there may be deserialization issues."
    end
end

## Utility functions for ZonedDateTime

# Regex parser that converts clauses like ":59.82-" to well formatted ":59.820-"
function _fixSubseconds(a)
    length(a) == 4 && return a[1:3]*".000"*a[4]
    frac = a[5:length(a)-1]
    frac = length(frac) > 3 ? frac[1:3] : frac*'0'^(3-length(frac))
    return a[1:4]*frac*a[length(a)]
end

function getStandardZDTString(stringTimestamp::String)
    # Additional check+fix for the ultraweird "2020-08-12T12:00Z"
    ts = replace(stringTimestamp, r"T(\d\d):(\d\d)(Z|z|\+|-)" => s"T\1:\2:00.000\3")

    # This is finding :59Z or :59.82-05:00 and fixing it to always have 3 subsecond digits.
    # Temporary fix until TimeZones.jl gets an upstream fix.
    return replace(ts, r":\d\d(\.\d+)?(Z|z|\+|-)" => _fixSubseconds)
end

# Corrects any `::ZonedDateTime` fields of T in corresponding `interm::Dict` as `dateformat"yyyy-mm-ddTHH:MM:SS.ssszzz"`
function standardizeZDTStrings!(T, interm::Dict)
    for (name, typ) in zip(fieldnames(T), T.types)
        if typ <: ZonedDateTime
            namestr = string(name)
            interm[namestr] = getStandardZDTString(interm[namestr])
        end
    end
    nothing
end

function string2ZonedDateTime(stringTimestamp)
    #   ss = split(stringTimestamp, r"(T[0-9.:]*?\K(?=[-+Zz]))|[\[\]]")
  ss = split(stringTimestamp, r"T[\d.:]{5,12}?\K(?=[-+Zz])")
  length(ss) != 2 && error("Misformed zoned timestamp string $stringTimestamp")
  ZonedDateTime(DateTime(ss[1]), TimeZone(ss[2]))
end

# variableType module.type string functions
function typeModuleName(variableType::InferenceVariable)
    io = IOBuffer()
    ioc = IOContext(io, :module=>DistributedFactorGraphs)
    show(ioc, typeof(variableType))
    return String(take!(io))
end

typeModuleName(varT::Type{<:InferenceVariable}) = typeModuleName(varT())

function getTypeFromSerializationModule(variableTypeString::String)
    try
        # split the type at last `.`
        split_st = split(variableTypeString, r"\.(?!.*\.)")
        #if module is specified look for the module in main, otherwise use Main
        if length(split_st) == 2
            m = getfield(Main, Symbol(split_st[1]))
        else
            m = Main
        end
        noparams = split(split_st[end], r"{")
        ret = if 1 < length(noparams)
            # fix #671, but does not work with specific module yet
            bidx = findfirst(r"{", split_st[end])[1]
            Core.eval(m, Base.Meta.parse("$(noparams[1])$(split_st[end][bidx:end])"))
            # eval(Base.Meta.parse("Main.$(noparams[1])$(split_st[end][bidx:end])"))
        else
            getfield(m, Symbol(split_st[end]))
        end

        return ret

    catch ex
        @error "Unable to deserialize soft type $(variableTypeString)"
        io = IOBuffer()
        showerror(io, ex, catch_backtrace())
        err = String(take!(io))
        @error(err)
    end
    nothing
end

##==============================================================================
## Variable Packing and unpacking
##==============================================================================
function packVariable(dfg::AbstractDFG, v::DFGVariable)
    props = Dict{String, Any}()
    props["label"] = string(v.label)
    props["timestamp"] = Dates.format(v.timestamp, "yyyy-mm-ddTHH:MM:SS.ssszzz")
    props["nstime"] = string(v.nstime.value)
    props["tags"] = JSON2.write(v.tags)
    props["ppeDict"] = JSON2.write(v.ppeDict)
    props["solverDataDict"] = JSON2.write(Dict(keys(v.solverDataDict) .=> map(vnd -> packVariableNodeData(dfg, vnd), values(v.solverDataDict))))
    props["smallData"] = JSON2.write(v.smallData)
    props["solvable"] = v.solvable
    props["variableType"] = typeModuleName(getVariableType(v))
    props["dataEntry"] = JSON2.write(Dict(keys(v.dataDict) .=> map(bde -> JSON.json(bde), values(v.dataDict))))
    props["dataEntryType"] = JSON2.write(Dict(keys(v.dataDict) .=> map(bde -> typeof(bde), values(v.dataDict))))
    props["_version"] = _getDFGVersion()
    return props::Dict{String, Any}
end

# returns a DFGVariable
function unpackVariable(dfg::G,
                        packedProps::Dict{String, Any};
                        unpackPPEs::Bool=true,
                        unpackSolverData::Bool=true,
                        unpackBigData::Bool=true) where G <: AbstractDFG
    @debug "Unpacking variable:\r\n$packedProps"
    # Version checking.
    _versionCheck(packedProps)
    label = Symbol(packedProps["label"])
    # Make sure that the timestamp is correctly formatted with subseconds
    packedProps["timestamp"] = getStandardZDTString(packedProps["timestamp"])
    # Parse it
    timestamp = ZonedDateTime(packedProps["timestamp"])
    nstime = Nanosecond(get(packedProps, "nstime", 0))
    # Supporting string serialization using packVariable and CGDFG serialization (Vector{String})
    if packedProps["tags"] isa String
        tags = JSON2.read(packedProps["tags"], Vector{Symbol})
    else
        tags = Symbol.(packedProps["tags"])
    end
    ppeDict = unpackPPEs ? JSON2.read(packedProps["ppeDict"], Dict{Symbol, MeanMaxPPE}) : Dict{Symbol, MeanMaxPPE}()
    smallData = JSON2.read(packedProps["smallData"], Dict{Symbol, SmallDataTypes})

    variableTypeString = if haskey(packedProps, "softtype")
        # TODO Deprecate, remove in v0.12
        @warn "Packed field `softtype` is deprecated and replaced with `variableType`"
        packedProps["softtype"]
    else
        packedProps["variableType"]
    end

    variableType = getTypeFromSerializationModule(variableTypeString)
    isnothing(variableType) && error("Cannot deserialize variableType '$variableTypeString' in variable '$label'")
    pointType = getPointType(variableType)

    if unpackSolverData
        packed = JSON2.read(packedProps["solverDataDict"], Dict{String, PackedVariableNodeData})
        solverData = Dict{Symbol, VariableNodeData{variableType, pointType}}(Symbol.(keys(packed)) .=> map(p -> unpackVariableNodeData(dfg, p), values(packed)))
    else
        solverData = Dict{Symbol, VariableNodeData{variableType, pointType}}()
    end
    # Rebuild DFGVariable using the first solver variableType in solverData
    # @info "dbg Serialization 171" variableType Symbol(packedProps["label"]) timestamp nstime ppeDict solverData smallData Dict{Symbol,AbstractDataEntry}() Ref(packedProps["solvable"])
    # variable = DFGVariable{variableType}(Symbol(packedProps["label"]), timestamp, nstime, Set(tags), ppeDict, solverData,  smallData, Dict{Symbol,AbstractDataEntry}(), Ref(packedProps["solvable"]))
    variable = DFGVariable( Symbol(packedProps["label"]),
                            variableType,
                            timestamp=timestamp,
                            nstime=nstime,
                            tags=Set{Symbol}(tags),
                            estimateDict=ppeDict,
                            solverDataDict=solverData,
                            smallData=smallData,
                            dataDict=Dict{Symbol,AbstractDataEntry}(),
                            solvable=packedProps["solvable"] )
    #

    # Now rehydrate complete DataEntry type.
    if unpackBigData
        #TODO Deprecate - for backward compatibility between v0.8 and v0.9, remove in v0.10
        if haskey(packedProps, "bigDataElemType")
            @warn "`bigDataElemType` is deprecate, please save data again with new version that uses `dataEntryType`"
            dataElemTypes = JSON2.read(packedProps["bigDataElemType"], Dict{Symbol, Symbol})
        else
            dataElemTypes = JSON2.read(packedProps["dataEntryType"], Dict{Symbol, Symbol})
            for (k,name) in dataElemTypes
                dataElemTypes[k] = Symbol(split(string(name), '.')[end])
            end
        end

        #TODO Deprecate - for backward compatibility between v0.8 and v0.9, remove in v0.10
        if haskey(packedProps, "bigData")
            @warn "`bigData` is deprecate, please save data again with new version"
            dataIntermed = JSON2.read(packedProps["bigData"], Dict{Symbol, String})
        else
            dataIntermed = JSON2.read(packedProps["dataEntry"], Dict{Symbol, String})
        end

        for (k,bdeInter) in dataIntermed
            interm = JSON.parse(bdeInter)
            objType = getfield(DistributedFactorGraphs, dataElemTypes[k])
            standardizeZDTStrings!(objType, interm)
            fullVal = Unmarshal.unmarshal(objType, interm)
            variable.dataDict[k] = fullVal
        end
    end

    return variable
end

# returns a PackedVariableNodeData
function packVariableNodeData(::G, d::VariableNodeData{T}) where {G <: AbstractDFG, T <: InferenceVariable}
  @debug "Dispatching conversion variable -> packed variable for type $(string(d.variableType))"
  # TODO change to Vector{Vector{Float64}} which can be directly packed by JSON
  castval = if 0 < length(d.val)
    precast = getCoordinates.(T, d.val)
    @cast castval[i,j] := precast[j][i]
    castval
  else
    zeros(1,0)
  end
  _val = castval[:]
#   castbw = if 0 < length(d.bw)
#     @cast castbw[i,j] := d.bw[j][i]
#     castbw
#   else
#     zeros(1,0)
#   end
#   _bw = castbw[:]
  return PackedVariableNodeData(_val, size(castval,1),
                                d.bw[:], size(d.bw,1),
                                d.BayesNetOutVertIDs,
                                d.dimIDs, d.dims, d.eliminated,
                                d.BayesNetVertID, d.separator,
                                typeModuleName(d.variableType),
                                d.initialized,
                                d.infoPerCoord,
                                d.ismargin,
                                d.dontmargin,
                                d.solveInProgress,
                                d.solvedCount,
                                d.solveKey)
end

function unpackVariableNodeData(dfg::G, d::PackedVariableNodeData) where G <: AbstractDFG
    @debug "Dispatching conversion packed variable -> variable for type $(string(d.variableType))"
    # Figuring out the variableType
    # TODO deprecated remove in v0.11 - for backward compatibility for saved variableTypes.
    ststring = string(split(d.variableType, "(")[1])
    T = getTypeFromSerializationModule(ststring)
    isnothing(T) && error("The variable doesn't seem to have a variableType. It needs to set up with an InferenceVariable from IIF. This will happen if you use DFG to add serialized variables directly and try use them. Please use IncrementalInference.addVariable().")

    r3 = d.dimval
    c3 = r3 > 0 ? floor(Int,length(d.vecval)/r3) : 0
    M3 = reshape(d.vecval,r3,c3)
    @cast val_[j][i] := M3[i,j]
    vals = Vector{getPointType(T)}(undef, length(val_))
    # vals = getPoint.(T, val_)
    for (i,v) in enumerate(val_)
      vals[i] = getPoint(T, v)
    end

    r4 = d.dimbw
    c4 = r4 > 0 ? floor(Int,length(d.vecbw)/r4) : 0
    BW = reshape(d.vecbw,r4,c4)

    #
    return VariableNodeData{T, getPointType(T)}(vals, BW, d.BayesNetOutVertIDs,
        d.dimIDs, d.dims, d.eliminated, d.BayesNetVertID, d.separator,
        T(), d.initialized, d.infoPerCoord, d.ismargin, d.dontmargin,
        d.solveInProgress, d.solvedCount, d.solveKey,
        Dict{Symbol,Threads.Condition}() )
end

##==============================================================================
## Factor Packing and unpacking
##==============================================================================

function packFactor(dfg::G, f::DFGFactor)::Dict{String, Any} where G <: AbstractDFG
    # Construct the properties to save
    props = Dict{String, Any}()
    props["label"] = string(f.label)
    props["timestamp"] = Dates.format(f.timestamp, "yyyy-mm-ddTHH:MM:SS.ssszzz")
    props["nstime"] = string(f.nstime.value)
    props["tags"] = JSON2.write(f.tags)
    # Pack the node data
    fnctype = getSolverData(f).fnc.usrfnc!
    try
        packtype = convertPackedType(fnctype)
        packed = convert(PackedFunctionNodeData{packtype}, getSolverData(f))
        props["data"] = JSON2.write(packed)
    catch ex
        io = IOBuffer()
        showerror(io, ex, catch_backtrace())
        err = String(take!(io))
        msg = "Error while packing '$(f.label)' as '$fnctype', please check the unpacking/packing converters for this factor - \r\n$err"
        error(msg)
    end
    # Include the type
    props["fnctype"] = String(_getname(fnctype))
    props["_variableOrderSymbols"] = JSON2.write(f._variableOrderSymbols)
    props["solvable"] = getSolvable(f)
    props["_version"] = _getDFGVersion()
    return props
end

function decodePackedType(::Type{T}, packeddata::GenericFunctionNodeData{PT}) where {T<:FactorOperationalMemory, PT}
    # usrtyp = convert(FunctorInferenceType, packeddata.fnc)
    # Also look at parentmodule
    usrtyp = convertStructType(PT)
    fulltype = DFG.FunctionNodeData{T{usrtyp}}
    factordata = convert(fulltype, packeddata)
    return factordata
end

function unpackFactor(dfg::G, packedProps::Dict{String, Any})::DFGFactor where G <: AbstractDFG
    # Version checking.
    _versionCheck(packedProps)

    label = packedProps["label"]
    # Make sure that the timestamp is correctly formatted with subseconds
    packedProps["timestamp"] = getStandardZDTString(packedProps["timestamp"])
    # Parse it
    timestamp = ZonedDateTime(packedProps["timestamp"])
    nstime = Nanosecond(get(packedProps, "nstime", 0))
    if packedProps["tags"] isa String
        tags = JSON2.read(packedProps["tags"], Vector{Symbol})
    else
        tags = Symbol.(packedProps["tags"])
        # If tags is empty we need to make sure it's a Vector{Symbol}
        if length(tags) == 0
            tags = Vector{Symbol}()
        end
    end
    data = packedProps["data"]
    datatype = packedProps["fnctype"]
    @debug "DECODING factor type = '$(datatype)' for factor '$label'"
    packtype = getTypeFromSerializationModule(dfg, Symbol("Packed"*datatype))

    # FIXME type instability from nothing to T
    packed = nothing
    fullFactorData = nothing
    try
        packed = JSON2.read(data, GenericFunctionNodeData{packtype})
        decodeType = getFactorOperationalMemoryType(dfg)
        fullFactorData = decodePackedType(decodeType, packed)
    catch ex
        io = IOBuffer()
        showerror(io, ex, catch_backtrace())
        err = String(take!(io))
        msg = "Error while unpacking '$label' as '$datatype', please check the unpacking/packing converters for this factor - \r\n$err"
        error(msg)
    end

    # Include the type
    if packedProps["_variableOrderSymbols"] isa String
        _variableOrderSymbols = JSON2.read(packedProps["_variableOrderSymbols"], Vector{Symbol})
    else
        _variableOrderSymbols = Symbol.(packedProps["_variableOrderSymbols"])
    end
    solvable = packedProps["solvable"]

    # Rebuild DFGFactor
    #TODO use constuctor to create factor
    factor = DFGFactor( Symbol(label),
                        timestamp,
                        nstime,
                        Set(tags),
                        fullFactorData,
                        solvable,
                        Tuple(_variableOrderSymbols))
    #

    # Note, once inserted, you still need to call rebuildFactorMetadata!
    return factor
end

##==============================================================================
## Serialization
##==============================================================================

"""
  $(SIGNATURES)
Get a type from the serialization module inside DFG.
"""
function getTypeFromSerializationModule(dfg::G, moduleType::Symbol) where G <: AbstractDFG
    st = nothing
    try
        st = getfield(Main, Symbol(moduleType))
    catch ex
        @error "Unable to deserialize soft type $(moduleType)"
        io = IOBuffer()
        showerror(io, ex, catch_backtrace())
        err = String(take!(io))
        @error(err)
    end
    return st
end


text:
# copied from https://github.com/JuliaStats/StatsModels.jl/blob/0e59b8ed8c6325175ef9425e89ecab04d87bf70f/src/statsmodel.jl#L1-L33

##############################################################################
#
# A macro for doing delegation
#
# This macro call
#
#     @delegate MyContainer.elems [:size, :length, :ndims, :endof]
#
# produces this block of expressions
#
#     size(a::MyContainer) = size(a.elems)
#     length(a::MyContainer) = length(a.elems)
#     ndims(a::MyContainer) = ndims(a.elems)
#     endof(a::MyContainer) = endof(a.elems)
#
##############################################################################

macro delegate(source, targets)
    typename = esc(source.args[1])
    sourceargs2 = source.args[2]
    fieldname = Expr(:quote, isa(sourceargs2, QuoteNode) ? sourceargs2.value : sourceargs2.args[1])
    funcnames = targets.args
    n = length(funcnames)
    result = quote begin end end
    for i in 1:n
        funcname = esc(funcnames[i])
        f = quote
            ($funcname)(a::($typename), args...; kwargs...) = ($funcname)(getfield(a, $fieldname), args...; kwargs...)
        end
        push!(result.args[2].args, f)
    end
    return result
end


text:
"""
    Assembly{TF, TP<:AbstractVector{<:AbstractVector{TF}},
        TC<:AbstractVector{<:Integer}, TE<:AbstractVector{Element{TF}}}

Composite type that defines an assembly of connected nonlinear beam elements for
analysis.

# Fields
 - `points`: Array of all beam element endpoints
 - `start`: Array containing point index where each beam element starts
 - `stop`: Array containing point index where each beam element stops
 - `elements`: Array of `Element`s
"""
struct Assembly{TF, TP<:AbstractVector{<:AbstractVector{TF}}, TC<:AbstractVector{<:Integer}, TE<:AbstractVector{Element{TF}}}
    points::TP
    start::TC
    stop::TC
    elements::TE
end
Base.eltype(::Assembly{TF, TP, TC, TE}) where {TF, TP, TC, TE} = TF

"""
    Assembly(points, start, stop; kwargs...)

Construct an assembly of connected nonlinear beam elements for analysis.  Beam lengths
and midpoints may be manually specified in case beam elements are curved rather than
straight.

# Arguments
 - `points`: Array of all beam element endpoints
 - `start`: Array containing point indices where each beam element starts
 - `stop`: Array containing point indices where each beam element stops

# Keyword Arguments
 - `stiffness`: Array of (6 x 6) stiffness matrices for each beam element,
        acts as an alternative to providing `compliance`
 - `compliance`: Array of (6 x 6) compliance matrices for each beam element,
        defaults to `zeros(6,6)` for each beam element
 - `mass`: Array of (6 x 6) mass matrices for each beam element,
        defaults to `zeros(6,6)` for each beam element
 - `frames`: Array of (3 x 3) tranformation matrices for each beam element.
        Transforms from the local undeformed beam frame to the global frame) and defaults
        to the identity matrix for each beam element
 - `lengths`: Array containing the length of each beam, defaults to the distance between
        beam endpoints
 - `midpoints`: Array containing the midpoint of each beam element, defaults to the average
        of the beam element endpoints
"""
function Assembly(points, start, stop;
    stiffness = nothing,
    compliance = nothing,
    mass = nothing,
    frames = nothing,
    lengths = norm.(points[stop] - points[start]),
    midpoints = (points[stop] + points[start])/2)

    nelem = length(start)

    if isnothing(compliance)
        if isnothing(stiffness)
            compliance = fill((@SMatrix zeros(6,6)), nelem)
        else
            compliance = [(@MMatrix zeros(eltype(eltype(stiffness)), 6,6)) for i=1:nelem] #can't use fill because it copies the reference. Need a different value for every i.
            for i = 1:nelem
                filled_cols = findall(vec(mapslices(col -> any(row -> !isapprox(row, 0), col), stiffness[i], dims = 1)))
                compliance[i][filled_cols,filled_cols] .= inv(Matrix(stiffness[i][filled_cols, filled_cols]))
            end
            compliance = SMatrix.(compliance)
        end
    end

    if isnothing(mass)
        mass = fill((@SMatrix zeros(6,6)), nelem)
    end

    if isnothing(frames)
        frames = fill(I3, nelem)
    end

    elements = Element.(lengths, midpoints, compliance, mass, frames)

    return Assembly(points, promote(start, stop)..., elements)
end

"""
    discretize_beam(L, start, discretization; frame, curvature)

Discretize a beam of length `L` located at `start` according to the discretization provided
in `discretization`

Return the lengths, endpoints, midpoints, and reference frame of the beam elements.

# Arguments
 - `L`: Beam length
 - `start`: Beam starting point
 - `discretization`: Number of beam elements, or the normalized endpoints of each beam
        element, with values ranging from 0 to 1.

# Keyword Arguments
 - `frame`: Reference frame at the start of the beam element, represented by a 3x3
         transformation matrix from the undeformed local frame to the body frame.
 - `curvature`: curvature vector
"""
function discretize_beam(L::Number, start::AbstractVector, discretization::Integer;
    frame=I3, curvature=(@SVector zeros(3)))

    return discretize_beam(L, start, range(0, 1, length=discretization+1);
        frame=frame, curvature=curvature)
end

function discretize_beam(L::Number, start::AbstractVector, discretization;
    frame=I3, curvature=(@SVector zeros(3)))

    r1 = SVector{3}(start)
    Cab = SMatrix{3,3}(frame)
    k = SVector{3}(curvature)

    # discretize beam
    sp = L*discretization
    sm = (sp[1:end-1] .+ sp[2:end])/2
    ΔL = sp[2:end] .- sp[1:end-1]

    # precompute some curvature quantities
    kkt = k*k'
    ktilde = tilde(k)
    kn = sqrt(k'*k)

    if iszero(curvature)
        triads = fill(Cab, length(sm))
        xp = [r1 + s*Cab*e1 for s in sp]
        xm = [r1 + s*Cab*e1 for s in sm]
    else
        triads = curve_triad.(Ref(Cab), Ref(kkt), Ref(ktilde), Ref(kn), sm)
        xp = curve_coordinates.(Ref(r1), Ref(Cab), Ref(kkt), Ref(ktilde), Ref(kn), sp)
        xm = curve_coordinates.(Ref(r1), Ref(Cab), Ref(kkt), Ref(ktilde), Ref(kn), sm)
    end

    return ΔL, xp, xm, triads
end

"""
    discretize_beam(start, stop, discretization; frame, curvature)

Discretize a beam from `start` to `stop` according to the discretization provided in
`discretization`.

Return the lengths, endpoints, midpoints, and reference frame of the beam elements.

# Arguments
 - `start`: Beam starting point
 - `stop`: Beam ending point
 - `discretization`: Number of beam elements, or the normalized endpoints of each beam
        element, with values ranging from 0 to 1.

# Keyword Arguments
 - `frame`: Reference frame at the start of the beam element, represented by a 3x3
        transformation matrix from the undeformed local frame to the body frame.
 - `curvature`: curvature vector
"""
function discretize_beam(start::AbstractVector, stop::AbstractVector,
    discretization::Integer; frame=I3, curvature=(@SVector zeros(3)))

    return discretize_beam(start, stop, range(0, 1, length=discretization+1);
        frame=frame, curvature=curvature)
end

function discretize_beam(start::AbstractVector, stop::AbstractVector, discretization;
    frame=I3, curvature=(@SVector zeros(3)))

    # compute curve length
    L = curve_length(start, stop, curvature)

    # discretize beam
    ΔL, xp, xm, triads = discretize_beam(L, start, discretization; frame=frame, curvature=curvature)

    return ΔL, xp, xm, triads
end

"""
    curve_length(start, stop, curvature)

Calculate the length of a curve given its endpoints and its curvature vector.
"""
function curve_length(start, stop, curvature)

    r1 = SVector{3}(start)
    r2 = SVector{3}(stop)
    k = SVector{3}(curvature)

    # precomputed quantities
    kn = sqrt(k'*k)
    r = r2 - r1
    rn2 = r'*r
    rn = sqrt(rn2)

    # compute curve length
    if iszero(kn) || (iszero(k[2]) && iszero(k[3]))
        # the beam is straight or twisted, but not curved
        ΔL = rn
    elseif iszero(k[1])
        # the beam is curved, but not twisted
        ΔL = 2*asin(kn*rn/2)/kn
    else
        # the beam is curved and twisted
        f = (L) -> (2*(kn^2-k[1]^2)*(1-cos(kn*L)) + k[1]^2*(kn*L)^2)/kn^4 - rn2
        ΔL = Roots.find_zero(f, (rn, rn*kn/k), Roots.Brent())
    end

    return ΔL
end

"""
    curve_triad(Cab, k, s)
    curve_triad(Cab, kkt, ktilde, kn, s)

Return the transformation matrix at `s` along the length of the beam given
the curvature vector `k` and the initial transformation matrix `Cab`.
"""
@inline curve_triad(Cab, k, s) = curve_triad(Cab, k*k', tilde(k), sqrt(k'*k), s)
@inline curve_triad(Cab, kkt, ktilde, kn, s) = SMatrix{3,3}(Cab*((I - kkt/kn^2)*cos(kn*s) +
    ktilde/kn*sin(kn*s) + kkt/kn^2))

"""
    curve_coordiantes(r, Cab, k, s)
    curve_coordinates(r, Cab, kkt, ktilde, kn, s)

Return the coordinates at `s` along the length of the beam given the starting
point `r`, initial transformation matrix `Cab`, and curvature vector `k`.
"""
@inline curve_coordinates(r, Cab, k, s) = curve_coordinates(r, Cab, k*k', tilde(k), sqrt(k'*k), s)
@inline curve_coordinates(r, Cab, kkt, ktilde, kn, s) = r + SVector{3}(Cab*((I/kn - kkt/kn^3)*sin(kn*s) + ktilde/kn^2*(1-cos(kn*s)) + kkt/kn^2*s)*e1)


text:
function test_woodbury()
    len = 500;
    ls_solver_julia = linear_solver_JULIA();
    solver = ls_solver_julia;
    mat = sparse(rand(len,len));
    initialize!(solver, mat, :unsymmetric);

    @test solver._SparseMatrix == mat

    inertia = ls_factor!(solver, len, 0)
    @test inertia == 1

    sol = 1.0*zeros(len);
    rhs = 1.0*ones(len);

    U = rand(len,2);
    #U[1,1] = 2.0;
    #U[2,2] = 2.0;
    V = copy(U)';

    W = woodbury_identity(solver, U, V);

    sol = ls_solve(W,rhs);

    err = norm((mat + U * V) * sol - rhs)
    @show err
    @fact err --> less_than(1e-10)
end

facts("helpers") do
  context("woodbury") do
      test_woodbury();
  end
end

max_forks_repo_name: ohinder/Everest


text:
#####
##### Shared IO functions
#####

export OutputStyle, SingleFile, FolderStructure

abstract type OutputStyle end

struct FolderStructure <: OutputStyle
  root_folder
  filename
end
FolderStructure() = FolderStructure(nothing, nothing)

struct SingleFile <: OutputStyle
  filename
end
SingleFile() = SingleFile(nothing)


text:
"""
    TwoSectionView{T<:Real} <: AbstractGraph{Int64}

Represents a 2-section view of a hypergraph `h`.
Note this is a view - changes to the original hypergraph will be automatically reflected in the view.

**Constructors**

TwoSectionView(::Hypergraph)

The 2-section view of a hypergraph is suitable for processing with the LightGraphs.jl package.
Several LightGraphs methods are provided for the compability.

"""
struct TwoSectionView{T<:Real} <: AbstractGraph{Int}
    h::Hypergraph{T}
end

LightGraphs.nv(t::TwoSectionView) = length(t.h.v2he)

LightGraphs.vertices(t::TwoSectionView) = Base.OneTo(nv(t))

function LightGraphs.ne(t::TwoSectionView)
    s = 0
    for x in t.h.he2v
        s += length(x) * (length(x) - 1)
    end
    div(s, 2)
end

"""
    LightGraphs.all_neighbors(t::TwoSectionView, v::Integer)

Returns N(v) (the vertex v is not included in N(v))
"""
function LightGraphs.all_neighbors(t::TwoSectionView, v::Integer)
    neighbors = Set{Int}()
    for he in keys(t.h.v2he[v])
        union!(neighbors, keys(t.h.he2v[he]))
    end
    delete!(neighbors, v) #remove v from its neighborhood
    collect(neighbors) #returns the corresponding array
end

function LightGraphs.has_edge(t::TwoSectionView, s, d)
    !isempty(intersect(keys(t.h.v2he[s]), keys(t.h.v2he[d])))
end

LightGraphs.outneighbors(t::TwoSectionView, v::Integer) =
    LightGraphs.all_neighbors(t::TwoSectionView, v)

LightGraphs.inneighbors(t::TwoSectionView, v::Integer) =
    LightGraphs.all_neighbors(t::TwoSectionView, v)

"""
    LightGraphs.SimpleGraph(t::TwoSectionView)

Creates a `LightGraphs.SimpleGraph` representation of a `TwoSectionView` t.

This creates a copy of the date. Note that the weights information is not stored
in the created `SimpleGraph`.
"""
function LightGraphs.SimpleGraph(t::TwoSectionView)
    g = SimpleGraph(nv(t))
    for v in LightGraphs.vertices(t)
        neighbors_v = LightGraphs.all_neighbors(t, v)
        for neighbor in neighbors_v
            add_edge!(g, v, neighbor)
        end
    end
    g
end

LightGraphs.is_directed(t::TwoSectionView) = false

"""
    shortest_path(t::TwoSectionView,source::Int, target::Int)

Finds a single shortest path in a graph `b` between vertices
`source` and `target`.
Note that if several paths of the same length exist, only one
will be returned.

"""
function shortest_path(t::TwoSectionView, source::Int, target::Int)
    checkbounds(t.h.v2he, source)
    checkbounds(t.h.v2he, target)
    dj = dijkstra_shortest_paths(t, source)
    enumerate_paths(dj)[target]
end


text:
@testset "Network" begin
  m = Network("Test")

  # Adding a Reaction with an unknown Species should throw an error
  @test_throws ErrorException m <= Reaction("bad",      1.0, "X --> 0")
  @test_throws ErrorException m <= Reaction("also bad", 1.0, "0 --> X")
  @test number_reactions(m) == 0

  # Adding a Species with negative counts should throw an error
  @test_throws ErrorException m <= Species("X", -1)

  # Adding a Species with non-negative counts is okay
  m <= Species("X", 0)
  m <= Species("Y", 100)

  @test number_species(m) == 2

  # Adding a Reaction with an unknown Species should throw an error
  @test_throws ErrorException m <= Reaction("should not work",  1.0, "X --> Z")
  @test_throws ErrorException m <= Reaction("should also fail", 1.0, "Z --> Y")
  @test number_reactions(m) == 0

  # Adding the same Species twice should change its initial value
  m <= Species("X", 50)

  @test number_species(m) == 2
  @test get_species(m, "X").population == 50
end


text:
module BEDFunc

export bedtools_intersect
export convert_methpipebed_to_bedfile
export write_bedfile
include(Pkg.dir("BEDFunc","src","bedtools_intersect.jl"))

export load_dataframe_from_url
export run_cmd
include(Pkg.dir("BEDFunc","src","misc.jl"))

export query
export BEDQuery
export bgzip_and_index
include(Pkg.dir("BEDFunc","src","tabix.jl"))

export parse_gtf_gene_features_to_df
include(Pkg.dir("BEDFunc","src","gtf.jl"))
end # module


text:
"""
        edas(decisionMat, weights, fns)

Apply EDAS (Evaluation based on Distance from Average Solution) for a given matrix and weights.

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n alterntives and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fns::Array{Function, 1}`: m-vector of functions to be applied on the columns.

# Description
edas() applies the EDAS method to rank n alterntives subject to m criteria which are supposed to be
either maximized or minimized.

# Output
- `::EDASResult`: EDASResult object that holds multiple outputs including scores, rankings, and best index.

# Examples
```julia-repl
julia> decmat = [5000 5 5300 450;
       4500 5 5000 400;
       4500 4 4700 400;
       4000 4 4200 400;
       5000 4 7100 500;
       5000 5 5400 450;
       5500 5 6200 500;
       5000 4 5800 450]
8×4 Array{Int64,2}:
 5000  5  5300  450
 4500  5  5000  400
 4500  4  4700  400
 4000  4  4200  400
 5000  4  7100  500
 5000  5  5400  450
 5500  5  6200  500
 5000  4  5800  450

julia> df = makeDecisionMatrix(decmat)
8×4 DataFrame
 Row │ Crt1     Crt2     Crt3     Crt4
     │ Float64  Float64  Float64  Float64
─────┼────────────────────────────────────
   1 │  5000.0      5.0   5300.0    450.0
   2 │  4500.0      5.0   5000.0    400.0
   3 │  4500.0      4.0   4700.0    400.0
   4 │  4000.0      4.0   4200.0    400.0
   5 │  5000.0      4.0   7100.0    500.0
   6 │  5000.0      5.0   5400.0    450.0
   7 │  5500.0      5.0   6200.0    500.0
   8 │  5000.0      4.0   5800.0    450.0

julia> weights = [0.25, 0.25, 0.25, 0.25];

julia> fns = [maximum, maximum, minimum, minimum];

julia> result = edas(df, weights, fns);

julia> result.scores
8-element Array{Float64,1}:
 0.7595941163602383
 0.8860162461306114
 0.6974721951442592
 0.739657763190231
 0.05908329207449442
 0.7318326305342363
 0.6416913873322523
 0.38519414045559647

julia> result.bestIndex
2
```

# References

Keshavarz Ghorabaee, M., Zavadskas, E. K., Olfat, L., & Turskis, Z. (2015). Multi-criteria inventory classification using a new method of evaluation based on distance from average solution (EDAS). Informatica, 26(3), 435-451.

Ulutaş, A. (2017). EDAS Yöntemi Kullanılarak Bir Tekstil Atölyesi İçin Dikiş Makinesi Seçimi. İşletme Araştırmaları Dergisi, 9(2), 169-183.

"""
function edas(decisionMat::DataFrame, weights::Array{Float64,1}, fns::Array{Function,1})::EDASResult

    # df = convert(Matrix, decisionMat)
    df = Matrix(decisionMat)

    row, col = size(df)

    w = unitize(weights)

        PDAMatrix = similar(df)
    NDAMatrix = similar(df)

    AV = zeros(Float64, col)

    @inbounds for i in 1:col
        AV[i] = mean(df[:, i])
        for j in 1:row
            if fns[i] == maximum
                PDAMatrix[j, i] = max(0, df[j, i] - AV[i]) ./ AV[i]
                NDAMatrix[j, i] = max(0, AV[i] - df[j, i]) ./ AV[i]
            elseif fns[i] == minimum
                PDAMatrix[j, i] = max(0, AV[i] - df[j, i]) ./ AV[i]
                NDAMatrix[j, i] = max(0, df[j, i] - AV[i]) ./ AV[i]
            end
        end
    end

    SP = zeros(Float64, row)
    SN = zeros(Float64, row)

    for i in 1:row
        SP[i] = w .* PDAMatrix[i, :] |> sum
        SN[i] = w .* NDAMatrix[i, :] |> sum
    end

    NSP = SP ./ maximum(SP)
    NSN = 1 .- SN ./ maximum(SN)

    scores = (NSP .+ NSN) ./ 2

    rankings = sortperm(scores)

    bestIndex = rankings |> last

    result = EDASResult(
        decisionMat,
        w,
        scores,
        rankings,
        bestIndex
    )

    return result
end

max_stars_repo_name: bahadirfyildirim/JMcDM

max_issues_repo_name: bahadirfyildirim/JMcDM

max_forks_repo_name: bahadirfyildirim/JMcDM


text:
function _add_dual_vars_in_dual_cones(
    dual_obj_affine_terms::Dict{VI,T},
    dual_model::MOI.ModelLike,
    primal_model::MOI.ModelLike,
    primal_dual_map::PrimalDualMap{T},
    dual_names::DualNames,
    ::Type{F},
    ::Type{S},
) where {T,F,S}
    for ci in MOI.get(primal_model, MOI.ListOfConstraintIndices{F,S}()) # Constraints of type {F, S}
        # If `F` not one of these two, we can skip the `in` check.
        if (F === MOI.VectorOfVariables || F === MOI.VariableIndex) &&
           haskey(primal_dual_map.constrained_var_dual, ci)
            continue
        end
        # Add dual variable to dual cone
        # Fill a dual objective dictionary
        # Fill the primal_con_dual_var dictionary
        ci_dual = add_dual_variable(
            dual_model,
            primal_model,
            dual_names,
            primal_dual_map.primal_con_dual_var,
            dual_obj_affine_terms,
            ci,
        )
        push_to_primal_con_dual_con!(
            primal_dual_map.primal_con_dual_con,
            ci,
            ci_dual,
        )
        push_to_primal_con_constants!(
            primal_model,
            primal_dual_map.primal_con_constants,
            ci,
        )
    end
    return
end

function add_dual_vars_in_dual_cones(
    dual_model::MOI.ModelLike,
    primal_model::MOI.ModelLike,
    primal_dual_map::PrimalDualMap{T},
    dual_names::DualNames,
    con_types::Vector{Tuple{Type,Type}},
) where {T}
    dual_obj_affine_terms = Dict{VI,T}()
    for (F, S) in con_types
        _add_dual_vars_in_dual_cones(
            dual_obj_affine_terms,
            dual_model,
            primal_model,
            primal_dual_map,
            dual_names,
            F,
            S,
        )
    end
    return dual_obj_affine_terms
end

# Utils for primal_con_constants dict
function push_to_primal_con_constants!(
    primal_model::MOI.ModelLike,
    primal_con_constants::Dict{CI,Vector{T}},
    ci::CI{F,S},
) where {T,F<:MOI.AbstractScalarFunction,S<:MOI.AbstractScalarSet}
    push!(primal_con_constants, ci => get_scalar_term(primal_model, ci))
    return
end

function push_to_primal_con_constants!(
    primal_model::MOI.ModelLike,
    primal_con_constants::Dict{CI,Vector{T}},
    ci::CI{F,S},
) where {T,F<:MOI.AbstractVectorFunction,S<:MOI.AbstractVectorSet}
    # No constants need to be passed to the DualOptimizer in this case,
    # because the VectorSet's do not have constants inside them that
    # will need to be used in result queries as some ScalarSet's might.
    # Hence, don't need to push zero to the dict.
    return
end

# Utils for primal_con_dual_con dict
function push_to_primal_con_dual_con!(
    primal_con_dual_con::Dict{CI,CI},
    ci::CI,
    ci_dual::CI,
)
    push!(primal_con_dual_con, ci => ci_dual)
    return
end

function push_to_primal_con_dual_con!(
    primal_con_dual_con::Dict{CI,CI},
    ci::CI,
    ci_dual::Nothing,
)
    return # Don't put in the dict a nothing value
end

# Utils for dual_obj_affine_terms dict
function push_to_dual_obj_aff_terms!(
    primal_model::MOI.ModelLike,
    dual_obj_affine_terms::Dict{VI,T},
    vi::VI,
    func::MOI.AbstractFunction,
    set::MOI.AbstractSet,
    i::Int,
) where {T}
    value = set_dot(i, set, T) * get_scalar_term(func, set, i)
    if !iszero(value) # If value is different than 0 add to the dictionary
        push!(dual_obj_affine_terms, vi => value)
    end
    return
end

function push_to_dual_obj_aff_terms!(
    ::MOI.ModelLike,
    ::Dict{VI},
    ::VI,
    ::MOI.VectorOfVariables,
    ::MOI.AbstractVectorSet,
    i::Int,
)
    return # It is zero so don't push to the dual_obj_affine_terms
end

function add_dual_variable(
    dual_model::MOI.ModelLike,
    primal_model::MOI.ModelLike,
    dual_names::DualNames,
    primal_con_dual_var::Dict{CI,Vector{VI}},
    dual_obj_affine_terms::Dict{VI,T},
    ci::CI{F,S},
) where {T,F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
    vis, con_index = add_dual_cone_constraint(dual_model, primal_model, ci)
    push!(primal_con_dual_var, ci => vis) # Add the map of the added dual variable to the relationated constraint
    # Get constraint name
    ci_name = MOI.get(primal_model, MOI.ConstraintName(), ci)
    # Add each vi to the dictionary
    func = get_function(primal_model, ci)
    set = get_set(primal_model, ci)
    for (i, vi) in enumerate(vis)
        push_to_dual_obj_aff_terms!(
            primal_model,
            dual_obj_affine_terms,
            vi,
            func,
            set,
            i,
        )
        if !is_empty(dual_names)
            set_dual_variable_name(
                dual_model,
                vi,
                i,
                ci_name,
                dual_names.dual_variable_name_prefix,
            )
        end
    end
    return con_index
end

function set_dual_variable_name(
    dual_model::MOI.ModelLike,
    vi::VI,
    i::Int,
    ci_name::String,
    prefix::String,
)
    isempty(ci_name) && return
    MOI.set(dual_model, MOI.VariableName(), vi, prefix * ci_name * "_$i")
    return
end

function add_primal_parameter_vars(
    dual_model::MOI.ModelLike,
    primal_model::MOI.ModelLike,
    primal_dual_map::PrimalDualMap{T},
    dual_names::DualNames,
    variable_parameters::Vector{VI},
    primal_objective,
    ignore_objective::Bool,
) where {T}
    # if objective is ignored we only need parameters that appear in the
    # quadratic objective
    if ignore_objective
        # only crossed terms (parameter times primal variable) of the objective
        # are required
        added = Set{VI}()
        for vec in values(primal_objective.quad_cross_parameters)
            for term in vec
                ind = term.variable
                if ind in added
                    # do nothing
                else
                    push!(added, ind)
                    vi = MOI.add_variable(dual_model)
                    push_to_primal_parameter!(
                        primal_dual_map.primal_parameter,
                        ind,
                        vi,
                    )
                    # set name
                    if !is_empty(dual_names)
                        vi_name = MOI.get(primal_model, MOI.VariableName(), ind)
                        set_parameter_variable_name(
                            dual_model,
                            vi,
                            vi_name,
                            dual_names,
                        )
                    end
                end
            end
        end
    elseif length(variable_parameters) > 0
        vis = MOI.add_variables(dual_model, length(variable_parameters))
        for i in eachindex(vis)
            push_to_primal_parameter!(
                primal_dual_map.primal_parameter,
                variable_parameters[i],
                vis[i],
            )
            if !is_empty(dual_names)
                vi_name = MOI.get(
                    primal_model,
                    MOI.VariableName(),
                    variable_parameters[i],
                )
                set_parameter_variable_name(
                    dual_model,
                    vis[i],
                    vi_name,
                    dual_names,
                )
            end
        end
    end
    return
end

# Save mapping between primal parameter and dual parameter
function push_to_primal_parameter!(
    primal_parameter::Dict{VI,VI},
    vi::VI,
    vi_dual::VI,
)
    push!(primal_parameter, vi => vi_dual)
    return
end

# Add name to parameter variable
function set_parameter_variable_name(
    dual_model::MOI.ModelLike,
    vi::VI,
    vi_name::String,
    dual_names,
)
    prefix =
        dual_names.parameter_name_prefix == "" ? "param_" :
        dual_names.parameter_name_prefix
    MOI.set(dual_model, MOI.VariableName(), vi, prefix * vi_name)
    return
end

function add_quadratic_slack_vars(
    dual_model::MOI.ModelLike,
    primal_model::MOI.ModelLike,
    primal_dual_map::PrimalDualMap{T},
    dual_names::DualNames,
    primal_objective,
) where {T}
    # only main quadratic terms (primal variable times primal variable)
    # of the objective are required
    added = Set{VI}()
    for term in primal_objective.obj.quadratic_terms
        for ind in [term.variable_1, term.variable_2]
            if ind in added
                #do nothing
            else
                push!(added, ind)
                vi = MOI.add_variable(dual_model)
                push_to_quad_slack!(
                    primal_dual_map.primal_var_dual_quad_slack,
                    ind,
                    vi,
                )
                # set name
                if !is_empty(dual_names)
                    vi_name = MOI.get(primal_model, MOI.VariableName(), ind)
                    set_quad_slack_name(dual_model, vi, vi_name, dual_names)
                end
            end
        end
    end
    return
end

# Save mapping between primal variable and dual quadratic slack
function push_to_quad_slack!(dual_quad_slack::Dict{VI,VI}, vi::VI, vi_dual::VI)
    push!(dual_quad_slack, vi => vi_dual)
    return
end

# set name for dual quadratic slack
function set_quad_slack_name(
    dual_model::MOI.ModelLike,
    vi::VI,
    vi_name::String,
    dual_names,
)
    prefix = if dual_names.quadratic_slack_name_prefix == ""
        "quadslack_"
    else
        dual_names.quadratic_slack_name_prefix
    end
    MOI.set(dual_model, MOI.VariableName(), vi, prefix * vi_name)
    return
end


text:
@sumProductRule(:node_type     => Gamma,
                :outbound_type => Message{Gamma},
                :inbound_types => (Nothing, Message{PointMass}, Message{PointMass}),
                :name          => SPGammaOutNPP)

@naiveVariationalRule(:node_type     => Gamma,
                      :outbound_type => Message{Gamma},
                      :inbound_types => (Nothing, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VBGammaOut)

@naiveVariationalRule(:node_type     => Gamma,
                      :outbound_type => Message{Function},
                      :inbound_types => (ProbabilityDistribution, Nothing, ProbabilityDistribution),
                      :name          => VBGammaA)

@naiveVariationalRule(:node_type     => Gamma,
                      :outbound_type => Message{Gamma},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, Nothing),
                      :name          => VBGammaB)