text:
###############################################################################
#
# Update carbon and water fluxes
#
# TODO: work more to make this function more customized
#
###############################################################################
"""
    layer_fluxes!(
                node::SPACMono{FT};
                updating::Bool = false
    ) where {FT<:AbstractFloat}

Run carbon, water, energy, and SIF fluxes for all canopy layers, given
- `node` [`SPACMono`](@ref) type struct
- `updating` If true, update cavitation history
"""
function layer_fluxes!(
            node::SPACMono{FT};
            updating::Bool = false
) where {FT<:AbstractFloat}
    # 0.1 unpack data
    @unpack angles, envirs, f_SL, ga, in_rad, leaves_rt, n_canopy, plant_ps, rt_con, rt_dim, soil_opt, stomata_model, wl_set = node;
    canopy_rt = node.canopy_rt;
    can_rad = node.can_rad;
    can_opt = node.can_opt;
    plant_hs = node.plant_hs;
    nSL = canopy_rt.nAzi * canopy_rt.nIncl;

    canopy_geometry!(canopy_rt, angles, can_opt, rt_con);
    canopy_matrices!(leaves_rt, can_opt);
    short_wave!(canopy_rt, can_opt, can_rad, in_rad, soil_opt, rt_con);
    canopy_fluxes!(canopy_rt, can_opt, can_rad, in_rad, soil_opt, leaves_rt,
                   wl_set, rt_con);

    #
    f_H₂O = 0;
    f_GPP = 0;
    f_NPP = 0;
    for i_can in 1:n_canopy
        iEN = envirs[i_can];
        iHS = plant_hs.leaves[i_can];
        iPS = plant_ps[i_can];
        iRT = n_canopy + 1 - i_can;

        iPS.T = can_rad.T_sun[i_can];
        update_leaf_TP!(iPS, iHS, iEN);
        temperature_effects!(iHS, iPS.T);

        # calculate the fraction of sunlit and shaded leaves
        f_view = (can_opt.Ps[iRT] + can_opt.Ps[iRT+1]) / 2;
        for iLF in 1:nSL
            iPS.APAR[iLF] = can_rad.absPAR_sunCab[(iRT-1)*nSL+iLF] * FT(1e6);
            iPS.LAIx[iLF] = f_view * f_SL[iLF];
        end
        iPS.APAR[end] = can_rad.absPAR_shadeCab[iRT] * FT(1e6);
        iPS.LAIx[end] = 1 - f_view;

        # iterate for N times until sum(iPS.Ag) does not change any more
        sum_ag_last = sum(iPS.Ag);
        count       = 0;
        while true
            # calculate the photosynthetic rates
            count += 1;
            gas_exchange!(iPS, iEN, GswDrive());
            if typeof(stomata_model) <: EmpiricalStomatalModel
                #
                #
                #
                #
                # also need to determine which β function is used
                #     to tune g1
                #     to tune Vcmax
                #
                #
                #
                #
                prognostic_gsw!(iPS, iEN, stomata_model, FT(1), FT(120));
            else
                prognostic_gsw!(iPS, iHS, iEN, stomata_model, FT(120));
            end

            # update flow and pressure profile (except for history)
            # TODO move this part to PlantHydraulics.jl
            gsw_control!(iPS, iEN);
            for i_can in 1:n_canopy
                iEN = envirs[i_can];
                iLF = plant_hs.leaves[i_can];
                iPS = plant_ps[i_can];
                iLF.flow = sum(iPS.g_lw .* iPS.LAIx) * (iPS.p_sat - iEN.p_H₂O) / iEN.P_AIR;
            end
            flow_profile!(plant_hs);
            pressure_profile!(plant_hs, SteadyStateMode(); update=false);

            # TODO make sure this is done in future SPAC module
            # update canopy layer p_ups, which will be passed to each leaf
            for _i_can in 1:n_canopy
                _iHS = plant_hs.leaves[_i_can];
                _iPS = plant_ps[_i_can];
                _iPS.p_ups = _iHS.p_ups;
            end;

            sum_ag_curr = sum(iPS.Ag);
            if (abs(sum_ag_curr - sum_ag_last) < 1e-4)
                break
            elseif count > 100000
                @info tinfo("layer_fluxes!: iterations exceed 100000 times");
                break
            else
                sum_ag_last = sum_ag_curr;
            end
        end

        # update the fluorescence quantum yield from leaf level calculation
        can_rad.ϕ_sun[:,:,iRT] .= reshape(view(iPS.φs,1:nSL), canopy_rt.nIncl,
                                          canopy_rt.nAzi);
        can_rad.ϕ_shade[iRT] = iPS.φs[end];

        # update the flow rates
        for iLF in 1:(nSL+1)
            f_GPP += iPS.Ag[iLF] * iPS.LAIx[iLF] * iPS.LA;
            f_NPP += iPS.An[iLF] * iPS.LAIx[iLF] * iPS.LA;
            f_H₂O += iPS.g_lw[iLF] * (iPS.p_sat - iEN.p_H₂O) / iEN.P_AIR * iPS.LAIx[iLF] * iPS.LA;
        end
    end

    # do SIF simulation
    SIF_fluxes!(leaves_rt, can_opt, can_rad, canopy_rt, soil_opt, wl_set,
                rt_con, rt_dim);

    # update flow profile and pressure history along the tree
    if updating
        for i_can in 1:n_canopy
            iEN = envirs[i_can];
            iLF = plant_hs.leaves[i_can];
            iPS = plant_ps[i_can];
            iLF.flow = sum(iPS.g_lw .* iPS.LAIx) * (iPS.p_sat - iEN.p_H₂O) / iEN.P_AIR;
        end
        flow_profile!(plant_hs);
        pressure_profile!(plant_hs, SteadyStateMode(); update=true);
    end

    # update the flows in SPACMono
    node.f_gpp = f_GPP / ga;
    node.f_npp = f_NPP / ga;
    node.f_H₂O = f_H₂O / ga;

    return nothing
end

"""
    layer_fluxes!(
                node::SPACMono{FT},
                Δt::FT;
                updating::Bool = false
    ) where {FT<:AbstractFloat}

Run carbon, water, energy, and SIF fluxes for all canopy layers, given
- `node` [`SPACMono`](@ref) type struct
- `Δt` Time step to forward in time
- `updating` If true, update cavitation history
"""
function layer_fluxes!(
            node::SPACMono{FT},
            Δt::FT;
            updating::Bool = false
) where {FT<:AbstractFloat}
    # 0.1 unpack data
    @unpack angles, envirs, f_SL, ga, in_rad, leaves_rt, n_canopy, plant_ps, rt_con, rt_dim, soil_opt, stomata_model, wl_set = node;
    canopy_rt = node.canopy_rt;
    can_rad = node.can_rad;
    can_opt = node.can_opt;
    plant_hs = node.plant_hs;
    nSL = canopy_rt.nAzi * canopy_rt.nIncl;

    canopy_geometry!(canopy_rt, angles, can_opt, rt_con);
    canopy_matrices!(leaves_rt, can_opt);
    short_wave!(canopy_rt, can_opt, can_rad, in_rad, soil_opt, rt_con);
    canopy_fluxes!(canopy_rt, can_opt, can_rad, in_rad, soil_opt, leaves_rt,
                   wl_set, rt_con);

    #
    f_H₂O = 0;
    f_GPP = 0;
    f_NPP = 0;
    for i_can in 1:n_canopy
        iEN = envirs[i_can];
        iHS = plant_hs.leaves[i_can];
        iPS = plant_ps[i_can];
        iRT = n_canopy + 1 - i_can;

        iPS.T = can_rad.T_sun[i_can];
        update_leaf_TP!(iPS, iHS, iEN);
        temperature_effects!(iHS, iPS.T);

        # calculate the fraction of sunlit and shaded leaves
        f_view = (can_opt.Ps[iRT] + can_opt.Ps[iRT+1]) / 2;
        for iLF in 1:nSL
            iPS.APAR[iLF] = can_rad.absPAR_sunCab[(iRT-1)*nSL+iLF] * FT(1e6);
            iPS.LAIx[iLF] = f_view * f_SL[iLF];
        end
        iPS.APAR[end] = can_rad.absPAR_shadeCab[iRT] * FT(1e6);
        iPS.LAIx[end] = 1 - f_view;

        # step forward with Δt
        gas_exchange!(iPS, iEN, GswDrive());
        if typeof(stomata_model) <: EmpiricalStomatalModel
            #
            #
            #
            #
            # also need to determine which β function is used
            #     to tune g1
            #     to tune Vcmax
            #
            #
            #
            #
            prognostic_gsw!(iPS, iEN, stomata_model, FT(1), Δt);
        else
            prognostic_gsw!(iPS, iHS, iEN, stomata_model, Δt);
        end

        # update flow and pressure profile (except for history)
        # TODO move this part to PlantHydraulics.jl
        gsw_control!(iPS, iEN);
        for i_can in 1:n_canopy
            iEN = envirs[i_can];
            iLF = plant_hs.leaves[i_can];
            iPS = plant_ps[i_can];
            iLF.flow = sum(iPS.g_lw .* iPS.LAIx) * (iPS.p_sat - iEN.p_H₂O) / iEN.P_AIR;
        end
        flow_profile!(plant_hs);
        pressure_profile!(plant_hs, SteadyStateMode(); update=false);

        # TODO make sure this is done in future SPAC module
        # update canopy layer p_ups, which will be passed to each leaf
        for _i_can in 1:n_canopy
            _iHS = plant_hs.leaves[_i_can];
            _iPS = plant_ps[_i_can];
            _iPS.p_ups = _iHS.p_ups;
        end;

        # update the fluorescence quantum yield from leaf level calculation
        can_rad.ϕ_sun[:,:,iRT] .= reshape(view(iPS.φs,1:nSL), canopy_rt.nIncl, canopy_rt.nAzi);
        can_rad.ϕ_shade[iRT] = iPS.φs[end];

        # update the flow rates
        for iLF in 1:(nSL+1)
            f_GPP += iPS.Ag[iLF] * iPS.LAIx[iLF] * iPS.LA;
            f_NPP += iPS.An[iLF] * iPS.LAIx[iLF] * iPS.LA;
            f_H₂O += iPS.g_lw[iLF] * (iPS.p_sat - iEN.p_H₂O) / iEN.P_AIR * iPS.LAIx[iLF] * iPS.LA;
        end
    end

    # do SIF simulation
    SIF_fluxes!(leaves_rt, can_opt, can_rad, canopy_rt, soil_opt, wl_set, rt_con, rt_dim);

    # update flow profile and pressure history along the tree
    if updating
        for i_can in 1:n_canopy
            iEN = envirs[i_can];
            iLF = plant_hs.leaves[i_can];
            iPS = plant_ps[i_can];
            iLF.flow = sum(iPS.g_lw .* iPS.LAIx) * (iPS.p_sat - iEN.p_H₂O) / iEN.P_AIR;
        end
        flow_profile!(plant_hs);
        pressure_profile!(plant_hs, SteadyStateMode(); update=true);
    end

    # update the flows in SPACMono
    node.f_gpp = f_GPP / ga;
    node.f_npp = f_NPP / ga;
    node.f_H₂O = f_H₂O / ga;

    return nothing
end


text:
using Documenter, AzSessions

makedocs(sitename="AzSessions", modules=[AzSessions])

deploydocs(
    repo = "github.com/ChevronETC/AzSessions.jl.git",
)

ext: jl


text:
@testset "Carriers" begin
    @field(F1) do
        I₀ = 1.0
        T = 2.0
        τ = 3.0
        σmax = 6.0
        ϕ = π/4
    end

    car1 = carrier(F1)
    @test car1 isa ElectricFields.FixedCarrier
    @test car1(0.4) ≈ sin(π*(0.4+1/4))

    @test string(car1) == "Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV); CEP = 0.25π"

    @test wavelength(car1) |> u"nm" == 2*austrip(1u"c")*u"bohr" |> u"nm"
    @test austrip(period(car1)) == 2.0
    @test austrip(frequency(car1)) == 0.5
    @test wavenumber(car1) |> u"cm^-1" == 1/(2*austrip(1u"c")*u"bohr") |> u"cm^-1"
    @test fundamental(car1) ≈ π
    @test photon_energy(car1) ≈ π
    @test phase(car1) == π/4

    @field(F2) do
        I₀ = 1.0
        T = 2.0
        τ = 3.0
        σmax = 6.0
        ϕ = π/4
        kind = :transverse
    end

    car2 = carrier(F2)
    @test car2 isa ElectricFields.LinearTransverseCarrier

    @test string(car2) == "LinearTransverseCarrier: Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV); CEP = 0.25π"

    @test car2(0.4) ≈ [0,0,sin(π*(0.4+1/4))]
    @test phase_shift(car2, -π/4)(0.4) ≈ [0,0,sin(π*0.4)]

    @test wavelength(car2) |> u"nm" == 2*austrip(1u"c")*u"bohr" |> u"nm"
    @test austrip(period(car2)) == 2.0
    @test austrip(frequency(car2)) == 0.5
    @test wavenumber(car2) |> u"cm^-1" == 1/(2*austrip(1u"c")*u"bohr") |> u"cm^-1"
    @test fundamental(car2) ≈ π
    @test photon_energy(car2) ≈ π
    @test phase(car2) == π/4

    for (ξ,label) in [(0, "0.00"),
                      (1, "1.00 (RCP)"),
                      (-1, "-1.00 (LCP)"),
                      (1/2, "0.50 (right)"),
                      (-1/2, "-0.50 (left)")]
        @field(F3) do
            I₀ = 1.0
            T = 2.0
            τ = 3.0
            σmax = 6.0
            ϕ = π/4
            ξ = ξ
        end
        car3 = carrier(F3)

        @test string(car3) == "Elliptical carrier with ξ = $(label) @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV); CEP = 0.25π"

        @test wavelength(car3) |> u"nm" == 2*austrip(1u"c")*u"bohr" |> u"nm"
        @test austrip(period(car3)) == 2.0
        @test austrip(frequency(car3)) == 0.5
        @test wavenumber(car3) |> u"cm^-1" == 1/(2*austrip(1u"c")*u"bohr") |> u"cm^-1"
        @test fundamental(car3) ≈ π
        @test photon_energy(car3) ≈ π
        @test phase(car3) == π/4
    end
end


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

using Documenter
using PkgTestSuite

makedocs(modules=[PkgTestSuite],
         format = :html,
         checkdocs = :all,
         sitename = "PkgTestSuite.jl",
         pages = ["index.md"]
        )


text:
using QuadraticAssignmentProblem

using ProgressMeter
using DataFrames
using CSV
using Dates
using Random

using Gurobi
qap_set_optimizer!(Gurobi.Optimizer)

iterations = 10

problem = (
	N = 10,
	objective = :max,
	matrix_type = :known,
	density = missing,
)

algorithms = [
	# (name = "opt",    method = qap_exact),
	# (name = "lin",    method = qap_linearization),
	(name = "lpr",    method = qap_lprounding),
	(name = "faq",    method = faq),
	(name = "rand",   method = qap_random),
	(name = "mms",    method = qap_mms),
	(name = "ns",     method = qap_ns),
	# (name = "star",   method = qap_starpacking),
	# (name = "dense",  method = qap_densemapping),
]

function compare(problem, algorithms, iterations)
	warmup_problem = merge(problem, (;N=5))
	compare_step(warmup_problem, algorithms)

	results = []
	@showprogress for i in 1:iterations
		append!(results, compare_step(problem, algorithms))
	end

	save_results(results)

	return results
end

function compare_step(problem, algorithms)
	A, B, opt = generate_example(problem)
	uid = rand(UInt)

	results = []
	for a in algorithms
		alg = a.method
		time = @elapsed P, _ = alg(A, B, problem.objective)
		obj = qap_objective(A, B, P)

		if (a.name == "opt") opt = obj end

		r = merge(problem, (;uid, method=a.name, time, obj, ratio=obj/opt))
		push!(results, r)
	end

	return results
end

function generate_example(problem)
	N = problem.N
	opt = missing
	if problem.matrix_type == :uniform
		A = uniform_matrix(N)
		B = uniform_matrix(N)
	elseif problem.matrix_type == :metric
		A = metric_matrix(N)
		B = metric_matrix(N)
	elseif problem.matrix_type == :undirected_graph
		A = zeroone_matrix(N, p=problem.density, selfedges=false, symmetric=true)
		B = zeroone_matrix(N, p=problem.density, selfedges=false, symmetric=true)
	elseif problem.matrix_type == :directed_graph
		A = zeroone_matrix(N, p=problem.density, selfedges=false, symmetric=false)
		B = zeroone_matrix(N, p=problem.density, selfedges=false, symmetric=false)
	elseif problem.matrix_type == :known
		A, B, perm = generate_qap(N)
		opt = qap_objective(A, B, perm)
	else
		error("Invalid matrix type: $(problem.matrix_type)")
	end

	return A, B, opt
end

function save_results(results)

	d = Dates.format(now(), "yyyymmddHHMM")
	d = isdir("results/$d/") ? d * "_" * randstring(8) : d

	dir = joinpath("./results/", d)
	mkpath(dir)

	metric_cols = [:time, :obj, :ratio]
	problem_cols = setdiff(keys(results[1]), metric_cols, [:uid])

	results_df_long = DataFrame(results)
	results_df_long |> CSV.write(joinpath(dir, "results_long.csv"))

	mean(xs) = sum(xs) / length(xs)
	results_mean_df = combine(groupby(results_df_long, problem_cols), [m => mean => m for m in metric_cols]...)
	results_mean_df |> CSV.write(joinpath(dir, "results_mean.csv"))

	results_df_wide = stack(results_df_long, metric_cols)
	transform!(results_df_wide, [:variable, :method] => ByRow((metric, method) -> metric*"_"*method) => :key)
	select!(results_df_wide, Not([:variable, :method]))
	results_df_wide = unstack(results_df_wide, :key, :value)
	results_df_wide |> CSV.write(joinpath(dir, "results_wide.csv"))

	return
end

if abspath(PROGRAM_FILE) == @__FILE__
	compare(problem, algorithms, iterations)
end


text:
function fib(max_val)
	res, val, i_m2 = 0, 0, 0
	i_m1 = 1

	while val < max_val
		val = i_m1 + i_m2
		i_m1, i_m2 = val, i_m1
		((val & 1) == 0) && (res += val)
	end
	println("The imperative answer to Euler Problem 2 is $res")
end

@time fib(4000000)

function fibarr(max_val, )
	arr = [1,1]
	# println(arr)
	while arr[end] < max_val
		arr = [arr, arr[end]+arr[end-1]]
	end
	arr
end

@time (farr = fibarr(4000000);
       res2 = +(filter(x -> x & 1 == 0, farr)...))

println("The slightly more idiomatic answer to Problem 3 is $res2")

max_stars_repo_path: Problem 2.jl

max_stars_repo_name: andrew-christianson/Polyglot-Euler

max_issues_repo_path: Problem 2.jl

max_issues_repo_name: andrew-christianson/Polyglot-Euler

max_forks_repo_path: Problem 2.jl

max_forks_repo_name: andrew-christianson/Polyglot-Euler


text:
#=

This file is part of Pseudospectra.jl.

Copyright (c) 2017 Ralph A. Smith

SPDX-License-Identifier: BSD-3-Clause
License-Filename: LICENSES/BSD-3-Clause_Eigtool
=#

"""
object to hold state for the GUI used by Pseudospectra
"""
abstract type GUIState end

mutable struct LevelDesc
    isunif::Bool
    full_levels
    first
    step
    last
    function LevelDesc(levs)
        if isa(levs,AbstractRange)
            levels = collect(levs)
        elseif isa(levs,Vector) && (eltype(levs) <: Real)
            levels = levs
        else
            throw(ArgumentError("level specification must be a range or real vector"))
        end
        isunif = ((length(levels) == 2)
                  || ((length(levels) > 2)
                      && (maximum(abs.(diff(diff(levels)))) < 10*eps())))
        if isunif
            first = levels[1]
            step = levels[2]-levels[1]
            last = levels[end]
            full_levels = nothing
        else
            full_levels = copy(levels)
            first = nothing; step = nothing; last = nothing
        end
        new(isunif,full_levels,first,step,last)
    end
end

"""
    Portrait

structure representing a spectral portrait; includes a mesh for `z=x+iy`,
values of the resolvent norm on the mesh, suitable contour levels,
and some other metadata.
"""
mutable struct Portrait
    x::Vector
    y::Vector
    Z::Matrix
    npts::Int
    ax::Vector
    levels::LevelDesc
    autolev::Bool
    proj_lev
    dims
    computed::Bool
    scale_equal::Bool
end

import Base: show

function show(io::IO,z::Portrait)
    k = isempty(z.Z) ? "missing" : "present"
    print(io,"npts: $(z.npts), ax: $(z.ax), autolev: $(z.autolev), "
          * "computed: $(z.computed), data $k, proj $(z.proj_lev), "
          * "scale_eq: $(z.scale_equal)\n",
          "levels: $(z.levels)")
end

"""
    ArpackOptions{T}(; nev=6, ncv=0, which=:LM,
                                           tol=zero(T),maxiter=300,
                                           v0=Vector{T}(0), have_v0=false,
                                           sigma=nothing)

constructs an object to manage the Arnoldi scheme; see the documentation
for `eigs` in the `Arpack` package for the meaning of fields.
"""
mutable struct ArpackOptions{T}
    # Control of iterative computations:
    nev::Int # nr. of eigenvalues for eigs() to search for
    ncv::Int # max. subspace size for eigs()
    which::Symbol # as for eigs()
    tol::Real # tolerance for eigs()
    maxiter::Int # bound for eigs()
    v0::Vector # initial Ritz vector for eigs()
    sigma # shift parameter for eigs()

    function ArpackOptions{T}(; nev=6, ncv=0, which=:LM,
                              tol=zero(T),maxiter=300,
                              v0=Vector{T}(),
                              sigma=nothing) where T
        if ncv==0
            ncv = max(20,2*nev+1)
        end
        new{T}(nev,ncv,which,tol,maxiter,v0,sigma)
    end
end

function Base.:(==)(l::ArpackOptions{T},r::ArpackOptions{S}) where {T,S}
    l === r && return true
    for f in fieldnames(ArpackOptions)
        (getfield(l,f) == getfield(r,f)) || return false
    end
    true
end
function Base.hash(l::ArpackOptions{T},h::UInt) where T
    h1 = hash(:ArpackOptions,h)
    for f in fieldnames(ArpackOptions)
        h1 = hash(getfield(l,f),h1)
    end
    h1
end

"""
Wrapper structure for Pseudospectra session data
"""
mutable struct PSAStruct
    matrix
    unitary_mtx
    input_matrix
    input_unitary_mtx
    ps_dict::Dict{Symbol,Any}
    zoom_list::Vector{Portrait}
    zoom_pos::Int
    # aside: why is there a :proj_lev entry in ps_dict?
    function PSAStruct(m1,u1,m1i,u1i,dict)
        new(m1,u1,m1i,u1i,dict,Vector{Portrait}(),0)
    end
end
#=
 A (probably incomplete) list of keys in ps_dict
 :Aisreal
 :isHessenberg
 :schur_mtx
 :schur_unitary_mtx
 :direct
 :ews
 :orig_ews
 :ew_estimates
 :arpack_opts
 :init_opts
 :init_direct
 :verbosity
 :fov - samples of the numerical range (field of values)
 :projection_on
 :proj_lev
 :proj_valid
 :proj_ews
 :proj_matrix - the Hessenberg matrix of coeffts from IRAM (ARPACK)
 :proj_unitary_mtx - the Krylov basis from IRAM (ARPACK)
 :matrix2 - second matrix (if QZ is used for rectangular A)

 :reentering - flag for logic (not in EigTool)
 :proj_axes
 :comp_proj_lev
 :mode_markers

 Used only in GUI code:
 :zpsradius
 :zpsabscissa
=#


text:
@default function contour(scene, kw_args)
    levels = to_float(levels)
    color = to_color(color)
    linewidth = to_float(1)
    fillrange = to_bool(fillrange)
end

function contour(scene::makie, x, y, z, attributes)
    attributes = contour_defaults(scene, attributes)

    if to_value(attributes[:fillrange])
        return heatmap(scene, x, y, z, attributes)
    else
        levels = round(Int, to_value(attributes[:levels]))
        T = eltype(z)

        contours = Contour.contours(T.(x), T.(y), z, levels)
        result = Point2f0[]
        colors = RGBA{Float32}[]
        col = to_value(attributes[:color])
        cols = if isa(col, AbstractVector)
            if length(col) != levels
                error("Please have one color per level. Found: $(length(col)) colors and $levels level")
            end
            col
        else
            repeated(col, levels)
        end
        for (color, c) in zip(cols, Contour.levels(contours))
            for elem in Contour.lines(c)
                append!(result, elem.vertices)
                push!(result, Point2f0(NaN32))
                append!(colors, fill(color, length(elem.vertices) + 1))
            end
        end
        attributes[:color] = colors
        return lines(scene, result, attributes)
    end
end


text:
function Rmat(n, l)
    if l == 0
        e = ones(n)
        R = myspdiagm(spdiagm(-1 => 0.5*e, 0 => e, 1 => 0.5*e), n)
        R[2, 1] = 1
    else
        nn = collect(0:n-1)
        e = ones(n)
        e1 = 0.5*(nn.+1)./(nn.+l)
        e2 = 0.5*(nn.+2*l.-1)./(nn.+l)
        R = myspdiagm(spdiagm(-1 => e1, 0 => e, 1 => e2[2:end]), n)
    end
    return R
end

myspdiagm(M, N) = size(M, 1) !== N ? (return M[1:N, 1:N]) : (return M)


text:
abstract type AbstractManifold{N} end

function get_steer_vec(manifold, doa, attitude)
    get_steer_vec(manifold, attitude * doa)
end

function get_steer_vec(manifold, doa::Spherical, attitude)
    get_steer_vec(manifold, sph2cart(doa), attitude)
end

get_num_ants(::AbstractManifold{N}) where N = N


text:
using Documenter
using Blosc2

makedocs(
    sitename = "Blosc2.jl",
    format = Documenter.HTML(),
    modules = [Blosc2],
    pages = Any[
        "Introduction" => "index.md",
        "API" => [
            "Main API" => "api.md"
            "Low Level API" => "lib.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/waralex/Blosc2.jl.git",
    target = "build"
)


text:
experimentdesign_insts(ext::MatSpecExt) = [
    [(d, ext)
    for d in vcat(3, 25:25:100, 200:100:400)] # includes compile run
    ]

insts = OrderedDict()
insts["nat"] = (nothing, vcat(
    experimentdesign_insts(MatNegGeom()),
    experimentdesign_insts(MatNegSqrt()),
    experimentdesign_insts(MatNegSqrtConj()),
    experimentdesign_insts(MatNegPower01(1/3)),
    experimentdesign_insts(MatNegPower01Conj(1/3)),
    experimentdesign_insts(MatPower12(1.5)),
    ))
insts["ext"] = (nothing, vcat(
    experimentdesign_insts(MatNegGeomEFExp()),
    experimentdesign_insts(MatNegSqrtEigOrd()),
    experimentdesign_insts(MatNegSqrtConjDirect()),
    experimentdesign_insts(MatNegPower01EigOrd(1/3)),
    experimentdesign_insts(MatNegPower01ConjEigOrd(1/3)),
    experimentdesign_insts(MatPower12EigOrd(1.5)),
    ))
return (ExperimentDesignJuMP, insts)


text:
using Base: product
function evaluateEnergy!(options::zefiroOptions,sizes::zefiroSizes,data::zefiroData,output::zefiroOutput)
    AEP = zeros(sizes.winds,sizes.turbines)
    CP = zeros(sizes.winds,sizes.turbines)

    for i in 1:sizes.winds,j in 1:sizes.turbines
        AEP[i,j] = sum(data.RFwindspeed[:,i].*data.powercurve[:,j])*year2hour
        CP[i,j] = AEP[i,j]*percent/data.turbinePR[j]/year2hour
    end
    output.AEP = AEP
    output.CP = CP
    nothing
end

max_stars_repo_name: Thiago-NovaesB/zefiro

max_issues_repo_name: Thiago-NovaesB/zefiro

max_forks_repo_name: Thiago-NovaesB/zefiro


text:
function edge_to_support(s)
  vals = Dict{Tuple{Int64, Int64}, Float64}()
  I = Vector{Int64}()
  J = Vector{Int64}()
  V = Vector{Float64}()
  for e in 1:ne(s)
      de = elementary_duals(Val{1},s, e)
      ev = volume(Val{1}, s, e)
      dv = sum([dual_volume(Val{1}, s, d) for d in de])
      for d in de
          dt = incident(s, d, :D_∂e0)
          append!(I, dt)
          append!(J, fill(e, length(dt)))
          append!(V, fill(1/(dv*ev), length(dt)))
      end
  end
  sparse(I,J,V)
end

function tri_to_support(s)
  vals = Dict{Tuple{Int64, Int64}, Float64}()
  I = Vector{Int64}()
  J = Vector{Int64}()
  V = Vector{Float64}()
  for t in 1:ntriangles(s)
      dt = vcat(incident(s, incident(s, triangle_center(s, t), :D_∂v0), :D_∂e1)...)
      tv = volume(Val{2}, s, t)
      append!(I, dt)
      append!(J, fill(t, length(dt)))
      append!(V, fill(1/tv#= * sign(Val{2}, s, t)=#, length(dt)))
  end
  sparse(I,J,V)
end

function support_to_tri(s)
  vals = Dict{Tuple{Int64, Int64}, Float64}()
  I = Vector{Int64}()
  J = Vector{Int64}()
  V = Vector{Float64}()
  for t in 1:nv(s)
      dt = elementary_duals(Val{0},s, t)
      for d in dt
          push!(I, t)
          push!(J, d)
          push!(V, 1)
      end
  end
  sparse(I,J,V)
end

diag_vols(s) = spdiagm([dual_volume(Val{2}, s, dt) for dt in 1:nparts(s, :DualTri)])

wedge_mat(::Type{Val{(1,1)}}, s) = 2.0 * (support_to_tri(s)*diag_vols(s)*edge_to_support(s))
wedge_mat(::Type{Val{(2,0)}}, s) = support_to_tri(s)*diag_vols(s)*tri_to_support(s)

wedge_edge(::Type{Val{(1,1)}}, s) = dual_boundary(Val{2}, s) * 2.0 * (support_to_tri(s)*diag_vols(s)*edge_to_support(s))
wedge_edge(::Type{Val{(2,0)}}, s) = dual_boundary(Val{2}, s) * support_to_tri(s)*diag_vols(s)*tri_to_support(s)

function pd_wedge!(x, ::Type{Val{(1,1)}}, s, α, β; wedge_t = Dict((1,1)=>wedge_mat(Val{(1,1)}, s)), caches=[zeros(nv(s)), zeros(ne(s))], kw...)
  broadcast!(*, caches[2], α, β)
  mul!(x, wedge_t[(1,1)], caches[2])
end

function pd_wedge(::Type{Val{(1,1)}}, s, α, β; wedge_t = Dict((1,1)=>wedge_mat(Val{(1,1)}, s)), caches=[zeros(nv(s)), zeros(ne(s))], kw...)
  broadcast!(*, caches[2], α, β)
  wedge_t[(1,1)] * caches[2]
end

function pd_wedge(::Type{Val{(2,0)}}, s, α, β; wedge_t = Dict((2,0)=>wedge_mat(Val{(2,0)},s)), kw...)
  wedge_t[(2,0)] * (α .* β)
end

function pd_wedge(::Type{Val{(0,2)}}, s, α, β; wedge_t = nothing, kw...)
  α .* β
end

function init_wedge_ops(s)
  (d_mat=Dict(:d₁=>d(Val{1}, s), :dual_d₀=>dual_derivative(Val{0}, s)),
   wedge_e=Dict((1,1)=>wedge_edge(Val{(1,1)},s), (2,0)=>wedge_edge(Val{(2,0)},s)),
   wedge_t=Dict((1,1)=>wedge_mat(Val{(1,1)}, s), (2,0)=>wedge_mat(Val{(2,0)},s)),
   caches=[zeros(nv(s)), zeros(ne(s)), zeros(ntriangles(s))])
end

vect(s, e) = (s[s[e,:∂v1], :point] - s[s[e,:∂v0], :point]) * sign(1, s, e)
vect(s, e::AbstractVector) = [vect(s, el) for el in e]
t_vects(s,t) = vect(s, triangle_edges(s,t)) .* ([1,-1,1] * sign(2,s,t))
function comp_support(sd)
  vects = []
  for t in 1:ntriangles(sd)
      inds = triangle_edges(sd, t)
      for i in 1:3
          push!(vects, (t, inds[i]))
      end
  end
  v2comp = Dict{Tuple{Int64, Int64}, Int64}()
  for (i, v) in enumerate(vects)
      v2comp[v] = i
  end
  v2comp
end
function changes(s, v2comp)
  orient_vals = [1,-1,1]
  I = Vector{Int64}()
  J = Vector{Int64}()
  V = Vector{Float64}()
  for t in 1:ntriangles(s)
      inds = triangle_edges(s, t)
      e_vects = t_vects(s,t)
      vals = zeros(1:3)
      for i in 1:3
          ns = [(i+1)%3 + 1, i%3+1]
          ort = e_vects[i] × (e_vects[i] × e_vects[ns[1]])
          n_ort = ort / norm(ort)
          append!(J, v2comp[(t,inds[i])])
          append!(I, inds[ns[1]])
          append!(V, dot(n_ort, e_vects[ns[1]]) * orient_vals[ns[1]] * sign(1, s, ns[1])* orient_vals[i]* sign(2,s,t) / 3.0)
          append!(J, v2comp[(t,inds[i])])
          append!(I, inds[ns[2]])
          append!(V, dot(n_ort, e_vects[ns[2]]) * orient_vals[ns[2]] * sign(1, s, ns[2])* orient_vals[i]* sign(2,s,t) / 3.0)
      end
  end
  sparse(I,J,V, ne(s), ntriangles(s)*3)
end
function edge2comp(s, v2comp)
  I = Vector{Int64}()
  J = Vector{Int64}()
  V = Vector{Float64}()
  for t in 1:ntriangles(sd)
      inds = triangle_edges(sd, t)
      for i in 1:3
          push!(I, v2comp[(t,inds[i])])
          push!(J, inds[i])
          push!(V, 1 / volume(Val{1}, s, inds[i]))
      end
  end
  sparse(I,J,V)
end
function tri2comp(s, v2comp)
  I = Vector{Int64}()
  J = Vector{Int64}()
  V = Vector{Float64}()
  for t in 1:ntriangles(sd)
      inds = triangle_edges(sd, t)
      for i in 1:3
          push!(I, v2comp[(t,inds[i])])
          push!(J, t)
          push!(V, 1)
      end
  end
  sparse(I,J,V)
end

function cp_2_1!(x, α, β, matrices)
  mul!(matrices[:α_cache], matrices[:t2c], α)
  mul!(matrices[:β_cache], matrices[:e2c], β)
  broadcast!(*, matrices[:β_cache], matrices[:α_cache], matrices[:β_cache])
  mul!(x, matrices[:cross], matrices[:β_cache])
  #x .= matrices[:cross] * ((matrices[:t2c]*α).*(matrices[:e2c]*β))
end

function avg_mat(::Type{Val{(0,1)}},s)
  I = Vector{Int64}()
  J = Vector{Int64}()
  V = Vector{Float64}()
  for e in 1:ne(s)
      append!(J, [s[e,:∂v0],s[e,:∂v1]])
      append!(I, [e,e])
      append!(V, [0.5, 0.5])
  end
  sparse(I,J,V)
end

max_stars_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_issues_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_forks_repo_name: AlgebraicJulia/DECAPODES-Benchmarks


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 728d5546-7a96-11eb-1757-93681c4e7955
begin
	import Pkg
    Pkg.activate(".")
    Pkg.add("PlutoUI")
	Pkg.add("Plots")

	using PlutoUI
	using Plots
end

# ╔═╡ 8dfd7e3e-80b8-11eb-0a8b-99c7a0102380
md"""# SIR-Models

A collection of variations of SIR Models:

* Standard SIR with birth & death
* Inzidenz dependent infection rates
* 6D-SIR with two subpopulations

First set the time horizon and the step size for the euler method for all simulations below.

Start $(@bind T0 NumberField(0:10_000, default=0)) |
Step size Δt $(@bind Δt NumberField(0.0001:0.0001:1.0,default=0.01)) |
Stopp $(@bind Tend NumberField(0:10_000, default=150))

"""

# ╔═╡ e6d00224-8bfe-11eb-1b62-79ccf754d3d7
round(Int,1/Δt+0.1)

# ╔═╡ 62645b8c-7b7a-11eb-1b35-a3957a0252fb
md"""
---
### Classical SIR-Model

The standard SIR-Model with constant birth and death, infection and recovery rates.

${\frac{dS}{dt} = \nu N - \beta \frac{SI}{N} - \mu S}$

${\frac{dI}{dt} = \beta \frac{SI}{N} - \gamma I - \mu I}$

${\frac{dR}{dt} = \gamma I - \mu R \phantom{-\beta \frac{SI}{N}}}$

"""

# ╔═╡ e64e9d62-7b7b-11eb-3620-a7a57b954fcc
md"""
recovery rate γ $(@bind γ Slider(0.0:0.01:1,show_value=true,default=0.15)) |
infection rate β $(@bind β Slider(0.0:0.01:1,show_value=true,default=0.7))

death rate μ $(@bind μ Slider(0.0:0.01:1,show_value=true)) |
birth rate ν $(@bind ν Slider(0.0:0.01:1,show_value=true))
"""

# ╔═╡ 81608c56-80b9-11eb-00f0-3b50905ec998
md"""
---
### Incidence dependent infection rate

Here the infection rate $\beta=\beta(i)$ is dependent on the incidence in a population of constant size.

${\frac{dS}{dt} = - \beta(i) \frac{SI}{N} \phantom{- \gamma I} }$

${\frac{dI}{dt} = \phantom{-} \beta(i) \frac{SI}{N} - \gamma I}$

${\frac{dR}{dt} = \phantom{-\beta(i) \frac{SI}{N}-} \gamma I }$

with the following infection rate functions

* Constant: $\beta(i)\equiv\beta$
* Step Function: $\beta(i) = \beta (1 - \mathbb{1}_{[I_{c},\infty)}(I) (1-\frac{1}{r_s}))$
* Logistic Decay: $\beta(i) = \frac{\beta}{1+\exp(r_l(I-I_c))}$
* Inverse Proportional $\beta(i) = \beta - \frac{\beta}{I_c} I$

where $I_c$ is the critical incidence at which the measurements stat and $r_s, r_f$ is the factor resp. rate at which the infection rate decreases.
"""

# ╔═╡ 5e81dd2c-80b9-11eb-364e-63fa4ebd27da
md"""
Choice of Incidence Function $(@bind Β_name Select(["LogisticDecay", "StepFunction","Constant","InverseProportional"])) |
Delay for application of measurements $(@bind d NumberField(0:100,default=0))
"""

# ╔═╡ 4f8dee16-80ba-11eb-1d0f-f3ecdc9ac671
if Β_name == "LogisticDecay"
	md"""
	maximum infection rate β $(@bind b Slider(0.0:0.01:1,show_value=true,default=0.7))

	Maximum decrease of infection rate at $(@bind high_inc NumberField(0:1000,default=50)) per 10^5 | Rate of decrease $(@bind r_ld Slider(0.0:0.001:0.1,show_value=true,default=0.01))
	"""
elseif Β_name == "StepFunction"
	md"""
	high infection rate β $(@bind b Slider(0.0:0.01:1,show_value=true,default=0.7))

	Reduce infection rate if incidence is high by factor of: $(@bind high_inc_f NumberField(0:0.1:100,default=2)) |
	High incidence at $(@bind high_inc NumberField(0:1000,default=50)) per 10^5
	"""
elseif Β_name == "Constant"
	md"""
	infection rate β $(@bind b Slider(0.0:0.01:1,show_value=true,default=0.7))
	"""
elseif Β_name == "InverseProportional"
	md"""
	maximum infection rate β $(@bind b Slider(0.0:0.01:1,show_value=true,default=0.7)) |
	High incidence at $(@bind high_inc NumberField(0:1000,default=50)) per 10^5
	"""
else
	md""
end

# ╔═╡ 7da7ede2-8bf0-11eb-0829-49e179733bb7
begin
	#Infection Rate Functions
	StepFunction(i) = i ≥ high_inc ? b/high_inc_f : b
	LogisticDecay(i) = b / (1+exp(r_ld*((i-high_inc))))
	Constant(i) = b
	InverseProportional(i) = max(b - (b/high_inc)*i,0)

	function_choice = Dict(
		"StepFunction" => StepFunction,
		"InverseProportional" => InverseProportional,
		"LogisticDecay" => LogisticDecay,
		"Constant" => Constant
	)
end;

# ╔═╡ 59879918-8bea-11eb-25f1-0908e2cdad4b
begin
	function incidence(F::Array{Array{Float64,1},1},t;ndays=7,scale=100_000)
		#starting time
		t₀ = max(t-ndays,1)
		I₀ = t₀ > 1 ? F[t-ndays-1][2] : 0
		#set the range
		Days = @view F[t₀:t]
		#count the newly infected individuals
		I = sum(getindex.(Days,2))-I₀
		#calculate the rescaled, averaged total population size
		N = (sum(sum.(Days))/ndays) / scale
		#return Incidence
		return I/N
	end
	incidence(F::Array{Float64,1},t₀,ndays) = sum(view(F,t₀+1:t₀+ndays+1)) - F[t₀]
	function incidence(F::Array{Float64,1},pop_size;ndays=7,scale=100_000)
		I = cumsum(view(F,1:ndays))
		F_extended = vcat(F,fill(F[end],ndays+1))
		for t₀ ∈ ndays+1:length(F)
			append!(I,incidence(F_extended,t₀,ndays))
		end
		return I .* scale/pop_size
	end
end;

# ╔═╡ 79a2aad2-7c2a-11eb-00b0-6bebd6658619
md"""
---
### Six dimensional SIR-Model
"""

# ╔═╡ 968af7ba-7c2a-11eb-0aa4-6d51510020b1
md"""
Devide the population in two subgroups with different infection rates within and between these groups.

${\frac{dS_1}{dt} = - \beta_{11} \frac{S_1I_1}{N} - \beta_{12} \frac{S_1I_2}{N} \phantom{-\beta \frac{SI}{N}}}$

${\frac{dI_1}{dt} = \phantom{-} \beta_{11} \frac{S_1I_1}{N} + \beta_{12} \frac{S_1I_2}{N} -   \gamma_1 I_1 }$

${\frac{dR_1}{dt} = \phantom{- \beta_{11} \frac{S_1I_1}{N} - \beta_{12} \frac{S_1I_2}{N} -}\gamma_1 I_1 }$

${\frac{dS_2}{dt} = - \beta_{22} \frac{S_2I_2}{N} - \beta_{21} \frac{S_2I_1}{N} \phantom{-\beta \frac{SI}{N}}}$

${\frac{dI_2}{dt} = \phantom{-} \beta_{22} \frac{S_2I_2}{N} + \beta_{21} \frac{S_2I_1}{N} -   \gamma_2 I_2 }$

${\frac{dR_2}{dt} = \phantom{- \beta_{11} \frac{S_1I_1}{N} - \beta_{12} \frac{S_1I_2}{N} -}\gamma_2 I_2 }$

Here $\beta_{11}$ is the infection rate within group one, $\beta_{12}$ is the infection rate from group one to group two, and so on.
"""

# ╔═╡ cbf8bba4-7c26-11eb-3301-07233740f82a
md"""

β₁₁ = $(@bind β11 Slider(0.0:0.01:1.0, show_value=true, default=0.75)) |
β₁₂ = $(@bind β12 Slider(0.0:0.01:1.0, show_value=true, default=0.5))

β₂₁ = $(@bind β21 Slider(0.0:0.01:1.0, show_value=true, default=0.5)) |
β₂₂ = $(@bind β22 Slider(0.0:0.01:1.0, show_value=true, default=0.75))

γ₁ = $(@bind γ1 Slider(0.0:0.01:1.0, show_value=true, default=0.25)) |
γ₂ = $(@bind γ2 Slider(0.0:0.01:1.0, show_value=true, default=0.25))

"""

# ╔═╡ 32f4bc7c-80b9-11eb-0111-074cc68e892f
md"---"

# ╔═╡ ecad4e78-7a96-11eb-2b54-e7f910e83835
begin
	function euler(f,t0,tn,Δt,x₀)
		T = t0:Δt:tn
		F = Vector{typeof(x₀)}(undef,length(T))
		F[1] = x₀
		for (n,t) ∈ enumerate(T[2:end])
			F[n+1] = F[n] .+ f(t,F[n]) .* Δt
		end
		return F
	end

	function euler2(f,t0,tn,Δt,x₀)
		T = t0:Δt:tn
		F = Vector{typeof(x₀)}(undef,length(T))
		F[1] = x₀
		for (n,t) ∈ enumerate(T[2:end])
			F[n+1] = F[n] .+ f(t,F,n) .* Δt
		end
		return F
	end

	reshape_result(F) = [[f[i] for f ∈ F] for i ∈ 1:length(F[1])]

	T = T0:Δt:Tend
end;

# ╔═╡ 6b37eb20-7b7a-11eb-07de-a1042efa6721
begin
	function SIR(t,x)
		N = sum(x)
		S, I, R = x[1], x[2], x[3]
		return [
			(ν*N - β*S*I/N - μ*S),	# dS/dt
			(β*S*I/N - γ*I - μ*I),	# dI/dt
			(γ*I - μ*R)				# dR/dt
		]
	end

	SIR_0 = [9999.9,0.1,0.0]

	F = euler(SIR,T[1],T[end],Δt,SIR_0)
	rF = reshape_result(F)
	N = sum.(F)
end;

# ╔═╡ 68f66db2-7b7c-11eb-138d-2f40b3bfbb4e
let
	p = plot(framestlye=:zerolines,size=(680,450))

	plot!(p,T,rF,label=["S" "I" "R"])

	p
end

# ╔═╡ a9f94f00-8be9-11eb-36f9-d32e2ab1fca7
begin
	function SIR_H(t,F,n)
		#current state
		x = F[n]
		#current sizes
		S, I, R = x[1], x[2], x[3]
		#incidence rate
		n₀ = max(1,n-delay)
		β = Β(incidence(F,n₀))
		return [
			(- β*S*I/N_H ),		# dS/dt
			(β*S*I/N_H - γ*I),	# dI/dt
			(γ*I)				# dR/dt
		]
	end

	#---

	Β = function_choice[Β_name]
	delay = round(Int,1/Δt)*d

	SIR_H_0 = [9999.9,0.1,0.0]
	N_H = sum(SIR_H_0)

	F_H = euler2(SIR_H,T[1],T[end],Δt,SIR_H_0)
	rF_H = reshape_result(F_H)
end;

# ╔═╡ 59382348-8bf3-11eb-0fcf-73857b994092
let
	p = plot(framestlye=:zerolines,size=(680,450))

	Inc = incidence(rF_H[2],N_H)

	plot!(p,T,Inc,label="Incidence per 100.000",color=:red,legend=:topright)
	p2 = twinx()
	plot!(p2,T,Β.(Inc), label="Infection rate",color=:orange,legend=:bottomright)

	p
end

# ╔═╡ 58e3898c-7c21-11eb-3364-9168fcedb22c
begin
	function SIR6D(t,x)
		N = sum(x)
		S1, S2 = x[1], x[4]
		I1, I2 = x[2], x[5]
		R1, R2 = x[3], x[6]

		return [
			(-β11*S1*I1/N - β21*S1*I2/N),			# dS₁/dt
			(β11*S1*I1/N + β21*S1*I2/N - γ1*I1),	# dI₁/dt
			(γ1*I1),								# dR₁/dt
			(-β22*S2*I2/N - β12*S2*I1/N),			# dS₂/dt
			(β22*S2*I2/N + β12*S2*I1/N - γ2*I2),	# dI₂/dt
			(γ2*I2),								# dR₂/dt
		]
	end

	SIR6D_0 = [
		4999.0,1.0,0.0,
		5000.0,0.0,0.0,
	]

	F6D = euler(SIR6D,T[1],T[end],Δt,SIR6D_0)
	rF6D = reshape_result(F6D)
	N6D = sum.(F6D)
end;

# ╔═╡ d1644026-7c27-11eb-17ae-67c414382cdc
let
	p = plot(framestlye=:zerolines,size=(680,450))

	plot!(p,T,rF6D,
		#label=["S₁" "I₁" "R₁" "S₂" "I₂" "R₂"],
		label="",
		alpha=0.5,
		lw=1,
		color=[:darkblue :darkred :darkgreen :lightblue :orangered :lightgreen]
	)

	plot!(p,T,[rF6D[i] .+ rF6D[i+3] for i ∈ 1:3],
		label=["S" "I" "R"],
		color=[:blue :red :green],
		lw=2
	)

	p
end

# ╔═╡ Cell order:
# ╟─8dfd7e3e-80b8-11eb-0a8b-99c7a0102380
# ╠═e6d00224-8bfe-11eb-1b62-79ccf754d3d7
# ╟─62645b8c-7b7a-11eb-1b35-a3957a0252fb
# ╠═6b37eb20-7b7a-11eb-07de-a1042efa6721
# ╟─e64e9d62-7b7b-11eb-3620-a7a57b954fcc
# ╠═68f66db2-7b7c-11eb-138d-2f40b3bfbb4e
# ╟─81608c56-80b9-11eb-00f0-3b50905ec998
# ╟─5e81dd2c-80b9-11eb-364e-63fa4ebd27da
# ╟─4f8dee16-80ba-11eb-1d0f-f3ecdc9ac671
# ╟─59382348-8bf3-11eb-0fcf-73857b994092
# ╠═a9f94f00-8be9-11eb-36f9-d32e2ab1fca7
# ╠═7da7ede2-8bf0-11eb-0829-49e179733bb7
# ╠═59879918-8bea-11eb-25f1-0908e2cdad4b
# ╟─79a2aad2-7c2a-11eb-00b0-6bebd6658619
# ╟─968af7ba-7c2a-11eb-0aa4-6d51510020b1
# ╠═58e3898c-7c21-11eb-3364-9168fcedb22c
# ╟─cbf8bba4-7c26-11eb-3301-07233740f82a
# ╟─d1644026-7c27-11eb-17ae-67c414382cdc
# ╟─32f4bc7c-80b9-11eb-0111-074cc68e892f
# ╠═728d5546-7a96-11eb-1757-93681c4e7955
# ╠═ecad4e78-7a96-11eb-2b54-e7f910e83835

max_stars_repo_name: roccminton/PlutoNotebooks

max_issues_repo_name: roccminton/PlutoNotebooks

max_forks_repo_name: roccminton/PlutoNotebooks


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

using Test, FEMBase, FEMPlates

# Create 1x1 element, assemble stiffness matrix.

X = Dict(1 => [0.0, 0.0], 2 => [1.0, 0.0],
         3 => [1.0, 1.0], 4 => [0.0, 1.0])
element = Element(Quad4, [1, 2, 3, 4])
update!(element, "geometry", X)
update!(element, "youngs modulus", 4880.0*6.0*3.0*4.0)
update!(element, "thickness", 1.0)
update!(element, "poissons ratio", 1/3)
problem = Problem(MindlinPlate, "test plate", 3)
add_elements!(problem, [element])
time = 0.0
assemble!(problem, time)
K = Matrix(problem.assembly.K)
#display(K)
fixed_dofs = [1, 2, 3, 10, 11, 12]
free_dofs = setdiff(1:12, fixed_dofs)
f = zeros(12)
f[[4, 7]] .= 10.0
u = zeros(12)
u[free_dofs] = K[free_dofs, free_dofs] \ f[free_dofs]
#display(u)
u_expected = [
0.0
0.0
0.0
0.000306921
-0.000249543
-0.0000458344
0.000306921
-0.000249543
0.0000458344
0.0
0.0
0.0
]
@test isapprox(u, u_expected,atol=1.0e-8)


text:
""
function objective_min_fuel_cost(pm::_PM.AbstractPowerModel)
    model = _PM.check_cost_models(pm)
    if model == 1
        return objective_min_pwl_fuel_cost(pm)
    elseif model == 2
        return objective_min_polynomial_fuel_cost(pm)
    else
        error("Only cost models of types 1 and 2 are supported at this time, given cost model type of $(model)")
    end

end

""
function objective_min_polynomial_fuel_cost(pm::_PM.AbstractPowerModel)
    order = _PM.calc_max_cost_index(pm.data)-1

    if order == 1
        return _objective_min_polynomial_fuel_cost_linear(pm)
    elseif order == 2
        return _objective_min_polynomial_fuel_cost_quadratic(pm)
    else
        error("cost model order of $(order) is not supported")
    end
end

function _objective_min_polynomial_fuel_cost_linear(pm::_PM.AbstractPowerModel)
    from_idx = Dict()
    for (n, nw_ref) in nws(pm)
        from_idx[n] = Dict(arc[1] => arc for arc in nw_ref[:arcs_from_dc])
    end

    return JuMP.@objective(pm.model, Min,
        sum(
            sum(   gen["cost"][1]*sum( var(pm, n, c, :pg, i) for c in conductor__PM.ids(pm, n))+
                   gen["cost"][2] for (i,gen) in nw_ref[:gen])
        for (n, nw_ref) in nws(pm))
    )
end
""
function _objective_min_polynomial_fuel_cost_quadratic(pm::_PM.AbstractPowerModel)
    from_idx = Dict()
    for (n, nw_ref) in nws(pm)
        from_idx[n] = Dict(arc[1] => arc for arc in nw_ref[:arcs_from_dc])
    end

    return JuMP.@objective(pm.model, Min,
        sum(
            sum(   gen["cost"][1]*sum( var(pm, n, c, :pg, i) for c in conductor__PM.ids(pm, n))^2 +
                   gen["cost"][2]*sum( var(pm, n, c, :pg, i) for c in conductor__PM.ids(pm, n))+
                   gen["cost"][3] for (i,gen) in nw_ref[:gen])
        for (n, nw_ref) in nws(pm))
    )
end
""
function objective_min_polynomial_fuel_cost(pm::_PM.AbstractConicModel)
    _PM.check_polynomial_cost_models(pm)
    pg_sqr = Dict()
    for (n, nw_ref) in _PM.nws(pm)
        for cnd in _PM.conductor__PM.ids(pm, n)
            pg_sqr = _PM.var(pm, n, cnd)[:pg_sqr] = JuMP.@variable(pm.model,
                [i in _PM.ids(pm, n, :gen)], base_name="$(n)_$(cnd)_pg_sqr",
                lower_bound = _PM.ref(pm, n, :gen, i, "pmin", cnd)^2,
                upper_bound = _PM.ref(pm, n, :gen, i, "pmax", cnd)^2
            )
            for (i, gen) in nw_ref[:gen]
                JuMP.@constraint(pm.model, [pg_sqr[i], var(pm, n, cnd, :pg, i)/sqrt(2), var(pm, n, cnd, :pg, i)/sqrt(2)] in JuMP.SecondOrderCone())
            end

        end
    end
    return JuMP.@objective(pm.model, Min,
        sum(
            sum( gen["cost"][1]*sum(_PM.var(pm, n, cnd, :pg_sqr, i) for cnd in _PM.conductor__PM.ids(pm, n)) +
                 gen["cost"][2]*sum(_PM.var(pm, n, cnd, :pg, i) for cnd in _PM.conductor__PM.ids(pm, n)) +
                 gen["cost"][3] for (i,gen) in nw_ref[:gen])
        for (n, nw_ref) in _PM.nws(pm))
    )
end
""
function objective_min_pwl_fuel_cost(pm::_PM.AbstractPowerModel)

    for (n, nw_ref) in _PM.nws(pm)
        pg_cost = _PM.var(pm, n)[:pg_cost] = JuMP.@variable(pm.model,
            [i in _PM.ids(pm, n, :gen)], base_name="$(n)_pg_cost"
        )

        # pwl cost
        gen_lines = _PM.get_lines(nw_ref[:gen])
        for (i, gen) in nw_ref[:gen]
            for line in gen_lines[i]
                JuMP.@constraint(pm.model, pg_cost[i] >= line["slope"]*sum(_PM.var(pm, n, cnd, :pg, i) for cnd in _PM.conductor__PM.ids(pm, n)) + line["intercept"])
            end
        end

    end

    return JuMP.@objective(pm.model, Min,
        sum(
            sum( _PM.var(pm, n,:pg_cost, i) for (i,gen) in nw_ref[:gen])
        for (n, nw_ref) in _PM.nws(pm))
        )
end

##################### TNEP Objective   ###################
function objective_min_cost(pm::_PM.AbstractPowerModel)
        gen_cost = Dict()
        for (n, nw_ref) in _PM.nws(pm)
            for (i,gen) in nw_ref[:gen]
                pg = _PM.var(pm, n, :pg, i)

                if length(gen["cost"]) == 1
                    gen_cost[(n,i)] = gen["cost"][1]
                elseif length(gen["cost"]) == 2
                    gen_cost[(n,i)] = gen["cost"][1]*pg + gen["cost"][2]
                elseif length(gen["cost"]) == 3
                    gen_cost[(n,i)] = gen["cost"][2]*pg + gen["cost"][3]
                else
                    gen_cost[(n,i)] = 0.0
                end
            end
        end

        return JuMP.@objective(pm.model, Min,
            sum(
                sum(conv["cost"]*_PM.var(pm, n, :conv_ne, i) for (i,conv) in nw_ref[:convdc_ne])
                +
                sum(branch["cost"]*_PM.var(pm, n, :branchdc_ne, i) for (i,branch) in nw_ref[:branchdc_ne])
                +
                sum( gen_cost[(n,i)] for (i,gen) in nw_ref[:gen] )
                for (n, nw_ref) in _PM.nws(pm)
                    )
        )
end

function objective_min_cost_acdc(pm::_PM.AbstractPowerModel)
        gen_cost = Dict()
        for (n, nw_ref) in _PM.nws(pm)
            for (i,gen) in nw_ref[:gen]
                pg = _PM.var(pm, n, :pg, i)

                if length(gen["cost"]) == 1
                    gen_cost[(n,i)] = gen["cost"][1]
                elseif length(gen["cost"]) == 2
                    gen_cost[(n,i)] = gen["cost"][1]*pg + gen["cost"][2]
                elseif length(gen["cost"]) == 3
                    gen_cost[(n,i)] = gen["cost"][2]*pg + gen["cost"][3]
                else
                    gen_cost[(n,i)] = 0.0
                end
            end
        end

        return JuMP.@objective(pm.model, Min,
            sum(
                sum(conv["cost"]*_PM.var(pm, n, :conv_ne, i) for (i,conv) in nw_ref[:convdc_ne])
                +
                sum(branch["construction_cost"]*_PM.var(pm, n, :branchdc_ne, i) for (i,branch) in nw_ref[:ne_branch])
                +
                sum(branch["cost"]*_PM.var(pm, n, :branchdc_ne, i) for (i,branch) in nw_ref[:branchdc_ne])
                +
                sum( gen_cost[(n,i)] for (i,gen) in nw_ref[:gen] )
                for (n, nw_ref) in _PM.nws(pm)
                    )
        )
end

# function objective_min_cost_OPF(pm::_PM.AbstractPowerModel)
#         base_nws = pm.setting["base_list"]
#         cont_nws = pm.setting["Cont_list"]
#         Total_sample = pm.setting["Total_sample"]
#         curt_gen = pm.setting["curtailed_gen"]
#         max_curt = pm.setting["max_curt"]
#
#         gen_cost = Dict()
#         FFR_cost = Dict()
#         FCR_cost = Dict()
#         fail_prob = Dict()
#         Gen_cost = _PM.var(pm)[:Gen_cost] = JuMP.@variable(pm.model, start = 0)
#         FFRReserves = _PM.var(pm)[:FFR_Reserves] = JuMP.@variable(pm.model, start = 0)
#         FCRReserves = _PM.var(pm)[:FCR_Reserves] = JuMP.@variable(pm.model, start = 0)
#         Curt = _PM.var(pm)[:curt] = JuMP.@variable(pm.model, start = 0)
#
#         sol_component_value_mod_wonw(pm,   :Gen_cost, Gen_cost)
#         sol_component_value_mod_wonw(pm,   :FFR_Reserves, FFRReserves)
#         sol_component_value_mod_wonw(pm,   :FCR_Reserves, FCRReserves)
#
#         for (n, nw_ref) in _PM.nws(pm)
#             for (r, reserves) in nw_ref[:reserves]
#                 FFR_cost[(n,r)] = reserves["Cf"]
#                 FCR_cost[(n,r)] = reserves["Cg"]
#             end
#             for (i,gen) in nw_ref[:gen]
#                 pg = _PM.var(pm, n, :pg, i)
#                 if length(gen["cost"]) == 1
#                     gen_cost[(n,i)] = gen["cost"][1]
#                 elseif length(gen["cost"]) == 2
#                     gen_cost[(n,i)] = gen["cost"][1]*pg + gen["cost"][2]
#                 elseif length(gen["cost"]) == 3
#                     gen_cost[(n,i)] = gen["cost"][2]*pg + gen["cost"][3]
#                 else
#                     gen_cost[(n,i)] = 0.0
#                 end
#             end
#         end
#
#         # display(base_nws)
#         # display(cont_nws)
#         Scale = 8760/Total_sample*5 # 5 for no. of gap years between two time steps
#         # multiperiod
#         JuMP.@constraint(pm.model, Gen_cost == sum(sum(Scale/10^6*gen_cost[(b,i)] for (i,gen) in _PM.nws(pm)[b][:gen]) for b in base_nws) )
#         JuMP.@constraint(pm.model, FFRReserves == sum(FFR_cost[(c,2)]*100*Scale/10^6* _PM.var(pm, c, :Pff, 2) for (b,c,br) in cont_nws) )
#         JuMP.@constraint(pm.model, FCRReserves ==  sum(FCR_cost[(c,2)]*100*Scale/10^6*_PM.var(pm, c, :Pgg, 2)  for (b,c,br) in cont_nws) )
#
#       curtailment = Dict()
#       capacity = Dict()
#       for b in base_nws
#           for c in curt_gen
#             curtailment[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c)) - _PM.var(pm, b, :pg, c)
#             capacity[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c))
#          end
#       end
#       JuMP.@constraint(pm.model, Curt == sum(sum(Scale*curtailment[(b,c)] for c in curt_gen) for b in base_nws) / sum(sum(Scale*capacity[(b,c)] for c in curt_gen) for b in base_nws) )
#       JuMP.@constraint(pm.model, Curt <= max_curt)
#       return  JuMP.@objective(pm.model, Min, Gen_cost + FFRReserves + FCRReserves)
#   end
#
#
  function sol_component_value_mod_wonw(pm::_PM.AbstractPowerModel, comp_name::Symbol, variables)
      _PM.sol(pm)[comp_name] = variables
  end

  function objective_min_cost_OPF(pm::_PM.AbstractPowerModel)
          base_nws = pm.setting["base_list"]
          cont_nws = pm.setting["Cont_list"]
          Total_sample = pm.setting["Total_sample"]
          curt_gen = pm.setting["curtailed_gen"]
          max_curt = pm.setting["max_curt"]
          weights = pm.setting["weights"]

          gen_cost = Dict()
          FFR_cost = Dict()
          FCR_cost = Dict()
          fail_prob = Dict()
          Gen_cost = _PM.var(pm)[:Gen_cost] = JuMP.@variable(pm.model, start = 0)
          FFRReserves = _PM.var(pm)[:FFR_Reserves] = JuMP.@variable(pm.model, start = 0)
          FCRReserves = _PM.var(pm)[:FCR_Reserves] = JuMP.@variable(pm.model, start = 0)
          Cont = _PM.var(pm)[:Cont]     = JuMP.@variable(pm.model, start = 0)
          Curt = _PM.var(pm)[:curt] = JuMP.@variable(pm.model, start = 0)

          sol_component_value_mod_wonw(pm,   :Gen_cost, Gen_cost)
          sol_component_value_mod_wonw(pm,   :FFR_Reserves, FFRReserves)
          sol_component_value_mod_wonw(pm,   :FCR_Reserves, FCRReserves)
          sol_component_value_mod_wonw(pm,   :Cont, Cont)

          for (n, nw_ref) in _PM.nws(pm)
              for (r, reserves) in nw_ref[:reserves]
                  FFR_cost[(n,r)] = reserves["Cf"]
                  FCR_cost[(n,r)] = reserves["Cg"]
              end
              for (i,gen) in nw_ref[:gen]
                  pg = _PM.var(pm, n, :pg, i)
                  if length(gen["cost"]) == 1
                      gen_cost[(n,i)] = gen["cost"][1]
                  elseif length(gen["cost"]) == 2
                      gen_cost[(n,i)] = gen["cost"][1]*pg + gen["cost"][2]
                  elseif length(gen["cost"]) == 3
                      gen_cost[(n,i)] = gen["cost"][2]*pg + gen["cost"][3]
                  else
                      gen_cost[(n,i)] = 0.0
                  end
              end
          end
          Gen_cost = _PM.var(pm)[:Gen_cost]
          # Reserves = _PM.var(pm)[:Reserves]
          Cont = _PM.var(pm)[:Cont]

          display("#####################################cost of conv and branch")

          display(base_nws)
          display(cont_nws)
          Scale = 8760/Total_sample*5 # 5 for no. of gap years between two time steps
          # multiperiod
          display(JuMP.@constraint(pm.model, Gen_cost == sum(sum(Scale/10^6*gen_cost[(b,i)] for (i,gen) in _PM.nws(pm)[b][:gen] *weights[b]) for b in base_nws) ))

          display(JuMP.@constraint(pm.model, FFRReserves == sum(weights[b]*FFR_cost[(c,2)]*100*Scale/10^6* _PM.var(pm, c, :Pff, 2) for (b,c,br) in cont_nws) ))

         display(JuMP.@constraint(pm.model, FCRReserves ==  sum(weights[b]*FCR_cost[(c,2)]*100*Scale/10^6*_PM.var(pm, c, :Pgg, 2)  for (b,c,br) in cont_nws) ) )

         display(JuMP.@constraint(pm.model, Cont == sum( sum( (Scale/10^6) *gen_cost[(c,i)]*_PM.ref(pm, b, :branchdc, br)["fail_prob"]
                      for (i,gen) in  _PM.nws(pm)[c][:gen])
                      for (b,c,br) in cont_nws) ) )

        curtailment = Dict()
        capacity = Dict()

        for b in base_nws
            for c in curt_gen
              curtailment[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c)) - _PM.var(pm, b, :pg, c)
              capacity[(b,c)] = JuMP.upper_bound(_PM.var(pm, b, :pg, c))
           end
       end
       display(curtailment)
       display(capacity)

        display(JuMP.@constraint(pm.model, Curt == sum(sum(Scale*curtailment[(b,c)] for c in curt_gen) for b in base_nws) / sum(sum(Scale*capacity[(b,c)] for c in curt_gen) for b in base_nws) ) )
        display(JuMP.@constraint(pm.model, Curt <= max_curt))
         # singleperiod
        #  display(JuMP.@constraint(pm.model, Gen_cost == sum(sum(8760*30/10^6*gen_cost[(b,i)] for (i,gen) in _PM.nws(pm)[b][:gen]) for b in base_nws) ))
        #
        #  display(JuMP.@constraint(pm.model, FFRReserves == sum(FFR_cost[(c,2)]*100*30*8760/10^6* _PM.var(pm, c, :Pff, 2) for (b,c,br) in cont_nws) ))
        #
        # display(JuMP.@constraint(pm.model, FCRReserves ==  sum(FCR_cost[(c,2)]*100*30*8760/10^6*_PM.var(pm, c, :Pgg, 2)  for (b,c,br) in cont_nws) ) )
        #
        # display( JuMP.@constraint(pm.model, Cont == sum( sum( (8760*30/10^6) *gen_cost[(c,i)]*_PM.ref(pm, b, :branchdc, br)["fail_prob"]
        #              for (i,gen) in  _PM.nws(pm)[c][:gen])
        #              for (b,c,br) in cont_nws) ) )

         #  display(JuMP.@constraint(pm.model, Gen_cost == sum(8760*30/10^6*gen_cost[(1,i)] for (i,gen) in _PM.nws(pm)[1][:gen]) ) )
         #  display(JuMP.@constraint(pm.model, FFRReserves   == sum(FFR_cost[(n,2)]*100*30*8760/10^6* _PM.var(pm, n, :Pff, 2) *nw_ref[:branchdc][n-1]["fail_prob"] for (n, nw_ref) in _PM.nws(pm) if n > 1) ))
         #
         # display(JuMP.@constraint(pm.model, FCRReserves   == sum(FCR_cost[(n,2)]*100*30*8760/10^6*_PM.var(pm, n, :Pgg, 2) *nw_ref[:branchdc][n-1]["fail_prob"]  for (n, nw_ref) in _PM.nws(pm) if n > 1) ))
         #
         #
         #  display(JuMP.@constraint(pm.model, Cont     == sum( sum( (8760*30/10^6) *gen_cost[(n,i)]*nw_ref[:branchdc][n-1]["fail_prob"]
         #              for (i,gen) in nw_ref[:gen])   for (n, nw_ref) in _PM.nws(pm) if n > 1) ) )

          return display( JuMP.@objective(pm.model, Min, Gen_cost + FFRReserves + FCRReserves))
  end


text:
using PETScBinaryIO
using Test
using SparseArrays

function idempotent_io(int_type, scalar_type, itype_check, stype_check)
    try
        itype, stype = if int_type == nothing
            PETScBinaryIO.parse_petsc_config()
        else
            int_type, scalar_type
        end
        A = rand(stype, 1000)
        B = sprand(stype, 1000,1000,0.1)
        C = rand(itype, 1000)
        writepetsc("vec.petsc", [A, B, C]; int_type=int_type, scalar_type=scalar_type)
        AA, BB, CC = readpetsc("vec.petsc"; int_type=int_type, scalar_type=scalar_type)
        AA == A && BB == B && CC == C && eltype(A) == stype_check && eltype(C) == itype_check
    finally
        if isfile("vec.petsc")
            rm("vec.petsc")
        end
    end
end

@testset "Test Suite" begin
    @testset "Idempotent IO" begin
        for itype in [Int32, Int64]
            for stype in [Float32, Float64]
                @testset "$itype $stype" begin
                    @test idempotent_io(itype, stype, itype, stype)
                end

                @testset "$itype $stype config" begin
                    withenv("PETSC_DIR" => pwd(), "PETSC_ARCH" => "arch-i$(sizeof(itype) * 8)-f$(sizeof(stype) * 8)") do
                        @test idempotent_io(nothing, nothing, itype, stype)
                    end
                end
            end
        end
    end

    @testset "Parse Config" begin
        for itype in [Int32, Int64]
            for stype in [Float32, Float64]
                @testset "$itype $stype config" begin
                    withenv("PETSC_DIR" => pwd(), "PETSC_ARCH" => "arch-i$(sizeof(itype) * 8)-f$(sizeof(stype) * 8)") do
                        @test PETScBinaryIO.parse_petsc_config() == (itype, stype)
                    end
                end
            end
        end
    end
end


text:
"""
    NumaAllocators.LibNUMA

NUMA Support for Linux

See https://github.com/numactl/numactl
"""
module LibNUMA

using ..NumaAllocators: AbstractNumaAllocator
import ArrayAllocators: AbstractArrayAllocator, nbytes, allocate
import ArrayAllocators: lineage_finalizer

using NUMA_jll

@static if ( VERSION >= v"1.6" && NUMA_jll.is_available() ) || isdefined(NUMA_jll, :libnuma)
    const libnuma = NUMA_jll.libnuma
else
    const libnuma = nothing
end

function numa_alloc_onnode(size, node)
    return ccall((:numa_alloc_onnode, libnuma), Ptr{Nothing}, (Csize_t, Cint), size, node)
end

function numa_alloc_local(size)
    return ccall((:numa_alloc_local, libnuma), Ptr{Nothing}, (Csize_t,), size)
end

function numa_free(arr::AbstractArray)
    numa_free(arr, sizeof(arr))
end

function numa_free(mem, size)
    return ccall((:numa_free, libnuma), Nothing, (Ptr{Nothing}, Csize_t), mem, size)
end

function numa_num_task_nodes()
    return ccall((:numa_numa_task_nodes, libnuma), Cint, ())
end

function numa_max_node()
    return ccall((:numa_max_node, libnuma), Cint, ())
end

function numa_node_of_cpu(cpu = sched_getcpu())
    return ccall((:numa_node_of_cpu, libnuma), Cint, (Cint,), cpu)
end

function sched_getcpu()
    return ccall(:sched_getcpu, Cint, ())
end

function wrap_numa(::Type{ArrayType}, ptr::Ptr{T}, dims) where {T, ArrayType <: AbstractArray{T}}
    if ptr == C_NULL
        throw(OutOfMemoryError())
    end
    arr = unsafe_wrap(ArrayType, ptr, dims; own = false)
    lineage_finalizer(numa_free, arr)
    return arr
end

abstract type AbstractLibNumaAllocator{B} <: AbstractNumaAllocator{B} end
Base.unsafe_wrap(::AbstractLibNumaAllocator, args...) = wrap_numa(args...)

"""
    LibNumaAllocator{B}

Allocate memory via `numa_alloc_onnode`.

See https://linux.die.net/man/3/numa
"""
struct LibNumaAllocator{B} <: AbstractLibNumaAllocator{B}
    node::Cint
end

function allocate(n::LibNumaAllocator, num_bytes)
    return numa_alloc_onnode(num_bytes, n.node)
end

end # module LibNUMA


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

# # Model Specifications

# ## Host Specifications
#
# The host may be in one of 4 states:
#
# - susceptible (S_H)
#
# - incubation period (E_H)
#
# - infectious (I_H)
#
# - recovered (R_H)
#
# ## Vector Specifications
#
# The vector may be in one of 3 states:
#
# - susceptible (S_V)
#
# - incubation period (E_V)
#
# - infectious (I_V)

# # Model Parameters
#
# σ_H, σ_V: the incubation rates for hosts & vectors (units: per time)
#
# μ_H, μ_V: the mortality rates for hosts & vectors (units: per time)
#
# λ: the clearance (or recovery) rate for hosts (units: per time)
#
# β: the infection rate (units: per capita per time)

# +
function F(du,u,p,t)
    S_H, E_H, I_H, R_H, S_V, E_V, I_V = u

    # host dynamics
    host_infection = (p.β*S_H*I_V)/p.N_H
    host_mortality = p.μ_H .* u[1:4] # include S_H, so easier to remove mortality
    host_births = sum(host_mortality)
    host_progression = p.σ_H*E_H
    recovery = p.λ*I_H

    du[1] = -host_infection + host_births
    du[2] = host_infection - host_progression
    du[3] = host_progression - recovery
    du[4] = recovery
    du[1:4] -= host_mortality

    # vector dynamics
    vec_infection = (p.β*S_V*I_H)/p.N_H
    vec_mortality = p.μ_V .* u[5:7] # include S_V, so easier to remove mortality
    vec_births = sum(vec_mortality)
    vec_progression = p.σ_V*E_V

    du[5] = -vec_infection + vec_births
    du[6] = vec_infection - vec_progression
    du[7] = vec_progression
    du[5:7] -= vec_mortality

end

# +
using DifferentialEquations
using IterableTables, DataFrames

u0 = [
    100.0,  0.0, 1.0, 0.0,
    10000.0, 0.0, 0.0
]
p = (
  μ_H=1/365, μ_V=1/30, σ_H=1/3, σ_V=1/7, λ=1/14,
  β=0.05, N_H = sum(u0[1:4])
)
tspan = (0.0, 365.0)
prob = ODEProblem(F, u0, tspan, p)
sol = @time solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8,saveat=range(0,stop = 365,length = 365*10+1))

df = DataFrame(sol')
rename!(df,:x1 => :S_H,:x2 => :E_H, :x3 => :I_H, :x4 => :R_H,
  :x5 => :S_V, :x6 => :E_V, :x7 => :I_V)
df.t = collect(range(0,stop = 365,length = 365*10+1));
# -

first(df,10)

using StatPlots

plot(df.t, [df.S_H,df.E_H,df.I_H,df.R_H],xlabel="t",ylabel="value",title="Host Dynamics",label=["S_H","E_H","I_H","R_H"])

plot(df.t,[df.S_V,df.E_V,df.I_V],xlabel="t",ylabel="value",title="Virus Dynamics",label=["S_V","E_V","I_V"])

max_stars_repo_name: epirecipes/epiRecipes

max_issues_repo_name: epirecipes/epiRecipes

max_forks_repo_name: epirecipes/epiRecipes


text:
using Rasters, Test, Dates, Plots

ga2 = Raster(ones(91) * (-25:15)', (X(0.0:4.0:360.0), Y(-25.0:1.0:15.0), ); name=:Test)
ga3 = Raster(rand(10, 41, 91), (Z(100:100:1000), Y(-20.0:1.0:20.0), X(0.0:4.0:360.0)))
ga4ti = Raster(
    rand(10, 41, 91, 4),
    (Z(100:100:1000), Y(-20.0:1.0:20.0), X(0.0:4.0:360.0), Ti(Date(2001):Year(1):Date(2004)))
)
ga4x = Raster(
    rand(10, 41, 91, 4),
    (Z(100:100:1000), Y(-20.0:1.0:20.0), X(0.0:4.0:360.0), X())
)

plot(ga2)
plot(ga3[Y(At(0.0))])
plot(ga3[X(At(180.0))])
# Line plot with Z on the vertical axis
plot(ga3[X(At(0.0)), Y(At(0.0))])
# DD fallback line plot with Z as key (not great really)
plot(ga4x[X(At(0.0)), Y(At(0.0))])
# DD fallback heatmap with Z on Y axis
heatmap(ga4x[X(At(0.0)), Y(At(0.0))])
# Cant plot 4d
@test_throws ErrorException plot(ga4x)
# 3d plot by NoLookupArray X dim
plot(ga4x[Y(1)])
# 3d plot by Ti dim
plot(ga4ti[Z(1)])
# Rasters handles filled contours
contourf(ga2)
# RasterStack plot
plot(RasterStack(ga2, ga3))

# DD fallback
contour(ga2)


text:
function saveDFG(dfg::G, folder::String) where G <: AbstractDFG
    variables = getVariables(dfg)
    factors = getFactors(dfg)
    varFolder = "$folder/variables"
    factorFolder = "$folder/factors"
    # Folder preparations
    if !isdir(folder)
        @info "Folder '$folder' doesn't exist, creating..."
        mkpath(folder)
    end
    !isdir(varFolder) && mkpath(varFolder)
    !isdir(factorFolder) && mkpath(factorFolder)
    # Clearing out the folders
    map(f -> rm("$varFolder/$f"), readdir(varFolder))
    map(f -> rm("$factorFolder/$f"), readdir(factorFolder))
    # Variables
    for v in variables
        vPacked = packVariable(dfg, v)
        io = open("$varFolder/$(v.label).json", "w")
        JSON2.write(io, vPacked)
        close(io)
    end
    # Factors
    for f in factors
        fPacked = packFactor(dfg, f)
        io = open("$folder/factors/$(f.label).json", "w")
        JSON2.write(io, fPacked)
        close(io)
    end
end

function loadDFG(folder::String, iifModule, dfgLoadInto::G) where G <: AbstractDFG
    variables = DFGVariable[]
    factors = DFGFactor[]
    varFolder = "$folder/variables"
    factorFolder = "$folder/factors"
    # Folder preparations
    !isdir(folder) && error("Can't load DFG graph - folder '$folder' doesn't exist")
    !isdir(varFolder) && error("Can't load DFG graph - folder '$folder' doesn't exist")
    !isdir(factorFolder) && error("Can't load DFG graph - folder '$folder' doesn't exist")

    varFiles = readdir(varFolder)
    factorFiles = readdir(factorFolder)
    for varFile in varFiles
        io = open("$varFolder/$varFile")
        packedData = JSON2.read(io, Dict{String, Any})
        push!(variables, unpackVariable(dfgLoadInto, packedData))
    end
    @info "Loaded $(length(variables)) variables - $(map(v->v.label, variables))"
    @info "Inserting variables into graph..."
    # Adding variables
    map(v->addVariable!(dfgLoadInto, v), variables)

    for factorFile in factorFiles
        io = open("$factorFolder/$factorFile")
        packedData = JSON2.read(io, Dict{String, Any})
        push!(factors, unpackFactor(dfgLoadInto, packedData, iifModule))
    end
    @info "Loaded $(length(variables)) factors - $(map(f->f.label, factors))"
    @info "Inserting factors into graph..."
    # # Adding factors
    map(f->addFactor!(dfgLoadInto, f._variableOrderSymbols, f), factors)

    # Finally, rebuild the CCW's for the factors to completely reinflate them
    @info "Rebuilding CCW's for the factors..."
    for factor in factors
        iifModule.rebuildFactorMetadata!(dfgLoadInto, factor)
    end

    # PATCH - To update the fncargvID for factors, it's being cleared somewhere in rebuildFactorMetadata.
    # TEMPORARY
    # TODO: Remove in future
    map(f->solverData(f).fncargvID = f._variableOrderSymbols, getFactors(dfgLoadInto))

    return dfgLoadInto
end


text:
"Utilities"
module Util

using Base: sym_in, merge_names, merge_types
using Spec

export applymany, ntranspose, Counter, reset!, increment!, UTuple, *ₛ

include("misc.jl")  # Miscellaneous
include("specs.jl") # Domain General specification tools

end


text:
### Cache
struct CodeCache
    dict::Dict{MethodInstance,Vector{CodeInstance}}
    callback::Function
    CodeCache(callback) = new(Dict{MethodInstance,Vector{CodeInstance}}(), callback)
end

function Base.show(io::IO, ::MIME"text/plain", cc::CodeCache)
    print(io, "CodeCache: ")
    for (mi, cis) in cc.dict
        println(io)
        print(io, "- ")
        show(io, mi)

        for ci in cis
            println(io)
            print(io, "  - ")
            print(io, (ci.min_world, ci.max_world))
        end
    end
end

function Core.Compiler.setindex!(cache::CodeCache, ci::CodeInstance, mi::MethodInstance)
    if !isdefined(mi, :callbacks)
        mi.callbacks = Any[cache.callback]
    else
        # Check if callback is present
        if all(cb -> cb !== cache.callback, mi.callbacks)
            push!(mi.callbacks, cache.callback)
        end
    end

    cis = get!(cache.dict, mi, CodeInstance[])
    push!(cis, ci)
end

### world view of the cache

using Core.Compiler: WorldView

function Core.Compiler.haskey(wvc::WorldView{CodeCache}, mi::MethodInstance)
    Core.Compiler.get(wvc, mi, nothing) !== nothing
end

function Core.Compiler.get(wvc::WorldView{CodeCache}, mi::MethodInstance, default)
    cache = wvc.cache
    for ci in get!(cache.dict, mi, CodeInstance[])
        if ci.min_world <= wvc.worlds.min_world && wvc.worlds.max_world <= ci.max_world
            # TODO: if (code && (code == jl_nothing || jl_ir_flag_inferred((jl_array_t*)code)))
            return ci
        end
    end

    return default
end

function Core.Compiler.getindex(wvc::WorldView{CodeCache}, mi::MethodInstance)
    r = Core.Compiler.get(wvc, mi, nothing)
    r === nothing && throw(KeyError(mi))
    return r::CodeInstance
end

Core.Compiler.setindex!(wvc::WorldView{CodeCache}, ci::CodeInstance, mi::MethodInstance) =
    Core.Compiler.setindex!(wvc.cache, ci, mi)

# invalidation
# invalidate_method_instance, but for our cache
function invalidate(cache::CodeCache, replaced::MethodInstance, max_world, depth)
    cis = get(cache.dict, replaced, nothing)
    if cis === nothing
        return
    end
    for ci in cis
        if ci.max_world == ~0 % Csize_t
            @assert ci.min_world - 1 <= max_world "attempting to set illogical constraints"
            ci.max_world = max_world
        end
        @assert ci.max_world <= max_world
    end

    # recurse to all backedges to update their valid range also
    if isdefined(replaced, :backedges)
        backedges = replaced.backedges
        # Don't touch/empty backedges `invalidate_method_instance` in C will do that later
        # replaced.backedges = Any[]

        for mi in backedges
            invalidate(cache, mi, max_world, depth + 1)
        end
    end
end


text:
module Semirings
using ..SuiteSparseGraphBLAS: isGxB, isGrB, TypedSemiring, AbstractSemiring, splitconstant
using ..libgb
export Semiring
function _semiringnames(name)
    if isGxB(name) || isGrB(name)
        simple = name[5:end]
    else
        simple = name
    end
    simple = replace(simple, "_SEMIRING" => "")
    containername = Symbol(simple, "_T")
    exportedname = Symbol(simple)
    return containername, exportedname
end

function Semiring(name)
    containername, exportedname = _semiringnames(name)
    structquote = quote
        struct $containername <: AbstractSemiring
            typedops::Dict{DataType, TypedSemiring}
            name::String
            $containername() = new(Dict{DataType, TypedSemiring}(), $name)
        end
    end
    @eval($structquote)
    constquote = quote
        const $exportedname = $containername()
        export $exportedname
    end
    @eval($constquote)
    return getproperty(Semirings, exportedname)
end

end

SemiringUnion = Union{AbstractSemiring, TypedSemiring}

#TODO: Rewrite
function _createsemirings()
    builtins = [
        "GxB_PLUS_FIRST",
        "GxB_TIMES_FIRST",
        "GxB_ANY_FIRST",
        "GxB_PLUS_SECOND",
        "GxB_TIMES_SECOND",
        "GxB_ANY_SECOND",
        "GxB_MIN_PAIR",
        "GxB_MAX_PAIR",
        "GxB_PLUS_PAIR",
        "GxB_TIMES_PAIR",
        "GxB_ANY_PAIR",
        "GxB_MIN_MIN",
        "GxB_TIMES_MIN",
        "GxB_ANY_MIN",
        "GxB_MAX_MAX",
        "GxB_PLUS_MAX",
        "GxB_TIMES_MAX",
        "GxB_ANY_MAX",
        "GxB_PLUS_PLUS",
        "GxB_TIMES_PLUS",
        "GxB_ANY_PLUS",
        "GxB_MIN_MINUS",
        "GxB_MAX_MINUS",
        "GxB_PLUS_MINUS",
        "GxB_TIMES_MINUS",
        "GxB_ANY_MINUS",
        "GxB_TIMES_TIMES",
        "GxB_ANY_TIMES",
        "GxB_MIN_DIV",
        "GxB_MAX_DIV",
        "GxB_PLUS_DIV",
        "GxB_TIMES_DIV",
        "GxB_ANY_DIV",
        "GxB_MIN_RDIV",
        "GxB_MAX_RDIV",
        "GxB_PLUS_RDIV",
        "GxB_TIMES_RDIV",
        "GxB_ANY_RDIV",
        "GxB_MIN_RMINUS",
        "GxB_MAX_RMINUS",
        "GxB_PLUS_RMINUS",
        "GxB_TIMES_RMINUS",
        "GxB_ANY_RMINUS",
        "GxB_MIN_ISEQ",
        "GxB_MAX_ISEQ",
        "GxB_PLUS_ISEQ",
        "GxB_TIMES_ISEQ",
        "GxB_ANY_ISEQ",
        "GxB_MIN_ISNE",
        "GxB_MAX_ISNE",
        "GxB_PLUS_ISNE",
        "GxB_TIMES_ISNE",
        "GxB_ANY_ISNE",
        "GxB_MIN_ISGT",
        "GxB_MAX_ISGT",
        "GxB_PLUS_ISGT",
        "GxB_TIMES_ISGT",
        "GxB_ANY_ISGT",
        "GxB_MIN_ISLT",
        "GxB_MAX_ISLT",
        "GxB_PLUS_ISLT",
        "GxB_TIMES_ISLT",
        "GxB_ANY_ISLT",
        "GxB_MIN_ISGE",
        "GxB_MAX_ISGE",
        "GxB_PLUS_ISGE",
        "GxB_TIMES_ISGE",
        "GxB_ANY_ISGE",
        "GxB_MIN_ISLE",
        "GxB_MAX_ISLE",
        "GxB_PLUS_ISLE",
        "GxB_TIMES_ISLE",
        "GxB_ANY_ISLE",
        "GxB_MIN_LOR",
        "GxB_MAX_LOR",
        "GxB_PLUS_LOR",
        "GxB_TIMES_LOR",
        "GxB_ANY_LOR",
        "GxB_MIN_LAND",
        "GxB_MAX_LAND",
        "GxB_PLUS_LAND",
        "GxB_TIMES_LAND",
        "GxB_ANY_LAND",
        "GxB_MIN_LXOR",
        "GxB_MAX_LXOR",
        "GxB_PLUS_LXOR",
        "GxB_TIMES_LXOR",
        "GxB_ANY_LXOR",
        "GxB_LOR_NE",
        "GxB_LOR_EQ",
        "GxB_LAND_EQ",
        "GxB_LXOR_EQ",
        "GxB_EQ_EQ",
        "GxB_ANY_EQ",
        "GxB_LAND_NE",
        "GxB_LXOR_NE",
        "GxB_EQ_NE",
        "GxB_ANY_NE",
        "GxB_LOR_GT",
        "GxB_LAND_GT",
        "GxB_LXOR_GT",
        "GxB_EQ_GT",
        "GxB_ANY_GT",
        "GxB_LOR_LT",
        "GxB_LAND_LT",
        "GxB_LXOR_LT",
        "GxB_EQ_LT",
        "GxB_ANY_LT",
        "GxB_LOR_GE",
        "GxB_LAND_GE",
        "GxB_LXOR_GE",
        "GxB_EQ_GE",
        "GxB_ANY_GE",
        "GxB_LOR_LE",
        "GxB_LAND_LE",
        "GxB_LXOR_LE",
        "GxB_EQ_LE",
        "GxB_ANY_LE",
        "GxB_LOR_FIRST",
        "GxB_LAND_FIRST",
        "GxB_LXOR_FIRST",
        "GxB_EQ_FIRST",
        "GxB_LOR_SECOND",
        "GxB_LAND_SECOND",
        "GxB_LXOR_SECOND",
        "GxB_EQ_SECOND",
        "GxB_LOR_PAIR",
        "GxB_LAND_PAIR",
        "GxB_LXOR_PAIR",
        "GxB_EQ_PAIR",
        "GxB_LOR_LOR",
        "GxB_LXOR_LOR",
        "GxB_EQ_LOR",
        "GxB_LAND_LAND",
        "GxB_EQ_LAND",
        "GxB_LOR_LXOR",
        "GxB_LAND_LXOR",
        "GxB_LXOR_LXOR",
        "GxB_EQ_LXOR",
        "GxB_BOR_BOR",
        "GxB_BOR_BAND",
        "GxB_BOR_BXOR",
        "GxB_BOR_BXNOR",
        "GxB_BAND_BOR",
        "GxB_BAND_BAND",
        "GxB_BAND_BXOR",
        "GxB_BAND_BXNOR",
        "GxB_BXOR_BOR",
        "GxB_BXOR_BAND",
        "GxB_BXOR_BXOR",
        "GxB_BXOR_BXNOR",
        "GxB_BXNOR_BOR",
        "GxB_BXNOR_BAND",
        "GxB_BXNOR_BXOR",
        "GxB_BXNOR_BXNOR",
        "GxB_MIN_FIRSTI",
        "GxB_MAX_FIRSTI",
        "GxB_ANY_FIRSTI",
        "GxB_PLUS_FIRSTI",
        "GxB_TIMES_FIRSTI",
        "GxB_MIN_FIRSTI1",
        "GxB_MAX_FIRSTI1",
        "GxB_ANY_FIRSTI1",
        "GxB_PLUS_FIRSTI1",
        "GxB_TIMES_FIRSTI1",
        "GxB_MIN_FIRSTJ",
        "GxB_MAX_FIRSTJ",
        "GxB_ANY_FIRSTJ",
        "GxB_PLUS_FIRSTJ",
        "GxB_TIMES_FIRSTJ",
        "GxB_MIN_FIRSTJ1",
        "GxB_MAX_FIRSTJ1",
        "GxB_ANY_FIRSTJ1",
        "GxB_PLUS_FIRSTJ1",
        "GxB_TIMES_FIRSTJ1",
        "GxB_MIN_SECONDI",
        "GxB_MAX_SECONDI",
        "GxB_ANY_SECONDI",
        "GxB_PLUS_SECONDI",
        "GxB_TIMES_SECONDI",
        "GxB_MIN_SECONDI1",
        "GxB_MAX_SECONDI1",
        "GxB_ANY_SECONDI1",
        "GxB_PLUS_SECONDI1",
        "GxB_TIMES_SECONDI1",
        "GxB_MIN_SECONDJ",
        "GxB_MAX_SECONDJ",
        "GxB_ANY_SECONDJ",
        "GxB_PLUS_SECONDJ",
        "GxB_TIMES_SECONDJ",
        "GxB_MIN_SECONDJ1",
        "GxB_MAX_SECONDJ1",
        "GxB_ANY_SECONDJ1",
        "GxB_PLUS_SECONDJ1",
        "GxB_TIMES_SECONDJ1",
        "GrB_PLUS_TIMES_SEMIRING",
        "GrB_PLUS_MIN_SEMIRING",
        "GrB_MIN_PLUS_SEMIRING",
        "GrB_MIN_TIMES_SEMIRING",
        "GrB_MIN_FIRST_SEMIRING",
        "GrB_MIN_SECOND_SEMIRING",
        "GrB_MIN_MAX_SEMIRING",
        "GrB_MAX_PLUS_SEMIRING",
        "GrB_MAX_TIMES_SEMIRING",
        "GrB_MAX_FIRST_SEMIRING",
        "GrB_MAX_SECOND_SEMIRING",
        "GrB_MAX_MIN_SEMIRING",
        "GrB_LOR_LAND_SEMIRING",
        "GrB_LAND_LOR_SEMIRING",
        "GrB_LXOR_LAND_SEMIRING",
        "GrB_LXNOR_LOR_SEMIRING",
    ]

    for name ∈ builtins
        Semiring(name)
    end
end

#Add typed ⊕ and ⊗ to semiring
function _addsemiring(rig::AbstractSemiring, add::TypedMonoid, mul::TypedBinaryOperator)
    rigref = Ref{TypedSemiring}()
    libgb.GrB_Semiring_new(rigref, add, mul)
    rig.typedops[xtype(add)] = TypedSemiring(rigref[])
    return nothing
end

#New semiring with typed ⊕ and ⊗
function Semirings.Semiring(name::String, add::TypedMonoid, mul::TypedBinaryOperator)
    rig = Semirings.Semiring(name)
    _addsemiring(rig, add, mul)
    return rig
end

#New semiring with all supported types in the intersection of validtypes(⊕) and validtypes(⊗)
function Semirings.Semiring(name::String, add::AbstractMonoid, mul::AbstractBinaryOp)
    tadd = validtypes(add)
    tmul = validtypes(mul)
    trig = intersect(tadd, tmul)
    rig = Semirings.Semiring(name)
    for t ∈ trig
        _addsemiring(rig, add[t], mul[t])
    end
    return rig
end

ztype(::TypedSemiring{X, Y, Z}) where {X, Y, Z} = Z
xtype(::TypedSemiring{X, Y, Z}) where {X, Y, Z} = X
ytype(::TypedSemiring{X, Y, Z}) where {X, Y, Z} = Y

function _load(rig::AbstractSemiring)
    booleans = ["GxB_LOR_FIRST",
        "GxB_LAND_FIRST",
        "GxB_LXOR_FIRST",
        "GxB_EQ_FIRST",
        "GxB_ANY_FIRST",
        "GxB_LOR_SECOND",
        "GxB_LAND_SECOND",
        "GxB_LXOR_SECOND",
        "GxB_EQ_SECOND",
        "GxB_ANY_SECOND",
        "GxB_LOR_PAIR",
        "GxB_LAND_PAIR",
        "GxB_LXOR_PAIR",
        "GxB_EQ_PAIR",
        "GxB_ANY_PAIR",
        "GxB_LOR_LOR",
        "GxB_LXOR_LOR",
        "GxB_EQ_LOR",
        "GxB_ANY_LOR",
        "GxB_LAND_LAND",
        "GxB_EQ_LAND",
        "GxB_ANY_LAND",
        "GxB_LOR_LXOR",
        "GxB_LAND_LXOR",
        "GxB_LXOR_LXOR",
        "GxB_EQ_LXOR",
        "GxB_ANY_LXOR",
        "GxB_LOR_EQ",
        "GxB_LAND_EQ",
        "GxB_LXOR_EQ",
        "GxB_EQ_EQ",
        "GxB_ANY_EQ",
        "GxB_LOR_GT",
        "GxB_LAND_GT",
        "GxB_LXOR_GT",
        "GxB_EQ_GT",
        "GxB_ANY_GT",
        "GxB_LOR_LT",
        "GxB_LAND_LT",
        "GxB_LXOR_LT",
        "GxB_EQ_LT",
        "GxB_ANY_LT",
        "GxB_LOR_GE",
        "GxB_LAND_GE",
        "GxB_LXOR_GE",
        "GxB_EQ_GE",
        "GxB_ANY_GE",
        "GxB_LOR_LE",
        "GxB_LAND_LE",
        "GxB_LXOR_LE",
        "GxB_EQ_LE",
        "GxB_ANY_LE",
        "GrB_LOR_LAND_SEMIRING",
        "GrB_LAND_LOR_SEMIRING",
        "GrB_LXOR_LAND_SEMIRING",
        "GrB_LXNOR_LOR_SEMIRING",
    ]

    integers = ["GxB_PLUS_FIRST",
        "GxB_TIMES_FIRST",
        "GxB_ANY_FIRST",
        "GxB_PLUS_SECOND",
        "GxB_TIMES_SECOND",
        "GxB_ANY_SECOND",
        "GxB_MIN_PAIR",
        "GxB_MAX_PAIR",
        "GxB_PLUS_PAIR",
        "GxB_TIMES_PAIR",
        "GxB_ANY_PAIR",
        "GxB_MIN_MIN",
        "GxB_TIMES_MIN",
        "GxB_ANY_MIN",
        "GxB_MAX_MAX",
        "GxB_PLUS_MAX",
        "GxB_TIMES_MAX",
        "GxB_ANY_MAX",
        "GxB_PLUS_PLUS",
        "GxB_TIMES_PLUS",
        "GxB_ANY_PLUS",
        "GxB_MIN_MINUS",
        "GxB_MAX_MINUS",
        "GxB_PLUS_MINUS",
        "GxB_TIMES_MINUS",
        "GxB_ANY_MINUS",
        "GxB_TIMES_TIMES",
        "GxB_ANY_TIMES",
        "GxB_MIN_DIV",
        "GxB_MAX_DIV",
        "GxB_PLUS_DIV",
        "GxB_TIMES_DIV",
        "GxB_ANY_DIV",
        "GxB_MIN_RDIV",
        "GxB_MAX_RDIV",
        "GxB_PLUS_RDIV",
        "GxB_TIMES_RDIV",
        "GxB_ANY_RDIV",
        "GxB_MIN_RMINUS",
        "GxB_MAX_RMINUS",
        "GxB_PLUS_RMINUS",
        "GxB_TIMES_RMINUS",
        "GxB_ANY_RMINUS",
        "GxB_MIN_ISEQ",
        "GxB_MAX_ISEQ",
        "GxB_PLUS_ISEQ",
        "GxB_TIMES_ISEQ",
        "GxB_ANY_ISEQ",
        "GxB_MIN_ISNE",
        "GxB_MAX_ISNE",
        "GxB_PLUS_ISNE",
        "GxB_TIMES_ISNE",
        "GxB_ANY_ISNE",
        "GxB_MIN_ISGT",
        "GxB_MAX_ISGT",
        "GxB_PLUS_ISGT",
        "GxB_TIMES_ISGT",
        "GxB_ANY_ISGT",
        "GxB_MIN_ISLT",
        "GxB_MAX_ISLT",
        "GxB_PLUS_ISLT",
        "GxB_TIMES_ISLT",
        "GxB_ANY_ISLT",
        "GxB_MIN_ISGE",
        "GxB_MAX_ISGE",
        "GxB_PLUS_ISGE",
        "GxB_TIMES_ISGE",
        "GxB_ANY_ISGE",
        "GxB_MIN_ISLE",
        "GxB_MAX_ISLE",
        "GxB_PLUS_ISLE",
        "GxB_TIMES_ISLE",
        "GxB_ANY_ISLE",
        "GxB_MIN_LOR",
        "GxB_MAX_LOR",
        "GxB_PLUS_LOR",
        "GxB_TIMES_LOR",
        "GxB_ANY_LOR",
        "GxB_MIN_LAND",
        "GxB_MAX_LAND",
        "GxB_PLUS_LAND",
        "GxB_TIMES_LAND",
        "GxB_ANY_LAND",
        "GxB_MIN_LXOR",
        "GxB_MAX_LXOR",
        "GxB_PLUS_LXOR",
        "GxB_TIMES_LXOR",
        "GxB_ANY_LXOR",
        "GxB_LOR_NE",
        "GxB_LOR_EQ",
        "GxB_LAND_EQ",
        "GxB_LXOR_EQ",
        "GxB_EQ_EQ",
        "GxB_ANY_EQ",
        "GxB_LAND_NE",
        "GxB_LXOR_NE",
        "GxB_EQ_NE",
        "GxB_ANY_NE",
        "GxB_LOR_GT",
        "GxB_LAND_GT",
        "GxB_LXOR_GT",
        "GxB_EQ_GT",
        "GxB_ANY_GT",
        "GxB_LOR_LT",
        "GxB_LAND_LT",
        "GxB_LXOR_LT",
        "GxB_EQ_LT",
        "GxB_ANY_LT",
        "GxB_LOR_GE",
        "GxB_LAND_GE",
        "GxB_LXOR_GE",
        "GxB_EQ_GE",
        "GxB_ANY_GE",
        "GxB_LOR_LE",
        "GxB_LAND_LE",
        "GxB_LXOR_LE",
        "GxB_EQ_LE",
        "GxB_ANY_LE",
        "GrB_PLUS_TIMES_SEMIRING",
        "GrB_PLUS_MIN_SEMIRING",
        "GrB_MIN_PLUS_SEMIRING",
        "GrB_MIN_TIMES_SEMIRING",
        "GrB_MIN_FIRST_SEMIRING",
        "GrB_MIN_SECOND_SEMIRING",
        "GrB_MIN_MAX_SEMIRING",
        "GrB_MAX_PLUS_SEMIRING",
        "GrB_MAX_TIMES_SEMIRING",
        "GrB_MAX_FIRST_SEMIRING",
        "GrB_MAX_SECOND_SEMIRING",
        "GrB_MAX_MIN_SEMIRING",
    ]

    unsignedintegers = ["GxB_PLUS_FIRST",
        "GxB_TIMES_FIRST",
        "GxB_ANY_FIRST",
        "GxB_PLUS_SECOND",
        "GxB_TIMES_SECOND",
        "GxB_ANY_SECOND",
        "GxB_MIN_PAIR",
        "GxB_MAX_PAIR",
        "GxB_PLUS_PAIR",
        "GxB_TIMES_PAIR",
        "GxB_ANY_PAIR",
        "GxB_MIN_MIN",
        "GxB_TIMES_MIN",
        "GxB_ANY_MIN",
        "GxB_MAX_MAX",
        "GxB_PLUS_MAX",
        "GxB_TIMES_MAX",
        "GxB_ANY_MAX",
        "GxB_PLUS_PLUS",
        "GxB_TIMES_PLUS",
        "GxB_ANY_PLUS",
        "GxB_MIN_MINUS",
        "GxB_MAX_MINUS",
        "GxB_PLUS_MINUS",
        "GxB_TIMES_MINUS",
        "GxB_ANY_MINUS",
        "GxB_TIMES_TIMES",
        "GxB_ANY_TIMES",
        "GxB_MIN_DIV",
        "GxB_MAX_DIV",
        "GxB_PLUS_DIV",
        "GxB_TIMES_DIV",
        "GxB_ANY_DIV",
        "GxB_MIN_RDIV",
        "GxB_MAX_RDIV",
        "GxB_PLUS_RDIV",
        "GxB_TIMES_RDIV",
        "GxB_ANY_RDIV",
        "GxB_MIN_RMINUS",
        "GxB_MAX_RMINUS",
        "GxB_PLUS_RMINUS",
        "GxB_TIMES_RMINUS",
        "GxB_ANY_RMINUS",
        "GxB_MIN_ISEQ",
        "GxB_MAX_ISEQ",
        "GxB_PLUS_ISEQ",
        "GxB_TIMES_ISEQ",
        "GxB_ANY_ISEQ",
        "GxB_MIN_ISNE",
        "GxB_MAX_ISNE",
        "GxB_PLUS_ISNE",
        "GxB_TIMES_ISNE",
        "GxB_ANY_ISNE",
        "GxB_MIN_ISGT",
        "GxB_MAX_ISGT",
        "GxB_PLUS_ISGT",
        "GxB_TIMES_ISGT",
        "GxB_ANY_ISGT",
        "GxB_MIN_ISLT",
        "GxB_MAX_ISLT",
        "GxB_PLUS_ISLT",
        "GxB_TIMES_ISLT",
        "GxB_ANY_ISLT",
        "GxB_MIN_ISGE",
        "GxB_MAX_ISGE",
        "GxB_PLUS_ISGE",
        "GxB_TIMES_ISGE",
        "GxB_ANY_ISGE",
        "GxB_MIN_ISLE",
        "GxB_MAX_ISLE",
        "GxB_PLUS_ISLE",
        "GxB_TIMES_ISLE",
        "GxB_ANY_ISLE",
        "GxB_MIN_LOR",
        "GxB_MAX_LOR",
        "GxB_PLUS_LOR",
        "GxB_TIMES_LOR",
        "GxB_ANY_LOR",
        "GxB_MIN_LAND",
        "GxB_MAX_LAND",
        "GxB_PLUS_LAND",
        "GxB_TIMES_LAND",
        "GxB_ANY_LAND",
        "GxB_MIN_LXOR",
        "GxB_MAX_LXOR",
        "GxB_PLUS_LXOR",
        "GxB_TIMES_LXOR",
        "GxB_ANY_LXOR",
        "GxB_LOR_NE",
        "GxB_LOR_EQ",
        "GxB_LAND_EQ",
        "GxB_LXOR_EQ",
        "GxB_EQ_EQ",
        "GxB_ANY_EQ",
        "GxB_LAND_NE",
        "GxB_LXOR_NE",
        "GxB_EQ_NE",
        "GxB_ANY_NE",
        "GxB_LOR_GT",
        "GxB_LAND_GT",
        "GxB_LXOR_GT",
        "GxB_EQ_GT",
        "GxB_ANY_GT",
        "GxB_LOR_LT",
        "GxB_LAND_LT",
        "GxB_LXOR_LT",
        "GxB_EQ_LT",
        "GxB_ANY_LT",
        "GxB_LOR_GE",
        "GxB_LAND_GE",
        "GxB_LXOR_GE",
        "GxB_EQ_GE",
        "GxB_ANY_GE",
        "GxB_LOR_LE",
        "GxB_LAND_LE",
        "GxB_LXOR_LE",
        "GxB_EQ_LE",
        "GxB_ANY_LE",
        "GxB_BOR_BOR",
        "GxB_BOR_BAND",
        "GxB_BOR_BXOR",
        "GxB_BOR_BXNOR",
        "GxB_BAND_BOR",
        "GxB_BAND_BAND",
        "GxB_BAND_BXOR",
        "GxB_BAND_BXNOR",
        "GxB_BXOR_BOR",
        "GxB_BXOR_BAND",
        "GxB_BXOR_BXOR",
        "GxB_BXOR_BXNOR",
        "GxB_BXNOR_BOR",
        "GxB_BXNOR_BAND",
        "GxB_BXNOR_BXOR",
        "GxB_BXNOR_BXNOR",
        "GrB_PLUS_TIMES_SEMIRING",
        "GrB_PLUS_MIN_SEMIRING",
        "GrB_MIN_PLUS_SEMIRING",
        "GrB_MIN_TIMES_SEMIRING",
        "GrB_MIN_FIRST_SEMIRING",
        "GrB_MIN_SECOND_SEMIRING",
        "GrB_MIN_MAX_SEMIRING",
        "GrB_MAX_PLUS_SEMIRING",
        "GrB_MAX_TIMES_SEMIRING",
        "GrB_MAX_FIRST_SEMIRING",
        "GrB_MAX_SECOND_SEMIRING",
        "GrB_MAX_MIN_SEMIRING",
    ]

    floats = ["GxB_PLUS_FIRST",
        "GxB_TIMES_FIRST",
        "GxB_ANY_FIRST",
        "GxB_PLUS_SECOND",
        "GxB_TIMES_SECOND",
        "GxB_ANY_SECOND",
        "GxB_MIN_PAIR",
        "GxB_MAX_PAIR",
        "GxB_PLUS_PAIR",
        "GxB_TIMES_PAIR",
        "GxB_ANY_PAIR",
        "GxB_MIN_MIN",
        "GxB_TIMES_MIN",
        "GxB_ANY_MIN",
        "GxB_MAX_MAX",
        "GxB_PLUS_MAX",
        "GxB_TIMES_MAX",
        "GxB_ANY_MAX",
        "GxB_PLUS_PLUS",
        "GxB_TIMES_PLUS",
        "GxB_ANY_PLUS",
        "GxB_MIN_MINUS",
        "GxB_MAX_MINUS",
        "GxB_PLUS_MINUS",
        "GxB_TIMES_MINUS",
        "GxB_ANY_MINUS",
        "GxB_TIMES_TIMES",
        "GxB_ANY_TIMES",
        "GxB_MIN_DIV",
        "GxB_MAX_DIV",
        "GxB_PLUS_DIV",
        "GxB_TIMES_DIV",
        "GxB_ANY_DIV",
        "GxB_MIN_RDIV",
        "GxB_MAX_RDIV",
        "GxB_PLUS_RDIV",
        "GxB_TIMES_RDIV",
        "GxB_ANY_RDIV",
        "GxB_MIN_RMINUS",
        "GxB_MAX_RMINUS",
        "GxB_PLUS_RMINUS",
        "GxB_TIMES_RMINUS",
        "GxB_ANY_RMINUS",
        "GxB_MIN_ISEQ",
        "GxB_MAX_ISEQ",
        "GxB_PLUS_ISEQ",
        "GxB_TIMES_ISEQ",
        "GxB_ANY_ISEQ",
        "GxB_MIN_ISNE",
        "GxB_MAX_ISNE",
        "GxB_PLUS_ISNE",
        "GxB_TIMES_ISNE",
        "GxB_ANY_ISNE",
        "GxB_MIN_ISGT",
        "GxB_MAX_ISGT",
        "GxB_PLUS_ISGT",
        "GxB_TIMES_ISGT",
        "GxB_ANY_ISGT",
        "GxB_MIN_ISLT",
        "GxB_MAX_ISLT",
        "GxB_PLUS_ISLT",
        "GxB_TIMES_ISLT",
        "GxB_ANY_ISLT",
        "GxB_MIN_ISGE",
        "GxB_MAX_ISGE",
        "GxB_PLUS_ISGE",
        "GxB_TIMES_ISGE",
        "GxB_ANY_ISGE",
        "GxB_MIN_ISLE",
        "GxB_MAX_ISLE",
        "GxB_PLUS_ISLE",
        "GxB_TIMES_ISLE",
        "GxB_ANY_ISLE",
        "GxB_MIN_LOR",
        "GxB_MAX_LOR",
        "GxB_PLUS_LOR",
        "GxB_TIMES_LOR",
        "GxB_ANY_LOR",
        "GxB_MIN_LAND",
        "GxB_MAX_LAND",
        "GxB_PLUS_LAND",
        "GxB_TIMES_LAND",
        "GxB_ANY_LAND",
        "GxB_MIN_LXOR",
        "GxB_MAX_LXOR",
        "GxB_PLUS_LXOR",
        "GxB_TIMES_LXOR",
        "GxB_ANY_LXOR",
        "GxB_LOR_NE",
        "GxB_LOR_EQ",
        "GxB_LAND_EQ",
        "GxB_LXOR_EQ",
        "GxB_EQ_EQ",
        "GxB_ANY_EQ",
        "GxB_LAND_NE",
        "GxB_LXOR_NE",
        "GxB_EQ_NE",
        "GxB_ANY_NE",
        "GxB_LOR_GT",
        "GxB_LAND_GT",
        "GxB_LXOR_GT",
        "GxB_EQ_GT",
        "GxB_ANY_GT",
        "GxB_LOR_LT",
        "GxB_LAND_LT",
        "GxB_LXOR_LT",
        "GxB_EQ_LT",
        "GxB_ANY_LT",
        "GxB_LOR_GE",
        "GxB_LAND_GE",
        "GxB_LXOR_GE",
        "GxB_EQ_GE",
        "GxB_ANY_GE",
        "GxB_LOR_LE",
        "GxB_LAND_LE",
        "GxB_LXOR_LE",
        "GxB_EQ_LE",
        "GxB_ANY_LE",
        "GrB_PLUS_TIMES_SEMIRING",
        "GrB_PLUS_MIN_SEMIRING",
        "GrB_MIN_PLUS_SEMIRING",
        "GrB_MIN_TIMES_SEMIRING",
        "GrB_MIN_FIRST_SEMIRING",
        "GrB_MIN_SECOND_SEMIRING",
        "GrB_MIN_MAX_SEMIRING",
        "GrB_MAX_PLUS_SEMIRING",
        "GrB_MAX_TIMES_SEMIRING",
        "GrB_MAX_FIRST_SEMIRING",
        "GrB_MAX_SECOND_SEMIRING",
        "GrB_MAX_MIN_SEMIRING",
    ]

    positionals = [
        "GxB_MIN_FIRSTI",
        "GxB_MAX_FIRSTI",
        "GxB_ANY_FIRSTI",
        "GxB_PLUS_FIRSTI",
        "GxB_TIMES_FIRSTI",
        "GxB_MIN_FIRSTI1",
        "GxB_MAX_FIRSTI1",
        "GxB_ANY_FIRSTI1",
        "GxB_PLUS_FIRSTI1",
        "GxB_TIMES_FIRSTI1",
        "GxB_MIN_FIRSTJ",
        "GxB_MAX_FIRSTJ",
        "GxB_ANY_FIRSTJ",
        "GxB_PLUS_FIRSTJ",
        "GxB_TIMES_FIRSTJ",
        "GxB_MIN_FIRSTJ1",
        "GxB_MAX_FIRSTJ1",
        "GxB_ANY_FIRSTJ1",
        "GxB_PLUS_FIRSTJ1",
        "GxB_TIMES_FIRSTJ1",
        "GxB_MIN_SECONDI",
        "GxB_MAX_SECONDI",
        "GxB_ANY_SECONDI",
        "GxB_PLUS_SECONDI",
        "GxB_TIMES_SECONDI",
        "GxB_MIN_SECONDI1",
        "GxB_MAX_SECONDI1",
        "GxB_ANY_SECONDI1",
        "GxB_PLUS_SECONDI1",
        "GxB_TIMES_SECONDI1",
        "GxB_MIN_SECONDJ",
        "GxB_MAX_SECONDJ",
        "GxB_ANY_SECONDJ",
        "GxB_PLUS_SECONDJ",
        "GxB_TIMES_SECONDJ",
        "GxB_MIN_SECONDJ1",
        "GxB_MAX_SECONDJ1",
        "GxB_ANY_SECONDJ1",
        "GxB_PLUS_SECONDJ1",
        "GxB_TIMES_SECONDJ1",
    ]

    complexes = [
        "GxB_PLUS_FIRST",
        "GxB_TIMES_FIRST",
        "GxB_ANY_FIRST",
        "GxB_PLUS_SECOND",
        "GxB_TIMES_SECOND",
        "GxB_ANY_SECOND",
        "GxB_PLUS_PAIR",
        "GxB_TIMES_PAIR",
        "GxB_ANY_PAIR",
        "GxB_PLUS_PLUS",
        "GxB_TIMES_PLUS",
        "GxB_ANY_PLUS",
        "GxB_PLUS_MINUS",
        "GxB_TIMES_MINUS",
        "GxB_ANY_MINUS",
        "GxB_PLUS_TIMES",
        "GxB_TIMES_TIMES",
        "GxB_ANY_TIMES",
        "GxB_PLUS_DIV",
        "GxB_TIMES_DIV",
        "GxB_ANY_DIV",
        "GxB_PLUS_RDIV",
        "GxB_TIMES_RDIV",
        "GxB_ANY_RDIV",
        "GxB_PLUS_RMINUS",
        "GxB_TIMES_RMINUS",
        "GxB_ANY_RMINUS",
    ]
    name = rig.name
    if name ∈ booleans
        rig.typedops[Bool] = TypedSemiring(load_global(name * "_BOOL", libgb.GrB_Semiring))
    end

    if name ∈ integers
        rig.typedops[Int8] = TypedSemiring(load_global(name * "_INT8", libgb.GrB_Semiring))
        rig.typedops[Int16] = TypedSemiring(load_global(name * "_INT16", libgb.GrB_Semiring))
        rig.typedops[Int32] = TypedSemiring(load_global(name * "_INT32", libgb.GrB_Semiring))
        rig.typedops[Int64] = TypedSemiring(load_global(name * "_INT64", libgb.GrB_Semiring))
    end

    if name ∈ unsignedintegers
        rig.typedops[UInt8] = TypedSemiring(load_global(name * "_UINT8", libgb.GrB_Semiring))
        rig.typedops[UInt16] = TypedSemiring(load_global(name * "_UINT16", libgb.GrB_Semiring))
        rig.typedops[UInt32] = TypedSemiring(load_global(name * "_UINT32", libgb.GrB_Semiring))
        rig.typedops[UInt64] = TypedSemiring(load_global(name * "_UINT64", libgb.GrB_Semiring))
    end

    if name ∈ floats
        rig.typedops[Float32] = TypedSemiring(load_global(name * "_FP32", libgb.GrB_Semiring))
        rig.typedops[Float64] = TypedSemiring(load_global(name * "_FP64", libgb.GrB_Semiring))
    end
    if name ∈ positionals
        rig.typedops[Any] = TypedSemiring(load_global(name * "_INT64", libgb.GrB_Semiring))
    end
    name = replace(name, "GrB_" => "GxB_")
    name = replace(name, "_SEMIRING" => "")
    if name ∈ complexes
        rig.typedops[ComplexF32] = TypedSemiring(load_global(name * "_FC32", libgb.GrB_Semiring))
        rig.typedops[ComplexF64] = TypedSemiring(load_global(name * "_FC64", libgb.GrB_Semiring))
    end
end


text:
@testset "AWSBatchNodeManager" begin
    @testset "AWS_BATCH_JOB_ID DNE" begin
        withenv(
            "AWS_BATCH_JOB_MAIN_NODE_INDEX" => 1
        ) do
            @test_throws ErrorException AWSBatchNodeManager()
        end
    end

    @testset "AWS_BATCH_JOB_MAIN_NODE_INDEX DNE" begin
        withenv(
            "AWS_BATCH_JOB_ID" => "job_id"
        ) do
            @test_throws ErrorException AWSBatchNodeManager()
        end
    end

    @testset "AWS_BATCH_JOB_NUM_NODES DNE" begin
        withenv(
            "AWS_BATCH_JOB_ID" => "job_id",
            "AWS_BATCH_JOB_NODE_INDEX" => 1,
            "AWS_BATCH_JOB_MAIN_NODE_INDEX" => 1,
        ) do
            @test_throws KeyError AWSBatchNodeManager()
        end
    end

    @testset "AWS_BATCH_JOB_NODE_INDEX != AWS_BATCH_JOB_MAIN_NODE_INDEX" begin
        withenv(
            "AWS_BATCH_JOB_ID" => "job_id",
            "AWS_BATCH_JOB_NODE_INDEX" => 0,
            "AWS_BATCH_JOB_MAIN_NODE_INDEX" => 1
        ) do
            @test_throws ErrorException AWSBatchNodeManager()
        end
    end

    @testset "AWS_BATCH_JOB_NUM_NODES -- String" begin
        withenv(
            "AWS_BATCH_JOB_ID" => "job_id",
            "AWS_BATCH_JOB_MAIN_NODE_INDEX" => 1,
            "AWS_BATCH_JOB_NODE_INDEX" => 1,
            "AWS_BATCH_JOB_NUM_NODES" => "foobar"
        ) do
            @test_throws ArgumentError AWSBatchNodeManager()
        end
    end

    @testset "AWS_BATCH_JOB_NUM_NODES -- Int" begin
        expected_workers = 9

        withenv(
            "AWS_BATCH_JOB_ID" => "job_id",
            "AWS_BATCH_JOB_MAIN_NODE_INDEX" => 1,
            "AWS_BATCH_JOB_NODE_INDEX" => 1,
            "AWS_BATCH_JOB_NUM_NODES" => expected_workers + 1  # Add one to account for the manager
        ) do
            result = AWSBatchNodeManager()

            @test result.num_workers == expected_workers
            @test result.timeout == AWSClusterManagers.AWS_BATCH_NODE_TIMEOUT
        end
    end

    @testset "AWS_BATCH_JOB_NUM_NODES -- Int as a String" begin
        expected_workers = 9

        withenv(
            "AWS_BATCH_JOB_ID" => "job_id",
            "AWS_BATCH_JOB_MAIN_NODE_INDEX" => 1,
            "AWS_BATCH_JOB_NODE_INDEX" => 1,
            "AWS_BATCH_JOB_NUM_NODES" => string(expected_workers + 1)  # Add one to account for the manager
        ) do
            result = AWSBatchNodeManager()

            @test result.num_workers == expected_workers
            @test result.timeout == AWSClusterManagers.AWS_BATCH_NODE_TIMEOUT
        end
    end

    @testset "AWS_BATCH_JOB_NUM_NODES -- Set timeout" begin
        expected_workers = 9
        expected_timeout = Second(5)

        withenv(
            "AWS_BATCH_JOB_ID" => "job_id",
            "AWS_BATCH_JOB_MAIN_NODE_INDEX" => 1,
            "AWS_BATCH_JOB_NODE_INDEX" => 1,
            "AWS_BATCH_JOB_NUM_NODES" => expected_workers + 1  # Add one to account for the manager
        ) do
            result = AWSBatchNodeManager(; timeout=expected_timeout)

            @test result.num_workers == expected_workers
            @test result.timeout == expected_timeout
        end
    end
end

@testset "parse_job_id" begin
    passing_cases = [
        string(repeat("a", 36), "#123"),
        string("0-", repeat("a", 34), "#123")
    ]

    @testset "passing case: $(case)" for case in passing_cases
        str = string("job_id:", case)

        @test AWSClusterManagers.parse_job_id(str) == case
    end

    failing_cases = [
        string(repeat("a", 35), "#123"),  # Fails because too short
        string(repeat("a", 37), "#123"),  # Fails because too long
        string(repeat("!", 36), "#123"),  # Fails because non `a-z0-9-` before #
        string(repeat("a", 36), "#"),  # Fails because no decimal after #
        string(repeat("a", 36), "#abc"),  # Fails because non-decimal after #
    ]

    @testset "failing case: $(case)" for case in failing_cases
        str = string("job_id:", case)

        @test_throws ErrorException AWSClusterManagers.parse_job_id(str)
    end
end

@testset "parse_cookie" begin
    @test AWSClusterManagers.parse_cookie("julia_cookie:1000") == "1000"
    @test AWSClusterManagers.parse_cookie("julia_cookie:1.5") == "1"
    @test AWSClusterManagers.parse_cookie("julia_cookie:foobar") == "foobar"

    @test_throws ErrorException AWSClusterManagers.parse_cookie("julia_cookie:-1")
    @test_throws ErrorException AWSClusterManagers.parse_cookie("foobar")
    @test_throws ErrorException AWSClusterManagers.parse_cookie("foobar:julia_cookie:1")
end

@testset "parse_worker_timeout" begin
    @test AWSClusterManagers.parse_worker_timeout("julia_worker_timeout:1000") == 1000
    @test AWSClusterManagers.parse_worker_timeout("julia_worker_timeout:1.5") == 1

    @test_throws ErrorException AWSClusterManagers.parse_worker_timeout("julia_worker_timeout:-1")
    @test_throws ErrorException AWSClusterManagers.parse_worker_timeout("julia_worker_timeout:foobar")
    @test_throws ErrorException AWSClusterManagers.parse_worker_timeout("foobar")
    @test_throws ErrorException AWSClusterManagers.parse_worker_timeout("foobar:julia_worker_timeout:1")
end


text:
"""
The module implements logic of permutation of input image
"""
module JWebImageDemo

using Images
using Random

export permute_image

"""
Calculate first and last points after centering by specific length
"""
center_begin_end(actual::Int, required::Int) =
    round(Int, (actual - required) / 2) .+ [1, required]

"""
Resize the image
"""
function resize(img::Matrix, maxSize::Int)
    sz = round.(Int, size(img) .* (maxSize / max(size(img)...)))
    σ = map((o, n) -> 0.05 * o / n, size(img), sz)
    kern = KernelFactors.gaussian(σ)   # from ImageFiltering
    return imresize(imfilter(img, kern, NA()), sz)
end

"""
Transform zone number into matrix indices
"""
zone_to_coords(ind::Vector{Int}, step::Int) =
    ind .* step |> list -> map(x -> range(x + 1, stop = x + step), list)

"""
Perform image permutation

# Arguments
- `img::Matrix`: the source image as a matrix of colors
- `num::Int`: the number of parts for splitting each axis
"""
function permute_image(img::Matrix, num = 5::Int)::Matrix
    l = min(size(img)...)
    x = center_begin_end(size(img, 2), l)
    y = center_begin_end(size(img, 1), l)

    data_c = img[first(y):last(y), first(x):last(x)]
    img_round = resize(data_c, round(Int, ceil(l / num) * num))

    img_out = similar(img_round)
    step = round(Int, size(img_round, 1) / num)

    indices = [[i, j] for i = 0:(num-1) for j = 0:(num-1)]
    for (ind, ind_base) in zip(Random.shuffle(indices), indices)
        r_x, r_y = zone_to_coords(ind, step)
        rb_x, rb_y = zone_to_coords(ind_base, step)
        img_out[r_x, r_y] = img_round[rb_x, rb_y]
    end
    return img_out
end

include("server.jl")

end # module