text:
function log_one_node(com, nvars, back_idx, step_nr)
    tree_log_node = TreeLogNode()
    tree_log_node.step_nr = step_nr
    tree_log_node.id = back_idx
    parent_idx = 0
    if length(com.backtrack_vec) > 0
        tree_log_node.status = com.backtrack_vec[back_idx].status
        tree_log_node.var_idx = com.backtrack_vec[back_idx].variable_idx
        tree_log_node.set_val = com.backtrack_vec[back_idx].pval
        tree_log_node.best_bound = com.backtrack_vec[back_idx].best_bound

        parent_idx = com.backtrack_vec[back_idx].parent_idx
    else # for initial solve
        tree_log_node.status = :Closed
        tree_log_node.var_idx = 0
        tree_log_node.set_val = 0
        tree_log_node.best_bound = com.best_bound
    end
    tree_log_node.var_changes = Dict{Int64,Vector{Tuple{Symbol, Int64, Int64, Int64}}}()
    tree_log_node.var_states = Dict{Int64,Vector{Int64}}()
    if tree_log_node.status == :Closed
        for var in com.search_space
            if length(var.changes[back_idx]) > 0
                tree_log_node.var_states[var.idx] = sort!(values(var))
                tree_log_node.var_changes[var.idx] = var.changes[back_idx]
            end
        end
    end
    tree_log_node.children = Vector{TreeLogNode}()
    if parent_idx > 0
        changed = false
        for (i,child) in enumerate(com.logs[parent_idx].children)
            if child.id == back_idx
                com.logs[parent_idx].children[i] = tree_log_node
                changed = true
                break
            end
        end
        if !changed
            push!(com.logs[parent_idx].children, tree_log_node)
        end
    end
    return tree_log_node
end

function bfs_list(start_node::CS.TreeLogNode)
    to_process = Vector{CS.TreeLogNode}()
    depths = Vector{Int}()
    nodes_list = Vector{CS.TreeLogNode}()
    num_children = Vector{Int}()

    push!(to_process, start_node)
    push!(depths, 0)
    node = deepcopy(start_node)
    push!(num_children, length(node.children))
    node.children = TreeLogNode[]
    push!(nodes_list, node)

    while !isempty(to_process)
        current_node = popfirst!(to_process)
        depth = popfirst!(depths)

        node = deepcopy(current_node)
        push!(num_children, length(node.children))
        node.children = TreeLogNode[]
        push!(nodes_list, node)

        for child_node in current_node.children
            push!(to_process, child_node)
            push!(depths, depth + 1)
        end
    end
    return nodes_list, num_children
end

function same_logs(log1, log2)
    nodes_list1, num_children1 = bfs_list(log1)
    nodes_list2, num_children2 = bfs_list(log2)
    if length(nodes_list1) != length(nodes_list2)
        println("Different length")
        return false
    end
    for i=1:length(nodes_list1)
        node1 = nodes_list1[i]
        node2 = nodes_list2[i]
        num_child1 = num_children1[i]
        num_child2 = num_children2[i]
        if num_child1 != num_child2
            println("Different number of children at: ", i)
            return false
        end

        if node1.id != node2.id || node1.status != node2.status || node1.best_bound != node2.best_bound ||
            node1.step_nr != node2.step_nr || node1.var_idx != node2.var_idx || node1.set_val != node2.set_val || node1.var_states != node2.var_states
            println("Not identical at i=", i)
            println("node1: ")
            println("id: ", node1.id, " status: ", node1.status, " best_bound: ", node1.best_bound, " step_nr: ", node1.step_nr)
            println("var_idx: ", node1.var_idx, " set_val: ", node1.set_val, " var_states: ", node1.var_states)
            println("node2: ")
            println("id: ", node2.id, " status: ", node2.status, " best_bound: ", node2.best_bound, " step_nr: ", node2.step_nr)
            println("var_idx: ", node2.var_idx, " set_val: ", node2.set_val, " var_states: ", node2.var_states)
            return false
        end
    end
    return true
end

function get_logs(com::CS.CoM)
    logs = Dict{Symbol, Any}()
    nvars = length(com.search_space)

    logs[:init] = Vector{Vector{Int64}}(undef, nvars)
    for var in com.init_search_space
        logs[:init][var.idx] = values(var)
    end

    logs[:tree] = com.logs[1]
    return logs
end

"""
    save_logs(com::CS.CoM, filepath)

Save the tree structure and some additional problem information in a json file `filepath`.
Can be only used if `keep_logs` is set to `true` in the [`solve!`](@ref) call.
"""
function save_logs(com::CS.CoM, filepath)
    logs = get_logs(com)
    write(filepath, JSON.json(logs))
end


text:
# Test layers and data/model movements on and off the GPU
# Add tests for layers and their gradients on the GPU
# Most of the forward passes should be fine being applied
# to bitstype objects, but this gives higher coverage for our use-cases
# Check that getting the gradients does not throw

# generic movement tests
@testset "Basic GPU Movement" begin
  @test gradient(x -> sum(gpu(x)), rand(3,3)) isa Tuple
  @test gradient(x -> sum(cpu(x)), gpu(rand(3,3))) isa Tuple
end

# TODO: These layers get into scalar indexing
# `AlphaDropout` throws a compilation error on GPUs,
# whereas, the rest are scalar indexing issues.
const BROKEN_LAYERS = [DepthwiseConv,
		                   AlphaDropout,
                       InstanceNorm,
                       GroupNorm]

function gradtest(name::String, layers::Vector, xs = nothing, args...)
  isnothing(xs) && error("Missing input to test the layers against.")
  @testset "$name GPU grad tests" begin
    for layer in layers
      @testset "$layer GPU grad test" begin
        l = gpu(layer(args...))
        xs = gpu(xs)
        if any(x -> isa(l, x), BROKEN_LAYERS)
          ps = Flux.params(l)
          @test_broken gradient(() -> sum(l(xs)), ps) isa Flux.Zygote.Grads
        else
          ps = Flux.params(l)
          @test gradient(() -> sum(l(xs)), ps) isa Flux.Zygote.Grads
          gs = gradient(() -> sum(l(xs)), ps)

          # Handle pooling layers
          if !isempty(ps)
            @test gs[first(ps)] isa Flux.CUDA.CuArray
          end
        end
      end
    end
  end
end

# Repeats from Conv, CrossCor

r = rand(Float32, 28, 28, 1, 1)
conv_layers = [Conv, ConvTranspose, CrossCor, DepthwiseConv]
gradtest("Conv", conv_layers, r, (2,2), 1=>3)

pooling_layers = [MaxPool, MeanPool]
gradtest("Pooling", pooling_layers, r, (2,2))

adaptive_pooling_layers = [AdaptiveMaxPool, AdaptiveMeanPool]
gradtest("AdaptivePooling", adaptive_pooling_layers, r, (7,7))

dropout_layers = [Dropout, AlphaDropout]
gradtest("Dropout", dropout_layers, r, 0.5f0)

norm_layers = [LayerNorm, BatchNorm]
gradtest("Normalising", norm_layers, rand(Float32, 28,28,3,1), 1)

instancenorm = [InstanceNorm]
gradtest("InstanceNorm", instancenorm, r, 1)

groupnorm = [GroupNorm]
gradtest("GroupNorm", groupnorm, rand(Float32, 28,28,3,1), 3, 1)

const stateless_layers = [Flux.normalise]

const stateless_layers_broadcasted = []

function stateless_gradtest(f, args...)
  @test gradient((args...) -> sum(f(args...)), args...)[1] isa CuArray
end

function stateless_gradtest_broadcasted(f, args...)
  @test gradient((args...) -> sum(f.(args...)), args...)[1] isa CuArray
end

@testset "Stateless GPU grad tests" begin
  x = gpu(rand(3,3))
  y = gpu(rand(3,3))

  for layer in stateless_layers
    if layer == Flux.normalise
      stateless_gradtest(x -> layer(x, dims=1), x)
    else
      stateless_gradtest(layer, x, y)
    end
  end

  for layer in stateless_layers_broadcasted
    stateless_gradtest_broadcasted(layer, x, y)
  end
end

ext: jl


text:
module PkgUtils

import FileWatching
import Pkg
import ..Pluto
import ..Pluto: Notebook, save_notebook, load_notebook, load_notebook_nobackup, withtoken, Token, readwrite, PkgCompat
import ..Pluto.PkgCompat: project_file, manifest_file

using Markdown

export activate_notebook

ensure_has_nbpkg(notebook::Notebook) = if notebook.nbpkg_ctx === nothing

    # TODO: update_save the notebook to init packages and stuff?
    error("""
    This notebook is not using Pluto's package manager. This means that either:
    1. The notebook contains Pkg.activate or Pkg.add calls, or
    2. The notebook was created before Pluto 0.15.

    Open the notebook using Pluto to get started.
    """)
else
    for f in [notebook |> project_file, notebook |> manifest_file]
        isfile(f) || touch(f)
    end
end

function assert_has_manifest(dir::String)
    @assert isdir(dir)
    if !isfile(dir |> project_file)
        error("The given directory does not contain a Project.toml file -- it is not a package environment.")
    end
    if !isfile(dir |> manifest_file)
        error("The given directory does not contain a Manifest.toml file. Use `Pkg.resolve` to generate a Manifest.toml from a Project.toml.")
    end
end

function nb_and_dir_environments_equal(notebook::Notebook, dir::String)
    try
        ensure_has_nbpkg(notebook)
        assert_has_manifest(dir)
        true
    catch
        false
    end && let
        read(notebook |> project_file) == read(dir |> project_file) &&
        read(notebook |> manifest_file) == read(dir |> manifest_file)
    end
end

function write_nb_to_dir(notebook::Notebook, dir::String)
    ensure_has_nbpkg(notebook)
    mkpath(dir)

    readwrite(notebook |> project_file, dir |> project_file)
    readwrite(notebook |> manifest_file, dir |> manifest_file)
end

function write_dir_to_nb(dir::String, notebook::Notebook)
    assert_has_manifest(dir)

    notebook.nbpkg_ctx = Pluto.PkgCompat.create_empty_ctx()

    readwrite(dir |> project_file, notebook |> project_file)
    readwrite(dir |> manifest_file, notebook |> manifest_file)

    save_notebook(notebook)
end

write_dir_to_nb(dir::String, notebook_path::String) = write_dir_to_nb(dir::String, load_notebook(notebook_path))

write_nb_to_dir(notebook_path::String, dir::String) = write_nb_to_dir(load_notebook(notebook_path), dir)
nb_and_dir_environments_equal(notebook_path::String, dir::String) = nb_and_dir_environments_equal(load_notebook(notebook_path), dir)

"""
```julia
reset_notebook_environment(notebook_path::String; keep_project::Bool=false, backup::Bool=true)
```

Remove the embedded `Project.toml` and `Manifest.toml` from a notebook file, modifying the file. If `keep_project` is true, only `Manifest.toml` will be deleted. A backup file is created by default.
"""
function reset_notebook_environment(path::String; kwargs...)
    Pluto.reset_nbpkg(
        load_notebook_nobackup(path);
        kwargs...
    )
end

"""
```julia
reset_notebook_environment(notebook_path::String; backup::Bool=true, level::Pkg.UpgradeLevel=Pkg.UPLEVEL_MAJOR)
```

Update the embedded `Project.toml` and `Manifest.toml` in a notebook file, modifying the file. A [`Pkg.UpgradeLevel`](@ref) can be passed to the `level` keyword argument. A backup file is created by default.
"""
function update_notebook_environment(path::String; kwargs...)
    Pluto.update_nbpkg(
        Pluto.ServerSession(),
        load_notebook_nobackup(path);
        kwargs...
    )
end

function activate_notebook_environment(path::String)
    notebook_ref = Ref(load_notebook(path))

    ensure_has_nbpkg(notebook_ref[])

    ourpath = joinpath(mktempdir(), basename(path))
    mkpath(ourpath)

    still_needed = Ref(true)

    save_token = Token()

    function maybe_update_nb()
        withtoken(save_token) do
            if still_needed[]
                if !nb_and_dir_environments_equal(notebook_ref[], ourpath)
                    write_dir_to_nb(ourpath, notebook_ref[])
                    println()
                    @info "Saved notebook package environment ✓"
                    println()
                end
            end
        end
    end

    function maybe_update_dir()
        withtoken(save_token) do
            if still_needed[]
                if !nb_and_dir_environments_equal(notebook_ref[], ourpath)
                    write_nb_to_dir(notebook_ref[], ourpath)
                    println()
                    @info "New notebook package environment written to directory ✓"
                    println()
                end
            end
        end
    end

    maybe_update_dir()

    atexit(maybe_update_nb)

    Base.ACTIVE_PROJECT[] = ourpath

    # WATCH DIR PROJECT FILE
    Pluto.@asynclog begin
        while Base.ACTIVE_PROJECT[] == ourpath
            FileWatching.watch_file(ourpath |> project_file)
            # @warn "DIR PROJECT UPDATED"
            sleep(.2)

            maybe_update_nb()
        end
        still_needed[] = false
    end

    # WATCH DIR MANIFEST FILE
    Pluto.@asynclog begin
        while Base.ACTIVE_PROJECT[] == ourpath
            FileWatching.watch_file(ourpath |> manifest_file)
            # @warn "DIR MANIFEST UPDATED"
            sleep(.5)

            maybe_update_nb()
        end
        still_needed[] = false
    end

    # WATCH NOTEBOOK FILE
    Pluto.@asynclog begin
        while Base.ACTIVE_PROJECT[] == ourpath
            FileWatching.watch_file(path)
            # @warn "NOTEBOOK FILE UPDATED"
            # we update the dir after a longer delay, because changes to the environment in the directory (written to the notebook) take precedence.
            sleep(1)

            notebook_ref[] = load_notebook(path)

            maybe_update_dir()
        end
        still_needed[] = false
    end

    # CHECK EVERY 5 SECONDS
    # Pluto.@asynclog begin
    #     while still_needed[]
    #         sleep(5)
    #         maybe_update_nb()
    #         # we update the dir second, because changes to the environment in the directory (written to the notebook) take precedence.
    #         maybe_update_dir()
    #     end
    # end

    println()
    """

    > Notebook environment activated!

    ## Step 1.
    _Press `]` to open the Pkg REPL._

    The notebook environment is currently active.

    ## Step 2.
    The notebook file and your REPL environment are now synced. This means that:
    1. Any changes you make in the REPL environment will be written to the notebook file. For example, you can `pkg> update` or `pkg> add SomePackage`, and the notebook file will update.
    2. Whenever the notebook file changes, the REPL environment will be updated from the notebook file.

    ## Step 3.
    When you are done, you can exit the notebook environment by deactivating it:

    ```
    pkg> activate
    ```
    """ |> Markdown.parse |> display
    println()

end

const activate_notebook = activate_notebook_environment

function testnb()
    t = tempname()

    readwrite(Pluto.project_relative_path("test","packages","nb.jl"), t)
    t
end

end


text:
module PandorasBoxLever

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/lever" Lever(x::Integer, y::Integer, active::Bool=false, flag::String="")

const placements = Ahorn.PlacementDict(
    "Lever (Pandora's Box)" => Ahorn.EntityPlacement(
        Lever
    )
)

sprite = "objects/pandorasBox/lever/lever0"

function getSprite(lever::Lever)
    return get(lever, "active", false) ? "objects/pandorasBox/lever/lever4" : "objects/pandorasBox/lever/lever0"
end

function Ahorn.selection(lever::Lever)
    x, y = Ahorn.position(lever)
    sprite = getSprite(lever)

    return Ahorn.getSpriteRectangle(sprite, x, y, jx=0.5, jy=1.0)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, lever::Lever, room::Maple.Room) = Ahorn.drawSprite(ctx, getSprite(lever), 0, 0, jx=0.5, jy=1.0)

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

# This is loaded via @require

using .AbstractPlotting: RGB, N0f8, Vec3f0, Scene, SceneSpace, Axis, meshscatter!, scatter!, lines!

export spacefilling, ballstick

colortype(c::Color) = RGB{N0f8}(c.r/255, c.g/255, c.b/255)

"""
    spacefilling(mol::UndirectGraph; radius="van der Waals", tform=identity, H::Bool=true)

Represent `mol` as a space-filling (Calotte) model in three dimensions. `mol` should have 3d atom positions represented
in Angstroms. (3D SDF files can be downloaded from sites such as PubChem.) The two supported options for `radius` are
`"van der Waals"` and `"covalent"`; the former are available only for main-group elements, and the latter are available for
all.

`tform(xyz)` optionally transforms the positions before plotting them. `H=false` causes hydrogens to be omitted from the plot.

This function requires that you load one of the backends of the Makie/WGLMakie/CairoMakie family.
"""
spacefilling(mol::UndirectedGraph; kwargs...) = spacefilling!(Scene(), mol; kwargs...)

function spacefilling!(scene::Scene, mol::UndirectedGraph; radius::String="van der Waals", tform=identity, H::Bool=true)
    syms = [a.symbol for a in nodeattrs(mol)]
    pos = reduce(hcat, [tform(atom.coords) for atom in nodeattrs(mol)])
    size(pos, 1) == 3 || error("this plots only in 3d")
    if !H
        keep = syms .!= :H
        syms, pos = syms[keep], pos[:,keep]
        atomidx = findall(keep)
    else
        atomidx = collect(1:length(syms))
    end
    isvdW = radius == "van der Waals"
    radii = map(atomidx, syms) do idx, sym
        an = ATOMSYMBOLMAP[string(sym)]
        if isvdW
            r = ATOM_VANDERWAALS_RADII[an]
        else
            r = ATOM_COVALENT_RADII[an]
            if !isa(r, Real)
                # Carbon and a few metals have multiple values
                if an == 6
                    nbrs = neighbors(mol, idx)
                    key = nbrs == 3 ? "Csp3" :
                        nbrs == 2 ? "Csp2" : "Csp"
                    r = r[key]
                else
                    # For metals, it's safest to choose the smallest radius
                    r = minimum(values(r))
                end
            end
        end
        r
    end
    return axisoff!(meshscatter!(scene, pos[1,:], pos[2,:], pos[3,:], color=colortype.(atomcolor(syms; setting=DRAW_SETTING3)), markersize=radii, markerspace=SceneSpace))
end

"""
    ballstick(mol::UndirectGraph; tform=identity, H::Bool=true, markersize=1)

Represent `mol` as a ball-and-stick model in three dimensions. `mol` should have 3d atom positions represented
in Angstroms. 3D SDF files can be downloaded from sites such as PubChem.

`tform(xyz)` optionally transforms the positions before plotting them. `H=false` causes hydrogens to be omitted from the plot.
`markersize` specifies the diameter of the balls, in Angstroms.

This function requires that you load one of the backends of the Makie/WGLMakie/CairoMakie family.
"""
ballstick(mol::UndirectedGraph; kwargs...) = ballstick!(Scene(), mol; kwargs...)

function ballstick!(scene::Scene, mol::UndirectedGraph; tform=identity, H::Bool=true, markersize=1)
    syms = [a.symbol for a in nodeattrs(mol)]
    pos = reduce(hcat, [tform(atom.coords) for atom in nodeattrs(mol)])
    size(pos, 1) == 3 || error("this plots only in 3d")
    bondinfo = map(edgesiter(mol), edgeattrs(mol)) do e, attr
        ([Vec3f0(pos[:,e[1]]), Vec3f0(pos[:,e[2]])], attr.order)
    end
    if !H
        keep = syms .!= :H
        syms, pos = syms[keep], pos[:,keep]
        bondinfo = [bi for (bi, e) in zip(bondinfo, edgesiter(mol)) if (keep[e[1]] & keep[e[2]])]
    end
    scene = scatter!(scene, pos[1,:], pos[2,:], pos[3,:], color=colortype.(atomcolor(syms; setting=DRAW_SETTING3)), markersize=markersize, markerspace=SceneSpace)
    for (l, w) in bondinfo
        lines!(scene, l, linewidth=10*w)
    end
    return axisoff!(scene)
end

axisoff!(scene::Scene) = (axisoff!(scene[Axis]); return scene)
function axisoff!(ax)
    # ax.showgrid[] = (false, false, false)   # disabled due to https://github.com/JuliaPlots/Makie.jl/issues/774
    ax.showaxis[] = (false, false, false)
    ax.showticks[] = (false, false, false)
    return ax
end


text:
"""this is used to download the data from biomd and the sbml-test-suite"""
module SBMLBioModelsRepository

const datadir = joinpath(@__DIR__, "..", "data")
const logdir = joinpath(@__DIR__, "..", "test", "logs")

using CSV, DataFrames, JSON3, JSONTables, Glob
using Base.Threads, Base.Iterators, Downloads

using ModelingToolkit, OrdinaryDiffEq, Sundials
using CSV, DataFrames, BenchmarkTools
using SBML, SBMLToolkit
using Plots

include("lower.jl")
include("biomd.jl")
include("suite.jl")

export jsonfn_to_df, curl_biomd_xmls, curl_biomd_metadata, read_biomd_metadf, get_biomd_fns, biomodels
export sbml_test_suite, get_sbml_suite_fns, setup_settings_txt, verify_case, verify_all, verify_plot
export datadir, logdir
export lower_one, lower_fns

end


text:
using Test

function alloc_con(con,z)
    p,w = length(con), TO.width(con)
    ∇c = zeros(p,w)
    allocs  = @allocated evaluate(con,z)
    allocs += @allocated jacobian!(∇c, con,z)
end

model = Dynamics.DubinsCar()
x,u = rand(model)
z = KnotPoint(x,u,0.1)
n,m = size(model)

# All static
A = @SMatrix rand(4,3)
b = @SVector rand(4)
∇c = zeros(4,3)
con = LinearConstraint{Inequality,State}(n,m,A,b)
evaluate(con, z)
@test length(con) == 4
@test TO.width(con) == 3
@test evaluate(con, z) == A*x - b
jacobian!(∇c, con, x)
@test ∇c == A
@test alloc_con(con,z) == 0
@test state_dim(con) == n
@test_throws MethodError control_dim(con)

# both dynamic
con = LinearConstraint{Inequality,State}(n,m,Matrix(A),Vector(b))
@test evaluate(con, z) == A*x - b
jacobian!(∇c, con, x)
@test ∇c == A
@test alloc_con(con,z) == 0

# mixed
con = LinearConstraint{Inequality,State}(n,m,A,Vector(b))
@test evaluate(con, z) == A*x - b
jacobian!(∇c, con, x)
@test ∇c == A
@test alloc_con(con,z) == 0

# wrong input size
@test_throws AssertionError LinearConstraint{Inequality,State}(m,m,A,b)
@test_throws AssertionError LinearConstraint{Inequality,Control}(n,m,A,b)
con = LinearConstraint{Inequality,Control}(n,n,A,b)
@test evaluate(con, x) == A*x - b
jacobian!(∇c, con, x)
@test ∇c == A

model = Dynamics.Quadrotor2{UnitQuaternion{Float64,VectorPart}}()
x,u = rand(model)
n,m = size(model)

# Goal Constraint
xf = @SVector rand(n)
con = GoalConstraint(xf)
size(con) == (13,13)
∇c = zeros(size(con)...)
@test evaluate(con, x) == x-xf
@test jacobian!(∇c, con, x) == I(n)

con = GoalConstraint(xf, 1:3)
∇c = zeros(size(con))
@test evaluate(con, x) == (x-xf)[1:3]
@test jacobian!(∇c, con, x) == Matrix(I,3,n)

noquat = collect(1:n)
noquat = deleteat!(noquat, 4:7)
noquat = SVector{9}(noquat)
con = GoalConstraint(xf, noquat)
∇c = zeros(size(con))
@test evaluate(con, x) == (x-xf)[noquat]
@test jacobian!(∇c, con, x) == I(n)[noquat,:]

# Sphere Constraint
xc = @SVector rand(4)
yc = @SVector rand(4)
zc = @SVector rand(4)
r = @SVector fill(0.1,4)
con = SphereConstraint(n, xc, yc, zc, r)
∇c = zeros(size(con))
@test state_dim(con) == n
@test evaluate(con, x) isa SVector{4}
@test norm(jacobian!(∇c, con, x)[:,4:end],Inf) == 0

# Norm Constraint
con = TO.NormConstraint{Equality,Control}(m,4.0)
@test TO.sense(con) == Equality
@test TO.contype(con) == Control
@test evaluate(con, u)[1] == norm(u)^2 - 4
@test_throws MethodError state_dim(con)
@test control_dim(con) == m
@test TO.check_dims(con,n,m)

con = TO.NormConstraint{Inequality,State}(n, 2.0)
@test TO.sense(con) == Inequality
@test TO.contype(con) == State
@test evaluate(con, x)[1] == norm(x)^2 - 2
@test_throws MethodError control_dim(con)
@test state_dim(con) == n
@test TO.check_dims(con,n,m)

# Bounds Constraint
x_max = @SVector fill(+1,n)
x_min = @SVector fill(-1,n)
u_max = @SVector fill(+2,m)
u_min = @SVector fill(-2,m)

bnd = BoundConstraint(n,m, x_max=x_max, x_min=x_min)
@test_throws ArgumentError BoundConstraint(n,m, x_max=x_min, x_min=x_max)
@test_nowarn BoundConstraint(n,m, x_max=1, x_min=x_min)
@test_nowarn BoundConstraint(n,m, x_max=x_max, x_min=-1.)
@test TO.is_bound(bnd)
@test TO.upper_bound(bnd) == [x_max; u_max*Inf]
@test TO.lower_bound(bnd) == [x_min; u_min*Inf]

# Variable Bound Constraint
N = 101
X_max = [copy(x_max) for k = 1:N]
bnd = VariableBoundConstraint(n,m,N, x_max=X_max)
@test state_dim(bnd) == n
@test control_dim(bnd) == m
vals = [@SVector zeros(n) for k = 1:N]
Z = Traj(x,u,0.1,N)
@test_nowarn evaluate!(vals, bnd, Z)
@test jacobian(bnd, Z[1]) == Matrix(I,n,n+m)

# Indexed Constraint
bnd = BoundConstraint(n,m, x_max=x_max, x_min=x_min)
con = TO.IndexedConstraint(2n,2m, bnd)
∇c = zeros(size(con))
x2 = repeat(x,2)
u2 = repeat(u,2)
z = KnotPoint(x,u,0.1)
z2 = KnotPoint(x2,u2,0.1)
@test evaluate(bnd, z) ≈ evaluate(con, z2)
∇c0 = zeros(size(bnd))
jacobian!(∇c0, bnd, z)
jacobian!(∇c, con, z2)
@test ∇c ≈ [∇c0 zeros(26, 17)]

@test TO.width(con) == 2(n+m)
@test TO.width(bnd) == n+m

con = TO.NormConstraint{Equality,Control}(m,4.0)
∇c0 = zeros(size(con))
ix = n .+ @SVector [i for i in 1:n]
iu = m .+ @SVector [i for i in 1:m]
idx = TO.IndexedConstraint(2n,2m, con, ix, iu)
∇c = zeros(size(idx))
@test TO.contype(idx) == Control
@test TO.sense(idx) == Equality
@test evaluate(idx, z2) == evaluate(con, z)
jacobian!(∇c0, con, z)
jacobian!(∇c, idx, z2)
@test ∇c == [zeros(1,m) ∇c0]

con = TO.NormConstraint{Equality,State}(n,4.0)
idx = TO.IndexedConstraint(2n,2m, con)
@test evaluate(idx, z2) == evaluate(con, z)
∇c0 = zeros(size(con))
∇c = zeros(size(idx))
jacobian!(∇c0, con, z)
jacobian!(∇c, idx, z2)
@test ∇c == [∇c0 zeros(1,n)]


text:
function approve!(repo::GitHub.Repo,
                  pr::GitHub.PullRequest,
                  pr_head_commit_sha_to_approve::String;
                  auth::GitHub.Authorization,
                  body::String="",
                  whoami)
    if pr.user.login == whoami
    else
        repo_full_name = full_name(repo)
        pr_number = number(pr)
        endpoint = "/repos/$(repo_full_name)/pulls/$(pr_number)/reviews"
        approving_review_body = """
                $body
                <!---
                approved_pr_head_commit_sha=\"$(pr_head_commit_sha_to_approve)\"
                --->
                """
        myparams = Dict("event" => "APPROVE",
                        "body" => approving_review_body)
        GitHub.gh_post_json(GitHub.DEFAULT_API,
                            endpoint;
                            auth=auth,
                            params = myparams)
    end
    return nothing
end

author_login(pull_request::GitHub.PullRequest) = pull_request.user.login

base_repo(pull_request::GitHub.PullRequest) = pull_request.base.repo

body(c::GitHub.Comment) = c.body

function created_at(pull_request::GitHub.PullRequest)
    result = time_is_already_in_utc(pull_request.created_at)
    return result
end

function delete_all_of_my_reviews!(repo::GitHub.Repo,
                                  pr::GitHub.PullRequest;
                                  auth::GitHub.Authorization,
                                  whoami::String)
    all_pr_reviews = get_all_pull_request_reviews(repo, pr; auth = auth)
    for rev in all_pr_reviews
        if reviewer_login(rev) == whoami
            delete_pr_review!(repo, pr, rev; auth = auth)
        end
    end
    return nothing
end

function delete_merged_branch!(repo::GitHub.Repo, pr::GitHub.PullRequest; auth::GitHub.Authorization)
    updated_pr = _get_updated_pull_request(pr; auth=auth)
    if is_merged(updated_pr)
        try
            head_branch = pull_request_head_branch(updated_pr)
            repo = head_branch.repo
            ref = "heads/$(head_branch.ref)"
            GitHub.delete_reference(repo, ref; auth=auth)
        catch ex
            showerror(stderr, ex)
            Base.show_backtrace(stderr, catch_backtrace())
            println(stderr)
        end
    end
    return nothing
end

function delete_pr_review!(repo::GitHub.Repo, pr::GitHub.PullRequest, r::GitHub.Review; auth::GitHub.Authorization)
    return nothing
    repo_full_name = full_name(repo)
    pr_number = number(pr)
    review_id = r.id
    endpoint = "/repos/$(repo_full_name)/pulls/$(pr_number)/reviews/$(review_id)"
    GitHub.gh_delete_json(GitHub.DEFAULT_API,
                          endpoint;
                          auth=auth)
    return nothing
end

full_name(repo::GitHub.Repo) = repo.full_name

function _get_updated_pull_request(pull_request::GitHub.PullRequest; auth::GitHub.Authorization)
    pr_base_repo = base_repo(pull_request)
    pr_number = number(pull_request)
    updated_pr = GitHub.pull_request(pr_base_repo, pr_number; auth=auth)
    return updated_pr
end

function get_all_pull_request_comments(repo::GitHub.Repo,
                                       pr::GitHub.PullRequest;
                                       auth::GitHub.Authorization)
    all_comments = Vector{GitHub.Comment}(undef, 0)
    myparams = Dict("per_page" => 100, "page" => 1)
    cs, page_data = GitHub.comments(repo, pr, :pr; auth=auth, params = myparams, page_limit = 100)
    append!(all_comments, cs)
    while haskey(page_data, "next")
        cs, page_data =  GitHub.comments(repo, pr, :pr; auth=auth, page_limit = 100, start_page = page_data["next"])
        append!(all_comments, cs)
    end
    unique!(all_comments)
    return all_comments
end

function get_all_pull_request_reviews(repo::GitHub.Repo,
                                      pr::GitHub.PullRequest;
                                      auth::GitHub.Authorization)
    all_reviews = Vector{GitHub.Review}(undef, 0)
    myparams = Dict("per_page" => 100, "page" => 1)
    revs, page_data = GitHub.reviews(repo, pr; auth=auth, params = myparams, page_limit = 100)
    append!(all_reviews, revs)
    while haskey(page_data, "next")
        revs, page_data = GitHub.reviews(repo, pr; auth=auth, page_limit = 100, start_page = page_data["next"])
        append!(all_reviews, revs)
    end
    unique!(all_reviews)
    return all_reviews
end

function get_all_pull_requests(repo::GitHub.Repo,
                               state::String;
                               auth::GitHub.Authorization)
    all_pull_requests = Vector{GitHub.PullRequest}(undef, 0)
    myparams = Dict("state" => state, "per_page" => 100, "page" => 1)
    prs, page_data = GitHub.pull_requests(repo; auth=auth, params = myparams, page_limit = 100)
    append!(all_pull_requests, prs)
    while haskey(page_data, "next")
        prs, page_data = GitHub.pull_requests(repo; auth=auth, page_limit = 100, start_page = page_data["next"])
        append!(all_pull_requests, prs)
    end
    unique!(all_pull_requests)
    return all_pull_requests
end

function get_changed_filenames(repo::GitHub.Repo, pull_request::GitHub.PullRequest; auth::GitHub.Authorization)
    files = GitHub.pull_request_files(repo, pull_request; auth=auth)
    n = length(files)
    filenames = Vector{String}(undef, n)
    for i = 1:n
        filenames[i] = files[i].filename
    end
    return filenames
end

is_merged(pull_request::GitHub.PullRequest) = pull_request.merged

function is_open(pull_request::GitHub.PullRequest)
    result = pr_state(pull_request) == "open"
    !result && @error("Pull request is not open")
    return result
end

function merge!(registry_repo::GitHub.Repo,
                pr::GitHub.PullRequest,
                approved_pr_head_sha::AbstractString;
                auth::GitHub.Authorization)
    pr = wait_pr_compute_mergeability(registry_repo, pr; auth = auth)
    _approved_pr_head_sha = convert(String, strip(approved_pr_head_sha))::String
    pr_number = number(pr)
    @info("Attempting to squash-merge pull request #$(pr_number)")
    @debug("sha = $(_approved_pr_head_sha)")
    @debug("pr.mergeable = $(pr.mergeable)")
    params = Dict("sha" => _approved_pr_head_sha, "merge_method" => "squash")
    try
        GitHub.merge_pull_request(registry_repo,
                                  pr_number;
                                  auth=auth,
                                  params=params)
    catch ex
        showerror(stderr, ex)
        Base.show_backtrace(stderr, catch_backtrace())
        println(stderr)
    end
    try
        delete_merged_branch!(registry_repo, pr; auth=auth)
    catch
    end
    return nothing
end

function wait_pr_compute_mergeability(repo::GitHub.Repo,
                                      pr::GitHub.PullRequest;
                                      auth::GitHub.Authorization)
    sleep(5)
    max_tries = 10
    num_tries = 0
    pr = GitHub.pull_request(repo, pr.number; auth = auth)
    while !(pr.mergeable isa Bool) && num_tries <= max_tries
        num_tries = num_tries + 1
        sleep(5)
        pr = GitHub.pull_request(repo, pr.number; auth = auth)
    end
    return pr
end

num_changed_files(pull_request::GitHub.PullRequest) = pull_request.changed_files

number(pull_request::GitHub.PullRequest) = pull_request.number

function post_comment!(repo::GitHub.Repo,
                       pr::GitHub.PullRequest,
                       body::String;
                       auth::GitHub.Authorization)
    myparams = Dict("body" => body)
    GitHub.create_comment(repo, pr, :pr; auth=auth, params = myparams)
    return nothing
end

pull_request_head_branch(pull_request::GitHub.PullRequest) = pull_request.head

pull_request_head_sha(pull_request::GitHub.PullRequest) = pull_request.head.sha

repo_url(repo::GitHub.Repo) = repo.html_url.uri

pr_state(pull_request::GitHub.PullRequest) = pull_request.state

review_state(r::GitHub.Review) = r.state

reviewer_login(r::GitHub.Review) = r.user.login

function time_since_pr_creation(pull_request::GitHub.PullRequest)
    _pr_created_at_utc = created_at(pull_request)
    _now_utc = now_utc()
    result = _now_utc - _pr_created_at_utc
    return result
end

title(pull_request::GitHub.PullRequest) = pull_request.title

function username(auth::GitHub.Authorization)
    user_information = GitHub.gh_get_json(GitHub.DEFAULT_API,
                                          "/user";
                                          auth = auth)
    return user_information["login"]::String
end


text:
#=
Code related with input output (IO) of .nc files directly to/from ClimArrays
utilizing the NCDatasets.jl package and a buttload of convenience code.
An initial version of parts of this code was taken from:
https://github.com/rafaqz/GeoData.jl
=#
using NCDatasets: NCDatasets, NCDataset
export NCDatasets, NCDataset
export nckeys, ncdetails, globalattr, ncsize
export ncread, ncwrite

dim_to_commonname(::Lat) = "lat"
dim_to_commonname(::Lon) = "lon"
dim_to_commonname(::Time) = "time"
dim_to_commonname(::Pre) = "level"
dim_to_commonname(D::Dim) = string(DimensionalData.name(D))
dim_to_commonname(::Coord) = "cell"

const POSSIBLE_CELL_NAMES = ("ncells", "cell", "rgrid", "grid")

"""
    nckeys(file::String)
Return all keys of the `.nc` file in `file`.
"""
function nckeys(path::String)
    NCDataset(path) do ds
        return keys(ds)
    end
end
nckeys(a::NCDataset) = keys(a)

"""
    ncdetails(file, io = stdout)
Print details about the `.nc` file in `file` on `io`.
"""
function ncdetails(file, io = stdout)
    NCDataset(file) do ds
        show(io, MIME"text/plain"(), ds)
    end
end
ncdetails(ds::NCDataset, io = stdout) = show(io, MIME"text/plain"(), ds)

"""
    ncsize(file, var)
Return the size of the variable of the `.nc` file without actually loading any data.
"""
function ncsize(file, var)
    NCDataset(file) do ds
        return size(ds[var])
    end
end

"""
    globalattr(file::String) → Dict
Return the global attributes of the .nc file.
"""
function globalattr(file::String)
    NCDataset(file) do ds
        return Dict(ds.attrib)
    end
end

#########################################################################
# Imports
#########################################################################
include("netcdf_read.jl")
include("netcdf_write.jl")


text:
# This file contains tests for the `@exception` macro

module TestException

using Exceptions
using SyntaxTree
using Test

# Print the header
println("\e[1;32mRUNNING\e[0m: exception.jl")

macro test_throws(exception::Expr, expr::Expr)
    return esc(
        quote
            try
                eval($(expr) |> linefilter!)
            catch e
                if e isa LoadError
                    @test sprint(showerror, e) == """
                    LoadError: $(sprint(showerror, $(exception)()))
                    in expression starting at <macrocall>:0"""
                elseif !(e isa $(exception))
                    rethrow()
                end
            end
        end
    )
end

@testset "@exception" begin

    @test_nowarn eval(quote @exception m0 end)
    @test_nowarn eval(quote @exception m1 file::String end)
    @test_nowarn eval(quote
        @exception m2 file::String s::Symbol context = begin
            docstring = "Doc" * docstring
        end
    end)
    @test_nowarn eval(quote @exception m3 file s::Symbol end)

    @test_throws(
        Exceptions.Internal.OnlyOneEquation,
        quote @exception(m, var=v1) end,
    )

    @test_throws(
        Exceptions.Internal.FieldsOnly,
        quote @exception(m, 1+1) end,
    )

end

# No arguments (macro)
@testset "@m0" begin

    @test_nowarn eval(quote @m0 e0 "Docstring" "ErrorMessage" end)

    @test_throws(
        Exceptions.Internal.DocstringIsNotAString,
        quote @m0 e 1+1 "ErrorMessage" end
    )

    @test_throws(
        Exceptions.Internal.ErrorMessageIsNotAString,
        quote @m0 e "ErrorMessage" 2+2 end
    )

end

# No arguments (exception)
@testset "e0" begin

    @test fieldnames(e0) == ()
    @test "$(@doc(e0))" == "Docstring\n"
    @test sprint(showerror, e0()) == "\n\nMain.TestException.e0:\nErrorMessage\n"

end

# One argument
@testset "@m1" begin

    @test_nowarn eval(quote @m1 e1 "Docstring" "ErrorMessage" end)

end

# One arguments
@testset "e1" begin

    @test fieldnames(e1) == (:file,)

end

# Two arguments and a context
@testset "@m2" begin

    @test_nowarn eval(quote @m2 e2 "string" "ErrorMessage" end)

end

# Two arguments and a context
@testset "e2" begin

    @test fieldnames(e2) == (:file, :s)
    @test "$(@doc(e2))" == "Docstring\n"

end

# Arguments with type `Any`
@testset "@m3" begin

    @test_nowarn eval(quote @m3 e3 "Docstring" "ErrorMessage" end)
    @test e3.types[1] == Any
    @test e3.types[2] == Symbol

end

end


text:
using Test, Random
using LinearAlgebra
using Arpack
using Polynomials

rng = MersenneTwister(123)
for n in [3, 5, 10, 30, 100, 1000]
    P = randn(rng, n, n); P = (P + P')/2;
    q = randn(rng, n)/10
    m = Int(floor(n/2));
    A = randn(rng, m, n)
    b = randn(rng, m)
    # C_ = randn(n, n); C_ = (C_ + C_')/2;
    # Make C positive definite
    # C_ -= 1.01*min(minimum(eigvals(C_)), 0)*I
    r = [1e-4 1e-2 1 10 100 1000]
    eye = Matrix{Float64}(I, n, n)

    P_ = zeros(0, 0)
    q_ = zeros(0)
    r_ = 0
    N = zeros(0, 0)
    x0 = zeros(0)
    for i = 1:length(r)
        if rand(rng) < .1
            # Make P positive definite
            shift = 1.1*abs(eigs(P, nev=1, which=:SR)[1][1])
            P = P + shift*I;
        elseif rand(rng) < .1
            # Make P negative definite
            shift = -1.1*abs(eigs(P, nev=1, which=:LR)[1][1])
            P = P + shift*I;
        end

        N = nullspace(A)
        x = A\b
        d = N*(N'*randn(n))
        P_ = Symmetric(N'*P*N)

        if norm(x) <= r[i]
            alpha = roots(Polynomial([norm(x)^2 - r[i]^2, 2*d'*x, norm(d)^2]))
            x += alpha[1]*d

            x0 = x - N*(N'*x)
            q_ = N'*(q + P*x0)

            r_ = sqrt(r[i]^2 - norm(x0)^2)

            X, info = trs(P, q, r[i], A, b, compute_local=true)
            x_g = X[:, 1]
            X_reduced, info_reduced = TRS.trs(P_, q_, r_, compute_local=true)
            x_g_reduced = X_reduced[:, 1]
            str = "Constrained TRS.trs - r:"*string(r[i])
            @testset "$str" begin
                @test info.λ[1] - info_reduced.λ[1] <= 1e-6*abs(info_reduced.λ[1])
                @test norm(x_g - N*x_g_reduced - x0) <= 1e-3*r[i]
                if length(info_reduced.λ) > 1
                    diff_lengths = length(info.λ) - length(info_reduced.λ)
                    if diff_lengths < 0
                        info.λ = [info.λ; NaN*ones(diff_lengths)]
                    elseif diff_lengths > 0
                        info_reduced.λ= [info_reduced.λ; NaN*ones(diff_lengths)]
                    end
                    println("Lagrange Mutlipliers")
                    show(stdout, "text/plain", [info.λ info_reduced.λ]); println()
                    println("Norm difference of the second solution")
                    println(norm(X[:, 2] - N*X_reduced[:, 2] - x0))
                    @test info.λ[2] - info_reduced.λ[2] <= 1e-6*abs(info_reduced.λ[2])
                    @test norm(X[:, 2] - N*X_reduced[:, 2] - x0) <= 1e-3*r[i]
                end
            end
        end
    end
    # hard case
    λ_min, v, _ = eigs(-P, nev=1, which=:LR)
    v = v/norm(v)
    q = (I - v*v')*q
    q_ = N'*(q + P*x0)
    X, info = TRS.trs(P, q, r[end], A, b, compute_local=true)
    x_g = X[:, 1]
    X_reduced, info_reduced = TRS.trs(P_, q_, r_, compute_local=true)
    x_g_reduced = X_reduced[:, 1]
    @testset "Constrainted trs - hard case" begin
        @test info.λ[1] - info_reduced.λ[1] <= 1e-6*abs(info_reduced.λ[1])
        @test norm(x_g - N*x_g_reduced - x0) <= 1e-3*r[end]
        if length(info_reduced.λ) > 1
            diff_lengths = length(info.λ) - length(info_reduced.λ)
            if diff_lengths < 0
                info.λ = [info.λ; NaN*ones(diff_lengths)]
            elseif diff_lengths > 0
                info_reduced.λ= [info_reduced.λ; NaN*ones(diff_lengths)]
            end
            println("Lagrange Mutlipliers")
            show(stdout, "text/plain", [info.λ info_reduced.λ]); println()
            println("Norm difference of the second solution")
            println(norm(X[:, 2] - N*X_reduced[:, 2] - x0))
            @test info.λ[2] - info_reduced.λ[2] <= 1e-6*abs(info_reduced.λ[2])
            @test norm(X[:, 2] - N*X_reduced[:, 2] - x0) <= 1e-3*r[end]
        end
    end
end

nothing


text:
## show.jl.  show and print methods for NamedArray
## (c) 2013-2016 David A. van Leeuwen

## Julia type that implements a drop-in replacement of Array with named dimensions.

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

## Displaying things in the REPL drives me crazy, there are all these things
## like print, show, display, writemime, etc., and I have no clue as to
## what gets called how in which circumstance.

import Base.print, Base.show, Base.summary, Base.display
import SparseArrays: nonzeros, getcolptr, rowvals, nnz

## fallback
function summary(n::NamedArray{T,N,AT}) where {T,N,AT}
    return Base.dims2string(size(n)) * string(" Named ", AT)
end

print(n::NamedArray) = print(n.array)

## This seems to be the essential function to overload for displaying in REPL:
Base.show(io::IO, ::MIME"text/plain", n::NamedArray) = show(io, n)

## ndims==1 is dispatched below
function show(io::IO, n::NamedArray)
    print(io, summary(n))
    s = size(n)
    limit = get(io, :limit, false)
    if ndims(n) == 0
        println(io)
        show(io, n.array[1])
    elseif ndims(n) == 2
        println(io)
        if limit
            maxnrow = displaysize(io)[1] - 7 # summary, header, dots, + 2 empty lines, + 2 Julia prompt...
            show(io, n, min(maxnrow, s[1]))
        else
            show(io, n, s[1])
        end
    else
        nlinesneeded = prod(s[3:end]) * (s[1] + 3) + 1
        if limit && nlinesneeded > displaysize(io)[1]
            maxnrow = clamp((displaysize(io)[1] - 3) ÷ (prod(s[3:end])) - 3, 3, s[1])
        else
            maxnrow = s[1]
        end
        maxrepeat = displaysize(io)[1] ÷ (maxnrow + 4)
        i = 1
        for idx in CartesianIndices(size(n)[3:end])
            if i > maxrepeat
                print(io, "\n⋮")
                break
            end
            cartnames = [string(strdimnames(n, 2+i), "=", strnames(n, 2+i)[ind]) for (i, ind) in enumerate(idx.I)]
            println(io, "\n")
            println(io, "[:, :, ", join(cartnames, ", "), "] =")
            show(io, n[:, :, idx], maxnrow)
            i += 1
        end
    end
end

#show(io::IO, x::NamedVector) = invoke(show, (IO, NamedArray), io, x)

function show(io::IO, v::NamedVector)
    println(io, summary(v))
    limit = get(io, :limit, false)
    if limit
        maxnrow = displaysize(io)[1] - 7
        show(io, v, min(maxnrow, length(v)))
    else
        show(io, v, length(v))
    end
end

## compute the ranges to be displayed, plus a total index comprising all ranges.
function compute_range(maxn, n)
    if maxn < n
        hn = div(maxn,2)
        r = (1:hn, n-hn+1:n)
    else
        r = (1:n,)
    end
    totr = vcat(map(collect, r)...)
    r, totr
end

leftalign(s, l) = rpad(s, l, " ")
rightalign(s, l) = lpad(s, l, " ")
sprint_colpart(width::Int, s::Vector) = join(map(s->lpad(s, width, " "), s), "  ")
function sprint_row(namewidth::Int, name, width::Int, names::Tuple; dots="…", sep=" │ ")
    s = string(leftalign(name, namewidth), sep, sprint_colpart(width, names[1]))
    if length(names)>1
        s = string(s, "  ", dots, "  ", sprint_colpart(width, names[2]))
    end
    s
end

## for 2D printing
function show(io::IO, n::NamedMatrix, maxnrow::Int)
    @assert ndims(n)==2
    nrow, ncol = size(n)
    limit = get(io, :limit, false)
    ## rows
    rowrange, totrowrange = compute_range(maxnrow, nrow)
    s = [sprint(show, n.array[i,j], context=:compact => true) for i=totrowrange, j=1:ncol]
    rowname, colname = strnames(n)
    strlen(x) = length(string(x))
    colwidth = max(maximum(map(length, s); init=0), maximum(map(strlen, colname); init=0))
    rownamewidth = max(maximum(map(strlen, rowname); init=0), sum(map(length, strdimnames(n)))+3)
    if limit
        maxncol = div(displaysize(io)[2] - rownamewidth - 4, colwidth+2) # dots, spaces between
    else
        maxncol = ncol
    end
    ## columns
    colrange, totcorange = compute_range(maxncol, ncol)
    ## header
    header = sprint_row(rownamewidth, rightalign(join(strdimnames(n), " ╲ "), rownamewidth),
                        colwidth, map(i->colname[i], colrange))
    println(io, header)
    print(io, "─"^(rownamewidth+1), "┼", "─"^(length(header)-rownamewidth-2))
    ## data
    l = 1
    for i in 1:length(rowrange)
        if i > 1
            vdots = map(i->["⋮" for i=1:length(i)], colrange)
            println(io)
            print(io, sprint_row(rownamewidth, "⋮", colwidth, vdots, dots="⋱", sep="   "))
        end
        r = rowrange[i]
        for j in 1:length(r)
            row = s[l,:]
            if (i == 1 && j == 1) || (i == length(rowrange) && j == length(r))
                dots = "…"
            else
                dots = " "
            end
            println(io)
            print(io, sprint_row(rownamewidth, rowname[totrowrange[l]], colwidth,
                                 map(r -> row[r], colrange), dots=dots))
            l += 1
        end
    end
end

## special case of sparse matrix, based on base/sparse/sparsematrix.c
function show(io::IO, n::NamedArray{T1,2,SparseMatrixCSC{T1,T2}}) where {T1,T2}
    S = n.array
    print(io, S.m, "×", S.n, " Named sparse matrix with ", nnz(S), " ", eltype(S), " nonzero entries", nnz(S) == 0 ? "" : ":")

    rownames, colnames = strnames(n)
    rowpad = maximum([length(s) for s in rownames])
    colpad = maximum([length(s) for s in colnames])
    sep = "\n\t"

    function format_line(r, col, rowpad, colpad)
        print(io, sep, '[', rpad(rownames[rowvals(S)[r]], rowpad), ", ", lpad(colnames[col], colpad), "]  =  ")
        if isassigned(nonzeros(S), r)
            Base.show(io, nonzeros(S)[r])
        else
            print(io, Base.undef_ref_str)
        end
    end

    function get_cols(from, to)
        idx = eltype(getcolptr(S))[]
        c = searchsortedlast(getcolptr(S), from)
        for i = from:to
            while i == getcolptr(S)[c+1]
                c +=1
            end
            push!(idx, c)
        end
        idx
    end

    maxnrows = displaysize(io)[1] - 4
    if !get(io, :limit, false) || maxnrows >= nnz(S)
        cols = get_cols(1, nnz(S))
        format_line.(1:nnz(S), cols, rowpad, colpad)
    else
        s1, e1 = 1, div(maxnrows - 1, 2) # -1 accounts for \vdots
        s2, e2 = nnz(S) - (maxnrows - 1 - e1) + 1, nnz(S)
        cols1 = get_cols(s1, e1)
        format_line.(s1:e1, cols1, rowpad, colpad)
        print(io, sep, lpad("⋮", rowpad + colpad + 7))
        cols2 = get_cols(s2, e2)
        format_line.(s2:e2, cols2, rowpad, colpad)
    end
    nothing
end

function show(io::IO, v::NamedVector, maxnrow::Int)
    nrow = size(v, 1)
    rownames = strnames(v,1)
    rowrange, totrowrange = compute_range(maxnrow, nrow)
    s = [sprint(show, v.array[i], context=:compact => true) for i=totrowrange]
    colwidth = maximum(map(length,s); init = 0)
    rownamewidth = max(maximum(map(length, rownames); init = 0), 1+length(strdimnames(v)[1]))
    ## header
    println(io, string(leftalign(strdimnames(v, 1), rownamewidth), " │ "))
    print(io, "─"^(rownamewidth+1), "┼", "─"^(colwidth+1))
    ## data
    l = 1
    for i in 1:length(rowrange)
        if i > 1
            vdots = ["⋮"]
            println(io)
            print(io, sprint_row(rownamewidth, "⋮", colwidth, (vdots,), sep="   "))
        end
        r = rowrange[i]
        for j in 1:length(r)
            row = s[l]
            println(io)
            print(io, sprint_row(rownamewidth, rownames[totrowrange[l]], colwidth, ([row],)))
            l += 1
        end
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Advisory reentrant lock
"""
    ReentrantLock()

Creates a re-entrant lock for synchronizing [`Task`](@ref)s.
The same task can acquire the lock as many times as required.
Each [`lock`](@ref) must be matched with an [`unlock`](@ref).
"""
mutable struct ReentrantLock <: AbstractLock
    locked_by::Union{Task, Nothing}
    cond_wait::GenericCondition{Threads.SpinLock}
    reentrancy_cnt::Int

    ReentrantLock() = new(nothing, GenericCondition{Threads.SpinLock}(), 0)
end

"""
    islocked(lock) -> Status (Boolean)

Check whether the `lock` is held by any task/thread.
This should not be used for synchronization (see instead [`trylock`](@ref)).
"""
function islocked(rl::ReentrantLock)
    return rl.reentrancy_cnt != 0
end

"""
    trylock(lock) -> Success (Boolean)

Acquire the lock if it is available,
and return `true` if successful.
If the lock is already locked by a different task/thread,
return `false`.

Each successful `trylock` must be matched by an [`unlock`](@ref).
"""
function trylock(rl::ReentrantLock)
    t = current_task()
    lock(rl.cond_wait)
    try
        if rl.reentrancy_cnt == 0
            rl.locked_by = t
            rl.reentrancy_cnt = 1
            return true
        elseif t == notnothing(rl.locked_by)
            rl.reentrancy_cnt += 1
            return true
        end
        return false
    finally
        unlock(rl.cond_wait)
    end
end

"""
    lock(lock)

Acquire the `lock` when it becomes available.
If the lock is already locked by a different task/thread,
wait for it to become available.

Each `lock` must be matched by an [`unlock`](@ref).
"""
function lock(rl::ReentrantLock)
    t = current_task()
    lock(rl.cond_wait)
    try
        while true
            if rl.reentrancy_cnt == 0
                rl.locked_by = t
                rl.reentrancy_cnt = 1
                return
            elseif t == notnothing(rl.locked_by)
                rl.reentrancy_cnt += 1
                return
            end
            wait(rl.cond_wait)
        end
    finally
        unlock(rl.cond_wait)
    end
end

"""
    unlock(lock)

Releases ownership of the `lock`.

If this is a recursive lock which has been acquired before, decrement an
internal counter and return immediately.
"""
function unlock(rl::ReentrantLock)
    t = current_task()
    rl.reentrancy_cnt == 0 && error("unlock count must match lock count")
    rl.locked_by == t || error("unlock from wrong thread")
    lock(rl.cond_wait)
    try
        rl.reentrancy_cnt -= 1
        if rl.reentrancy_cnt == 0
            rl.locked_by = nothing
            notify(rl.cond_wait)
        end
    finally
        unlock(rl.cond_wait)
    end
    return
end

function unlockall(rl::ReentrantLock)
    t = current_task()
    n = rl.reentrancy_cnt
    rl.locked_by == t || error("unlock from wrong thread")
    n == 0 && error("unlock count must match lock count")
    lock(rl.cond_wait)
    try
        rl.reentrancy_cnt = 0
        rl.locked_by = nothing
        notify(rl.cond_wait)
    finally
        unlock(rl.cond_wait)
    end
    return n
end

function relockall(rl::ReentrantLock, n::Int)
    t = current_task()
    lock(rl)
    n1 = rl.reentrancy_cnt
    rl.reentrancy_cnt = n
    n1 == 1 || concurrency_violation()
    return
end

function lock(f, l::AbstractLock)
    lock(l)
    try
        return f()
    finally
        unlock(l)
    end
end

function trylock(f, l::AbstractLock)
    if trylock(l)
        try
            return f()
        finally
            unlock(l)
        end
    end
    return false
end

@eval Threads begin
    """
        Threads.Condition([lock])

    A thread-safe version of [`Base.Condition`](@ref).

    !!! compat "Julia 1.2"
        This functionality requires at least Julia 1.2.
    """
    const Condition = Base.GenericCondition{Base.ReentrantLock}

    """
    Special note for [`Threads.Condition`](@ref):

    The caller must be holding the [`lock`](@ref) that owns `c` before calling this method.
    The calling task will be blocked until some other task wakes it,
    usually by calling [`notify`](@ref)` on the same Condition object.
    The lock will be atomically released when blocking (even if it was locked recursively),
    and will be reacquired before returning.
    """
    wait(c::Condition)
end

const ThreadSynchronizer = GenericCondition{Threads.SpinLock}

"""
    Semaphore(sem_size)

Create a counting semaphore that allows at most `sem_size`
acquires to be in use at any time.
Each acquire must be matched with a release.
"""
mutable struct Semaphore
    sem_size::Int
    curr_cnt::Int
    cond_wait::Threads.Condition
    Semaphore(sem_size) = sem_size > 0 ? new(sem_size, 0, Threads.Condition()) : throw(ArgumentError("Semaphore size must be > 0"))
end

"""
    acquire(s::Semaphore)

Wait for one of the `sem_size` permits to be available,
blocking until one can be acquired.
"""
function acquire(s::Semaphore)
    lock(s.cond_wait)
    try
        while s.curr_cnt >= s.sem_size
            wait(s.cond_wait)
        end
        s.curr_cnt = s.curr_cnt + 1
    finally
        unlock(s.cond_wait)
    end
    return
end

"""
    release(s::Semaphore)

Return one permit to the pool,
possibly allowing another task to acquire it
and resume execution.
"""
function release(s::Semaphore)
    lock(s.cond_wait)
    try
        s.curr_cnt > 0 || error("release count must match acquire count")
        s.curr_cnt -= 1
        notify(s.cond_wait; all=false)
    finally
        unlock(s.cond_wait)
    end
    return
end

"""
    Event()

Create a level-triggered event source. Tasks that call [`wait`](@ref) on an
`Event` are suspended and queued until `notify` is called on the `Event`.
After `notify` is called, the `Event` remains in a signaled state and
tasks will no longer block when waiting for it.

!!! compat "Julia 1.1"
    This functionality requires at least Julia 1.1.
"""
mutable struct Event
    notify::Threads.Condition
    set::Bool
    Event() = new(Threads.Condition(), false)
end

function wait(e::Event)
    e.set && return
    lock(e.notify)
    try
        while !e.set
            wait(e.notify)
        end
    finally
        unlock(e.notify)
    end
    nothing
end

function notify(e::Event)
    lock(e.notify)
    try
        if !e.set
            e.set = true
            notify(e.notify)
        end
    finally
        unlock(e.notify)
    end
    nothing
end

@eval Threads begin
    import .Base: Event
    export Event
end

max_stars_repo_name: musm/julia

max_issues_repo_name: musm/julia

max_forks_repo_name: musm/julia


text:
using MetabolicFluxAnalysis
using OpenCL

# current tests..
fluxi_1_f = open("testdata/test_01.fluxi")
fluxi_1_s = readstring(fluxi_1_f)
fluxi1 = MetabolicFluxAnalysis.parseFluxi( fluxi_1_s )

# create an evaluator for the forward simulation
ilee_jul = MetabolicFluxAnalysis.create_julia_ile_evaluator( fluxi1.net , fluxi1.substrate_configs["SC_010"] ,fluxi1.measurements["MassDistribution_01"] )

# create opencl evaluator for the fw simulation:
ocl_devices = OpenCL.cl.devices()
ocl_device  = ocl_devices[1]
ocl_context =  MetabolicFluxAnalysis.OCLContext(ocl_device)
ilee_ocl = MetabolicFluxAnalysis.create_opencl_ile_evaluator( fluxi1.net , fluxi1.substrate_configs["SC_010"] ,fluxi1.measurements["MassDistribution_01"] , ocl_context )
print("\nOpenCL Device: $(ocl_device)\n")

# create flux polytope:

fp1 = MetabolicFluxAnalysis.get_flux_polyhedron( fluxi1.net , fluxi1.linear_constraints["ConstraintsExchange"] , max_val=200 )

# parametrize the flux polytope
rfd_parametrization = MetabolicFluxAnalysis.get_rounded_fulldimensional_parametrization( fp1 ; remove_red_constr=1)

@assert MetabolicFluxAnalysis.isinside( rfd_parametrization.p_opt , rfd_parametrization.cheby_p_opt[:,:] )

#                    v1  v2     v3   v4   v5   v6      Ain  Eout  Fout
reference_flux = [  100.; 110 ; 50 ; 20 ; 20 ; 80   ;  100 ; 60 ; 80 ]

ilee_jul.f_fwsim( ilee_jul.substrate_mds , reference_flux[:,[1]])
ilee_jul.f_fwsim( ilee_jul.substrate_mds , reference_flux[:,[1;1;1;1]])
ilee_ocl.f_fwsim( ilee_jul.substrate_mds , reference_flux[:,[1]])
ilee_ocl.f_fwsim( ilee_jul.substrate_mds , reference_flux[:,[1;1;1;1]])

flux_test_a = reference_flux[:,repmat([1],100000)]

tic(); results = ilee_jul.f_fwsim( ilee_jul.substrate_mds , flux_test_a ); toc()
tic(); results = ilee_ocl.f_fwsim( ilee_jul.substrate_mds , flux_test_a ); toc()

md_xx  = MetabolicFluxAnalysis.eval_flux2mds( ilee_jul , repmat(reference_flux[:,:],1,100) )
sse_xx = MetabolicFluxAnalysis.eval_flux2sse( ilee_jul , repmat(reference_flux[:,:],1,100) )

print(ilee_jul.fwsim.code)
print(ilee_ocl.fwsim.code)


text:
# Synchronization
using CEnum

export MemoryFlags, threadgroup_barrier, simdgroup_barrier

"""
    MemoryFlags

Flags to set the memory synchronization behavior of threadgroup\\_barrier and simdgroup\\_barrier.

Possible values:

    None: Set barriers to only act as an execution barrier and not apply a memory fence.

    Device: Ensure the GPU correctly orders the memory operations to device memory
            for threads in the threadgroup or simdgroup.

    ThreadGroup: Ensure the GPU correctly orders the memory operations to threadgroup
            memory for threads in a threadgroup or simdgroup.

    Texture: Ensure the GPU correctly orders the memory operations to texture memory for
            threads in a threadgroup or simdgroup for a texture with the read_write access qualifier.

    ThreadGroup_ImgBlock: Ensure the GPU correctly orders the memory operations to threadgroup imageblock memory
            for threads in a threadgroup or simdgroup.
"""
@cenum MemoryFlags::UInt32 begin
    MemoryFlagNone                  = 0
    MemoryFlagDevice                = 1
    MemoryFlagThreadGroup           = 2
    MemoryFlagTexture               = 4
    MemoryFlagThreadGroup_ImgBlock  = 8
end

"""
    threadgroup_barrier(flag::MemoryFlags=MemoryFlagNone)

Synchronize all threads in a threadgroup.

Possible flags that affect the memory synchronization behavior are found in [`MemoryFlags`](@ref)
"""
@inline threadgroup_barrier(flag::MemoryFlags=MemoryFlagNone) = ccall("extern air.wg.barrier", llvmcall, Cvoid, (Cuint, Cuint, ), flag, UInt32(1))

"""
    simdgroup_barrier(flag::MemoryFlags=MemoryFlagNone)

Synchronize all threads in a SIMD-group.

Possible flags that affect the memory synchronization behavior are found in [`MemoryFlags`](@ref)
"""
@inline simdgroup_barrier(flag::MemoryFlags=MemoryFlagNone) = ccall("extern air.simdgroup.barrier", llvmcall, Cvoid, (Cuint, Cuint, ), flag, UInt32(1))


text:
# ---
# title: 754. Reach a Number
# id: problem754
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/reach-a-number/description/>
# hidden: true
# ---
#
# You are standing at position `0` on an infinite number line. There is a goal
# at position `target`.
#
# On each move, you can either go left or right. During the _n_ -th move
# (starting from 1), you take _n_ steps.
#
# Return the minimum number of steps required to reach the destination.
#
# **Example 1:**
#
#
#
#     Input: target = 3
#     Output: 2
#     Explanation:
#     On the first move we step from 0 to 1.
#     On the second step we step from 1 to 3.
#
#
# **Example 2:**
#
#
#
#     Input: target = 2
#     Output: 3
#     Explanation:
#     On the first move we step from 0 to 1.
#     On the second move we step  from 1 to -1.
#     On the third move we step from -1 to 2.
#
#
# **Note:**
#
# * `target` will be a non-zero integer in the range `[-10^9, 10^9]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
##########################################
# Master file for running all test cases #
##########################################
using Zygote, ReverseDiff, Memoization, Turing; turnprogress(false)
using Pkg
using Random
using Test

include("test_utils/AllUtils.jl")

# Begin testing.
@testset "Turing" begin
    @testset "core" begin
        include("core/ad.jl")
        include("core/container.jl")
    end

    test_adbackends = if VERSION >= v"1.2"
        [:forwarddiff, :tracker, :reversediff]
    else
        [:forwarddiff, :tracker]
    end
    Turing.setrdcache(false)
    for adbackend in test_adbackends
        Turing.setadbackend(adbackend)
        @testset "inference: $adbackend" begin
            @testset "samplers" begin
                include("inference/gibbs.jl")
                include("inference/hmc.jl")
                include("inference/is.jl")
                include("inference/mh.jl")
                include("inference/ess.jl")
                include("inference/emcee.jl")
                include("inference/AdvancedSMC.jl")
                include("inference/Inference.jl")
                include("contrib/inference/dynamichmc.jl")
            end
        end

        @testset "variational algorithms : $adbackend" begin
            include("variational/advi.jl")
        end

        @testset "modes" begin
            include("modes/ModeEstimation.jl")
        end
    end
    @testset "variational optimisers" begin
        include("variational/optimisers.jl")
    end

    Turing.setadbackend(:forwarddiff)
    @testset "stdlib" begin
        include("stdlib/distributions.jl")
        include("stdlib/RandomMeasures.jl")
    end

    @testset "utilities" begin
        # include("utilities/stan-interface.jl")
        include("inference/utilities.jl")
    end
end


text:
# * Field creation
# ** Parameter calculation

@doc raw"""
    calc_params!(field_params)

This function performs the calculation of different quantities from
the information provided.
"""
function calc_params!(field_params::Dict{Symbol,Any})
    test_field_parameters(field_params, [:λ, :T, :f, :ν, :ω, :ħω])
    test_field_parameters(field_params, [:I₀, :E₀, :Uₚ, :A₀])

    for k in keys(field_params)
        field_params[k] = get_unitful_quantity(field_params, k)
    end

    @namespace!(field_params) do
        if λ || T
            if λ
                T = λ/u"c"
            elseif T
                λ = T*u"c"
            end
            ν = 1/λ
            f = 1/T
            ω = 2π*u"rad"*f
            ħω = u"ħ"*ω
        else # ∝ Frequency specified
            if f || ν
                if f
                    ν = f/u"c"
                elseif ν
                    f = ν*u"c"
                end
                ω = 2π*u"rad"*f
            elseif ω || ħω
                if ω
                    ħω = u"ħ"*ω
                elseif ħω
                    ω = ħω/u"ħ"
                end
                f = ω/(2π*u"rad")
                ν = f/u"c"
            end
            T = 1/f
            λ = 1/ν
        end

        if I₀ || Uₚ
            if Uₚ
                I₀ = Uₚ / (2*u"q"^2/(u"c"*u"ε0"*u"me")) * 4ω^2
            end
            E₀ = √(2I₀/(u"ε0"*u"c"))
            A₀ = E₀/ω
        elseif E₀ || A₀
            if E₀
                A₀ = E₀/ω
            else
                E₀ = ω*A₀
            end
            I₀ = u"ε0"*u"c"/2*E₀^2
        end
        if !Uₚ
            Uₚ = 2*u"q"^2/(u"c"*u"ε0"*u"me") * I₀/4ω^2
        end
    end

    if get(field_params, :kind, :linear) == :transverse
        field_params[:R] = compute_rotation(get(field_params, :rotation, nothing))
    end

    field_params
end

# ** Frontend macro

function make_field(field_params::Dict{Symbol,Any})
    is_transverse = (:ξ ∈ keys(field_params) ||
                     :rotation ∈ keys(field_params) ||
                     get(field_params, :carrier, :fixed) ∈ [:linear, :elliptical, :circular])
    kind = get(field_params, :kind, is_transverse ? :transverse : :linear)

    if kind == :linear
        calc_params!(field_params)

        # Maybe these two blocks can be implicitly deduced from the passed
        # parameters? E.g. if a chirp parameter is given, the carrier type
        # should autmatically be resolved as ChirpedCarrier. Similarly, if
        # ramp and flat are given, a trapezoidal pulse is requested.

        carrier_sym = get(field_params, :carrier,
                          :q ∉ keys(field_params) ? :fixed : :harmonic)
        carrier_sym ∉ keys(carrier_types) &&
            error("Unknown carrier type $(carrier_sym), valid choices are $(keys(carrier_types))")
        :q ∈ keys(field_params) && carrier_sym != :harmonic &&
            error("Invalid carrier type, $(carrier_sym), for field with harmonic components")
        carrier = carrier_types[carrier_sym](field_params)

        env_sym = get(field_params, :env, :gauss)
        env_sym ∉ keys(envelope_types) &&
            error("Unknown envelope type $(env_sym), valid choices are $(keys(envelope_types))")
        env = envelope_types[env_sym](field_params, carrier)

        :ξ in keys(field_params) || :rotation in keys(field_params) &&
            error("Elliptical (transverse) fields are not linear!")

        LinearField(carrier, env, field_params)
    elseif kind == :transverse
        field_params[:kind] = :transverse
        calc_params!(field_params)

        carrier_sym = get(field_params, :carrier, :ξ ∉ keys(field_params) ? :linear : :elliptical)
        carrier_sym ∉ keys(carrier_types) &&
            error("Unknown carrier type $(carrier_sym), valid choices are $(keys(carrier_types))")
        :ξ ∈ keys(field_params) && carrier_sym != :elliptical &&
            error("Invalid carrier type, $(carrier_sym), for field with elliptical polarization")
        carrier = carrier_types[carrier_sym](field_params)

        env_sym = get(field_params, :env, :gauss)
        env_sym ∉ keys(envelope_types) &&
            error("Unknown envelope type $(env_sym), valid choices are $(keys(envelope_types))")
        env = envelope_types[env_sym](field_params, carrier)

        TransverseField(carrier, env, field_params)
    elseif kind ∈ keys(field_types)
        field_types[kind](field_params)
    else
        throw(ArgumentError("Unknown field kind $(kind)"))
    end
end

"""
    @field(name) do ... end

Frontend macro for creating an electric field and storing it in the
variable `name`.

# Example

```julia
julia> @field(F) do
           I₀ = 2.0
           T = 2.0
           σ = 3.0
           Tmax = 3.0
       end
Linearly polarized field with
  - I₀ = 2.0000e+00 au = 7.0188904e16 W cm⁻² =>
    - E₀ = 1.4142e+00 au = 727.2178 GV m⁻¹
    - A₀ = 0.4502 au
  – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
  – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
```
"""
macro field(spec, var)
    spec.head == :-> ||
        error("Expected a block with parameters for definition of the field")
    block = spec.args[2]
    block.head == :block ||
        error("Expected a block with parameters for definition of the field")

    field_params = gensym()
    local tree = walk(spec, field_params, true)
    quote
        $field_params = Dict{Symbol,Any}()
        $(tree)()
        $(esc(var)) = make_field($field_params)
    end
end

export @field


text:
using Test
using LAMMPS

LMP() do lmp
    @test LAMMPS.version(lmp) >= 0
end

LMP(["-screen", "none"]) do lmp
    @test LAMMPS.version(lmp) >= 0
    command(lmp, "clear")

    @test_throws ErrorException command(lmp, "nonsense")
end


text:
immutable WrappedMap{T}<:AbstractLinearMap{T}
    lmap::Union{AbstractMatrix{T},AbstractLinearMap{T}}
    _isreal::Bool
    _issym::Bool
    _ishermitian::Bool
    _isposdef::Bool
    WrappedMap(A::Union{AbstractMatrix{T},AbstractLinearMap{T}};isreal::Bool=Base.isreal(A),issym::Bool=Base.issym(A),ishermitian::Bool=Base.ishermitian(A),isposdef::Bool=Base.isposdef(A))=new(A,isreal,issym,ishermitian,isposdef)
end
WrappedMap{T}(A::Union{AbstractMatrix{T},AbstractLinearMap{T}};isreal::Bool=Base.isreal(A),issym::Bool=Base.issym(A),ishermitian::Bool=Base.ishermitian(A),isposdef::Bool=Base.isposdef(A))=WrappedMap{T}(A;isreal=isreal,issym=issym,ishermitian=ishermitian,isposdef=isposdef)

# properties
Base.size(A::WrappedMap,n)=(n==1 || n==2 ? size(A.lmap,n) : error("AbstractLinearMap objects have only 2 dimensions"))
Base.size(A::WrappedMap)=size(A.lmap)
Base.isreal(A::WrappedMap)=A._isreal
Base.issym(A::WrappedMap)=A._issym
Base.ishermitian(A::WrappedMap)=A._ishermitian
Base.isposdef(A::WrappedMap)=A._isposdef

# comparison
==(A::WrappedMap,B::WrappedMap)=(A.lmap==B.lmap && isreal(A)==isreal(B) && issym(A)==issym(B) && ishermitian(A)==ishermitian(B) && isposdef(A)==isposdef(B))

# multiplication with vector
Base.A_mul_B!(y::AbstractVector,A::WrappedMap,x::AbstractVector)=Base.A_mul_B!(y,A.lmap,x)
*(A::WrappedMap,x::AbstractVector)=*(A.lmap,x)

Base.At_mul_B!(y::AbstractVector,A::WrappedMap,x::AbstractVector)=Base.At_mul_B!(y,A.lmap,x)
Base.At_mul_B(A::WrappedMap,x::AbstractVector)=Base.At_mul_B(A.lmap,x)

Base.Ac_mul_B!(y::AbstractVector,A::WrappedMap,x::AbstractVector)=Base.Ac_mul_B!(y,A.lmap,x)
Base.Ac_mul_B(A::WrappedMap,x::AbstractVector)=Base.Ac_mul_B(A.lmap,x)

# combine AbstractLinearMap and Matrix objects: linear combinations and map composition
+(A1::AbstractLinearMap,A2::AbstractMatrix)=+(A1,WrappedMap(A2))
+(A1::AbstractMatrix,A2::AbstractLinearMap)=+(WrappedMap(A1),A2)

*(A1::AbstractLinearMap,A2::AbstractMatrix)=*(A1,WrappedMap(A2))
*(A1::AbstractMatrix,A2::AbstractLinearMap)=*(WrappedMap(A1),A2)


text:
using ApproxFun, OrdinaryDiffEq, Sundials, BenchmarkTools, DiffEqOperators
using DiffEqDevTools
using LinearAlgebra
using Plots; gr()

S = Fourier()
n = 512
x = points(S, n)
D2 = Derivative(S,2)[1:n,1:n]
D  = (Derivative(S) → S)[1:n,1:n]
T = ApproxFun.plan_transform(S, n)
Ti = ApproxFun.plan_itransform(S, n)

D4 = Derivative(S,4)[1:n,1:n]
û₀ = T*(cos.(x./16).*(1 .+ sin.(x./2.04)))
tmp=similar(û₀)
q = (D,T,Ti,tmp,similar(tmp),similar(tmp))
function kuramoto_sivashinsky(dû,û,q,t)
    D,T,Ti,tmp,u,uc = q
    mul!(u, D, û)
    mul!(tmp, Ti, u)
    mul!(u, Ti, û)
    @. tmp=tmp*u
    mul!(u,T, tmp)
    #mul!(uc, D2, û)
    @. dû = - u
end

prob = SplitODEProblem(DiffEqArrayOperator(-Diagonal(D4+D2)), kuramoto_sivashinsky, û₀, (0.0,5.0), q)
sol  = solve(prob,RadauIIA5(autodiff=false); reltol=1e-14,abstol=1e-14)
test_sol = TestSolution(sol)

tslices=[0.0 1.0 2.0 3.0 5.0]
ys=hcat((Ti*sol(t) for t in tslices)...)
labels=["t=$t" for t in tslices]
plot(x,ys,label=labels)

diag_linsolve=LinSolveFactorize(W->let tmp = tmp
    for i in 1:size(W, 1)
        tmp[i] = W[i, i]
    end
    Diagonal(tmp)
end)

abstols = 0.1 .^ (5:8)
reltols = 0.1 .^ (1:4)
multipliers =  0.5 .^ (0:3)
setups = [Dict(:alg => IMEXEuler(linsolve=diag_linsolve), :dts => 1e-3 * multipliers),
          Dict(:alg => CNAB2(linsolve=diag_linsolve), :dts => 5e-3 * multipliers),
          Dict(:alg => CNLF2(linsolve=diag_linsolve), :dts => 5e-3 * multipliers),
          Dict(:alg => SBDF2(linsolve=diag_linsolve), :dts => 1e-3 * multipliers)]
labels = ["IMEXEuler" "CNAB2" "CNLF2" "SBDF2"]
@time wp1 = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true,names=labels,
                            numruns=5,seconds=5,
                            save_everystop=false,appxsol=test_sol,maxiters=Int(1e5));

plot(wp1,label=labels,markershape=:auto,title="IMEX methods, diagonal linsolve, low order")

abstols = 0.1 .^ (5:8) # all fixed dt methods so these don't matter much
reltols = 0.1 .^ (1:4)
multipliers = 0.5 .^ (0:3)
setups = [Dict(:alg => NorsettEuler(), :dts => 1e-3 * multipliers),
          Dict(:alg => ETDRK2(), :dts => 1e-2 * multipliers)]
labels = hcat("NorsettEuler",
              "ETDRK2 (caching)")
@time wp2 = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp2, label=labels, markershape=:auto, title="ExpRK methods, low order")

abstols = 0.1 .^ (5:8) # all fixed dt methods so these don't matter much
reltols = 0.1 .^ (1:4)
multipliers = 0.5 .^ (0:3)
setups = [Dict(:alg => CNAB2(linsolve=diag_linsolve), :dts => 5e-3 * multipliers)
          Dict(:alg => ETDRK2(), :dts => 1e-2 * multipliers)]
labels = ["CNAB2 (diagonal linsolve)" "ETDRK2"]
@time wp3 = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp3, label=labels, markershape=:auto, title="Between family, low orders")

abstols = 0.1 .^ (7:13)
reltols = 0.1 .^ (4:10)
setups = [Dict(:alg => ARKODE(Sundials.Implicit(), order=3, linear_solver=:Band, jac_upper=1, jac_lower=1)),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=4, linear_solver=:Band, jac_upper=1, jac_lower=1)),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:Band, jac_upper=1, jac_lower=1))]
labels = hcat("ARKODE3", "ARKODE4", "ARKODE5")
@time wp4 = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp4, label=labels, markershape=:auto, title="IMEX methods, band linsolve, medium order")

abstols = 0.1 .^ (7:11) # all fixed dt methods so these don't matter much
reltols = 0.1 .^ (4:8)
multipliers = 0.5 .^ (0:4)
setups = [Dict(:alg => ETDRK3(), :dts => 1e-2 * multipliers),
          Dict(:alg => ETDRK4(), :dts => 1e-2 * multipliers),
          Dict(:alg => HochOst4(), :dts => 1e-2 * multipliers)]
labels = hcat("ETDRK3 (caching)", "ETDRK4 (caching)",
              "HochOst4 (caching)")
@time wp5 = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp5, label=labels, markershape=:auto, title="ExpRK methods, medium order")

abstols = 0.1 .^ (7:11)
reltols = 0.1 .^ (4:8)
multipliers = 0.5 .^ (0:4)
setups = [Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:Band, jac_upper=1, jac_lower=1)),
          Dict(:alg => ETDRK3(), :dts => 1e-2 * multipliers),
          Dict(:alg => ETDRK4(), :dts => 1e-2 * multipliers)]
labels = hcat("ARKODE (nondiagonal linsolve)", "ETDRK3 ()", "ETDRK4 ()")
@time wp6 = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp6, label=labels, markershape=:auto, title="Between family, medium order")

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
# # Single degree-of-freedom

using ReachabilityAnalysis, StructuralDynamicsODESolvers

# ### Equations of motion

# Struct that holds a problem describing an harmonic oscillator with frequency ω:
# x''(t) + ω^2 x(t) = 0
#
# solution x(t) = Acos(ωt + B), v(t) = -ωAsin(ωt + B)
# x(0) = Acos(B)
# v(0) = -ωAsin(B)
#
# special case: if x(0) = A, v(0) = 0 => x(t) = A cos ωt
struct SDOF{T, ST}
    ivp::T
    X0::ST
    ω::Float64
    Ampl::Float64
    T::Float64
end

amplitude(p::SDOF) = p.Ampl
frequency(p::SDOF) = p.ω
period(p::SDOF) = p.T
MathematicalSystems.initial_state(p::SDOF) = p.X0
MathematicalSystems.InitialValueProblem(p::SDOF) = p.ivp
MathematicalSystems.state_matrix(p::SDOF) = state_matrix(p.ivp)

function analytic_solution(p::SDOF; A=p.Ampl, B=0.0, x0=nothing, v0=nothing)
    @assert p.X0 isa Singleton "the analytic solution requires a singleton initial condition"

    ω = p.ω
    if !isnothing(x0) && !isnothing(v0)
        A = sqrt(x0^2 + v0^2 / ω^2)
        B = atan(-v0 / (ω*x0))
    end
    return t -> A * cos(ω * t + B)
end

function analytic_derivative(p::SDOF; A=p.Ampl, B=0.0, x0=nothing, v0=nothing)
    @assert p.X0 isa Singleton "the analytic solution requires a singleton initial condition"

    ω = p.ω
    if !isnothing(x0) && !isnothing(v0)
        A = sqrt(x0^2 + v0^2 / ω^2)
        B = atan(-v0 / (ω*x0))
    end
     return t -> -ω * A * sin(ω * t + B)
end

function InitialValueProblem_quad(p::SDOF)
    M = hcat(1.0)
    ω = p.ω
    K = hcat(ω^2)
    C = zeros(1, 1)
    R = zeros(1)
    S = SecondOrderAffineContinuousSystem(M, C, K, R)
    @assert p.X0 isa Singleton
    x0 = element(p.X0)
    U0 = [x0[1]]
    U0dot = [x0[2]]
    return IVP(S, (U0, U0dot))
end

function sdof(; T=0.5,     # period
                Ampl=1.0,  # amplitude
                X0 = Singleton([Ampl, 0.0])) # initial condition
    ## frequency
    ω = 2π / T

    ## cast as a first-order system:
    ## x' = v
    ## v' = -ω^2 * x
    A = [ 0.0     1.0;
         -ω^2     0.0]

    prob = @ivp(X' = AX, X(0) ∈ X0)
    return SDOF(prob, X0, ω, Ampl, T)
end

# ================================
# Funciones para calcular PE y AD
# ================================

function consecutivos(indis)
    diffdx1 = diff(indis);
    lista_indx = [];

    j = 1
    while j <= length(indis)
        inds = [indis[j]]
        while (j < length(indis)) && (diffdx1[j]==1)
            j+=1
            push!(inds, indis[j])
        end
        j+=1
        push!(lista_indx, inds)
    end
    return lista_indx
end

function period_evolution2(T, tf, α; X0 = Singleton([1.0, 0.0]), polardirs = PolarDirections(200))

    ω = 2π / T

    ## x' = v
    ## v' = -ω^2 * x
    A = [0.0     1.0;
         -ω^2    0.0]

    prob = @ivp(X' = AX, X(0) ∈ X0)
    guard = HPolyhedron([HalfSpace([ 0.0,  1.0], 0.0),    # v ≤ 0
                         HalfSpace([ 0.0, -1.0], 0.0),    # v ≥ 0
                         HalfSpace([-1.0,  0.0], 0.0)])   # x ≥ 0

    alg(αi) = VREP(δ=αi*T, approx_model=StepIntersect(setops=:concrete))
    ## alg(αi) = GLGM06(δ=αi*T))
    ## alg(αi) = LGG09(δ=αi*T, template=polardirs))

    conjuntos = []
    for i in 1:length(α)
        sol = solve(prob, tspan=(0.0, tf), alg=alg(α[i]))
        ind_i = findall(R -> !isdisjoint(set(R), guard), array(sol)) #TODO sol[2:end]
        push!(conjuntos, [sol[x] for x in consecutivos(ind_i[2:end])])
    end

    return conjuntos
end

function AD_numerica(xb, N, npers)
    ts_barrido = Int(N); #ultimos N periodos
    return (1-(maximum(xb[end-ts_barrido:end]))^(1/(Int(npers)-N)))*100;
end

function compute_amplitude_decay(T, N, numpers, α)
    tf = numpers*T
    conjs = period_evolution2(T, tf, α)

    ## intervalos de tiempo en los que se da una vuelta
    Ts = [[tspan(R) for R in c] for c in conjs]

    ## divido entre el numero de vueltas hechas para tener los periodos
    periodos = [[R/i for (i, R) in enumerate(c)] for c in Ts]

    amplitudes = [[overapproximate(project(overapproximate(ConvexHullArray(set.(R)), Hyperrectangle), [1]), Interval).dat for R in c] for c in conjs]
    amps_min = [inf.(a1) for a1 in amplitudes]
    amps_min_num = AD_numerica.(amps_min, N, numpers)

    return amps_min_num
end

function compute_period_elongation(T, N, numpers, α)

    tf = numpers*T

    conjs = period_evolution2(T, tf, α)

    Ts = [[tspan(R) for R in c] for c in conjs] #intervalos de tiempo en los que se da una vuelta
    periodos = [[R/i for (i, R) in enumerate(c)] for c in Ts]

    intersecciones = [];
    for i in 1:length(periodos)
        z = [];
        for j in 1:length(periodos[i])
            push!(z, reduce(∩, periodos[i][1:j]))
        end
        push!(intersecciones, z)
    end

    ## ¿Todos los periodos limite contienen al T del sistema?
    @assert all([T ∈ interss[end] for interss in intersecciones])

    x = [Singleton([alphai]) × Interval((Ai[end]-T)/T*100) for (alphai, Ai) in zip(α, intersecciones)];
    y1 = UnionSetArray([ConvexHull(x[k], x[k+1]) for k in 1:length(x)-1]);

    PE_RFEM_sup = ([sup(Ai[end]) for Ai in intersecciones].-T)./T*100;
    ##@. PE_RFEM_sup = (PE_RFEM_sup/T - 1) * 100

    PE_RFEM_inf = ([inf(Ai[end]) for Ai in intersecciones].-T)./T*100;

    return PE_RFEM_inf, PE_RFEM_sup
end

max_stars_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_issues_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_forks_repo_name: JuliaReach/SetPropagation-FEM-Examples


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

# What is included in this file:
# - the `fmi3ComponentState`-enum which mirrors the internal FMU state (state-machine, not the system state)
# - the `FMU3Component`-struct
# - the `FMU3`-struct
# - string/enum-converters for FMI-attribute-structs (e.g. `fmi3StatusToString`, ...)

# this is a custom type to catch the internal state of the component
@enum fmi3ComponentState begin
    # ToDo
    fmi3ComponentToDo
end

"""
Source: FMISpec3.0, Version D5ef1c1:: 2.2.1. Header Files and Naming of Functions

The mutable struct represents a pointer to an FMU specific data structure that contains the information needed to process the model equations or to process the co-simulation of the model/subsystem represented by the FMU.
"""
mutable struct FMU3Component
    compAddr::Ptr{Nothing}
    fmu
    state

    loggingOn::Bool
    instanceName::String
    continuousStatesChanged::fmi3Boolean

    t::fmi3Float64             # the system time
    #x::Array{fmi3Float64, 1}   # the system states (or sometimes u)
    #ẋ::Array{fmi3Float64, 1}   # the system state derivative (or sometimes u̇)
    #u::Array{fmi3Float64, 1}   # the system inputs
    #y::Array{fmi3Float64, 1}   # the system outputs
    #p::Array{fmi3Float64, 1}   # the system parameters
    z::Array{fmi3Float64, 1}   # the system event indicators
    z_prev::Union{Nothing, Array{fmi3Float64, 1}}   # the last system event indicators

    realValues::Dict

    x_vrs::Array{fmi3ValueReference, 1}   # the system state value references
    ẋ_vrs::Array{fmi3ValueReference, 1}   # the system state derivative value references
    u_vrs::Array{fmi3ValueReference, 1}   # the system input value references
    y_vrs::Array{fmi3ValueReference, 1}   # the system output value references
    p_vrs::Array{fmi3ValueReference, 1}   # the system parameter value references

    # FMIFlux

    jac_dxy_x::Matrix{fmi3Float64}
    jac_dxy_u::Matrix{fmi3Float64}
    jac_x::Array{fmi3Float64}
    jac_t::fmi3Float64
    skipNextDoStep::Bool    # allows skipping the next `fmi2DoStep` like it is not called

    # custom

    rootsFound::Array{fmi3Int32}
    stateEvent::fmi3Boolean
    timeEvent::fmi3Boolean
    stepEvent::fmi3Boolean

    function fmi3Component(compAddr, fmu)
        inst = new()
        inst.compAddr = compAddr
        inst.fmu = fmu
        inst.t = 0.0

        inst.z_prev = nothing

        inst.realValues = Dict()
        inst.x_vrs = Array{fmi3ValueReference, 1}()
        inst.ẋ_vrs = Array{fmi3ValueReference, 1}()
        inst.u_vrs = Array{fmi3ValueReference, 1}()
        inst.y_vrs = Array{fmi3ValueReference, 1}()
        inst.p_vrs = Array{fmi3ValueReference, 1}()

        # initialize further variables
        inst.jac_x = Array{fmi3Float64, 1}()
        inst.jac_t = -1.0
        inst.jac_dxy_x = zeros(fmi3Float64, 0, 0)
        inst.jac_dxy_u = zeros(fmi3Float64, 0, 0)
        inst.skipNextDoStep = false

        return inst
    end
end

""" Overload the Base.show() function for custom printing of the FMU2Component"""
Base.show(io::IO, fmu::FMU3Component) = print(io,
"FMU2:         $(fmu.fmu)
State:         $(fmu.state)
Logging:       $(fmu.loggingOn)
Instance name: $(fmu.instanceName)
System time:   $(fmu.t)
Values:        $(fmu.realValues)"
)

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.1. Header Files and Naming of Functions

The mutable struct representing an FMU in the FMI 3.0 Standard.
Also contains the paths to the FMU and ZIP folder as well als all the FMI 3.0 function pointers
"""
mutable struct FMU3 <: FMU
    modelName::String
    instanceName::String
    fmuResourceLocation::String

    modelDescription::fmi3ModelDescription

    type::fmi3Type
    instanceEnvironment::fmi3InstanceEnvironment
    components::Array # {fmi3Component}

    # c-functions
    cInstantiateModelExchange::Ptr{Cvoid}
    cInstantiateCoSimulation::Ptr{Cvoid}
    cInstantiateScheduledExecution::Ptr{Cvoid}

    cGetVersion::Ptr{Cvoid}
    cFreeInstance::Ptr{Cvoid}
    cSetDebugLogging::Ptr{Cvoid}
    cEnterConfigurationMode::Ptr{Cvoid}
    cExitConfigurationMode::Ptr{Cvoid}
    cEnterInitializationMode::Ptr{Cvoid}
    cExitInitializationMode::Ptr{Cvoid}
    cTerminate::Ptr{Cvoid}
    cReset::Ptr{Cvoid}
    cGetFloat32::Ptr{Cvoid}
    cSetFloat32::Ptr{Cvoid}
    cGetFloat64::Ptr{Cvoid}
    cSetFloat64::Ptr{Cvoid}
    cGetInt8::Ptr{Cvoid}
    cSetInt8::Ptr{Cvoid}
    cGetUInt8::Ptr{Cvoid}
    cSetUInt8::Ptr{Cvoid}
    cGetInt16::Ptr{Cvoid}
    cSetInt16::Ptr{Cvoid}
    cGetUInt16::Ptr{Cvoid}
    cSetUInt16::Ptr{Cvoid}
    cGetInt32::Ptr{Cvoid}
    cSetInt32::Ptr{Cvoid}
    cGetUInt32::Ptr{Cvoid}
    cSetUInt32::Ptr{Cvoid}
    cGetInt64::Ptr{Cvoid}
    cSetInt64::Ptr{Cvoid}
    cGetUInt64::Ptr{Cvoid}
    cSetUInt64::Ptr{Cvoid}
    cGetBoolean::Ptr{Cvoid}
    cSetBoolean::Ptr{Cvoid}
    cGetString::Ptr{Cvoid}
    cSetString::Ptr{Cvoid}
    cGetBinary::Ptr{Cvoid}
    cSetBinary::Ptr{Cvoid}
    cGetFMUState::Ptr{Cvoid}
    cSetFMUState::Ptr{Cvoid}
    cFreeFMUState::Ptr{Cvoid}
    cSerializedFMUStateSize::Ptr{Cvoid}
    cSerializeFMUState::Ptr{Cvoid}
    cDeSerializeFMUState::Ptr{Cvoid}
    cGetDirectionalDerivative::Ptr{Cvoid}
    cGetAdjointDerivative::Ptr{Cvoid}
    cEvaluateDiscreteStates::Ptr{Cvoid}
    cGetNumberOfVariableDependencies::Ptr{Cvoid}
    cGetVariableDependencies::Ptr{Cvoid}

    # Co Simulation function calls
    cGetOutputDerivatives::Ptr{Cvoid}
    cEnterStepMode::Ptr{Cvoid}
    cDoStep::Ptr{Cvoid}

    # Model Exchange function calls
    cGetNumberOfContinuousStates::Ptr{Cvoid}
    cGetNumberOfEventIndicators::Ptr{Cvoid}
    cGetContinuousStates::Ptr{Cvoid}
    cGetNominalsOfContinuousStates::Ptr{Cvoid}
    cEnterContinuousTimeMode::Ptr{Cvoid}
    cSetTime::Ptr{Cvoid}
    cSetContinuousStates::Ptr{Cvoid}
    cGetContinuousStateDerivatives::Ptr{Cvoid}
    cGetEventIndicators::Ptr{Cvoid}
    cCompletedIntegratorStep::Ptr{Cvoid}
    cEnterEventMode::Ptr{Cvoid}
    cUpdateDiscreteStates::Ptr{Cvoid}

    # Scheduled Execution function calls
    cSetIntervalDecimal::Ptr{Cvoid}
    cSetIntervalFraction::Ptr{Cvoid}
    cGetIntervalDecimal::Ptr{Cvoid}
    cGetIntervalFraction::Ptr{Cvoid}
    cGetShiftDecimal::Ptr{Cvoid}
    cGetShiftFraction::Ptr{Cvoid}
    cActivateModelPartition::Ptr{Cvoid}

    # paths of ziped and unziped FMU folders
    path::String
    zipPath::String

    # c-libraries
    libHandle::Ptr{Nothing}

    # START: experimental section (to FMIFlux.jl)
    dependencies::Matrix

    # linearization jacobians
    A::Matrix{fmi2Real}
    B::Matrix{fmi2Real}
    C::Matrix{fmi2Real}
    D::Matrix{fmi2Real}

    # END: experimental section

    # Constructor
    function FMU3()
        inst = new()
        inst.components = []
        return inst
    end
end

""" Overload the Base.show() function for custom printing of the FMU3"""
Base.show(io::IO, fmu::FMU3) = print(io,
"Model name:       $(fmu.modelName)
Instance name:     $(fmu.instanceName)
Model description: $(fmu.modelDescription)
Type:              $(fmu.type)
Components:        $(fmu.components)"
)

"""
Formats the fmi3Status into a String.
"""
function fmi3StatusToString(status::fmi3Status)
    if status == fmi3StatusOK
        return "OK"
    elseif status == fmi3StatusWarning
        return "Warning"
    elseif status == fmi3StatusDiscard
        return "Discard"
    elseif status == fmi3StatusError
        return "Error"
    elseif status == fmi3StatusFatal
        return "Fatal"
    else
        return "Unknown"
    end
end

"""
Formats the fmi3Status into a Integer.
"""
function fmi3StatusToString(status::Integer)
    if status == fmi3StatusOK
        return "OK"
    elseif status == fmi3StatusWarning
        return "Warning"
    elseif status == fmi3StatusDiscard
        return "Discard"
    elseif status == fmi3StatusError
        return "Error"
    elseif status == fmi3StatusFatal
        return "Fatal"
    else
        return "Unknown"
    end
end

# ToDo: Equivalent of fmi2CausalityToString, fmi2StringToCausality, ...


text:
is_directed(::AdjacencyMap{D}) where {D} = D
is_directed(::Type{<:AdjacencyMap{D}}) where {D} = D

###
### has_edge
###
has_edge(g::AdjacencyMap, e::AbstractEdge) = has_edge(g, src(e), dst(e))
function has_edge(g::AdjMap, i::Integer, j::Integer)
    if i < j
        if j <= nvertices(g)
            unsafe_has_edge(getcolptr(g), rowvals(g), i, j)
        else
            return false
        end
    else ## j <= i
        if i <= nvertices(g)
            unsafe_has_edge(getcolptr(g), rowvals(g), j, i)
        else
            return false
        end
    end
end
function has_edge(g::AdjDiMap, i::Integer, j::Integer)
    N = nvertices(g)
    if !((1 <= i <= N) & (1 <= j <= N))
        return false
    else
        unsafe_has_edge(getcolptr(g), rowvals(g), i, j)
    end
end
function unsafe_has_edge(cptr, rvals, i, j)
    coljlastk = @inbounds(cptr[j + 1] - 1)
    searchk = searchsortedfirst(rvals, i, @inbounds(cptr[j]), coljlastk, Base.Order.Forward)
    return searchk <= coljlastk && @inbounds(rvals[searchk]) == i
end

###
### add_edge!
###
add_edge!(g::AdjacencyMap, i::Int, j::Int) = add_edge!(g, i, j, no_metadata)
function add_edge!(g::AdjMap, i::Int, j::Int, val) where {T,M}
    N = nvertices(g)
    if !((1 <= i <= N) & (1 <= j <= N))
        return false
    elseif i < j
        unsafe_add_edge!(g, i, j, val)
    else
        unsafe_add_edge!(g, j, i, val)
    end
end
function add_edge!(g::AdjDiMap, i::Int, j::Int, val) where {T,M}
    N = nvertices(g)
    if !((1 <= i <= N) & (1 <= j <= N))
        return false
    else
        return unsafe_add_edge!(g, i, j, val)
    end
end
function unsafe_add_edge!(g, i, j, val)
    coljlastk = @inbounds(getcolptr(g)[j + 1] - 1)
    searchk = searchsortedfirst(rowvals(g), i, @inbounds(getcolptr(g)[j]), coljlastk, Base.Order.Forward)
    if searchk <= coljlastk && rowvals(g)[searchk] == i
        _setindex!(metadata(g), val, searchk)
        return true
    else
        N = nvertices(g)
        nz = getcolptr(g)[N + 1]
        # if nnz(A) < length(rowval/nzval): no need to grow rowval and preserve values
        _insert!(rowvals(g), searchk, i, nz)
        _insert!(metadata(g), searchk, val, nz)
        @simd for i in (j + 1):(N + 1)
            @inbounds getcolptr(g)[i] += 1
        end
    end
    return true
end
# insert item at position pos, shifting only from pos+1 to nz
_insert!(m, pos, ::NoMetadata, nz) = error()
_insert!(::NoMetadata, pos, ::NoMetadata, nz) = nothing
function _insert!(m, pos, val, nz)
    if nz > length(m)
        insert!(m, pos, val)
    else # nz < length(v)
        Base.unsafe_copyto!(m, pos + 1, m, pos, nz - pos)
        @inbounds(setindex!(m, val, pos))
        m
    end
end
_insert!(m::ListOfMeta, pos, val, nz) = _insert!(parent(m), pos, val, nz)

_setindex!(nzvals, val, i) = @inbounds(setindex!(nzvals, val, i))
_setindex!(m::ListOfMeta, val, i) = _setindex!(parent(m), val, i)
_setindex!(::NoMetadata, ::NoMetadata, i) = nothing
_setindex!(::NoMetadata, val, i) = nothing
_setindex!(nzvals, ::NoMetadata, i) = nothing

###
### rem_edge!
###
rem_edge!(g::AdjDiMap, i::Integer, j::Integer) = unsafe_rem_edge!(g, i, j)
function rem_edge!(g::AdjMap, i::Integer, j::Integer)
    if i < j
        return unsafe_rem_edge!(g, i, j)
    else
        return unsafe_rem_edge!(g, j, i)
    end
end
function unsafe_rem_edge!(g::AdjacencyMap, i::Integer, j::Integer)
    N = nvertices(g)
    if !((1 <= i <= N) && (1 <= j <= N))
        return false
    else
        cptr = getcolptr(g)
        rval = rowvals(g)
        coljlastk = @inbounds(cptr[j + 1]) - 1
        searchk = searchsortedfirst(rval, i, @inbounds(cptr[j]), coljlastk, Base.Order.Forward)
        if searchk > coljlastk || @inbounds(rval[searchk]) != i
            return false
        else
            Awritepos = 1
            oldAcolptrAj = 1
            @inbounds for Aj in OneTo(N)
                for Ak in oldAcolptrAj:(cptr[Aj + 1]-1)
                    if Ak != searchk && Awritepos != Ak
                        # If this element should be kept, rewrite in new position
                        _move_value!(rval, Awritepos, Ak)
                        _move_value!(metadata(g), Awritepos, Ak)
                    end
                    Awritepos += 1
                end
                oldAcolptrAj = Acolptr[Aj + 1]
                Acolptr[Aj + 1] = Awritepos
            end

            # Trim A's storage if necessary
            Annz = Acolptr[end] - 1
            resize!(rval, Annz)
            resize!(m, Annz)
            return true
        end
    end
end

_move_value!(m, d, s) = @inbounds(setindex!(m, m[s], d))
_move_value!(::NoMetadata, d, s) = nothing

function SparseArrays.fkeep!(A::AdjacencyMap, f)
    An = size(A, 2)
    Acolptr = getcolptr(A)
    Arowval = rowvals(A)
    Anzval = metadata(A)

    # Sweep through columns, rewriting kept elements in their new positions
    # and updating the column pointers accordingly as we go.
    Awritepos = 1
    oldAcolptrAj = 1
    @inbounds for Aj in 1:An
        for Ak in oldAcolptrAj:(Acolptr[Aj+1]-1)
            Ai = Arowval[Ak]
            Ax = Anzval[Ak]
            # If this element should be kept, rewrite in new position
            if f(Ai, Aj, Ax)
                if Awritepos != Ak
                    _move_value!(rval, Awritepos, Ak)
                    _move_value!(Anzval, Awritepos, Ak)
                end
                Awritepos += 1
            end
        end
        oldAcolptrAj = Acolptr[Aj+1]
        Acolptr[Aj+1] = Awritepos
    end

    # Trim A's storage if necessary
    Annz = Acolptr[end] - 1
    resize!(Arowval, Annz)
    resize!(Anzval, Annz)

    return A
end

###
### add_vertices!
###
function add_vertices!(g::AdjacencyMap, n::Integer)
    if can_change_size(g)
        return _unsafe_add_vertices!(g, max(0, n))
    else
        return false
    end
end

function _unsafe_add_vertices!(g::AdjacencyMap, n)
    cptr = getcolptr(g)
    c = last(cptr)
    for _ in OneTo(n)
        push!(cptr, c)
    end
    return true
end

function add_vertex!(g::AdjacencyMap)
    cptr = getcolptr(g)
    push!(cptr, last(cptr))
    return true
end

###
### add_vertices!
###
# TODO LightGraphs.rem_vertex!(m::AdjMap, n::Integer)
# TODO LightGraphs.rem_vertices!(m::AdjMap, n::AbstractVector)
# TODO check this a bunch
function rem_vertex!(g::AdjacencyMap, n::Integer)
    fkeep!(g, (i, j, x) -> (i != n && j != n))
    deleteat!(getcolptr(g), n)
end


text:
@testset "All Tests" begin
  include("core/RUNME.jl")
  @test_skip "Stdlib Tests"
  @test_skip "External Tests"
  @test_skip "Behavior Tests"
  @test_skip "Perf Tests"
end


text:
import ..SamplingConstraints: SamplingConstraint

"""
    ConstrainedResampling(constraint::NTuple{N, Union{SamplingConstraint, Vector{<:SamplingConstraint}}})

Indicates that resampling should be performed with constraints on the uncertain values.
"""
struct ConstrainedResampling{N} <: AbstractUncertainDataResampling where N
    constraints::NTuple{N, Union{SamplingConstraint, Vector{<:SamplingConstraint}}}
end

ConstrainedResampling(constraints::Union{SamplingConstraint, Vector{<:SamplingConstraint}}...) =
    ConstrainedResampling(constraints)

ConstrainedResampling(constraints::Vector{Union{SamplingConstraint, Vector{<:SamplingConstraint}}}) =
    ConstrainedResampling((constraints...))

Base.length(resampling::ConstrainedResampling) = length(resampling.constraints)
Base.getindex(resampling::ConstrainedResampling, i) = resampling.constraints[i]
Base.firstindex(resampling::ConstrainedResampling) = 1
Base.lastindex(resampling::ConstrainedResampling) = length(resampling)

function summarise(resampling::ConstrainedResampling)
    _type = typeof(resampling)
    constraint_types = [typeof(c) for c in resampling.constraints]

    strs = ["$constraint" for  constraint in constraint_types]
    return "$_type" * "(" * join(strs, ", ") * ")"
end

Base.show(io::IO, resampling::ConstrainedResampling) = print(io, summarise(resampling))

export ConstrainedResampling