text:
include("../utils.jl")

function quick_sort!(A, start, stop)
    if start < stop
        mid = partition(A, start, stop)
        quick_sort!(A, start, mid - 1)
        quick_sort!(A, mid + 1, stop)
    end
end

function partition(A, start, stop)
    pivot = A[stop]
    i = start - 1
    for j in start:(stop-1)
        if A[j] ≤ pivot
            i += 1
            exchange!(A, i, j)
        end
    end

    exchange!(A, i+1, stop)

    return i+1
end

function randomized_partition!(A, start, stop)
    pivot = rand(start:stop)
    exchange!(A, pivot, stop)
    return partition(A, start, stop)
end

function randomized_quick_sort!(A, start, stop)
    if start < stop
        mid = randomized_partition!(A, start, stop)
        randomized_quick_sort!(A, start, mid - 1)
        randomized_quick_sort!(A, mid + 1, stop)
    end
end

max_stars_repo_name: blallen/CodePractice

max_issues_repo_name: blallen/CodePractice

max_forks_repo_name: blallen/CodePractice


text:
## Least squares
ei,ej = copy.(eachcol(Int.(readdlm("undir-graph-14.edges"))))
xy = readdlm("undir-graph-14.xy")'
n = max(maximum(ei),maximum(ej))
A = sparse(ei,ej,1,n,n)
A = A+A'
D = Diagonal(vec(sum(A,dims=1)))
m = length(ei)
B = sparse(1:m,ei,1,m,n) - sparse(1:m,ej,1,m,n)
##
idmat = sparse(1.0I,n,n)
known = [4,9]
unknown = Set(1:n)
for v in known
  delete!(unknown, v)
end
unknown = collect(unknown)
IK = idmat[known,:]
IU = idmat[unknown,:]
labels = ones(length(known))
##
M = [B; IK; sqrt(0.1)*IU]
b = [zeros(m); labels; zeros(size(IU,1))]
##
x = M\b
##
pyplot(size=(300,300))
p = graphplot(ei, ej, x =xy[1,:], y=xy[2,:],
  markercolor=:white, markerstrokecolor=:black,
  markersize=14, linecolor=3, linealpha=1.0, linewidth=1.0,
  marker_z = x,
  axis_buffer=0.02, background=:transparent,dpi=300)
for i=1:size(A,1)
  annotate!(xy[1,i],xy[2,i],text("$i",9,:black))
end
plot!()
savefig("least-squares-undir.pdf")

ext: jl


text:
using Roots, Dierckx, LsqFit, CoolProp, PyPlot, JLD
import JSON
include("../functions/rel_perms.jl")
include("../functions/frac_flow_funcs.jl")

# read the data files
L_core = 0.05
D_core = 0.03
poros  = 0.256
A_core  = π*D_core^2/4
PV_core = L_core*A_core*poros

swi     = 0.334
soi     = 1-swi
T_c     = 55.6
mu_oil = 0.35e-3
pv_inj_R_exp = [0.0, 0.029, 0.152, 0.274, 1.195, 2.177, 3.527, 4.939, 6.289, 7.700, 9.051]
R_oil_exp = [0.0, 0.043, 0.236, 0.430, 0.534, 0.564, 0.584, 0.593, 0.607, 0.610, 0.610]
R_oil_exp_int = Spline1D(pv_inj_R_exp, R_oil_exp, k=1)
pv_inj_dp = [0.08, 0.18, 0.24, 0.43, 0.61, 0.86, 1.43, 2.20, 3.02, 3.96, 4.90,
5.84, 6.80, 7.73, 8.73]
dp_exp = [0.001098, 0.486813, 0.53956, 0.486813, 0.434065, 0.381318, 0.334065, 0.302197,
0.271428, 0.261538, 0.247252, 0.23956, 0.231868, 0.224175, 0.221978]*20*1e5

pv_inj_exp = pv_inj_dp[:] # same as the time steps for the dp
R_oil = R_oil_exp_int(pv_inj_exp) # just to have both dp and R for the same time steps

pv_inj = pv_inj_exp[end]
injection_rate = 0.001/(3600*24) # m^3/day
recovery_final = R_oil_exp[end]

so_final = (1-recovery_final)*soi # final oil saturation

t_sec = pv_inj_exp*PV_core/injection_rate

# conversion of the input data to the right unit
u_inj   = injection_rate/A_core # Darcy velocity
T_K = T_c + 273.15
perm_ave = 0.001e-12 # m^2 permeability
mu_water = PropsSI("viscosity", "T", T_K, "P", 5e5, "H2O")

# define the objective functions
# param = [sor, swc, kro0, krw0, no, nw]
function error_calc(param, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp, dp_exp=0.0)
  (sor, swc, kro0, krw0, no, nw) = param
  (xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R)=
  frac_flow_wf(muw=muw, muo=muo, ut=ut, phi=phi,
    k=k, swc=swc, sor=sor, kro0=kro0, no=no, krw0=krw0,
    nw=nw, sw0=sw0, sw_inj=sw_inj, L=L, pv_inj=pv_inj)
  R_int = Spline1D(t, R, k=1)
  dp_int = Spline1D(t, p_inj, k=1)
  error_vals_R = R_int(t_exp) - R_exp
  error_vals = error_vals_R
  if dp_exp != 0.0
      error_vals_dp = (dp_int(t_exp) - dp_exp)/maximum(dp_exp) # normalize
      error_vals = [error_vals_R; error_vals_dp;]
  end
  return error_vals
end

function error_calc_lsq(param, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp, dp_exp=0.0)
  (sor, swc, kro0, krw0, no, nw) = param
  (xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R)=
  frac_flow_wf(muw=muw, muo=muo, ut=ut, phi=phi,
    k=k, swc=swc, sor=sor, kro0=kro0, no=no, krw0=krw0,
    nw=nw, sw0=sw0, sw_inj=sw_inj, L=L, pv_inj=pv_inj)
  R_int = Spline1D(t, R, k=1)
  dp_int = Spline1D(t, p_inj, k=1)
  error_vals_R = R_int(t_exp)
  error_vals = error_vals_R
  if dp_exp != 0.0
      error_vals_dp = dp_int(t_exp) # normalize
      error_vals = [error_vals_R; error_vals_dp;]
  end
  return error_vals
end

model = (t_exp, param) -> error_calc_lsq(param, mu_water, mu_oil, u_inj, poros, perm_ave, swi, 1.0, L_core, pv_inj, t_sec, R_oil, dp_exp)

function error_calculation(param, param_ind, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj,
  t_exp, R_exp, dp_exp; grad_calc=true)
  eps1 = 1e-8
  error_vals = error_calc(param, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp, dp_exp)
  if grad_calc
    jac  = zeros(length(error_vals), length(param_ind))
  else
    jac = zeros(0)
  end
  if grad_calc
    for j in eachindex(param_ind)
      param2 = copy(param)
      param2[param_ind[j]]+=eps1
      error_val2 = error_calc(param2, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp, dp_exp)
      jac[:, j] = (error_val2 - error_vals) / eps1
    end
  end
  return error_vals, jac
end

function objective_function(param::Vector{Float64}, grad::Vector{Float64}, param_all,
  param_ind, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp, dp_exp; w=1.0)
  param2 = copy(param_all)
  param2[param_ind] = param
  grad_calc = false
  if length(grad)>0
    grad_calc = true
  end
  error_val, jac = error_calculation(param2, param_ind, muw, muo, ut, phi, k, sw0,
    sw_inj, L, pv_inj, t_exp, R_exp, dp_exp, grad_calc = grad_calc)
  obj_func_value = sum(abs2, w.*error_val)
  if length(grad)>0
      grad[:] = 2.0*sum(w.*error_val.*jac, 1)
  end
  return obj_func_value
end

function plot_results(param, muw, muo, ut, phi, k, sw0, sw_inj, L, pv_inj, t_exp, R_exp, dp_exp)
  (sor, swc, kro0, krw0, no, nw) = param
  (xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R)=
  frac_flow_wf(muw=muw, muo=muo, ut=ut, phi=phi,
    k=k, swc=swc, sor=sor, kro0=kro0, no=no, krw0=krw0,
    nw=nw, sw0=sw0, sw_inj=sw_inj, L=L, pv_inj=pv_inj)
  R_int = Spline1D(t, R, k=1)
  dp_int = Spline1D(t, p_inj, k=1)
  error_vals = R_int(t_exp) - R_exp
  t_plot = linspace(0.0, maximum(t_exp), 100)
  figure(1)
  plot(t_plot, R_int(t_plot), t_exp, R_exp, "o")
  xlabel("time [s]")
  ylabel("Recovery factor")
  figure(2)
  plot(t_plot, dp_int(t_plot))
  if dp_exp !=0.0
      plot(t_exp, dp_exp, "o")
  end
  xlabel("time [s]")
  ylabel("pressure drop [Pa]")
end

plot_quick = param -> plot_results(param, mu_water, mu_oil, u_inj, poros, perm_ave, swi, 1.0, L_core, pv_inj, t_sec, R_oil, dp_exp)
# Optimization problem:

# parameter values
# [sor, swc, kro0, krw0, no, nw]
param_all = [0.2, swi/3, 0.95, 0.6, 2.0, 2.0]
param_ind = [1, 2,3,4,5,6]
w = ones(2*length(R_oil))
# w[end]=2
# w[end-1]=2
# w[7] = 3
obj_fun = (param, grad) -> objective_function(param::Vector{Float64}, grad::Vector{Float64}, param_all,
  param_ind, mu_water, mu_oil, u_inj, poros, perm_ave, swi, 1.0, L_core, pv_inj, t_sec, R_oil, dp_exp, w=w)

# x_init = [0.2, 0.8, 0.4, 2, 2]
# x_lb = [0.1, 0.1, 0.1, 0.5, 0.5]
# x_ub = [0.5, 1.0, 1.0, 5.0, 5.0]

x_init = copy(param_all)
x_lb = [0.05, 0.05, 0.1, 0.1, 1.0, 1.0]
x_ub = [0.3,0.3, 1.0, 1.0, 4.0, 4.0]

fit = curve_fit(model, [t_sec;t_sec], [R_oil;dp_exp], w, x_init, lower = x_lb, upper = x_ub)

plot_quick(fit.param)

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
mutable struct PrimalData{T, TV1<:AbstractVector{T}, TV2, TM<:AbstractMatrix{T}}
    σ::TV1
    α::TV1 # Lower move limit
    β::TV1 # Upper move limit
    p0::TV1
    q0::TV1
    p::TM
    q::TM
    ρ::TV2
    r::TV2
    r0::T
    x0::TV1
    x::TV1 # Optimal value for x in dual iteration
    x1::TV1 # Optimal value for x in previous outer iteration
    x2::TV1 # Optimal value for x in previous outer iteration
    f_x::T # Function value at current iteration
    f_x_previous::T
    g::TV2 # Inequality values at current iteration
    ∇f_x::TV1 # Function gradient at current iteration
    ∇g::TM # Inequality gradients [var, ineq] at current iteration
end
GPUUtils.whichdevice(p::PrimalData) = whichdevice(p.x)

struct XUpdater{T, TV, TPD<:PrimalData{T, TV}}
    pd::TPD
end
function (xu::XUpdater)(λ::DualSolution{:CPU})
    @unpack x, p0, p, q0, q, σ, x1, α, β, ρ = xu.pd
    λρ = dot(λ.cpu, ρ)
    tmap!(x, 1:length(x)) do j
        λpj = @matdot(λ.cpu, p, j)
        λqj = @matdot(λ.cpu, q, j)
        getxj(λρ, λpj, λqj, j, p0, q0, σ, x1, α, β)
    end
    return
end
function (xu::XUpdater{T, TV})(λ::DualSolution{:GPU}) where {T, TV <: CuVector{T}}
    @unpack x, p0, p, q0, q, σ, x1, α, β, ρ = xu.pd
    λρ = dot(λ.cpu, ρ)
    togpu!(λ)
    λgpu = λ.gpu
    args = (x, λgpu, λρ, p, q, p0, q0, σ, x1, α, β)
    callkernel(dev, xupdater_kernel, args)
    CUDAdrv.synchronize(ctx)
    return
end
function xupdater_kernel(x, λgpu, λρ, p, q, p0, q0, σ, x1, α, β)
    j = @thread_global_index()
	offset = @total_threads()
    while j <= length(x)
        λpj = @matdot(λgpu, p, j)
        λqj = @matdot(λgpu, q, j)
        x[j] = getxj(λρ, λpj, λqj, j, p0, q0, σ, x1, α, β)
        j += offset
    end
    return
end

function getxj(λρ, λpj, λqj, j, p0, q0, σ, x1, α, β)
    lj1 = p0[j] + λpj + λρ*σ[j]/4
    lj2 = q0[j] + λqj + λρ*σ[j]/4

    αj, βj = α[j], β[j]
    Lj, Uj = minus_plus(x1[j], σ[j])

    Ujαj = Uj - αj
    αjLj = αj - Lj
    ljαj = lj1/Ujαj^2 - lj2/αjLj^2

    Ujβj = Uj - βj
    βjLj = βj - Lj
    ljβj = lj1/Ujβj^2 - lj2/βjLj^2

    fpj = sqrt(lj1)
    fqj = sqrt(lj2)
    xj = (fpj * Lj + fqj * Uj) / (fpj + fqj)
    xj = ifelse(ljαj >= 0, αj, ifelse(ljβj <= 0, βj, xj))

    return xj
end

# Primal problem functions
struct ConvexApproxGradUpdater{T, TV, TPD<:PrimalData{T}, TM <: Model{T, TV}}
    pd::TPD
    m::TM
end

function (gu::ConvexApproxGradUpdater{T, TV})() where {T, TV <: AbstractVector}
    @unpack pd, m = gu
    @unpack f_x, g, r, x, σ, x1, p0, q0, ∇f_x, p, q, ρ, ∇g = pd
    n = dim(m)
    r0 = f_x
    r0 -= tmapreduce(+, 1:n, init = zero(T)) do j
        getgradelement(x, σ, x1, p0, q0, ∇f_x, j)
    end
    for i in 1:length(constraints(m))
        r[i] = g[i]
        r[i] -= tmapreduce(+, 1:n, init = zero(T)) do j
            getgradelement(x, σ, p, q, ρ[i], ∇g, (j, i))
        end
    end
    pd.r0 = r0
end

function (gu::ConvexApproxGradUpdater{T, TV})() where {T, TV <: CuVector}
    @unpack pd, m = gu
    @unpack f_x, g, r, x, σ, x1, p0, q0, ∇f_x, p, q, ρ, ∇g = pd
    n = dim(m)
    r0 = compute_r0(f_x, x, σ, x1, p0, q0, ∇f_x)
    update_r!(r, g, x, σ, p, q, ρ, ∇g)
    pd.r0 = r0
end

function compute_r0(r0, x::AbstractVector{T}, σ, x1, p0, q0, ∇f_x, ::Val{blocksize} = Val(80), ::Val{threads} = Val(256)) where {T, blocksize, threads}
    result = similar(x, T, (blocksize,))
    args = (x, σ, x1, p0, q0, ∇f_x, result, Val(threads))
    @cuda blocks = blocksize threads = threads gradupdater_kernel1(args...)
    CUDAnative.synchronize()
    r0 -= sum(Array(result))
    return r0
end

function gradupdater_kernel1(x::AbstractVector{T}, σ, x1, p0, q0, ∇f_x, result, ::Val{LMEM}) where {T, LMEM}
    @mapreduce_block(j, length(x), +, T, LMEM, result, begin
        getgradelement(x, σ, x1, p0, q0, ∇f_x, j)
    end)

    return
end

function update_r!(r, g, x::AbstractVector{T}, σ, p, q, ρ, ∇g, ::Val{blocksize} = Val(80), ::Val{threads} = Val(256)) where {T, blocksize, threads}
    result = similar(x, T, (blocksize,))
    for i in 1:length(r)
        r[i] = g[i]
        args = (x, σ, p, q, i, ρ[i], ∇g, result, Val(threads))
        @cuda blocks = blocksize threads = threads gradupdater_kernel2(args...)
        CUDAnative.synchronize()
        r[i] -= sum(Array(result))
    end
    return
end

function gradupdater_kernel2(x::AbstractVector{T}, σ, p, q, i, ρi, ∇g, result, ::Val{LMEM}) where {T, LMEM}
    @mapreduce_block(j, length(x), +, T, LMEM, result, begin
        getgradelement(x, σ, p, q, ρi, ∇g, (j, i))
    end)

    return
end

function getgradelement(x::AbstractVector{T}, σ, x1, p0, q0, ∇f_x, j::Int) where {T}
    xj = x[j]
    σj = σ[j]
    Lj, Uj = minus_plus(xj, σj) # x == x1
    ∇fj = ∇f_x[j]
    abs2σj∇fj = abs2(σj)*∇fj
    (p0j, q0j) = ifelse(∇fj > 0, (abs2σj∇fj, zero(T)), (zero(T), -abs2σj∇fj))
    p0[j], q0[j] = p0j, q0j
    return (p0[j] + q0[j])/σj
end
function getgradelement(x::AbstractVector{T}, σ, p, q, ρi, ∇g, ji::Tuple) where {T}
    j, i = ji
    σj = σ[j]
    xj = x[j]
    Lj, Uj = minus_plus(xj, σj) # x == x1
    ∇gj = ∇g[j,i]
    abs2σj∇gj = abs2(σj)*∇gj
    (pji, qji) = ifelse(∇gj > 0, (abs2σj∇gj, zero(T)), (zero(T), -abs2σj∇gj))
    p[j,i], q[j,i] = pji, qji
    Δ = ρi*σj/4
    return (pji + qji + 2Δ)/σj
end

struct VariableBoundsUpdater{T, TV, TPD <: PrimalData{T}, TModel <: Model{T, TV}}
    pd::TPD
    m::TModel
    μ::T
end
function (bu::VariableBoundsUpdater{T, TV})() where {T, TV <: AbstractVector}
    @unpack pd, m, μ = bu
    @unpack α, β, x, σ = pd
    n = dim(m)
    αβ = StructArray{Tuple{T, T}}((α, β))
    tmap!(αβ, 1:n) do j
        xj = x[j]
        Lj, Uj = minus_plus(xj, σ[j]) # x == x1 here
        αj = max(Lj + μ * (xj - Lj), min(m, j))
        βj = min(Uj - μ * (Uj - xj), max(m, j))
        αj, βj
    end
    return
end

function (bu::VariableBoundsUpdater{T, TV})() where {T, TV <: CuVector}
    @unpack pd, m, μ = bu
    @unpack α, β, x, σ = pd
    @unpack box_max, box_min = m
    n = dim(m)
    args = (α, β, σ, x, box_max, box_min)
    callkernel(dev, bounds_kernel, args)
    CUDAdrv.synchronize(ctx)
    return
end

function bounds_kernel(α, β, σ, x, box_max, box_min)
    j = @thread_global_index()
    offset = @total_threads()
    while j <= length(σ)
        xj = x[j]
        Lj, Uj = minus_plus(xj, σ[j]) # x == x1 here
        αj = max(Lj + μ * (xj - Lj), box_min[j])
        βj = min(Uj - μ * (Uj - xj), box_max[j])
        α[j] = αj
        β[j] = βj
        j += offset
    end
    return
end

struct AsymptotesUpdater{T, TV <: AbstractVector{T}, TModel <: Model{T}}
    m::TModel
    σ::TV
    x::TV
    x1::TV
    x2::TV
    s_init::T
    s_incr::T
    s_decr::T
end

function (au::AsymptotesUpdater{T, TV})(k::Iteration) where {T, TV <: AbstractVector}
    @unpack σ, m, s_init, x, x1, x2, s_incr, s_decr = au
    if k.i == 1 || k.i == 2
        tmap!(σ, 1:dim(m)) do j
            s_init * (max(m, j) - min(m, j))
        end
    else
        tmap!(σ, 1:dim(m)) do j
            σj = σ[j]
            xj = x[j]
            x1j = x1[j]
            x2j = x2[j]
            d = ifelse((xj == x1j || x1j == x2j),
                σj, ifelse(xor(xj > x1j, x1j > x2j),
                σj * s_decr, σj * s_incr))
            diff = max(m, j) - min(m, j)
            _min = T(0.01)*diff
            _max = 10diff
            ifelse(d <= _min, _min, ifelse(d >= _max, _max, d))
        end
    end

    return
end

function (au::AsymptotesUpdater{T, TV})(k::Iteration) where {T, TV <: CuVector}
    @unpack σ, m, s_init, s_incr, s_decr, x, x1, x2 = au
    @unpack box_max, box_min = m

    if k.i == 1 || k.i == 2
        args = (σ, s_init, box_max, box_min)
        callkernel(dev, asymptotes_kernel1, args)
        CUDAdrv.synchronize(ctx)
    else
        args = (σ, x, x1, x2, box_max, box_min, s_incr, s_decr)
        callkernel(dev, asymptotes_kernel2, args)
        CUDAdrv.synchronize(ctx)
    end
end

function asymptotes_kernel1(σ, s_init, box_max, box_min)
    j = @thread_global_index()
    offset = @total_threads()
    while j <= length(σ)
        σ[j] = s_init * (box_max[j] - box_min[j])
        j += offset
    end
    return
end

function asymptotes_kernel2(σ::AbstractVector{T}, x, x1, x2, box_max, box_min, s_incr, s_decr) where {T}
    j = @thread_global_index()
    offset = @total_threads()
    while j <= length(σ)
        σj = σ[j]
        xj = x[j]
        x1j = x1[j]
        x2j = x2[j]
        d = ifelse((xj == x1j || x1j == x2j),
            σj, ifelse(xor(xj > x1j, x1j > x2j),
            σj * s_decr, σj * s_incr))
        diff = box_max[j] - box_min[j]
        _min = T(0.01)*diff
        _max = 10diff
        σ[j] = ifelse(d <= _min, _min, ifelse(d >= _max, _max, d))
        j += offset
    end
    return
end


text:
#=
Replication of Artificial Intelligence, Algorithmic Pricing, and Collusion
    by: Calvano, Calzolari, Denicolò (2020)
    at: https://www.aeaweb.org/articles?id=10.1257/aer.20190623
Code
    author: Matteo Courthoud
    date: 07/05/2021
    git: https://github.com/matteocourthoud
    myself: https://matteocourthoud.github.io/
=#

include("input/model.jl")
include("input/qlearning.jl")

# Init algorithm
game = model.init_game();

# Compute equilibrium
game_equilibrium = qlearning.simulate_game(game);


text:
# Test the pure Rayleigh case
@testset "pure_rayleigh" begin

    println("Testing pure rayleigh RT...")

    ϵ = 1e-4

    parameters = vSmartMOM.parameters_from_yaml("test_parameters/PureRayleighParameters.yaml")
    model = model_from_parameters(parameters)

    R_model = rt_run(model)

    R_test = [1.6376193549450170E-002  2.1561609277437792E-003  5.5387198966105483E-019 0.0 ;
              1.2818490637667997E-002  3.7825131902305691E-004  2.0685375352086151E-019 0.0 ;
              1.0839231562441865E-002 -3.0496438891765027E-005  4.2614463687432547E-021 0.0 ;
              9.4104223210835589E-003  2.9434545102941607E-004 -1.4472678739285088E-019 0.0 ;
              8.2598796504639114E-003  1.1181829243926966E-003  0.0000000000000000      0.0 ;
              7.3749583918955438E-003  2.3298093802174317E-003  0.0000000000000000      0.0 ;
              6.9126577200309845E-003  3.8960774035191144E-003  0.0000000000000000      0.0 ;
              7.2827593896516067E-003  5.9139825670394474E-003  0.0000000000000000      0.0 ;
              9.6437995452080226E-003  8.8885549319859258E-003  0.0000000000000000      0.0 ]

    @test isapprox(R_model[:,:,1], R_test, atol=ϵ)

end

@testset "compare against 6SV1" begin

    include("benchmarks/6SV1_R_trues.jl")
    R_modeled_all = zeros(6, 3, 3, 16);
    R_deltas_all = zeros(6, 3, 3, 16);
    parameters = parameters_from_yaml("benchmarks/6SV1_1.yaml");

    function test_against_6SV1(case_i, azs, szas, λ, τ, ρ)
        for sza_i in 1:length(szas)
            for az_i in 1:length(azs)
                R_true = R_trues[case_i]
                R_modeled = view(R_modeled_all, case_i, :, :, :)
                R_deltas = view(R_deltas_all, case_i, :, :, :)
                parameters.spec_bands = [1e7/λ (1e7/λ + 1)]
                parameters.vaz = repeat([azs[az_i]], 16)
                parameters.sza = szas[sza_i]
                parameters.brdf = [RadiativeTransfer.vSmartMOM.LambertianSurfaceScalar(ρ * π)]
                model = model_from_parameters(parameters);
                model.τ_rayl[1] .= τ
                R_modeled[sza_i, az_i, :] = vSmartMOM.rt_run(model, i_band=1)[:,1,1] / model.quad_points.μ₀
                R_deltas[sza_i, az_i, :] = abs.(R_true[sza_i][az_i] - R_modeled[sza_i, az_i, :]) ./ R_true[sza_i][az_i]
            end
        end
        return maximum(R_deltas_all[case_i,:,:,:])
    end

    ϵ = 0.006

    @test test_against_6SV1(1, [180, 90, 0], [23.0739, 53.1301, 78.4630], 530, 0.1, 0.0) < ϵ
    @test test_against_6SV1(2, [180, 90, 0], [0.0001, 36.8699, 66.4218], 530, 0.1, 0.25) < ϵ
    @test test_against_6SV1(3, [180, 90, 0], [0.0001, 36.8699, 66.4218], 440, 0.25, 0.0) < ϵ
    @test test_against_6SV1(4, [180, 90, 0], [23.0739, 53.1301, 78.4630], 440, 0.25, 0.25) < ϵ
    @test test_against_6SV1(5, [180, 90, 0], [23.0739, 53.1301, 78.4630], 360, 0.50, 0.0) < ϵ
    @test test_against_6SV1(6, [180, 90, 0], [0.0001, 36.8699, 66.4218], 360, 0.50, 0.25) < ϵ

end

@testset "compare against natraj paper" begin

    I_modeled_all = zeros(7, 16);
    I_deltas_all  = zeros(7, 16);
    Q_modeled_all = zeros(7, 16);
    Q_deltas_all  = zeros(7, 16);
    U_modeled_all = zeros(7, 16);
    U_deltas_all  = zeros(7, 16);

    μ = [0.02, 0.06, 0.10, 0.16, 0.20, 0.28, 0.32, 0.40, 0.52, 0.64, 0.72, 0.84, 0.92, 0.96, 0.98, 1.00]
    ϕs = collect(0.0:30.0:180.0)
    τ = 0.5

    include("benchmarks/natraj_trues.jl")
    parameters = parameters_from_yaml("benchmarks/natraj.yaml");

    parameters.spec_bands = [1e7/360.0 (1e7/360.0 + 1)]
    parameters.vza = acosd.(μ)
    parameters.sza = acosd.(0.2)

    for ϕ_i in 1:length(ϕs)
        parameters.vaz = repeat([ϕs[ϕ_i]], 16)
        model = model_from_parameters(parameters);
        model.τ_rayl[1] .= τ

        R = vSmartMOM.rt_run(model, i_band=1)
        @show size(R)

        I_modeled_all[ϕ_i,:] = R[:,1,1]
        Q_modeled_all[ϕ_i,:] = R[:,2,1]
        U_modeled_all[ϕ_i,:] = R[:,3,1]

        I_deltas_all[ϕ_i,:] = abs.(I_trues[:,ϕ_i] - I_modeled_all[ϕ_i,:]) ./ I_trues[:,ϕ_i]
        Q_deltas_all[ϕ_i,:] = abs.(Q_trues[:,ϕ_i] - Q_modeled_all[ϕ_i,:]) ./ Q_trues[:,ϕ_i]
        U_deltas_all[ϕ_i,:] = abs.(U_trues[:,ϕ_i] - U_modeled_all[ϕ_i,:]) ./ U_trues[:,ϕ_i]
    end

    ϵ = 0.008

    @test maximum(I_deltas_all) < 0.002
    @test maximum(Q_deltas_all[findall(i -> i >= 0.01, Q_modeled_all)]) < ϵ
    @test maximum(filter(!isnan, U_deltas_all[findall(i -> i >= 0.01, U_modeled_all)])) < ϵ

end

# # Test the O2 case
# @testset "O2 RT" begin

#     println("Testing RT with only O₂...")

#     parameters = vSmartMOM.parameters_from_yaml("helper/O2Parameters.yaml")

# end

# # Test 3 bands, 3 profiles case (O₂, CO₂, H₂O)
# @testset "pure_rayleigh" begin

#     println("Testing 3 bands, 3 profiles case (O₂, CO₂, H₂O)")

#     parameters = vSmartMOM.parameters_from_yaml("helper/ThreeBandsParameters.yaml")

# end


text:
#!/cvmfs/soft.computecanada.ca/easybuild/software/2017/avx512/Compiler/gcc7.3/julia/1.1.0/bin/julia
#SBATCH --time=00:05:00 # Running time of hours
#SBATCH --ntasks=4
#SBATCH --account=def-whitem

using Logging, Distributed

include("/home/mkschleg/mkschleg/ClusterManagers.jl/src/ClusterManagers.jl")

ClusterManagers.addprocs_slurm(4; exeflags=["--project=.", "--color=yes"], job_file_loc="test_loc")

@sync begin
    @async @sync for job_id in collect(1:100) @spawn begin
        println("Hello World $(job_id)")
    end
    end
end


text:
using Biomechanics, Test, Random

# include("timenormalize.jl")
@testset "Step metrics" begin include("stepmetrics.jl") end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(pretty_table)),NamedTuple{(:alignment, :crop), Tuple{Symbol, Symbol}},typeof(pretty_table),Vector{Any}})   # time: 0.05833952
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Vector{Float16}})   # time: 0.049024533
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Vector{Int16}})   # time: 0.045383085
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Vector{UInt16}})   # time: 0.043031204
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Matrix{Int8}})   # time: 0.03944839
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),AnsiTextCell})   # time: 0.03538433
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Vector{String}})   # time: 0.034589678
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Vector{UInt64}})   # time: 0.034464628
    Base.precompile(Tuple{typeof(_create_printing_recipe),Display,Int64,Int64,Int64,Int64,Int64,Vector{Int64},Vector{Int64},UnitRange{Int64},Vector{Int64},Vector{Int64},Int64,Int64,Symbol,Bool,Bool,Symbol})   # time: 0.033365928
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Vector{Float32}})   # time: 0.03203411
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Vector{Int8}})   # time: 0.029632116
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Vector{Int32}})   # time: 0.029488461
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Matrix{Int64}})   # time: 0.028275946
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Vector{UInt32}})   # time: 0.028247943
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Vector{UInt8}})   # time: 0.027401254
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Matrix{Float32}})   # time: 0.026978565
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Vector{Bool}})   # time: 0.026965057
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Vector{Float64}})   # time: 0.025837606
    Base.precompile(Tuple{typeof(_compute_row_fill_vectors),UnitRange{Int64},Int64,Symbol})   # time: 0.025818532
    Base.precompile(Tuple{typeof(_print_table_data),IOContext{IOBuffer},Display,Any,Matrix{Vector{String}},UnitRange{Int64},Vector{Int64},Int64,Vector{Int64},Vector{Int64},Vector{NTuple{4, Int64}},Vector{Int64},Vector{Symbol},NTuple{4, Char},Dict{Tuple{Int64, Int64}, Symbol},Symbol,Int64,Ref{Any},Bool,Bool,Bool,TextFormat,Crayon,Crayon,Crayon})   # time: 0.023430934
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),UInt32})   # time: 0.022836197
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Vector{Int64}})   # time: 0.02276013
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Matrix{Float64}})   # time: 0.02257849
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Symbol})   # time: 0.022379683
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Tuple{Int64, Int64}})   # time: 0.019844497
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Matrix{UInt64}})   # time: 0.016211173
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Vector{Char}})   # time: 0.01579448
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Matrix{Float16}})   # time: 0.015664909
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Matrix{UInt8}})   # time: 0.015432802
    Base.precompile(Tuple{typeof(_apply_alignment_anchor_regex!),Matrix{Vector{String}},Vector{Int64},UnitRange{Int64},UnitRange{Int64},Int64,Symbol,Dict{Int64, Symbol},Dict{Int64, Vector{Regex}},Dict{Tuple{Int64, Int64}, Symbol},Vector{Int64},Vector{Int64},Vector{Int64}})   # time: 0.015359414
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Matrix{UInt16}})   # time: 0.015224181
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Matrix{Char}})   # time: 0.015171672
    Base.precompile(Tuple{typeof(_create_printing_recipe),Display,Int64,Int64,Int64,Int64,Int64,Vector{Int64},Vector{Int64},Vector{Int64},Vector{Int64},Vector{Int64},Int64,Int64,Symbol,Bool,Bool,Symbol})   # time: 0.014925075
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Matrix{UInt32}})   # time: 0.014733393
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Matrix{Int16}})   # time: 0.01392427
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Matrix{String}})   # time: 0.013914018
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Matrix{Bool}})   # time: 0.013874153
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Matrix{Int32}})   # time: 0.01370403
    Base.precompile(Tuple{typeof(_apply_alignment_anchor_regex!),Matrix{Vector{String}},Vector{Int64},Vector{Int64},UnitRange{Int64},Int64,Symbol,Dict{Int64, Symbol},Dict{Int64, Vector{Regex}},Dict{Tuple{Int64, Int64}, Symbol},Vector{Int64},Vector{Int64},Vector{Int64}})   # time: 0.01366833
    Base.precompile(Tuple{typeof(_apply_alignment_anchor_regex!),Matrix{Vector{String}},Vector{Int64},UnitRange{Int64},Vector{Int64},Int64,Symbol,Dict{Int64, Symbol},Dict{Int64, Vector{Regex}},Dict{Tuple{Int64, Int64}, Symbol},Vector{Int64},Vector{Int64},Vector{Int64}})   # time: 0.013031783
    Base.precompile(Tuple{typeof(_print_table_data),IOContext{IOBuffer},Display,Any,Matrix{Vector{String}},UnitRange{Int64},UnitRange{Int64},Int64,Vector{Int64},Vector{Int64},Vector{NTuple{4, Int64}},Vector{Int64},Vector{Symbol},NTuple{4, Char},Dict{Tuple{Int64, Int64}, Symbol},Symbol,Int64,Ref{Any},Bool,Bool,Bool,TextFormat,Crayon,Crayon,Crayon})   # time: 0.012541843
    Base.precompile(Tuple{typeof(_fill_matrix_data!),Matrix{String},Matrix{Vector{String}},Vector{Int64},Vector{Int64},Vector{Int64},Vector{Int64},Vector{Int64},Int64,Any,Any,Ref{Any},Display,Bool,Bool,Vector{Int64},Bool,Bool,Bool,Bool,Vector{Int64},Vector{Int64},Bool,Val{:print},Symbol})   # time: 0.01242532
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Int64})   # time: 0.012378445
    Base.precompile(Tuple{typeof(_apply_alignment_anchor_regex!),Matrix{Vector{String}},Vector{Int64},Vector{Int64},Vector{Int64},Int64,Symbol,Dict{Int64, Symbol},Dict{Int64, Vector{Regex}},Dict{Tuple{Int64, Int64}, Symbol},Vector{Int64},Vector{Int64},Vector{Int64}})   # time: 0.01235328
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Int16})   # time: 0.011694558
    Base.precompile(Tuple{typeof(_fill_matrix_data!),Matrix{String},Matrix{Vector{String}},Vector{Int64},Vector{Int64},UnitRange{Int64},Vector{Int64},Vector{Int64},Int64,Any,Any,Ref{Any},Display,Bool,Bool,Vector{Int64},Bool,Bool,Bool,Bool,Vector{Int64},Vector{Int64},Bool,Val{:print},Symbol})   # time: 0.011301746
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),UInt64})   # time: 0.010704046
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Float32})   # time: 0.010583799
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Int32})   # time: 0.010526446
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Char})   # time: 0.010324799
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Float16})   # time: 0.009823338
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Int8})   # time: 0.009639627
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Float64})   # time: 0.009383389
    Base.precompile(Tuple{typeof(_print_table_data),IOContext{IOBuffer},Display,Any,Matrix{Vector{String}},Vector{Int64},UnitRange{Int64},Int64,Vector{Int64},Vector{Int64},Vector{NTuple{4, Int64}},Vector{Int64},Vector{Symbol},NTuple{4, Char},Dict{Tuple{Int64, Int64}, Symbol},Symbol,Int64,Ref{Any},Bool,Bool,Bool,TextFormat,Crayon,Crayon,Crayon})   # time: 0.009361974
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Bool})   # time: 0.008951217
    Base.precompile(Tuple{typeof(_print_table_header!),IOContext{IOBuffer},Display,Any,Matrix{String},UnitRange{Int64},UnitRange{Int64},Int64,Int64,Vector{Int64},Vector{Int64},Vector{Symbol},Vector{Symbol},Ref{Any},Bool,Bool,TextFormat,Crayon,Vector{Crayon},Vector{Crayon},Crayon,Crayon})   # time: 0.008807218
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),UInt16})   # time: 0.008730905
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),String})   # time: 0.008627564
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),UInt8})   # time: 0.00835835
    Base.precompile(Tuple{typeof(_print_table_data),IOContext{IOBuffer},Display,Any,Matrix{Vector{String}},Vector{Int64},Vector{Int64},Int64,Vector{Int64},Vector{Int64},Vector{NTuple{4, Int64}},Vector{Int64},Vector{Symbol},NTuple{4, Char},Dict{Tuple{Int64, Int64}, Symbol},Symbol,Int64,Ref{Any},Bool,Bool,Bool,TextFormat,Crayon,Crayon,Crayon})   # time: 0.007898714
    Base.precompile(Tuple{typeof(_print_table_header!),IOContext{IOBuffer},Display,Any,Matrix{String},UnitRange{Int64},Vector{Int64},Int64,Int64,Vector{Int64},Vector{Int64},Vector{Symbol},Vector{Symbol},Ref{Any},Bool,Bool,TextFormat,Crayon,Vector{Crayon},Vector{Crayon},Crayon,Crayon})   # time: 0.007784846
    Base.precompile(Tuple{typeof(_print_table_header!),IOContext{IOBuffer},Display,Any,Matrix{String},Vector{Int64},UnitRange{Int64},Int64,Int64,Vector{Int64},Vector{Int64},Vector{Symbol},Vector{Symbol},Ref{Any},Bool,Bool,TextFormat,Crayon,Vector{Crayon},Vector{Crayon},Crayon,Crayon})   # time: 0.007097781
    Base.precompile(Tuple{typeof(_print_table_header!),IOContext{IOBuffer},Display,Any,Matrix{String},Vector{Int64},Vector{Int64},Int64,Int64,Vector{Int64},Vector{Int64},Vector{Symbol},Vector{Symbol},Ref{Any},Bool,Bool,TextFormat,Crayon,Vector{Crayon},Vector{Crayon},Crayon,Crayon})   # time: 0.006893446
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_first_line_only, :column_width, :compact_printing, :linebreaks, :renderer), Tuple{Bool, Bool, Int64, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),String})   # time: 0.006176668
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),URLTextCell})   # time: 0.004810641
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),Symbol})   # time: 0.004649655
    Base.precompile(Tuple{Core.kwftype(typeof(_parse_cell_text)),NamedTuple{(:autowrap, :cell_data_type, :cell_first_line_only, :column_width, :compact_printing, :has_color, :limit_printing, :linebreaks, :renderer), Tuple{Bool, DataType, Bool, Int64, Bool, Bool, Bool, Bool, Val{:print}}},typeof(_parse_cell_text),String})   # time: 0.004569251
    Base.precompile(Tuple{typeof(_compute_cell_alignment_override),Vector{AnsiTextCell},UnitRange{Int64},UnitRange{Int64},Int64,Int64,Int64,Base.RefValue{Any}})   # time: 0.004113377
    Base.precompile(Tuple{typeof(_compute_row_fill_vectors),Vector{Int64},Int64,Symbol})   # time: 0.00348641
    Base.precompile(Tuple{typeof(_compute_cell_alignment_override),Matrix{Any},UnitRange{Int64},Vector{Int64},Int64,Int64,Int64,Base.RefValue{Any}})   # time: 0.003172167
    Base.precompile(Tuple{typeof(_compute_cell_alignment_override),ColumnTable,UnitRange{Int64},UnitRange{Int64},Int64,Int64,Int64,Base.RefValue{Any}})   # time: 0.003055841
    Base.precompile(Tuple{typeof(_compute_cell_alignment_override),Matrix{Any},UnitRange{Int64},UnitRange{Int64},Int64,Int64,Int64,Base.RefValue{Any}})   # time: 0.002971007
    Base.precompile(Tuple{typeof(_compute_cell_alignment_override),Vector{Any},UnitRange{Int64},UnitRange{Int64},Int64,Int64,Int64,Base.RefValue{Any}})   # time: 0.00282743
    Base.precompile(Tuple{typeof(_compute_cell_alignment_override),Matrix{Any},Vector{Int64},UnitRange{Int64},Int64,Int64,Int64,Base.RefValue{Any}})   # time: 0.00279592
    Base.precompile(Tuple{typeof(_compute_cell_alignment_override),Matrix{Any},Vector{Int64},Vector{Int64},Int64,Int64,Int64,Base.RefValue{Any}})   # time: 0.002695775
end


text:
using LinearAlgebra
using FileIO
using Makie

const TOLERANCE = 1e-7
N = 51
ϵ = 1e-7

"""
get_rotation_axis(i, n)

Calculates a unit quaternion corresponding to the axis of a 3D rotation with
the given initial and end unit vectors.
"""
function get_rotation_axis(i, n)
    u = normalize(cross(i, n))
    # The angle of rotation
    θ = acos(dot(i, n))
    u, θ/2
end

"""
build_surface(scene, points, color; transparency, shading)

Builds a surface with the given scene, points, color, transparency
and shading.
"""
function build_surface(
    scene,
    points,
    color;
    transparency = true,
    shading = false,
)
    surface!(
        scene,
        @lift($points[:, :, 1]),
        @lift($points[:, :, 2]),
        @lift($points[:, :, 3]),
        color = color,
        transparency = transparency,
        shading = shading,
    )
end

f(x, y) = (-x .* exp.(-x .^ 2 .- (y') .^ 2)) .* 4
f(A) = begin
    x, y = eachcol(A)
    (-x .* exp.(-x .^ 2 .- y .^ 2)) .* 4
end
d(Φ, A) = begin
    rows, cols = size(A)
    D = Matrix{Float64}(undef, rows, cols)
    E = Matrix{Float64}(I, cols, cols) .* ϵ
    for (i, v) in enumerate(eachrow(A))
        P = repeat(v, 1, cols)
        P′ = P + E
        D[i, :] .= (f(transpose(P′)) - f(transpose(P))) ./ ϵ
    end
    D
end
ξΦ(Φ, ξ, A) = begin
    derivative = d(Φ, A)
    D = similar(derivative)
    for (i, v) in enumerate(eachrow(derivative))
        D[i, :] .= v .* ξ
    end
    D
end

lspace = range(-2, stop = 2, length = N)
points = f(lspace, lspace)

universe = surface(lspace,
                   lspace,
                   f(lspace, lspace),
                   colormap = :rainbow,
                   show_axis = false,
                   resolution = (360, 360),
                   backgroundcolor = :royalblue)
xm, ym, zm = minimum(scene_limits(universe))

slu, olu = textslider(-2.0:0.01:2.0, "u", start = -0.5)
slv, olv = textslider(-2.0:0.01:2.0, "v", start = -0.09)
slξ, olξ = textslider(0.0:0.01:2pi, "ξ", start = 0)

tangent_tail(x) =
    [Point3f0(x[i, 1], x[i, 2], f(x[i, 1], x[i, 2])) for i = 1:size(x, 1)]
tangent_head(x, f, ξ) = [
    Point3f0(
        ξ[1] * ϵ,
        ξ[2] * ϵ,
        f(x[i, 1] + ξ[1] * ϵ, x[i, 2] + ξ[2] * ϵ) - f(x[i, 1], x[i, 2]),
    ) ./ ϵ for i = 1:size(x, 1)
]
stiff_tail(x) = [Point3f0(x[i, 1], x[i, 2], zm) for i = 1:size(x, 1)]
stiff_head(x, f, ξ) = begin
    derivative = ξΦ(f, ξ, x)
    [
        Point3f0(derivative[i, 1], derivative[i, 2], 0)
        for i = 1:size(x, 1)
    ]
end
normal(x, f) = begin
    tx = tangent_head(x, f, [1; 0])
    ty = tangent_head(x, f, [0; 1])
    [normalize(cross(tx[i], ty[i])) for i = 1:size(x, 1)]
end
plane(x, f, ξ) = begin
    t = tangent_tail(x)
    tx = tangent_head(x, f, [ξ[1]*0.1+1; 0])
    ty = tangent_head(x, f, [0; ξ[2]*0.1+1])
    surf(i) = begin
        p = Array{Float64}(undef, 2, 2, 3)
        p[1, 1, :] = normalize(tx[i] + ty[i] .+ 1e-5 * rand()) + t[i]
        p[1, 2, :] = normalize(-tx[i] + ty[i] .+ 1e-5 * rand()) + t[i]
        p[2, 1, :] = normalize(tx[i] - ty[i] .+ 1e-5 * rand()) + t[i]
        p[2, 2, :] = normalize(-tx[i] - ty[i] .+ 1e-5 * rand()) + t[i]
        p
    end
    [surf(i) for i = 1:size(x, 1)]
end
number = 1
P = @lift([$olu $olv])
tangent_arrowtail = @lift(tangent_tail($P))
tangent_arrowhead = @lift(tangent_head($P, f, [cos($olξ); sin($olξ)]))
yconst_tangent_arrowhead = @lift(tangent_head($P, f, [cos($olξ); 0]))
xconst_tangent_arrowhead = @lift(tangent_head($P, f, [0; sin($olξ)]))
normal_arrowhead = @lift(normal($P, f))
stiff_arrowtail = @lift(stiff_tail($P))
ξarrowhead = @lift(stiff_head($P, f, [cos($olξ); sin($olξ)]))
yconst_stiff_arrowhead = @lift(stiff_head($P, f, [cos($olξ); 0]))
xconst_stiff_arrowhead = @lift(stiff_head($P, f, [0; sin($olξ)]))

arrowmesh = load("test/data/boquarrow.obj")
arrow = mesh!(universe, arrowmesh, color = :gold)[end]
scale!(arrow, 0.02, 0.02, 0.02)
on(olu) do val
    i = [0; 1; 0]
    n = normalize(to_value(normal_arrowhead)[1])
    u, θ = get_rotation_axis(i, n)
    rotate!(arrow, Quaternionf0((sin(θ).*u)..., cos(θ)))
    x = to_value(olu)
    y = to_value(olv)
    z = f(x, y)
    translate!(arrow, Vec3f0(x, y, z))
end

on(olv) do val
    i = [0; 1; 0]
    n = normalize(to_value(normal_arrowhead)[1])
    u, θ = get_rotation_axis(i, n)
    #rotate!(arrow, Vec3f0(u...), θ)
    rotate!(arrow, Quaternionf0((sin(θ).*u)..., cos(θ)))
    x = to_value(olu)
    y = to_value(olv)
    z = f(x, y)
    translate!(arrow, Vec3f0(x, y, z))
end

# Instantiate a horizontal box for holding the visuals and the controls
#box = vbox(slu, slv, slξ)
#scene = hbox(universe, box, parent = Scene(resolution = (1920, 1080)))
contour!(
    universe,
    lspace,
    lspace,
    f(lspace, lspace),
    levels = 15,
    linewidth = 2,
    transformation = (:xy, zm),
    colormap = :cinferno,
)
wireframe!(
    universe,
    lspace,
    lspace,
    f(lspace, lspace),
    overdraw = true,
    transparency = true,
    color = (:black, 0.1),
)
image = load("gallery/porta.jpg")
for i in 1:number
    build_surface(universe,
                  @lift(plane($P, f, [cos($olξ); sin($olξ)])[i]),
                  image)
end
arrows!(
    universe,
    tangent_arrowtail,
    tangent_arrowhead,
    arrowsize = 0.1,
    linecolor = (:white, 2.0),
    linewidth = 3,
)
arrows!(
    universe,
    tangent_arrowtail,
    xconst_tangent_arrowhead,
    arrowsize = 0.1,
    linecolor = (:red, 2.0),
    linewidth = 3,
)
arrows!(
    universe,
    tangent_arrowtail,
    yconst_tangent_arrowhead,
    arrowsize = 0.1,
    linecolor = (:green, 2.0),
    linewidth = 3,
)

arrows!(
    universe,
    stiff_arrowtail,
    ξarrowhead,
    arrowsize = 0.1,
    linecolor = (:white, 2.0),
    linewidth = 3,
)
arrows!(
    universe,
    stiff_arrowtail,
    xconst_stiff_arrowhead,
    arrowsize = 0.1,
    linecolor = (:red, 2.0),
    linewidth = 3
)
arrows!(
    universe,
    stiff_arrowtail,
    yconst_stiff_arrowhead,
    arrowsize = 0.1,
    linecolor = (:green, 2.0),
    linewidth = 3,
)
center!(universe) # center the Scene on the display
# update eye position
eye_position, lookat = Vec3f0(3, 3, 3), Vec3f0(0)
update_cam!(universe, eye_position, lookat)
universe.center = false # prevent scene from recentering on display
#Makie.save("gallery/surfaces.jpg", scene)
display(universe)

record(universe, "gallery/surfaces.gif") do io
    frames = 90
    for i = 1:frames
        progress = i / frames
        olu[] = cos(progress * 2pi + pi/2)
        olv[] = sin(progress * 2pi + pi/2)
        olξ[] = progress * 2pi # animate scene
        eye_position = Vec3f0(1 - sin(progress * 2pi),
                              1 - sin(progress * 2pi),
                              1 + sin(progress * 2pi)) .* 3 .+ 0.00001
        update_cam!(universe, eye_position, to_value(tangent_arrowtail)[1])
        recordframe!(io) # record a new frame
    end
end


text:
module DAMMmodel

include("DAMM.jl")
export DAMM

end


text:
push!(LOAD_PATH,"../src/")
using Documenter, Literate, TrialDocs

#=
If Literate is included then the following must be added to the .travis.yml file
julia --project=docs/ -e 'using Pkg; Pkg.instantiate();
                    Pkg.add("Documenter");
                    Pkg.add("Literate");
                    Pkg.add("Plots");
                    Pkg.develop(PackageSpec(path=pwd()))'
specifically the lines Pkg.add("Documenter"); Pkg.add("Literate"); Pkg.add("Plots");
=#
# doesn't necessarily need to be generated on the fly
const examples_directory = pwd() * "/docs/example/"
const output_directory = pwd() * "/docs/src/generated/"
# read the examples in the directory
examples = readdir(examples_directory)

for example in examples
    example_filepath = examples_directory * example
    Literate.markdown(example_filepath, output_directory, documenter = true)
end

makedocs(
    modules = [TrialDocs],
    pages = [
        "Home" => "index.md",
        "Contribution Guide" => [
            "Iterative Solvers" => "iterative_solvers/contribution_guide.md",
            "Time Stepping" => "time_steppers/contribution_guide.md",
        ],
        "Examples" => [
            "first example" => "generated/literate_example.md"
        ],
        "Function Index" => "function_index.md"
    ],
    sitename = "TrialDocs"
)

deploydocs(repo = "github.com/sandreza/TrialDocs.git")

max_stars_repo_name: sandreza/TrialDocs

max_issues_repo_name: sandreza/TrialDocs

max_forks_repo_name: sandreza/TrialDocs


text:
@with_kw mutable struct AsynchronousData{T <: AbstractFloat}
    timestamp::Int = 1
    incumbent::Int = 1
    active::Int = 0
end

"""
    AsynchronousExecution

Functor object for using synchronous execution in an L-shaped algorithm (assuming multiple Julia cores are available). Create by supplying a [`Asynchronous`](@ref) object through `execution` in the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

"""
struct AsynchronousExecution{H <: AbstractFeasibilityHandler,
                             T <: AbstractFloat,
                             A <: AbstractVector} <: AbstractLShapedExecution
    data::AsynchronousData{T}
    subobjectives::Vector{A}
    model_objectives::Vector{A}
    finished::Vector{Int}
    subworkers::Vector{SubWorker{H,T}}
    decisions::Vector{DecisionChannel}
    work::Vector{Work}
    finalize::Vector{Work}
    metadata::Vector{MetaData}
    iterates::RemoteIterates{A}
    cutqueue::CutQueue{T}
    active_workers::Vector{Future}
    triggered::Vector{Bool}
    added::Vector{Bool}
    max_active::Int
    κ::T

    function AsynchronousExecution(structure::VerticalStructure{2, 1, <:Tuple{DistributedScenarioProblems}},
                                   max_active::Int, κ::T,
                                   ::Type{F}, ::Type{T}, ::Type{A}) where {F <: AbstractFeasibility,
                                                                           T <: AbstractFloat,
                                                                           A <: AbstractVector}
        H = HandlerType(F)
        return new{H,T,A}(AsynchronousData{T}(),
                          Vector{A}(),
                          Vector{A}(),
                          Vector{Int}(),
                          Vector{SubWorker{H,T}}(undef, nworkers()),
                          scenarioproblems(structure).decisions,
                          Vector{Work}(undef, nworkers()),
                          Vector{Work}(undef, nworkers()),
                          Vector{MetaData}(undef, nworkers()),
                          RemoteChannel(() -> IterateChannel(Dict{Int,A}())),
                          RemoteChannel(() -> Channel{QCut{T}}(max_active * nworkers() * num_scenarios(structure))),
                          Vector{Future}(undef, nworkers()),
                          Vector{Bool}(),
                          Vector{Bool}(),
                          max_active,
                          κ)
    end
end

function initialize_subproblems!(execution::AsynchronousExecution,
                                 scenarioproblems::DistributedScenarioProblems)
    load_subproblems!(execution.subworkers, scenarioproblems, execution.decisions)
    return nothing
end

function finish_initilization!(lshaped::AbstractLShaped, execution::AsynchronousExecution)
    # Load initial decision
    put!(execution.iterates, 1, copy(lshaped.x))
    # Prepare memory
    push!(execution.subobjectives, fill(1e10, num_thetas(lshaped)))
    push!(execution.model_objectives, fill(-1e10, num_thetas(lshaped)))
    push!(execution.finished, 0)
    push!(execution.triggered, false)
    push!(execution.added, false)
    push!(lshaped.Q_history, Inf)
    push!(lshaped.θ_history, -Inf)
    log_regularization!(lshaped)
    # Prepare work channels
    for w in workers()
        execution.work[w-1] = RemoteChannel(() -> Channel{Int}(execution.max_active+1), w)
        execution.finalize[w-1] = RemoteChannel(() -> Channel{Int}(1), w)
        execution.metadata[w-1] = RemoteChannel(() -> MetaChannel(), w)
        put!(execution.work[w-1], 1)
        put!(execution.metadata[w-1], 1, :gap, Inf)
    end
    return nothing
end

function start_workers!(lshaped::AbstractLShaped, execution::AsynchronousExecution)
    for w in workers()
        worker_aggregator = remote_aggregator(lshaped.aggregation, scenarioproblems(lshaped.structure), w)
        execution.active_workers[w-1] = remotecall(work_on_subproblems!,
                                                   w,
                                                   execution.subworkers[w-1],
                                                   execution.decisions[w-1],
                                                   execution.work[w-1],
                                                   execution.finalize[w-1],
                                                   execution.cutqueue,
                                                   execution.iterates,
                                                   execution.metadata[w-1],
                                                   worker_aggregator)
    end
    return nothing
end

function close_workers!(::AbstractLShaped, execution::AsynchronousExecution)
    t = execution.data.timestamp
    map((w, aw)->!isready(aw) && put!(w, t), execution.finalize, execution.active_workers)
    map((w, aw)->!isready(aw) && put!(w, -1), execution.work, execution.active_workers)
    map(wait, execution.active_workers)
    return nothing
end

function resolve_subproblems!(::AbstractLShaped, ::AsynchronousExecution)
    return nothing
end

function current_decision(lshaped::AbstractLShaped, execution::AsynchronousExecution)
    t = timestamp(lshaped)
    return fetch(execution.iterates, t)
end

function timestamp(::AbstractLShaped, execution::AsynchronousExecution)
    return execution.data.timestamp
end

function incumbent_decision(::AbstractLShaped, t::Integer, regularization::AbstractRegularization, execution::AsynchronousExecution)
    if t > 1
        return fetch(execution.iterates, regularization.incumbents[t])
    else
        return map(regularization.ξ) do ξᵢ
            return ξᵢ.value
        end
    end
end

function incumbent_decision(lshaped::AbstractLShaped, t::Integer, ::NoRegularization, execution::AsynchronousExecution)
    lshaped.x
end

function incumbent_objective(::AbstractLShaped, t::Integer, regularization::AbstractRegularization, ::AsynchronousExecution)
    return t > 1 ? regularization.Q̃_history[regularization.incumbents[t]] : regularization.data.Q̃
end

function incumbent_trustregion(::AbstractLShaped, t::Integer, rd::RegularizedDecomposition, ::AsynchronousExecution)
    return rd.σ_history[t]
end

function incumbent_trustregion(::AbstractLShaped, t::Integer, tr::TrustRegion, ::AsynchronousExecution)
    return tr.Δ_history[t]
end

function readd_cuts!(lshaped::AbstractLShaped, consolidation::Consolidation, execution::AsynchronousExecution)
    for i in eachindex(consolidation.cuts)
        for cut in consolidation.cuts[i]
            add_cut!(lshaped, cut, execution.model_objectives[i], execution.subobjectives[i], fetch(execution.iterates, i), check = false)
        end
        for cut in consolidation.feasibility_cuts[i]
            add_cut!(lshaped, cut, execution.model_objectives[i], execution.subobjectives[i], Inf)
        end
    end
    return nothing
end

function subobjectives(lshaped::AbstractLShaped, execution::AsynchronousExecution)
    t = timestamp(lshaped)
    return execution.subobjectives[t]
end

function set_subobjectives(lshaped::AbstractLShaped, Qs::AbstractVector, execution::AsynchronousExecution)
    t = timestamp(lshaped)
    execution.subobjectives[t] .= Qs
    return nothing
end

function model_objectives(lshaped::AbstractLShaped, execution::AsynchronousExecution)
    t = timestamp(lshaped)
    θs = t > 1 ? execution.model_objectives[t-1] : fill(-1e10, num_thetas(lshaped))
    return θs
end

function set_model_objectives(lshaped::AbstractLShaped, θs::AbstractVector, execution::AsynchronousExecution)
    t = timestamp(lshaped)
    ids = active_model_objectives(lshaped)
    execution.model_objectives[t][ids] .= θs[ids]
    return nothing
end

function restore_subproblems!(::AbstractLShaped, execution::AsynchronousExecution)
    restore_subproblems!(execution.subworkers)
    return nothing
end

function solve_master!(lshaped::AbstractLShaped, execution::AsynchronousExecution)
    try
        MOI.optimize!(lshaped.master)
    catch
        status = MOI.get(lshaped.master, MOI.TerminationStatus())
        # Master problem could not be solved for some reason.
        @unpack Q,θ = lshaped.data
        gap = abs(θ-Q)/(abs(Q)+1e-10)
        # Always print this warning
        @warn "Master problem could not be solved, solver returned status $status. The following relative tolerance was reached: $(@sprintf("%.1e",gap)). Aborting procedure."
        rethrow(err)
    end
    status = MOI.get(lshaped.master, MOI.TerminationStatus())
    if status == MOI.INFEASIBLE
        # Asynchronicity can sometimes yield an infeasible master. If so, try to continue
        if execution.max_active > 1 || execution.κ < 1.0
            # Ensure that there are still cuts to consider
            if isready(execution.cutqueue)
                # Return false OPTIMAL to continue procedure
                return MOI.OPTIMAL
            end
        end
    end
    # Otherwise, return status as usual
    return status
end

function iterate!(lshaped::AbstractLShaped, execution::AsynchronousExecution{H,T}) where {H <: AbstractFeasibilityHandler, T <: AbstractFloat}
    wait(execution.cutqueue)
    while isready(execution.cutqueue)
        new_iterate = false
        # Add new cuts from subworkers
        t::Int, cut::SparseHyperPlane{T} = take!(execution.cutqueue)
        execution.data.timestamp = t
        # Break if any subproblem is infeasible or unbounded
        if infeasible(cut)
            @warn "Stochastic program is not second-stage feasible at the current decision. Rerun procedure with feasibility_cuts = true to use feasibility cuts."
            # Early termination log
            log!(lshaped; status = MOI.INFEASIBLE)
            return MOI.INFEASIBLE
        end
        if !bounded(cut)
            # Early termination log
            log!(lshaped; status = MOI.DUAL_INFEASIBLE)
            return MOI.DUAL_INFEASIBLE
        end
        # Otherwise, add new cut to master and update bookkeeping
        execution.added[t] |= add_cut!(lshaped, cut)
        execution.finished[t] += num_subproblems(cut)
        # Asynchronicity parameter should be 1 first iteration
        κ = t > 1 ? execution.κ : 1.0
        if execution.finished[t] == num_thetas(lshaped)
            # All work from iteration t complete
            lshaped.Q_history[t] = current_objective_value(lshaped)
            lshaped.data.Q = lshaped.Q_history[t]
            lshaped.data.θ = t > 1 ? lshaped.θ_history[t-1] : -1e10
            # Update incumbent (if applicable)
            take_step!(lshaped)
            # Early optimality check if using level sets
            if lshaped.regularization isa LevelSet && check_optimality(lshaped)
                # Resolve subproblems with optimal vector
                lshaped.x .= decision(lshaped)
                t = lshaped.data.iterations
                put!(execution.iterates, t+1, copy(lshaped.x))
                for w in workers()
                    if !isready(execution.active_workers[w-1])
                        put!(execution.work[w-1], t+1)
                        put!(execution.metadata[w-1], t+1, :gap, gap(lshaped))
                    end
                end
                execution.data.timestamp = t + 1
                # Optimal, final log
                log!(lshaped; optimal = true)
                return MOI.OPTIMAL
            end
            # Optimal if no cuts were added and current decision
            # is not different from incumbent projection target
            if !execution.added[t] && norm(incumbent_decision(lshaped, t, lshaped.regularization, execution) - lshaped.x) <= sqrt(eps())
                # Optimal, final log
                log!(lshaped, t; optimal = true)
                return MOI.OPTIMAL
            end
            # Consolidate (if applicable)
            consolidate!(lshaped, lshaped.consolidation)
            # Decrease number of active iterations
            execution.data.active = max(0, execution.data.active - 1)
            # Determine if new iterate should be generated
            new_iterate |= !(lshaped.regularization isa NoRegularization && execution.triggered[t])
            new_iterate |= t == 1
        elseif execution.finished[t] >= κ * num_thetas(lshaped) &&
               !execution.triggered[t] &&
               execution.data.active < execution.max_active &&
               execution.added[t]
            execution.triggered[t] = true
            new_iterate = true
        else
            new_iterate = false
        end
        # Generate new candidate decision
        if new_iterate
            t = lshaped.data.iterations
            execution.data.timestamp = t
            # Resolve master
            status = solve_master!(lshaped)
            if !(status ∈ AcceptableTermination)
                # Early termination log
                log!(lshaped; status = status)
                return status
            end
            # Update master solution
            update_solution!(lshaped)
            θ = calculate_estimate(lshaped)
            lshaped.data.θ = θ
            lshaped.θ_history[t] = θ
            # Check if optimal
            if check_optimality(lshaped)
                # Optimal, final log
                log!(lshaped, t)
                return MOI.OPTIMAL
            end
            # Log progress at current timestamp
            log_regularization!(lshaped, t)
            # Update workers
            put!(execution.iterates, t+1, copy(lshaped.x))
            for w in workers()
                if !isready(execution.active_workers[w-1])
                    put!(execution.work[w-1], t+1)
                    put!(execution.metadata[w-1], t+1, :gap, gap(lshaped))
                end
            end
            # New active iteration
            execution.data.active += 1
            # Prepare memory for next iteration
            push!(execution.subobjectives, fill(1e10, num_thetas(lshaped)))
            push!(execution.model_objectives, fill(-1e10, num_thetas(lshaped)))
            push!(execution.finished, 0)
            push!(execution.triggered, false)
            push!(execution.added, false)
            if lshaped.consolidation isa Consolidation
                allocate!(lshaped.consolidation)
            end
            # Log progress
            log!(lshaped)
            lshaped.θ_history[t+1] = -Inf
        end
    end
    # Dont return a status as procedure should continue
    return nothing
end

# API
# ------------------------------------------------------------
function (execution::Asynchronous)(structure::VerticalStructure{2, 1, <:Tuple{DistributedScenarioProblems}},
                                   ::Type{F}, ::Type{T}, ::Type{A}) where {F <: AbstractFeasibility,
                                                                           T <: AbstractFloat,
                                                                           A <: AbstractVector}
    return AsynchronousExecution(structure, execution.max_active, execution.κ, F, T, A)
end

function str(::Asynchronous)
    return "Asynchronous "
end


text:
include("../Utils/PolymeshReaderUtils.jl")

function getpolymeshfiles(polymeshdir)
    polymeshfiles = Dict()
    polymeshfiles["points"] = joinpath(polymeshdir, "points")
    polymeshfiles["faces"] = joinpath(polymeshdir, "faces")
    polymeshfiles["boundary"] = joinpath(polymeshdir, "boundary")
    polymeshfiles["owner"] = joinpath(polymeshdir, "owner")
    polymeshfiles["neighbour"] = joinpath(polymeshdir, "neighbour")
    return polymeshfiles
end

function readpoints(pointsfile)
    points = []
    lines, npoints, linestart = getfileitems(pointsfile)
    for i in 1:npoints
        push!(points, stringtofloatlist(lines[i + linestart - 1]))
    end
    return points
end

function readfaces(facesfile)
    faces = []
    lines, nfaces, linestart = getfileitems(facesfile)
    for i in 1:nfaces
        push!(faces, stringtointlist(lines[i + linestart - 1][2:end]))
        faces[i] = faces[i] .+ 1
    end
    return faces
end

function readboundary(boundaryfile)
    boundary = []
    lines, nboundaries, linestart = getfileitems(boundaryfile)
    j = linestart
    i = 0
    while true
        if lines[j] == ")"
            break
        else
            if lines[j + 1] == "{" # face detected
                i += 1
                push!(boundary, Dict())
                boundary[i]["name"] = lines[j]
                k = j + 2
                while true
                    if lines[k] == "}" # end of a particular face
                        break
                    else
                        key, value = linetokeyvalue(lines[k])
                        if key == "nFaces"
                            value = parse(Int64, value)
                        elseif key == "startFace"
                            value = parse(Int64, value) + 1
                        end
                        boundary[i][key] = value
                    k += 1
                    end
                end
            end
        end
        j += 1
    end
    return boundary
end

function readowner(ownerfile)
    owner = []
    lines, nfaces, linestart = getfileitems(ownerfile)
    for i in 1:nfaces
        push!(owner, parse(Int64, lines[i + linestart - 1]))
        owner[i] = owner[i] .+ 1
    end
    return owner
end

function readneighbour(neighbourfile)
    neighbour = []
    lines, nifaces, linestart = getfileitems(neighbourfile)
    for i in 1:nifaces
        push!(neighbour, parse(Int64, lines[i + linestart - 1]))
        neighbour[i] = neighbour[i] .+ 1
    end
    return neighbour
end

function readpolymesh(polymeshfiles)
    polymesh = Dict()
    polymesh["points"] = readpoints(polymeshfiles["points"])
    polymesh["faces"] = readfaces(polymeshfiles["faces"])
    polymesh["boundary"] = readboundary(polymeshfiles["boundary"])
    polymesh["owner"] = readowner(polymeshfiles["owner"])
    polymesh["neighbour"] = readneighbour(polymeshfiles["neighbour"])
    return polymesh
end