text:
function get_quadruped(;
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=0.8,
    spring=0.0,
    damper=0.0,
    contact_feet=true,
    contact_body=true,
    limits=true,
    spring_offset=true,
    path=joinpath(@__DIR__, "../deps/quadruped.urdf"),
    joint_limits=[[-0.5, -0.5, -2.5,],
                  [ 0.5,  1.5, -1.0,]],
    T=Float64)

    mech = Mechanism(path, true, T,
        gravity=gravity,
        timestep=timestep,
        spring=spring,
        damper=damper)

    # Adding springs and dampers
    for joint in mech.joints[2:end]
        joint.damper = true
        joint.spring = true
        joint.translational.spring=spring
        joint.translational.damper=damper
        joint.rotational.spring=spring
        joint.rotational.damper=damper
    end

    # joint limits
    joints = deepcopy(mech.joints)
    if limits
        for group in [:FR, :FL, :RR, :RL]
            hip_joint = get_joint(mech, Symbol(group, :_hip_joint))
            joints[hip_joint.id] = add_limits(mech, hip_joint,
                rot_limits=[SVector{1}(joint_limits[1][1]), SVector{1}(joint_limits[2][1])])

            thigh_joint = get_joint(mech, Symbol(group, :_thigh_joint))
            joints[thigh_joint.id] = add_limits(mech, thigh_joint,
                rot_limits=[SVector{1}(joint_limits[1][2]), SVector{1}(joint_limits[2][2])])

            calf_joint = get_joint(mech, Symbol(group, :_calf_joint))
            joints[calf_joint.id] = add_limits(mech, calf_joint,
                rot_limits=[SVector{1}(joint_limits[1][3]), SVector{1}(joint_limits[2][3])])
        end
        mech = Mechanism(Origin{T}(), [mech.bodies...], [joints...],
            gravity=gravity,
            timestep=timestep,
            spring=spring,
            damper=damper)
    end

    origin = Origin{T}()
    bodies = mech.bodies
    eqs = mech.joints

    contacts = ContactConstraint{T}[]
    if contact_feet
        # Foot contact
        normal = [0.0; 0.0; 1.0]
        foot_contact = [-0.006; 0.0; -0.092]
        foot_contact_radius = 0.021
        hip_contact = [0.0; 0.05; 0.0]
        hip_contact_radius = 0.05
        elbow_contactR = [-0.005; -0.023; -0.16]
        elbow_contactL = [-0.005; 0.023; -0.16]
        elbow_contact_radius = 0.023

        foot_contacts1 = contact_constraint(get_body(mech,:FR_calf), normal;
            friction_coefficient=friction_coefficient,
            contact_origin=foot_contact,
            contact_radius=foot_contact_radius,
            name=:FR_contact)
        foot_contacts2 = contact_constraint(get_body(mech,:FL_calf), normal;
            friction_coefficient=friction_coefficient,
            contact_origin=foot_contact,
            contact_radius=foot_contact_radius,
            name=:FL_contact)
        foot_contacts3 = contact_constraint(get_body(mech,:RR_calf), normal;
            friction_coefficient=friction_coefficient,
            contact_origin=foot_contact,
            contact_radius=foot_contact_radius,
            name=:RR_contact)
        foot_contacts4 = contact_constraint(get_body(mech,:RL_calf), normal;
            friction_coefficient=friction_coefficient,
            contact_origin=foot_contact,
            contact_radius=foot_contact_radius,
            name=:RL_contact)
        contacts = [contacts...; foot_contacts1; foot_contacts2; foot_contacts3; foot_contacts4]

        if contact_body
            elbow_contacts1 = contact_constraint(get_body(mech,:FR_thigh), normal;
                friction_coefficient=friction_coefficient,
                contact_origin=elbow_contactR,
                contact_radius=elbow_contact_radius,
                name=:FR_hip_contact)
            elbow_contacts2 = contact_constraint(get_body(mech,:FL_thigh), normal;
                friction_coefficient=friction_coefficient,
                contact_origin=elbow_contactL,
                contact_radius=elbow_contact_radius,
                name=:FL_hip_contact)
            elbow_contacts3 = contact_constraint(get_body(mech,:RR_thigh), normal;
                friction_coefficient=friction_coefficient,
                contact_origin=elbow_contactR,
                contact_radius=elbow_contact_radius,
                name=:RR_hip_contact)
            elbow_contacts4 = contact_constraint(get_body(mech,:RL_thigh), normal;
                friction_coefficient=friction_coefficient,
                contact_origin=elbow_contactL,
                contact_radius=elbow_contact_radius,
                name=:RL_hip_contact)
            push!(contacts, elbow_contacts1, elbow_contacts2, elbow_contacts3, elbow_contacts4)
            hip_contacts1 = contact_constraint(get_body(mech,:FR_hip), normal;
                friction_coefficient=friction_coefficient,
                contact_origin=-hip_contact,
                contact_radius=hip_contact_radius,
                name=:FR_hip_contact)
            hip_contacts2 = contact_constraint(get_body(mech,:FL_hip), normal;
                friction_coefficient=friction_coefficient,
                contact_origin=+hip_contact,
                contact_radius=hip_contact_radius,
                name=:FL_hip_contact)
            hip_contacts3 = contact_constraint(get_body(mech,:RR_hip), normal;
                friction_coefficient=friction_coefficient,
                contact_origin=-hip_contact,
                contact_radius=hip_contact_radius,
                name=:RR_hip_contact)
            hip_contacts4 = contact_constraint(get_body(mech,:RL_hip), normal;
                friction_coefficient=friction_coefficient,
                contact_origin=+hip_contact,
                contact_radius=hip_contact_radius,
                name=:RL_hip_contact)
            push!(contacts, hip_contacts1, hip_contacts2, hip_contacts3, hip_contacts4)
        end

        set_minimal_coordinates!(mech, get_joint(mech, :floating_base), [0.0; 0.0; 0.32; 0.0; 0.0; 0.0])

        mech = Mechanism(origin, bodies, eqs, contacts,
            gravity=gravity,
            timestep=timestep,
            spring=spring,
            damper=damper)
    end

    # Spring Offset
    if spring_offset
        θ_thigh = 0.9
        θ_calf = -1.425
        get_node(mech, :FR_hip_joint).rotational.spring_offset=0.0*sones(1)
        get_node(mech, :FL_hip_joint).rotational.spring_offset=0.0*sones(1)
        get_node(mech, :RR_hip_joint).rotational.spring_offset=0.0*sones(1)
        get_node(mech, :RL_hip_joint).rotational.spring_offset=0.0*sones(1)

        get_node(mech, :FR_thigh_joint).rotational.spring_offset=θ_thigh*sones(1)
        get_node(mech, :FL_thigh_joint).rotational.spring_offset=θ_thigh*sones(1)
        get_node(mech, :RR_thigh_joint).rotational.spring_offset=θ_thigh*sones(1)
        get_node(mech, :RL_thigh_joint).rotational.spring_offset=θ_thigh*sones(1)

        get_node(mech, :FR_calf_joint).rotational.spring_offset=θ_calf*sones(1)
        get_node(mech, :FL_calf_joint).rotational.spring_offset=θ_calf*sones(1)
        get_node(mech, :RR_calf_joint).rotational.spring_offset=θ_calf*sones(1)
        get_node(mech, :RL_calf_joint).rotational.spring_offset=θ_calf*sones(1)
    end
    return mech
end

function initialize_quadruped!(mechanism::Mechanism{T};
    body_position=[0.0, 0.0, 0.0],
    body_orientation=[0.0, 0.0, 0.0],
    body_linear_velocity=[0.0, 0.0, 0.0],
    body_angular_velocity=[0.0, 0.0, 0.0],
    link_angle=0.95) where T
    #TODO: add more leg angles

    body_position += [0.0, 0.0, 0.32]
    zero_velocity!(mechanism)
    set_minimal_coordinates!(mechanism, get_joint(mechanism, :floating_base),
        [body_position; body_orientation])

    set_minimal_coordinates!(mechanism, get_joint(mechanism, :FR_hip_joint), [0.0])
    set_minimal_coordinates!(mechanism, get_joint(mechanism, :FR_thigh_joint), [link_angle])
    set_minimal_coordinates!(mechanism, get_joint(mechanism, :FR_calf_joint), [-1.5 * link_angle])

    set_minimal_coordinates!(mechanism, get_joint(mechanism, :FL_hip_joint), [0.0])
    set_minimal_coordinates!(mechanism, get_joint(mechanism, :FL_thigh_joint), [link_angle * 0.9])
    set_minimal_coordinates!(mechanism, get_joint(mechanism, :FL_calf_joint), [-1.5 * link_angle])

    set_minimal_coordinates!(mechanism, get_joint(mechanism, :RR_hip_joint), [0.0])
    set_minimal_coordinates!(mechanism, get_joint(mechanism, :RR_thigh_joint), [link_angle * 0.9])
    set_minimal_coordinates!(mechanism, get_joint(mechanism, :RR_calf_joint), [-1.5 * link_angle])

    set_minimal_coordinates!(mechanism, get_joint(mechanism, :RL_hip_joint), [0.0])
    set_minimal_coordinates!(mechanism, get_joint(mechanism, :RL_thigh_joint), [link_angle])
    set_minimal_coordinates!(mechanism, get_joint(mechanism, :RL_calf_joint), [-1.5 * link_angle])

    set_minimal_velocities!(mechanism, get_joint(mechanism, :floating_base),
        [body_linear_velocity; body_angular_velocity])

    set_minimal_velocities!(mechanism, get_joint(mechanism, :FR_hip_joint), [0.0])
    set_minimal_velocities!(mechanism, get_joint(mechanism, :FR_thigh_joint), [0.0])
    set_minimal_velocities!(mechanism, get_joint(mechanism, :FR_calf_joint), [0.0])

    set_minimal_velocities!(mechanism, get_joint(mechanism, :FL_hip_joint), [0.0])
    set_minimal_velocities!(mechanism, get_joint(mechanism, :FL_thigh_joint), [0.0])
    set_minimal_velocities!(mechanism, get_joint(mechanism, :FL_calf_joint), [0.0])

    set_minimal_velocities!(mechanism, get_joint(mechanism, :RR_hip_joint), [0.0])
    set_minimal_velocities!(mechanism, get_joint(mechanism, :RR_thigh_joint), [0.0])
    set_minimal_velocities!(mechanism, get_joint(mechanism, :RR_calf_joint), [0.0])

    set_minimal_velocities!(mechanism, get_joint(mechanism, :RL_hip_joint), [0.0])
    set_minimal_velocities!(mechanism, get_joint(mechanism, :RL_thigh_joint), [0.0])
    set_minimal_velocities!(mechanism, get_joint(mechanism, :RL_calf_joint), [0.0])
end


text:
# Will be evaluated _inside_ the workspace process.

# Pluto does most things on process 1 (the server), and it uses little workspace processes to evaluate notebook code in.
# These baby processes don't import Pluto, they only import this module. Functions from this module are called by WorkspaceManager.jl, using Distributed

# So when reading this file, pretend that you are living in process 2, and you are communicating with Pluto's server, who lives in process 1.

module PlutoRunner

using Markdown
import Markdown: html, htmlinline, LaTeX, withtag, htmlesc
import Distributed
import Base64
import REPL.REPLCompletions: completions, complete_path, completion_text, Completion, ModuleCompletion
import Base: show, istextmime
import UUIDs: UUID
import Logging

export @bind

MimedOutput = Tuple{Union{String,Vector{UInt8}}, MIME}

###
# WORKSPACE MANAGER
###

#Will be set to the latest workspace module
"The current workspace where your variables live. See [`move_vars`](@ref)."
current_module = Main

function set_current_module(newname)
    # Revise.jl support
    if isdefined(current_module, :Revise) &&
        isdefined(current_module.Revise, :revise) && current_module.Revise.revise isa Function &&
        isdefined(current_module.Revise, :revision_queue) && current_module.Revise.revision_queue isa AbstractSet

        if !isempty(current_module.Revise.revision_queue) # to avoid the sleep(0.01) in revise()
            current_module.Revise.revise()
        end
    end

    global iocontext = IOContext(iocontext, :module => current_module)
    global iocontext_compact = IOContext(iocontext_compact, :module => current_module)

    global current_module = getfield(Main, newname)
end

const cell_results = Dict{UUID, WeakRef}()

function formatted_result_of(id::UUID, ends_with_semicolon::Bool)::NamedTuple{(:output_formatted, :errored, :interrupted, :runtime),Tuple{MimedOutput,Bool,Bool,Union{UInt64, Missing}}}
    ans = cell_results[id].value
    errored = ans isa CapturedException
    output_formatted = (!ends_with_semicolon || errored) ? format_output(ans) : ("", MIME"text/plain"())
    (output_formatted = output_formatted, errored = errored, interrupted = false, runtime = Main.runtime)
end

"""
Move some of the globals over from one workspace to another. This is how Pluto "deletes" globals - it doesn't, it just executes your new code in a new module where those globals are not defined.

Notebook code does run in `Main` - it runs in workspace modules. Every time that you run cells, a new module is created, called `Main.workspace123` with `123` an increasing number.

The trick boils down to two things:
1. When we create a new workspace module, we move over some of the global from the old workspace. (But not the ones that we want to 'delete'!)
2. If a function used to be defined, but now we want to delete it, then we go through the method table of that function and snoop out all methods that we defined by us, and not by another package. This is how we reverse extending external functions. For example, if you run a cell with `Base.sqrt(s::String) = "the square root of" * s`, and then delete that cell, then you can still call `sqrt(1)` but `sqrt("one")` will err. Cool right!
"""
function move_vars(old_workspace_name::Symbol, new_workspace_name::Symbol, vars_to_delete::Set{Symbol}, funcs_to_delete::Set{Tuple{UUID,Vector{Symbol}}}, module_imports_to_move::Set{Expr})
    old_workspace = getfield(Main, old_workspace_name)
    new_workspace = getfield(Main, new_workspace_name)

    for expr in module_imports_to_move
        try
            Core.eval(new_workspace, expr)
        catch; end # TODO catch specificallly
    end

    # TODO: delete
    Core.eval(new_workspace, :(import ..($(old_workspace_name))))

    old_names = names(old_workspace, all=true, imported=true)

    funcs_with_no_methods_left = filter(funcs_to_delete) do f
        !try_delete_toplevel_methods(old_workspace, f)
    end
    name_symbols_of_funcs_with_no_methods_left = last.(last.(funcs_with_no_methods_left))
    for symbol in old_names
        if (symbol ∈ vars_to_delete) || (symbol ∈ name_symbols_of_funcs_with_no_methods_left)
            # var will be redefined - unreference the value so that GC can snoop it

            # free memory for other variables
            # & delete methods created in the old module:
            # for example, the old module might extend an imported function:
            # `import Base: show; show(io::IO, x::Flower) = print(io, "🌷")`
            # when you delete/change this cell, you want this extension to disappear.
            if isdefined(old_workspace, symbol)
                # try_delete_toplevel_methods(old_workspace, symbol)

                try
                    # it could be that `symbol ∈ vars_to_move`, but the _value_ has already been moved to the new reference in `new_module`.
                    # so clearing the value of this reference does not affect the reference in `new_workspace`.
                    Core.eval(old_workspace, :($(symbol) = nothing))
                catch; end # sometimes impossible, eg. when $symbol was constant
            end
        else
            # var will not be redefined in the new workspace, move it over
            if !(symbol == :eval || symbol == :include || string(symbol)[1] == '#' || startswith(string(symbol), "workspace"))
                try
                    val = getfield(old_workspace, symbol)

                    # Expose the variable in the scope of `new_workspace`
                    Core.eval(new_workspace, :(import ..($(old_workspace_name)).$(symbol)))
                catch ex
                    if !(ex isa UndefVarError)
                        @warn "Failed to move variable $(symbol) to new workspace:"
                        showerror(stderr, ex, stacktrace(catch_backtrace()))
                    end
                end
            end
        end
    end
end

"Return whether the `method` was defined inside this notebook, and not in external code."
isfromcell(method::Method, cell_id::UUID) = endswith(String(method.file), string(cell_id))

"Delete all methods of `f` that were defined in this notebook, and leave the ones defined in other packages, base, etc. ✂

Return whether the function has any methods left after deletion."
function delete_toplevel_methods(f::Function, cell_id::UUID)::Bool
    # we can delete methods of functions!
    # instead of deleting all methods, we only delete methods that were defined in this notebook. This is necessary when the notebook code extends a function from remote code
    methods_table = typeof(f).name.mt
    deleted_sigs = Set{Type}()
    Base.visit(methods_table) do method # iterates through all methods of `f`, including overridden ones
        if isfromcell(method, cell_id) && getfield(method, deleted_world) == alive_world_val
            Base.delete_method(method)
            push!(deleted_sigs, method.sig)
        end
    end

    # if `f` is an extension to an external function, and we defined a method that overrides a method, for example,
    # we define `Base.isodd(n::Integer) = rand(Bool)`, which overrides the existing method `Base.isodd(n::Integer)`
    # calling `Base.delete_method` on this method won't bring back the old method, because our new method still exists in the method table, and it has a world age which is newer than the original. (our method has a deleted_world value set, which disables it)
    #
    # To solve this, we iterate again, and _re-enable any methods that were hidden in this way_, by adding them again to the method table with an even newer`primary_world`.
    if !isempty(deleted_sigs)
        to_insert = Method[]
        Base.visit(methods_table) do method
            if !isfromcell(method, cell_id) && method.sig ∈ deleted_sigs
                push!(to_insert, method)
            end
        end
        # separate loop to avoid visiting the recently added method
        for method in Iterators.reverse(to_insert)
            setfield!(method, primary_world, one(typeof(alive_world_val))) # `1` will tell Julia to increment the world counter and set it as this function's world
            ccall(:jl_method_table_insert, Cvoid, (Any, Any, Ptr{Cvoid}), methods_table, method, C_NULL) # i dont like doing this either!
        end
    end
    return !isempty(methods(f).ms)
end

# function try_delete_toplevel_methods(workspace::Module, name::Symbol)
#     try_delete_toplevel_methods(workspace, [name])
# end

function try_delete_toplevel_methods(workspace::Module, (cell_id, name_parts)::Tuple{UUID,Vector{Symbol}})::Bool
    try
        val = workspace
        for name in name_parts
            val = getfield(val, name)
        end
        try
            (val isa Function) && delete_toplevel_methods(val, cell_id)
        catch ex
            @warn "Failed to delete methods for $(name_parts)"
            showerror(stderr, ex, stacktrace(catch_backtrace()))
            false
        end
    catch
        false
    end
end

# these deal with some inconsistencies in Julia's internal (undocumented!) variable names
const primary_world = filter(in(fieldnames(Method)), [:primary_world, :min_world]) |> first # Julia v1.3 and v1.0 resp.
const deleted_world = filter(in(fieldnames(Method)), [:deleted_world, :max_world]) |> first # Julia v1.3 and v1.0 resp.
const alive_world_val = getfield(methods(Base.sqrt).ms[1], deleted_world) # typemax(UInt) in Julia v1.3, Int(-1) in Julia 1.0

###
# FORMATTING
###

"Because even showerror can error... 👀"
function try_showerror(io::IO, e, args...)
    try
        showerror(io, e, args...)
    catch show_ex
        print(io, "\nFailed to show error:\n\n")
        try_showerror(io, show_ex, stacktrace(catch_backtrace()))
    end
end

# We add a method for the Markdown -> HTML conversion that takes a LaTeX chunk from the Markdown tree and adds our custom span
function htmlinline(io::IO, x::LaTeX)
    withtag(io, :span, :class => "tex") do
        print(io, '$')
        htmlesc(io, x.formula)
        print(io, '$')
    end
end

# this one for block equations: (double $$)
function html(io::IO, x::LaTeX)
    withtag(io, :p, :class => "tex") do
        print(io, '$', '$')
        htmlesc(io, x.formula)
        print(io, '$', '$')
    end
end

"The `IOContext` used for converting arbitrary objects to pretty strings."
iocontext = IOContext(stdout, :color => false, :limit => true, :displaysize => (18, 88))
iocontext_compact = IOContext(iocontext, :compact => true)

const imagemimes = [MIME"image/svg+xml"(), MIME"image/png"(), MIME"image/jpg"(), MIME"image/jpeg"(), MIME"image/bmp"(), MIME"image/gif"()]
# in order of coolness
# text/plain always matches
"""
The MIMEs that Pluto supports, in order of how much I like them.

`text/plain` should always match - the difference between `show(::IO, ::MIME"text/plain", x)` and `show(::IO, x)` is an unsolved mystery.
"""
const allmimes = [MIME"application/vnd.pluto.tree+xml"(); MIME"text/html"(); imagemimes; MIME"text/latex"(); MIME"text/plain"()]

"""
Format `val` using the richest possible output, return formatted string and used MIME type.

See [`allmimes`](@ref) for the ordered list of supported MIME types.
"""
function format_output(@nospecialize(val))::MimedOutput
    try
        result, mime = sprint_withreturned(show_richest, val; context=iocontext)
        if mime ∈ imagemimes
            result, mime
        else
            String(result), mime
        end
    catch ex
        title = ErrorException("Failed to show value: \n" * sprint(try_showerror, ex))
        bt = stacktrace(catch_backtrace())
        format_output(CapturedException(title, bt))
    end
end

format_output(val::Nothing)::MimedOutput = "", MIME"text/plain"()

function format_output(val::CapturedException)::MimedOutput
    ## We hide the part of the stacktrace that belongs to Pluto's evalling of user code.
    stack = [s for (s,_) in val.processed_bt]

    for _ in 1:2
        until = findfirst(b -> b.func == :eval, reverse(stack))
        stack = until === nothing ? stack : stack[1:(length(stack) - until)]
    end

    pretty = map(stack[1:end]) do s
        Dict(
            :call => pretty_stackcall(s, s.linfo),
            :inlined => s.inlined,
            :file => basename(String(s.file)),
            :line => s.line,
        )
    end
    sprint(json, Dict(:msg => sprint(try_showerror, val.ex), :stacktrace => pretty)), MIME"application/vnd.pluto.stacktrace+json"()
end

# from the Julia source code:
function pretty_stackcall(frame::Base.StackFrame, linfo::Nothing)
    if frame.func isa Symbol
        String(frame.func)
    else
        repr(frame.func)
    end
end

function pretty_stackcall(frame::Base.StackFrame, linfo::Core.CodeInfo)
    "top-level scope"
end

function pretty_stackcall(frame::Base.StackFrame, linfo::Core.MethodInstance)
    if linfo.def isa Method
        sprint(Base.show_tuple_as_call, linfo.def.name, linfo.specTypes)
    else
        sprint(Base.show, linfo)
    end
end

"Like `Base.sprint`, but return a `(String, Any)` tuple containing function output as the second entry."
function sprint_withreturned(f::Function, args...; context=nothing, sizehint::Integer=0)
    s = IOBuffer(sizehint=sizehint)
    val = if context !== nothing
        f(IOContext(s, context), args...)
    else
        f(s, args...)
    end
    resize!(s.data, s.size), val
end

"Super important thing don't change."
struct 🥔 end
const struct_showmethod = which(show, (IO, 🥔))
const struct_showmethod_mime = which(show, (IO, MIME"text/plain", 🥔))

"""
Like two-argument `Base.show`, except:
1. the richest MIME type available to Pluto will be used
2. the used MIME type is returned
3. 'raw' data (e.g. image data) is always base64 encoded, with base64 header. This will change when/if we switch to a binary message format

With `onlyhtml=true`, the returned MIME type will always be MIME"text/html", and other MIME types are converted to this type. For example, an image with MIME type MIME"image/png" defined will display as:
```
<img src="data:image/png;base64,ahsdf87hf278hwh7823hr..." >
```
instead of (`onlyhtml=false`)
```
data:image/png;base64,ahsdf87hf278hwh7823hr...
```
"""
function show_richest(io::IO, @nospecialize(x); onlyhtml::Bool=false)::MIME
    mime = Iterators.filter(m -> Base.invokelatest(showable, m, x), allmimes) |> first
    t = typeof(x)

    # types that have no specialized show methods (their fallback is text/plain) are displayed using Pluto's interactive tree viewer.
    # this is how we check whether this display method is appropriate:
    isstruct = try
        mime isa MIME"text/plain" &&
        t isa DataType &&
        # there are two ways to override the plaintext show method:
        which(show, (IO, MIME"text/plain", t)) === struct_showmethod_mime &&
        which(show, (IO, t)) === struct_showmethod
    catch
        false
    end

    if isstruct
        show_struct(io, x)
        return MIME"application/vnd.pluto.tree+xml"()
    end

    if mime ∈ imagemimes
        if onlyhtml
            # if only html output is accepted, we need to base64 encode the result and use it as image source.
            enc_pipe = Base64.Base64EncodePipe(io)
            io_64 = IOContext(enc_pipe, iocontext)

            print(io, "<img src=\"data:", mime, ";base64,")
            show(io_64, mime, x)
            close(enc_pipe)
            print(io, "\">")
            return MIME"text/html"()
        else
            show(io, mime, x)
            return mime
        end
    else
        if onlyhtml || mime isa MIME"text/latex"
            # see onlyhtml description in docstring
            if mime isa MIME"text/plain"
                withtag(io, :pre) do
                    htmlesc(io, repr(mime, x; context=iocontext_compact))
                end
            elseif mime isa MIME"text/latex"
                # Wrapping with `\text{}` allows for LaTeXStrings with mixed text/math
                texed = repr(mime, x)
                html(io, Markdown.LaTeX("\\text{$texed}"))
            else
                show(io, mime, x)
            end
            return MIME"text/html"()
        else
            # the classic:
            show(io, mime, x)
            return mime
        end
    end
end

###
# TREE VIEWER
###

# We invent our own MIME _because we can_ but don't use it somewhere else because it might change :)

const tree_display_limit = 50

function show_array_row(io::IO, pair::Tuple)
    i, element = pair
    print(io, "<r><k>", i, "</k><v>")
    show_richest(io, element; onlyhtml=true)
    print(io, "</v></r>")
end

function show_array_elements(io::IO, indices::AbstractVector{<:Integer}, x::AbstractArray{<:Any, 1})
    for i in indices
        if isassigned(x, i)
            show_array_row(io, (i, x[i]))
        else
            show_array_row(io, (i, Text(Base.undef_ref_str)))
        end
    end
end

function show_dict_row(io::IO, pair::Union{Pair,Tuple})
    k, element = pair
    print(io, "<r><k>")
    if pair isa Pair
        show_richest(io, k; onlyhtml=true)
    else
        # this is an entry of a NamedTuple, the first element of the Tuple is a Symbol, which we want to print as `x` instead of `:x`
        print(io, k)
    end
    print(io, "</k><v>")
    show_richest(io, element; onlyhtml=true)
    print(io, "</v></r>")
end

istextmime(::MIME"application/vnd.pluto.tree+xml") = true

function array_prefix(io, x::Array{<:Any, 1})
    print(io, eltype(x))
end
function array_prefix(io, x)
    original = sprint(Base.showarg, x, false)
    print(io, lstrip(original, ':'))
    print(io, ": ")
end

Base.showable(::MIME"application/vnd.pluto.tree+xml", x::AbstractRange) = false

function show(io::IO, ::MIME"application/vnd.pluto.tree+xml", x::AbstractArray{<:Any, 1})
    print(io, """<jltree class="collapsed" onclick="onjltreeclick(this, event)">""")
    array_prefix(io, x)
    print(io, "<jlarray>")
    indices = eachindex(x)

    if length(x) <= tree_display_limit
        show_array_elements(io, indices, x)
    else
        firsti = firstindex(x)
        from_end = tree_display_limit > 20 ? 10 : 1

        show_array_elements(io, indices[firsti:firsti-1+tree_display_limit-from_end], x)

        print(io, "<r><more></more></r>")

        show_array_elements(io, indices[end+1-from_end:end], x)
    end

    print(io, "</jlarray>")
    print(io, "</jltree>")
end

function show(io::IO, ::MIME"application/vnd.pluto.tree+xml", x::Tuple)
    print(io, """<jltree class="collapsed" onclick="onjltreeclick(this, event)">""")
    print(io, """<jlarray class="Tuple">""")
    show_array_row.([io], zip(eachindex(x), x))
    print(io, "</jlarray>")
    print(io, "</jltree>")
end

function show(io::IO, ::MIME"application/vnd.pluto.tree+xml", x::AbstractDict{<:Any, <:Any})
    print(io, """<jltree class="collapsed" onclick="onjltreeclick(this, event)">""")
    print(io, typeof(x) |> trynameof)
    print(io, "<jldict>")
    row_index = 1
    for pair in x
        show_dict_row(io, pair)
        if row_index == tree_display_limit
            print(io, "<r><more></more></r>")
            break
        end
        row_index += 1
    end

    print(io, "</jldict>")
    print(io, "</jltree>")
end

function show(io::IO, ::MIME"application/vnd.pluto.tree+xml", x::NamedTuple)
    print(io, """<jltree class="collapsed" onclick="onjltreeclick(this, event)">""")
    print(io, """<jldict class="NamedTuple">""")
    show_dict_row.([io], zip(eachindex(x), x))
    print(io, "</jldict>")
    print(io, "</jltree>")
end

function show(io::IO, ::MIME"application/vnd.pluto.tree+xml", x::Pair)
    print(io, """<jlpair>""")
    show_dict_row(io, x)
    print(io, "</jlpair>")
end

# Based on Julia source code, but HTML-ified
function show_struct(io::IO, @nospecialize(x))
    t = typeof(x)
    nf = nfields(x)
    nb = sizeof(x)
    if nf != 0 || nb == 0
        print(io, """<jltree class="collapsed" onclick="onjltreeclick(this, event)">""")
        show(io, t)
        print(io, "<jlstruct>")

        if !Base.show_circular(io, x)
            recur_io = IOContext(io, Pair{Symbol,Any}(:SHOWN_SET, x),
                                 Pair{Symbol,Any}(:typeinfo, Any))
            for i in 1:nf
                f = fieldname(t, i)
                if !isdefined(x, f)
                    print(io, "<r>", Base.undef_ref_str, "</r>")
                else
                    show_array_row(recur_io, (f, getfield(x, i)))
                end
            end
        end

        print(io, "</jlstruct>")
        print(io, "</jltree>")
    else
        Base.show_default(io, x)
    end
end

trynameof(x::DataType) = nameof(x)
trynameof(x::Any) = Symbol()

###
# JSON SERIALIZER
###

# We define a minimal JSON serializer here so that the notebook process does not need to depend on JSON.jl
# Performance is about 0.5-1.0x JSON.jl, but that's okay since it is only used for special output types like stack traces
# Not designed/tested for use outside of Pluto

struct ReplacePipe <: IO
	outstream::IO
end

# to get these character codes:
# [c => UInt8(c) for c in "\"\\/\b\f\n\r\t"]
# we can do this escaping per-byte because UTF-8 is backwards compatible with ASCII, i.e. these special characters are never part of a UTF-8 encoded character other than the ASCII characters they represent. Cool!
function Base.write(rp::ReplacePipe, x::UInt8)
	if x == 0x22 || x== 0x5c || x== 0x2f # https://www.json.org/json-en.html
        write(rp.outstream, '\\')
        write(rp.outstream, x)
    elseif x < 0x10 # ish
        write(rp.outstream, escape_string(String([Char(x)]))) # the Julia escaping 'happens' to coincide with what we want
    else
        write(rp.outstream, x)
    end
end
function sanitize_pipe(func::Function, outstream::IO, args...)
	func(ReplacePipe(outstream), args...)
end

function json(io, arr::AbstractArray)
    write(io, '[')
    len = length(arr)
    for (i, val) in enumerate(arr)
        json(io, val)
        (i != len) && write(io, ',')
    end
    write(io, ']')
end

function json(io, d::Dict{Symbol, T}) where T
    write(io, '{')
    len = length(d)
    for (i, val) in enumerate(d)
        write(io, '"', val.first, '"', ':')
        json(io, val.second)
        (i != len) && write(io, ',')
    end
    write(io, '}')
end

function json(io, str::T) where T<:AbstractString
    write(io, '"')
    sanitize_pipe(write, io, str)
    write(io, '"')
end

function json(io, val::Any)
    show(io, val)
end

###
# REPL THINGS
###

function completion_priority((s, description, exported))
	c = first(s)
	if islowercase(c)
		1 - 10exported
	elseif isuppercase(c)
		2 - 10exported
	else
		3 - 10exported
	end
end

completed_object_description(x::Function) = "Function"
completed_object_description(x::Number) = "Number"
completed_object_description(x::AbstractString) = "String"
completed_object_description(x::Module) = "Module"
completed_object_description(x::AbstractArray) = "Array"
completed_object_description(x::Any) = "Any"

completion_description(c::ModuleCompletion) = try
    completed_object_description(getfield(c.parent, Symbol(c.mod)))
catch
    nothing
end
completion_description(::Completion) = nothing

function completions_exported(cs::Vector{<:Completion})
    completed_modules = Set(c.parent for c in cs if c isa ModuleCompletion)
    completed_modules_exports = Dict(m => string.(names(m, all=false, imported=true)) for m in completed_modules)

    map(cs) do c
        if c isa ModuleCompletion
            c.mod ∈ completed_modules_exports[c.parent]
        else

            true
        end
    end
end

"You say Linear, I say Algebra!"
function completion_fetcher(query, pos, workspace::Module=current_module)
    results, loc, found = completions(query, pos, workspace)

    texts = completion_text.(results)
    descriptions = completion_description.(results)
    exported = completions_exported(results)

    smooshed_together = zip(texts, descriptions, exported)

    final = sort(collect(smooshed_together); alg=MergeSort, by=completion_priority)
    (final, loc, found)
end

# Based on /base/docs/bindings.jl from Julia source code
function binding_from(x::Expr, workspace::Module=current_module)
    if x.head == :macrocall
        Docs.Binding(workspace, x.args[1])
    elseif x.head == :.
        Docs.Binding(Core.eval(workspace, x.args[1]), x.args[2].value)
    else
        error("Invalid @var syntax `$x`.")
    end
end
binding_from(s::Symbol, workspace::Module=current_module) = Docs.Binding(workspace, s)
binding_from(r::GlobalRef, workspace::Module=current_module) = Docs.Binding(r.mod, r.name)
binding_from(other, workspace::Module=current_module) = error("Invalid @var syntax `$other`.")

"You say doc_fetch, I say You say doc_fetch, I say You say doc_fetch, I say You say doc_fetch, I say ...!!!!"
function doc_fetcher(query, workspace::Module=current_module)
    try
        binding = binding_from(Meta.parse(query), workspace)::Docs.Binding
        (repr(MIME"text/html"(), Docs.doc(binding)), :👍)
    catch ex
        (nothing, :👎)
    end
end

###
# BONDS
###

"""
_“The name is Bond, James Bond.”_

Wraps around an `element` and not much else. When you `show` a `Bond` with the `text/html` MIME type, you will get:

```html
<bond def="\$(bond.defines)">
\$(repr(MIME"text/html"(), bond.element))
</bond>
```

For example, `Bond(html"<input type=range>", :x)` becomes:

```html
<bond def="x">
<input type=range>
</bond>
```

The actual reactive-interactive functionality is not done in Julia - it is handled by the Pluto front-end (JavaScript), which searches cell output for `<bond>` elements, and attaches event listeners to them. Put on your slippers and have a look at the JS code to learn more.
"""
struct Bond
    element::Any
    defines::Symbol
    Bond(element, defines::Symbol) = showable(MIME"text/html"(), element) ? new(element, defines) : error("""Can only bind to html-showable objects, ie types T for which show(io, ::MIME"text/html", x::T) is defined.""")
end

import Base: show
function show(io::IO, ::MIME"text/html", bond::Bond)
    withtag(io, :bond, :def => bond.defines) do
        show(io, MIME"text/html"(), bond.element)
    end
end

"""
    `@bind symbol element`

Return the HTML `element`, and use its latest JavaScript value as the definition of `symbol`.

# Example

```julia
@bind x html"<input type=range>"
```
and in another cell:
```julia
x^2
```

The first cell will show a slider as the cell's output, ranging from 0 until 100.
The second cell will show the square of `x`, and is updated in real-time as the slider is moved.
"""
macro bind(def, element)
	if def isa Symbol
		quote
			local el = $(esc(element))
            global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
			PlutoRunner.Bond(el, $(Meta.quot(def)))
		end
	else
		:(throw(ArgumentError("""\nMacro example usage: \n\n\t@bind my_number html"<input type='range'>"\n\n""")))
	end
end

"""
Will be inserted in saved notebooks that use the @bind macro, make sure that they still contain legal syntax when executed as a vanilla Julia script. Overloading `Base.get` for custom UI objects gives bound variables a sensible value.
"""
const fake_bind = """macro bind(def, element)
    quote
        local el = \$(esc(element))
        global \$(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end"""

###
# LOGGING
###

const log_channel = Channel{Any}(10)
const old_logger = Ref{Any}(nothing)

struct PlutoLogger <: Logging.AbstractLogger
    stream
end

function Logging.shouldlog(::PlutoLogger, level, _module, _...)
    # Accept logs
    # - From the user's workspace module
    # - Info level and above for other modules
    _module === current_module || convert(Logging.LogLevel, level) >= Logging.Info
end
Logging.min_enabled_level(::PlutoLogger) = Logging.Debug
Logging.catch_exceptions(::PlutoLogger) = false
function Logging.handle_message(::PlutoLogger, level, msg, _module, group, id, file, line; kwargs...)
    try
        put!(log_channel, (
            level=string(level),
            msg=(msg isa String) ? msg : repr(msg),
            group=group,
            # id=id,
            file=file,
            line=line,
            kwargs=Dict((k=>repr(v) for (k,v) in kwargs)...),
        ))
        # also print to console
        Logging.handle_message(old_logger[], level, msg, _module, group, id, file, line; kwargs...)
    catch e
        println(stderr, "Failed to relay log from PlutoRunner")
        showerror(stderr, e, stacktrace(catch_backtrace()))
    end
end

# we put this in __init__ to fix a world age problem
function __init__()
    if Distributed.myid() != 1
        old_logger[] = Logging.global_logger()
        Logging.global_logger(PlutoLogger(nothing))
    end
end

end

ext: jl


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Model of basic streams
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt and Rafael de P. Soares
#* $Id$
#*---------------------------------------------------------------------
type sink
	sink()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of chemical components",
				:Lower=>1
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Component Mol Weight",
				:Protected=>true
			)),(NComp)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet Stream",
				:PosX=>0,
				:PosY=>0.5308,
				:Protected=>true,
				:Symbol=>"_{in}"
			)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Vapourization fraction",
				:Hidden=>true
			)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Fraction",
				:Hidden=>true
			)),(NComp)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Vapour Molar Fraction",
				:Hidden=>true
			)),(NComp)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Mass Fraction"
			)),(NComp)),
			molweight(Dict{Symbol,Any}(
				:Brief=>"Average Mol Weight"
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Molar Volume"
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Stream Mass Density"
			)),
			dens_mol(Dict{Symbol,Any}(
				:Brief=>"Stream Molar Density"
			)),
			flow_mass(Dict{Symbol,Any}(
				:Brief=>"Stream Mass Flow"
			)),
			flow_vol(Dict{Symbol,Any}(
				:Brief=>"Volumetric Flow"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Temperature in  C",
				:Lower=>-250,
				:Upper=>5000
			)),
			[
				:([v, x, y] = PP.FlashPH(Inlet.P, Inlet.h, Inlet.z)),
				:(Mw = sum(M*Inlet.z)),
				:(rhom * vm = 1),
				:(rhom * Mw = rho),
				:(Fw = Mw*Inlet.F),
				:(vm = (1-v)*PP.LiquidVolume(Inlet.T, Inlet.P, x) + v*PP.VapourVolume(Inlet.T,Inlet.P,y)),
				:(Fvol = Inlet.F*vm),
				:(zmass = M*Inlet.z / Mw),
				:(T_Cdeg = Inlet.T/"K" - 273.15),
			],
			[
				"Flash Calculation","Average Molecular Weight","Molar Density","Mass or Molar Density","Flow Mass","Molar Volume","Volumetric Flow","Mass Fraction","Temperature in  C",
			],
			[:PP,:NComp,:M,],
			[:Inlet,:v,:x,:y,:zmass,:Mw,:vm,:rho,:rhom,:Fw,:Fvol,:T_Cdeg,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	M::Array{molweight}
	Inlet::stream
	v::fraction
	x::Array{fraction}
	y::Array{fraction}
	zmass::Array{fraction}
	Mw::molweight
	vm::volume_mol
	rho::dens_mass
	rhom::dens_mol
	Fw::flow_mass
	Fvol::flow_vol
	T_Cdeg::DanaReal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export sink
function set(in::sink)
	M = PP.MolecularWeight()

end
function setEquationFlow(in::sink)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
end
function atributes(in::sink,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Sink"
	fields[:Brief]="Material stream sink"
	fields[:Info]="
	This model should be used for boundary streams when additional
	information about the stream is desired.

	Some of the additional informations calculated by this models are:
	 * Mass density
	 * Mass flow
	 * Mass compostions
	 * Specific volume
	 * Vapour fraction
	 * Volumetric flow
	 * Liquid and Vapour compositions
	"
	drive!(fields,_)
	return fields
end
sink(_::Dict{Symbol,Any})=begin
	newModel=sink()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
import Splittings: UniformMesh

@testset "UniformMesh" begin

mesh = UniformMesh( 0, 1, 11; endpoint=true)

@test mesh.start   == 0.0
@test mesh.step    == 0.1
@test mesh.stop    == 1.0
@test mesh.length  == 11

mesh = UniformMesh( 0, 1, 10; endpoint=false)

@test mesh.start   == 0.0
@test mesh.step    == 0.1
@test mesh.stop    == 1.0
@test mesh.length  == 10

end


text:
let
    #= line 1 =#
    var"|main|b_main"::Any = 1::Any
    var"|main|c_x"::Any = 1::Any
    var"|main|d_unit"::Type = Type::Type
    var"|main|f_println"::Function = println::Function
    var"|main|g_eq"::Function = (x->begin
                    #= none:1 =#
                    #= none:3 =#
                    y->begin
                            #= none:3 =#
                            #= none:5 =#
                            x == y
                        end
                end)::Function
    (var"|main|f_println"::Function)(if ((var"|main|g_eq"::Function)(1::Any)::Function)(2::Any)::Bool
                    1::Any
                else
                    2::Any
                end::Any)::Any
    var"|main|h_id"::Function = function (var"|main|i_x",)
                #= main.mlfs:23 =#
                var"|main|i_x"::Any
            end::Function
    var"|main|j_choose"::Function = function (var"|main|k_x",)
                #= main.mlfs:26 =#
                function (var"|main|l_y",)
                        #= main.mlfs:26 =#
                        var"|main|k_x"::Any
                    end::Function
            end::Function
    var"|main|m_choose_id"::Function = (var"|main|j_choose"::Function)(var"|main|h_id"::Function)::Function
    var"|main|n_x"::Function = (var"|main|m_choose_id"::Function)(var"|main|h_id"::Function)::Function
    var"|main|o_Nat"::Type = Type::Type
    var"|main|q_S"::Function = (x::UInt->begin
                    #= none:1 =#
                    #= none:3 =#
                    1 + x
                end)::Function
    var"|main|r_Z"::Any = UInt(0)::Any
    var"|main|s_matchNat"::Function = (f1->begin
                    #= none:1 =#
                    #= none:3 =#
                    f2->begin
                            #= none:3 =#
                            #= none:5 =#
                            x->begin
                                    #= none:5 =#
                                    #= none:7 =#
                                    if x == 0
                                        #= none:8 =#
                                        f2(x)
                                    else
                                        #= none:10 =#
                                        f1(x)
                                    end
                                end
                        end
                end)::Function
    (var"|main|f_println"::Function)((var"|main|n_x"::Function)("hellow world!"::String)::String)::Any
    (var"|main|f_println"::Function)((((var"|main|s_matchNat"::Function)(function (var"|main|t_x",)
                                    #= main.mlfs:50 =#
                                    "isS"::String
                                end::Function)::Function)(function (var"|main|u_x",)
                                #= main.mlfs:50 =#
                                "isZ"::String
                            end::Function)::Function)((var"|main|q_S"::Function)(var"|main|r_Z"::Any)::Any)::String)::Any
    var"|main|"::Tuple = ()
    main::Tuple = (("eq", "i64", "main", "i32", "str", "S", "i16", "choose", "bool", "matchNat", "x", "Z", "unit", "field", "module", "char", "println", "Nat", "fieldnames", "choose_id", "i8", "id", "Type"), (var"|main|g_eq", Type, var"|main|b_main", Type, Type, var"|main|q_S", Type, var"|main|j_choose", Type, var"|main|s_matchNat", var"|main|n_x", var"|main|r_Z", Type, Type, Type, Type, var"|main|f_println", Type, Type, var"|main|m_choose_id", Type, var"|main|h_id", Type))
    0
end


text:
using Documenter, NLSolvers

makedocs(
    doctest = false,
    sitename = "NLSolvers.jl"
)

deploydocs(
    deps = Deps.pip("pygments", "mkdocs", "python-markdown-math", "mkdocs-windmill"),
    repo = "github.com/JuliaNLSolvers/NLSolvers.jl.git",
)


text:
@testset "load_file_line_by_line" begin
    file_txt = ["1","2","","a","b c","d e f",""]

    @test load_file_line_by_line("data/file.txt") == file_txt
    @test_throws ErrorException("File a does not exist.") load_file_line_by_line("a")
end


text:
function _precompile_tullio()
    precompile(extern_prod, (Vector{Float64}, Vector{Float64},))
    precompile(pairwise_mul!, (Matrix{Float64}, Matrix{Float64},))
    precompile(pairwise_diff, (Vector{Float64}, Vector{Float64},))
    precompile(prod_reduction, (Matrix{Float64}, Vector{Float64}, Vector{Float64},))
    precompile(prod_reduction, (Matrix{Float64}, Vector{Float64},))
end

max_stars_repo_name: edwinb-ai/LeastSquaresSVM

max_issues_repo_name: edwinb-ai/Elysivm

max_forks_repo_name: edwinb-ai/LeastSquaresSVM


text:
get_aggregation_topology_accessor(::Type{Area}) = get_area
get_aggregation_topology_accessor(::Type{LoadZone}) = get_load_zone

set_load_zone!(bus::Bus, load_zone::LoadZone) = bus.load_zone = load_zone
set_area!(bus::Bus, area::Area) = bus.area = area

"""
    removes the aggregation topology in a Bus
"""
_remove_aggregration_topology!(bus::Bus, ::LoadZone) = bus.load_zone = nothing
_remove_aggregration_topology!(bus::Bus, ::Area) = bus.area = nothing

function set_dynamic_injector!(
    static_injector::T,
    dynamic_injector::U,
) where {
    T <: Union{
        GenericBattery,
        HydroDispatch,
        HydroEnergyReservoir,
        InterruptibleLoad,
        PowerLoad,
        RenewableDispatch,
        RenewableFix,
        ThermalStandard,
    },
    U <: Union{Nothing, DynamicInjection},
}
    current_dynamic_injector = get_dynamic_injector(static_injector)
    if !isnothing(current_dynamic_injector) && !isnothing(dynamic_injector)
        throw(ArgumentError("cannot assign a dynamic injector on a device that already has one"))
    end

    # All of these types implement this field.
    static_injector.dynamic_injector = dynamic_injector
end

function set_static_injector!(
    dynamic_injector::T,
    static_injector::U,
) where {
    T <: Union{DynamicGenerator, DynamicInverter},
    U <: Union{Nothing, StaticInjection},
}
    current_static_injector = get_static_injector(dynamic_injector)
    if !isnothing(current_static_injector) && !isnothing(static_injector)
        throw(ArgumentError("cannot assign a static injector on a device that already has one"))
    end

    # All of these types implement this field.
    dynamic_injector.static_injector = static_injector
end


text:
mutable struct IterateChannel{A <: AbstractArray} <: AbstractChannel{A}
    decisions::Dict{Int,A}
    cond_take::Condition
    IterateChannel(decisions::Dict{Int,A}) where A <: AbstractArray = new{A}(decisions, Condition())
end

function put!(channel::IterateChannel, t, x)
    channel.decisions[t] = copy(x)
    notify(channel.cond_take)
    return channel
end

function take!(channel::IterateChannel, t)
    x = fetch(channel, t)
    delete!(channel.decisions, t)
    return x
end

isready(channel::IterateChannel) = length(channel.decisions) > 1
isready(channel::IterateChannel, t) = haskey(channel.decisions, t)

function fetch(channel::IterateChannel, t)
    wait(channel, t)
    return channel.decisions[t]
end

function wait(channel::IterateChannel, t)
    while !isready(channel, t)
        wait(channel.cond_take)
    end
end

RemoteIterates{A} = RemoteChannel{IterateChannel{A}}

mutable struct MetaChannel <: AbstractChannel{Any}
    metadata::Dict{Tuple{Int,Symbol},Any}
    cond_take::Condition
    MetaChannel() = new(Dict{Tuple{Int,Symbol},Any}(), Condition())
end

function put!(channel::MetaChannel, t, key, x)
    channel.metadata[(t,key)] = copy(x)
    notify(channel.cond_take)
    return channel
end

function take!(channel::MetaChannel, t, key)
    x = fetch(channel, t, key)
    delete!(channel.metadata, (t,k))
    return x
end

isready(channel::MetaChannel) = length(channel.metadata) > 1
isready(channel::MetaChannel, t, key) = haskey(channel.metadata, (t,key))

function fetch(channel::MetaChannel, t, key)
    wait(channel, t, key)
    return channel.metadata[(t,key)]
end

function wait(channel::MetaChannel, t, key)
    while !isready(channel, t, key)
        wait(channel.cond_take)
    end
end

MetaData = RemoteChannel{MetaChannel}


text:
# Trace for Profiling
# ===================

mutable struct Trace
    starttime::Nullable{DateTime}
    snapshots::Vector{Tuple{UInt64,Stats}}

    function Trace()
        return new(nothing, [])
    end
end

function takesnapshot!(trace, stats)
    t = time_ns()
    if isempty(trace.snapshots)
        trace.starttime = now()
    end
    push!(trace.snapshots, (t, deepcopy(stats)))
    return trace
end

function elapsedtime(trace)
    if isempty(trace.snapshots)
        return typemax(UInt64)
    else
        return time_ns() - trace.snapshots[end][1]
    end
end

function Base.print(io::IO, trace::Trace)
    join(io,
         (:snapshot,
          :datetime,
          :calltime_ns,
          :selecttime_ns,
          :remotecalltime_ns,
          :waitargstime_ns), '\t')
    if isempty(trace.snapshots)
        return
    end
    starttime = get(trace.starttime)
    starttime_ns, _ = trace.snapshots[1]
    for i in 1:endof(trace.snapshots)
        println(io)
        time_ns, stats = trace.snapshots[i]
        datetime = starttime + Dates.Nanosecond(time_ns - starttime_ns)
        join(io,
             (i,
              datetime,
              stats.calltime_ns,
              stats.selecttime_ns,
              stats.remotecalltime_ns,
              stats.waitargstime_ns), '\t')
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SeqGaussSim(var₁=>param₁, var₂=>param₂, ...)

A sequential Gaussian simulation solver.

## Parameters

* `variogram` - Variogram model (default to `GaussianVariogram()`)
* `mean`      - Simple Kriging mean
* `degree`    - Universal Kriging degree
* `drifts`    - External Drift Kriging drift functions

Latter options override former options. For example, by specifying
`drifts`, the user is telling the algorithm to ignore `degree` and
`mean`. If no option is specified, Ordinary Kriging is used by
default with the `variogram` only.

* `neighborhood` - Neighborhood on which to search neighbors
* `maxneighbors` - Maximum number of neighbors (default to 10)
* `path`         - Simulation path (default to `RandomPath`)

For each location in the simulation `path`, a maximum number of
neighbors `maxneighbors` is used to fit a Gaussian distribution.
The neighbors are searched according to a `neighborhood`.
"""
@simsolver SeqGaussSim begin
  @param variogram = GaussianVariogram()
  @param mean = nothing
  @param degree = nothing
  @param drifts = nothing
  @param neighborhood
  @param minneighbors = 1
  @param maxneighbors = 10
  @param path = nothing
end

function preprocess(problem::SimulationProblem, solver::SeqGaussSim)
  # retrieve problem info
  pdomain = domain(problem)

  params = []
  for (var, V) in variables(problem)
    # get user parameters
    if var ∈ keys(solver.params)
      varparams = solver.params[var]
    else
      varparams = SeqGaussSimParam()
    end

    # determine which Kriging variant to use
    if varparams.drifts ≠ nothing
      estimator = ExternalDriftKriging(varparams.variogram, varparams.drifts)
    elseif varparams.degree ≠ nothing
      estimator = UniversalKriging(varparams.variogram, varparams.degree, ndims(pdomain))
    elseif varparams.mean ≠ nothing
      estimator = SimpleKriging(varparams.variogram, varparams.mean)
    else
      estimator = OrdinaryKriging(varparams.variogram)
    end

    # determine marginal distribution
    marginal = Normal()

    # determine simulation path
    path = varparams.path ≠ nothing ? varparams.path : LinearPath(pdomain)

    # equivalent parameters for SeqSim solver
    param = (estimator=estimator,
             neighborhood=varparams.neighborhood,
             minneighbors=varparams.minneighbors,
             maxneighbors=varparams.maxneighbors,
             marginal=marginal, path=path)

    push!(params, var => param)
  end

  preprocess(problem, SeqSim(params...))
end

solve_single(problem::SimulationProblem, var::Symbol, solver::SeqGaussSim, preproc) =
  solve_single(problem, var, SeqSim(), preproc)


text:
#=
Functions for temperature dependence of biological rates :

- NoEffectTemperature : temperature is not included
- ExtendedEppley : Extended (hump-shaped) Eppley function
- ExponentialBA : Exponential Boltzmann-Arrhénius function
- ExtendedBA : Extended (hump-shaped) Boltzmann-Arrhénius function
- Gaussian : Gaussian, or inverse Gaussian, function

These functions take as argument :
- rate_affected : the type of biological rate (metabolic, growth, attack rate or handling time) specified with a key (symbol)
- parameters_tuple : optionnal, a tuple with the parameters corresponding to the specified temperature dependence function

If no parameters are specified, default parameters are provided. However, as the parameters for temperature dependence functions can be quite specific,
we encourage to specify parameters corresponding to your study system.
All parameters can be changed or just a few of them.

The functions are called in model parameters as, for instance :

growthrate = NoEffectTemperature(:growthrate) -> no effect of temperature for growth rate
metabolicrate = ExtendedEppley(:metabolicrate) -> Extended Eppley function for metabolic rate
attackrate = ExponentialBA(:attackrate, parameters_tuple = parameters) -> Exponental BA function for attack rate, parameters are specified in the tuple 'parameters'

According to the 'rate_affected' parameter, each function calls the specified temperature dependence function defined in temperature_dependence_functions.jl
For instance : when metabolicrate = ExtendedEppley(:metabolicrate), the function extended_eppley_x from temperature_dependence_functions.jl is used to define the metabolic rate

The functions return either a value, a vector or a matrix for the value of the biological rate, depending on which biological rate is chosen.
=#

"""
**No effect of temperature on biological rates**

    NoEffectTemperature(::Symbol; parameters_tuple...)

Calculates biological rates independently of temperature.

Called by default in `model_parameters` to define temperature independent biological rates. Default
parameters values are provided, but can be overwritten through the keyword argument `parameters_tuple`
(see example). See the complete documentation for the default values and more details.

Arguments:

- rate_affected : the biological rate that the function should calculate. May be :growth, :metabolism, :handlingtime, :attackrate
- parameters_tuple : a named tuple specifying parameters values. Note that if you provide a one-element tuple, it should end with a comma to avoid being treated as a vector.

# Example :
```julia-repl
# Metabolic rate with default parameters (this is the default behavior of model_parameters)
julia> p = model_parameters(A, metabolicrate = NoEffectTemperature(:x))
# Growth rate, changing the default parameters
julia> p = model_parameters(A, growthrate = NoEffectTemperature(:r, parameters_tuple = (r = 1.5,)))
```
"""
function NoEffectTemperature(rate_affected::Symbol; parameters_tuple...)
    if rate_affected ∉ [:growth, :r, :metabolism, :x, :handlingtime, :attackrate]
        error("rate_affected should be either :growth (alternatively :r or :growthrate), :metabolism (alternatively :x or :metabolicrate),
        :handlingtime or :attackrate")
    end
    if rate_affected ∈ [:growth, :r, :growthrate]
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            no_effect_r(passed_temp_parameters = pt)
        else
            no_effect_r()
        end
    elseif rate_affected ∈ [:metabolism, :x, :metabolicrate]
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            no_effect_x(passed_temp_parameters = pt)
        else
            no_effect_x()
        end
    elseif rate_affected == :handlingtime
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            no_effect_handlingt(passed_temp_parameters = pt)
        else
            no_effect_handlingt()
        end
    elseif rate_affected == :attackrate
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            no_effect_attackr(passed_temp_parameters = pt)
        else
            no_effect_attackr()
        end
    end
end

"""
**Extended Eppley**

    ExtendedEppley(::Symbol, parameters_tuple...)

Calculates biological rates using an Extended Eppley function.

Can be called in `model_parameters` (see example) to define temperature dependent biological rates. The Extended Eppley
function is only compatible with metabolic and growth rate and cannot be used to calculate handling time or attack rate.
Default parameters values are provided, but can be overwritten through the keyword argument `parameters_tuple`. See the
complete documentation for the default values, their sources and more details.

Arguments:
- rate_affected : the biological rate that the function should calculate. May be :growth, :metabolism
- parameters_tuple : a named tuple specifying parameters values. Note that if you provide a one-element tuple, it should end with a comma to avoid being treated as a vector.

# Example
```julia-repl
# Metabolic rate with default parameters values
p = model_parameters(A, metabolicrate = ExtendedEppley(:x))
# Metabolic rate, changing some of the default parameters values
p = model_parameters(A, metabolicrate = ExtendedEppley(:x,
                                                       parameters_tuple = (maxrate_0_producer = 0.7,
                                                                           maxrate_0_vertebrate = 0.9))

```
"""
function ExtendedEppley(rate_affected::Symbol; parameters_tuple...)
    if rate_affected ∉ [:growth, :growthrate, :r, :metabolism, :x, :metabolicrate] ; error("rate_affected should be either :growth (alternatively :r or :growthrate) or :metabolism (alternatively :x or :metabolicrate)") ; end
    if rate_affected ∈ [:growth, :r, :growthrate]
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            extended_eppley_r(passed_temp_parameters = pt)
        else
            extended_eppley_r()
        end
    elseif rate_affected ∈ [:metabolism, :x, :metabolicrate]
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            extended_eppley_x(passed_temp_parameters = pt)
        else
            extended_eppley_x()
        end
    end
end

"""
**Exponential Boltzmann-Arrhenius**

    ExponentialBA(::Symbol, parameters_tuple...)

Calculates biological rates using an Exponential Boltzmann Arrhenius function.

Can be called in `model_parameters` (see example) to define temperature dependent biological rates. Default
parameters values are provided, but can be overwritten through the keyword argument `parameters_tuple`. See
the complete documentation for the default values, their sources and more details.

Arguments:
- rate_affected : the biological rate that the function should calculate. May be :growth, :metabolism, :handlingtime, :attackrate
- parameters_tuple : a named tuple specifying parameters values. Note that if you provide a one-element tuple, it should end with a comma to avoid being treated as a vector.

# Example
```julia-repl
# Handling time with default parameters values
p = model_parameters(A, handlingtime = ExponentialBA(:handlingtime))
# Metabolic rate, changing some of the default parameters values
p = model_parameters(A, metabolicrate = ExponentialBA(:x, parameters_tuple = (T0 = 300.0,))
```
"""
function ExponentialBA(rate_affected::Symbol; parameters_tuple...)
    if rate_affected ∉ [:growth, :r, :growthrate, :metabolism, :x, :metabolicrate, :handlingtime, :attackrate]
        error("rate_affected should be either :growth (alternatively :r or :growthrate), :metabolism (alternatively :x or :metabolicrate), :attackrate or :handlingtime")
    end
    if rate_affected ∈ [:growth, :r, :growthrate]
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            exponential_BA_r(passed_temp_parameters = pt)
        else
            exponential_BA_r()
        end
    elseif rate_affected ∈ [:metabolism, :x, :metabolicrate]
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            exponential_BA_x(passed_temp_parameters = pt)
        else
            exponential_BA_x()
        end
    elseif rate_affected == :attackrate
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            exponential_BA_attackr(passed_temp_parameters = pt)
        else
            exponential_BA_attackr()
        end
    elseif rate_affected == :handlingtime
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            exponential_BA_handlingt(passed_temp_parameters = pt)
        else
            exponential_BA_handlingt()
        end
    end
end

"""
**Extended  Boltzmann-Arrhenius**

    ExtendedBA(::Symbol, parameters_tuple...)

Calculates biological rates using an Extended Boltzmann Arrhenius function.

Can be called in `model_parameters` (see example) to define temperature dependent biological rates. Only compatible
with metabolic, growth and attack rates. Default parameters values are provided, but can be overwritten through the
keyword argument `parameters_tuple`. See the complete documentation for the default values, their sources and more details.

Arguments:
- rate_affected : the biological rate that the function should calculate. May be :growth, :metabolism, :attackrate
- parameters_tuple : a named tuple specifying parameters values. Note that is you provide a one-element tuple, it should end with a comma to avoid being treated as a vector.

# Example
```julia-repl
# Attack rate with default parameters values
p = model_parameters(A, attackrate = ExtendedBA(:attackrate))
# Metabolic rate, changing some of the default parameters values
p = model_parameters(A, metabolicrate = ExtendedBA(:x, parameters_tuple = (T_opt = 290.0, activation_energy = 0.47))
```
"""
function ExtendedBA(rate_affected::Symbol; parameters_tuple...)
    if rate_affected ∉ [:growth, :r, :growthrate, :metabolism, :x, :metabolicrate, :attackrate]
        error("rate_affected should be either :growth (alternatively :r or :growthrate), :metabolism (alternatively :x or :metabolicrate), or :attackrate")
    end
    if rate_affected ∈ [:growth, :r, :growthrate]
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            extended_BA_r(passed_temp_parameters = pt)
        else
            extended_BA_r()
        end
    elseif rate_affected ∈ [:metabolism, :x, :metabolicrate]
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            extended_BA_x(passed_temp_parameters = pt)
        else
            extended_BA_x()
        end
    elseif rate_affected == :attackrate
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            extended_BA_attackr(passed_temp_parameters = pt)
        else
            extended_BA_attackr()
        end
    end
end

"""
**Gaussian**

    Gaussian(::Symbol, parameters_tuple...)

Calculates biological rates using a Gaussian function.

Can be called in `model_parameters` (see example) to define temperature dependent biological rates. Default parameters
values are provided, but can be overwritten through the keyword argument `parameters_tuple`. See the complete documentation
for the default values, their sources and more details.

Arguments:
- rate_affected : the biological rate that the function should calculate. May be :growth, :metabolism, :handlingtime, :attackrate
- parameters_tuple : a named tuple specifying parameters values. Note that is you provide a one-element tuple, it should end with a comma to avoid being treated as a vector.

# Example
```julia-repl
# Growth rate with default parameters values
p = model_parameters(A, growthrate = Gaussian(:growthrate))
# Metabolic rate, changing some of the default parameters values
p = model_parameters(A, metabolicrate = ExtendedBA(:x, parameters_tuple = (range_invertebrate = 25.0, ))
```
"""
function Gaussian(rate_affected::Symbol; parameters_tuple...)
    if rate_affected ∉ [:growth, :r, :growthrate, :metabolism, :x, :metabolicrate, :handlingtime, :attackrate]
        error("rate_affected should be either :growth (alternatively :r or :growthrate), :metabolism (alternatively :x or :metabolicrate) or :functionalresponse")
    end
    if rate_affected ∈ [:growth, :r, :growthrate]
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            gaussian_r(passed_temp_parameters = pt)
        else
            gaussian_r()
        end
    elseif rate_affected ∈ [:metabolism, :x, :metabolicrate]
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            gaussian_x(passed_temp_parameters = pt)
        else
            gaussian_x()
        end
    elseif rate_affected ∈ [:attackrate]
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            gaussian_attackr(passed_temp_parameters = pt)
        else
            gaussian_attackr()
        end
    elseif rate_affected ∈ [:handlingtime]
        if length(parameters_tuple) != 0
            pt = parameters_tuple[:parameters_tuple]
            gaussian_handlingt(passed_temp_parameters = pt)
        else
            gaussian_handlingt()
        end
    end
end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
import UncertainData: UncertainDataset, UncertainIndexDataset, UncertainValueDataset, resample

randomshuffle(x::UncertainDataset) = randomshuffle(resample(x))
randomphases(x::UncertainDataset) = randomphases(resample(x))
randomamplitudes(x::UncertainDataset) = randomamplitudes(resample(x))
aaft(x::UncertainDataset) = aaft(resample(x))
iaaft(x::UncertainDataset, args...; kwargs...) = iaaft(resample(x), args...; kwargs...)

randomshuffle(x::UncertainValueDataset) = randomshuffle(resample(x))
randomphases(x::UncertainValueDataset) = randomphases(resample(x))
randomamplitudes(x::UncertainValueDataset) = randomamplitudes(resample(x))
aaft(x::UncertainValueDataset) = aaft(resample(x))
iaaft(x::UncertainValueDataset, args...; kwargs...) = iaaft(resample(x), args...; kwargs...)

randomshuffle(x::UncertainIndexDataset) = randomshuffle(resample(x))
randomphases(x::UncertainIndexDataset) = randomphases(resample(x))
randomamplitudes(x::UncertainIndexDataset) = randomamplitudes(resample(x))
aaft(x::UncertainIndexDataset) = aaft(resample(x))
iaaft(x::UncertainIndexDataset, args...; kwargs...) = iaaft(resample(x), args...; kwargs...)


text:
function test_twodnavierstokes_lambdipole(n, dt, dev::Device=CPU(); L=2π, Ue=1, Re=L/20, ν=0.0, nν=1, ti=L/Ue*0.01, nm=3)
  nt = round(Int, ti/dt)
  prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, dt=dt, stepper="FilteredRK4")
  zeta₀ = lambdipole(Ue, Re, prob.grid)
  TwoDNavierStokes.set_zeta!(prob, zeta₀)

  xzeta = zeros(nm) # centroid of abs(zeta)
  Ue_m = zeros(nm)  # measured dipole speed
  x, y = gridpoints(prob.grid)
  zeta = prob.vars.zeta

  for i = 1:nm # step forward
    stepforward!(prob, nt)
    TwoDNavierStokes.updatevars!(prob)
    xzeta[i] = mean(@. abs(zeta)*x) / mean(abs.(zeta))
    if i > 1
      Ue_m[i] = (xzeta[i]-xzeta[i-1]) / ((nt-1)*dt)
    end
  end
  isapprox(Ue, mean(Ue_m[2:end]), rtol=rtol_lambdipole)
end

function test_twodnavierstokes_stochasticforcing_energybudget(dev::Device=CPU(); n=256, L=2π, dt=0.005, ν=1e-7, nν=2, μ=1e-1, nμ=0, tf=0.1/μ)
  nt = round(Int, tf/dt)

  # Forcing parameters
  kf, dkf = 12.0, 2.0
  ε = 0.1

  grid = TwoDGrid(dev, n, L)
  x, y = gridpoints(grid)

  Kr = ArrayType(dev)([CUDA.@allowscalar grid.kr[i] for i=1:grid.nkr, j=1:grid.nl])

  forcing_spectrum = ArrayType(dev)(zero(grid.Krsq))
  @. forcing_spectrum = exp(-(sqrt(grid.Krsq) - kf)^2 / (2 * dkf^2))
  @. forcing_spectrum = ifelse(grid.Krsq < 2^2, 0, forcing_spectrum)
  @. forcing_spectrum = ifelse(grid.Krsq > 20^2, 0, forcing_spectrum)
  @. forcing_spectrum = ifelse(Kr < 2π/L, 0, forcing_spectrum)
  ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly)
  forcing_spectrum .= ε / ε0 * forcing_spectrum

  Random.seed!(1234)

  function calcF!(Fh, sol, t, clock, vars, params, grid)
    eta = ArrayType(dev)(exp.(2π * im * rand(Float64, size(sol))) / sqrt(clock.dt))
    CUDA.@allowscalar eta[1, 1] = 0.0
    @. Fh = eta * sqrt(forcing_spectrum)
    return nothing
  end

  prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt,
   stepper="RK4", calcF=calcF!, stochastic=true)

  TwoDNavierStokes.set_zeta!(prob, 0*x)

  E = Diagnostic(TwoDNavierStokes.energy,             prob, nsteps=nt)
  D = Diagnostic(TwoDNavierStokes.energy_dissipation, prob, nsteps=nt)
  R = Diagnostic(TwoDNavierStokes.energy_drag,        prob, nsteps=nt)
  W = Diagnostic(TwoDNavierStokes.energy_work,        prob, nsteps=nt)
  diags = [E, D, W, R]

  stepforward!(prob, diags, nt)

  TwoDNavierStokes.updatevars!(prob)

  dEdt_numerical = (E[2:E.i] - E[1:E.i-1]) / prob.clock.dt

  # If the Ito interpretation was used for the work
  # then we need to add the drift term
  # dEdt_computed = W[2:E.i] + ε - D[1:E.i-1] - R[1:E.i-1]      # Ito
  dEdt_computed = W[2:E.i] - D[1:E.i-1] - R[1:E.i-1]        # Stratonovich

  return isapprox(dEdt_numerical, dEdt_computed, rtol = 1e-3)
end

function test_twodnavierstokes_stochasticforcing_enstrophybudget(dev::Device=CPU(); n=256, L=2π, dt=0.005, ν=1e-7, nν=2, μ=1e-1, nμ=0, tf=0.1/μ)
  nt = round(Int, tf/dt)

  # Forcing parameters
  kf, dkf = 12.0, 2.0
  εᶻ = 0.1

  grid = TwoDGrid(dev, n, L)
  x, y = gridpoints(grid)

  Kr = ArrayType(dev)([CUDA.@allowscalar grid.kr[i] for i=1:grid.nkr, j=1:grid.nl])

  forcing_spectrum = ArrayType(dev)(zero(grid.Krsq))
  @. forcing_spectrum = exp(-(sqrt(grid.Krsq) - kf)^2 / (2 * dkf^2))
  @. forcing_spectrum = ifelse(grid.Krsq < 2^2, 0, forcing_spectrum)
  @. forcing_spectrum = ifelse(grid.Krsq > 20^2, 0, forcing_spectrum)
  @. forcing_spectrum = ifelse(Kr < 2π/L, 0, forcing_spectrum)
  εᶻ0 = parsevalsum(forcing_spectrum / 2, grid) / (grid.Lx * grid.Ly)
  forcing_spectrum .= εᶻ / εᶻ0 * forcing_spectrum

  Random.seed!(1234)

  function calcF!(Fh, sol, t, cl, v, p, g)
    eta = ArrayType(dev)(exp.(2π * im * rand(Float64, size(sol))) / sqrt(cl.dt))
    CUDA.@allowscalar eta[1, 1] = 0.0
    @. Fh = eta * sqrt(forcing_spectrum)
    nothing
  end

  prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt,
   stepper="RK4", calcF=calcF!, stochastic=true)

  TwoDNavierStokes.set_zeta!(prob, 0*x)
  Z = Diagnostic(TwoDNavierStokes.enstrophy,             prob, nsteps=nt)
  D = Diagnostic(TwoDNavierStokes.enstrophy_dissipation, prob, nsteps=nt)
  R = Diagnostic(TwoDNavierStokes.enstrophy_drag,        prob, nsteps=nt)
  W = Diagnostic(TwoDNavierStokes.enstrophy_work,        prob, nsteps=nt)
  diags = [Z, D, W, R]

  stepforward!(prob, diags, nt)

  TwoDNavierStokes.updatevars!(prob)

  dZdt_numerical = (Z[2:Z.i] - Z[1:Z.i-1]) / prob.clock.dt

  # If the Ito interpretation was used for the work
  # then we need to add the drift term
  # dZdt_computed = W[2:Z.i] + εᶻ - D[1:Z.i-1] - R[1:Z.i-1]      # Ito
  dZdt_computed = W[2:Z.i] - D[1:Z.i-1] - R[1:Z.i-1]        # Stratonovich

  return isapprox(dZdt_numerical, dZdt_computed, rtol = 1e-3)
end

function test_twodnavierstokes_deterministicforcing_energybudget(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1, nμ=0)
  n, L  = 256, 2π
  ν, nν = 1e-7, 2
  μ, nμ = 1e-1, 0
  dt, tf = 0.005, 0.1/μ
  nt = round(Int, tf/dt)

  grid = TwoDGrid(dev, n, L)
  x, y = gridpoints(grid)

  # Forcing = 0.01cos(4x)cos(5y)cos(2t)
  f = @. 0.01cos(4x) * cos(5y)
  fh = rfft(f)
  function calcF!(Fh, sol, t, clock, vars, params, grid)
    @. Fh = fh * cos(2t)
    return nothing
  end

  prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt,
   stepper="RK4", calcF=calcF!, stochastic=false)

  TwoDNavierStokes.set_zeta!(prob, 0*x)

  E = Diagnostic(TwoDNavierStokes.energy,      prob, nsteps=nt)
  D = Diagnostic(TwoDNavierStokes.energy_dissipation, prob, nsteps=nt)
  R = Diagnostic(TwoDNavierStokes.energy_drag,        prob, nsteps=nt)
  W = Diagnostic(TwoDNavierStokes.energy_work,        prob, nsteps=nt)
  diags = [E, D, W, R]

  stepforward!(prob, diags, nt)

  TwoDNavierStokes.updatevars!(prob)

  dEdt_numerical = (E[3:E.i] - E[1:E.i-2]) / (2 * prob.clock.dt)
  dEdt_computed  = W[2:E.i-1] - D[2:E.i-1] - R[2:E.i-1]

  return isapprox(dEdt_numerical, dEdt_computed, rtol = 1e-4)
end

function test_twodnavierstokes_deterministicforcing_enstrophybudget(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1, nμ=0)
  n, L  = 256, 2π
  ν, nν = 1e-7, 2
  μ, nμ = 1e-1, 0
  dt, tf = 0.005, 0.1/μ
  nt = round(Int, tf/dt)

  grid = TwoDGrid(dev, n, L)
  x, y = gridpoints(grid)

  # Forcing = 0.01cos(4x)cos(5y)cos(2t)
  f = @. 0.01cos(4x) * cos(5y)
  fh = rfft(f)
  function calcF!(Fh, sol, t, clock, vars, params, grid)
    @. Fh = fh * cos(2t)
    return nothing
  end

  prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt,
   stepper="RK4", calcF=calcF!, stochastic=false)

  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

  TwoDNavierStokes.set_zeta!(prob, 0*x)

  Z = Diagnostic(TwoDNavierStokes.enstrophy,      prob, nsteps=nt)
  D = Diagnostic(TwoDNavierStokes.enstrophy_dissipation, prob, nsteps=nt)
  R = Diagnostic(TwoDNavierStokes.enstrophy_drag,        prob, nsteps=nt)
  W = Diagnostic(TwoDNavierStokes.enstrophy_work,        prob, nsteps=nt)
  diags = [Z, D, W, R]

  stepforward!(prob, diags, nt)

  TwoDNavierStokes.updatevars!(prob)

  dZdt_numerical = (Z[3:Z.i] - Z[1:Z.i-2]) / (2 * prob.clock.dt)
  dZdt_computed  = W[2:Z.i-1] - D[2:Z.i-1] - R[2:Z.i-1]

  return isapprox(dZdt_numerical, dZdt_computed, rtol = 1e-4)
end

"""
    testnonlinearterms(dt, stepper; kwargs...)

Tests the advection term in the twodnavierstokes module by timestepping a
test problem with timestep dt and timestepper identified by the string stepper.
The test problem is derived by picking a solution ζf (with associated
streamfunction ψf) for which the advection term J(ψf, ζf) is non-zero. Next, a
forcing Ff is derived according to Ff = ∂ζf/∂t + J(ψf, ζf) - ν∇²ζf. One solution
to the vorticity equation forced by this Ff is then ζf. (This solution may not
be realized, at least at long times, if it is unstable.)
"""
function test_twodnavierstokes_advection(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=1e-2, nν=1, μ=0.0, nμ=0)
  n, L  = 128, 2π
  ν, nν = 1e-2, 1
  μ, nμ = 0.0, 0
  tf = 1.0
  nt = round(Int, tf/dt)

  grid = TwoDGrid(dev, n, L)
  x, y = gridpoints(grid)

   psif = @.   sin(2x)*cos(2y) +  2sin(x)*cos(3y)
  zetaf = @. -8sin(2x)*cos(2y) - 20sin(x)*cos(3y)

  Ff = @. -(
    ν*( 64sin(2x)*cos(2y) + 200sin(x)*cos(3y) )
    + 8*( cos(x)*cos(3y)*sin(2x)*sin(2y) - 3cos(2x)*cos(2y)*sin(x)*sin(3y) )
  )

  Ffh = rfft(Ff)

  # Forcing
  function calcF!(Fh, sol, t, cl, v, p, g)
    Fh .= Ffh
    return nothing
  end

  prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt, stepper=stepper, calcF=calcF!, stochastic=false)
  sol, cl, p, v, g = prob.sol, prob.clock, prob.params, prob.vars, prob.grid
  TwoDNavierStokes.set_zeta!(prob, zetaf)

  stepforward!(prob, nt)
  TwoDNavierStokes.updatevars!(prob)

  isapprox(prob.vars.zeta, zetaf, rtol=rtol_twodnavierstokes)
end

function test_twodnavierstokes_energyenstrophy(dev::Device=CPU())
  nx, Lx  = 128, 2π
  ny, Ly  = 128, 3π

  grid = TwoDGrid(dev, nx, Lx, ny, Ly)
  x, y = gridpoints(grid)

  k₀, l₀ = 2π/grid.Lx, 2π/grid.Ly # fundamental wavenumbers
   ψ₀ = @. sin(2k₀*x)*cos(2l₀*y) + 2sin(k₀*x)*cos(3l₀*y)
   ζ₀ = @. -((2k₀)^2+(2l₀)^2)*sin(2k₀*x)*cos(2l₀*y) - (k₀^2+(3l₀)^2)*2sin(k₀*x)*cos(3l₀*y)

  energy_calc = 29/9
  enstrophy_calc = 2701/162

  prob = TwoDNavierStokes.Problem(dev; nx=nx, Lx=Lx, ny=ny, Ly=Ly, stepper="ForwardEuler")

  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid;

  TwoDNavierStokes.set_zeta!(prob, ζ₀)
  TwoDNavierStokes.updatevars!(prob)

  energyζ₀ = TwoDNavierStokes.energy(prob)
  enstrophyζ₀ = TwoDNavierStokes.enstrophy(prob)

  params = TwoDNavierStokes.Params(p.ν, p.nν)

  (isapprox(energyζ₀, energy_calc, rtol=rtol_twodnavierstokes) &&
   isapprox(enstrophyζ₀, enstrophy_calc, rtol=rtol_twodnavierstokes) &&
   TwoDNavierStokes.addforcing!(prob.timestepper.N, sol, cl.t, cl, v, p, g)==nothing && p == params)
end

function test_twodnavierstokes_problemtype(dev, T)
  prob = TwoDNavierStokes.Problem(dev; T=T)

  A = ArrayType(dev)

  (typeof(prob.sol)<:A{Complex{T},2} && typeof(prob.grid.Lx)==T && eltype(prob.grid.x)==T && typeof(prob.vars.u)<:A{T,2})
end


text:
using Revise
using DiffPrecTest
using Statistics, StatsBase, LinearAlgebra
using SparseArrays
using ProximalBase, CoordinateDescent, CovSel
using Random, Distributions
using JLD

#using Plots

pArr = [13, 200]
elemArr = [(5,5), (3, 2), (50, 25), (21, 20), (30, 30)]
n = 300
est      = Array{Any}(undef, 5)   # number of methods

out = Vector{DiffPrecResultNormal}(undef, 1000)

ip = 1
iElem = 2

p = pArr[ip]
Random.seed!(1234)

# generate model
Ωx = Matrix{Float64}(I, p, p)
for l=1:p-1
    Ωx[l  , l+1] = 0.6
    Ωx[l+1, l  ] = 0.6
end
for l=1:p-2
    Ωx[l  , l+2] = 0.3
    Ωx[l+2, l  ] = 0.3
end
Δ = zeros(p, p)
# generate Delta
for j=1:p
    #i = 5 + (j-1)
    i = j
    Δ[i, i] = rand(Uniform(0.1, 0.2))
end
for j=1:p-1
    #i = 5 + (j-1)
    i = j
    v = rand(Uniform(0.2, 0.5))
    Δ[i  , i+1] = v
    Δ[i+1, i  ] = v
end
Ωy = Ωx - Δ
d = Vector{Float64}(undef, p)
rand!(Uniform(0.5, 2.5), d)
d .= sqrt.(d)
D = Diagonal(d)
Σx = inv(Symmetric(D * Ωx * D))
Σy = inv(Symmetric(D * Ωy * D))

tΔ = D * Δ * D

dist_X = MvNormal(convert(Matrix, Σx))
dist_Y = MvNormal(convert(Matrix, Σy))

ri, ci = elemArr[iElem]
indE = sub2indLowerTriangular(p, ri, ci)

indOracle = Vector{Int64}()
for ci=1:p
  for ri=ci:p
    global indOracle

    if !iszero(tΔ[ri, ci])
      push!(indOracle, sub2indLowerTriangular(p, ri, ci))
    end
  end
end

Random.seed!(1234)
X = rand(dist_X, n)'
Y = rand(dist_Y, n)'

# Sx = Symmetric(cov(X))
# Sy = Symmetric(cov(Y))

# C = DiffPrecTest.skron(Sx, Sy)
# b = DiffPrecTest.svec(Sy) - DiffPrecTest.svec(Sx)

# f = CDQuadraticLoss(C, b)

@time res1, _ = DiffPrecTest.estimate(ReducedNormal(), X, Y, indE);
@time res2, _ = DiffPrecTest.estimate(SymmetricNormal(), X, Y, 3, 2);

for rep = 1:1000
  global out
  global indOracle
  global indE

  # generate data
  # Random.seed!(1234 + rep)
  X = rand(dist_X, n)'
  Y = rand(dist_Y, n)'

  #@time out[rep], _, _, _, _, indS = DiffPrecTest.estimate(ReducedNormal(), X, Y, indE)
  #@time out[rep], _, _, _, _, indS = DiffPrecTest.estimate(SymmetricNormal(), X, Y, 3, 2)

  @time out[rep] = DiffPrecTest.estimate(ReducedOracleNormal(), Symmetric(cov(X)), Symmetric(cov(Y)), X, Y, indE, indOracle)
end

computeSimulationResult(out, tΔ[elemArr[iElem]...])

true_v = tΔ[elemArr[iElem]...]
histogram( [(out[i].p - true_v) / out[i].std  for i = 1:1000] )

# @time est[2] = DiffPrecTest.estimate(SeparateNormal(), X, Y, indE)
# @time est[3] = DiffPrecTest.estimate(SymmetricOracleBoot(), X, Y, indS)
# @time est[4] = DiffPrecTest.estimate(SymmetricOracleNormal(), Symmetric(cov(X)), n, Symmetric(cov(Y)), n, indOracle)
# @time est[5] = DiffPrecTest.estimate(SymmetricOracleBoot(), X, Y, indOracle)

######################
# test variance

Random.seed!(1234+5)
X = rand(dist_X, n)'
Y = rand(dist_Y, n)'

Sx = Symmetric(cov(X))
Sy = Symmetric(cov(Y))

indS = indOracle

C = DiffPrecTest.skron(Sx, Sy, indS)
b = DiffPrecTest.svec(Sy, indS) - DiffPrecTest.svec(Sx, indS)

tmp = zeros(length(indS))
tmp[1] = 1.
ω = C \ tmp
Δ = C \ b

@time v = variance(ReducedOracleNormal(), Sx, Sy, X, Y, ω, indS, Δ, indS)

@time variance(DiffPrecTest.createVarianceReducedEstim(ω, indS, Δ, indS), X, Y)


text:
# 1D diffusion problem

# Packages and inclusions
using ModelingToolkit,DiffEqOperators,LinearAlgebra,Test,OrdinaryDiffEq
using ModelingToolkit: Differential

# Beam Equation
@test_broken begin
    @parameters x, t
    @variables u(..)
    Dt = Differential(t)
    Dtt = Differential(t)^2
    Dx = Differential(x)
    Dxx = Differential(x)^2
    Dx3 = Differential(x)^3
    Dx4 = Differential(x)^4

    g = -9.81
    EI = 1
    mu = 1
    L = 10.0
    dx = 0.4

    eq = Dtt(u(t,x)) ~ -mu*EI*Dx4(u(t,x)) + mu*g

    bcs = [u(0, x) ~ 0,
           u(t,0) ~ 0,
           Dx(u(t,0)) ~ 0,
           Dxx(u(t, L)) ~ 0,
           Dx3(u(t, L)) ~ 0]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,1.0),
               x ∈ IntervalDomain(0.0,L)]

    pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])
    discretization = MOLFiniteDifference([x=>dx],t, centered_order=4)
    prob = discretize(pdesys,discretization)
end

# Beam Equation with Velocity
@test_broken begin
    @parameters x, t
    @variables u(..), v(..)
    Dt = Differential(t)
    Dtt = Differential(t)^2
    Dx = Differential(x)
    Dxx = Differential(x)^2
    Dx3 = Differential(x)^3
    Dx4 = Differential(x)^4

    g = -9.81
    EI = 1
    mu = 1
    L = 10.0
    dx = 0.4

    eqs = [v(t, x) ~ Dt(u(t,x)),
           Dt(v(t,x)) ~ -mu*EI*Dx4(u(t,x)) + mu*g]

    bcs = [u(0, x) ~ 0,
           v(0, x) ~ 0,
           u(t,0) ~ 0,
           v(t,0) ~ 0,
           Dxx(u(t, L)) ~ 0,
           Dx3(u(t, L)) ~ 0]

    # Space and time domains
    domains = [t ∈ IntervalDomain(0.0,1.0),
               x ∈ IntervalDomain(0.0,L)]

    pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t,x)])
    discretization = MOLFiniteDifference([x=>dx],t, centered_order=4)
    prob = discretize(pdesys,discretization)
end

@testset "Kuramoto–Sivashinsky equation" begin
    @parameters x, t
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dx2 = Differential(x)^2
    Dx3 = Differential(x)^3
    Dx4 = Differential(x)^4

    α = 1
    β = 4
    γ = 1
    eq = Dt(u(x,t)) ~ -u(x,t)*Dx(u(x,t)) - α*Dx2(u(x,t)) - β*Dx3(u(x,t)) - γ*Dx4(u(x,t))

    u_analytic(x,t;z = -x/2+t) = 11 + 15*tanh(z) -15*tanh(z)^2 - 15*tanh(z)^3
    du(x,t;z = -x/2+t) = 15/2*(tanh(z) + 1)*(3*tanh(z) - 1)*sech(z)^2

    bcs = [u(x,0) ~ u_analytic(x,0),
           u(-10,t) ~ u_analytic(-10,t),
           u(10,t) ~ u_analytic(10,t),
           Dx(u(-10,t)) ~ du(-10,t),
           Dx(u(10,t)) ~ du(10,t)]

    # Space and time domains
    domains = [x ∈ IntervalDomain(-10.0,10.0),
               t ∈ IntervalDomain(0.0,1.0)]
    # Discretization
    dx = 0.4; dt = 0.2

    discretization = MOLFiniteDifference([x=>dx],t;centered_order=4,grid_align=center_align)
    pdesys = PDESystem(eq,bcs,domains,[x,t],[u(x,t)])
    prob = discretize(pdesys,discretization)

    sol = solve(prob,Tsit5(),saveat=0.1,dt=dt)

    @test sol.retcode == :Success

    xs = domains[1].domain.lower+dx+dx:dx:domains[1].domain.upper-dx-dx
    ts = sol.t

    u_predict = sol.u
    u_real = [[u_analytic(x, t) for x in xs] for t in ts]
    u_diff = u_real - u_predict
    @test_broken u_diff[:] ≈ zeros(length(u_diff)) atol=0.01;
    #plot(xs, u_diff)
end


text:
for line in split(strip(readstring("REQUIRE")), '\n')[2:end]
    name = split(line)[1]
    println(name, "\t", Pkg.installed(name))
end


text:
# This file is part of Kpax3. License is MIT.

function test_mcmc_clusterweight()
  n = 10
  k = 2

  priorR = Kpax3.EwensPitman(0.5, -0.1)
  @test isapprox(Kpax3.clusterweight(n, priorR), log(n - 0.5), atol=ε)

  priorR = Kpax3.EwensPitman(0.5, 0.0)
  @test isapprox(Kpax3.clusterweight(n, priorR), log(n - 0.5), atol=ε)

  priorR = Kpax3.EwensPitman(0, 2.0)
  @test isapprox(Kpax3.clusterweight(n, priorR), log(n), atol=ε)

  priorR = Kpax3.EwensPitman(-2, 5)
  @test isapprox(Kpax3.clusterweight(n, priorR), log(n + 2.0), atol=ε)

  priorR = Kpax3.EwensPitman(0.5, -0.1)
  @test isapprox(Kpax3.clusterweight(n, k, priorR), log(0.5 * k - 0.1), atol=ε)

  priorR = Kpax3.EwensPitman(0.5, 0.0)
  @test isapprox(Kpax3.clusterweight(n, k, priorR), log(0.5 * k), atol=ε)

  priorR = Kpax3.EwensPitman(0, 2.0)
  @test isapprox(Kpax3.clusterweight(n, k, priorR), log(2), atol=ε)

  priorR = Kpax3.EwensPitman(-2, 5)
  @test isapprox(Kpax3.clusterweight(n, k, priorR), log(-2.0 * (k - 5)), atol=ε)

  nothing
end

test_mcmc_clusterweight()


text:
#--------------------------------------------------------------------
# DNSS.jl
# Soham 03-2022
# Simulate a scalar field in fixed Minkowski spacetime
#--------------------------------------------------------------------
using NLsolve, Printf, DoubleFloats

function psi(u::Real, v::Real)::Real
    t = v + u
    r = v - u
    if r == 0
        return cos(k * t) # Explicitly regularized for r = 0
    else
        return cos(k * t) * (sin(k* r)  / (k * r))
    end
end

function computeUboundary(PS::ProductSpace{S1, S2})::NTuple{1, Field{S2}} where {S1, S2}
    ϕ0 = extractUboundary(Field(PS, psi), :incoming)
    return (ϕ0, )
end

function computeVboundary(PS::ProductSpace{S1, S2})::NTuple{1, Field{S1}} where {S1, S2}
    ϕ0 = extractVboundary(Field(PS, psi), :incoming)
    return (ϕ0, )
end

function initialguess(PS::ProductSpace{S1, S2})::NTuple{1, Field{ProductSpace{S1, S2}}} where {S1, S2}
    ϕ0 = Field(PS, (u,v)->1)
    return (ϕ0, )
end

function computePatch(PS::ProductSpace{S1, S2}, ubnd::NTuple{1, Field{S2}}, vbnd::NTuple{1, Field{S1}})::NTuple{1, Field{ProductSpace{S1, S2}}} where {S1, S2}
    B = incomingboundary(PS)
    I = identity(PS)
    DU, DV = derivative(PS)
    ∂U = combineUVboundary(ubnd, vbnd, :incoming)

    function offaxis(U::NTuple{1, Field{S}})::NTuple{1, Field{S}} where {S}
        (ϕ, ) = U
        η  = Field(PS, (u,v)->v-u)
        a  = Field(PS, (u,v)->2)
        F3 = DU*(DV*ϕ) + (1/η)*(DU*η)*(DV*ϕ) + (1/η)*(DV*η)*(DU*ϕ)
        return (F3, )
    end

    function onaxis(U::NTuple{1, Field{S}})::NTuple{1, Field{S}} where {S}
        (ϕ, ) = U
        A3 = DU*(DV*ϕ)
        return (A3, )
    end

    function boundary(U::NTuple{1, Field{S}})::NTuple{1, Field{S}} where {S}
        return U - ∂U
    end

    function residual(U::NTuple{1, Field{S}})::NTuple{1, Field{S}} where {S}
        F = offaxis(U)
        enforceregularity!(F, onaxis(U))
        enforcebc!(F, boundary(U))
        return F
    end

    function f!(fvec::Array{T,1}, x::Array{T,1}) where {T}
        fvec[:] = reshapeFromTuple(residual(reshapeToTuple(PS, 1, x)))
    end

    s = nlsolve(f!, reshapeFromTuple(initialguess(PS)); method=solver, autodiff=:forward,
                show_trace=debug >= 3, ftol=1e-12, iterations=100)
    k = reshapeToTuple(PS, 1, s.zero)

    if debug >= 2
        println(" Converged? = ", converged(s))
    end
    return k
end

function deltapsi(U::NTuple{1, Field})::NTuple{1, Field} where {S}
    psinum   = project(first(U), prolongate(first(U).space))
    psiexact = Field(psinum.space, psi)
    return (psinum - psiexact, )
end

function pconv(range::StepRange, nh::Int=1)
    println("Testing p convergence")
    n_ = collect(range)
    l_ = zeros(size(n_))
    for index in CartesianIndices(n_)
        np = n_[index]
        params = Parameters((np, np), (nh,nh), urange, vrange, 1)
        q = distribute(params, computePatch, computeUboundary, computeVboundary)
        l_[index] = rmse(extract(map(deltapsi, q), 1))
        @printf("  n = %i, rmse = %e\n", n_[index], l_[index])
    end
    return (n_, l_)
end

function hconv(range::StepRange, np::Int=4)
    println("Testing h convergence")
    n_ = collect(range)
    l_ = zeros(size(n_))
    for index in CartesianIndices(n_)
        nh = n_[index] = 2^n_[index]
        params = Parameters((np, np), (nh, nh), urange, vrange, 1)
        l_[index] = rmse(extract(map(deltapsi, distribute(params, computePatch, computeUboundary, computeVboundary)), 1))
        l0 = index[1] > 1 ? l_[index[1] - 1] : 1
        @printf("  n = %3i, rmse[2h] / rmse[h] = %e\n", n_[index], l0 / l_[index])
    end
    return (n_, l_)
end

#-----------------------------------------
# Setup simulation grid parameters,
#-----------------------------------------
npoints  = (18, 18)
npatches = (1, 1)
urange   = (0.0, 4.0)
vrange   = (0.0, 4.0)
params   = Parameters(npoints, npatches, urange, vrange, 1)
solver   = :trust_region
debug    = 3
k        = 0.5

# Check the solution
# It is unstable, and the instability is most promiment near the axis and grows with time.
# However, these are with partiy conditions on the axis. If we impose the other condition on the axis, the solution
# looks more stable.
# TODO: Check mode fall-off?
S = distribute(params, computePatch, computeUboundary, computeVboundary, debug)
save("output/data/scalarwave_3D/k_04_np_18_nh_1/scalarwave_3D_psi", extract(S, 1))
save("output/data/scalarwave_3D/k_04_np_18_nh_1/scalarwave_3D_deltapsi", extract(map(deltapsi, S), 1))

#-----------------------------------------
# Now test h-p convergence of the error
#-----------------------------------------
# np, ep = pconv(2:2:26, 12)
# nh, eh = hconv(0:1:9, 4)

#-----------------------------------------
# Save all the data for solution and convergence
#-----------------------------------------
# save("output/data/scalar_wave_spherical_symmetry/scalar_waave_spherical_symmetry_off_axis_pconv", np, ep)
# save("output/data/scalar_wave_spherical_symmetry/scalar_waave_spherical_symmetry_off_axis_hconv", nh, eh)


text:
# TODO:
#   Figure out ARPACK's calling:
#     Symmetric Eigenvectors/values
#     Generalised Symmetric Eigenvectors/values
#   Created weighted svd for MCA

function translate!{T<:AbstractFloat}(A::Array{T}, b::T)
    @inbounds for i = 1:length(A)
        A[i] += b
    end
    A
end
translate!{T<:AbstractFloat}(b::T, A::Array{T}) = translate!(A, b)

function translate!{T<:AbstractFloat}(b::Vector{T}, A::Matrix{T})
    (n = size(A,1)) == length(b) || throw(DimensionMismatch("first dimension of A does not match length of b"))
    @inbounds for j = 1:size(A,2), i = 1:n
        A[i,j] += b[i]
    end
    A
end

function translate!{T<:AbstractFloat}(A::Matrix{T}, b::Vector{T})
    (n = size(A,2)) == length(b) || throw(DimensionMismatch("second dimension of A does not match length of b"))
    @inbounds for j = 1:n, i = 1:size(A,1)
        A[i,j] += b[j]
    end
    A
end

# A := A + ϵI
function perturb!{T<:AbstractFloat}(A::Matrix{T}, ϵ::T)
    (n = size(A,1)) == size(A,2) || throw(DimensionMismatch("Matrix A must be square."))
    @inbounds for i = 1:n
        A[i,i] += ϵ
    end
    A
end

# A := (1-α)A + αβ
function regularize!{T<:AbstractFloat}(A::Matrix{T}, α::T, β::T)
    (n = size(A,1)) == size(A,2) || throw(DimensionMismatch("Matrix A must be square."))
    @inbounds for j = 1:n
        for i = 1:n
            A[i,j] *= (one(T)-α)
        end
        A[j,j] += α*β
    end
    A
end

# Symmetrize the lower half of matrix S using the upper half of S
function syml!(S::Matrix)
    (p = size(S,1)) == size(S,2) || throw(ArgumentError("S ∈ ℝ$(p)×$(size(S, 2)) must be square"))
    @inbounds for j = 1:(p - 1), i = (j + 1):p
        S[i, j] = S[j, i]
    end
    S
end
syml(S::Matrix) = syml!(copy(S))

# Symmetrize the upper off-diagonal of matrix S using the lower half of S
function symu!(S::Matrix)
    (p = size(S,1)) == size(S,2) || throw(ArgumentError("S ∈ ℝ$(p)×$(size(S, 2)) must be square"))
    @inbounds for j = 2:p, i = 1:(j-1)
        S[i,j] = S[j,i]
    end
    S
end
symu(S::Matrix) = symu!(copy(S))

#========================
 Import Divide & Conquer BLAS Algorithm
========================#

const liblapack = Base.liblapack_name

import Base.blasfunc

import Base.LinAlg: BlasFloat, BlasInt, LAPACKException,
    DimensionMismatch, SingularException, PosDefException, chkstride1, chksquare

typealias BlasChar Char

#Generic LAPACK error handlers
macro assertargsok() #Handle only negative info codes - use only if positive info code is useful!
    :(info[1]<0 && throw(ArgumentError("invalid argument #$(-info[1]) to LAPACK call")))
end
macro lapackerror() #Handle all nonzero info codes
    :(info[1]>0 ? throw(LAPACKException(info[1])) : @assertargsok )
end

macro assertnonsingular()
    :(info[1]>0 && throw(SingularException(info[1])))
end
macro assertposdef()
    :(info[1]>0 && throw(PosDefException(info[1])))
end

#Check that upper/lower (for special matrices) is correctly specified
macro chkuplo()
    :((uplo=='U' || uplo=='L') || throw(ArgumentError("""invalid uplo = $uplo
Valid choices are 'U' (upper) or 'L' (lower).""")))
end

for (syevd, elty) in
    ((:dsyevd_, :Float64),
     (:ssyevd_, :Float32))
    @eval begin
        #       SUBROUTINE dsyevd( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, LIWORK, INFO )
        # *     .. Scalar Arguments ..
        #       CHARACTER          JOBZ, UPLO
        #       INTEGER            INFO, LDA, LIWORK, LWORK, N
        # *     .. Array Arguments ..
        #       INTEGER            IWORK( * )
        #       DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * )
        function syevd!(jobz::BlasChar, uplo::BlasChar, A::StridedMatrix{$elty})
            chkstride1(A)
            n = chksquare(A)
            W     = similar(A, $elty, n)
            work  = Array($elty, 1)
            lwork = convert(BlasInt, -1)
            iwork  = Array(BlasInt, 1)
            liwork = convert(BlasInt, -1)
            info  = Array(BlasInt, 1)
            for i in 1:2
                ccall(($(blasfunc(syevd)), liblapack), Void,
                     (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
                      Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
                      &jobz, &uplo, &n, A, &max(1,stride(A,2)), W, work, &lwork, iwork, &liwork, info)
                @lapackerror
                if lwork < 0
                    lwork = convert(BlasInt, real(work[1]))
                    work = Array($elty, lwork)
                    liwork = iwork[1]
                    iwork = Array(BlasInt, liwork)
                end
            end
            jobz=='V' ? (W, A) : W
        end
    end
end

#========================
  Eigen Solvers
========================#

# Singular Value Decomposition
#   Vectors in increasing order
function components_svd!{T<:AbstractFloat}(X::Matrix{T})
    _U, D, Vᵀ = LAPACK.gesdd!('S', X)
    for i = 1:length(D)
        D[i] = D[i]^2
    end
    (transpose(Vᵀ), D)
end

# Divide and conquer singular value decomposition
function components_eig!{T<:AbstractFloat}(S::Matrix{T})
    (p = size(S,1)) == size(S,2) || throw(DimensionMismatch("Matrix A must be square."))
    D, V = syevd!('V', 'U', S)  # S = VDVᵀ
    (V[:,end:-1:1], D[end:-1:1])
end

#tol::T = eps(T)*maximum(size(S_x))*maximum(S_x)
function components_geig!{T<:AbstractFloat}(S_m::Matrix{T}, S_x::Matrix{T})
    (p = size(S_x, 1)) == size(S_x, 2) || throw(DimensionMismatch("Covariance matrix for X must be square."))
    size(S_m, 2) == size(S_m, 2)       || throw(DimensionMismatch("Covariance matrix for M must be square."))
    p == size(S_m, 2)                  || throw(DimensionMismatch("Covariance matrices for X and M must be of the same order."))
    D, V, _U = LAPACK.sygvd!(1, 'V', 'U', S_m, S_x)
    (V[:,end:-1:1], D[end:-1:1])
end

# MCA - WIP weighted svd
#=
function wsvd!{T<:AbstractFloat}(X::Matrix{T}, Wu::Matrix{T}, Wv::Matrix{T})
    (n = size(X, 1)) == size(Wu, 1) == size(Wu, 2) || throw(DimensionMismatch("The order of Wu must match the number of rows of X."))
    (m = size(X, 2)) == size(Wv, 1) == size(Wv, 2) || throw(DimensionMismatch("The order of Wv must match the number of columns of X."))
    Λu, Qu = LAPACK.syev!('V', 'U', Wu)  # QΛQᵀ = Wu
    Λv, Qv = LAPACK.syev!('V', 'U', Wv)  # QΛQᵀ = Wv
    U, D, Vᵀ = LAPACK.gesdd!('S', S_u * X * S_v)
end
=#