text:
export TR_generic

# One dimensionnal trust region algorithm.
# We approximate our function h by a 2nd degree polynomial (Taylor) q.
# If q is a good approximation of h then we move towards the minimizer of q
# Otherwise we reduce our Trust region interval (Δ) until we have a good
# approximation.
# q(d) = h(t) + h'(t) * d + 0.5 * h''(t) * d²
#
# Specific key word arguments
# eps1 :: Float64. Below that treshold we don't move and reduce Δ.
#                  Treshold to determines if q is a bad approximation of h
# eps2 :: Float64. Over that treshold we augment the size of Δ
#                  Treshold to determines if q is a very good approximation of h
# red :: Float64. Rate at which we reduce Δ if we have a bad approximation
# aug :: Float64. Rate at which we augment Δ if we have a bad approximation
# Δ :: Float64. Size of the trust region at the begining of the algorithm.

function TR_generic(h :: AbstractNLPModel,
                    nlpstop :: AbstractStopping;
                    verbose :: Bool = true,
                    eps1 :: Float64 = 0.2,
                    eps2 :: Float64 = 0.8,
                    red:: Float64 = 0.5,
                    aug :: Float64 = 2.0,
                    Δ :: Float64 = 1.0,
                    direction :: Symbol = :Nwt)
    Δ = [Δ]
    (length(h.meta.x0) > 1) && warn("Not a 1-D problem ")
    t = h.meta.x0; iter = 0; # Establish starting point and iteration counter.
    fₖ = obj(h, t); f₀ = copy(fₖ)
    gₖ = grad(h, t); g₀ = copy(gₖ)

    # H denotes the (approximation of the) second derivative
    if direction == :Sec || direction == :SecA
      H = [1.0]
    else
      H = hess(h, t)
    end

    OK = update_and_start!(nlpstop, x = t, fx = fₖ, gx = gₖ, g0 = g₀, Hx = H)

    verbose && @printf(" iter  t         gₖ          Δ        pred         ared\n")
    verbose && @printf(" %4d %7.2e  %7.2e  %7.2e \n", iter, t[1], gₖ[1], Δ[1])

    # We loop until we have a minimizer or we have reached the maximum number of
    # iterations.
    while !OK
        dN = -gₖ ./ H;                 # (approximation of the) Newton direction
        d = TR_step_computation(H, gₖ, dN, Δ) # determines the right direction
                                              # depending on our trust region

        # Numerical reduction computation
        ftestTR = obj(h, t + d)        # Value of h and h' at t + d
        gtestTR = grad(h, t + d)

        # We check to see if we have a good approximation of h using the ratio
        # of the actual reduction and the predicted reduction.
        (pred, ared, ratio) =
            pred_ared_computation(gₖ, fₖ, H, d, ftestTR, gtestTR)

        if (ratio .< eps1) &&  (true in (abs.(d) .== Δ))
            # Bad approximation of h. Reduction of the size of the trust region
            Δ = red * Δ
            add_stop!(nlpstop.meta)
            # printstyled("on a reduit delta \n", color = :red)
        else
            # Good approximation we move towards the minimizer of q
            (t, fₖ, gₖ, H) =
                step_computation(direction, h, t, d, gₖ, fₖ, ftestTR, gtestTR)

            if ratio > eps2
                # Very good approximation of h. We increase our trust region.
                Δ = aug * Δ
            end
            OK = update_and_stop!(nlpstop, x = t, fx = fₖ, gx = gₖ, Hx = H)
            # @show OK
        end

        iter += 1
        verbose && @printf(" %4d %7.2e  %7.2e  %7.2e %7.2e %7.2e\n",
                           iter, t[1], gₖ[1], Δ[1], pred, ared)
    end

    optimal = nlpstop.meta.optimal
    # @show optimal
    return optimal, nlpstop
end


text:
module PkgFetcher

using Tar
using Downloads
using Downloads.Curl: Multi

using Gzip_jll: gzip_path
using zrl_jll: zrle_path, zrld_path

const gzip = `$gzip_path`
const gzcat = `$gzip_path -c -d`
const zrle = `$zrle_path`
const zrld = `$zrld_path`

include("BSPatch.jl")
using .BSPatch

function download(curl::Multi, url::String)
    path, io = mktemp()
    # TODO: add various Pkg headers
    request = Request(io, url, Pair{String,String}[])
    try
        response = Downloader.get(request, curl) # TODO: hook up progress
        response.status == 200 ||
            error("HTTP $(response.status) response")
    catch error
        @warn "download error" url error
        close(io)
        rm(path, force=true)
        return
    end
    close(io)
    return path
end

skeleton_file(path::AbstractString) = "$path.skel"

struct Resource
    hash::String
    path::String
    # TODO: add artifacts after package download
    # post::Function
end

Resource(hash::AbstractString) = Resource(hash, tempname())

abstract type Source end

const Resources = Dict{Resource,<:AbstractVector{<:Source}}

struct Tarball <: Source
    url::String
end

struct Patch <: Source
    url::String
    old::String
end

struct GitRepo <: Source
    url::String
    get::Bool
end

fetch(resources::Resources) = fetch(Multi(), resources)

function fetch(curl::Multi, resources::Resources)
    @sync for (resource, sources) in resources
        @async for source in sources
            tarball = fetch(curl, source)
            tarball === nothing && continue
            hash = Tar.tree_hash(`$zcat $tarball`)
            if hash != resource.hash
                @warn "hash mismatch" resource source hash
                rm(tarball, force=true)
                continue
            end
            skeleton = skeleton_file(resource.path)
            chmod(resource.path, 0o700, recursive=true)
            rm(resource.path, force=true, recursive=true)
            Tar.extract(`$zcat $tarball`, resource.path; skeleton)
            # mv(tarball, resource.path)
            break # success!
        end
    end
end

fetch(curl::Multi, source::Tarball) = download(curl, source.url)

function fetch(curl::Multi, source::Patch)
    patch = download(curl, source.url)
    patch === nothing && return
    try
        mktemp() do old_tarball, old_tar
            skeleton = skeleton_file(source.old)
            if isfile(skeleton)
                Tar.create(source.old, old_tar; skeleton)
            else
                # no skeleton: assume everything goes in the tarball
                # might fail or give wrong data, but we'll catch it
                Tar.create(source.old, old_tar)
            end
            close(old_tar)
            new_tarball = tempname()
            old_tar_zrl = `$zrle $old_tarball`
            new_tar_zrl = pipeline(zrld, gzip, new_tarball)
            try
                bspatch(old_tar_zrl, new_tar_zrl, patch)
            catch error
                @warn "could not apply patch" error
                rm(new_tarball, force=true)
                return
            end
            new_tarball
        end
    finally
        rm(patch, force=true)
    end
end

const SYSTEM_GIT = Ref{Union{Nothing,Bool}}(nothing)

function system_git()
    SYSTEM_GIT[] isa Bool && return SYSTEM_GIT[]
    SYSTEM_GIT[] = Sys.which("git") !== nothing
end

function fetch(curl::Multi, source::GitRepo)
    if system_git()
    else # use LibGit2 instead
    end
end

end # module

ext: jl


text:
# Arrays

glorot_uniform(dims...) = glorot_uniform(STD_REAL_PREC, dims...)
glorot_normal(dims...) = glorot_normal(STD_REAL_PREC, dims...)
rescaled_normal(scale::Real, dims...) = rescaled_normal(STD_REAL_PREC, scale, dims...)

glorot_uniform(T::Type, dims...) = (rand(T, dims...) .- T(0.5)) .* sqrt(T(24)/sum(dims))
glorot_normal(T::Type, dims...) = randn(T, dims...) .* sqrt(T(2)/sum(dims))
rescaled_normal(T::Type, scale::Real, dims::Integer...) = randn(T, dims...) .* T(scale) .* sqrt(T(24)/sum(dims))


text:
function replace_zero(x, tol = 1e-4)
    return isapprox(x, 0.0, atol = tol) ? 0.0 : x
end

function replace_inf(x, lim = 1e5)
    if x == Inf
        return lim
    elseif x == -Inf
        return -lim
    else
        return x
    end
end

function plot_subsets(sampled_subsets, factor_title = " "; columns = [:D, :DELB])
    upscale    = 2
    small_font = Plots.font("sans-serif", 10.0 * upscale)
    large_font = Plots.font("sans-serif", 14.0 * upscale)
    default(titlefont  = large_font,
            guidefont  = large_font,
            tickfont   = small_font,
            legendfont = small_font)
    default(size = (700 * upscale, 900 * upscale))
    default(dpi = 300)

    plotly()

    subplots = []

    global_max_d = 0.0

    sampled_subsets = deepcopy(sampled_subsets)

    for subset in sampled_subsets
        global_max_d = max(global_max_d, subset[:D]...)
    end

    global_max_x = global_max_d == 0.0 ? string(global_max_d) : string("1e",
                                                                       floor(log(10, global_max_d)))

    for subset in sampled_subsets
        for column in columns
            subset[column] = replace_inf(subset[column])
            subset[column] = replace_zero.(subset[column])
        end

        max_d = max(subset[:D]...)
        max_x = max_d == 0.0 ? string(max_d) : string("1e",
                                                      floor(log(10, max_d)))

        push!(subplots,
              histogram(Array(subset[:D]), labels = "Designs",
                        normalize = :probability,
                        title  = string("D-Optimality for ", subset[:Length][1],
                                        " Experiments, ", size(subset, 1),
                                        " Samples and ", factor_title, " Factors"),
                        xticks = ([0, max_d, global_max_d],
                                  ["0", max_x, global_max_x]),
                        xlims = (0.0, global_max_d),
                        color  = :lightblue),
              histogram(Array(subset[:DELB]), labels = "Designs",
                        normalize = :probability,
                        title = string("D-Efficiency for ", subset[:Length][1],
                                        " Experiments, ", size(subset, 1),
                                        " Samples and ", factor_title, " Factors"),
                        xlims = (0.0, 1.0),
                        color = :lightgreen))
    end

    plot(subplots..., legend = false,
         layout = (length(sampled_subsets), 2))
end


text:
using BinDeps

@BinDeps.setup

libmongoc = library_dependency(
    "libmongoc",
    aliases = ["libmongoc", "libmongoc-1.0"]
    )

provides(Sources, Dict(
    URI("http://github.com/mongodb/mongo-c-driver/releases/download/1.0.0/mongo-c-driver-1.0.0.tar.gz") => libmongoc
    ))

provides(
    BuildProcess,
    Autotools(libtarget = "libmongoc-1.0.la"),
    libmongoc,
    os = :Unix
    )

@osx_only begin
    using Homebrew
    provides(Homebrew.HB, Dict("mongo-c" => libmongoc))
end

@BinDeps.install Dict(:libmongoc => :libmongoc)


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

using FEMBase
using FEMBase: is_field_problem, is_boundary_problem
using FEMBase: get_parent_field_name, get_global_solution
using FEMBase: get_assembly, get_elements
using Test

import FEMBase: get_unknown_field_dimension, get_unknown_field_name
import FEMBase: get_formulation_type, assemble_elements!

mutable struct P1 <: FieldProblem
    A :: Bool
end

function P1()
    return P1(true)
end

mutable struct P2 <: BoundaryProblem
    formulation :: Symbol
end

function P2()
    return P2(:incremental)
end

mutable struct P3 <: FieldProblem
end

get_unknown_field_name(p::P1) = "P1"
get_formulation_type(p::Problem{P2}) = p.properties.formulation

@testset "test get_unknown_field_name" begin
    p = Problem(P3, "P3", 1)
    @test get_unknown_field_name(p) == "N/A"
end

@testset "test creating new problems" begin
    p1 = Problem(P1, "P1", 1)
    p2 = Problem(P2, "P2", 2, "p")
    @test get_formulation_type(p1) == :incremental
    @test get_unknown_field_name(p2) == "lambda"
    @test isempty(get_assembly(p1))
    @test is_field_problem(p1)
    @test !is_field_problem(p2)
    @test !is_boundary_problem(p1)
    @test is_boundary_problem(p2)
    @test p1.properties.A
    update!(p1.properties, "A" => "false")
    @test !p1.properties.A
    el = Element(Seg2, [1, 2])
    add_elements!(p1, [el])
    add_elements!(p2, [el])
    initialize!(p1, el, 0.0)
    initialize!(p2, el, 0.0)
    #@test haskey(el, "P1")
    #@test haskey(el, "P2")

    A = Assembly()
    A2 = Assembly()
    u = ones(2)
    la = ones(2)
    update!(p1, A, u, la)
    update!(p1, A, get_elements(p1), 0.0)
    @test length(A.u) == length(u)
    @test length(A.la) == length(la)
    @test size(get_global_solution(p1, A)[1]) == (2,)

    update!(p2, A2, ones(4), ones(4))
    update!(p2, A2, get_elements(p2), 0.0)
    u, la = get_global_solution(p2, A2)
    o = Vector{Float64}[ones(2), ones(2)]
    @test isapprox(u, o)
    @test isapprox(la, o)
    update!(p2, A2, ones(4), ones(4))
    u, la = get_global_solution(p2, A2)
    @test isapprox(u, 2*o)
    @test isapprox(la, o)
    p2.properties.formulation = :forwarddiff
    update!(p2, A2, ones(4), ones(4))
    u, la = get_global_solution(p2, A2)
    @test isapprox(u, 3*o)
    @test isapprox(la, 2*o)
    p2.properties.formulation = :total
    update!(p2, A2, ones(4), ones(4))
    u, la = get_global_solution(p2, A2)
    @test isapprox(u, o)
    @test isapprox(la, o)
    p2.properties.formulation = :nottell
    @test_throws Exception update!(p2, A2, ones(4), ones(4))

    @test length(p1) == 1
    update!(p1, "geometry", Dict(1 => [0.0, 0.0], 2 => [1.0, 0.0]))
    @test haskey(el, "geometry")
    @debug "p1 geometry field" X = p1("geometry", 0.0)
    @debug "e1 geometry field" X = el("geometry", 0.0)
    @test isapprox(p1("geometry", 0.0)[1], [0.0, 0.0])
    @test get_parent_field_name(p2) == "p"
    @test get_gdofs(p1, el) == [1, 2]

    p3 = Problem(P2, "P3", 1, "p")
    as = get_assembly(p3)
    e1 = Element(Seg2, [1, 2])
    e2 = Element(Seg2, [2, 3])
    e3 = Element(Seg2, Int[])
    update!(e1, "f", [1.0, 1.0])
    update!(e2, "f", [2.0, 2.0])
    push!(p3, e1)
    push!(p3, [e2], [e3])
    initialize!(p3, e1, 0.0)
    initialize!(p3, e2, 0.0)
    p3.fields["f"] = field(1.0)
    @test haskey(p3, "f")
    @test isapprox(p3["f"].data, 1.0)
    f = p3("f", 0.0)
    @test_throws Exception get_gdofs(problem, e3)
end

@testset "test initializing new problems" begin
    p1 = Problem(P1, "P1", 2)
    p2 = Problem(P2, "P2", 2, "P1")
    @test get_unknown_field_name(p1) == "P1"
    @test get_unknown_field_name(p2) == "lambda"
    @test get_parent_field_name(p2) == "P1"
    e1 = Element(Seg2, [1, 2])
    e2 = Element(Seg2, [1, 2])
    add_elements!(p1, [e1])
    add_elements!(p2, [e2])
    initialize!(p1, 0.0)
    initialize!(p2, 0.0)
    @debug "e1 keys" keys(e1.fields)
    @test isapprox(e1("P1", (0.0,), 0.0), [0.0, 0.0])
    @test isapprox(e2("P1", (0.0,), 0.0), [0.0, 0.0])
    @test isapprox(e2("lambda", (0.0,), 0.0), [0.0, 0.0])
end

function assemble_elements!(problem::Problem{P1},
                            assembly::Assembly,
                            elements::Vector{Element{E}},
                            time::Float64) where E

    @debug "assemble_elements!" eltype=E
    bi = BasisInfo(E)
    ndofs = length(E)
    Ke = zeros(ndofs, ndofs)
    K = assembly.K

    for element in elements
        fill!(Ke, 0.0)
        for ip in get_integration_points(element)
            J, detJ, N, dN = element_info!(bi, element, ip, time)
            c = element("coefficient", ip, time)
            Ke += ip.weight * c*dN'*dN * detJ
        end
        gdofs = get_gdofs(problem, element)
        add!(K, gdofs, gdofs, Ke)
    end

    return nothing

end

@testset "test assemble field problem" begin
    el1 = Element(Quad4, [1, 2, 3, 4])
    el2 = Element(Tri3, [3, 2, 5])
    X = Dict(1 => [0.0, 0.0],
             2 => [1.0, 0.0],
             3 => [1.0, 1.0],
             4 => [0.0, 1.0],
             5 => [2.0, 1.0])
    elements = [el1, el2]
    update!(elements, "geometry", X)
    update!(elements, "coefficient", 6.0)

    problem = Problem(P1, "test problem", 1)
    add_elements!(problem, elements)
    time = 0.0
    assemble!(problem, time)
    K_expected = [
                  4.0 -1.0 -2.0 -1.0  0.0
                 -1.0  7.0 -4.0 -2.0  0.0
                 -2.0 -4.0 10.0 -1.0 -3.0
                 -1.0 -2.0 -1.0  4.0  0.0
                  0.0  0.0 -3.0  0.0  3.0]
    @test isapprox(problem.assembly.K, K_expected)
end

mutable struct DirBC <: BoundaryProblem
end

function assemble_elements!(problem::Problem{DirBC},
                            assembly::Assembly,
                            elements::Vector{Element{E}},
                            time::Float64) where E

    name = get_parent_field_name(problem)
    dim = get_unknown_field_dimension(problem)

    data = Dict{Int,Float64}()
    for element in elements
        for i=1:dim
            haskey(element, "$name $dim") || continue
            gdofs = get_gdofs(problem, element)
            ldofs = gdofs[i:dim:end]
            xis = get_reference_element_coordinates(E)
            for (ldof, xi) in zip(ldofs, xis)
                data[ldof] = interpolate(element, "$name $dim", xi, time)
            end
        end
    end

    for (dof, val) in data
        add!(assembly.C1, dof, dof, 1.0)
        add!(assembly.C2, dof, dof, 1.0)
        add!(assembly.g, dof, val)
    end

    return nothing

end

@testset "test assemble boundary problem" begin
    el1 = Element(Seg2, [1, 2])
    el2 = Element(Seg2, [2, 3])
    X = Dict(1 => [0.0, 0.0],
             2 => [1.0, 0.0],
             3 => [1.0, 1.0],
             4 => [0.0, 1.0],
             5 => [2.0, 1.0])
    elements = [el1, el2]
    update!(elements, "geometry", X)
    T0 = Dict(1 => 0.5, 2 => 0.0, 3 => 1.0)
    update!(elements, "temperature 1", T0)
    problem = Problem(DirBC, "boundary", 1, "temperature")
    add_elements!(problem, elements)
    time = 0.0
    assemble!(problem, time)
    C1_expected = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
    g_expected = [0.5, 0.0, 1.0]
    @test isapprox(problem.assembly.C1, C1_expected)
    @test isapprox(problem.assembly.g, g_expected)
end

@testset "test deprecated assembly procedure" begin
    problem = Problem(P3, "P3", 1)
    elements = [Element(Seg2, [1, 2])]
    add_elements!(problem, elements)
    assemble!(problem)
    assemble!(problem.assembly, problem, first(elements), 0.0)
    assemble!(problem.assembly, problem, problem.elements, 0.0)
    assemble_elements!(problem, problem.assembly, elements, 0.0)
end

@testset "test set and get global dofs for element" begin
    element = Element(Seg2, [1, 2])
    problem = Problem(P3, "P3", 2)
    @test get_gdofs(problem, element) == [1, 2, 3, 4]
    set_gdofs!(problem, element, [2, 3, 4, 5])
    @test get_gdofs(problem, element) == [2, 3, 4, 5]
    element2 = Element(Seg2, Int[])
    @test_throws ErrorException get_gdofs(problem, element2)
end


text:
using Visualize
using Visualize: orthographicprojection, perspectiveprojection, lookat, normalmesh, JLCanvas, Area, JLRasterizer, Framebuffer
using FileIO, Images, GeometryTypes
include(Visualize.dir("src", "julia","mesh.jl"))

catmesh = normalmesh(load(Pkg.dir("GLVisualize", "assets", "cat.obj")))
catmesh = Base.view(reinterpret(MeshRender.Vertex, catmesh.parent), catmesh.indexes[1])
proj = perspectiveprojection(42f0, 1f0, 0.1f0, 10f0)
view = lookat(Vec3f0(1.5), Vec3f0(0.5), Vec3f0(0, -1, 0))
solid_color = Vec3f0(0.0, 0.7, 1.0)

light = MeshRender.Light(
    Vec3f0(10),
    Vec3f0(0.0),
    Vec3f0(0.9),
    1.0f0,
    Vec3f0(0.2),
    1.0f0
)
shading = MeshRender.Shading(
    Vec3f0(0.0),
    Vec3f0(0.3),
    8.0f0
)

program_args = (light, shading, (solid_color, proj, view))
canvas = JLCanvas(Area => SimpleRectangle(0, 0, 500, 500))

draw = JLRasterizer(
    catmesh, program_args,
    MeshRender.vertexshader, MeshRender.fragmentshader,
)

@time draw(canvas, catmesh, program_args)

save("test.png", clamp01.(canvas[Framebuffer][1]))


text:
n = 10 # rand(1:10)
@test matrixdepot("tridiag", [1,1,1], [1,1,1,1], [1,1,1]) == matrixdepot("tridiag", 4, 1,1,1)
@test matrixdepot("tridiag", Float64, n) == matrixdepot("tridiag", n)

A = matrixdepot("tridiag", n)
@test isposdef(Array(A))
@test issymmetric(Array(A))
println("'tridiag' passed test...")


text:
#
#=
Examples for developing the JAC tools
=#
println("A)  Tests & development of the electronic structure part.")
println("---------------------------------------------------------")
println("Aa) Test of the radial integration in CI computations (Randolf).")
println("Ab) Read in a full multiplet representation from files compatible to Ratip2012, Grasp2018, ....")
println("Ac) Test of several radial e-e potentials.")
println("Ad) Tests towards a SCF field: B-spline primitives and one-particle spectra in a local potential.")
println("Ae) Test of the CI part for an internally generated neon multiplet without Breit interaction.")
println("Af) Test of the frequency-independent Breit interaction for an internally generated neon multiplet.")
println("Ag) Test of the QED model corrections to the level structure of atoms and ions.")
println("Ah) Generate and normalize continuum orbitals in a local potential.")
println("Ai) Test of the jj-LS transformation of levels from a given multiplet.")
println("An) Test of a mean-field basis and configuration-interaction (CI) expansion.")
println("Ap) Test of restricted-active-space (RAS) expansion.")
println("Aq) Test of the Green(function) expansion.")
println("As) Test of spin-angular coefficients.")
println("At) Test of parallel computating with Julia.")
println("")
println("B)  Tests & development of atomic amplitudes")
println("--------------------------------------------")
println("Ba) Tests of the dipole, em multipole and momentum-transfer amplitudes.")
println("Bb) Tests of the parity non-conservation, Schiff moment and anapole moment amplitudes.")
println("")
println("C)   Tests & development of atomic properties")
println("---------------------------------------------")
println("Ca)  Test of the Einstein module with ASF from an internally generated multiplet.")
println("Cb)  Test of the Hfs module for HFS A,B parameters and hyperfine representation with ASF from an internally generated multiplet.")
println("Cc)  Test of the IsotopeShift module with ASF from an internally generated multiplet.")
println("Cd)  Test of the FormFactor module with ASF from an internally generated multiplet.")
println("Ce)  Test of the PlasmaShift module with ASF from an internally generated multiplet.")
println("Cfg) Test of the LandeZeeman module with ASF from an internally generated multiplet.")
println("Ch)  Test of the AlphaVariation module with ASF from an internally generated multiplet.")
println("Ci)  Test of the DecayYields module with ASF from an internally generated multiplet.")
println("Ck)  Test of the MultipolePolarizibility module with ASF from an internally generated multiplet.")
println("")
println("D)   Tests & development of atomic processes")
println("--------------------------------------------")
println("Da)  Test of the PhotoEmission module with ASF from an internally generated initial- and final-state multiplet.")
println("Db)  Test of the PhotoExcitation module with ASF from an internally generated initial- and final-state multiplet.")
println("Dc)  Test of the Photoionization module with ASF from an internally generated initial- and final-state multiplet.")
println("Dcp) Test of the Photoionization and PlasmaShift module with ASF from an internally generated initial- and final-state multiplet.")
println("Dd)  Test of the PhotoRecombination module with ASF from an internally generated initial- and final-state multiplet.")
println("De)  Test of the Auger module with ASF from an internally generated initial- and final-state multiplet.")
println("Dep) Test of the Auger and PlasmaShift modules with ASF from an internally generated initial- and final-state multiplet.")
println("Df)  Test of the Dielectronic module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
println("Dg)  Test of the PhotoExcitationFluores module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
println("Dh)  Test of the PhotoExcitationAutoIon module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
println("Di)  Test of the RayleighCompton module with ASF from an internally generated initial- and final-state multiplets.")
println("Dj)  Test of the MultiPhotonDeExcitation module with ASF from an internally generated initial- and final-state multiplet.")
println("Dk)  Test of the CoulombExcitation module with ASF from an internally generated initial- and final-state multiplet.")
println("Dl)  Test of the PhotoIonizationFluores module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
println("Dm)  Test of the PhotoIonizationAutoion module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
println("Dn)  Test of the ImpactExcitation module with ASF from an internally generated initial- and final-state multiplet.")
println("Dp)  Test of the ImpactExcitationAutoion module with ASF from an internally generated initial-, intermediate and final-state multiplet.")
println("Dq)  Test of the RadiativeAuger module with ASF from an internally generated initial- and final-state multiplet.")
println("Dr)  Test of the MultiPhotonIonization module with ASF from an internally generated initial- and final-state multiplet.")
println("Ds)  Test of the MultiPhotonDoubleIon module with ASF from an internally generated initial- and final-state multiplet.")
println("Dt)  Test of the InternalConversion module with ASF from an internally generated initial- and final-state multiplet.")
println("Du) Test of the CoulombIonization module with ASF from an internally generated initial- and final-state multiplet.")
println("Dv) Test of the PairAnnihilation1Photon.jl module with ASF from an internally generated initial- and final-state multiplet.")
println("")
println("E)  Tests & development of atomic cascades")
println("------------------------------------------")
println("Ea) Three-step cascade after 1s-3p photo-excitation of Si^+: Configuration model.")
println("Eb) Simulation of the four-step cascade after 1s-3p photo-excitation of Si^+: Configuration model.")
println("Ec) Three-step cascade computations and simulations for the decay of the neon 1s^-1 3p hole states: AverageSCA model.")
println("Ed) Photoabsorption of Ne^+: AverageSCA model.")
println("")
println("F)  Tests & development of symbolic evaluations of Racah expressions")
println("--------------------------------------------------------------------")
println("Fa) Symbolic evaluation by means of special values & recursion relations of the Wigner n-j symbols.")
println("Fb) Symbolic simplification of Racah expressions by means of sum rules.")
println("Fc) Symbolic simplification of recoupling coefficients by means of sum rules.")
println("Fd) Symbolic simplification of Racah expressions including spherical harmonics and Wigner rotation matrices.")
println("Fe) Symbolic simplification of spherical tensor operators, matrix elements and spherical amplitudes.")
println("")
println("G)  Tests & development of semiempirical computations")
println("-----------------------------------------------------")
println("Ga) Tests of the semi-empirical computation of ADK rates.")
println("")
println("H)  Tests & development of atomic response computations")
println("-----------------------------------------------------")
println("Ha) Tests of the HighHarmonic module to calculate a high-harmonic spectrum in single-electron approximation.")
println("Hb) Tests of the StrongField module to calculate energy and momentum distributions of photoelectron in SFA.")


text:
## Sensor model that returns position and velocities of other traffic participants
# Gaussian Noise
# FP/FN
# Late detection?
# Position noise depends on distance

"""
Sensor model that returns position and velocities of other traffic participants
Gaussian Noise
FP/FN
Late detection?
Position noise depends on distance
"""
@with_kw struct GaussianSensor{P <: AbstractNoiseModel, V <: AbstractNoiseModel, R<:AbstractRNG} <: AbstractSensor
    pos_noise::P = LinearNoise()
    vel_noise::V = LinearNoise()
    false_positive_rate::Float64 = 0.1
    false_negative_rate::Float64 = 0.1
    rng::R = Random.GLOBAL_RNG
end

function measure(sensor::GaussianSensor, ego::Entity, scene::EntityScene{S,D,I}, roadway::Roadway, obstacles::Vector{ConvexPolygon}) where {S,D,I}
    obs = Vector{Entity{VehicleState, D, I}}()
    sizehint!(obs, 10) #XXX
    for veh in scene
        if veh.id == ego.id
            continue
        end
        zveh = measure(sensor, ego, veh, roadway, obstacles, sensor.rng)
        if zveh != nothing
            push!(obs, Entity(zveh, veh.def, veh.id))
        end
    end
    if isempty(obs)
        false_car = car_false_positive(sensor.false_positive_rate, roadway, scene, obstacles, sensor.rng)
        false_car == nothing ? nothing : push!(obs, false_car)
        false_ped = ped_false_positive(sensor.false_positive_rate, roadway, scene, obstacles, sensor.rng)
        false_ped == nothing ? nothing : push!(obs, false_ped)
    end
    return obs
end

function measure(sensor::GaussianSensor, ego::Entity, veh::Entity, roadway::Roadway, obstacles::Vector{ConvexPolygon}, rng::AbstractRNG)
    occluded = occlusion_checker(ego.state, veh.state, obstacles)
    FN = rand(rng) < sensor.false_negative_rate
    if occluded || FN
        return nothing
    else
        dist = sqrt((posg(ego.state).x - posg(veh.state).x)^2 + (posg(ego.state).y - posg(veh.state).y)^2)
        pos_std = noise(sensor.pos_noise, dist)
        vel_std = noise(sensor.vel_noise, dist)
        z_pos = VecSE2(posg(veh.state).x + pos_std*randn(rng),
                       posg(veh.state).y + pos_std*randn(rng),
                       posg(veh.state).θ)
        z_v = veh.state.v + vel_std*randn(rng)
        return VehicleState(z_pos, roadway, z_v)
    end
end

function obs_weight(sensor::GaussianSensor, ego::VehicleState, obs::Union{Nothing, VehicleState}, veh::Union{Nothing, VehicleState}, obstacles::Vector{ConvexPolygon})
    weight = 1.0
    visible = veh != nothing && !occlusion_checker(ego, veh, obstacles)
    if !visible
        if obs == nothing
            return weight*(1. - sensor.false_positive_rate)
        else
            return weight*sensor.false_positive_rate # should handle false positive
        end
    else
        if obs == nothing
            return sensor.false_positive_rate
        else
            dist = sqrt((ego.posG.x - veh.posG.x)^2 + (ego.posG.y - veh.posG.y)^2)
            pos_std = noise(sensor.pos_noise, dist)
            vel_std = noise(sensor.vel_noise, dist)
            weight *= (1 - sensor.false_positive_rate)
            weight *= 1/(sqrt(2*pi*pos_std*vel_std))*exp(-1/pos_std*(veh.posG.x - obs.posG.x)^2 -1/pos_std*(veh.posG.y - obs.posG.y)^2 - 1/vel_std*(veh.v - obs.v)^2)
            return weight
        end
    end
end

@with_kw struct GaussianSensorOverlay{E<:Entity}
    sensor::GaussianSensor = GaussianSensor()
    o::Vector{E} = Vector{E}()
    color::Colorant = RGBA(0.976, 0.592, 0.122, 0.7) # orange
end

function AutomotiveVisualization.add_renderable!(rendermodel::RenderModel, overlay::GaussianSensorOverlay, scene::Scene, roadway::R) where R
    for veh in overlay.o
        p = posg(veh.state)
        # add_instruction!(rendermodel, render_vehicle, (p.x, p.y, p.θ, veh.def.length, veh.def.width, overlay.color))
        add_instruction!(rendermodel, render_vehicle,
                            (p.x, p.y, p.θ,
                            veh.def.length, veh.def.width, overlay.color,ARGB(1.,1.,1.,0.),ARGB(1.,1.,1.,0.)))
    end
    return rendermodel
end


text:
@defcomp ocean begin
    # Sea-level rise in meters
    sea = Variable(index=[time])

    lifesea = Parameter(default = 500)

    seas = Parameter(default = 2)

    delaysea = Variable()

    # Temperature incrase in C°
    temp = Parameter(index=[time])

    function run_timestep(p, v, d, t)

        if is_first(t)
            # Delay in sea-level rise
            v.delaysea = 1.0 / p.lifesea
            v.sea[t] = 0.0
        else
            # Calculate sea level rise
            ds = v.delaysea * p.seas * p.temp[t] - v.delaysea * v.sea[t - 1]

            v.sea[t] = v.sea[t - 1] + ds
        end
    end
end


text:
const PlusMinus  = "±"
const PosInfinity = "+∞"
const NegInfinity = "-∞"
const Empty      = "∅"
const Full       = "▣"
const Exactly    = "⌁"
const Inexactly  = "~"

# select characters corresponding to boundries as specified (lo,hi)
#   true for Open, false for Closed
#                          ClCl            ClOp            OpCl          OpOp

#const Delimiters = [ ( "⟨" , "⟩" ) , ( "⟨" , "⟫" ),  ( "⟪" , "⟩" ), ( "⟪" , "⟫" )  ];
#const Delimiters =  [ ( "⎰" , "⎱" ) , ( "⎰" , "⎰" ) ,  ( "⎱" , "⎱" ), ( "⎱" , "⎰" )  ];
const Delimiters =  [ ( "⟨" , "⟩" ) , ( "⟨" , "⦒" ) ,  ( "⦑" , "⟩" ), ( "⦑" , "⦒" )  ];

@inline delimiters(loIsOpen::Bool, hiIsOpen::Bool) =
      Delimiters[ one(Int8)+(reinterpret(Int8,loIsOpen)<<1)+hiIsOpen ]

Base.string(::Type{ClCl}) = "ClCl"
Base.string(::Type{ClOp}) = "ClOp"
Base.string(::Type{OpCl}) = "OpCl"
Base.string(::Type{OpOp}) = "OpOp"

function show{G<:Grasp,R<:Real}(io::IO, x::Rvl{G,R})
    tiestr = ", "
    Gstr = string(G)
    s = if (x.lo != x.hi)
            string(Gstr,"(", x.lo, tiestr, x.hi, ")")
        elseif G==ClCl
            string(x.lo)
        else
            string(Gstr,"(", x.lo, ")")
        end
    print(io, s)
end

#=
function show{G<:Grasp,R<:Real}(io::IO, x::Rvl{G,R})
    tiestr = ", "
    loDelim,hiDelim = delimiters(boundaries(G)...)
    s = if (x.lo != x.hi)
            string(loDelim, x.lo, tiestr, x.hi, hiDelim)
        elseif G==ClCl
            string(x.lo)
        else
            string(loDelim, x.lo, hiDelim)
        end
    print(io, s)
end
=#

function showcompact{G<:Grasp,R<:Real}(io::IO, x::Rvl{G,R})
    tiestr = ", "
    Gstr = string(G)
    lo = @sprintf("%7.5g", x.lo); lo = strip(lo)
    hi = @sprintf("%7.5g", x.hi); hi = strip(hi)
    s = if (x.lo != x.hi)
            string(Gstr,"(", lo, tiestr, hi, ")")
        elseif G==ClCl
            string(x.lo)
        else
            string(Gstr,"(", lo, ")")
        end
    print(io, s)
end

function showmidrad{G<:Grasp,R<:Real}(io::IO, x::Rvl{G,R})
    md = @sprintf("%7.5g", mid(x)); md = strip(md)
    rd = @sprintf("%7.5g", rad(x)); rd = strip(rd)
    s = string(md, PlusMinus, rd)
    print(io, s)
end

showmidrad{G<:Grasp,R<:Real}(x::Rvl{G,R}) = showmidrad(STDOUT,x)

function showmidpointradius{G<:Grasp,R<:Real}(io::IO, x::Rvl{G,R})
    tiesym = (Q==INEXACT) ? Inexactly : Exactly
    md = mid(x)
    rd = rad(x)
    s  = string(md, PlusMinus, rd)
    print(io, s)
end

showmidpointradius{G<:Grasp,R<:Real}(x::Rvl{G,R}) = showmidpointradius(STDOUT,x)


text:
# TODO: Needs testing
# Taken from Clustering.jl/src/randindex.jl)

struct Contingency
  t::Matrix{Int}
  x::Vector{Int}
  y::Vector{Int}
end

function Contingency(M::Dict{Tuple{Int, Int}, Int})
  x, y = zip(keys(M)...)
  ux = sort(unique(x))
  uy = sort(unique(y))
  t = zeros(Int, length(ux), length(uy))
  for i in ux
    for j in uy
      idx_i = findfirst(index -> index == i, ux)
      idx_j = findfirst(index -> index == j, uy)
      if haskey(M, (i, j))
        t[idx_i, idx_j] = M[i, j] + 0
      end
    end
  end
  return Contingency(t, ux, uy)
end

function table(x::AbstractVector, y::AbstractVector)
  n = length(x)
  @assert length(y) == n

  range_x = minimum(x):maximum(x)
  range_y = minimum(y):maximum(y)
  C = StatsBase.counts(x, y, (range_x, range_y))

  M = Dict{Tuple{Int, Int}, Int}()
  for i in 1:length(range_x)
    for j in 1:length(range_y)
      c_ij = C[i, j]
      if c_ij > 0
        M[(range_x[i], range_y[j])] = C[i, j]
      end
    end
  end

  return M
end

ctable(x, y) = Contingency(table(x, y))

function pretty(C::Contingency)
  M = zeros(Union{Int, Missing}, length(C.x) + 1, length(C.y) + 1)
  M[2:end, 2:end] .= C.t
  M[2:end, 1] = C.x
  M[1, 2:end] = C.y
  M[1, 1] = missing
  M
end

#= TEST
x = [1,3,5,2,4,5,5,1,1]
y = [2,4,6,1,1,1,1,1,2]
C = ctable(x, y)

C.t
C.x
C.y
pretty(C)
=#

"""
Adjusted rand index ∈ [0, 1] between two sets
of cluster labels (x and y)
"""
function ari(x::AbstractVector, y::AbstractVector)
  n = length(x)
  @assert length(y) == n

  C = ctable(x, y)
  c = C.t

  n = round(Int,sum(c))
  nis = sum(abs2, sum(c, dims=2))        # sum of squares of sums of rows
  njs = sum(abs2, sum(c, dims=1))        # sum of squares of sums of columns

  t1 = binomial(n,2)            # total number of pairs of entities
  t2 = sum(c.^2)                # sum over rows & columnns of nij^2
  t3 = .5*(nis+njs)

  # Expected index (for adjustment)
  nc = (n*(n^2+1)-(n+1)*nis-(n+1)*njs+2*(nis*njs)/n)/(2*(n-1))

  A = t1+t2-t3;        # no. agreements
  D = -t2+t3;          # no. disagreements

  if t1 == nc
      # avoid division by zero; if k=1, define Rand = 0
      ARI = 0
  else
      # adjusted Rand - Hubert & Arabie 1985
      ARI = (A-nc)/(t1-nc)
  end

  return ARI
end


text:
using CUDA
using KernelAbstractions: @kernel, @index, CUDADevice
using Oceananigans.Architectures: device, GPU
using Oceananigans.Utils: work_layout

function set!(Φ::NamedTuple; kwargs...)
    for (fldname, value) in kwargs
        ϕ = getproperty(Φ, fldname)
        set!(ϕ, value)
    end
    return nothing
end

set!(u::AbstractField, v::Number) = @. u.data.parent = v

set!(u::AbstractField{X, Y, Z, A}, v::AbstractField{X, Y, Z, A}) where {X, Y, Z, A} =
    @. u.data.parent = v.data.parent

# Niceties
const AbstractCPUField =
    AbstractField{X, Y, Z, A, G} where {X, Y, Z, A<:OffsetArray{T, D, <:Array} where {T, D}, G}

const AbstractReducedCPUField =
    AbstractReducedField{X, Y, Z, A, G} where {X, Y, Z, A<:OffsetArray{T, D, <:Array} where {T, D}, G}

"Set the CPU field `u` to the array `v`."
function set!(u::AbstractCPUField, v::Array)

    Sx, Sy, Sz = size(u)
    for k in 1:Sz, j in 1:Sy, i in 1:Sx
        u[i, j, k] = v[i, j, k]
    end

    return nothing
end

""" Returns an AbstractReducedField on the CPU. """
function similar_cpu_field(u::AbstractReducedField)
    FieldType = typeof(u).name.wrapper
    return FieldType(location(u), CPU(), u.grid; dims=u.dims)
end

""" Set the CPU field `u` data to the function `f(x, y, z)`. """
set!(u::AbstractCPUField, f::Function) = interior(u) .= f.(nodes(u; reshape=true)...)

#####
##### set! for fields on the GPU
#####

@hascuda begin
    const AbstractGPUField =
        AbstractField{X, Y, Z, A, G} where {X, Y, Z, A<:OffsetArray{T, D, <:CuArray} where {T, D}, G}

    const AbstractReducedGPUField =
        AbstractReducedField{X, Y, Z, A, G} where {X, Y, Z, A<:OffsetArray{T, D, <:CuArray} where {T, D}, G}

    """ Returns a field on the CPU with `nothing` boundary conditions. """
    function similar_cpu_field(u)
        FieldType = typeof(u).name.wrapper
        return FieldType(location(u), CPU(), u.grid, nothing)
    end

    """ Set the GPU field `u` to the array `v`. """
    function set!(u::AbstractGPUField, v::Array)
        v_field = similar_cpu_field(u)

        set!(v_field, v)
        set!(u, v_field)

        return nothing
    end

    """ Set the GPU field `u` to the CuArray `v`. """
    function set!(u::AbstractGPUField, v::CuArray)

        launch!(GPU(), u.grid, :xyz, _set_gpu!, u.data, v, u.grid,
                include_right_boundaries=true, location=location(u))

        return nothing
    end

    @kernel function _set_gpu!(u, v, grid)
        i, j, k = @index(Global, NTuple)
        @inbounds u[i, j, k] = v[i, j, k]
    end

    """ Set the CPU field `u` data to the GPU field data of `v`. """
    set!(u::AbstractCPUField, v::AbstractGPUField) = u.data.parent .= Array(v.data.parent)

    """ Set the GPU field `u` data to the CPU field data of `v`. """
    set!(u::AbstractGPUField, v::AbstractCPUField) = copyto!(u.data.parent, v.data.parent)

    """ Set the GPU field `u` data to the function `f(x, y, z)`. """
    function set!(u::AbstractGPUField, f::Function)
        # Create a temporary field with bcs = nothing.
        v_field = similar_cpu_field(u)

        set!(v_field, f)
        set!(u, v_field)

        return nothing
    end
end


text:
export AbstractSequence, Sequence, sequence, generator

abstract type AbstractSequence{T,D} end

struct Generator{D}
    generator::D
end

(G::Generator)(i) = G.generator(i)
Base.getindex(G::Generator, i) = G(i)

mutable struct Sequence{T,D} <: AbstractSequence{T,D}
    sequence::Vector{T}
    generator::Generator{D}
end

Sequence(generator::D) where D<:Base.Callable = Sequence([generator(1)],Generator(generator))

sequence(S::AbstractSequence) = S.sequence
generator(S::AbstractSequence) = S.generator

Base.eltype(::AbstractSequence{T}) where T = T
Base.size(::AbstractSequence) = (Inf,)
Base.length(::AbstractSequence) = Inf

Base.getindex(S::AbstractSequence,i::Integer) = (generate!(S,i); sequence(S)[i])
Base.getindex(S::AbstractSequence,ir::AbstractRange) = (generate!(S,maximum(ir)); sequence(S)[ir])

max_stars_repo_name: MikaelSlevinsky/SequenceTransformations

max_issues_repo_name: MikaelSlevinsky/SequenceTransformations

max_forks_repo_name: MikaelSlevinsky/SequenceTransformations


text:
# Scheduling of processes in GeMM

"""
    simulate!(world, timesteps)

This is the central function of the model with the main event loop. It defines
the scheduling for all submodels and output functions.
"""
function simulate!(world::Array{Patch,1}, timesteps::Int=1000, timeoffset::Int = 0)
    simlog("Starting simulation.")
    checkviability!(world)
    for t in (timeoffset + 1):(timeoffset + timesteps)
        simlog("UPDATE $t")
        # ecological processes are outsourced to specialised methods below
        if setting("mode") == "default"
            defaultexperiment(world)
        elseif setting("mode") == "invasion"
            invasionexperiment(world, t)
        elseif setting("mode") == "zosterops"
            zosteropsexperiment(world)
        else
            simlog("Mode setting not recognised: $(setting("mode"))", 'e')
        end
        if iszero(mod(t, setting("outfreq")))
            writedata(world, t)
        end
        if iszero(mod(t, setting("fastaoutfreq"))) && (setting("fasta") != "off")
            writefasta(world, t)
        end
    end
end

"""
    defaultexperiment(world)

The standard annual update procedure, designed primarily for plant communities.
"""
function defaultexperiment(world::Array{Patch,1})
    establish!(world)
    survive!(world)
    grow!(world)
    compete!(world)
    reproduce!(world)
    if setting("mutate")
        mutate!(world)
    end
    disperse!(world)
    changehabitat!(world) # model output
end

"""
    invasionexperiment(world)

The annual update procedure for the invasion experiments.
"""
function invasionexperiment(world::Array{Patch,1}, t::Int)
    establish!(world)
    survive!(world)
    grow!(world)
    compete!(world)
    reproduce!(world)
    if 0 < setting("burn-in") < t
        disturb!(world)
        invade!(world)
    end
    disperse!(world)
end

"""
    zosteropsexperiment(world)

The annual update procedure for the Zosterops experiments, this time for bird populations.
"""
function zosteropsexperiment(world::Array{Patch,1})
    # The first four processes are patch-internal and can therefore be parallelised
    # Note: multithreading requires calling Julia with the -p parameter (and -t?)
    #FIXME we're not entirely certain this works as advertised - better not use it
    # until further tested
    Threads.@threads for patch in world
        establish!(patch)
        survive!(patch)
        zreproduce!(patch)
        if setting("mutate")
            mutate!(patch)
        end
    end
    zdisperse!(world)
end

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
@testset "2-area, 11-bus, 4-generator benchmark system Parsing " begin
    test_dir = mktempdir()
    ger4_raw_file = joinpath(PSSE_RAW_DIR, "Benchmark_4ger_33_2015.RAW")
    ger4_dyr_file = joinpath(PSSE_DYR_DIR, "Benchmark_4ger_33_2015.dyr")
    sys = PSB.build_system(PSB.PSSETestSystems, "psse_Benchmark_4ger_33_2015_sys")
    @test_throws PSY.DataFormatError System(ger4_raw_file, ger4_raw_file)
    @test_throws PSY.DataFormatError System(ger4_dyr_file, ger4_raw_file)
    dyn_injectors = get_components(DynamicInjection, sys)
    @test length(dyn_injectors) == 4
    for g in dyn_injectors
        m = PSY.get_machine(g)
        @test typeof(m) <: RoundRotorExponential
    end
    path = joinpath(test_dir, "test_dyn_system_serialization.json")
    to_json(sys, path)
    parsed_sys = System(path)
    dyn_injectors = get_components(DynamicInjection, parsed_sys)
    @test length(dyn_injectors) == 4
    for g in dyn_injectors
        m = PSY.get_machine(g)
        @test typeof(m) <: RoundRotorExponential
    end
end

@testset "GENCLS dyr parsing" begin
    threebus_raw_file = joinpath(PSSE_TEST_DIR, "ThreeBusNetwork.raw")
    gencls_dyr_file = joinpath(PSSE_TEST_DIR, "TestGENCLS.dyr")
    nogencls_dyr_file = joinpath(PSSE_TEST_DIR, "Test-NoCLS.dyr")
    sys_gencls = PSB.build_system(PSB.PSSETestSystems, "psse_3bus_gen_cls_sys")
    sys_nogencls = PSB.build_system(PSB.PSSETestSystems, "psse_3bus_no_cls_sys")

    #Check that generator at bus 102 (H = 0) is a Source, and not ThermalStandard.
    @test isnothing(get_component(ThermalStandard, sys_gencls, "generator-102-1"))
    @test typeof(get_component(StaticInjection, sys_gencls, "generator-102-1")) <: Source
    #Check that generator at bus 102 (not provided in dyr data) is a Source.
    @test isnothing(get_component(ThermalStandard, sys_nogencls, "generator-102-1"))
    @test typeof(get_component(StaticInjection, sys_nogencls, "generator-102-1")) <: Source

    #Check Parameters of GENROE raw file
    genroe_gen = get_component(DynamicInjection, sys_gencls, "generator-101-1")
    m = PSY.get_machine(genroe_gen)
    s = PSY.get_shaft(genroe_gen)
    @test typeof(m) <: RoundRotorExponential

    @test get_R(m) == 0.0
    @test get_Td0_p(m) == 8.0
    @test get_Td0_pp(m) == 0.03
    @test get_Tq0_p(m) == 0.4
    @test get_Tq0_pp(m) == 0.05
    @test get_Xd(m) == 1.8
    @test get_Xq(m) == 1.7
    @test get_Xd_p(m) == 0.3
    @test get_Xd_pp(m) == 0.25
    @test get_Xl(m) == 0.2
    @test get_H(s) == 6.5
    @test get_D(s) == 0.0

    #Check Parameters of other GENCLS file
    gencls_gen = get_component(DynamicInjection, sys_gencls, "generator-103-1")
    m = PSY.get_machine(gencls_gen)
    s = PSY.get_shaft(gencls_gen)
    @test typeof(m) <: BaseMachine

    @test get_R(m) == 0.0
    @test get_Xd_p(m) == 0.2
    @test get_H(s) == 3.1
    @test get_D(s) == 2.0

    #Check Parameters of other GENCLS file
    gencls_gen = get_component(DynamicInjection, sys_nogencls, "generator-103-1")
    m = PSY.get_machine(gencls_gen)
    s = PSY.get_shaft(gencls_gen)
    @test typeof(m) <: BaseMachine

    @test get_R(m) == 0.0
    @test get_Xd_p(m) == 0.2
    @test get_H(s) == 3.1
    @test get_D(s) == 2.0
end

@testset "2000-Bus Parsing" begin
    test_dir = mktempdir()
    texas2000_raw_file = joinpath(TAMU_DIR, "ACTIVSg2000.RAW")
    texas2000_dyr_file = joinpath(PSSE_DYR_DIR, "ACTIVSg2000_dynamics.dyr")
    sys = PSB.build_system(PSB.PSSETestSystems, "psse_ACTIVSg2000_sys")
    for g in get_components(ThermalStandard, sys)
        @test !isnothing(get_dynamic_injector(g))
    end
    path = joinpath(test_dir, "test_dyn_system_serialization_2000.json")
    to_json(sys, path)
    parsed_sys = System(path)
    dyn_injectors = get_components(DynamicInjection, parsed_sys)
    @test length(dyn_injectors) == 435
    for g in get_components(ThermalStandard, parsed_sys)
        @test !isnothing(get_dynamic_injector(g))
    end
end


text:
# Automatically generated using Clang.jl

const rocfft_version_major = 0
const rocfft_version_minor = 9
const rocfft_version_patch = 7
const rocfft_plan_t = Cvoid
const rocfft_plan = Ptr{rocfft_plan_t}
const rocfft_plan_description_t = Cvoid
const rocfft_plan_description = Ptr{rocfft_plan_description_t}
const rocfft_execution_info_t = Cvoid
const rocfft_execution_info = Ptr{rocfft_execution_info_t}

@cenum rocfft_status_e::UInt32 begin
    rocfft_status_success = 0
    rocfft_status_failure = 1
    rocfft_status_invalid_arg_value = 2
    rocfft_status_invalid_dimensions = 3
    rocfft_status_invalid_array_type = 4
    rocfft_status_invalid_strides = 5
    rocfft_status_invalid_distance = 6
    rocfft_status_invalid_offset = 7
end

const rocfft_status = rocfft_status_e

@cenum rocfft_transform_type_e::UInt32 begin
    rocfft_transform_type_complex_forward = 0
    rocfft_transform_type_complex_inverse = 1
    rocfft_transform_type_real_forward = 2
    rocfft_transform_type_real_inverse = 3
end

const rocfft_transform_type = rocfft_transform_type_e

@cenum rocfft_precision_e::UInt32 begin
    rocfft_precision_single = 0
    rocfft_precision_double = 1
end

const rocfft_precision = rocfft_precision_e

@cenum rocfft_result_placement_e::UInt32 begin
    rocfft_placement_inplace = 0
    rocfft_placement_notinplace = 1
end

const rocfft_result_placement = rocfft_result_placement_e

@cenum rocfft_array_type_e::UInt32 begin
    rocfft_array_type_complex_interleaved = 0
    rocfft_array_type_complex_planar = 1
    rocfft_array_type_real = 2
    rocfft_array_type_hermitian_interleaved = 3
    rocfft_array_type_hermitian_planar = 4
    rocfft_array_type_unset = 5
end

const rocfft_array_type = rocfft_array_type_e

@cenum rocfft_execution_mode_e::UInt32 begin
    rocfft_exec_mode_nonblocking = 0
    rocfft_exec_mode_nonblocking_with_flush = 1
    rocfft_exec_mode_blocking = 2
end

const rocfft_execution_mode = rocfft_execution_mode_e

@cenum rocfft_layer_mode_::UInt32 begin
    rocfft_layer_mode_none = 0
    rocfft_layer_mode_log_trace = 1
    rocfft_layer_mode_log_bench = 2
    rocfft_layer_mode_log_profile = 4
end

const rocfft_layer_mode = rocfft_layer_mode_


text:
#=
sweeping:
- Julia version:
- Author: qd4314
- Date: 2019-03-19
=#

function sweeping(
    KS::KineticSolver,
    psi::Array{Float64,3},
    Source::Array{Float64,3},
    ScatteringKernel::Array{Float64,2},
)
    rhs = zeros(KS.nquadpoints, KS.nx, KS.ny)
    io = open("mytestfile.txt", "a")
    #println("")
    #println("sweeping is called")
    # loop over pseudo-time
    for k = 1:KS.nSweepingSteps
        # setup rhs
        print(" ")
        println(k)
        @inbounds for j in KS.rangex, i in KS.rangey
            rhs[:, i, j] =
                Source[:, i, j] +
                KS.convolutionmagnitude .* KS.kernelTimesWeightsSparse * (psi[:, i, j])
        end

        # check residual
        psiOld = deepcopy(psi)
        # perform sweep
        solveFluxSweeping!(KS, psi, rhs)

        # check residual
        #println("-> Sweeping norm: ",norm(psi-psiOld)*KS.contractionEstimate)
        if norm(psi - psiOld) * KS.contractionEstimate < KS.epsSweeping
            writedlm(io, k)
            close(io)
            return psi
        end
    end
    return psi
end

function sweepingWithKernels(
    KS::KineticSolver,
    psi::Array{Float64,3},
    Source::Array{Float64,3},
    asKernel::Array{Float64,2},
    physKernel::Array{Float64,2},
)
    rhs = zeros(KS.nquadpoints, KS.nx, KS.ny)
    io = open("mytestfile.txt", "a")
    #println("")
    #println("sweeping is called")
    # loop over pseudo-time
    for k = 1:KS.nSweepingSteps
        # setup rhs
        print(" ")
        println(k)
        @inbounds for j in KS.rangex, i in KS.rangey
            rhs[:, i, j] =
                Source[:, i, j] +
                KS.SigmaS[i, j] .* physKernel * psi[:, i, j] +
                KS.convolutionmagnitude .* asKernel * psi[:, i, j]
        end

        # check residual
        psiOld = deepcopy(psi)
        # perform sweep
        solveFluxSweeping!(KS, psi, rhs)
        #solveFluxSweepingFirstOrder!(KS,psi,rhs);

        # check residual
        println(
            "-> Sweeping norm: ",
            norm(psi - psiOld),
            " ; constraction: ",
            KS.contractionEstimate,
        )
        if norm(psi - psiOld) * KS.contractionEstimate < KS.epsSweeping
            writedlm(io, k)
            close(io)
            return psi
        end
    end
    return psi
end

function solveFluxSweepingBiggerStencil!(
    KS::KineticSolver,
    psi::Array{Float64,3},
    rhs::Array{Float64,3},
)
    # computes the numerical flux over cell boundaries for each ordinate
    # for faster computation, we split the iteration over quadrature points
    # into four different blocks: North West, Nort East, Sout West, South East
    # this corresponds to the direction the ordinates point to
    idxPosPos = findall((KS.Q.pointsxyz[:, 1] .>= 0.0) .& (KS.Q.pointsxyz[:, 2] .>= 0.0))
    idxPosNeg = findall((KS.Q.pointsxyz[:, 1] .>= 0.0) .& (KS.Q.pointsxyz[:, 2] .< 0.0))
    idxNegPos = findall((KS.Q.pointsxyz[:, 1] .< 0.0) .& (KS.Q.pointsxyz[:, 2] .>= 0.0))
    idxNegNeg = findall((KS.Q.pointsxyz[:, 1] .< 0.0) .& (KS.Q.pointsxyz[:, 2] .< 0.0))

    # PosPos
    @inbounds for j in KS.rangex, i in KS.rangey, q in idxPosPos
        s0 = psi[q, i, j-3]
        s1 = psi[q, i, j-2]
        s2 = psi[q, i, j-1]

        fluxExplY = KS.Q.pointsxyz[q, 2] ./ KS.dy .* (-s2 + 1 / 4 * s1)
        fluxExplY =
            fluxExplY - KS.Q.pointsxyz[q, 2] ./ KS.dy .* ((7 / 4) * s2 - s1 + 1 / 4 * s0)
        fluxImplY = KS.Q.pointsxyz[q, 2] ./ KS.dy .* (7 / 4)

        s0 = psi[q, i-3, j]
        s1 = psi[q, i-2, j]
        s2 = psi[q, i-1, j]

        fluxExplX = KS.Q.pointsxyz[q, 1] ./ KS.dx .* (-s2 + 1 / 4 * s1)
        fluxExplX =
            fluxExplX - KS.Q.pointsxyz[q, 1] ./ KS.dx .* ((7 / 4) * s2 - s1 + 1 / 4 * s0)
        fluxImplX = KS.Q.pointsxyz[q, 1] ./ KS.dx .* (7 / 4)

        fluxImplicit = fluxImplY + fluxImplX

        flux = fluxExplX + fluxExplY

        psi[q, i, j] =
            KS.dt .* (-flux + rhs[q, i, j]) /
            (1 + KS.dt .* (KS.SigmaT[i, j] + fluxImplicit))
    end
    #PosNeg
    @inbounds for j in Iterators.reverse(KS.rangex), i in KS.rangey, q in idxPosNeg
        s0 = psi[q, i, j+1]
        s1 = psi[q, i, j+2]
        s2 = psi[q, i, j+3]

        fluxExplY = KS.Q.pointsxyz[q, 2] ./ KS.dy .* ((7 / 4) * s0 - s1 + (1 / 4) * s2)
        fluxExplY = fluxExplY - KS.Q.pointsxyz[q, 2] ./ KS.dy .* (-s0 + (1 / 4) * s1)
        fluxImplY = -KS.Q.pointsxyz[q, 2] ./ KS.dy .* (7 / 4)

        s0 = psi[q, i-3, j]
        s1 = psi[q, i-2, j]
        s2 = psi[q, i-1, j]

        fluxExplX = KS.Q.pointsxyz[q, 1] ./ KS.dx .* (-s2 + 1 / 4 * s1)
        fluxExplX =
            fluxExplX - KS.Q.pointsxyz[q, 1] ./ KS.dx .* ((7 / 4) * s2 - s1 + 1 / 4 * s0)
        fluxImplX = KS.Q.pointsxyz[q, 1] ./ KS.dx .* (7 / 4)

        fluxImplicit = fluxImplY + fluxImplX

        flux = fluxExplX + fluxExplY

        psi[q, i, j] =
            KS.dt .* (-flux + rhs[q, i, j]) /
            (1 + KS.dt .* (KS.SigmaT[i, j] + fluxImplicit))
    end

    # NegPos
    @inbounds for j in KS.rangex, i in Iterators.reverse(KS.rangey), q in idxNegPos
        s0 = psi[q, i, j-3]
        s1 = psi[q, i, j-2]
        s2 = psi[q, i, j-1]

        fluxExplY = KS.Q.pointsxyz[q, 2] ./ KS.dy .* (-s2 + 1 / 4 * s1)
        fluxExplY =
            fluxExplY - KS.Q.pointsxyz[q, 2] ./ KS.dy .* ((7 / 4) * s2 - s1 + 1 / 4 * s0)
        fluxImplY = KS.Q.pointsxyz[q, 2] ./ KS.dy .* (7 / 4)

        s0 = psi[q, i+1, j]
        s1 = psi[q, i+2, j]
        s2 = psi[q, i+3, j]

        fluxExplX = KS.Q.pointsxyz[q, 1] ./ KS.dx .* ((7 / 4) * s0 - s1 + (1 / 4) * s2)
        fluxExplX = fluxExplX - KS.Q.pointsxyz[q, 1] ./ KS.dx .* (-s0 + 1 / 4 * s1)
        fluxImplX = -KS.Q.pointsxyz[q, 1] ./ KS.dx .* (7 / 4)

        fluxImplicit = fluxImplY + fluxImplX

        flux = fluxExplX + fluxExplY

        psi[q, i, j] =
            KS.dt .* (-flux + rhs[q, i, j]) /
            (1 + KS.dt .* (KS.SigmaT[i, j] + fluxImplicit))
    end

    # NegNeg
    @inbounds for j in Iterators.reverse(KS.rangex),
        i in Iterators.reverse(KS.rangey),
        q in idxNegNeg

        s0 = psi[q, i, j+1]
        s1 = psi[q, i, j+2]
        s2 = psi[q, i, j+3]

        fluxExplY = KS.Q.pointsxyz[q, 2] ./ KS.dy .* ((7 / 4) * s0 - s1 + (1 / 4) * s2)
        fluxExplY = fluxExplY - KS.Q.pointsxyz[q, 2] ./ KS.dy .* (-s0 + (1 / 4) * s1)
        fluxImplY = -KS.Q.pointsxyz[q, 2] ./ KS.dy .* (7 / 4)

        s0 = psi[q, i+1, j]
        s1 = psi[q, i+2, j]
        s2 = psi[q, i+3, j]

        fluxExplX = KS.Q.pointsxyz[q, 1] ./ KS.dx .* ((7 / 4) * s0 - s1 + (1 / 4) * s2)
        fluxExplX = fluxExplX - KS.Q.pointsxyz[q, 1] ./ KS.dx .* (-s0 + 1 / 4 * s1)
        fluxImplX = -KS.Q.pointsxyz[q, 1] ./ KS.dx .* (7 / 4)

        fluxImplicit = fluxImplY + fluxImplX

        flux = fluxExplX + fluxExplY

        psi[q, i, j] =
            KS.dt .* (-flux + rhs[q, i, j]) /
            (1 + KS.dt .* (KS.SigmaT[i, j] + fluxImplicit))
    end
end

function solveFluxSweeping!(KS::KineticSolver, psi::Array{Float64,3}, rhs::Array{Float64,3})
    # computes the numerical flux over cell boundaries for each ordinate
    # for faster computation, we split the iteration over quadrature points
    # into four different blocks: North West, Nort East, Sout West, South East
    # this corresponds to the direction the ordinates point to
    idxPosPos = findall((KS.Q.pointsxyz[:, 1] .>= 0.0) .& (KS.Q.pointsxyz[:, 2] .>= 0.0))
    idxPosNeg = findall((KS.Q.pointsxyz[:, 1] .>= 0.0) .& (KS.Q.pointsxyz[:, 2] .< 0.0))
    idxNegPos = findall((KS.Q.pointsxyz[:, 1] .< 0.0) .& (KS.Q.pointsxyz[:, 2] .>= 0.0))
    idxNegNeg = findall((KS.Q.pointsxyz[:, 1] .< 0.0) .& (KS.Q.pointsxyz[:, 2] .< 0.0))

    # PosPos
    @inbounds for j in KS.rangex, i in KS.rangey, q in idxPosPos
        s2 = psi[q, i, j-2]
        s1 = psi[q, i, j-1]

        fluxExplY = KS.Q.pointsxyz[q, 2] ./ KS.dy .* (-2 * s1 + 1 / 2 * s2)
        fluxImplY = KS.Q.pointsxyz[q, 2] ./ KS.dy .* (3 / 2)

        s2 = psi[q, i-2, j]
        s1 = psi[q, i-1, j]

        fluxExplX = KS.Q.pointsxyz[q, 1] ./ KS.dx .* (-2 * s1 + 1 / 2 * s2)
        fluxImplX = KS.Q.pointsxyz[q, 1] ./ KS.dx .* (3 / 2)

        fluxImplicit = fluxImplY + fluxImplX

        flux = fluxExplX + fluxExplY

        psi[q, i, j] =
            KS.dt .* (-flux + rhs[q, i, j]) /
            (1 + KS.dt .* (KS.SigmaT[i, j] + fluxImplicit))
    end
    #PosNeg
    @inbounds for j in Iterators.reverse(KS.rangex), i in KS.rangey, q in idxPosNeg
        s1 = psi[q, i, j+1]
        s2 = psi[q, i, j+2]

        fluxExplY = KS.Q.pointsxyz[q, 2] ./ KS.dy .* (2 * s1 - (1 / 2) * s2)
        fluxImplY = -KS.Q.pointsxyz[q, 2] ./ KS.dy .* (3 / 2)

        s2 = psi[q, i-2, j]
        s1 = psi[q, i-1, j]

        fluxExplX = KS.Q.pointsxyz[q, 1] ./ KS.dx .* (-2 * s1 + 1 / 2 * s2)
        fluxImplX = KS.Q.pointsxyz[q, 1] ./ KS.dx .* (3 / 2)

        fluxImplicit = fluxImplY + fluxImplX

        flux = fluxExplX + fluxExplY

        psi[q, i, j] =
            KS.dt .* (-flux + rhs[q, i, j]) /
            (1 + KS.dt .* (KS.SigmaT[i, j] + fluxImplicit))
    end

    # NegPos
    @inbounds for j in KS.rangex, i in Iterators.reverse(KS.rangey), q in idxNegPos
        s2 = psi[q, i, j-2]
        s1 = psi[q, i, j-1]

        fluxExplY = KS.Q.pointsxyz[q, 2] ./ KS.dy .* (-2 * s1 + 1 / 2 * s2)
        fluxImplY = KS.Q.pointsxyz[q, 2] ./ KS.dy .* (3 / 2)

        s1 = psi[q, i+1, j]
        s2 = psi[q, i+2, j]

        fluxExplX = KS.Q.pointsxyz[q, 1] ./ KS.dx .* (2 * s1 - (1 / 2) * s2)
        fluxImplX = -KS.Q.pointsxyz[q, 1] ./ KS.dx .* (3 / 2)

        fluxImplicit = fluxImplY + fluxImplX

        flux = fluxExplX + fluxExplY

        psi[q, i, j] =
            KS.dt .* (-flux + rhs[q, i, j]) /
            (1 + KS.dt .* (KS.SigmaT[i, j] + fluxImplicit))
    end

    # NegNeg
    @inbounds for j in Iterators.reverse(KS.rangex),
        i in Iterators.reverse(KS.rangey),
        q in idxNegNeg

        s1 = psi[q, i, j+1]
        s2 = psi[q, i, j+2]

        fluxExplY = KS.Q.pointsxyz[q, 2] ./ KS.dy .* (2 * s1 - (1 / 2) * s2)
        fluxImplY = -KS.Q.pointsxyz[q, 2] ./ KS.dy .* (3 / 2)

        s1 = psi[q, i+1, j]
        s2 = psi[q, i+2, j]

        fluxExplX = KS.Q.pointsxyz[q, 1] ./ KS.dx .* (2 * s1 - (1 / 2) * s2)
        fluxImplX = -KS.Q.pointsxyz[q, 1] ./ KS.dx .* (3 / 2)

        fluxImplicit = fluxImplY + fluxImplX

        flux = fluxExplX + fluxExplY

        psi[q, i, j] =
            KS.dt .* (-flux + rhs[q, i, j]) /
            (1 + KS.dt .* (KS.SigmaT[i, j] + fluxImplicit))
    end
end

function solveFluxSweepingFirstOrder!(
    KS::KineticSolver,
    psi::Array{Float64,3},
    rhs::Array{Float64,3},
)
    # computes the numerical flux over cell boundaries for each ordinate
    # for faster computation, we split the iteration over quadrature points
    # into four different blocks: North West, Nort East, Sout West, South East
    # this corresponds to the direction the ordinates point to
    idxPosPos = findall((KS.Q.pointsxyz[:, 1] .>= 0.0) .& (KS.Q.pointsxyz[:, 2] .>= 0.0))
    idxPosNeg = findall((KS.Q.pointsxyz[:, 1] .>= 0.0) .& (KS.Q.pointsxyz[:, 2] .< 0.0))
    idxNegPos = findall((KS.Q.pointsxyz[:, 1] .< 0.0) .& (KS.Q.pointsxyz[:, 2] .>= 0.0))
    idxNegNeg = findall((KS.Q.pointsxyz[:, 1] .< 0.0) .& (KS.Q.pointsxyz[:, 2] .< 0.0))

    # PosPos
    @inbounds for j in KS.rangex, i in KS.rangey, q in idxPosPos
        s1 = psi[q, i, j-2]
        s2 = psi[q, i, j-1]
        s3 = psi[q, i, j]
        s4 = psi[q, i, j+1]
        northflux = s3#+0.5 .*slopefit(s2,s3,s4)
        southflux = s2#+0.5 .*slopefit(s1,s2,s3)

        s1 = psi[q, i-2, j]
        s2 = psi[q, i-1, j]
        s3 = psi[q, i, j]
        s4 = psi[q, i+1, j]
        eastflux = s3#+0.5 .*slopefit(s2,s3,s4)
        westflux = s2#+0.5 .*slopefit(s1,s2,s3)

        fluxImpllicit = KS.Q.pointsxyz[q, 1] ./ KS.dx + KS.Q.pointsxyz[q, 2] ./ KS.dy

        flux =
            KS.Q.pointsxyz[q, 1] ./ KS.dx .* (-westflux) +
            KS.Q.pointsxyz[q, 2] ./ KS.dy .* (-southflux)
        psi[q, i, j] =
            KS.dt .* (-flux + rhs[q, i, j]) /
            (1 + KS.dt .* (KS.SigmaT[i, j] + fluxImpllicit))
    end
    #PosNeg
    @inbounds for j in Iterators.reverse(KS.rangex), i in KS.rangey, q in idxPosNeg
        s1 = psi[q, i, j-1]
        s2 = psi[q, i, j]
        s3 = psi[q, i, j+1]
        s4 = psi[q, i, j+2]
        northflux = s3#-0.5 .* slopefit(s2,s3,s4)
        southflux = s2#-0.5 .*slopefit(s1,s2,s3)

        s1 = psi[q, i-2, j]
        s2 = psi[q, i-1, j]
        s3 = psi[q, i, j]
        s4 = psi[q, i+1, j]
        eastflux = s3#+0.5 .*slopefit(s2,s3,s4)
        westflux = s2#+0.5 .*slopefit(s1,s2,s3)

        fluxImpllicit =
            KS.Q.pointsxyz[q, 1] ./ KS.dx + KS.Q.pointsxyz[q, 2] ./ KS.dy .* (-1)
        flux =
            KS.Q.pointsxyz[q, 1] ./ KS.dx .* (-westflux) +
            KS.Q.pointsxyz[q, 2] ./ KS.dy .* (northflux)
        psi[q, i, j] =
            KS.dt .* (-flux + rhs[q, i, j]) /
            (1 + KS.dt .* (KS.SigmaT[i, j] + fluxImpllicit))
    end

    # NegPos
    @inbounds for j in KS.rangex, i in Iterators.reverse(KS.rangey), q in idxNegPos
        s1 = psi[q, i, j-2]
        s2 = psi[q, i, j-1]
        s3 = psi[q, i, j]
        s4 = psi[q, i, j+1]
        northflux = s3#+0.5 .*slopefit(s2,s3,s4)
        southflux = s2#+0.5 .*slopefit(s1,s2,s3)

        s1 = psi[q, i-1, j]
        s2 = psi[q, i, j]
        s3 = psi[q, i+1, j]
        s4 = psi[q, i+2, j]
        eastflux = s3#-0.5 .*slopefit(s2,s3,s4)
        westflux = s2#-0.5 .*slopefit(s1,s2,s3)

        fluxImpllicit =
            KS.Q.pointsxyz[q, 1] ./ KS.dx .* (-1) + KS.Q.pointsxyz[q, 2] ./ KS.dy
        flux =
            KS.Q.pointsxyz[q, 1] ./ KS.dx .* (eastflux) +
            KS.Q.pointsxyz[q, 2] ./ KS.dy .* (-southflux)
        psi[q, i, j] =
            KS.dt .* (-flux + rhs[q, i, j]) /
            (1 + KS.dt .* (KS.SigmaT[i, j] + fluxImpllicit))
    end

    # NegNeg
    @inbounds for j in Iterators.reverse(KS.rangex),
        i in Iterators.reverse(KS.rangey),
        q in idxNegNeg

        s1 = psi[q, i, j-1]
        s2 = psi[q, i, j]
        s3 = psi[q, i, j+1]
        s4 = psi[q, i, j+2]
        northflux = s3#-0.5 .* slopefit(s2,s3,s4)
        southflux = s2#-0.5 .* slopefit(s1,s2,s3)

        s1 = psi[q, i-1, j]
        s2 = psi[q, i, j]
        s3 = psi[q, i+1, j]
        s4 = psi[q, i+2, j]
        eastflux = s3#-0.5 .* slopefit(s2,s3,s4)
        westflux = s2#-0.5 .* slopefit(s1,s2,s3)

        fluxImpllicit =
            KS.Q.pointsxyz[q, 1] ./ KS.dx .* (-1) + KS.Q.pointsxyz[q, 2] ./ KS.dy .* (-1)
        flux =
            KS.Q.pointsxyz[q, 1] ./ KS.dx .* (eastflux) +
            KS.Q.pointsxyz[q, 2] ./ KS.dy .* (northflux)
        psi[q, i, j] =
            KS.dt .* (-flux + rhs[q, i, j]) /
            (1 + KS.dt .* (KS.SigmaT[i, j] + fluxImpllicit))
    end
end

max_stars_repo_name: vavrines/CSN

max_issues_repo_name: vavrines/CSN

max_forks_repo_name: vavrines/CSN


text:
using PlotlyJS

function area1()
    trace1 = scatter(;x=1:4, y=[0, 2, 3, 5], fill="tozeroy")
    trace2 = scatter(;x=1:4, y=[3, 5, 1, 7], fill="tonexty")
    plot([trace1, trace2])
end

function area2()
    function _stacked_area!(traces)
        for (i, tr) in enumerate(traces[2:end])
            for j in 1:min(length(traces[i]["y"]), length(tr["y"]))
                tr["y"][j] += traces[i]["y"][j]
            end
        end
        traces
    end

    traces = [scatter(;x=1:3, y=[2, 1, 4], fill="tozeroy"),
              scatter(;x=1:3, y=[1, 1, 2], fill="tonexty"),
              scatter(;x=1:3, y=[3, 0, 2], fill="tonexty")]
    _stacked_area!(traces)

    plot(traces, Layout(title="stacked and filled line chart"))
end

function area3()
    trace1 = scatter(;x=1:4, y=[0, 2, 3, 5], fill="tozeroy", mode="none")
    trace2 = scatter(;x=1:4, y=[3, 5, 1, 7], fill="tonexty", mode="none")
    plot([trace1, trace2],
         Layout(title="Overlaid Chart Without Boundary Lines"))
end


text:
module Muffins

include("FC.jl")
using .FC

include("Half.jl")
using .Half

include("Int.jl")
using .Int

include("Mid.jl")
using .Mid

include("EBM.jl")
using .EBM

include("HBM.jl")
using .HBM

include("Gap.jl")
using .Gap

include("FindProc.jl")
using .FindProc

include("Matrix.jl")
using .Matrix

export muffins

# Solves muffin problem for m muffins and s students -- Work in progress
# Authors: Antara Hebbar and Jason Liu
function muffins(m::Int64, s::Int64; output::Int64 = 1)
    alphas = [
        fc(m, s, output = 0),
        half(m, s, output = 0),
        int(m, s, output = 0),
        mid(m, s, output = 0),
        ebm(m, s, output = 0),
        hbm(m, s, output = 0),
        gap(m, s, output = 0),
    ]
    alpha = minimum(alphas)
    disp = output > 0

    method = ""
    if alpha == alphas[1]
        fc(m, s, output = output)
        method = "Floor-Ceiling Theorem"
    elseif alpha == alphas[2]
        half(m, s, output = output)
        method = "Half Method"
    elseif alpha == alphas[3]
        int(m, s, output = output)
        method = "Interval Method"
    elseif alpha == alphas[4]
        mid(m, s, output = output)
        method = "Midpoint Method"
    elseif alpha == alphas[5]
        ebm(m, s, output = output)
        method = "Easy Buddy Match"
    elseif alpha == alphas[6]
        hbm(m, s, output = output)
        method = "Hard Buddy Match"
    else
        gap(m, s, output = output)
        method = "Gap Method"
    end
    disp && println("\nOptimal α derived by $method ⮥")

    procedures = findproc(m, s, alpha, output = 0)[2]
    if procedures == Nothing
        disp &&
            println("\nFindProc could not confirm optimal α to be a lower bound\nmuffins($m,$s) failed")
        for i in vcat(numerator(alpha)-1:-1:denominator(alpha)//3, denominator(alpha) // 3)
            lowerB = i // denominator(alpha)
            if findproc(m, s, lowerB, output = 0)[2] != Nothing
                findproc(m, s, lowerB, output = output)
                disp &&
                    println("\nFindProc found $(numerator(lowerB))/$(denominator(lowerB)) to be a lower bound for muffins($m,$s) ⮥")
                return lowerB
            end
        end
        return 1 // 3
    end

    findproc(m, s, alpha, output = output)
    disp && println("\nFindProc found optimal α to be a lower bound as well ⮥")
    disp && println("\nmuffins($m,$s) = $(numerator(alpha))/$(denominator(alpha))")
    alpha
end

end


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#alkashi-log.jl

include("../unum.jl")
using Unums

#implements an alkashi-like logarithm algorithm for Float64 as a testground for
#doing it in unums.

#returns a Int16 "1" if it's negative, "0" if it's positive.
function signof(x::FloatingPoint)
  T = typeof(x)
  if (T == Float64)
    return ((reinterpret(UInt64, x) & 0x8000_0000_0000_0000) != 0) ? 1 : 0
  elseif (T == Float32)
    return ((reinterpret(UInt32, x) & 0x8000_0000) != 0) ? 1 : 0
  end
end

function exponentof(x::FloatingPoint)
  T = typeof(x)
  if (T == Float64)
    return Int64((reinterpret(Int64, x) & 0x7FF0_0000_0000_0000) >> 52 - 1023)
  elseif (T == Float32)
    return Int64((reinterpret(Int32, x) & 0x7F80_0000) >> 23 - 127)
  end
end

function castfrac(x::FloatingPoint)
  T = typeof(x)
  if (T == Float64)
    return (reinterpret(UInt64, x) & 0x000F_FFFF_FFFF_FFFF) << 12
  elseif (T == Float32)
    return UInt64(reinterpret(UInt32, x) & 0x007F_FFFF) << 41
  end
end

function maskfor(T::Type)
  if (T == Float64)
    return 0xFFFF_FFFF_FFFF_FF00  #56 digits suffices for Float64 (52 digits)
  elseif (T == Float32)
    return 0xFFFF_FFE0_0000_0000  #27 digits of precision seems to suffice for a Float32 (23 digits)
  end
end

o16 = one(UInt16)
z16 = zero(UInt16)
__clz_array=[0x0004,0x0003,0x0002,0x0002, o16, o16, o16, o16, z16,z16,z16,z16,z16,z16,z16,z16]
function leading_zeros(n)
  (n == 0) && return 64
  res::UInt16 = 0
  #use the binary search method
  (n & 0xFFFF_FFFF_0000_0000 == 0) && (n <<= 32; res += 0x0020)
  (n & 0xFFFF_0000_0000_0000 == 0) && (n <<= 16; res += 0x0010)
  (n & 0xFF00_0000_0000_0000 == 0) && (n <<= 8;  res += 0x0008)
  (n & 0xF000_0000_0000_0000 == 0) && (n <<= 4;  res += 0x0004)
  res + __clz_array[(n >> 60) + 1]
end

#simple fused - multiply - add.  Assumes num1 has hidden bits "carry" and num2
#has no hidden bits.
rm = 0x0000_0000_FFFF_FFFF
function sfma(carry, num1, num2)
  (fracprod, _) = Unums.__chunk_mult(num1, num2)
  (_carry, fracprod) = Unums.__carried_add(carry, num1, fracprod)
  ((carry & 0x1) != 0) && ((_carry, fracprod) = Unums.__carried_add(_carry, num2, fracprod))
  ((carry & 0x2) != 0) && ((_carry, fracprod) = Unums.__carried_add(_carry, lsh(num2, 1), fracprod))
  (_carry, fracprod)
end

#performs a simple multiply, Assumes that number 1 has a hidden bit of exactly one
#and number 2 has a hidden bit of exactly zero
#(1 + a)(0 + b) = b + ab
function smult(a::UInt64, b::UInt64)

  (fraction, _) = Unums.__chunk_mult(a, b)
  carry = one(UInt64)

  #only perform the respective adds if the *opposing* thing is not subnormal.
  ((carry, fraction) = Unums.__carried_add(carry, fraction, b))

  #carry may be as high as three!  So we must shift as necessary.
  (fraction, shift, is_ubit) = Unums.__shift_after_add(carry, fraction, _)
  fraction << 1
end

function reassemble(sign::UInt64, ev::UInt64, fv::UInt64)
  number = (fv >> 12) | ((ev + 1023) << 52) | (sign << 63)

  println(bits(number))

  reinterpret(Float64, number)
end

const t32 = 0xFFFF_FFFF_0000_0000
const l32 = 0x0000_0000_FFFF_FFFF

function wsqr(a::UInt64)
  top = (a & t32) >> 32
  bot = (a & l32)

  carry = (a & 0x8000_0000_0000_0000 != 0) ? 2 : 1
  res = a << 1
  res2 = res + (top * top)
  (res2 < res) && (carry += 1)
  res3 = res2 + (bot * top) >> 31
  (res3 < res2) && (carry += 1)
  (carry, res3)
end

function exlg(x::FloatingPoint)
  #exact floating point with the goldschmidt algorithm.
  T = typeof(x)

  isnan(x) && return (nan(T), false)

  x <= 0 && return nan(T)

  #calculate the exponent.
  exp_f::Int64 = exponentof(x) + (issubnormal(x) ? 1 : 0)

  #figure the decimals.
  fraction::UInt64 = castfrac(x)

  if (issubnormal(x))
<<<<<<< HEAD:oddsandends/alkashi-log.jl
    shift::UInt64 = clz(fraction) + 1
=======
    shift::UInt64 = leading_zeros(fraction) + 1
>>>>>>> directory-reorg:src/oddsandends/alkashi-log.jl
    fraction = fraction << shift
    exp_f -= shift
  end

  sign::UInt64 = 0
  if (exp_f < 0)
    exp_f = -exp_f - 1
    sign = 1
  end
<<<<<<< HEAD:oddsandends/alkashi-log.jl
  lz = clz(UInt64(exp_f))
=======
  lz = leading_zeros(UInt64(exp_f))
>>>>>>> directory-reorg:src/oddsandends/alkashi-log.jl
  resexp = 63 - lz
  #add the exponent part onto the result fraction.
  resfrac = UInt64(exp_f << (lz + 1))

  bitsofar = zero(UInt64)
  #do the alkashi-type algorithm.
  for idx = 1:52
    (c, q) = wsqr(fraction)

    if (c == 1)
      #do nothing.  Merely set x to y.
      fraction = q
    else
      bitsofar |= (0x8000_0000_0000_0000 >> (idx - 1))
      fraction = q >> 1
      (c == 3) && (fraction |= (0x8000_0000_0000_0000))
    end
  end

  iint::UInt64 = (sign == 0) ? bitsofar : -bitsofar
  resfrac |= iint >> resexp
  #resexp = 0
  #resfrac = 0
  #reassemble the value into the requisite floating point
  reassemble(sign, resexp, resfrac)
end

#one-time testing
#v = (rand(Int64) & 0xFFF0_0000_0000_0000) | (rand(Int64) & 0x0000_0000_0FFF_FFFF)
v = rand(Int64)
x = abs(reinterpret(Float64, v))
z = exlg(x)
a = log2(x)

println("input :    $x")
println("answer:    $(log2(x))")
println("calculate: $z")

println("xbits:     ", bits(x))
println("abits:     ", bits(a))
println("zbits:     ", bits(z))

println("diff:      ", reinterpret(UInt64, a) - reinterpret(UInt64, z))

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using Random
using LinearAlgebra

export LTMADS

"""
    LTMADS()

Return an empty LTMADS object.

LTMADS is a poll stage that creates a set of directions based
on a semi-randomly generated lower triangular matrix. This randomness
means that several runs of the algorithm may be needed to find a minimum.
"""
mutable struct LTMADS{T} <: AbstractPoll
    b::Dict{T,Vector{T}}
    i::Dict{T,Int}
    maximal_basis::Bool
    LTMADS(;kwargs...) = LTMADS{Float64}(;kwargs...)
    function LTMADS{T}(;maximal_basis=true) where T
        g = new()
        g.b = Dict{T, Vector{T}}()
        g.i = Dict{T, Int}()
        g.maximal_basis=maximal_basis
        return g
    end
end

"""
    GenerateDirections(p::DSProblem{T}, DG::LTMADS{T})::Vector{Vector{T}}

Generates columns and forms a basis matrix for direction generation.
"""
function GenerateDirections(p::AbstractProblem, DG::LTMADS{T})::Matrix{T} where T
    B = LT_basis_generation(p.mesh, p.N, DG)
    Dₖ = form_basis_matrix(p.N, B, DG.maximal_basis)

    return Dₖ
end

function form_basis_matrix(N::Int, B::Matrix{T}, max_basis::Bool) where T
    max_basis && return [B -B]

    d = zeros(T, N)
    for (i,_) in enumerate(d)
        d[i] = -sum(B[i,:])
    end

    return [B d]
end

function LT_basis_generation(m::Mesh, N::Int, DG::LTMADS{T}) where T
    b, i = b_l_generation(DG.b, DG.i, m.l, N)

    L = L_generation(N, m.l)

    B = B_generation(N, i, b, L)

    B′ = B′_generation(B, N)

    return B′
end

function B′_generation(B, N; perm=shuffle(1:N))
    B′ = zeros(N,N)
    for (i,e) in enumerate(eachcol(B))
        B′[:,perm[i]] = e
    end
    return B′
end

function b_l_generation(b::Dict{T,Vector{T}}, i::Dict{T,Int}, l::Int, N::Int
                       ) where T
    if !haskey(b, l)
        i[l] = rand(1:N)
        b[l] = zeros(T, N)

        for j in 1:N
            if j == i[l]
                b[l][j] = rand([-2^l, 2^l])
            else
                b[l][j] = rand(-2^l+1:2^l-1)
            end
        end
    end
    return b[l], i[l]
end

function L_generation(N, l)
    L = zeros(N-1,N-1)

    for i=1:N-1, j=1:N-1
        if j==i
            L[i,j] = rand([2^l, -2^l])
        elseif j < i
            L[i,j] = rand(1-2^l:-1+2^l)
        end
    end

    return L
end

function B_generation(N, i, b, L; perm=shuffle(setdiff(1:N, i)))
    B = zeros(N,N-1)
    for (i,e) in enumerate(eachrow(L))
        B[perm[i],:] = e
    end
    B = [B b]
    return B
end


text:
using Mocking

create_listener_1_patch = @patch function create_listener(a...; b...)
    return Dict{String, Any}("Listeners" => Any[Dict{String, Any}("ListenerArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:listener/app/my-load-balancer/50dc6c495c0c9188/f2f7dc8efc522ab2", "DefaultActions" => Any[Dict{String, Any}("Type" => "forward", "TargetGroupArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:targetgroup/my-targets/73e2d6bc24d8a067")], "Protocol" => "HTTP", "LoadBalancerArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:loadbalancer/app/my-load-balancer/50dc6c495c0c9188", "Port" => 80)])
end

create_listener_2_patch = @patch function create_listener(a...; b...)
    return Dict{String, Any}("Listeners" => Any[Dict{String, Any}("ListenerArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:listener/app/my-load-balancer/50dc6c495c0c9188/f2f7dc8efc522ab2", "SslPolicy" => "ELBSecurityPolicy-2015-05", "Certificates" => Any[Dict{String, Any}("CertificateArn" => "arn:aws:iam::123456789012:server-certificate/my-server-cert")], "DefaultActions" => Any[Dict{String, Any}("Type" => "forward", "TargetGroupArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:targetgroup/my-targets/73e2d6bc24d8a067")], "Protocol" => "HTTPS", "LoadBalancerArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:loadbalancer/app/my-load-balancer/50dc6c495c0c9188", "Port" => 443)])
end

create_load_balancer_1_patch = @patch function create_load_balancer(a...; b...)
    return Dict{String, Any}("LoadBalancers" => Any[Dict{String, Any}("AvailabilityZones" => Any[Dict{String, Any}("ZoneName" => "us-west-2a", "SubnetId" => "subnet-8360a9e7"), Dict{String, Any}("ZoneName" => "us-west-2b", "SubnetId" => "subnet-b7d581c0")], "CreatedTime" => "2016-03-25T21:26:12.920Z", "SecurityGroups" => Any["sg-5943793c"], "VpcId" => "vpc-3ac0fb5f", "State" => Dict{String, Any}("Code" => "provisioning"), "Scheme" => "internet-facing", "CanonicalHostedZoneId" => "Z2P70J7EXAMPLE", "DNSName" => "my-load-balancer-424835706.us-west-2.elb.amazonaws.com", "LoadBalancerName" => "my-load-balancer", "Type" => "application", "LoadBalancerArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:loadbalancer/app/my-load-balancer/50dc6c495c0c9188")])
end

create_load_balancer_2_patch = @patch function create_load_balancer(a...; b...)
    return Dict{String, Any}("LoadBalancers" => Any[Dict{String, Any}("AvailabilityZones" => Any[Dict{String, Any}("ZoneName" => "us-west-2a", "SubnetId" => "subnet-8360a9e7"), Dict{String, Any}("ZoneName" => "us-west-2b", "SubnetId" => "subnet-b7d581c0")], "CreatedTime" => "2016-03-25T21:29:48.850Z", "SecurityGroups" => Any["sg-5943793c"], "VpcId" => "vpc-3ac0fb5f", "State" => Dict{String, Any}("Code" => "provisioning"), "Scheme" => "internal", "CanonicalHostedZoneId" => "Z2P70J7EXAMPLE", "DNSName" => "internal-my-internal-load-balancer-1529930873.us-west-2.elb.amazonaws.com", "LoadBalancerName" => "my-internal-load-balancer", "Type" => "application", "LoadBalancerArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:loadbalancer/app/my-internal-load-balancer/5b49b8d4303115c2")])
end

create_rule_1_patch = @patch function create_rule(a...; b...)
    return Dict{String, Any}("Rules" => Any[Dict{String, Any}("Actions" => Any[Dict{String, Any}("Type" => "forward", "TargetGroupArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:targetgroup/my-targets/73e2d6bc24d8a067")], "Conditions" => Any[Dict{String, Any}("Values" => Any["/img/*"], "Field" => "path-pattern")], "IsDefault" => false, "Priority" => "10", "RuleArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:listener-rule/app/my-load-balancer/50dc6c495c0c9188/f2f7dc8efc522ab2/9683b2d02a6cabee")])
end

create_target_group_1_patch = @patch function create_target_group(a...; b...)
    return Dict{String, Any}("TargetGroups" => Any[Dict{String, Any}("HealthCheckProtocol" => "HTTP", "UnhealthyThresholdCount" => 2, "TargetGroupArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:targetgroup/my-targets/73e2d6bc24d8a067", "Protocol" => "HTTP", "TargetGroupName" => "my-targets", "HealthCheckPath" => "/", "HealthCheckTimeoutSeconds" => 5, "HealthyThresholdCount" => 5, "VpcId" => "vpc-3ac0fb5f", "Matcher" => Dict{String, Any}("HttpCode" => "200"), "HealthCheckPort" => "traffic-port", "Port" => 80, "HealthCheckIntervalSeconds" => 30)])
end

describe_listeners_1_patch = @patch function describe_listeners(a...; b...)
    return Dict{String, Any}("Listeners" => Any[Dict{String, Any}("ListenerArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:listener/app/my-load-balancer/50dc6c495c0c9188/f2f7dc8efc522ab2", "DefaultActions" => Any[Dict{String, Any}("Type" => "forward", "TargetGroupArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:targetgroup/my-targets/73e2d6bc24d8a067")], "Protocol" => "HTTP", "LoadBalancerArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:loadbalancer/app/my-load-balancer/50dc6c495c0c9188", "Port" => 80)])
end

describe_load_balancer_attributes_1_patch = @patch function describe_load_balancer_attributes(a...; b...)
    return Dict{String, Any}("Attributes" => Any[Dict{String, Any}("Value" => "false", "Key" => "access_logs.s3.enabled"), Dict{String, Any}("Value" => "60", "Key" => "idle_timeout.timeout_seconds"), Dict{String, Any}("Value" => "", "Key" => "access_logs.s3.prefix"), Dict{String, Any}("Value" => "false", "Key" => "deletion_protection.enabled"), Dict{String, Any}("Value" => "", "Key" => "access_logs.s3.bucket")])
end

describe_load_balancers_1_patch = @patch function describe_load_balancers(a...; b...)
    return Dict{String, Any}("LoadBalancers" => Any[Dict{String, Any}("AvailabilityZones" => Any[Dict{String, Any}("ZoneName" => "us-west-2a", "SubnetId" => "subnet-8360a9e7"), Dict{String, Any}("ZoneName" => "us-west-2b", "SubnetId" => "subnet-b7d581c0")], "CreatedTime" => "2016-03-25T21:26:12.920Z", "SecurityGroups" => Any["sg-5943793c"], "VpcId" => "vpc-3ac0fb5f", "State" => Dict{String, Any}("Code" => "active"), "Scheme" => "internet-facing", "CanonicalHostedZoneId" => "Z2P70J7EXAMPLE", "DNSName" => "my-load-balancer-424835706.us-west-2.elb.amazonaws.com", "LoadBalancerName" => "my-load-balancer", "Type" => "application", "LoadBalancerArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:loadbalancer/app/my-load-balancer/50dc6c495c0c9188")])
end

describe_rules_1_patch = @patch function describe_rules(a...; b...)
    return Dict{String, Any}("Rules" => Any[Dict{String, Any}("Actions" => Any[Dict{String, Any}("Type" => "forward", "TargetGroupArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:targetgroup/my-targets/73e2d6bc24d8a067")], "Conditions" => Any[Dict{String, Any}("Values" => Any["/img/*"], "Field" => "path-pattern")], "IsDefault" => false, "Priority" => "10", "RuleArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:listener-rule/app/my-load-balancer/50dc6c495c0c9188/f2f7dc8efc522ab2/9683b2d02a6cabee")])
end

describe_sslpolicies_1_patch = @patch function describe_sslpolicies(a...; b...)
    return Dict{String, Any}("SslPolicies" => Any[Dict{String, Any}("Ciphers" => Any[Dict{String, Any}("Priority" => 1, "Name" => "ECDHE-ECDSA-AES128-GCM-SHA256"), Dict{String, Any}("Priority" => 2, "Name" => "ECDHE-RSA-AES128-GCM-SHA256"), Dict{String, Any}("Priority" => 3, "Name" => "ECDHE-ECDSA-AES128-SHA256"), Dict{String, Any}("Priority" => 4, "Name" => "ECDHE-RSA-AES128-SHA256"), Dict{String, Any}("Priority" => 5, "Name" => "ECDHE-ECDSA-AES128-SHA"), Dict{String, Any}("Priority" => 6, "Name" => "ECDHE-RSA-AES128-SHA"), Dict{String, Any}("Priority" => 7, "Name" => "DHE-RSA-AES128-SHA"), Dict{String, Any}("Priority" => 8, "Name" => "ECDHE-ECDSA-AES256-GCM-SHA384"), Dict{String, Any}("Priority" => 9, "Name" => "ECDHE-RSA-AES256-GCM-SHA384"), Dict{String, Any}("Priority" => 10, "Name" => "ECDHE-ECDSA-AES256-SHA384"), Dict{String, Any}("Priority" => 11, "Name" => "ECDHE-RSA-AES256-SHA384"), Dict{String, Any}("Priority" => 12, "Name" => "ECDHE-RSA-AES256-SHA"), Dict{String, Any}("Priority" => 13, "Name" => "ECDHE-ECDSA-AES256-SHA"), Dict{String, Any}("Priority" => 14, "Name" => "AES128-GCM-SHA256"), Dict{String, Any}("Priority" => 15, "Name" => "AES128-SHA256"), Dict{String, Any}("Priority" => 16, "Name" => "AES128-SHA"), Dict{String, Any}("Priority" => 17, "Name" => "AES256-GCM-SHA384"), Dict{String, Any}("Priority" => 18, "Name" => "AES256-SHA256"), Dict{String, Any}("Priority" => 19, "Name" => "AES256-SHA")], "SslProtocols" => Any["TLSv1", "TLSv1.1", "TLSv1.2"], "Name" => "ELBSecurityPolicy-2015-05")])
end

describe_tags_1_patch = @patch function describe_tags(a...; b...)
    return Dict{String, Any}("TagDescriptions" => Any[Dict{String, Any}("ResourceArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:loadbalancer/app/my-load-balancer/50dc6c495c0c9188", "Tags" => Any[Dict{String, Any}("Value" => "lima", "Key" => "project"), Dict{String, Any}("Value" => "digital-media", "Key" => "department")])])
end

describe_target_group_attributes_1_patch = @patch function describe_target_group_attributes(a...; b...)
    return Dict{String, Any}("Attributes" => Any[Dict{String, Any}("Value" => "false", "Key" => "stickiness.enabled"), Dict{String, Any}("Value" => "300", "Key" => "deregistration_delay.timeout_seconds"), Dict{String, Any}("Value" => "lb_cookie", "Key" => "stickiness.type"), Dict{String, Any}("Value" => "86400", "Key" => "stickiness.lb_cookie.duration_seconds")])
end

describe_target_groups_1_patch = @patch function describe_target_groups(a...; b...)
    return Dict{String, Any}("TargetGroups" => Any[Dict{String, Any}("HealthCheckProtocol" => "HTTP", "UnhealthyThresholdCount" => 2, "TargetGroupArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:targetgroup/my-targets/73e2d6bc24d8a067", "Protocol" => "HTTP", "TargetGroupName" => "my-targets", "HealthCheckPath" => "/", "HealthCheckTimeoutSeconds" => 5, "HealthyThresholdCount" => 5, "VpcId" => "vpc-3ac0fb5f", "Matcher" => Dict{String, Any}("HttpCode" => "200"), "HealthCheckPort" => "traffic-port", "Port" => 80, "LoadBalancerArns" => Any["arn:aws:elasticloadbalancing:us-west-2:123456789012:loadbalancer/app/my-load-balancer/50dc6c495c0c9188"], "HealthCheckIntervalSeconds" => 30)])
end

describe_target_health_1_patch = @patch function describe_target_health(a...; b...)
    return Dict{String, Any}("TargetHealthDescriptions" => Any[Dict{String, Any}("TargetHealth" => Dict{String, Any}("Description" => "Given target group is not configured to receive traffic from ELB", "Reason" => "Target.NotInUse", "State" => "unused"), "Target" => Dict{String, Any}("Id" => "i-0f76fade", "Port" => 80)), Dict{String, Any}("TargetHealth" => Dict{String, Any}("State" => "healthy"), "HealthCheckPort" => "80", "Target" => Dict{String, Any}("Id" => "i-0f76fade", "Port" => 80))])
end

describe_target_health_2_patch = @patch function describe_target_health(a...; b...)
    return Dict{String, Any}("TargetHealthDescriptions" => Any[Dict{String, Any}("TargetHealth" => Dict{String, Any}("State" => "healthy"), "HealthCheckPort" => "80", "Target" => Dict{String, Any}("Id" => "i-0f76fade", "Port" => 80))])
end

modify_listener_1_patch = @patch function modify_listener(a...; b...)
    return Dict{String, Any}("Listeners" => Any[Dict{String, Any}("ListenerArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:listener/app/my-load-balancer/50dc6c495c0c9188/f2f7dc8efc522ab2", "DefaultActions" => Any[Dict{String, Any}("Type" => "forward", "TargetGroupArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:targetgroup/my-new-targets/2453ed029918f21f")], "Protocol" => "HTTP", "LoadBalancerArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:loadbalancer/app/my-load-balancer/50dc6c495c0c9188", "Port" => 80)])
end

modify_listener_2_patch = @patch function modify_listener(a...; b...)
    return Dict{String, Any}("Listeners" => Any[Dict{String, Any}("ListenerArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:listener/app/my-load-balancer/50dc6c495c0c9188/0467ef3c8400ae65", "SslPolicy" => "ELBSecurityPolicy-2015-05", "Certificates" => Any[Dict{String, Any}("CertificateArn" => "arn:aws:iam::123456789012:server-certificate/my-new-server-cert")], "DefaultActions" => Any[Dict{String, Any}("Type" => "forward", "TargetGroupArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:targetgroup/my-targets/73e2d6bc24d8a067")], "Protocol" => "HTTPS", "LoadBalancerArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:loadbalancer/app/my-load-balancer/50dc6c495c0c9188", "Port" => 443)])
end

modify_load_balancer_attributes_1_patch = @patch function modify_load_balancer_attributes(a...; b...)
    return Dict{String, Any}("Attributes" => Any[Dict{String, Any}("Value" => "true", "Key" => "deletion_protection.enabled"), Dict{String, Any}("Value" => "false", "Key" => "access_logs.s3.enabled"), Dict{String, Any}("Value" => "60", "Key" => "idle_timeout.timeout_seconds"), Dict{String, Any}("Value" => "", "Key" => "access_logs.s3.prefix"), Dict{String, Any}("Value" => "", "Key" => "access_logs.s3.bucket")])
end

modify_load_balancer_attributes_2_patch = @patch function modify_load_balancer_attributes(a...; b...)
    return Dict{String, Any}("Attributes" => Any[Dict{String, Any}("Value" => "30", "Key" => "idle_timeout.timeout_seconds"), Dict{String, Any}("Value" => "false", "Key" => "access_logs.s3.enabled"), Dict{String, Any}("Value" => "", "Key" => "access_logs.s3.prefix"), Dict{String, Any}("Value" => "true", "Key" => "deletion_protection.enabled"), Dict{String, Any}("Value" => "", "Key" => "access_logs.s3.bucket")])
end

modify_load_balancer_attributes_3_patch = @patch function modify_load_balancer_attributes(a...; b...)
    return Dict{String, Any}("Attributes" => Any[Dict{String, Any}("Value" => "true", "Key" => "access_logs.s3.enabled"), Dict{String, Any}("Value" => "my-load-balancer-logs", "Key" => "access_logs.s3.bucket"), Dict{String, Any}("Value" => "myapp", "Key" => "access_logs.s3.prefix"), Dict{String, Any}("Value" => "60", "Key" => "idle_timeout.timeout_seconds"), Dict{String, Any}("Value" => "false", "Key" => "deletion_protection.enabled")])
end

modify_rule_1_patch = @patch function modify_rule(a...; b...)
    return Dict{String, Any}("Rules" => Any[Dict{String, Any}("Actions" => Any[Dict{String, Any}("Type" => "forward", "TargetGroupArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:targetgroup/my-targets/73e2d6bc24d8a067")], "Conditions" => Any[Dict{String, Any}("Values" => Any["/images/*"], "Field" => "path-pattern")], "IsDefault" => false, "Priority" => "10", "RuleArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:listener-rule/app/my-load-balancer/50dc6c495c0c9188/f2f7dc8efc522ab2/9683b2d02a6cabee")])
end

modify_target_group_1_patch = @patch function modify_target_group(a...; b...)
    return Dict{String, Any}("TargetGroups" => Any[Dict{String, Any}("HealthCheckProtocol" => "HTTPS", "UnhealthyThresholdCount" => 2, "TargetGroupArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:targetgroup/my-https-targets/2453ed029918f21f", "Protocol" => "HTTPS", "TargetGroupName" => "my-https-targets", "HealthCheckTimeoutSeconds" => 5, "HealthyThresholdCount" => 5, "VpcId" => "vpc-3ac0fb5f", "Matcher" => Dict{String, Any}("HttpCode" => "200"), "HealthCheckPort" => "443", "Port" => 443, "LoadBalancerArns" => Any["arn:aws:elasticloadbalancing:us-west-2:123456789012:loadbalancer/app/my-load-balancer/50dc6c495c0c9188"], "HealthCheckIntervalSeconds" => 30)])
end

modify_target_group_attributes_1_patch = @patch function modify_target_group_attributes(a...; b...)
    return Dict{String, Any}("Attributes" => Any[Dict{String, Any}("Value" => "false", "Key" => "stickiness.enabled"), Dict{String, Any}("Value" => "600", "Key" => "deregistration_delay.timeout_seconds"), Dict{String, Any}("Value" => "lb_cookie", "Key" => "stickiness.type"), Dict{String, Any}("Value" => "86400", "Key" => "stickiness.lb_cookie.duration_seconds")])
end

set_rule_priorities_1_patch = @patch function set_rule_priorities(a...; b...)
    return Dict{String, Any}("Rules" => Any[Dict{String, Any}("Actions" => Any[Dict{String, Any}("Type" => "forward", "TargetGroupArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:targetgroup/my-targets/73e2d6bc24d8a067")], "Conditions" => Any[Dict{String, Any}("Values" => Any["/img/*"], "Field" => "path-pattern")], "IsDefault" => false, "Priority" => "5", "RuleArn" => "arn:aws:elasticloadbalancing:us-west-2:123456789012:listener-rule/app/my-load-balancer/50dc6c495c0c9188/f2f7dc8efc522ab2/1291d13826f405c3")])
end

set_security_groups_1_patch = @patch function set_security_groups(a...; b...)
    return Dict{String, Any}("SecurityGroupIds" => Any["sg-5943793c"])
end

set_subnets_1_patch = @patch function set_subnets(a...; b...)
    return Dict{String, Any}("AvailabilityZones" => Any[Dict{String, Any}("ZoneName" => "us-west-2a", "SubnetId" => "subnet-8360a9e7"), Dict{String, Any}("ZoneName" => "us-west-2b", "SubnetId" => "subnet-b7d581c0")])
end


text:
# Autogenerated wrapper script for Arb_jll for aarch64-linux-gnu
export libarb

using FLINT_jll
using GMP_jll
using MPFR_jll
JLLWrappers.@generate_wrapper_header("Arb")
JLLWrappers.@declare_library_product(libarb, "libarb.so.2")
function __init__()
    JLLWrappers.@generate_init_header(FLINT_jll, GMP_jll, MPFR_jll)
    JLLWrappers.@init_library_product(
        libarb,
        "lib/libarb.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
### CONCRETE TYPE: DIRECT PROX EVALUATION
# prox! is computed using a Cholesky factorization of Q + I/gamma.
# The factorization is cached and recomputed whenever gamma changes

using LinearAlgebra
using SparseArrays
using SuiteSparse

mutable struct QuadraticDirect{R <: Real, M <: AbstractMatrix{R}, V <: AbstractVector{R}, F <: Factorization} <: Quadratic
    Q::M
    q::V
    gamma::R
    temp::V
    fact::F
    function QuadraticDirect{R, M, V, F}(Q::M, q::V) where {R <: Real, M <: AbstractMatrix{R}, V <: AbstractVector{R}, F <: Factorization}
        if size(Q, 1) != size(Q, 2) || length(q) != size(Q, 2)
            error("Q must be squared and q must be compatible with Q")
        end
        new(Q, q, -1, similar(q))
    end
end

function QuadraticDirect(Q::M, q::V) where {R <: Real, I <: Integer, M <: SparseMatrixCSC{R, I}, V <: AbstractVector{R}}
    QuadraticDirect{R, M, V, SuiteSparse.CHOLMOD.Factor{R}}(Q, q)
end

function QuadraticDirect(Q::M, q::V) where {R <: Real, M <: DenseMatrix{R}, V <: AbstractVector{R}}
    QuadraticDirect{R, M, V, Cholesky{R, M}}(Q, q)
end

function (f::QuadraticDirect{R, M, V, F})(x::AbstractArray{R}) where {R, M, V, F}
    mul!(f.temp, f.Q, x)
    return 0.5*dot(x, f.temp) + dot(x, f.q)
end

function prox!(y::AbstractArray{R}, f::QuadraticDirect{R, M, V, F}, x::AbstractArray{R}, gamma::R=R(1)) where {R, M, V, F <: Cholesky}
    if gamma != f.gamma
        factor_step!(f, gamma)
    end
    y .= x./gamma
    y .-= f.q
    # Qy = U'Uy = b, therefore y = U\(U'\b)
    LAPACK.trtrs!('U', 'C', 'N', f.fact.factors, y)
    LAPACK.trtrs!('U', 'N', 'N', f.fact.factors, y)
    mul!(f.temp, f.Q, y)
    fy = 0.5*dot(y, f.temp) + dot(y, f.q)
    return fy
end

function prox!(y::AbstractArray{R}, f::QuadraticDirect{R, M, V, F}, x::AbstractArray{R}, gamma::R=R(1)) where {R, M, V, F <: SuiteSparse.CHOLMOD.Factor}
    if gamma != f.gamma
        factor_step!(f, gamma)
    end
    f.temp .= x./gamma
    f.temp .-= f.q
    y .= f.fact\f.temp
    mul!(f.temp, f.Q, y)
    fy = 0.5*dot(y, f.temp) + dot(y, f.q)
    return fy
end

function factor_step!(f::QuadraticDirect{R, M, V, F}, gamma::R) where {R, M <: DenseMatrix{R}, V, F}
    f.gamma = gamma
    f.fact = cholesky(f.Q + I/gamma)
end

function factor_step!(f::QuadraticDirect{R, M, V, F}, gamma::R) where {R, I, M <: SparseMatrixCSC{R, I}, V, F}
    f.gamma = gamma
    f.fact = ldlt(f.Q; shift = 1/gamma)
end

function gradient!(y::AbstractArray{R}, f::QuadraticDirect{R, M, V, F}, x::AbstractArray{R}) where {R, M, V, F}
    mul!(y, f.Q, x)
    y .+= f.q
    return 0.5*(dot(x, y) + dot(x, f.q))
end

function prox_naive(f::QuadraticDirect, x::AbstractArray{R}, gamma::R=one(R)) where R
    y = (gamma*f.Q + I)\(x - gamma*f.q)
    fy = 0.5*dot(y, f.Q*y) + dot(y, f.q)
    return y, fy
end