text:
# # Knapsack Problem - Single Expansion
# This first example demonstrates how a basic JuDGE model can be set up, with
# a 7-node tree, and a single expansion.
using JuMP, JuDGE, GLPK

# First we will define our tree with a depth of 2 and a degree of 2.
mytree = narytree(2,2)

# ## Data

# Now we specify data for each node of the tree. Here we define the data using
# functions, but you could use dictionaries, as shown in subsequent examples.
function invest_cost(node)
   if node == get_node(mytree,[1])
	  180.0
   elseif node == get_node(mytree,[1,1])
	  50.0
   elseif node == get_node(mytree,[1,2])
	  60.0
   elseif node == get_node(mytree,[1,1,1])
	  40.0
   elseif node == get_node(mytree,[1,1,2])
	  60.0
   elseif node == get_node(mytree,[1,2,1])
	  10.0
   elseif node == get_node(mytree,[1,2,2])
	  10.0
   end
end

function item_volume(node)
   if node == get_node(mytree,[1])
	  [6, 2, 1, 1, 1]
   elseif node == get_node(mytree,[1,1])
	  [8, 2, 2, 2, 1]
   elseif node == get_node(mytree,[1,2])
	  [8, 1, 1, 1, 3]
   elseif node == get_node(mytree,[1,1,1])
	  [4, 4, 3, 1, 2]
   elseif node == get_node(mytree,[1,1,2])
	  [1, 3, 1, 1, 2]
   elseif node == get_node(mytree,[1,2,1])
	  [7, 3, 1, 1, 1]
   elseif node == get_node(mytree,[1,2,2])
	  [2, 5, 2, 1, 2]
   end
end

function item_reward(node)
   if node == get_node(mytree,[1])
	  [60, 20, 10, 15, 10]
   elseif node == get_node(mytree,[1,1])
	  [8, 10, 20, 20, 10]
   elseif node == get_node(mytree,[1,2])
	  [8, 10, 15, 10, 30]
   elseif node == get_node(mytree,[1,1,1])
	  [40, 40, 35, 10, 20]
   elseif node == get_node(mytree,[1,1,2])
	  [15, 35, 15, 15, 20]
   elseif node == get_node(mytree,[1,2,1])
	  [70, 30, 15, 15, 10]
   elseif node == get_node(mytree,[1,2,2])
	  [25, 50, 25, 15, 20]
   end
end

# ## Subproblem definitions

# The JuDGE subproblems are defined through a function that takes a node as its
# single argument, returning a JuMP model. We use GLPK as our optimizer.
JuDGE_SP_Solver = optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 0, "mip_gap" => 0.0)
function sub_problems(node)
   model = Model(JuDGE_SP_Solver)
   @expansion(model, bag)
   @expansioncosts(model, bag*invest_cost(node))
   @variable(model, y[1:5], Bin)
   @constraint(model, BagExtension, sum(y[i]*item_volume(node)[i] for i in 1:5) <= 3 + 4 * bag)
   @objective(model, Min, sum(-item_reward(node)[i] * y[i] for i in 1:5))
   model
end

# Note that there are a few JuDGE-specific macros used to define a JuDGE subproblem

# This defines our single expansion variable
@expansion(model, bag)

# This defines the cost of buying this bag at particular node
@expansioncosts(model, bag*invest_cost(node))

# JuDGE subproblem objective is defined using the standard @objective macro:
@objective(model, Min, sum(-item_reward(node)[i] * y[i] for i in 1:5))

# ## Defining and solving the JuDGE model
# The JuDGEModel is defined based on a tree, probability distribution,
# sub_problems, and an optimizer using an interior point method.
JuDGE_MP_Solver = optimizer_with_attributes((method=GLPK.INTERIOR) -> GLPK.Optimizer(),
							"msg_lev" => 0, "mip_gap" => 0.0)
judy = JuDGEModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_MP_Solver)

# Solve the model
JuDGE.solve(judy)

# ## Displaying the output
# Print the objective, and optimal expansions
println("Objective: "*string(judy.bounds.UB))
JuDGE.print_expansions(judy,onlynonzero=false)

# Re-solve the subproblems and print the objective
println("Re-solved Objective: " * string(resolve_subproblems(judy)))

# ## Deterministic equivalent
# Here we set up and solve the deterministic equivalent, this uses the same structure as a `JuDGEModel`.
# We set up our GLPK solver to show its progress.
JuDGE_DE_Solver = optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 2, "mip_gap" => 0.0)
deteq = DetEqModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_DE_Solver)
JuDGE.solve(deteq)
println("Deterministic Equivalent Objective: " * string(objective_value(deteq.problem)))


text:
include("settings.jl")
include("SGSolver.jl")
include("DLRSolver.jl")
include("plotting.jl")

using NPZ;

close("all");

s = Settings();

###########################
# run conservative solver #
###########################
s.NCons = 2; # number of conserved basis functions
s.iCons = [1 1; 1 2]; # indices of conserved basis funtions
solver = DLRSolver(s);

PhiQuad = Array(solver.basis.PhiQuad');
PhiQuadCons = Array(solver.basis.PhiQuadCons');
PhiQuadW = Array(solver.basis.PhiQuadW');
fXi = 0.25;

X,S,W,uCons = FactorIC(solver);

@time tEnd, Xt,St,Wt,uCons = SolveSplitUnconventionalIntegrator(solver,X,S,W,uCons);

# compute moments
uQ = Xt*St*Wt'*PhiQuad + uCons*PhiQuadCons;
uDLRCons = uQ*Array(solver.basis.PhiQuadWFull')*fXi;

uDLRCons = Array(uDLRCons')

##################################
# run stochastic-Galerkin solver #
##################################
s.NCons = 0;
s.iCons = 0;

solver = Solver(s);

@time tEnd, uSG = Solve(solver);

###########################
# run conservative solver #
###########################
s.NCons = 2; # number of conserved basis functions
s.iCons = [1 1; 1 2]; # indices of conserved basis funtions
solver = DLRSolver(s);

@time tEnd, uDLRConsN = SolveNaiveSplitUnconventionalIntegrator(solver);

uDLRConsN = Array(uDLRConsN')

#########################
##### Plot solution #####
#########################
plotSolution = Plotting(s,solver.basis,solver.q,s.tEnd);

PlotExpectedValue(plotSolution,uSG,uDLRCons,uDLRConsN,"noFilter","Figure1");

ext: jl


text:
# =========================================================================
#    CORE FUNCTIONS AND MECHANISMS
# =========================================================================
function bnb_tolerance(tolerance::Float64,timelimit::Int64)
    feas = Model(get_optimizer())
    #set_silent(feas)

    set_optimizer_attribute(feas, "Cuts", 0)
    set_optimizer_attribute(feas, "Threads", 1)
    set_optimizer_attribute(feas, "MIPGap", tolerance)
    set_optimizer_attribute(feas, "TimeLimit", timelimit)

    @variable(feas, o[i = V(), j = V(), k = K(), t = T()] >= 0, Int)
    @variable(feas, p[i = V(), j = V(), k = K(), t = T()] >= 0, Int)
    @variable(feas, I[i = V(), t = vcat(first(T()) - 1, T())], Int)

    @constraint(feas, [i = V(), t = T()],
        I[i,t-1] - I[i,t] - d(i,t) +
        sum(sum(o[j,i,k,t] for j in V()) - sum(o[i,j,k,t] for j in V()) for k in K()) == 0
    ) #Inventory Balance

    @constraint(feas, [i = V(), t = T()],
        V(i).MIN <= I[i,t] <= V(i).MAX
    ) #Inventory Max Min

    @constraint(feas, [i = V()],
        I[i,first(T()) - 1] == V(i).START
    ) #Inventory Start

    @constraint(feas, [i = V(), j = V(), k = K(), t = T()],
        o[i,j,k,t] <= K(k).Q[i,j] * p[i,j,k,t]
    ) #quantity to trips proportion from ij

    @constraint(feas, [i = V(), j = V(), k = K(), t = T()],
        p[i,j,k,t] <= K(k).lim[i,j]
    ) #limit of trips from ij

    @objective(feas, Min,
        sum(
            K(k).f[i,j] * o[i,j,k,t] + K(k).g[i,j] * p[i,j,k,t]
            for i in V(), j in V(), k in K(), t in T()
        ) +
        sum(
            V(i).h * I[i,t]
            for i in V(), t in T()
        )
    )

    optimize!(feas)

    ori_sol = Vector{delivery}()
    for i in V(), j in V(), k in K(), t in T()
        if value(p[i,j,k,t]) > 0
            new_delivery = delivery(i,j,value(o[i,j,k,t]),value(p[i,j,k,t]),k,t)
            push!(ori_sol, new_delivery)
        end
    end #COLLECT SOLUTION

    return ori_sol
end #ORIGINAL FORMULATION OF IRP

function relaxed()
    feas = Model(get_optimizer())
    set_silent(feas)
    set_optimizer_attribute(feas, "Threads", 1)

    @variable(feas, o[i = V(), j = V(), k = K(), t = T()] >= 0) #RELAXED DELIVERY
    @variable(feas, p[i = V(), j = V(), k = K(), t = T()] >= 0) #RELAXED DELIVERY
    @variable(feas, I[i = V(), t = vcat(first(T()) - 1, T())], Int) #RELAXED INVENTORY

    @constraint(feas, [i = V(), t = T()],
        I[i,t-1] - I[i,t] - d(i,t) +
        sum(sum(o[j,i,k,t] for j in V()) - sum(o[i,j,k,t] for j in V()) for k in K()) == 0
    ) #Inventory Balance

    @constraint(feas, [i = V(), t = T()],
        V(i).MIN <= I[i,t] <= V(i).MAX
    ) #Inventory Max Min

    @constraint(feas, [i = V()],
        I[i,first(T()) - 1] == V(i).START
    ) #Inventory Start

    @constraint(feas, [i = V(), j = V(), k = K(), t = T()],
        o[i,j,k,t] <= K(k).Q[i,j] * p[i,j,k,t]
    ) #quantity to trips proportion from ij

    @constraint(feas, [i = V(), j = V(), k = K(), t = T()],
        p[i,j,k,t] <= K(k).lim[i,j]
    ) #limit of trips from ij

    @objective(feas, Min,
        sum(
            K(k).f[i,j] * o[i,j,k,t] + K(k).g[i,j] * p[i,j,k,t]
            for i in V(), j in V(), k in K(), t in T()
        ) +
        sum(
            V(i).h * I[i,t]
            for i in V(), t in T()
        )
    )

    optimize!(feas)

    relaxed_sol = Vector{delivery}()
    for i in V(), j in V(), k in K(), t in T()
        if value(p[i,j,k,t]) > 0
            new_delivery = delivery(i,j,value(o[i,j,k,t]),value(p[i,j,k,t]),k,t)
            push!(relaxed_sol, new_delivery)
        end
    end #COLLECT SOLUTION

    return relaxed_sol
end #SAME AS ORI BUT NO INTEGER DEFINITION

function cluster()
    Q_cluster = JuMP.Containers.DenseAxisArray{Float64}(undef,V(),V())
    f_cluster = JuMP.Containers.DenseAxisArray{Float64}(undef,V(),V())
    g_cluster = JuMP.Containers.DenseAxisArray{Float64}(undef,V(),V())
    lim_cluster = JuMP.Containers.DenseAxisArray{Float64}(undef,V(),V())

    Q_cluster .= 0
    f_cluster .= 0
    g_cluster .= 0
    lim_cluster .= 0

    for i in V(), j in V()
        active_mode = [k for k in K() if K(k).lim[i,j] > 0]

        if length(active_mode) >= 1
            weight_vec = Dict{String,Float64}()
            total_weight = sum(K(m).lim[i,j] for m in active_mode)
            for m in active_mode
                weight_vec[m] = K(m).lim[i,j]/total_weight
            end #vektor bobot proporsi terhadap total biaya peti

            Q_cluster[i,j] = ceil(sum(weight_vec[m] * K(m).Q[i,j] for m in active_mode))
            f_cluster[i,j] = sum(weight_vec[m] * K(m).f[i,j] for m in active_mode)
            g_cluster[i,j] = sum(weight_vec[m] * K(m).g[i,j] for m in active_mode)
            lim_cluster[i,j] = sum(K(m).lim[i,j] for m in active_mode)
        end #PROSES SEKIRANYA ADA KENDARAAN DI SEGMEN
    end

    return veh(Q_cluster,f_cluster,g_cluster,lim_cluster)
end #OUTPUTS NEW AGGREGATED VEHICLE

function reduced(cluster::veh)
    feas = Model(get_optimizer())
    set_silent(feas)

    set_optimizer_attribute(feas, "Cuts", 0)
    set_optimizer_attribute(feas, "Threads", 1)

    @variable(feas, o[i = V(), j = V(), t = T()] >= 0) #RELAXED DELIVERY
    @variable(feas, p[i = V(), j = V(), t = T()] >= 0) #RELAXED DELIVERY
    @variable(feas, I[i = V(), t = vcat(first(T()) - 1, T())], Int) #RELAXED INVENTORY

    @constraint(feas, [i = V(), t = T()],
        I[i,t-1] - I[i,t] - d(i,t) +
        sum(o[j,i,t] for j in V()) - sum(o[i,j,t] for j in V()) == 0
    ) #Inventory Balance

    @constraint(feas, [i = V(), t = T()],
        V(i).MIN <= I[i,t] <= V(i).MAX
    ) #Inventory Max Min

    @constraint(feas, [i = V()],
        I[i,first(T()) - 1] == V(i).START
    ) #Inventory Start

    @constraint(feas, [i = V(), j = V(), t = T()],
        o[i,j,t] <= cluster.Q[i,j] * p[i,j,t]
    ) #quantity to trips proportion from ij

    @constraint(feas, [i = V(), j = V(), t = T()],
        p[i,j,t] <= cluster.lim[i,j]
    ) #limit of trips from ij

    @objective(feas, Min,
        sum(
            cluster.f[i,j] * o[i,j,t] + cluster.g[i,j] * p[i,j,t]
            for i in V(), j in V(), t in T()
        ) +
        sum(
            V(i).h * I[i,t]
            for i in V(), t in T()
        )
    )

    optimize!(feas)

    deliveries = Vector{delivery}()
    for i in V(), j in V(), t in T()
        if value(p[i,j,t]) > 0
            dis = Model(get_optimizer())
            set_silent(dis)

            set_optimizer_attribute(feas, "Cuts", 0)
            set_optimizer_attribute(feas, "Threads", 1)

            @variable(dis, load[k = K()] >= 0, Int)
            @variable(dis, trip[k = K()] >= 0, Int)

            @constraint(dis, sum(load) >= value(o[i,j,t]))
            @constraint(dis, [k = K()], load[k] <= K(k).Q[i,j] * trip[k])
            @constraint(dis, [k = K()], trip[k] <= K(k).lim[i,j])

            @objective(dis, Min,
                sum(K(k).f[i,j] * load[k] + K(k).g[i,j] * trip[k] for k in K())
            ) #sum of all transport mode

            optimize!(dis)

            for k in K()
                if value(trip[k]) > 0
                    new_delivery = delivery(i,j,value(load[k]),value(trip[k]),k,t)
                    push!(deliveries,new_delivery)
                end
            end
        end
    end

    return deliveries
end #COMPUTE AGGREGATED SOLUTION

function directRounding(relaxed_sol::Vector{delivery})
    deliveries = Vector{delivery}()

    for r in relaxed_sol
        if r.trip > 0
            new_delivery = delivery(r.asal,r.tujuan,r.load,ceil(r.trip),r.k,r.t)
            push!(deliveries, new_delivery)
        end
    end

    return deliveries
end #ALTERNATIVE FOR FEASIBLE SOL

function optimalRounding(relaxed_sol::Vector{delivery})
    deliveries = Vector{delivery}()
    unique_asal = unique([p.asal for p in relaxed_sol])
    unique_tujuan = unique([p.tujuan for p in relaxed_sol])
    unique_t = unique([p.t for p in relaxed_sol])

    for i in unique_asal, j in unique_tujuan, t in unique_t
        active_delivery = filter(p -> p.asal == i && p.tujuan == j && p.t == t,relaxed_sol)
        aggregate = sum([r.load for r in active_delivery]) #collect aggregate delivery

        if aggregate > 0
            dis = Model(get_optimizer())
            set_silent(dis)

            set_optimizer_attribute(dis, "Cuts", 0)
            set_optimizer_attribute(dis, "Threads", 1)

            @variable(dis, o[k = K()] >= 0, Int)
            @variable(dis, p[k = K()] >= 0, Int)

            @constraint(dis, sum(o) == aggregate)
            @constraint(dis, [k = K()], o[k] <= K(k).Q[i,j] * p[k])
            @constraint(dis, [k = K()], p[k] <= K(k).lim[i,j])

            @objective(dis, Min, sum(K(k).f[i,j] * o[k] + K(k).g[i,j] * p[k] for k in K()))

            optimize!(dis)

            for k in K()
                if value(p[k]) > 0
                    new_delivery = delivery(i,j,value(o[k]),value(p[k]),k,t)
                    push!(deliveries,new_delivery)
                end
            end
        end
    end

    return deliveries
end #ALTERNATIVE FOR FEASIBLE SOL

function optimizeVolume(deliveries::Vector{delivery})
    OV = Model(get_optimizer())
    #set_silent(OV)
    set_optimizer_attribute(OV, "Threads", 1)
    set_optimizer_attribute(OV, "Cuts", 0)
    set_optimizer_attribute(OV, "TimeLimit", length(V()))

    @variable(OV, o[r = deliveries])
    @variable(OV, p[r = deliveries], Int) #bisa nambahin bisa ngurangin
    @variable(OV, I[i = V(), t = vcat(first(T()) - 1, T())], Int)

    @constraint(OV, [i = V(), t = T()],
        I[i,t-1] - I[i,t] - d(i,t) +
        sum(r.load + o[r] for r in filter(p -> p.tujuan == i && p.t == t,deliveries)) -
        sum(r.load + o[r] for r in filter(p -> p.asal == i && p.t == t,deliveries)) == 0
    ) #Inventory Balance

    @constraint(OV, [i = V(), t = T()],
        V(i).MIN <= I[i,t] <= V(i).MAX
    ) #Inventory Max Min

    @constraint(OV, [i = V()],
        I[i,first(T()) - 1] == V(i).START
    ) #Inventory Start

    @constraint(OV, [r = deliveries],
        (r.load + o[r]) <= K(r.k).Q[r.asal,r.tujuan] * (r.trip + p[r])
    ) #quantity to trips proportion from ij

    @constraint(OV, [r=deliveries],
        (r.load + o[r]) >= 0
    )

    @constraint(OV, [r = deliveries],
        (r.trip + p[r]) <= K(r.k).lim[r.asal,r.tujuan]
    )

    @constraint(OV, [r=deliveries],
        (r.trip + p[r]) >= 0
    )

    @objective(OV, Min,
        sum(
            (r.load + o[r]) * K(r.k).f[r.asal,r.tujuan] +
            (r.trip + p[r]) * K(r.k).g[r.asal,r.tujuan]
            for r in deliveries
        ) +
        sum(
            V(i).h * I[i,t]
            for i in V(), t in T()
        )
    )

    optimize!(OV)

    replacement = Vector{delivery}()
    for r in deliveries
        new_load = round(value(o[r]) + r.load)
        new_trip = round(value(p[r]) + r.trip) #revise trip

        if new_trip > 0
            new_delivery = delivery(r.asal, r.tujuan, new_load, new_trip, r.k, r.t)
            push!(replacement,new_delivery)
        end
    end #REVISE TRIP NUMBER

    return replacement
end #reoptimize the volume from trips

function iterOV(deliveries::Vector{delivery})
    mem = length(deliveries)
    deliveries = optimizeVolume(deliveries)

    while mem != length(deliveries)
        mem = length(deliveries)
        deliveries = optimizeVolume(deliveries)
    end

    return deliveries
end #iterate optimize volume until no improvement

function deliveryCost(deliveries)
    if isa(deliveries,Vector{delivery})
        total_cost = 0

        for r in deliveries
            indv_cost = r.load * K(r.k).f[r.asal,r.tujuan] + r.trip * K(r.k).g[r.asal,r.tujuan]
            total_cost += indv_cost
        end

        return total_cost
    else
        indv_cost = deliveries.load * K(deliveries.k).f[deliveries.asal,deliveries.tujuan] + deliveries.trip * K(deliveries.k).g[deliveries.asal,deliveries.tujuan]

        return indv_cost
    end
end #calculate cost of deliveries

function inventoryCost(deliveries::Vector{delivery})
    total_cost = 0
    I = stockLevels(deliveries)
    total_cost = sum(V(i).h * I[i,t] for i in V(), t in T())

    return total_cost
end #calculate inventory costs implied by deliveries

totalCost(deliveries) = deliveryCost(deliveries) + inventoryCost(deliveries)

function stockLevels(deliveries::Vector{delivery})
    I = JuMP.Containers.DenseAxisArray{Int64}(undef, V(), vcat(first(T()) - 1, T()))
    for i in V()
        I[i,first(T()) - 1] = V(i).START
    end #fill starting inventory level

    for i in V(), t in T()
        IN = filter(p -> p.tujuan == i && p.t == t,deliveries)
        OUT = filter(p -> p.asal == i && p.t == t,deliveries)

        I[i,t] = I[i,t-1] + sum([r.load for r in IN]) - sum([r.load for r in OUT]) - d(i,t)
    end

    return I
end #compute stock levels

function feasibility(deliveries::Vector{delivery})
    I = stockLevels(deliveries)
    infeasibilities = Vector{NamedTuple}()

    for i in V(), t in T()
        if I[i,t] > V(i).MAX || I[i,t] < V(i).MIN
            push!(infeasibilities,(i=i,t=t))
        end
    end

    if !isempty(infeasibilities)
        println("infeasible")
    end

    return infeasibilities
end #track infeasibilities in inventory level

function random_removal(deliveries::Vector{delivery})
    destroyed = deepcopy(deliveries) #make independent duplicate
    removed = Vector{delivery}()

    vertex_set = deepcopy(V()) #vertex set to remove
    percent_destroy = rand(truncated(Normal(0, 1), 0, 0.5))
    while length(removed)/(length(removed)+length(destroyed)) < percent_destroy
        key = rand(vertex_set) #random combination of random length
        deleteat!(vertex_set, findall(x -> x == key, vertex_set)) #remove from candidate

        to_remove = filter(p -> p.tujuan == key, destroyed)
        for r in to_remove
            deleteat!(destroyed, findall(x -> x == r, destroyed))
        end #remove all

        append!(removed,to_remove)
    end

    return destroyed,removed
end #removal based on cost

function worst_removal(deliveries::Vector{delivery})
    destroyed = deepcopy(deliveries) #make independent duplicate
    removed = Vector{delivery}()

    total_cost = deliveryCost(deliveries) #patokan percentage

    cost_group = DataFrame(c=String[],cost=Float64[],delivs=Vector{delivery}[])
    for c in V()
        deli_to_c = filter(p -> p.tujuan == c,destroyed)
        append!(cost_group,DataFrame(c=c,cost=deliveryCost(deli_to_c),delivs=[deli_to_c]))
    end
    sort!(cost_group,2;rev = true)

    percent_destroy = rand(truncated(Normal(0, 1), 0, 0.5))
    while deliveryCost(removed)/total_cost < percent_destroy
        left_cost = sum(cost_group.cost)
        p = [r.cost/left_cost for r in eachrow(cost_group)]
        roulette = rand()

        fitness_sum = 0
        chosen = 0
        for f in 1:length(p)
            fitness_sum += p[f]
            if fitness_sum > roulette
                chosen = f
                break
            end
        end

        to_remove = cost_group.delivs[chosen]
        append!(removed,to_remove)
        delete!(cost_group,chosen)
    end

    for r in removed
        deleteat!(destroyed, findall(x -> x == r, destroyed))
    end #remove all

    return destroyed,removed
end #removal based on random proportion

function repair(destroyed::Vector{delivery},removed::Vector{delivery})
    feas = Model(get_optimizer())
    set_silent(feas)

    set_optimizer_attribute(feas, "Threads", 1)
    set_optimizer_attribute(feas, "Cuts", 0)
    set_optimizer_attribute(feas, "TimeLimit", log(length(V())))

    @variable(feas, o[r = removed] >= 0, start = round(r.load))
    @variable(feas, p[r = removed] >= 0, Int)
    @variable(feas, I[i = V(), t = vcat(first(T()) - 1, T())], Int)

    @constraint(feas, balance[i = V(), t = T()],
        I[i,t-1] - I[i,t] - d(i,t) +
        sum(m.load for m in filter(p -> p.tujuan == i && p.t == t, destroyed)) -
        sum(m.load for m in filter(p -> p.asal == i && p.t == t, destroyed)) +
        sum(o[r] for r in filter(p -> p.tujuan == i && p.t == t, removed)) -
        sum(o[r] for r in filter(p -> p.asal == i && p.t == t, removed)) == 0
    ) #inventory balance

    @constraint(feas, [i = V(), t = T()],
        V(i).MIN <= I[i,t] <= V(i).MAX
    ) #Inventory Max Min

    @constraint(feas, [i = V()],
        I[i,first(T()) - 1] == V(i).START
    ) #Inventory Start

    @constraint(feas, [r = removed],
        o[r] <= K(r.k).Q[r.asal,r.tujuan] * p[r]
    ) #load to trips proportion from ij

    @constraint(feas, [r = removed],
        p[r] <= K(r.k).lim[r.asal,r.tujuan]
    ) #limit of trips from ij

    @objective(feas, Min,
        sum(
            K(r.k).g[r.asal,r.tujuan] * p[r] + K(r.k).f[r.asal,r.tujuan] * o[r]
            for r in removed
        ) +
        sum(
            V(i).h * I[i,t]
            for i in V(), t in T()
        )
    )

    optimize!(feas)

    new_deliveries = Vector{delivery}()
    for r in removed
        if value(p[r]) > 0
            new_delivery = delivery(
                r.asal, r.tujuan, round(value(o[r])), round(value(p[r])), r.k, r.t
            )
            push!(new_deliveries,new_delivery)
        end
    end #EXTRACT ADDITIONAL SOLUTION

    append!(destroyed,new_deliveries)

    return destroyed
end #make the solution feasible, if not optimal

function lns(x::Vector{delivery};method::String,stop_crit::String,stop_val::Int64)
    x_best = deepcopy(x)
    counter = 0

    limit = stop_val
    while counter < limit #&& totalCost(x_best) > 3.475e10
        if method == "worst"
            destroyed,removed = worst_removal(x)
            x_temp = repair(destroyed,removed)
        elseif method == "random"
            destroyed,removed = random_removal(x)
            x_temp = repair(destroyed,removed)
        end

        if isempty(feasibility(x_temp))
            if totalCost(x_temp) < totalCost(x)
                x = deepcopy(x_temp)

                if totalCost(x_temp) < totalCost(x_best)
                    x_best = deepcopy(x_temp)
                    println("new best cost: $(totalCost(x_best)) after $counter iter")

                    if stop_crit == "stagnation"
                        counter = 0 #resest counter
                    end
                else
                    if stop_crit == "stagnation"
                        counter += 1 #no improvement
                    end #println("counter: $counter")
                end
            else
                if stop_crit == "stagnation"
                    counter += 1 #no improvement
                end #println("counter: $counter")
            end
        end

        if stop_crit == "steps"
            counter += 1
        end

        println("counter: $counter")
    end

    return x_best
end

function stockouts(deliveries)
    I = stockLevels(deliveries)

    for i in V(), t in T()
        if I[i,t] >= V(i).MIN
            I[i,t] = 0
        else
            I[i,t] = V(i).MIN - I[i,t]
        end
    end

    return I
end

function overstocks(deliveries)
    I = stockLevels(deliveries)

    for i in V(), t in T()
        if I[i,t] <= V(i).MAX
            I[i,t] = 0
        else
            I[i,t] = I[i,t] - V(i).MAX
        end
    end

    return I
end

function available_supply(deliveries)
    I = stockLevels(deliveries)

    for i in V(), t in T()
        if I[i,t] > V(i).MIN
            I[i,t] = I[i,t] - V(i).MIN
        else
            I[i,t] = 0
        end
    end

    return I
end

function available_vaccancy(deliveries)
    I = stockLevels(deliveries)

    for i in V(), t in T()
        if I[i,t] < V(i).MAX
            I[i,t] = V(i).MAX - I[i,t]
        else
            I[i,t] = 0
        end
    end

    return I
end


text:
function check_chronology!(sim::Simulation, key::Pair, sync::Synchronize)
    source_problem = get_problems(sim)[key.first]
    destination_problem = get_problems(sim)[key.second]
    source_problem_horizon = get_horizon(source_problem)
    destination_problem_horizon = get_horizon(destination_problem)
    sequence = get_sequence(sim)
    source_problem_interval = get_interval(sequence, key.first)
    destination_problem_interval = get_interval(sequence, key.second)

    source_problem_resolution = get_resolution(source_problem)
    @debug source_problem_resolution, destination_problem_interval
    # How many times the second problem executes per solution retireved from the source_problem.
    # E.g. source_problem_resolution = 1 Hour, destination_problem_interval = 5 minutes => 12 executions per solution
    destination_problem_executions_per_solution =
        Int(source_problem_resolution / destination_problem_interval)
    # Number of periods in the horizon that will be synchronized between the source_problem and the destination_problem
    source_problem_sync = sync.periods

    if source_problem_sync > source_problem_horizon
        throw(
            IS.ConflictingInputsError(
                "The lookahead length $(source_problem_horizon) in problem is insufficient to syncronize with $(source_problem_sync) feedforward periods",
            ),
        )
    end

    if (source_problem_sync % destination_problem_executions_per_solution) != 0
        throw(
            IS.ConflictingInputsError(
                "The current configuration implies $(source_problem_sync / destination_problem_executions_per_solution) executions of $(key.second) per execution of $(key.first). The number of Synchronize periods $(sync.periods) in problem $(key.first) needs to be a mutiple of the number of problem $(key.second) execution for every problem $(key.first) interval.",
            ),
        )
    end

    return
end

function check_chronology!(sim::Simulation, key::Pair, ::Consecutive)
    source_problem = get_problems(sim)[key.first]
    destination_problem = get_problems(sim)[key.second]
    source_problem_horizon = get_horizon(source_problem)
    destination_problem_horizon = get_horizon(destination_problem)
    if source_problem_horizon != source_problem_interval
        @warn(
            "Consecutive Chronology Requires the same interval and horizon, the parameter horizon = $(source_problem_horizon) in problem $(key.first) will be replaced with $(source_problem_interval). If this is not the desired behviour consider changing your chronology to RecedingHorizon"
        )
    end
    get_sequence(sim).horizons[key.first] = get_interval(sim, key.first)
    return
end

check_chronology!(sim::Simulation, key::Pair, ::RecedingHorizon) = nothing
check_chronology!(sim::Simulation, key::Pair, ::FullHorizon) = nothing
# TODO: Add missing check
check_chronology!(sim::Simulation, key::Pair, ::Range) = nothing

function check_chronology!(
    sim::Simulation,
    key::Pair,
    ::T,
) where {T <: FeedForwardChronology}
    error("Chronology $(T) not implemented")
    return
end

############################ FeedForward Definitions ########################################

struct UpperBoundFF <: AbstractAffectFeedForward
    variable_source_problem::Symbol
    affected_variables::Vector{Symbol}
    cache::Union{Nothing, Type{<:AbstractCache}}
    function UpperBoundFF(
        variable_source_problem::AbstractString,
        affected_variables::Vector{<:AbstractString},
        cache::Union{Nothing, Type{<:AbstractCache}},
    )
        new(Symbol(variable_source_problem), Symbol.(affected_variables), cache)
    end
end

function UpperBoundFF(; variable_source_problem, affected_variables)
    return UpperBoundFF(variable_source_problem, affected_variables, nothing)
end

get_variable_source_problem(p::UpperBoundFF) = p.variable_source_problem

struct RangeFF <: AbstractAffectFeedForward
    variable_source_problem_ub::Symbol
    variable_source_problem_lb::Symbol
    affected_variables::Vector{Symbol}
    cache::Union{Nothing, Type{<:AbstractCache}}
    function RangeFF(
        variable_source_problem_ub::AbstractString,
        variable_source_problem_lb::AbstractString,
        affected_variables::Vector{<:AbstractString},
        cache::Union{Nothing, Type{<:AbstractCache}},
    )
        new(
            Symbol(variable_source_problem_ub),
            Symbol(variable_source_problem_lb),
            Symbol.(affected_variables),
            cache,
        )
    end
end

function RangeFF(;
    variable_source_problem_ub,
    variable_source_problem_lb,
    affected_variables,
)
    return RangeFF(
        variable_source_problem_ub,
        variable_source_problem_lb,
        affected_variables,
        nothing,
    )
end

get_bounds_source_problem(p::RangeFF) =
    (p.variable_source_problem_lb, p.variable_source_problem_lb)

struct SemiContinuousFF <: AbstractAffectFeedForward
    binary_source_problem::Symbol
    affected_variables::Vector{Symbol}
    cache::Union{Nothing, Type{<:AbstractCache}}
    function SemiContinuousFF(
        binary_source_problem::AbstractString,
        affected_variables::Vector{<:AbstractString},
        cache::Union{Nothing, Type{<:AbstractCache}},
    )
        new(Symbol(binary_source_problem), Symbol.(affected_variables), cache)
    end
end

function SemiContinuousFF(; binary_source_problem, affected_variables)
    return SemiContinuousFF(binary_source_problem, affected_variables, nothing)
end

get_binary_source_problem(p::SemiContinuousFF) = p.binary_source_problem
get_affected_variables(p::AbstractAffectFeedForward) = p.affected_variables

struct IntegralLimitFF <: AbstractAffectFeedForward
    variable_source_problem::Symbol
    affected_variables::Vector{Symbol}
    cache::Union{Nothing, Type{<:AbstractCache}}
    function IntegralLimitFF(
        variable_source_problem::AbstractString,
        affected_variables::Vector{<:AbstractString},
        cache::Union{Nothing, Type{<:AbstractCache}},
    )
        new(Symbol(variable_source_problem), Symbol.(affected_variables), cache)
    end
end

function IntegralLimitFF(; variable_source_problem, affected_variables)
    return IntegralLimitFF(variable_source_problem, affected_variables, nothing)
end

get_variable_source_problem(p::IntegralLimitFF) = p.variable_source_problem

struct ParameterFF <: AbstractAffectFeedForward
    variable_source_problem::Symbol
    affected_parameters::Any
    function ParameterFF(variable_source_problem::AbstractString, affected_parameters)
        new(Symbol(variable_source_problem), affected_parameters)
    end
end

function ParameterFF(; variable_source_problem, affected_parameters)
    return ParameterFF(variable_source_problem, affected_parameters)
end

####################### Feed Forward Affects ###############################################

@doc raw"""
        ub_ff(optimization_container::OptimizationContainer,
              cons_name::Symbol,
              constraint_infos::Vector{DeviceRangeConstraintInfo},
              param_reference::UpdateRef,
              var_name::Symbol)

Constructs a parametrized upper bound constraint to implement feedforward from other models.
The Parameters are initialized using the uppper boundary values of the provided variables.

# Constraints
``` variable[var_name, t] <= param_reference[var_name] ```

# LaTeX

`` x \leq param^{max}``

# Arguments
* optimization_container::OptimizationContainer : the optimization_container model built in PowerSimulations
* cons_name::Symbol : name of the constraint
* param_reference : Reference to the PJ.ParameterRef used to determine the upperbound
* var_name::Symbol : the name of the continuous variable
"""
function ub_ff(
    optimization_container::OptimizationContainer,
    cons_name::Symbol,
    constraint_infos::Vector{DeviceRangeConstraintInfo},
    param_reference::UpdateRef,
    var_name::Symbol,
)
    time_steps = model_time_steps(optimization_container)
    ub_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "ub")
    variable = get_variable(optimization_container, var_name)

    axes = JuMP.axes(variable)
    set_name = axes[1]
    @assert axes[2] == time_steps
    container = add_param_container!(optimization_container, param_reference, set_name)
    param_ub = get_parameter_array(container)
    multiplier_ub = get_multiplier_array(container)
    con_ub = add_cons_container!(optimization_container, ub_name, set_name, time_steps)

    for constraint_info in constraint_infos
        name = get_component_name(constraint_info)
        value = JuMP.upper_bound(variable[name, 1])
        param_ub[name] = add_parameter(optimization_container.JuMPmodel, value)
        # default set to 1.0, as this implementation doesn't use multiplier
        multiplier_ub[name] = 1.0
        for t in time_steps
            expression_ub = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
            for val in constraint_info.additional_terms_ub
                JuMP.add_to_expression!(expression_ub, variable[name, t])
            end
            con_ub[name, t] = JuMP.@constraint(
                optimization_container.JuMPmodel,
                expression_ub <= param_ub[name] * multiplier_ub[name]
            )
        end
    end
    return
end

@doc raw"""
        range_ff(optimization_container::OptimizationContainer,
                        cons_name::Symbol,
                        param_reference::NTuple{2, UpdateRef},
                        var_name::Symbol)

Constructs min/max range parametrized constraint from device variable to include feedforward.

# Constraints

``` param_reference[1][var_name] <= variable[var_name, t] ```
``` variable[var_name, t] <= param_reference[2][var_name] ```

where r in range_data.

# LaTeX

`` param^{min} \leq x ``
`` x \leq param^{max}``

# Arguments
* optimization_container::OptimizationContainer : the optimization_container model built in PowerSimulations
* param_reference::NTuple{2, UpdateRef} : Tuple with the lower bound and upper bound parameter reference
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the continuous variable
"""
function range_ff(
    optimization_container::OptimizationContainer,
    cons_name::Symbol,
    constraint_infos::Vector{DeviceRangeConstraintInfo},
    param_reference::NTuple{2, UpdateRef},
    var_name::Symbol,
)
    time_steps = model_time_steps(optimization_container)
    ub_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "ub")
    lb_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "lb")

    variable = get_variable(optimization_container, var_name)
    # Used to make sure the names are consistent between the variable and the infos
    axes = JuMP.axes(variable)
    set_name = axes[1]
    @assert axes[2] == time_steps

    # Create containers for the constraints
    container_lb =
        add_param_container!(optimization_container, param_reference[1], set_name)
    param_lb = get_parameter_array(container_lb)
    multiplier_lb = get_multiplier_array(container_lb)
    container_ub =
        add_param_container!(optimization_container, param_reference[2], set_name)
    param_ub = get_parameter_array(container_ub)
    multiplier_ub = get_multiplier_array(container_ub)
    # Create containers for the parameters
    con_lb = add_cons_container!(optimization_container, lb_name, set_name, time_steps)
    con_ub = add_cons_container!(optimization_container, ub_name, set_name, time_steps)

    for constraint_info in constraint_infos
        name = get_component_name(constraint_info)
        param_lb[name] = add_parameter(
            optimization_container.JuMPmodel,
            JuMP.lower_bound(variable[name, 1]),
        )
        param_ub[name] = add_parameter(
            optimization_container.JuMPmodel,
            JuMP.upper_bound(variable[name, 1]),
        )
        # default set to 1.0, as this implementation doesn't use multiplier
        multiplier_ub[name] = 1.0
        multiplier_lb[name] = 1.0
        for t in time_steps
            expression_ub = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
            for val in constraint_info.additional_terms_ub
                JuMP.add_to_expression!(expression_ub, variable[name, t])
            end
            expression_lb = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
            for val in constraint_info.additional_terms_lb
                JuMP.add_to_expression!(expression_lb, variable[name, t], -1.0)
            end
            con_ub[name, t] = JuMP.@constraint(
                optimization_container.JuMPmodel,
                expression_ub <= param_ub[name] * multiplier_ub[name]
            )
            con_lb[name, t] = JuMP.@constraint(
                optimization_container.JuMPmodel,
                expression_lb >= param_lb[name] * multiplier_lb[name]
            )
        end
    end

    return
end

@doc raw"""
            semicontinuousrange_ff(optimization_container::OptimizationContainer,
                                    cons_name::Symbol,
                                    var_name::Symbol,
                                    param_reference::UpdateRef)

Constructs min/max range constraint from device variable with parameter setting.

# Constraints
If device min = 0:

``` variable[var_name, t] <= r[2].max*param_reference[var_name] ```

Otherwise:

``` variable[var_name, t] <= r[2].max*param_reference[var_name] ```

``` variable[var_name, t] >= r[2].min*param_reference[var_name] ```

where r in range_data.

# LaTeX

`` 0.0 \leq x^{var} \leq r^{max} x^{param}, \text{ for } r^{min} = 0 ``

`` r^{min} x^{param} \leq x^{var} \leq r^{min} x^{param}, \text{ otherwise } ``

# Arguments
* optimization_container::OptimizationContainer : the optimization_container model built in PowerSimulations
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the continuous variable
* param_reference::UpdateRef : UpdateRef of the parameter
"""
function semicontinuousrange_ff(
    optimization_container::OptimizationContainer,
    cons_name::Symbol,
    constraint_infos::Vector{DeviceRangeConstraintInfo},
    param_reference::UpdateRef,
    var_name::Symbol,
)
    time_steps = model_time_steps(optimization_container)
    ub_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "ub")
    lb_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "lb")
    variable = get_variable(optimization_container, var_name)
    # Used to make sure the names are consistent between the variable and the infos
    axes = JuMP.axes(variable)
    set_name = [get_component_name(ci) for ci in constraint_infos]
    @assert axes[2] == time_steps
    container = add_param_container!(optimization_container, param_reference, set_name)
    multiplier = get_multiplier_array(container)
    param = get_parameter_array(container)
    con_ub = add_cons_container!(optimization_container, ub_name, set_name, time_steps)
    con_lb = add_cons_container!(optimization_container, lb_name, set_name, time_steps)

    for constraint_info in constraint_infos
        name = get_component_name(constraint_info)
        ub_value = JuMP.upper_bound(variable[name, 1])
        lb_value = JuMP.lower_bound(variable[name, 1])
        @debug "SemiContinuousFF" name ub_value lb_value
        # default set to 1.0, as this implementation doesn't use multiplier
        multiplier[name] = 1.0
        param[name] = add_parameter(optimization_container.JuMPmodel, 1.0)
        for t in time_steps
            expression_ub = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
            for val in constraint_info.additional_terms_ub
                JuMP.add_to_expression!(
                    expression_ub,
                    get_variable(optimization_container, val)[name, t],
                )
            end
            expression_lb = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
            for val in constraint_info.additional_terms_lb
                JuMP.add_to_expression!(
                    expression_lb,
                    get_variable(optimization_container, val)[name, t],
                    -1.0,
                )
            end
            mul_ub = ub_value * multiplier[name]
            mul_lb = lb_value * multiplier[name]
            con_ub[name, t] = JuMP.@constraint(
                optimization_container.JuMPmodel,
                expression_ub <= mul_ub * param[name]
            )
            con_lb[name, t] = JuMP.@constraint(
                optimization_container.JuMPmodel,
                expression_lb >= mul_lb * param[name]
            )
        end
    end

    # If the variable was a lower bound != 0, not removing the LB can cause infeasibilities
    for v in variable
        if JuMP.has_lower_bound(v)
            @debug "lb reset" v
            JuMP.set_lower_bound(v, 0.0)
        end
    end

    return
end

@doc raw"""
        integral_limit_ff(optimization_container::OptimizationContainer,
                        cons_name::Symbol,
                        param_reference::UpdateRef,
                        var_name::Symbol)

Constructs a parametrized integral limit constraint to implement feedforward from other models.
The Parameters are initialized using the upper boundary values of the provided variables.

# Constraints
``` sum(variable[var_name, t] for t in time_steps)/length(time_steps) <= param_reference[var_name] ```

# LaTeX

`` \sum_{t} x \leq param^{max}``
`` \sum_{t} x * DeltaT_lower \leq param^{max} * DeltaT_upper ``
    `` P_LL - P_max * ON_upper <= 0.0 ``
    `` P_LL - P_min * ON_upper >= 0.0 ``

# Arguments
* optimization_container::OptimizationContainer : the optimization_container model built in PowerSimulations
* cons_name::Symbol : name of the constraint
* param_reference : Reference to the PJ.ParameterRef used to determine the upperbound
* var_name::Symbol : the name of the continuous variable
"""
function integral_limit_ff(
    optimization_container::OptimizationContainer,
    cons_name::Symbol,
    param_reference::UpdateRef,
    var_name::Symbol,
)
    time_steps = model_time_steps(optimization_container)
    ub_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "integral_limit")
    variable = get_variable(optimization_container, var_name)

    axes = JuMP.axes(variable)
    set_name = axes[1]

    @assert axes[2] == time_steps
    container_ub = add_param_container!(optimization_container, param_reference, set_name)
    param_ub = get_parameter_array(container_ub)
    multiplier_ub = get_multiplier_array(container_ub)
    con_ub = add_cons_container!(optimization_container, ub_name, set_name)

    for name in axes[1]
        value = JuMP.upper_bound(variable[name, 1])
        param_ub[name] = add_parameter(optimization_container.JuMPmodel, value)
        # default set to 1.0, as this implementation doesn't use multiplier
        multiplier_ub[name] = 1.0
        con_ub[name] = JuMP.@constraint(
            optimization_container.JuMPmodel,
            sum(variable[name, t] for t in time_steps) / length(time_steps) <=
            param_ub[name] * multiplier_ub[name]
        )
    end
end

########################## FeedForward Constraints #########################################
function feedforward!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, <:AbstractDeviceFormulation},
    ff_model::Nothing,
) where {T <: PSY.Component}
    return
end

function feedforward!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, <:AbstractDeviceFormulation},
    ff_model::UpperBoundFF,
) where {T <: PSY.StaticInjection}
    constraint_infos = Vector{DeviceRangeConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        limits = PSY.get_active_power_limits(d)
        constraint_info = DeviceRangeConstraintInfo(name, limits)
        add_device_services!(constraint_info, d, model)
        constraint_infos[ix] = constraint_info
    end
    for prefix in get_affected_variables(ff_model)
        var_name = make_variable_name(prefix, T)
        parameter_ref = UpdateRef{JuMP.VariableRef}(var_name)
        ub_ff(
            optimization_container,
            constraint_name(FEEDFORWARD_UB, T),
            constraint_infos,
            parameter_ref,
            var_name,
        )
    end
end

function feedforward!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, <:AbstractDeviceFormulation},
    ff_model::SemiContinuousFF,
) where {T <: PSY.StaticInjection}
    bin_var = make_variable_name(get_binary_source_problem(ff_model), T)
    parameter_ref = UpdateRef{JuMP.VariableRef}(bin_var)
    constraint_infos = Vector{DeviceRangeConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        limits = PSY.get_active_power_limits(d)
        constraint_info = DeviceRangeConstraintInfo(name, limits)
        add_device_services!(constraint_info, d, model)
        constraint_infos[ix] = constraint_info
    end
    for prefix in get_affected_variables(ff_model)
        var_name = make_variable_name(prefix, T)
        semicontinuousrange_ff(
            optimization_container,
            make_constraint_name(FEEDFORWARD_BIN, T),
            constraint_infos,
            parameter_ref,
            var_name,
        )
    end
end

function feedforward!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, <:AbstractDeviceFormulation},
    ff_model::IntegralLimitFF,
) where {T <: PSY.StaticInjection}
    for prefix in get_affected_variables(ff_model)
        var_name = make_variable_name(prefix, T)
        parameter_ref = UpdateRef{JuMP.VariableRef}(var_name)
        integral_limit_ff(
            optimization_container,
            make_constraint_name(FEEDFORWARD_INTEGRAL_LIMIT, T),
            parameter_ref,
            var_name,
        )
    end
end

######################### FeedForward Variables Updating #####################################
# This makes the choice in which variable to get from the results.
function get_problem_variable(
    chron::RecedingHorizon,
    problems::Pair{OperationsProblem{T}, OperationsProblem{U}},
    device_name::AbstractString,
    var_ref::UpdateRef,
) where {T, U <: AbstractOperationsProblem}
    variable =
        get_variable(problems.first.internal.optimization_container, var_ref.access_ref)
    step = axes(variable)[2][chron.periods]
    var = variable[device_name, step]
    if JuMP.is_binary(var)
        return round(JuMP.value(var))
    else
        return JuMP.value(var)
    end
end

function get_problem_variable(
    ::Consecutive,
    problems::Pair{OperationsProblem{T}, OperationsProblem{U}},
    device_name::String,
    var_ref::UpdateRef,
) where {T, U <: AbstractOperationsProblem}
    variable =
        get_variable(problems.first.internal.optimization_container, var_ref.access_ref)
    step = axes(variable)[2][get_end_of_interval_step(problems.first)]
    var = variable[device_name, step]
    if JuMP.is_binary(var)
        return round(JuMP.value(var))
    else
        return JuMP.value(var)
    end
end

function get_problem_variable(
    chron::Synchronize,
    problems::Pair{OperationsProblem{T}, OperationsProblem{U}},
    device_name::String,
    var_ref::UpdateRef,
) where {T, U <: AbstractOperationsProblem}
    variable =
        get_variable(problems.first.internal.optimization_container, var_ref.access_ref)
    e_count = get_execution_count(problems.second)
    wait_count = get_execution_wait_count(get_trigger(chron))
    index = (floor(e_count / wait_count) + 1)
    step = axes(variable)[2][Int(index)]
    var = variable[device_name, step]
    if JuMP.is_binary(var)
        return round(JuMP.value(var))
    else
        return JuMP.value(var)
    end
end

function get_problem_variable(
    ::FullHorizon,
    problems::Pair{OperationsProblem{T}, OperationsProblem{U}},
    device_name::String,
    var_ref::UpdateRef,
) where {T, U <: AbstractOperationsProblem}
    variable =
        get_variable(problems.first.internal.optimization_container, var_ref.access_ref)
    vars = variable[device_name, :]
    if JuMP.is_binary(first(vars))
        return round.(JuMP.value(vars))
    else
        return JuMP.value.(vars)
    end
end

function get_problem_variable(
    chron::Range,
    problems::Pair{OperationsProblem{T}, OperationsProblem{U}},
    device_name::String,
    var_ref::UpdateRef,
) where {T, U <: AbstractOperationsProblem}
    variable =
        get_variable(problems.first.internal.optimization_container, var_ref.access_ref)
    vars = variable[device_name, chron.range]
    if JuMP.is_binary(first(vars))
        return round.(JuMP.value(vars))
    else
        return JuMP.value.(vars)
    end
end

function feedforward_update!(
    destination_problem::OperationsProblem,
    source_problem::OperationsProblem,
    chronology::FeedForwardChronology,
    param_reference::UpdateRef{JuMP.VariableRef},
    param_array::JuMPParamArray,
    current_time::Dates.DateTime,
)
    trigger = get_trigger(chronology)
    if trigger_update(trigger)
        for device_name in axes(param_array)[1]
            var_value = get_problem_variable(
                chronology,
                (source_problem => destination_problem),
                device_name,
                param_reference,
            )
            previous_value = PJ.value(param_array[device_name])
            PJ.set_value(param_array[device_name], var_value)
            IS.@record :simulation FeedForwardUpdateEvent(
                "FeedForward",
                current_time,
                param_reference,
                device_name,
                var_value,
                previous_value,
                destination_problem,
                source_problem,
            )
        end
        reset_trigger_count!(trigger)
    end
    update_count!(trigger)
    return
end

function attach_feedforward(model, ff::AbstractAffectFeedForward)
    model.feedforward = ff
    return
end


text:
export read_meta, read_meta!

@doc """
    read_meta

Generic wrapper for reading metadata (i.e., instrument parameters, responses).

    S = read_meta(fmt, filestr [, keywords])

Read metadata in file format `fmt` matching file pattern `filestr` into a new
SeisData object `S`.

    read_meta!(S, fmt, filestr [, keywords])

Read metadata in file format `fmt` matching file pattern `filestr` into existing
SeisData object `S`.

### Supported File Formats
| Format                    | String          |
| :---                      | :---            |
| Dataless SEED             | dataless        |
| FDSN Station XML          | sxml            |
| SACPZ                     | sacpz           |
| SEED RESP                 | resp            |

### Keywords
|KW      | Used By  | Type      | Default   | Meaning                         |
|:---    |:---      |:---       |:---       |:---                             |
| msr    | sxml     | Bool      | false     | read full MultiStageResp?       |
| s      | all      | TimeSpec  |           | Start time                      |
| t      | all      | TimeSpec  |           | Termination (end) time          |
| units  | resp     | Bool      | false     | fill in MultiStageResp units?   |
|        | dataless |           |           |                                 |
| v      | all      | Int64     | 0         | verbosity                       |

### Notes
1. Unlike `read_data`, `read_meta` can't use `guess` for files of unknown type.
The reason is that most metadata formats are ASCII-based; generally only XML
files have reliable tests for uniqueness.

See also: SeisIO.KW, get_data, read_data
""" read_meta!
function read_meta!(S::GphysData, fmt::String, filestr::String;
  msr     ::Bool      = false                     ,  # read as MultiStageResp?
  s       ::TimeSpec  = "0001-01-01T00:00:00"     ,  # Start
  t       ::TimeSpec  = "9999-12-31T12:59:59"     ,  # End or Length (s)
  units   ::Bool      = false                     ,  # fill in units of CoeffResp stages?
  v       ::Int64     = KW.v                      ,  # verbosity level
  )

  one_file = safe_isfile(filestr)

  if fmt == "dataless"
    if one_file
      append!(S, read_dataless(filestr, s=s, t=t, v=v, units=units))
    else
      files = ls(filestr)
      for fname in files
        append!(S, read_dataless(fname, s=s, t=t, v=v, units=units))
      end
    end

  elseif fmt == "resp"
    read_seed_resp!(S, filestr, units=units)

  elseif fmt == "sacpz"
    if one_file
      read_sacpz!(S, filestr)
    else
      files = ls(filestr)
      for fname in files
        read_sacpz!(S, fname)
      end
    end

  elseif fmt == "sxml"
    if one_file
      append!(S, read_sxml(filestr, s=s, t=t, v=v, msr=msr))
    else
      files = ls(filestr)
      for fname in files
        append!(S, read_sxml(fname, s=s, t=t, v=v, msr=msr))
      end
    end

  else
    error("Unknown file format!")
  end

  return nothing
end

@doc (@doc read_meta!)
function read_meta(fmt::String, filestr::String;
  msr     ::Bool      = false                     ,  # read as MultiStageResp?
  s       ::TimeSpec  = "0001-01-01T00:00:00"     ,  # Start
  t       ::TimeSpec  = "9999-12-31T12:59:59"     ,  # End or Length (s)
  units   ::Bool      = false                     ,  # fill in units of CoeffResp stages?
  v       ::Int64     = KW.v                      ,  # verbosity level
  )

  S = SeisData()
  read_meta!(S, fmt, filestr, msr=msr, s=s, t=t, units=units, v=v)
  return S
end


text:
struct FunctionField{F} <: Field
  f::F
end

function_field(f::Function) =  FunctionField(f)

function field_cache(f::FunctionField,x)
  nx = length(x)
  Te = eltype(x)
  c = zeros(return_type(f.f,Te),nx)
  CachedArray(c)
end

function evaluate_field!(c,f::FunctionField,x)
  nx = length(x)
  setsize!(c,(nx,))
  for i in eachindex(x)
    c[i] = f.f(x[i])
  end
  c
end

function field_gradient(f::FunctionField)
  FunctionField(gradient(f.f))
end


text:
using CooperationPolarization
using Test

@testset "CooperationPolarization.jl" begin

    # setup.jl: random_sets
    @test random_sets(2, 1, 1, 2).affiliations == [1, 2]
    @test random_sets(8, 3, 2, 2).affiliations == [1, 1, 1, 1, 2, 2, 2, 2]
    @test random_sets(8, 3, 2).affiliations == random_sets(8, 3, 2, 2).affiliations

    # setup.jl: vectobase3
    @test vectobase3([0; 0; 0]) == 0
    @test vectobase3([0; 1; 1]) == 4
    @test vectobase3([0;-1;-1]) == -4

end


text:
using StochasticDiffEq, Test, Random
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_2Dlinear
Random.seed!(100)
prob = prob_sde_2Dlinear

## Solve and plot
sol1 =solve(prob,RKMil(),abstol=1,reltol=0,adaptive=true)
err1 = sol1.errors[:final]

sol1 =solve(prob,RKMil(),abstol=1e-3,reltol=1e-3,adaptive=true)
err12 = sol1.errors[:final]
@test err12 < err1

sol1 =solve(prob,RKMilGeneral(),abstol=1,reltol=0,adaptive=true)
err11 = sol1.errors[:final]

sol1 =solve(prob,RKMilGeneral(),adaptive=true)
err112 = sol1.errors[:final]
@test err112 < err11

sol2 =solve(prob,SRI(),dt=1/2^(4),abstol=1,reltol=0)
err2 = sol2.errors[:final]

println("1e-1")
sol3 =solve(prob,SRI(),dt=1/2^(4),abstol=1e-1,reltol=0)
err3 = sol3.errors[:final]

println("1e-2")
sol4 =solve(prob,SRI(),dt=1/2^(4),abstol=1e-2,reltol=0)
err4 = sol4.errors[:final]
@test err4 < err2
println("1e-3")
sol5 =solve(prob,SRI(),dt=1/2^(4),abstol=1e-3,reltol=0)
err5 = sol5.errors[:final]
@test err5 < err3
println("""
Final error for the solutions were:
          $err1
          $err11
          $err2
          $err3
          $err4
          $err5
          """)

sol6 =solve(prob,SRI())

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq


text:
abstract type AbstractLPForm{T} <: MOI.ModelLike end

# FIXME Taken from Polyhedra.jl, we should maybe move this to MOIU
structural_nonzero_indices(a::SparseArrays.SparseVector) = SparseArrays.nonzeroinds(a)
structural_nonzero_indices(a::AbstractVector) = eachindex(a)

function _dot_terms(a::AbstractVector{T}) where T
    return MOI.ScalarAffineTerm{T}[
        MOI.ScalarAffineTerm{T}(a[i], MOI.VariableIndex(i))
        for i in structural_nonzero_indices(a) if !iszero(a[i])
    ]
end
function _dot(a::AbstractVector{T}) where T
    return MOI.ScalarAffineFunction(_dot_terms(a), zero(T))
end

function MOI.get(model::AbstractLPForm, ::MOI.ListOfModelAttributesSet)
    list = MOI.AbstractModelAttribute[]
    push!(list, MOI.ObjectiveSense())
    if model.sense != MOI.FEASIBILITY_SENSE
        push!(list, MOI.ObjectiveFunction{MOI.get(model, MOI.ObjectiveFunctionType())}())
    end
    return list
end
MOI.get(model::AbstractLPForm, ::MOI.ObjectiveSense) = model.sense
function MOI.get(model::AbstractLPForm{T},
                 ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}) where T
    return _dot(model.c)
end
function MOI.get(model::AbstractLPForm{T}, ::MOI.ObjectiveFunctionType) where T
    return MOI.ScalarAffineFunction{T}
end

MOI.get(::AbstractLPForm, ::MOI.ListOfVariableAttributesSet) = MOI.AbstractVariableAttribute[]
MOI.get(model::AbstractLPForm, ::MOI.NumberOfVariables) = length(model.c)
function MOI.get(model::AbstractLPForm, ::MOI.ListOfVariableIndices)
    # TODO return `collect` with MOI v0.9.15 (see https://github.com/jump-dev/MathOptInterface.jl/pull/1110)
    return collect(MOIU.LazyMap{MOI.VariableIndex}(
        i -> MOI.VariableIndex(i), # FIXME `LazyMap` needs a `Function` so cannot just give `MOI.VariableIndex`
        1:MOI.get(model, MOI.NumberOfVariables())
    ))
end

MOI.get(::AbstractLPForm, ::MOI.ListOfConstraintAttributesSet) = MOI.AbstractConstraintAttribute[]
function MOI.get(model::AbstractLPForm{T}, ::MOI.ConstraintFunction,
                 ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}}) where T
    return _dot(model.A[ci.value, :])
end

"""
    LPStandardForm{T, AT<:AbstractMatrix{T}, VT <: AbstractVector{T}} <: AbstractLPForm{T}

Represents a problem of the form:
```
sense ⟨c, x⟩
s.t.  A x == b
      x ≥ 0
```
"""
struct LPStandardForm{T, AT<:AbstractMatrix{T}, VT <: AbstractVector{T}} <: AbstractLPForm{T}
    sense::MOI.OptimizationSense
    c::VT
    A::AT
    b::VT
end

function MOI.get(::LPStandardForm{T}, ::MOI.ListOfConstraintTypesPresent) where T
    return [(MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
            (MOI.VectorOfVariables, MOI.Nonnegatives)]
end
const EQ{T} = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}
function MOI.get(
    model::LPStandardForm{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}
) where T
    # TODO return `collect` with MOI v0.9.15 (see https://github.com/jump-dev/MathOptInterface.jl/pull/1110)
    return collect(MOIU.LazyMap{EQ{T}}(
        i -> EQ{T}(i), # FIXME `LazyMap` needs a `Function` so cannot just give `EQ{T}`
        1:size(model.A, 1)
    ))
end
function MOI.get(model::LPStandardForm, ::MOI.ConstraintSet,
                 ci::EQ)
    return MOI.EqualTo(model.b[ci.value])
end
const NONNEG = MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}
function MOI.get(
    ::LPStandardForm,
    ::MOI.ListOfConstraintIndices{
        MOI.VectorOfVariables, MOI.Nonnegatives}
)
    return [NONNEG(1)]
end
function MOI.get(model::LPStandardForm, ::MOI.ConstraintFunction,
                 ci::NONNEG)
    return MOI.VectorOfVariables(MOI.get(model, MOI.ListOfVariableIndices()))
end
function MOI.get(model::LPStandardForm, ::MOI.ConstraintSet,
                 ci::NONNEG)
    return MOI.Nonnegatives(MOI.get(model, MOI.NumberOfVariables()))
end

"""
    LPGeometricForm{T, AT<:AbstractMatrix{T}, VT <: AbstractVector{T}} <: AbstractLPForm{T}

Represents a linear problem of the form:
```
sense ⟨c, x⟩
s.t.  Ax <= b
```
"""
struct LPGeometricForm{T, AT<:AbstractMatrix{T}, VT <: AbstractVector{T}} <: AbstractLPForm{T}
    sense::MOI.OptimizationSense
    c::VT
    A::AT
    b::VT
end

function MOI.get(::LPGeometricForm{T}, ::MOI.ListOfConstraintTypesPresent) where T
    return [(MOI.ScalarAffineFunction{T}, MOI.LessThan{T})]
end
const LT{T} = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}
function MOI.get(
    model::LPGeometricForm{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}
) where T
    # FIXME `copy_constraint` needs a `Vector`
    return collect(MOIU.LazyMap{LT{T}}(
        i -> LT{T}(i), # FIXME `LazyMap` needs a `Function` so cannot just give `LT{T}`
        1:size(model.A, 1))
    )
end
function MOI.get(model::LPGeometricForm, ::MOI.ConstraintSet,
                 ci::LT)
    return MOI.LessThan(model.b[ci.value])
end

abstract type LPMixedForm{T} <: AbstractLPForm{T} end

function MOI.get(model::LPMixedForm{T}, ::MOI.ListOfConstraintTypesPresent) where T
    list = Tuple{DataType, DataType}[]
    for S in [MOI.EqualTo{T}, MOI.Interval{T}, MOI.GreaterThan{T}, MOI.LessThan{T}]
        for F in [MOI.VariableIndex, MOI.ScalarAffineFunction{T}]
            if !iszero(MOI.get(model, MOI.NumberOfConstraints{F, S}()))
                push!(list, (F, S))
            end
        end
    end
    return list
end

struct BoundSense <: MOI.AbstractVariableAttribute end
function MOI.get(model::LPMixedForm, ::BoundSense, vi::MOI.VariableIndex)
    return _bound_sense(model.v_lb[vi.value], model.v_ub[vi.value])
end

const LinearBounds{T} = Union{MOI.EqualTo{T}, MOI.Interval{T}, MOI.GreaterThan{T}, MOI.LessThan{T}}

function MOI.get(model::LPMixedForm{T}, ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, S}) where {T, S <: LinearBounds{T}}
    s = _sense(S)
    return count(1:size(model.A, 1)) do i
        _constraint_bound_sense(model, i) == s
    end
end
const AFF{T, S} = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, S}
function MOI.get(model::LPMixedForm{T}, ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, S}) where {T, S <: LinearBounds{T}}
    s = _sense(S)
    # FIXME `copy_constraint` needs a `Vector`
    return collect(MOIU.LazyMap{AFF{T, S}}(
        i -> AFF{T, S}(i), # FIXME `LazyMap` needs a `Function` so cannot just give `AFF{T, S}`
        Base.Iterators.Filter(1:size(model.A, 1)) do i
            _constraint_bound_sense(model, i) == s
        end
    ))
end

const VBOUND{S} = MOI.ConstraintIndex{MOI.VariableIndex, S}
function MOI.get(model::LPMixedForm{T}, ::MOI.NumberOfConstraints{MOI.VariableIndex, S}) where {T, S <: LinearBounds{T}}
    s = _sense(S)
    return count(MOI.get(model, MOI.ListOfVariableIndices())) do vi
        MOI.get(model, BoundSense(), vi) == s
    end
end
function MOI.get(model::LPMixedForm{T}, ::MOI.ListOfConstraintIndices{MOI.VariableIndex, S}) where {T, S <: LinearBounds{T}}
    s = _sense(S)
    return collect(MOIU.LazyMap{VBOUND{S}}(
        Base.Iterators.Filter(MOI.get(model, MOI.ListOfVariableIndices())) do vi
            MOI.get(model, BoundSense(), vi) == s
        end
    ) do vi
        VBOUND{S}(vi.value)
    end)
end
function MOI.get(::LPMixedForm, ::MOI.ConstraintFunction, ci::VBOUND)
    return MOI.VariableIndex(ci.value)
end
function MOI.get(model::LPMixedForm, ::MOI.ConstraintSet, ci::VBOUND)
    return _bound_set(model.v_lb[ci.value], model.v_ub[ci.value])
end

"""
    LPForm{T, AT<:AbstractMatrix{T}, VT <: AbstractVector{T}}

Represents a problem of the form:
```
sense ⟨c, x⟩
s.t.  c_lb <= Ax <= c_ub
      v_lb <=  x <= v_ub
```
"""
struct LPForm{T, AT<:AbstractMatrix{T}, VT <: AbstractVector{T}} <: LPMixedForm{T} #, V<:AbstractVector{T} #, M<:AbstractMatrix{T}}
    sense::MOI.OptimizationSense
    c::VT
    A::AT
    c_lb::VT
    c_ub::VT
    v_lb::VT
    v_ub::VT
end

function _constraint_bound_sense(model::LPForm, i)
    return _bound_sense(model.c_lb[i], model.c_ub[i])
end
function MOI.get(model::LPForm, ::MOI.ConstraintSet, ci::AFF)
    return _bound_set(model.c_lb[ci.value], model.c_ub[ci.value])
end

"""
    LPSolverForm{T, AT<:AbstractMatrix{T}, VT<:AbstractVector{T}}

Represents a problem of the form:
```
sense ⟨c, x⟩
s.t. Ax senses b
     v_lb <=  x <= v_ub
```
"""
struct LPSolverForm{T, AT<:AbstractMatrix{T}, VT<:AbstractVector{T}} <: LPMixedForm{T}
    sense::MOI.OptimizationSense
    c::VT
    A::AT
    b::VT
    senses::Vector{ConstraintSense}
    v_lb::VT
    v_ub::VT
end

function _constraint_bound_sense(model::LPSolverForm, i)
    return model.senses[i]
end
function MOI.get(model::LPSolverForm, ::MOI.ConstraintSet, ci::AFF)
    s = model.senses[ci.value]
    β = model.b[ci.value]
    if s == GREATER_THAN
        return MOI.GreaterThan(β)
    elseif s == LESS_THAN
        return MOI.LessThan(β)
    else
        s == EQUAL_TO || error("Invalid $s")
        return MOI.EqualTo(β)
    end
end

"""
    MILP{T, LP<:AbstractLPForm{T}}

A mixed-integer problem represented by a linear problem of type `LP`
and a vector indicating each `VariableType`.
"""
struct MILP{T, LP<:AbstractLPForm{T}}
    lp::LP
    variable_type::Vector{VariableType}
end


text:
struct ConstantBlackModel
    vol::Float64
    r::Float64
    q::Float64
end

struct TSBlackModel{S}
    surface::S
    r::Float64
    q::Float64
end

struct LocalVolatilityModel{S}
    surface::S
    r::Float64
    q::Float64
end


text:
@testset "wclogp" begin

@testset "wclogp" begin
    # Wildman and Crippen 1999, Table.2
    mol1 = smilestomol("C=1C=CC=C(OC)C=1O")
    @test wclogptype(mol1) == [
        :C18, :C18, :C18, :C18, :C23, :O4, :C3, :C23, :O2
    ]
    @test wclogphydrogentype(mol1) == [
        :H1, :H1, :H1, :H1, :undef, :undef, :H1, :undef, :H2
    ]
    @test wclogp(mol1) == 1.40

    mol2 = smilestomol("C1=CC=CC=C1C2=CC=CC=N2")
    @test wclogptype(mol2) == [
        :C18, :C18, :C18, :C18, :C18, :C20, :C20,
        :C18, :C18, :C18, :C18, :N11
    ]
    @test wclogphydrogentype(mol2) == [
        :H1, :H1, :H1, :H1, :H1, :undef, :undef,
        :H1, :H1, :H1, :H1, :undef
    ]
    @test wclogp(mol2) == 2.75

    # Test molecules
    TESTMOL_DIR = joinpath(dirname(@__FILE__), "..", "assets", "test")
    aliphatic = sdftomol(open(joinpath(TESTMOL_DIR, "wctype_C_alip.mol")))
    @test wclogptype(aliphatic) == [
        :C2, :C1, :C7, :C4, :N2, :C3, :N1, :C7, :C5, :O9,
        :C27, :Me1, :C27, :C27, :C27, :undef, :C6, :C26, :C21, :C21,
        :N12, :C18, :C21, :C21, :C11, :O3, :C12, :C1, :C1, :C1,
        :C9, :C10, :C1
    ]

    aromatic = sdftomol(open(joinpath(TESTMOL_DIR, "wctype_C_arom.mol")))
    @test wclogptype(aromatic) == [
        :C21, :C23, :C22, :C24, :C19, :C19, :C14, :C20, :C25, :O1,
        :O8, :O2, :N3, :S1, :C20, :C15, :C16, :C17, :C13, :C18,
        :Cl, :Br, :I, :C8, :P, :F
    ]

    nitrogen = sdftomol(open(joinpath(TESTMOL_DIR, "wctype_N.mol")))
    @test wclogptype(nitrogen) == [
        :N11, :C21, :C22, :C18, :C22, :C22, :N4, :C4, :N8, :C4,
        :N6, :N7, :N5, :C3, :N13, :C3, :C3, :C3, :C3, :N6,
        :N10, :C7, :N9, :N14, :N14
    ]

    oands = sdftomol(open(joinpath(TESTMOL_DIR, "wctype_OS.mol")))
    @test wclogptype(oands) == [
        :S2, :C4, :N13, :P, :C4, :O6, :O5, :O7, :C5, :O9,
        :O3, :C5, :O4, :C23, :S3, :C23, :C24, :C21, :O2, :O11,
        :C5, :O12, :O10, :O4, :S1, :O2, :O2
    ]
end

end # wclogp


text:
using Documenter, ScoreDrivenModels

makedocs(
    modules = [ScoreDrivenModels],
    doctest  = false,
    clean    = true,
    format   = Documenter.HTML(mathengine = Documenter.MathJax()),
    sitename = "ScoreDrivenModels.jl",
    authors  = "Guilherme Bodin and Raphael Saavedra",
    pages   = [
        "Home" => "index.md",
        "manual.md",
        "examples.md"
    ]
)

deploydocs(
    repo = "github.com/LAMPSPUC/ScoreDrivenModels.jl.git",
)


text:
using OpenQuantumTools, OrdinaryDiffEq, Plots
# define the Hamiltonian
H = DenseHamiltonian([(s)->1.0], [σz], unit=:ħ)
# define the initial state
u0 = PauliVec[1][1]*PauliVec[1][1]'
# define the Lindblad operator
# the rate and Lindblad operator can also be time-dependent functions
lind = Lindblad(0.1, σz)
# combine them into an Annealing object
annealing = Annealing(H, u0, interactions = InteractionSet(lind))

# define total annealing/evolution time
tf = 10
# solve the Lindblad equation
sol = solve_lindblad(annealing, 10, alg=Tsit5());

t_axis = range(0, 10, length=100)
bloch_vector = []
for t in t_axis
    # matrix_decompose projects a matrix onto a list of basis elements
    push!(bloch_vector, 2*real.(matrix_decompose(sol(t), [σx, σy, σz])))
end

off_diag = []
for t in t_axis
    push!(off_diag, abs(sol(t)[1,2]))
end

plot(t_axis, [c[1] for c in bloch_vector], label="X", linewidth=2)
plot!(t_axis, [c[2] for c in bloch_vector], label="Y", linewidth=2)
plot!(t_axis, [c[3] for c in bloch_vector], label="Z", linewidth=2)
xlabel!("t (ns)")
ylabel!("Bloch Vector")

plot(t_axis, off_diag, linewidth=2, label="ME")
plot!(t_axis, 0.5*exp.(-0.2*t_axis), linestyle=:dash, linewidth=3, label="Analytical")
xlabel!("t (ns)")
ylabel!("|ρ₀₁(t)|")

# For the quantum trajectories method, the u0 supplied to `Annealing` must be
# a state vector.
# We will show how to replace it with a pure state ensemble later
# in this example
u0 = PauliVec[1][1]
lind = Lindblad(0.1, σz)
annealing = Annealing(H, u0, interactions = InteractionSet(lind))
tf = 10
prob = build_ensembles(annealing, tf, :lindblad)
# We ran each trajectory serially for the sake of simplicity. The user is encouraged to try parallel algorithms.
sol = solve(prob, Tsit5(), EnsembleSerial(), trajectories=1000, saveat=range(0,tf,length=100))

vec_norm = []
# this is the index of the trajectory you want to look at
idx = 2
for v in sol[idx].u
    push!(vec_norm, norm(v))
end
plot(sol[idx].t, vec_norm, linewidth=2, label="", xlabel="t (ns)", ylabel="‖ψ̃(t)‖")

t_axis = range(0,tf,length=100)
dataset = []
for t in t_axis
    temp = []
    for so in sol
        v = normalize(so(t))
        push!(temp, real.(v'*σx*v))
    end
    push!(dataset, temp)
end

x_mean = []
x_sem = []
for data in dataset
    t_mean = sum(data)/1000
    t_sem = sqrt(sum((x)->(x-t_mean)^2, data))/1000
    push!(x_mean, t_mean)
    push!(x_sem, t_sem)
end

scatter(t_axis, x_mean, marker=:d, yerror=2 * x_sem, label="1000 trajectories", markersize=6, ylabel="<X>", xlabel="t (ns)")
plot!(t_axis, [c[1] for c in bloch_vector], linewidth=2, label="direct solver")

dataset = []
for t in t_axis
    temp = []
    for so in sol[1:100]
        v = normalize(so(t))
        push!(temp, real.(v'*σx*v))
    end
    push!(dataset, temp)
end

x_mean = []
x_sem = []
for data in dataset
    t_mean = sum(data)/100
    t_sem = sqrt(sum((x)->(x-t_mean)^2, data))/100
    push!(x_mean, t_mean)
    push!(x_sem, t_sem)
end

plot!(t_axis, x_mean, ribbon=2*x_sem, label="100 trajectories", markersize=6, ylabel="<X>", xlabel="t (ns)")

# PuliVec[1][1] is the plus state and PauliVec[1][2] is the minus state
# The first argument is a list of corresponding probabilities of the
# pure states in the second argument.
E = EᵨEnsemble([0.7, 0.3], [PauliVec[1][1], PauliVec[1][2]])

function prob_func(prob,i,repeat)
  prob.u0 .= sample_state_vector(E)
  prob
end

u0 = PauliVec[1][1]
lind = Lindblad(0.1, σz)
annealing = Annealing(H, u0, interactions = InteractionSet(lind))
tf = 10
prob = build_ensembles(annealing, tf, :lindblad, prob_func=prob_func)
sol = solve(prob, Tsit5(), EnsembleSerial(), trajectories=2000, saveat=range(0,tf,length=100))

initial_state_counter = zeros(2)
for so in sol
    if so.prob.u0 == PauliVec[1][1]
        initial_state_counter[1] += 1
    else
        initial_state_counter[2] += 1
    end
end
bar([0,1],initial_state_counter, label="", ylabel="Frequency", xticks=([0, 1], ["|+⟩","|-⟩"]))

t_axis = range(0,tf,length=100)
dataset = []
for t in t_axis
    temp = []
    for so in sol
        v = normalize(so(t))
        push!(temp, real.(v'*σx*v))
    end
    push!(dataset, temp)
end

x_mean = []
x_sem = []
for data in dataset
    t_mean = sum(data)/2000
    t_sem = sqrt(sum((x)->(x-t_mean)^2, data))/2000
    push!(x_mean, t_mean)
    push!(x_sem, t_sem)
end

# define Hamiltoian
H = DenseHamiltonian([(s)->1.0], [σz], unit=:ħ)
# define initial state
u0 = 0.7*PauliVec[1][1]*PauliVec[1][1]'+0.3*PauliVec[1][2]*PauliVec[1][2]'
# define Lindblad operator
lind = Lindblad(0.1, σz)
# combine them into an Annealing object
annealing = Annealing(H, u0, interactions = InteractionSet(lind))

# define total annealing/evolution time
tf = 10
# solve the Lindblad equation
sol = solve_lindblad(annealing, 10, alg=Tsit5());

t_axis = range(0, 10, length=100)
x_vector = []
for t in t_axis
    push!(x_vector, 2*real.(matrix_decompose(sol(t), [σx])))
end

scatter(t_axis, x_mean, marker=:d, yerror=2*x_sem, label="2000 trajectories", markersize=6, ylabel="<X>", xlabel="t (ns)")

plot!(t_axis, [c[1] for c in x_vector], linewidth=2, label="direct solver")

using HOQSTTutorials
HOQSTTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
ENV["GKSwstype"]="100" #src

# # Metaprogramming
#
# ## The ParticleGroup example

import Sobol
using Plots, LinearAlgebra

#     ParticleGroup{D,V}(n_particles, charge, mass)
#
# - `D` : number of dimension in physical space
# - `V` : number of dimension in phase space
# - `n` : number of particles
#
struct ParticleGroup{D,V}

    n_particles       :: Int64
    data              :: Array{Float64, 2}

    function ParticleGroup{D,V}(n) where {D, V}

        data = zeros( Float64, (D+V, n))
        new( n, data)
    end
end

#md # ---

#
# Set position of ith particle of p to x
#
@generated function set_x!( p :: ParticleGroup{D,V}, i, x :: Float64 ) where {D, V}

    :(p.data[1, i] = x)

end

#md # --

#
# Set position of ith particle of p to x when x is a vector
#
@generated function set_x!( p :: ParticleGroup{D,V}, i, x :: Vector{Float64} ) where {D, V}

     :(for j in 1:$D p.data[j, i] = x[j] end)

end

#md # ---

#
# Set velocity of ith particle of p to v
#
@generated function set_v!( p :: ParticleGroup{D,V}, i, v :: Float64 ) where {D, V}

    :(p.data[$D+1, i] = v)

end

#md # --

#
# Set velocity of ith particle of p to v
#
@generated function set_v!( p :: ParticleGroup{D,V}, i, v :: Vector{Float64} ) where {D, V}

    :(for j in 1:$V p.data[$D+j, i] = v[j] end)

end

#md # ---

#
# Get position of ith particle of p
#
@generated function get_x( p :: ParticleGroup{D,V}, i ) where {D, V}

     :(p.data[1:$D, i])

end

#
# Get velocity of ith particle of p
#
@generated function get_v( p :: ParticleGroup{D,V}, i ) where {D, V}

     :(p.data[$D+1:$D+$V, i])
end

#md # ---

#
# Sampling from a probability distribution to initialize Landau damping
#
function landau_sampling!( pg :: ParticleGroup{1,2}, alpha, kx )

    function newton(r)
        x0, x1 = 0.0, 1.0
        r *= 2π / kx
        while (abs(x1-x0) > 1e-12)
            p = x0 + alpha * sin( kx * x0) / kx
            f = 1 + alpha * cos( kx * x0)
            x0, x1 = x1, x0 - (p - r) / f
        end
        x1
    end

    s = Sobol.SobolSeq(2)
    nbpart = pg.n_particles

    for i=1:nbpart

        v = sqrt(-2 * log( (i-0.5)/nbpart))
        r1, r2 = Sobol.next!(s)
        θ = r1 * 2π
        set_x!(pg,  i, newton(r2))
        set_v!(pg,  i, [v * cos(θ), v * sin(θ)])
    end

end

#nb #+
#md # ---

n_particles = 10000
pg = ParticleGroup{1,2}( n_particles)
alpha, kx = 0.1, 0.5
landau_sampling!(pg, alpha, kx)

#md # --

#nb #+
xp = vcat([get_x(pg, i) for i in 1:pg.n_particles]...)
vp = vcat([get_v(pg, i) for i in 1:pg.n_particles]'...)

#md # ---
#nb #+
pp = plot(layout=(3,1))
histogram!(pp[1,1], xp, normalize=true, bins = 100, lab=:x)
plot!(pp[1,1], x -> (1+alpha*cos(kx*x))/(2π/kx), 0., 2π/kx, lab="")
histogram!(pp[2,1], vp[:,1], normalize=true, bins = 100, lab=:vx)
plot!(pp[2,1], v -> exp( - v^2 / 2) * 4 / π^2 , -6, 6, lab="")
histogram!(pp[3,1], vp[:,2], normalize=true, bins = 100, lab=:vy)
plot!(pp[3,1], v -> exp( - v^2 / 2) * 4 / π^2 , -6, 6, lab="")
#md savefig("particles.svg"); nothing # hide

#md # ![](particles.svg)

#md # ---

#nb #+
histogram2d(vp[:,1], vp[:,2], normalize=true, bins=100)
#md savefig("hist2d.svg"); nothing # hide

#md # ![](hist2d.svg)

@test true #src

#md # Next: [GEMPIC](/Numkin2019/05/build/index.html)


text:
using CellListMap
using StaticArrays
import Random

#
# In this test we compute the complete neighbor list of particles, meaning all the pairs
# that are within the cutoff distance. This function is implemented in the
#
# CellListMap.neighborlist
#
# routine. Here we illustrate how is is this implementation, to facilitate customization.
#
function neighborlist(;N=100_000,parallel=true,x=nothing)

    # Number of particles, sides and cutoff
    sides = @SVector [250,250,250]
    cutoff = 10.
    box = Box(sides, cutoff)

    # Particle positions
    Random.seed!(321)
    if x === nothing
        x = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
    end

    # Initialize auxiliary linked lists
    cl = CellList(x, box, parallel=parallel)

    # Function to be evalulated for each pair: push pair if d<cutoff
    function push_pair!(i, j, d2, pairs, cutoff)
        d = sqrt(d2)
        if d < cutoff
            push!(pairs, (i, j, d))
        end
        return pairs
    end

    # Reduction function
    function reduce_pairs(pairs, pairs_threaded)
        for i in 1:length(pairs_threaded)
            append!(pairs, pairs_threaded[i])
        end
        return pairs
    end

    # Initialize empty output
    pairs = Tuple{Int,Int,Float64}[]

    # Run pairwise computation
    map_pairwise!(
        (x, y, i, j, d2, pairs) -> push_pair!(i, j, d2, pairs, cutoff),
        pairs,box,cl,
        reduce=reduce_pairs,
        parallel=parallel
    )

    return pairs
end


text:
struct DistributedMultiFieldFEFunction{A,B,C} <: GridapType
  field_fe_fun::A
  part_fe_fun::B
  free_values::C
  function DistributedMultiFieldFEFunction(
    field_fe_fun::AbstractVector{<:DistributedSingleFieldFEFunction},
    part_fe_fun::AbstractPData{<:MultiFieldFEFunction},
    free_values::AbstractVector)
    A = typeof(field_fe_fun)
    B = typeof(part_fe_fun)
    C = typeof(free_values)
    new{A,B,C}(field_fe_fun,part_fe_fun,free_values)
  end
end

function FESpaces.get_free_dof_values(uh::DistributedMultiFieldFEFunction)
  uh.free_values
end

local_views(a::DistributedMultiFieldFEFunction) = a.part_fe_fun
MultiField.num_fields(m::DistributedMultiFieldFEFunction) = length(m.field_fe_fun)
Base.iterate(m::DistributedMultiFieldFEFunction) = iterate(m.field_fe_fun)
Base.iterate(m::DistributedMultiFieldFEFunction,state) = iterate(m.field_fe_fun,state)
Base.getindex(m::DistributedMultiFieldFEFunction,field_id::Integer) = m.field_fe_fun[field_id]

struct DistributedMultiFieldFESpace{A,B,C,D} <: DistributedFESpace
  field_fe_space::A
  part_fe_space::B
  gids::C
  vector_type::Type{D}
  function DistributedMultiFieldFESpace(
    field_fe_space::AbstractVector{<:DistributedSingleFieldFESpace},
    part_fe_space::AbstractPData{<:MultiFieldFESpace},
    gids::PRange,
    vector_type::Type{D}) where D
    A = typeof(field_fe_space)
    B = typeof(part_fe_space)
    C = typeof(gids)
    new{A,B,C,D}(field_fe_space,part_fe_space,gids,vector_type)
  end
end

local_views(a::DistributedMultiFieldFESpace) = a.part_fe_space
MultiField.num_fields(m::DistributedMultiFieldFESpace) = length(m.field_fe_space)
Base.iterate(m::DistributedMultiFieldFESpace) = iterate(m.field_fe_space)
Base.iterate(m::DistributedMultiFieldFESpace,state) = iterate(m.field_fe_space,state)
Base.getindex(m::DistributedMultiFieldFESpace,field_id::Integer) = m.field_fe_space[field_id]

function FESpaces.get_vector_type(fs::DistributedMultiFieldFESpace)
  fs.vector_type
end

function FESpaces.get_free_dof_ids(fs::DistributedMultiFieldFESpace)
  fs.gids
end

function MultiField.restrict_to_field(
  f::DistributedMultiFieldFESpace,free_values::PVector,field::Integer)
  values = map_parts(f.part_fe_space,free_values.values) do u,x
    restrict_to_field(u,x,field)
  end
  gids = f.field_fe_space[field].gids
  PVector(values,gids)
end

function FESpaces.FEFunction(
  f::DistributedMultiFieldFESpace,x::AbstractVector,isconsistent=false)
  free_values = change_ghost(x,f.gids)
  # This will cause also the single-field components to be consistent
  local_vals = consistent_local_views(free_values,f.gids,isconsistent)
  part_fe_fun = map_parts(FEFunction,f.part_fe_space,local_vals)
  field_fe_fun = DistributedSingleFieldFEFunction[]
  for i in 1:num_fields(f)
    free_values_i = restrict_to_field(f,free_values,i)
    fe_space_i = f.field_fe_space[i]
    fe_fun_i = FEFunction(fe_space_i,free_values_i,true)
    push!(field_fe_fun,fe_fun_i)
  end
  DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values)
end

function FESpaces.EvaluationFunction(
  f::DistributedMultiFieldFESpace,x::AbstractVector,isconsistent=false)
  free_values = change_ghost(x,f.gids)
  # This will cause also the single-field components to be consistent
  local_vals = consistent_local_views(free_values,f.gids,false)
  part_fe_fun = map_parts(EvaluationFunction,f.part_fe_space,local_vals)
  field_fe_fun = DistributedSingleFieldFEFunction[]
  for i in 1:num_fields(f)
    free_values_i = restrict_to_field(f,free_values,i)
    fe_space_i = f.field_fe_space[i]
    fe_fun_i = EvaluationFunction(fe_space_i,free_values_i)
    push!(field_fe_fun,fe_fun_i)
  end
  DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values)
end

function FESpaces.interpolate(objects,fe::DistributedMultiFieldFESpace)
  free_values = zero_free_values(fe)
  interpolate!(objects,free_values,fe)
end

function FESpaces.interpolate!(objects,free_values::AbstractVector,fe::DistributedMultiFieldFESpace)
  local_vals = consistent_local_views(free_values,fe.gids,true)
  part_fe_fun = map_parts(local_vals,local_views(fe)) do x,f
    interpolate!(objects,x,f)
  end
  field_fe_fun = DistributedSingleFieldFEFunction[]
  for i in 1:num_fields(fe)
    free_values_i = restrict_to_field(fe,free_values,i)
    fe_space_i = fe.field_fe_space[i]
    fe_fun_i = FEFunction(fe_space_i,free_values_i)
    push!(field_fe_fun,fe_fun_i)
  end
  DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values)
end

struct DistributedMultiFieldFEBasis{A,B} <: GridapType
  field_fe_basis::A
  part_fe_basis::B
  function DistributedMultiFieldFEBasis(
    field_fe_basis::AbstractVector{<:DistributedCellField},
    part_fe_basis::AbstractPData{<:MultiFieldCellField})
    A = typeof(field_fe_basis)
    B = typeof(part_fe_basis)
    new{A,B}(field_fe_basis,part_fe_basis)
  end
end

local_views(a::DistributedMultiFieldFEBasis) = a.part_fe_basis
MultiField.num_fields(m::DistributedMultiFieldFEBasis) = length(m.field_fe_basis)
Base.iterate(m::DistributedMultiFieldFEBasis) = iterate(m.field_fe_basis)
Base.iterate(m::DistributedMultiFieldFEBasis,state) = iterate(m.field_fe_basis,state)
Base.getindex(m::DistributedMultiFieldFEBasis,field_id::Integer) = m.field_fe_basis[field_id]

function FESpaces.get_fe_basis(f::DistributedMultiFieldFESpace)
  part_mbasis = map_parts(get_fe_basis,f.part_fe_space)
  field_fe_basis = DistributedCellField[]
  for i in 1:num_fields(f)
    basis_i = map_parts(b->b[i],part_mbasis)
    bi = DistributedCellField(basis_i)
    push!(field_fe_basis,bi)
  end
  DistributedMultiFieldFEBasis(field_fe_basis,part_mbasis)
end

function FESpaces.get_trial_fe_basis(f::DistributedMultiFieldFESpace)
  part_mbasis = map_parts(get_trial_fe_basis,f.part_fe_space)
  field_fe_basis = DistributedCellField[]
  for i in 1:num_fields(f)
    basis_i = map_parts(b->b[i],part_mbasis)
    bi = DistributedCellField(basis_i)
    push!(field_fe_basis,bi)
  end
  DistributedMultiFieldFEBasis(field_fe_basis,part_mbasis)
end

function FESpaces.TrialFESpace(objects,a::DistributedMultiFieldFESpace)
  TrialFESpace(a,objects)
end

function FESpaces.TrialFESpace(a::DistributedMultiFieldFESpace,objects)
  f_dspace_test = a.field_fe_space
  f_dspace = map( arg -> TrialFESpace(arg[1],arg[2]), zip(f_dspace_test,objects) )
  f_p_space = map(local_views,f_dspace)
  v(x...) = collect(x)
  p_f_space = map_parts(v,f_p_space...)
  p_mspace = map_parts(MultiFieldFESpace,p_f_space)
  gids = a.gids
  vector_type = a.vector_type
  DistributedMultiFieldFESpace(f_dspace,p_mspace,gids,vector_type)
end

# Factory

function MultiField.MultiFieldFESpace(
  f_dspace::Vector{<:DistributedSingleFieldFESpace})
  f_p_space = map(local_views,f_dspace)
  v(x...) = collect(x)
  p_f_space = map_parts(v,f_p_space...)
  p_mspace = map_parts(MultiFieldFESpace,p_f_space)
  gids = generate_multi_field_gids(f_dspace,p_mspace)
  vector_type = _find_vector_type(p_mspace,gids)
  DistributedMultiFieldFESpace(f_dspace,p_mspace,gids,vector_type)
end

function generate_multi_field_gids(
  f_dspace::Vector{<:DistributedSingleFieldFESpace},
  p_mspace::AbstractPData{<:MultiFieldFESpace})

  p_lids = map_parts(mspace->collect(get_free_dof_ids(mspace)),p_mspace)
  p_1lid_lid = map_parts(p_mspace,p_lids) do mspace, lids
    restrict_to_field(mspace,lids,1)
  end
  f_p_flid_lid = [p_1lid_lid]
  for f in 2:length(f_dspace)
    p_flid_lid = map_parts(p_mspace,p_lids) do mspace, lids
      restrict_to_field(mspace,lids,f)
    end
    push!(f_p_flid_lid,p_flid_lid)
  end
  f_frange = map(get_free_dof_ids,f_dspace)
  gids = generate_multi_field_gids(f_p_flid_lid,f_frange)
end

function generate_multi_field_gids(
  f_p_flid_lid::AbstractVector{<:AbstractPData{<:AbstractVector}},
  f_frange::AbstractVector{<:PRange})

  f_p_fiset = map(local_views,f_frange)

  v(x...) = collect(x)
  p_f_fiset = map_parts(v,f_p_fiset...)
  p_f_flid_lid = map_parts(v,f_p_flid_lid...)

  # Find the first gid of the multifield space in each part
  ngids = sum(map(num_gids,f_frange))
  p_noids = map_parts(f_fiset->sum(map(num_oids,f_fiset)),p_f_fiset)
  p_part = get_part_ids(p_noids)
  p_firstgid = xscan(+,p_noids,init=1)

  # Distributed gids to owned dofs
  p_lid_gid, p_lid_part = map_parts(
    p_f_flid_lid, p_f_fiset, p_firstgid) do f_flid_lid, f_fiset, firstgid
    nlids = sum(map(length,f_flid_lid))
    lid_gid = zeros(Int,nlids)
    lid_part = zeros(Int32,nlids)
    nf = length(f_fiset)
    gid = firstgid
    for f in 1:nf
      fiset = f_fiset[f]
      flid_lid = f_flid_lid[f]
      part = fiset.part
      for foid in 1:num_oids(fiset)
        flid = fiset.oid_to_lid[foid]
        lid = flid_lid[flid]
        lid_part[lid] = part
        lid_gid[lid] = gid
        gid += 1
      end
    end
    lid_gid,lid_part
  end

  # Now we need to propagate to ghost
  # to this end we use the already available
  # communicators in each of the single fields
  # We cannot use the cell wise dof like in the old version
  # since each field can be defined on an independent mesh.
  f_aux = map(frange->PVector{Int}(undef,frange),f_frange)
  propagate_to_ghost_multifield!(p_lid_gid,f_aux,f_p_flid_lid,f_p_fiset)
  propagate_to_ghost_multifield!(p_lid_part,f_aux,f_p_flid_lid,f_p_fiset)

  # Setup IndexSet
  p_iset = map_parts(IndexSet,p_part,p_lid_gid,p_lid_part)

  # Merge neighbors
  function merge_neigs(f_neigs)
    dict = Dict{Int32,Int32}()
    for f in length(f_neigs)
      for neig in f_neigs[f]
        dict[neig] = neig
      end
    end
    collect(keys(dict))
  end
  f_p_parts_snd = map(i->i.exchanger.parts_snd,f_frange)
  f_p_parts_rcv = map(i->i.exchanger.parts_rcv,f_frange)
  p_f_parts_snd = map_parts(v,f_p_parts_snd...)
  p_f_parts_rcv = map_parts(v,f_p_parts_rcv...)
  p_neigs_snd = map_parts(merge_neigs,p_f_parts_snd)
  p_neigs_rcv = map_parts(merge_neigs,p_f_parts_rcv)

  # Setup exchanger
  exchanger = Exchanger(p_iset,p_neigs_snd,p_neigs_rcv)

  # Setup the range
  ran = PRange(ngids,p_iset,exchanger)

  ran
end

function propagate_to_ghost_multifield!(
  p_lid_gid,f_gids,f_p_flid_lid,f_p_fiset)
  # Loop over fields
  nf = length(f_gids)
  for f in 1:nf
    # Write data into owned in single-field buffer
    gids = f_gids[f]
    p_flid_gid = gids.values
    p_flid_lid = f_p_flid_lid[f]
    p_fiset = f_p_fiset[f]
    map_parts(
      p_flid_gid,p_flid_lid,p_lid_gid,p_fiset) do flid_gid,flid_lid,lid_gid,fiset
      for foid in 1:num_oids(fiset)
        flid = fiset.oid_to_lid[foid]
        lid = flid_lid[flid]
        flid_gid[flid] = lid_gid[lid]
      end
    end
    # move to ghost
    exchange!(gids)
    # write again into multifield array on ghost ids
    map_parts(
      p_flid_gid,p_flid_lid,p_lid_gid,p_fiset) do flid_gid,flid_lid,lid_gid,fiset
      for fhid in 1:num_hids(fiset)
        flid = fiset.hid_to_lid[fhid]
        lid = flid_lid[flid]
        lid_gid[lid] = flid_gid[flid]
      end
    end
  end
end


text:
function parser(m::JuMP.Model)
    map::Map = Map()

    variable_list = JuMP.all_variables(m)
    @assert all(v -> JuMP.has_lower_bound(v) && JuMP.has_upper_bound(v), variable_list) "All variables must be bounded"

    for var in variable_list
        lower_bound = JuMP.lower_bound(var)
        upper_bound = JuMP.upper_bound(var)

        scale = upper_bound - lower_bound + 1

        num_created_city::Int64 = Int(ceil(log2(scale)))

        print(num_created_city)
    end

end