text:
module Scratch
import Base: UUID
using Dates

export with_scratch_directory, scratch_dir, get_scratch!, delete_scratch!, clear_scratchspaces!, @get_scratch!

const SCRATCH_DIR_OVERRIDE = Ref{Union{String,Nothing}}(nothing)
"""
    with_scratch_directory(f::Function, scratch_dir::String)

Helper function to allow temporarily changing the scratch space directory.  When this is
set, no other directory will be searched for spaces, and new spaces will be created
within this directory.  Similarly, removing a scratch space will only effect the given
scratch directory.
"""
function with_scratch_directory(f::Function, scratch_dir::String)
    try
        SCRATCH_DIR_OVERRIDE[] = scratch_dir
        f()
    finally
        SCRATCH_DIR_OVERRIDE[] = nothing
    end
end

"""
    scratch_dir(args...)

Returns a path within the current depot's `scratchspaces` directory.  This location can
be overridden via `with_scratch_directory()`.
"""
function scratch_dir(args...)
    if SCRATCH_DIR_OVERRIDE[] === nothing
        return abspath(first(Base.DEPOT_PATH), "scratchspaces", args...)
    else
        # If we've been given an override, use _only_ that directory.
        return abspath(SCRATCH_DIR_OVERRIDE[], args...)
    end
end

const uuid_re = r"uuid\s*=\s*(?i)\"([0-9a-z]{8}-[0-9a-z]{4}-[0-9a-z]{4}-[0-9a-z]{4}-[0-9a-z]{12})\""

find_uuid(uuid::UUID) = uuid
find_uuid(mod::Module) = find_uuid(Base.PkgId(parent_pkg).uuid)
function find_uuid(::Nothing)
    # Try and see if the current project has a UUID
    project = Base.active_project()
    if project !== nothing && isfile(project)
        str = read(project, String)
        if (m = match(uuid_re, str); m !== nothing)
            return UUID(m[1])
        end
    end
    # If we still haven't found a UUID, fall back to the "global namespace"
    return UUID(UInt128(0))
end

"""
    scratch_path(pkg_uuid, key)

Common utility function to return the path of a scratch space, keyed by the given
parameters.  Users should use `get_scratch!()` for most user-facing usage.
"""
function scratch_path(pkg_uuid::UUID, key::AbstractString)
    return scratch_dir(string(pkg_uuid), key)
end

# Session-based space access time tracker
## Should perhaps keep track of find_project_file(UUID) instead
## but since you can only load a package once per Julia session,
## and since these timers are reset for every session, keeping
## track of the calling UUID should be good enough.
const scratch_access_timers = Dict{Tuple{UUID,String},Float64}()
"""
    track_scratch_access(pkg_uuid, scratch_path)

We need to keep track of who is using which spaces, so we know when it is advisable to
remove them during a GC.  We do this by attributing accesses of spaces to `Project.toml`
files in much the same way that package versions themselves are logged upon install, only
instead of having the project information implicitly available, we must rescue it out
from the currently-active Pkg Env.  If we cannot do that, it is because someone is doing
something weird like opening a space for a Pkg UUID that is not loadable, which we will
simply not track; that space will be reaped after the appropriate time in an orphanage.

If `pkg_uuid` is explicitly set to `nothing`, this space is treated as belonging to the
current project, or if that does not exist, the default global project located at
`Base.load_path_expand("@v#.#")`.

While package and artifact access tracking can be done at `add()`/`instantiate()` time,
we must do it at access time for spaces, as we have no declarative list of spaces that
a package may or may not access throughout its lifetime.  To avoid building up a
ludicrously large number of accesses through programs that e.g. call `get_scratch!()` in a
loop, we only write out usage information for each space once per day at most.
"""
function track_scratch_access(pkg_uuid::UUID, scratch_path::AbstractString)
    # Don't write this out more than once per day within the same Julia session.
    curr_time = time()
    if get(scratch_access_timers, (pkg_uuid, scratch_path), 0.0) >= curr_time - 60*60*24
        return
    end

    function find_project_file(pkg_uuid::UUID)
        # The simplest case (`pkg_uuid` == UUID(0)) simply attributes the space to
        # the active project, and if that does not exist, the  global depot environment,
        # which will never cause the space to be GC'ed because it has been removed,
        # as long as the global environment within the depot itself is intact.
        if pkg_uuid === UUID(UInt128(0))
            p = Base.active_project()
            if p !== nothing
                return p
            end
            return Base.load_path_expand("@v#.#")
        end

        # Otherwise, we attempt to find the source location of the package identified
        # by `pkg_uuid`, then find its owning `Project.toml`:
        for (p, m) in Base.loaded_modules
            if p.uuid == pkg_uuid
                source_path = Base.pathof(m)
                if source_path !== nothing
                    return Base.current_project(dirname(source_path))
                end
            end
        end

        # Finally, make one last desperate attempt and check if the
        # active project has our UUID
        if pkg_uuid === find_uuid(nothing)
            p = Base.active_project()
            if p !== nothing
                return p
            end
        end

        # If we couldn't find anything to attribute the space to, return `nothing`.
        return nothing
    end

    # We must decide which manifest to attribute this space to.
    project_file = find_project_file(pkg_uuid)

    # If we couldn't find one, skip out.
    if project_file === nothing || !ispath(project_file)
        return
    end

    # We manually format some simple TOML entries so that we don't have
    # to depend on the whole TOML writer stdlib.
    toml_entry = string(
        "[[\"", escape_string(abspath(scratch_path)), "\"]]\n",
        "time = ", string(now()), "Z\n",
        "parent_projects = [\"", escape_string(abspath(project_file)), "\"]\n",
    )
    usage_file = usage_toml()
    mkpath(dirname(usage_file))
    open(usage_file, append=true) do io
        write(io, toml_entry)
    end

    # Record that we did, in fact, write out the space access time
    scratch_access_timers[(pkg_uuid, scratch_path)] = curr_time
end

usage_toml() = joinpath(first(Base.DEPOT_PATH), "logs", "scratch_usage.toml")

# We clear the access timers from every entry referencing this path
# even if the calling package might not match. This is safer,
# since it only means that we might print out some extra entries
# to scratch_usage.toml instead of missing to record some usage.
function prune_timers!(path)
    for k in keys(scratch_access_timers)
        _, recorded_path = k
        if path == recorded_path
            delete!(scratch_access_timers, k)
        end
    end
    return nothing
end

"""
    get_scratch!(parent_pkg = nothing, key::AbstractString, calling_pkg = parent_pkg)

Returns the path to (or creates) a space.

If `parent_pkg` is given (either as a `UUID` or as a `Module`), the scratch space is
namespaced with that package's UUID, so that it will not conflict with any other space
with the same name but a different parent package UUID.  The space's lifecycle is tied
to the calling package, allowing the space to be garbage collected if all versions of the
package that used it have been removed.  By default, `parent_pkg` and `calling_pkg` are
the same, however in rare cases a package may become dependent on a scratch space that is
namespaced within another package, in such cases they should identify themselves as the
`calling_pkg` so that the scratch space's lifecycle is tied to that calling package.

If `parent_pkg` is not defined, or is a `Module` without a root UUID (e.g. `Main`,
`Base`, an anonymous module, etc...) the created scratch space is namespaced within the
global environment for the current version of Julia.

Scratch spaces are removed if all calling projects that have accessed them are removed.
As an example, if a scratch space is used by two versions of the same package but not a
newer version, when the two older versions are removed the scratch space may be garbage
collected.  See `Pkg.gc()` and `track_scratch_access()` for more details.
"""
function get_scratch!(parent_pkg::Union{Module,UUID,Nothing}, key::AbstractString,
                      calling_pkg::Union{Module,UUID,Nothing} = parent_pkg)
    # Verify that the key is valid (only needed here at construction time)
    if match(r"^[a-zA-Z0-9-\.]+$", key) === nothing
        throw(ArgumentError(
            "invalid key \"$key\": keys may only include a-z, A-Z, 0-9, - and ."
            ))
    end
    parent_pkg = find_uuid(parent_pkg)
    calling_pkg = find_uuid(calling_pkg)
    # Calculate the path and create the containing folder
    path = scratch_path(parent_pkg, key)
    mkpath(path)

    # We need to keep track of who is using which spaces, so we track usage in a log
    track_scratch_access(calling_pkg, path)
    return path
end
get_scratch!(key::AbstractString) = get_scratch!(nothing, key)

"""
    delete_scratch!(parent_pkg, key)

Explicitly deletes a scratch space created through `get_scratch!()`.
"""
function delete_scratch!(parent_pkg::Union{Module,UUID,Nothing}, key::AbstractString, )
    parent_pkg = find_uuid(parent_pkg)
    path = scratch_path(parent_pkg, key)
    rm(path; force=true, recursive=true)
    prune_timers!(path)
    return nothing
end
delete_scratch!(key::AbstractString) = delete_scratch!(nothing, key)

"""
    clear_scratchspaces!()

Delete all scratch spaces in the current depot.
"""
function clear_scratchspaces!()
    rm(scratch_dir(); force=true, recursive=true)
    empty!(scratch_access_timers)
    return nothing
end

"""
    clear_scratchspaces!(parent_pkg::Union{Module,UUID})

Delete all scratch spaces for the given package.
"""
function clear_scratchspaces!(parent_pkg::Union{Module,UUID,Nothing})
    parent_pkg = find_uuid(parent_pkg)
    if parent_pkg === UUID(UInt128(0))
        # TODO: Why not make this a way to clear the global scratchspace ??
        throw(ArgumentError("Cannot find owning package for module"))
    end
    parent_prefix = scratch_dir(string(parent_pkg))
    # First prune the access timers from all references to paths belonging to this namespace
    for (_, path) in keys(scratch_access_timers)
        if startswith(path, parent_prefix)
            prune_timers!(path)
        end
    end
    # Next, remove the whole namespace
    rm(parent_prefix; force=true, recursive=true)
    return nothing
end

"""
    @get_scratch!(key)

Convenience macro that gets/creates a scratch space with the given key and parented to
the package the calling module belongs to.  If the calling module does not belong to a
package, (e.g. it is `Main`, `Base`, an anonymous module, etc...) the UUID will be taken
to be `nothing`, creating a global scratchspace.
"""
macro get_scratch!(key)
    # Note that if someone uses this in the REPL, it will return `nothing`, and thereby
    # create a global scratch space.
    uuid = Base.PkgId(__module__).uuid
    return quote
        get_scratch!($(esc(uuid)), $(esc(key)), $(esc(uuid)))
    end
end

end # module Scratch


text:
path = dirname(@__FILE__)

m = AnSchorfheide()

@testset "Test prepare forecast inputs" begin
    # Throw error if input_type=:subset but no subset_inds provided
    @test_throws ErrorException DSGE.prepare_forecast_inputs!(m, :subset, :none, [:histobs, :forecastobs])
    # Throw error if trying to compute shockdec under alternative policy rule
    m <= Setting(:alternative_policy, AltPolicy(:taylor93, eqcond, solve))
    @test_throws ErrorException DSGE.prepare_forecast_inputs!(m, :mode, :none, [:shockdec])
    # Go back to normal policy rule
    m <= Setting(:alternative_policy, AltPolicy(:historical, eqcond, solve))
    # Test with df passed in
    if haskey(ENV, "FRED_API_KEY") || isfile(joinpath(homedir(), ".freddatarc"))
        df_in = load_data(m)
        global output_vars, df_out = DSGE.prepare_forecast_inputs!(m, :mode, :none, [:histobs, :forecastobs], df = df_in)
        @test output_vars == [:histobs, :forecastobs, :bddforecastobs]
        @test isapprox(df_to_matrix(m, df_in), df_to_matrix(m, df_out), atol = 0., nans = true)
        # Test without df passed in
        global output_vars, df_out = DSGE.prepare_forecast_inputs!(m, :mode, :none, [:histobs, :forecastobs])
        @test output_vars == [:histobs, :forecastobs, :bddforecastobs]
        @test isapprox(df_to_matrix(m, df_in), df_to_matrix(m, df_out), atol = 0., nans = true)
    else
        @warn "Skipping prepare forecast inputs test because FRED_API_KEY not present"
    end
end

@testset "Test load draws" begin
    global overrides = forecast_input_file_overrides(m)
    global overrides[:mode] = "$path/../reference/paramsmode_.h5"
    global overrides[:mode_draw_shocks] = "$path/../reference/paramsmode_.h5"

    # Init
    @test typeof(load_draws(m, :init, verbose = :none)) == Vector{Float64}
    @test typeof(load_draws(m, :init_draw_shocks, verbose = :none)) == Vector{Float64}
    @test typeof(load_draws(m, :prior, verbose = :none)) == Matrix{Float64}
    @test typeof(load_draws(m, :init_draw_shocks, 1:1, verbose = :none))==Vector{Vector{Float64}}
    @test typeof(load_draws(m, :prior, 1:1, verbose = :none))==Vector{Vector{Float64}}

    # MH
    m <= Setting(:sampling_method, :MH)
    global overrides[:full] = "$path/../reference/mhsave_.h5"
    # load_draws without blocks
    @test typeof(load_draws(m, :mode, verbose = :none)) == Vector{Float64}
    @test typeof(load_draws(m, :mode_draw_shocks, verbose = :none)) == Vector{Float64}
    @test typeof(load_draws(m, :full, verbose = :none)) == Matrix{Float64}
    @test typeof(load_draws(m, :subset, subset_inds = 1:10, verbose = :none)) == Matrix{Float64}

    # load_draws with blocks
    @test typeof(load_draws(m, :full, 1:1, verbose = :none))==Vector{Vector{Float64}}
    @test typeof(load_draws(m, :mode_draw_shocks, 1:1, verbose = :none))==Vector{Vector{Float64}}
       # can't block with empty/invalid block indices
    @test_throws ErrorException typeof(load_draws(m, :full, 1:0, verbose = :none))==Vector{Vector{Float64}}
       # can't block with only single draw
    @test_throws ErrorException load_draws(m, :mode, 1:1, verbose = :none)
    @test_throws ErrorException load_draws(m, :init, 1:1, verbose = :none)

    # SMC
    m <= Setting(:sampling_method, :SMC)
    global overrides[:full] = "$path/../reference/smcsave_.h5"
    # load draws without blocks
    @test typeof(load_draws(m, :mode, verbose = :none, use_highest_posterior_value = true)) == Vector{Float64}
    @test typeof(load_draws(m, :mode_draw_shocks, verbose = :none)) == Vector{Float64}
    @test typeof(load_draws(m, :full, verbose = :none)) == Matrix{Float64}
    @test typeof(load_draws(m, :subset, subset_inds = 1:10, verbose = :none, use_highest_posterior_value = true)) == Matrix{Float64}

    # load draws with blocks
    @test typeof(load_draws(m, :full, 1:1, verbose = :none))==Vector{Vector{Float64}}
    @test_throws ErrorException typeof(load_draws(m, :full, 1:0, verbose = :none))==Vector{Vector{Float64}}
    @test typeof(load_draws(m, :mode_draw_shocks, 1:1, verbose = :none))==Vector{Vector{Float64}}

    # If give nonsense sampler, should throw error
    m <= Setting(:sampling_method, :marco)
    @test_throws ErrorException load_draws(m, :mode, verbose = :none)

    # Test it can load the csminwel mode when using SMC
    delete!(overrides, :mode) # make sure no modal override
    m <= Setting(:sampling_method, :SMC)
    m <= Setting(:saveroot, joinpath(dirname(@__FILE__), "../reference"))
    m <= Setting(:data_vintage, "")
    @test typeof(load_draws(m, :mode, verbose = :none,
                            use_highest_posterior_value = true)) == Vector{Float64}
    @test typeof(load_draws(m, :mode, verbose = :none,
                            use_highest_posterior_value = false)) == Vector{Float64}
    m <= Setting(:fix, "true", true, "fix", "") # Add tail to estimation file name
    @test typeof(load_draws(m, :mode, verbose = :none, # to test we default to the cloud again.
                            use_highest_posterior_value = false)) == Vector{Float64}
end

ext: jl


text:
#=
    an excitation tensor has 4 legs (1,2),(3,4)
    the first and the last are virtual, the second is physical, the third is the utility leg
=#
@with_kw struct QuasiparticleAnsatz <: Algorithm
    toler::Float64 = 1e-10;
    krylovdim::Int = 30;
end

include("excitransfers.jl")

"""
    excitations(H::Hamiltonian, alg::QuasiparticleAnsatz, args...; kwargs...)

Compute the first excited states and their energy gap above a groundstate.

This is an implementation of the algorithm found [here](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.111.080401).
"""
function excitations end

################################################################################
#                           Infinite Excitations                               #
################################################################################
"""
    excitations(H, alg, V₀::InfiniteQP, lenvs, renvs; num, solver)

Optimise the infinite quasiparticle state ```V₀```.

# Arguments
- `V₀::InfiniteQP`: initial quasiparticle state.
- `lenvs=environments(V₀.left_gs, H; solver=solver)`: left environment of the groundstate.
- `renvs=(V₀.trivial ? lenvs : environments(V₀.right_gs, H; solver=solver))`: right environment of the groundstate.
- `num=1`: number of excited states to compute.
- `solver=Defaults.solver`: the algorithm to compute the environments.
"""
function excitations(
    H::Hamiltonian, alg::QuasiparticleAnsatz, V₀::InfiniteQP, lenvs, renvs;
    num=1, solver=Defaults.solver
)
    qp_envs(V) = environments(V, H, lenvs, renvs, solver=solver)
    H_eff = @closure(V->effective_excitation_hamiltonian(H, V, qp_envs(V)))

    Es, Vs, convhist = eigsolve(
        H_eff, V₀, num, :SR, tol=alg.toler, krylovdim=alg.krylovdim
    )

    convhist.converged < num &&
        @warn "Quasiparticle didn't converge: $(convhist.normres)"

    return Es, Vs
end
function excitations(
    H::Hamiltonian, alg::QuasiparticleAnsatz, V₀::InfiniteQP, lenvs;
    num=1, solver=Defaults.solver
)
    # Infer `renvs` in function body as it depends on `solver`.
    renvs = V₀.trivial ? lenvs : environments(V₀.right_gs, H, solver=solver)
    return excitations(H, alg, V₀, lenvs, renvs; num, solver)
end
function excitations(
    H::Hamiltonian, alg::QuasiparticleAnsatz, V₀::InfiniteQP;
    num=1, solver=Defaults.solver
)
    # Infer `lenvs` in function body as it depends on `solver`.
    lenvs = environments(V₀.left_gs, H, solver=solver)
    return excitations(H, alg, V₀, lenvs; num, solver)
end

"""
    excitations(H, alg, p::Real, lmps, lenvs, rmps, renvs; sector, kwargs...)

Create and optimise an infinite quasiparticle state with momentum ```p```.
"""
function excitations(
    H::Hamiltonian, alg::QuasiparticleAnsatz, p::Real,
    lmps::InfiniteMPS, lenvs=environments(lmps, H),
    rmps::InfiniteMPS=lmps, renvs=lmps===rmps ? lenvs : environments(rmps, H);
    sector=first(sectors(oneunit(virtualspace(lmps, 1)))), num=1,
    solver=Defaults.solver
)
    V₀ = LeftGaugedQP(rand, lmps, rmps; sector=sector, momentum=p)
    return excitations(H, alg, V₀, lenvs, renvs; num=num, solver=solver)
end
function excitations(
    H, alg::QuasiparticleAnsatz, momenta, lmps, lenvs=environments(lmps, H),
    rmps=lmps, renvs=lmps===rmps ? lenvs : environments(rmps, H);
    verbose=Defaults.verbose, num=1, solver=Defaults.solver, sector=first(sectors(oneunit(virtualspace(lmps, 1))))
)
    tasks = map(momenta) do p
        @Threads.spawn begin
            (E, V) = excitations(H, alg, p, lmps, lenvs, rmps, renvs; num=num, solver=solver, sector=sector)
            verbose && @info "Found excitations for p = $(p)"
            (E, V)
        end
    end

    fetched = fetch.(tasks);

    Ep = permutedims(reduce(hcat, map(x->x[1][1:num], fetched)));
    Bp = permutedims(reduce(hcat, map(x->x[2][1:num], fetched)));

    return Ep, Bp
end

################################################################################
#                           Finite Excitations                                 #
################################################################################
"""
    excitations(H, alg, V₀::FiniteQP[, lenvs[, renvs]]; num)

Optimise the finite quasiparticle state ```V₀```.

# Arguments
- `V₀::FiniteQP`: initial quasiparticle state.
- `lenvs=environments(V₀.left_gs, H)`: left environment of the groundstate.
- `renvs=(V₀.trivial ? lenvs : environments(V₀.right_gs, H))`: right environment of the groundstate.
- `num=1`: number of excited states to compute.
"""
function excitations(
    H::Hamiltonian, alg::QuasiparticleAnsatz,
    V₀::FiniteQP, lenvs=environments(V₀.left_gs, H),
    renvs=V₀.trivial ? lenvs : environments(V₀.right_gs, H);
    num=1
)
    qp_envs(V) = environments(V, H, lenvs, renvs)
    H_eff = @closure(V->effective_excitation_hamiltonian(H, V, qp_envs(V)))

    Es, Vs, convhist = eigsolve(
        H_eff, V₀, num, :SR, tol=alg.toler, krylovdim=alg.krylovdim
    )

    convhist.converged < num &&
        @warn "Quasiparticle didn't converge: $(convhist.normres)"

    return Es, Vs
end

"""
    excitations(H, alg, [lmps, [lenvs, [rmps, [renvs]]]]; sector, kwargs...)

Create and optimise a finite quasiparticle state.
"""
function excitations(
    H::Hamiltonian, alg::QuasiparticleAnsatz,
    lmps::FiniteMPS, lenvs=environments(lmps, H),
    rmps::FiniteMPS=lmps, renvs=lmps===rmps ? lenvs : environments(rmps, H)environmentsleft;
    sector=first(sectors(oneunit(virtualspace(lmps, 1)))), num=1
)
    V₀ = LeftGaugedQP(rand, lmps, rmps; sector=sector);
    return excitations(H, alg, V₀, lenvs, renvs; num=num)
end

################################################################################
#                           Statmech Excitations                               #
################################################################################

function excitations(
    H::MPOMultiline,alg::QuasiparticleAnsatz,V₀::Multiline{<:InfiniteQP},lenvs,renvs;
    num = 1, solver = Defaults.solver)

    qp_envs(V) = environments(V, H, lenvs, renvs, solver=solver)
    function H_eff(oV)
        V = Multiline(oV.vecs);
        RecursiveVec(effective_excitation_hamiltonian(H, V, qp_envs(V)).data.data)
    end

    Es, Vs, convhist = eigsolve(
        H_eff, RecursiveVec(V₀.data.data), num, :LM, tol=alg.toler, krylovdim=alg.krylovdim
    )

    convhist.converged < num &&
        @warn "Quasiparticle didn't converge: $(convhist.normres)"

    return Es, map(x->Multiline(x.vecs),Vs)
end

function excitations(
    H::MPOMultiline, alg::QuasiparticleAnsatz, V₀::Multiline{<:InfiniteQP}, lenvs;
    num=1, solver=Defaults.solver
)
    # Infer `renvs` in function body as it depends on `solver`.
    renvs = V₀.trivial ? lenvs : environments(V₀.right_gs, H, solver=solver)
    return excitations(H, alg, V₀, lenvs, renvs; num, solver)
end
function excitations(
    H::MPOMultiline, alg::QuasiparticleAnsatz, V₀::Multiline{<:InfiniteQP};
    num=1, solver=Defaults.solver
)
    # Infer `lenvs` in function body as it depends on `solver`.
    lenvs = environments(V₀.left_gs, H, solver=solver)
    return excitations(H, alg, V₀, lenvs; num, solver)
end

function excitations(
    H::InfiniteMPO, alg::QuasiparticleAnsatz, p::Real,
    lmps, lenvs=environments(lmps, H),
    rmps=lmps, renvs=lmps===rmps ? lenvs : environments(rmps, H);
    sector=first(sectors(oneunit(virtualspace(lmps, 1)))), num=1,
    solver=Defaults.solver
)
    multiline_lmps = convert(MPSMultiline,lmps);
    if lmps === rmps
        excitations(convert(MPOMultiline,H),alg,p,multiline_lmps,lenvs,multiline_lmps,
            lenvs;sector,num,solver);
    else
        excitations(convert(MPOMultiline,H),alg,p,multiline_lmps,lenvs,convert(MPSMultiline,rmps),
            renvs;sector,num,solver);
    end
end

function excitations(
    H::MPOMultiline, alg::QuasiparticleAnsatz, p::Real,
    lmps::MPSMultiline, lenvs=environments(lmps, H),
    rmps=lmps, renvs=lmps===rmps ? lenvs : environments(rmps, H);
    sector=first(sectors(oneunit(virtualspace(lmps, 1)))), num=1,
    solver=Defaults.solver
)
    V₀ = Multiline(map(1:size(lmps,1)) do row
        LeftGaugedQP(rand,lmps[row],rmps[row];sector,momentum=p)
    end)

    return excitations(H, alg, V₀, lenvs, renvs; num=num, solver=solver)
end

################################################################################
#                                H_eff                                         #
################################################################################

function effective_excitation_hamiltonian(
    H::MPOHamiltonian, exci::QP,envs=environments(exci, H)
)
    Bs = [exci[i] for i in 1:length(exci)];
    toret = similar(exci);

    #do necessary contractions
    for i = 1:length(exci)
        T = zero(Bs[i]);

        for (j,k) in keys(H,i)
            @plansor T[-1 -2;-3 -4] +=    leftenv(envs.lenvs,i,exci.left_gs)[j][-1 5;4]*
                                                Bs[i][4 2;-3 1]*
                                                H[i,j,k][5 -2;2 3]*
                                                rightenv(envs.renvs,i,exci.right_gs)[k][1 3;-4]

            # <B|H|B>-<H>
            en = @plansor    conj(exci.left_gs.AC[i][2 6;4])*
                            leftenv(envs.lenvs,i,exci.left_gs)[j][2 5;3]*
                            exci.left_gs.AC[i][3 7;1]*
                            H[i,j,k][5 6;7 8]*
                            rightenv(envs.lenvs,i,exci.left_gs)[k][1 8;4]

            T -= Bs[i]*en
            if i>1 || exci isa InfiniteQP
                @plansor T[-1 -2;-3 -4] +=    envs.lBs[i-1][j][-1 4;-3 5]*
                                                    exci.right_gs.AR[i][5 2;1]*
                                                    H[i,j,k][4 -2;2 3]*
                                                    rightenv(envs.renvs,i,exci.right_gs)[k][1 3;-4]
            end
            if i<length(exci.left_gs) || exci isa InfiniteQP
                @plansor T[-1 -2;-3 -4] +=    leftenv(envs.lenvs,i,exci.left_gs)[j][-1 2;1]*
                                                    exci.left_gs.AL[i][1 3;4]*
                                                    H[i,j,k][2 -2;3 5]*
                                                    envs.rBs[i+1][k][4 5;-3 -4]
            end
        end

        toret[i] = T;

    end

    return toret

end

function effective_excitation_hamiltonian(
    H::MPOMultiline, exci::Multiline{<:InfiniteQP},envs=environments(exci, H)
)

    toreturn = Multiline(similar.(exci.data));

    left_gs = exci.left_gs;
    right_gs = exci.right_gs;

    for row in 1:size(H,1)
        Bs = [exci[row][i] for i in 1:size(H,2)];
        for col in 1:size(H,2)
            en = @plansor    conj(left_gs.AC[row,col][2 6;4])*
                            leftenv(envs.lenvs,row,col,left_gs)[2 5;3]*
                            left_gs.AC[row+1,col][3 7;1]*
                            H[row,col][5 6;7 8]*
                            rightenv(envs.lenvs,row,col,left_gs)[1 8;4]

            @plansor T[-1 -2;-3 -4] := leftenv(envs.lenvs,row,col,left_gs)[-1 5;4]*
                                                Bs[col][4 2;-3 1]*
                                                H[row,col][5 -2;2 3]*
                                                rightenv(envs.renvs,row,col,right_gs)[1 3;-4]

            @plansor T[-1 -2;-3 -4] +=    envs.lBs[row][col-1][-1 4;-3 5]*
                                                right_gs.AR[row,col][5 2;1]*
                                                H[row,col][4 -2;2 3]*
                                                rightenv(envs.renvs,row,col,right_gs)[1 3;-4]

            @plansor T[-1 -2;-3 -4] +=    leftenv(envs.lenvs,row,col,left_gs)[-1 2;1]*
                                                left_gs.AL[row,col][1 3;4]*
                                                H[row,col][2 -2;3 5]*
                                                envs.rBs[row][col+1][4 5;-3 -4]

            toreturn[row+1][col] = T/en;
        end
    end

    return toreturn;
end


text:
export BrickSilentStepperSpeedRamping
struct BrickSilentStepperSpeedRamping
    acceleration::Integer
    deacceleration::Integer
end

export BrickSilentStepperStepConfiguration
struct BrickSilentStepperStepConfiguration
    step_resolution::Integer
    interpolation::Bool
end

export BrickSilentStepperBasicConfiguration
struct BrickSilentStepperBasicConfiguration
    standstill_current::Integer
    motor_run_current::Integer
    standstill_delay_time::Integer
    power_down_time::Integer
    stealth_threshold::Integer
    coolstep_threshold::Integer
    classic_threshold::Integer
    high_velocity_chopper_mode::Bool
end

export BrickSilentStepperSpreadcycleConfiguration
struct BrickSilentStepperSpreadcycleConfiguration
    slow_decay_duration::Integer
    enable_random_slow_decay::Bool
    fast_decay_duration::Integer
    hysteresis_start_value::Integer
    hysteresis_end_value::Integer
    sine_wave_offset::Integer
    chopper_mode::Integer
    comparator_blank_time::Integer
    fast_decay_without_comparator::Bool
end

export BrickSilentStepperStealthConfiguration
struct BrickSilentStepperStealthConfiguration
    enable_stealth::Bool
    amplitude::Integer
    gradient::Integer
    enable_autoscale::Bool
    force_symmetric::Bool
    freewheel_mode::Integer
end

export BrickSilentStepperCoolstepConfiguration
struct BrickSilentStepperCoolstepConfiguration
    minimum_stallguard_value::Integer
    maximum_stallguard_value::Integer
    current_up_step_width::Integer
    current_down_step_width::Integer
    minimum_current::Integer
    stallguard_threshold_value::Integer
    stallguard_mode::Integer
end

export BrickSilentStepperMiscConfiguration
struct BrickSilentStepperMiscConfiguration
    disable_short_to_ground_protection::Bool
    synchronize_phase_frequency::Integer
end

export BrickSilentStepperDriverStatus
struct BrickSilentStepperDriverStatus
    open_load::Integer
    short_to_ground::Integer
    over_temperature::Integer
    motor_stalled::Bool
    actual_motor_current::Integer
    full_step_active::Bool
    stallguard_result::Integer
    stealth_voltage_amplitude::Integer
end

export BrickSilentStepperAllData
struct BrickSilentStepperAllData
    current_velocity::Integer
    current_position::Integer
    remaining_steps::Integer
    stack_voltage::Integer
    external_voltage::Integer
    current_consumption::Integer
end

export BrickSilentStepperSPITFPBaudrateConfig
struct BrickSilentStepperSPITFPBaudrateConfig
    enable_dynamic_baudrate::Bool
    minimum_dynamic_baudrate::Integer
end

export BrickSilentStepperSPITFPErrorCount
struct BrickSilentStepperSPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickSilentStepperProtocol1BrickletName
struct BrickSilentStepperProtocol1BrickletName
    protocol_version::Integer
    firmware_version::Vector{Integer}
    name::String
end

export BrickSilentStepperIdentity
struct BrickSilentStepperIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickSilentStepper
"""
Silently drives one bipolar stepper motor with up to 46V and 1.6A per phase
"""
mutable struct BrickSilentStepper <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickSilentStepper(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.brick_silent_stepper")
        deviceInternal = package.BrickSilentStepper(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_max_velocity
"""
    $(SIGNATURES)

Sets the maximum velocity of the stepper motor.
This function does *not* start the motor, it merely sets the maximum
velocity the stepper motor is accelerated to. To get the motor running use
either :func:`Set Target Position`, :func:`Set Steps`, :func:`Drive Forward` or
:func:`Drive Backward`.
"""
function set_max_velocity(device::BrickSilentStepper, velocity)
    device.deviceInternal.set_max_velocity(velocity)
end

export get_max_velocity
"""
    $(SIGNATURES)

Returns the velocity as set by :func:`Set Max Velocity`.
"""
function get_max_velocity(device::BrickSilentStepper)
    return device.deviceInternal.get_max_velocity()
end

export get_current_velocity
"""
    $(SIGNATURES)

Returns the *current* velocity of the stepper motor.
"""
function get_current_velocity(device::BrickSilentStepper)
    return device.deviceInternal.get_current_velocity()
end

export set_speed_ramping
"""
    $(SIGNATURES)

Sets the acceleration and deacceleration of the stepper motor.
An acceleration of 1000 means, that
every second the velocity is increased by 1000 *steps/s*.

For example: If the current velocity is 0 and you want to accelerate to a
velocity of 8000 *steps/s* in 10 seconds, you should set an acceleration
of 800 *steps/s²*.

An acceleration/deacceleration of 0 means instantaneous
acceleration/deacceleration (not recommended)
"""
function set_speed_ramping(device::BrickSilentStepper, acceleration, deacceleration)
    device.deviceInternal.set_speed_ramping(acceleration, deacceleration)
end

export get_speed_ramping
"""
    $(SIGNATURES)

Returns the acceleration and deacceleration as set by
:func:`Set Speed Ramping`.
"""
function get_speed_ramping(device::BrickSilentStepper)
    return device.deviceInternal.get_speed_ramping()
end

export full_brake
"""
    $(SIGNATURES)

Executes an active full brake.

.. warning::
 This function is for emergency purposes,
 where an immediate brake is necessary. Depending on the current velocity and
 the strength of the motor, a full brake can be quite violent.

Call :func:`Stop` if you just want to stop the motor.
"""
function full_brake(device::BrickSilentStepper)
    device.deviceInternal.full_brake()
end

export set_current_position
"""
    $(SIGNATURES)

Sets the current steps of the internal step counter. This can be used to
set the current position to 0 when some kind of starting position
is reached (e.g. when a CNC machine reaches a corner).
"""
function set_current_position(device::BrickSilentStepper, position)
    device.deviceInternal.set_current_position(position)
end

export get_current_position
"""
    $(SIGNATURES)

Returns the current position of the stepper motor in steps. On startup
the position is 0. The steps are counted with all possible driving
functions (:func:`Set Target Position`, :func:`Set Steps`, :func:`Drive Forward` or
:func:`Drive Backward`). It also is possible to reset the steps to 0 or
set them to any other desired value with :func:`Set Current Position`.
"""
function get_current_position(device::BrickSilentStepper)
    return device.deviceInternal.get_current_position()
end

export set_target_position
"""
    $(SIGNATURES)

Sets the target position of the stepper motor in steps. For example,
if the current position of the motor is 500 and :func:`Set Target Position` is
called with 1000, the stepper motor will drive 500 steps forward. It will
use the velocity, acceleration and deacceleration as set by
:func:`Set Max Velocity` and :func:`Set Speed Ramping`.

A call of :func:`Set Target Position` with the parameter *x* is equivalent to
a call of :func:`Set Steps` with the parameter
(*x* - :func:`Get Current Position`).
"""
function set_target_position(device::BrickSilentStepper, position)
    device.deviceInternal.set_target_position(position)
end

export get_target_position
"""
    $(SIGNATURES)

Returns the last target position as set by :func:`Set Target Position`.
"""
function get_target_position(device::BrickSilentStepper)
    return device.deviceInternal.get_target_position()
end

export set_steps
"""
    $(SIGNATURES)

Sets the number of steps the stepper motor should run. Positive values
will drive the motor forward and negative values backward.
The velocity, acceleration and deacceleration as set by
:func:`Set Max Velocity` and :func:`Set Speed Ramping` will be used.
"""
function set_steps(device::BrickSilentStepper, steps)
    device.deviceInternal.set_steps(steps)
end

export get_steps
"""
    $(SIGNATURES)

Returns the last steps as set by :func:`Set Steps`.
"""
function get_steps(device::BrickSilentStepper)
    return device.deviceInternal.get_steps()
end

export get_remaining_steps
"""
    $(SIGNATURES)

Returns the remaining steps of the last call of :func:`Set Steps`.
For example, if :func:`Set Steps` is called with 2000 and
:func:`Get Remaining Steps` is called after the motor has run for 500 steps,
it will return 1500.
"""
function get_remaining_steps(device::BrickSilentStepper)
    return device.deviceInternal.get_remaining_steps()
end

export set_step_configuration
"""
    $(SIGNATURES)

Sets the step resolution from full-step up to 1/256-step.

If interpolation is turned on, the Silent Stepper Brick will always interpolate
your step inputs as 1/256-step. If you use full-step mode with interpolation, each
step will generate 256 1/256 steps.

For maximum torque use full-step without interpolation. For maximum resolution use
1/256-step. Turn interpolation on to make the Stepper driving less noisy.

If you often change the speed with high acceleration you should turn the
interpolation off.
"""
function set_step_configuration(device::BrickSilentStepper, step_resolution, interpolation)
    device.deviceInternal.set_step_configuration(step_resolution, interpolation)
end

export get_step_configuration
"""
    $(SIGNATURES)

Returns the step mode as set by :func:`Set Step Configuration`.
"""
function get_step_configuration(device::BrickSilentStepper)
    return device.deviceInternal.get_step_configuration()
end

export drive_forward
"""
    $(SIGNATURES)

Drives the stepper motor forward until :func:`Drive Backward` or
:func:`Stop` is called. The velocity, acceleration and deacceleration as
set by :func:`Set Max Velocity` and :func:`Set Speed Ramping` will be used.
"""
function drive_forward(device::BrickSilentStepper)
    device.deviceInternal.drive_forward()
end

export drive_backward
"""
    $(SIGNATURES)

Drives the stepper motor backward until :func:`Drive Forward` or
:func:`Stop` is triggered. The velocity, acceleration and deacceleration as
set by :func:`Set Max Velocity` and :func:`Set Speed Ramping` will be used.
"""
function drive_backward(device::BrickSilentStepper)
    device.deviceInternal.drive_backward()
end

export stop
"""
    $(SIGNATURES)

Stops the stepper motor with the deacceleration as set by
:func:`Set Speed Ramping`.
"""
function stop(device::BrickSilentStepper)
    device.deviceInternal.stop()
end

export get_stack_input_voltage
"""
    $(SIGNATURES)

Returns the stack input voltage. The stack input voltage is the
voltage that is supplied via the stack, i.e. it is given by a
Step-Down or Step-Up Power Supply.
"""
function get_stack_input_voltage(device::BrickSilentStepper)
    return device.deviceInternal.get_stack_input_voltage()
end

export get_external_input_voltage
"""
    $(SIGNATURES)

Returns the external input voltage. The external input voltage is
given via the black power input connector on the Silent Stepper Brick.

If there is an external input voltage and a stack input voltage, the motor
will be driven by the external input voltage. If there is only a stack
voltage present, the motor will be driven by this voltage.

.. warning::
 This means, if you have a high stack voltage and a low external voltage,
 the motor will be driven with the low external voltage. If you then remove
 the external connection, it will immediately be driven by the high
 stack voltage
"""
function get_external_input_voltage(device::BrickSilentStepper)
    return device.deviceInternal.get_external_input_voltage()
end

export set_motor_current
"""
    $(SIGNATURES)

Sets the current with which the motor will be driven.

.. warning::
 Do not set this value above the specifications of your stepper motor.
 Otherwise it may damage your motor.
"""
function set_motor_current(device::BrickSilentStepper, current)
    device.deviceInternal.set_motor_current(current)
end

export get_motor_current
"""
    $(SIGNATURES)

Returns the current as set by :func:`Set Motor Current`.
"""
function get_motor_current(device::BrickSilentStepper)
    return device.deviceInternal.get_motor_current()
end

export enable
"""
    $(SIGNATURES)

Enables the driver chip. The driver parameters can be configured (maximum velocity,
acceleration, etc) before it is enabled.
"""
function enable(device::BrickSilentStepper)
    device.deviceInternal.enable()
end

export disable
"""
    $(SIGNATURES)

Disables the driver chip. The configurations are kept (maximum velocity,
acceleration, etc) but the motor is not driven until it is enabled again.

.. warning::
 Disabling the driver chip while the motor is still turning can damage the
 driver chip. The motor should be stopped calling :func:`Stop` function
 before disabling the motor power. The :func:`Stop` function will **not**
 wait until the motor is actually stopped. You have to explicitly wait for the
 appropriate time after calling the :func:`Stop` function before calling
 the :func:`Disable` function.
"""
function disable(device::BrickSilentStepper)
    device.deviceInternal.disable()
end

export is_enabled
"""
    $(SIGNATURES)

Returns *true* if the driver chip is enabled, *false* otherwise.
"""
function is_enabled(device::BrickSilentStepper)
    return device.deviceInternal.is_enabled()
end

export set_basic_configuration
"""
    $(SIGNATURES)

Sets the basic configuration parameters for the different modes (Stealth, Coolstep, Classic).

* Standstill Current: This value can be used to lower the current during stand still. This might
  be reasonable to reduce the heating of the motor and the Brick. When the motor is in standstill
  the configured motor phase current will be driven until the configured
  Power Down Time is elapsed. After that the phase current will be reduced to the standstill
  current. The elapsed time for this reduction can be configured with the Standstill Delay Time.
  The maximum allowed value is the configured maximum motor current
  (see :func:`Set Motor Current`).

* Motor Run Current: The value sets the motor current when the motor is running.
  Use a value of at least one half of the global maximum motor current for a good
  microstep performance. The maximum allowed value is the current
  motor current. The API maps the entered value to 1/32 ... 32/32 of the maximum
  motor current. This value should be used to change the motor current during motor movement,
  whereas the global maximum motor current should not be changed while the motor is moving
  (see :func:`Set Motor Current`).

* Standstill Delay Time: Controls the duration for motor power down after a motion
  as soon as standstill is detected and the Power Down Time is expired. A high Standstill Delay
  Time results in a smooth transition that avoids motor jerk during power down.

* Power Down Time: Sets the delay time after a stand still.

* Stealth Threshold: Sets the upper threshold for Stealth mode.
  If the velocity of the motor goes above this value, Stealth mode is turned
  off. Otherwise it is turned on. In Stealth mode the torque declines with high speed.

* Coolstep Threshold: Sets the lower threshold for Coolstep mode.
  The Coolstep Threshold needs to be above the Stealth Threshold.

* Classic Threshold: Sets the lower threshold for classic mode.
  In classic mode the stepper becomes more noisy, but the torque is maximized.

* High Velocity Chopper Mode: If High Velocity Chopper Mode is enabled, the stepper control
  is optimized to run the stepper motors at high velocities.

If you want to use all three thresholds make sure that
Stealth Threshold < Coolstep Threshold < Classic Threshold.
"""
function set_basic_configuration(device::BrickSilentStepper, standstill_current, motor_run_current, standstill_delay_time, power_down_time, stealth_threshold, coolstep_threshold, classic_threshold, high_velocity_chopper_mode)
    device.deviceInternal.set_basic_configuration(standstill_current, motor_run_current, standstill_delay_time, power_down_time, stealth_threshold, coolstep_threshold, classic_threshold, high_velocity_chopper_mode)
end

export get_basic_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Basic Configuration`.
"""
function get_basic_configuration(device::BrickSilentStepper)
    return device.deviceInternal.get_basic_configuration()
end

export set_spreadcycle_configuration
"""
    $(SIGNATURES)

Note: If you don't know what any of this means you can very likely keep all of
the values as default!

Sets the Spreadcycle configuration parameters. Spreadcycle is a chopper algorithm which actively
controls the motor current flow. More information can be found in the TMC2130 datasheet on page
47 (7 spreadCycle and Classic Chopper).

* Slow Decay Duration: Controls duration of off time setting of slow decay phase.
  0 = driver disabled, all bridges off. Use 1 only with Comparator Blank time >= 2.

* Enable Random Slow Decay: Set to false to fix chopper off time as set by Slow Decay Duration.
  If you set it to true, Decay Duration is randomly modulated.

* Fast Decay Duration: Sets the fast decay duration. This parameters is
  only used if the Chopper Mode is set to Fast Decay.

* Hysteresis Start Value: Sets the hysteresis start value. This parameter is
  only used if the Chopper Mode is set to Spread Cycle.

* Hysteresis End Value: Sets the hysteresis end value. This parameter is
  only used if the Chopper Mode is set to Spread Cycle.

* Sine Wave Offset: Sets the sine wave offset. This parameters is
  only used if the Chopper Mode is set to Fast Decay. 1/512 of the value becomes added to the absolute
  value of the sine wave.

* Chopper Mode: 0 = Spread Cycle, 1 = Fast Decay.

* Comparator Blank Time: Sets the blank time of the comparator. Available values are

  * 0 = 16 clocks,
  * 1 = 24 clocks,
  * 2 = 36 clocks and
  * 3 = 54 clocks.

  A value of 1 or 2 is recommended for most applications.

* Fast Decay Without Comparator: If set to true the current comparator usage for termination of the
  fast decay cycle is disabled.
"""
function set_spreadcycle_configuration(device::BrickSilentStepper, slow_decay_duration, enable_random_slow_decay, fast_decay_duration, hysteresis_start_value, hysteresis_end_value, sine_wave_offset, chopper_mode, comparator_blank_time, fast_decay_without_comparator)
    device.deviceInternal.set_spreadcycle_configuration(slow_decay_duration, enable_random_slow_decay, fast_decay_duration, hysteresis_start_value, hysteresis_end_value, sine_wave_offset, chopper_mode, comparator_blank_time, fast_decay_without_comparator)
end

export get_spreadcycle_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Basic Configuration`.
"""
function get_spreadcycle_configuration(device::BrickSilentStepper)
    return device.deviceInternal.get_spreadcycle_configuration()
end

export set_stealth_configuration
"""
    $(SIGNATURES)

Note: If you don't know what any of this means you can very likely keep all of
the values as default!

Sets the configuration relevant for Stealth mode.

* Enable Stealth: If set to true the stealth mode is enabled, if set to false the
  stealth mode is disabled, even if the speed is below the threshold set in :func:`Set Basic Configuration`.

* Amplitude: If autoscale is disabled, the PWM amplitude is scaled by this value. If autoscale is enabled,
  this value defines the maximum PWM amplitude change per half wave.

* Gradient: If autoscale is disabled, the PWM gradient is scaled by this value. If autoscale is enabled,
  this value defines the maximum PWM gradient. With autoscale a value above 64 is recommended,
  otherwise the regulation might not be able to measure the current.

* Enable Autoscale: If set to true, automatic current control is used. Otherwise the user defined
  amplitude and gradient are used.

* Force Symmetric: If true, A symmetric PWM cycle is enforced. Otherwise the PWM value may change within each
  PWM cycle.

* Freewheel Mode: The freewheel mode defines the behavior in stand still if the Standstill Current
  (see :func:`Set Basic Configuration`) is set to 0.
"""
function set_stealth_configuration(device::BrickSilentStepper, enable_stealth, amplitude, gradient, enable_autoscale, force_symmetric, freewheel_mode)
    device.deviceInternal.set_stealth_configuration(enable_stealth, amplitude, gradient, enable_autoscale, force_symmetric, freewheel_mode)
end

export get_stealth_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Stealth Configuration`.
"""
function get_stealth_configuration(device::BrickSilentStepper)
    return device.deviceInternal.get_stealth_configuration()
end

export set_coolstep_configuration
"""
    $(SIGNATURES)

Note: If you don't know what any of this means you can very likely keep all of
the values as default!

Sets the configuration relevant for Coolstep.

* Minimum Stallguard Value: If the Stallguard result falls below this value*32, the motor current
  is increased to reduce motor load angle. A value of 0 turns Coolstep off.

* Maximum Stallguard Value: If the Stallguard result goes above
  (Min Stallguard Value + Max Stallguard Value + 1) * 32, the motor current is decreased to save
  energy.

* Current Up Step Width: Sets the up step increment per Stallguard value. The value range is 0-3,
  corresponding to the increments 1, 2, 4 and 8.

* Current Down Step Width: Sets the down step decrement per Stallguard value. The value range is 0-3,
  corresponding to the decrements 1, 2, 8 and 16.

* Minimum Current: Sets the minimum current for Coolstep current control. You can choose between
  half and quarter of the run current.

* Stallguard Threshold Value: Sets the level for stall output (see :func:`Get Driver Status`).
  A lower value gives a higher sensitivity. You have to find a suitable value for your
  motor by trial and error, 0 works for most motors.

* Stallguard Mode: Set to 0 for standard resolution or 1 for filtered mode. In filtered mode the Stallguard
  signal will be updated every four full-steps.
"""
function set_coolstep_configuration(device::BrickSilentStepper, minimum_stallguard_value, maximum_stallguard_value, current_up_step_width, current_down_step_width, minimum_current, stallguard_threshold_value, stallguard_mode)
    device.deviceInternal.set_coolstep_configuration(minimum_stallguard_value, maximum_stallguard_value, current_up_step_width, current_down_step_width, minimum_current, stallguard_threshold_value, stallguard_mode)
end

export get_coolstep_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Coolstep Configuration`.
"""
function get_coolstep_configuration(device::BrickSilentStepper)
    return device.deviceInternal.get_coolstep_configuration()
end

export set_misc_configuration
"""
    $(SIGNATURES)

Note: If you don't know what any of this means you can very likely keep all of
the values as default!

Sets miscellaneous configuration parameters.

* Disable Short To Ground Protection: Set to false to enable short to ground protection, otherwise
  it is disabled.

* Synchronize Phase Frequency: With this parameter you can synchronize the chopper for both phases
  of a two phase motor to avoid the occurrence of a beat. The value range is 0-15. If set to 0,
  the synchronization is turned off. Otherwise the synchronization is done through the formula
  f_sync = f_clk/(value*64). In Classic Mode the synchronization is automatically switched off.
  f_clk is 12.8MHz.
"""
function set_misc_configuration(device::BrickSilentStepper, disable_short_to_ground_protection, synchronize_phase_frequency)
    device.deviceInternal.set_misc_configuration(disable_short_to_ground_protection, synchronize_phase_frequency)
end

export get_misc_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Misc Configuration`.
"""
function get_misc_configuration(device::BrickSilentStepper)
    return device.deviceInternal.get_misc_configuration()
end

export get_driver_status
"""
    $(SIGNATURES)

Returns the current driver status.

* Open Load: Indicates if an open load is present on phase A, B or both. This could mean that there is a problem
  with the wiring of the motor. False detection can occur in fast motion as well as during stand still.

* Short To Ground: Indicates if a short to ground is present on phase A, B or both. If this is detected the driver
  automatically becomes disabled and stays disabled until it is enabled again manually.

* Over Temperature: The over temperature indicator switches to "Warning" if the driver IC warms up. The warning flag
  is expected during long duration stepper uses. If the temperature limit is reached the indicator switches
  to "Limit". In this case the driver becomes disabled until it cools down again.

* Motor Stalled: Is true if a motor stall was detected.

* Actual Motor Current: Indicates the actual current control scaling as used in Coolstep mode.
  It represents a multiplier of 1/32 to 32/32 of the
  ``Motor Run Current`` as set by :func:`Set Basic Configuration`. Example: If a ``Motor Run Current``
  of 1000mA was set and the returned value is 15, the ``Actual Motor Current`` is 16/32*1000mA = 500mA.

* Stallguard Result: Indicates the load of the motor. A lower value signals a higher load. Per trial and error
  you can find out which value corresponds to a suitable torque for the velocity used in your application.
  After that you can use this threshold value to find out if a motor stall becomes probable and react on it (e.g.
  decrease velocity).
  During stand still this value can not be used for stall detection, it shows the chopper on-time for motor coil A.

* Stealth Voltage Amplitude: Shows the actual PWM scaling. In Stealth mode it can be used to detect motor load and
  stall if autoscale is enabled (see :func:`Set Stealth Configuration`).
"""
function get_driver_status(device::BrickSilentStepper)
    return device.deviceInternal.get_driver_status()
end

export set_minimum_voltage
"""
    $(SIGNATURES)

Sets the minimum voltage, below which the :cb:`Under Voltage` callback
is triggered. The minimum possible value that works with the Silent Stepper
Brick is 8V.
You can use this function to detect the discharge of a battery that is used
to drive the stepper motor. If you have a fixed power supply, you likely do
not need this functionality.
"""
function set_minimum_voltage(device::BrickSilentStepper, voltage)
    device.deviceInternal.set_minimum_voltage(voltage)
end

export get_minimum_voltage
"""
    $(SIGNATURES)

Returns the minimum voltage as set by :func:`Set Minimum Voltage`.
"""
function get_minimum_voltage(device::BrickSilentStepper)
    return device.deviceInternal.get_minimum_voltage()
end

export set_time_base
"""
    $(SIGNATURES)

Sets the time base of the velocity and the acceleration of the Silent Stepper
Brick.

For example, if you want to make one step every 1.5 seconds, you can set
the time base to 15 and the velocity to 10. Now the velocity is
10steps/15s = 1steps/1.5s.
"""
function set_time_base(device::BrickSilentStepper, time_base)
    device.deviceInternal.set_time_base(time_base)
end

export get_time_base
"""
    $(SIGNATURES)

Returns the time base as set by :func:`Set Time Base`.
"""
function get_time_base(device::BrickSilentStepper)
    return device.deviceInternal.get_time_base()
end

export get_all_data
"""
    $(SIGNATURES)

Returns the following parameters: The current velocity,
the current position, the remaining steps, the stack voltage, the external
voltage and the current consumption of the stepper motor.

The current consumption is calculated by multiplying the ``Actual Motor Current``
value (see :func:`Set Basic Configuration`) with the ``Motor Run Current``
(see :func:`Get Driver Status`). This is an internal calculation of the
driver, not an independent external measurement.

The current consumption calculation was broken up to firmware 2.0.1, it is fixed
since firmware 2.0.2.

There is also a callback for this function, see :cb:`All Data` callback.
"""
function get_all_data(device::BrickSilentStepper)
    return device.deviceInternal.get_all_data()
end

export set_all_data_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`All Data` callback is triggered
periodically. A value of 0 turns the callback off.
"""
function set_all_data_period(device::BrickSilentStepper, period)
    device.deviceInternal.set_all_data_period(period)
end

export get_all_data_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set All Data Period`.
"""
function get_all_data_period(device::BrickSilentStepper)
    return device.deviceInternal.get_all_data_period()
end

export set_spitfp_baudrate_config
"""
    $(SIGNATURES)

The SPITF protocol can be used with a dynamic baudrate. If the dynamic baudrate is
enabled, the Brick will try to adapt the baudrate for the communication
between Bricks and Bricklets according to the amount of data that is transferred.

The baudrate will be increased exponentially if lots of data is sent/received and
decreased linearly if little data is sent/received.

This lowers the baudrate in applications where little data is transferred (e.g.
a weather station) and increases the robustness. If there is lots of data to transfer
(e.g. Thermal Imaging Bricklet) it automatically increases the baudrate as needed.

In cases where some data has to transferred as fast as possible every few seconds
(e.g. RS485 Bricklet with a high baudrate but small payload) you may want to turn
the dynamic baudrate off to get the highest possible performance.

The maximum value of the baudrate can be set per port with the function
:func:`Set SPITFP Baudrate`. If the dynamic baudrate is disabled, the baudrate
as set by :func:`Set SPITFP Baudrate` will be used statically.

.. versionadded:: 2.0.4\$nbsp;(Firmware)
"""
function set_spitfp_baudrate_config(device::BrickSilentStepper, enable_dynamic_baudrate, minimum_dynamic_baudrate)
    device.deviceInternal.set_spitfp_baudrate_config(enable_dynamic_baudrate, minimum_dynamic_baudrate)
end

export get_spitfp_baudrate_config
"""
    $(SIGNATURES)

Returns the baudrate config, see :func:`Set SPITFP Baudrate Config`.

.. versionadded:: 2.0.4\$nbsp;(Firmware)
"""
function get_spitfp_baudrate_config(device::BrickSilentStepper)
    return device.deviceInternal.get_spitfp_baudrate_config()
end

export get_send_timeout_count
"""
    $(SIGNATURES)

Returns the timeout count for the different communication methods.

The methods 0-2 are available for all Bricks, 3-7 only for Master Bricks.

This function is mostly used for debugging during development, in normal operation
the counters should nearly always stay at 0.
"""
function get_send_timeout_count(device::BrickSilentStepper, communication_method)
    return device.deviceInternal.get_send_timeout_count(communication_method)
end

export set_spitfp_baudrate
"""
    $(SIGNATURES)

Sets the baudrate for a specific Bricklet port.

If you want to increase the throughput of Bricklets you can increase
the baudrate. If you get a high error count because of high
interference (see :func:`Get SPITFP Error Count`) you can decrease the
baudrate.

If the dynamic baudrate feature is enabled, the baudrate set by this
function corresponds to the maximum baudrate (see :func:`Set SPITFP Baudrate Config`).

Regulatory testing is done with the default baudrate. If CE compatibility
or similar is necessary in your applications we recommend to not change
the baudrate.
"""
function set_spitfp_baudrate(device::BrickSilentStepper, bricklet_port, baudrate)
    device.deviceInternal.set_spitfp_baudrate(bricklet_port, baudrate)
end

export get_spitfp_baudrate
"""
    $(SIGNATURES)

Returns the baudrate for a given Bricklet port, see :func:`Set SPITFP Baudrate`.
"""
function get_spitfp_baudrate(device::BrickSilentStepper, bricklet_port)
    return device.deviceInternal.get_spitfp_baudrate(bricklet_port)
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Brick side. All
Bricklets have a similar function that returns the errors on the Bricklet side.
"""
function get_spitfp_error_count(device::BrickSilentStepper, bricklet_port)
    return device.deviceInternal.get_spitfp_error_count(bricklet_port)
end

export enable_status_led
"""
    $(SIGNATURES)

Enables the status LED.

The status LED is the blue LED next to the USB connector. If enabled is is
on and it flickers if data is transfered. If disabled it is always off.

The default state is enabled.
"""
function enable_status_led(device::BrickSilentStepper)
    device.deviceInternal.enable_status_led()
end

export disable_status_led
"""
    $(SIGNATURES)

Disables the status LED.

The status LED is the blue LED next to the USB connector. If enabled is is
on and it flickers if data is transfered. If disabled it is always off.

The default state is enabled.
"""
function disable_status_led(device::BrickSilentStepper)
    device.deviceInternal.disable_status_led()
end

export is_status_led_enabled
"""
    $(SIGNATURES)

Returns *true* if the status LED is enabled, *false* otherwise.
"""
function is_status_led_enabled(device::BrickSilentStepper)
    return device.deviceInternal.is_status_led_enabled()
end

export get_protocol1_bricklet_name
"""
    $(SIGNATURES)

Returns the firmware and protocol version and the name of the Bricklet for a
given port.

This functions sole purpose is to allow automatic flashing of v1.x.y Bricklet
plugins.
"""
function get_protocol1_bricklet_name(device::BrickSilentStepper, port)
    return device.deviceInternal.get_protocol1_bricklet_name(port)
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has an
accuracy of ±15%. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickSilentStepper)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Brick. Calling this function
on a Brick inside of a stack will reset the whole stack.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickSilentStepper)
    device.deviceInternal.reset()
end

export write_bricklet_plugin
"""
    $(SIGNATURES)

Writes 32 bytes of firmware to the bricklet attached at the given port.
The bytes are written to the position offset * 32.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_bricklet_plugin(device::BrickSilentStepper, port, offset, chunk)
    device.deviceInternal.write_bricklet_plugin(port, offset, chunk)
end

export read_bricklet_plugin
"""
    $(SIGNATURES)

Reads 32 bytes of firmware from the bricklet attached at the given port.
The bytes are read starting at the position offset * 32.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function read_bricklet_plugin(device::BrickSilentStepper, port, offset)
    return device.deviceInternal.read_bricklet_plugin(port, offset)
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Brick is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position is the position in the stack from '0' (bottom) to '8' (top).

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickSilentStepper)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickSilentStepper, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
#trying out VCFTools.jl and GeneticVariation.jl
using Revise
using GeneticVariation
using Random
using VCFTools
using BenchmarkTools

cd("/Users/biona001/.julia/dev/VCFTools/test")
vcffile = "test.08Jun17.d8b.vcf"
reader = VCF.Reader(openvcf(vcffile, "r"))
record = read(reader)

record = VCF.Record("20\t14370\trs6054257\tG\tA\t29\tPASS\tNS=3;DP=14;AF=0.5;DB;H2\tGT:GQ:DP:HQ\t0|0:48:1:51,51\t1|0:48:8:51,51")
record = VCF.Record("20\t14370\t.\tG\tA\t29\tPASS;fdsa\t.\tDS\t1.99\t1.01")
record = VCF.Record("20\t14370\t.\tG\tA\t29\tPASS\t.\tGT\t1/0\t0/0")
record.genotype

#test convert_ht
using Revise
using GeneticVariation
using VCFTools

cd("/Users/biona001/.julia/dev/VCFTools/test")
vcffile = "test.08Jun17.d8b.vcf"
H = convert_ht(Float32, vcffile)

# byte representation mapping
String([0x09]) # '\t'
String([0x3a]) # ':'
String([0x30]) # '0'
String([0x31]) # '1'
String([0x2e]) # '.'
String([0x7c]) # '|'
String([0x2f]) # '/'

# test filter
using Revise
using GeneticVariation
using Random
using VCFTools
using BenchmarkTools

#import data
cd("/Users/biona001/.julia/dev/VCFTools/test")
vcffile = "test.08Jun17.d8b.vcf.gz"
samples = nsamples(vcffile)
records = nrecords(vcffile)
reader = VCF.Reader(openvcf(vcffile, "r"))
record = read(reader)

# define masks
Random.seed!(123)
record_mask = trues(records)
# sample_mask = bitrand(samples)
sample_mask = trues(samples)
sample_mask[1] = false
sample_mask[end] = false

des = "filtered.test.08Jun17.d8b.vcf.gz"
VCFTools.filter(vcffile, record_mask, sample_mask, des=des)
reader = VCF.Reader(openvcf(des, "r"))
record = read(reader)
record.genotype

nsamples(des)

X = convert_gt(Float32, vcffile)
X_filter = convert_gt(Float32, "filtered." * vcffile)

@benchmark VCFTools.filter(vcffile, record_mask, sample_mask)
#new: 149.763 ms, 214.48 MiB
#old: 134.655 ms, 177.42 MiB

des = "filter." * vcffile
reader = VCF.Reader(openvcf(vcffile, "r"))
writer = VCF.Writer(openvcf(des, "w"), filter_header(reader, sample_mask))

reader = VCF.Reader(openvcf(vcffile, "r"))
record = read(reader)

# test masking

using Revise
using GeneticVariation
using Random
using VCFTools
using BenchmarkTools

cd("/Users/biona001/.julia/dev/VCFTools/test")
vcffile = "test.08Jun17.d8b.vcf"
samples = nsamples(vcffile)
records = nrecords(vcffile)
masks = bitrand(records, samples)
mask_gt(vcffile, masks)
@benchmark mask_gt(vcffile, masks) #82.221 ms, 108.69 MiB, 800187 alloc

# test convert_ds
using Revise
using GeneticVariation
using Random
using VCFTools
using BenchmarkTools

cd("/Users/biona001/Benjamin_Folder/UCLA/research/2nd_project/benchmarks/AFRPed/minimac4")
vcffile = "minimac4_result.dose.vcf"
dosage = convert_ds(Float64, vcffile)

# test overwrite_alt_ref_allele
using Revise
using GeneticVariation
using Random
using VCFTools
using BenchmarkTools

function overwrite_alt_ref_allele(
    src::AbstractString,
    des::AbstractString,
    ref::String = "A",
    alt::String = "T"
    )

    # create input and output VCF files
    reader = VCF.Reader(openvcf(src, "r"))
    writer = VCF.Writer(openvcf(des, "w"), VCF.header(reader))

    # define ref/alt allele
    dict = Dict([("A", 0x41), ("T", 0x54), ("C", 0x43), ("G", 0x47), ("N", 0x4e)])
    haskey(dict, ref) || error("ref allele must be A, T, C, G, or N.")
    haskey(dict, alt) || error("ref allele must be A, T, C, G, or N.")
    ref, alt = dict[ref], dict[alt]
    for record in reader
        record.data[record.alt[1]] .= alt
        record.data[record.ref[1]] = ref
        write(writer, record)
    end
    flush(writer); close(reader); close(writer)
end

src = "full_test.vcf"
reader = VCF.Reader(openvcf(src, "r"))
record = read(reader)
# overwrite_alt_ref_allele("test.vcf")
overwrite_alt_ref_allele("test.vcf", "test_result.vcf")

# Simulate genotype matrix using real haplotypes
# 1000 Genomes Phase 1, chr 22 data downloaded from: https://genome.sph.umich.edu/wiki/Minimac4
using Revise
using MendelImpute
using VCFTools
using BenchmarkTools

cd("/Users/biona001/.julia/dev/MendelImpute/simulation/compare4/VCF_Files")

#2184 haplotypes each with 365644 SNPs -> 99MB of RAM
@time H = convert_ht(Bool, "ALL.chr22.phase1_v3.snps_indels_svs.genotypes.all.noSingleton.vcf.gz")

#simulate 100 samples each with 365644 SNPs
X = simulate_genotypes(H', 100)


text:
mutable struct SimulationParameters
  endtime::Float64
  filenamestub::String
  time::Float64
  iteration::Int
end

SimulationParameters(endtime, filenamestub) = SimulationParameters(endtime, filenamestub, 0.0, 0)

@concrete struct Simulation{F<:AbstractField}
  plasma::Plasma
  field::F
  timeintegrator
  diagnosticdumpevery::Int64
  endtime::Float64
  filenamestub::String
  time::Ref{Float64}
  iteration::Ref{Int}
  diagnosticdumpcounter::Ref{Int}
end

function Simulation(plasma::Plasma, field::AbstractField, timeintegrator::AbstractTimeIntegrator;
  diagnosticdumpevery=1, endtime=10.0, filenamestub=nothing)
  isnothing(filenamestub) && (filenamestub = "$(typeof(eltype(plasma[1])))-$(typeof(field))-$(timeintegrator)-")
  return Simulation(plasma, field, timeintegrator, diagnosticdumpevery, endtime, filenamestub, Ref(0.0), Ref(0), Ref(-1))
end
filenamestub(sim::Simulation) = sim.filenamestub
filename(filepathstub::String, i::Int) = filepathstub * lpad("$i", 6, "0") * ".jld2"
filename(sim::Simulation) = filename(sim.filenamestub, diagnosticdumpcounter(sim))
timestep(sim::Simulation) = timestep(sim.timeintegrator)
iteration(sim::Simulation) = sim.iteration[]
diagnosticdumpcounter(sim::Simulation) = sim.diagnosticdumpcounter[]
time(sim::Simulation) = sim.time[]
iterate!(sim::Simulation) = (sim.iteration[] += 1)
function save(sim::Simulation)
  sim.diagnosticdumpcounter[] += 1
  @save "$(filename(sim))" sim
end
function load(filepathstub::String, diagnosticdumpcount::Int)
  @load "$(filename(filepathstub, diagnosticdumpcount))" sim
  return sim
end

function init!(sim::Simulation)
  sim.diagnosticdumpcounter[] = -1
  zero!(sim.field)
  deposit!(sim.field, sim.plasma)
  solve!(sim.field)
  return sim
end
function run!(sim::Simulation)
  save(sim) # save initial conditions as file 0
  while time(sim) <= sim.endtime
    runtimestep!(sim)
  end
end

function runtimestep!(sim::Simulation)
  integrate!(sim.plasma, sim.field, sim.timeintegrator)
  sim.iteration[] += 1
  sim.time[] += timestep(sim.timeintegrator)
  if mod(iteration(sim), sim.diagnosticdumpevery) == 0
    save(sim)
  end
  return sim
end


text:
"""
    abstract type Transformable
This is a bit of a weird name, but all scenes and plots are transformable,
so that's what they all have in common. This might be better expressed as traits.
"""
abstract type Transformable end

abstract type AbstractPlot{Typ} <: Transformable end
abstract type AbstractScene <: Transformable end
abstract type ScenePlot{Typ} <: AbstractPlot{Typ} end

const SceneLike = Union{AbstractScene, ScenePlot}
# const Attributes = Dict{Symbol, Any}

abstract type AbstractCamera end

# placeholder if no camera is present
struct EmptyCamera <: AbstractCamera end

@enum RaymarchAlgorithm begin
    IsoValue
    Absorption
    MaximumIntensityProjection
    AbsorptionRGBA
    IndexedAbsorptionRGBA
end

const RealVector{T} = AbstractVector{T} where T <: Number

const Node = Observable

const Rect{N, T} = HyperRectangle{N, T}
const Rect2D{T} = HyperRectangle{2, T}
const FRect2D = Rect2D{Float32}

"A generic, three dimensional rectangle."
const Rect3D{T} = Rect{3, T}

"An float valued, three dimensional rectangle."
const FRect3D = Rect3D{Float32}

"An integer valued, three dimensional rectangle."
const IRect3D = Rect3D{Int}

"An integer valued, two dimensional rectangle."
const IRect2D = Rect2D{Int}

const Point2d{T} = NTuple{2, T}
const Vec2d{T} = NTuple{2, T}
const VecTypes{N, T} = Union{StaticVector{N, T}, NTuple{N, T}}
const NVec{N} = Union{StaticVector{N}, NTuple{N, Any}}
const RGBAf0 = RGBA{Float32}
const RGBf0 = RGB{Float32}

abstract type AbstractScreen <: AbstractDisplay end

"""
    IRect(x, y, w, h)

Creates a two dimensional rectangle of integer dimensions, at origin (x, y)
and with width w and height h
"""
function IRect(x, y, w, h)
    HyperRectangle{2, Int}(Vec(round(Int, x), round(Int, y)), Vec(round(Int, w), round(Int, h)))
end

"""
    IRect(xy::VecTypes, w, h)

Creates a two dimensional rectangle of integer dimensions, with origin
at vector xy, and with width w and height h
"""
function IRect(xy::VecTypes, w, h)
    IRect(xy[1], xy[2], w, h)
end

"""
    IRect(x, y, wh::VecTypes)

Creates a two dimensional rectangle of integer dimensions, with origin
at (x, y), and with width and height as the respective components of vector wh
"""
function IRect(x, y, wh::VecTypes)
    IRect(x, y, wh[1], wh[2])
end

"""
    IRect(xy::VecTypes, wh::VecTypes)

Creates a two dimensional rectangle of integer dimensions, with origin
at vector xy, and with width and height as the respective components of vector wh
"""
function IRect(xy::VecTypes, wh::VecTypes)
    IRect(xy[1], xy[2], wh[1], wh[2])
end

"""
    IRect(xy::NamedTuple{(:x, :y)}, wh::NamedTuple{(:width, :height)})

This takes two named tuples and constructs an integer valued rectangle with them.
"""
function IRect(xy::NamedTuple{(:x, :y)}, wh::NamedTuple{(:width, :height)})
    IRect(xy.x, xy.y, wh.width, wh.height)
end

function positive_widths(rect::HyperRectangle{N, T}) where {N, T}
    mini, maxi = minimum(rect), maximum(rect)
    realmin = min.(mini, maxi)
    realmax = max.(mini, maxi)
    HyperRectangle{N, T}(realmin, realmax .- realmin)
end

"""
    FRect(x, y, w, h)

Creates a two dimensional rectangle, at origin (x, y)
and with width w and height h.  Formally defined as the
Cartesian product of the intervals (x, y) and (w, h).
"""
function FRect(x, y, w, h)
    HyperRectangle{2, Float32}(Vec2f0(x, y), Vec2f0(w, h))
end
function FRect(r::SimpleRectangle)
    FRect(r.x, r.y, r.w, r.h)
end
function FRect(r::Rect)
    FRect(minimum(r), widths(r))
end
function FRect(xy::VecTypes, w, h)
    FRect(xy[1], xy[2], w, h)
end
function FRect(x, y, wh::VecTypes)
    FRect(x, y, wh[1], wh[2])
end
function FRect(xy::VecTypes, wh::VecTypes)
    FRect(xy[1], xy[2], wh[1], wh[2])
end

function FRect3D(x::Tuple{Tuple{<: Number, <: Number}, Tuple{<: Number, <: Number}})
    FRect3D(Vec3f0(x[1]..., 0), Vec3f0(x[2]..., 0))
end
function FRect3D(x::Tuple{Tuple{<: Number, <: Number, <: Number}, Tuple{<: Number, <: Number, <: Number}})
    FRect3D(Vec3f0(x[1]...), Vec3f0(x[2]...))
end

function FRect3D(x::Rect2D)
    FRect3D(Vec3f0(minimum(x)..., 0), Vec3f0(widths(x)..., 0.0))
end
# For now, we use Reactive.Signal as our Node type. This might change in the future
const Node = Observable

include("interaction/iodevices.jl")

"""
This struct provides accessible `Observable`s to monitor the events
associated with a Scene.

## Fields
$(TYPEDFIELDS)
"""
struct Events
    """
    The area of the window in pixels, as an [`IRect2D`](@ref).
    """
    window_area::Node{IRect2D}
    """
    The DPI resolution of the window, as a `Float64`.
    """
    window_dpi::Node{Float64}
    """
    The state of the window (open => true, closed => false).
    """
    window_open::Node{Bool}

    """
    The pressed mouse buttons.
    Updates when a mouse button is pressed.

    See also [`ispressed`](@ref).
    """
    mousebuttons::Node{Set{Mouse.Button}}
    """
    The position of the mouse as a [`Point2`](@ref).
    Updates whenever the mouse moves.
    """
    mouseposition::Node{Point2d{Float64}}
    """
    The
    """
    mousedrag::Node{Mouse.DragEnum}
    """
    The direction of scroll
    """
    scroll::Node{Vec2d{Float64}}

    """
    See also [`ispressed`](@ref).
    """
    keyboardbuttons::Node{Set{Keyboard.Button}}

    unicode_input::Node{Vector{Char}}
    dropped_files::Node{Vector{String}}
    """
    Whether the Scene window is in focus or not.
    """
    hasfocus::Node{Bool}
    entered_window::Node{Bool}
end

function Events()
    Events(
        node(:window_area, IRect(0, 0, 0, 0)),
        node(:window_dpi, 100.0),
        node(:window_open, false),

        node(:mousebuttons, Set{Mouse.Button}()),
        node(:mouseposition, (0.0, 0.0)),
        node(:mousedrag, Mouse.notpressed),
        node(:scroll, (0.0, 0.0)),

        node(:keyboardbuttons, Set{Keyboard.Button}()),

        node(:unicode_input, Char[]),
        node(:dropped_files, String[]),
        node(:hasfocus, false),
        node(:entered_window, false),
    )
end

"""
"""
mutable struct Camera
    view::Node{Mat4f0}
    projection::Node{Mat4f0}
    projectionview::Node{Mat4f0}
    resolution::Node{Vec2f0}
    eyeposition::Node{Vec3f0}
    steering_nodes::Vector{Any}
end

"""
Holds the transformations for Scenes.

## Fields
$(TYPEDFIELDS)
"""
struct Transformation <: Transformable
    parent::RefValue{Transformable}
    translation::Node{Vec3f0}
    scale::Node{Vec3f0}
    rotation::Node{Quaternionf0}
    model::Node{Mat4f0}
    flip::Node{NTuple{3, Bool}}
    align::Node{Vec2f0}
    data_func::Node{Any}
    function Transformation(translation, scale, rotation, model, flip, align, data_func)
        return new(
            RefValue{Transformable}(),
            translation, scale, rotation, model, flip, align, data_func
        )
    end
end

struct Attributes
    attributes::Dict{Symbol, Node}
end
Base.broadcastable(x::AbstractScene) = Ref(x)
Base.broadcastable(x::AbstractPlot) = Ref(x)
Base.broadcastable(x::Attributes) = Ref(x)

value_convert(x::Observables.AbstractObservable) = Observables.observe(x)
value_convert(@nospecialize(x)) = x

function value_convert(x::NTuple{N, Union{Any, Observables.AbstractObservable}}) where N
    # with an observable at first place,
    # lift correctly dispatches to create a new Node that updates on any
    lift(Node(nothing), x...) do _, args...
        args
    end
end

value_convert(x::NamedTuple) = Attributes(x)

node_pairs(pair::Union{Pair, Tuple{Any, Any}}) = (pair[1] => to_node(Any, value_convert(pair[2]), pair[1]))
node_pairs(pairs) = (node_pairs(pair) for pair in pairs)
Base.convert(::Type{<: Node}, x) = to_node(Any, x)
Base.convert(::Type{T}, x::T) where T <: Node = x
Base.convert(::Type{Node{T}}, x::Node) where T = to_node(T, x)

Attributes(; kw_args...) = Attributes(Dict{Symbol, Node}(node_pairs(kw_args)))
Attributes(pairs::Pair...) = Attributes(Dict{Symbol, Node}(node_pairs(pairs)))
Attributes(pairs::AbstractVector) = Attributes(Dict{Symbol, Node}(node_pairs.(pairs)))
Attributes(pairs::Iterators.Pairs) = Attributes(collect(pairs))
Attributes(nt::NamedTuple) = Attributes(; nt...)
attributes(x::Attributes) = getfield(x, :attributes)
Base.keys(x::Attributes) = keys(x.attributes)
Base.values(x::Attributes) = values(x.attributes)
function Base.iterate(x::Attributes, state...)
    s = iterate(keys(x), state...)
    s === nothing && return nothing
    return (s[1] => x[s[1]], s[2])
end

function Base.copy(attributes::Attributes)
    result = Attributes()
    for (k, v) in attributes
        # We need to create a new Signal to have a real copy
        result[k] = copy(v)
    end
    return result
end
Base.filter(f, x::Attributes) = Attributes(filter(f, attributes(x)))
Base.empty!(x::Attributes) = (empty!(attributes(x)); x)
Base.length(x::Attributes) = length(attributes(x))

function Base.merge!(target::Attributes, args::Attributes...)
    for elem in args
        merge_attributes!(target, elem)
    end
    return target
end
Base.merge(target::Attributes, args::Attributes...) = merge!(copy(target), args...)

@generated hasfield(x::T, ::Val{key}) where {T, key} = :($(key in fieldnames(T)))

@inline function Base.getproperty(x::T, key::Symbol) where T <: Union{Attributes, Transformable}
    if hasfield(x, Val(key))
        getfield(x, key)
    else
        getindex(x, key)
    end
end
@inline function Base.setproperty!(x::T, key::Symbol, value) where T <: Union{Attributes, Transformable}
    if hasfield(x, Val(key))
        setfield!(x, key, value)
    else
        setindex!(x, value, key)
    end
end

function getindex(x::Attributes, key::Symbol)
    x = attributes(x)[key]
    x[] isa Attributes ? x[] : x
end

function setindex!(x::Attributes, value, key::Symbol)
    if haskey(x, key)
        x.attributes[key][] = value
    else
        x.attributes[key] = to_node(Any, value, key)
    end
end

function setindex!(x::Attributes, value::Node, key::Symbol)
    if haskey(x, key)
        # error("You're trying to update an attribute node with a new node. This is not supported right now.
        # You can do this manually like this:
        # lift(val-> attributes[$key] = val, node::$(typeof(value)))
        # ")
        return x.attributes[key] = to_node(Any, value)
    else
        #TODO make this error. Attributes should be sort of immutable
        return x.attributes[key] = to_node(Any, value)
    end
    return x
end

function Base.show(io::IO,::MIME"text/plain", attr::Attributes)

    d = Dict()
    for p in pairs(attr.attributes)
        d[p.first] = to_value(p.second)
    end
    show(IOContext(io, :limit => false), MIME"text/plain"(), d)

end

Base.show(io::IO, attr::Attributes) = show(io, MIME"text/plain"(), attr)

struct Combined{Typ, T} <: ScenePlot{Typ}
    parent::SceneLike
    transformation::Transformation
    attributes::Attributes
    input_args::Tuple
    converted::Tuple
    plots::Vector{AbstractPlot}
end

theme(x::AbstractPlot) = x.attributes
isvisible(x) = haskey(x, :visible) && to_value(x[:visible])

#dict interface
const AttributeOrPlot = Union{AbstractPlot, Attributes}
Base.pop!(x::AttributeOrPlot, args...) = pop!(x.attributes, args...)
haskey(x::AttributeOrPlot, key) = haskey(x.attributes, key)
delete!(x::AttributeOrPlot, key) = delete!(x.attributes, key)
function get!(f::Function, x::AttributeOrPlot, key::Symbol)
    if haskey(x, key)
        return x[key]
    else
        val = f()
        x[key] = val
        return x[key]
    end
end

get!(x::AttributeOrPlot, key::Symbol, default) = get!(()-> default, x, key)
get(f::Function, x::AttributeOrPlot, key::Symbol) = haskey(x, key) ? x[key] : f()
get(x::AttributeOrPlot, key::Symbol, default) = get(()-> default, x, key)

# This is a bit confusing, since for a plot it returns the attribute from the arguments
# and not a plot for integer indexing. But, we want to treat plots as "atomic"
# so from an interface point of view, one should assume that a plot doesn't contain subplots
# Combined plots break this assumption in some way, but the way to look at it is,
# that the plots contained in a Combined plot are not subplots, but _are_ actually
# the plot itself.
getindex(plot::AbstractPlot, idx::Integer) = plot.converted[idx]
getindex(plot::AbstractPlot, idx::UnitRange{<:Integer}) = plot.converted[idx]
setindex!(plot::AbstractPlot, value, idx::Integer) = (plot.input_args[idx][] = value)
Base.length(plot::AbstractPlot) = length(plot.converted)

function getindex(x::AbstractPlot, key::Symbol)
    argnames = argument_names(typeof(x), length(x.converted))
    idx = findfirst(isequal(key), argnames)
    if idx == nothing
        return x.attributes[key]
    else
        x.converted[idx]
    end
end
function getindex(x::AttributeOrPlot, key::Symbol, key2::Symbol, rest::Symbol...)
    dict = to_value(x[key])
    dict isa Attributes || error("Trying to access $(typeof(dict)) with multiple keys: $key, $key2, $(rest)")
    dict[key2, rest...]
end

function setindex!(x::AttributeOrPlot, value, key::Symbol, key2::Symbol, rest::Symbol...)
    dict = to_value(x[key])
    dict isa Attributes || error("Trying to access $(typeof(dict)) with multiple keys: $key, $key2, $(rest)")
    dict[key2, rest...] = value
end

function setindex!(x::AbstractPlot, value, key::Symbol)
    argnames = argument_names(typeof(x), length(x.converted))
    idx = findfirst(isequal(key), argnames)
    if idx == nothing && haskey(x.attributes, key)
        return x.attributes[key][] = value
    elseif !haskey(x.attributes, key)
        x.attributes[key] = to_node(value)
    else
        return setindex!(x.converted[idx], value)
    end
end

function setindex!(x::AbstractPlot, value::Node, key::Symbol)
    argnames = argument_names(typeof(x), length(x.converted))
    idx = findfirst(isequal(key), argnames)
    if idx == nothing
        if haskey(x, key)
            # error("You're trying to update an attribute node with a new node. This is not supported right now.
            # You can do this manually like this:
            # lift(val-> attributes[$key] = val, node::$(typeof(value)))
            # ")
            return x.attributes[key] = value
        else
            return x.attributes[key] = value
        end
    else
        return setindex!(x.converted[idx], value)
    end
end
parent(x::AbstractPlot) = x.parent

"""
Remove `combined` from the current parent, and add it to a new subscene of the
parent scene. Returns the new parent.
"""
function detach!(x::Combined)
    p1 = parent(x)
    filter!(p-> p != x, p1.plots) # remove from parent

    p2 = parent_scene(x) # first scene that parents this node

    sub = Scene(p2, pixelarea(p2)) # subscene
    push!(x.plots, x)

    sub
end

function func2string(func::F) where F <: Function
    string(F.name.mt.name)
end

plotkey(::Type{<: AbstractPlot{Typ}}) where Typ = Symbol(lowercase(func2string(Typ)))
plotkey(::T) where T <: AbstractPlot = plotkey(T)

plotfunc(::Type{<: AbstractPlot{Func}}) where Func = Func
plotfunc(::T) where T <: AbstractPlot = plotfunc(T)
plotfunc(f::Function) = f

func2type(x::T) where T = func2type(T)
func2type(x::Type{<: AbstractPlot}) = x
func2type(f::Function) = Combined{f}

"""
Billboard attribute to always have a primitive face the camera.
Can be used for rotation.
"""
struct Billboard end

"""
Type to indicate that an attribute will get calculated automatically
"""
struct Automatic end

"""
Singleton instance to indicate that an attribute will get calculated automatically
"""
const automatic = Automatic()

const Vecf0{N} = Vec{N, Float32}
const Pointf0{N} = Point{N, Float32}
export Vecf0, Pointf0
const NativeFont = Vector{Ptr{FreeType.FT_FaceRec}}

"""
`PlotSpec{P<:AbstractPlot}(args...; kwargs...)`

Object encoding positional arguments (`args`), a `NamedTuple` of attributes (`kwargs`)
as well as plot type `P` of a basic plot.
"""
struct PlotSpec{P<:AbstractPlot}
    args::Tuple
    kwargs::NamedTuple
    PlotSpec{P}(args...; kwargs...) where {P<:AbstractPlot} = new{P}(args, values(kwargs))
end

PlotSpec(args...; kwargs...) = PlotSpec{Combined{Any}}(args...; kwargs...)

Base.getindex(p::PlotSpec, i::Int) = getindex(p.args, i)
Base.getindex(p::PlotSpec, i::Symbol) = getproperty(p.kwargs, i)

to_plotspec(::Type{P}, args; kwargs...) where {P} =
    PlotSpec{P}(args...; kwargs...)

to_plotspec(::Type{P}, p::PlotSpec{S}; kwargs...) where {P, S} =
    PlotSpec{plottype(P, S)}(p.args...; p.kwargs..., kwargs...)

plottype(::PlotSpec{P}) where {P} = P


text:
@render_program FXAAProgram
@render_program CompositingProgram

struct FinalRenderer <: AbstractRenderSystem end

function Overseer.requested_components(::FinalRenderer)
    return (Canvas, FullscreenVao, FXAAProgram, CompositingProgram, IOTarget)
end

function Overseer.prepare(::FinalRenderer, dio::Diorama)
    e = Entity(dio[DioEntity], 1)
    if isempty(dio[CompositingProgram])
        dio[e] = CompositingProgram(Program(compositing_shaders()))
    end
    if isempty(dio[FXAAProgram])
        dio[e] = FXAAProgram(Program(fxaa_shaders()))
    end
end

function Overseer.update(::FinalRenderer, m::AbstractLedger)
    compositing_program = singleton(m, CompositingProgram)
    fxaa_program        = singleton(m, FXAAProgram)
    canvas              = singleton(m, Canvas)
    vao                 = singleton(m, FullscreenVao)
    iofbo               = singleton(m, IOTarget)
    # bind(iofbo)
    bind(canvas)
    clear!(canvas)
    draw(canvas)

    glDisable(GL_DEPTH_TEST)
    bind(fxaa_program)
    gluniform(fxaa_program, :u_colorTexture, 0, color_attachment(iofbo.target, 1))
    gluniform(fxaa_program, :u_texelStep, Vec2f0((1 ./ size(iofbo))...))

    gluniform(fxaa_program, :u_showEdges, false)
    gluniform(fxaa_program, :u_fxaaOn, true)
    gluniform(fxaa_program, :u_lumaThreshold, 0.3f0)
    gluniform(fxaa_program, :u_mulReduce, 1 / 8.0f0)
    gluniform(fxaa_program, :u_minReduce, 1 / 128.0f0)
    gluniform(fxaa_program, :u_maxSpan, 8.0f0)
    # gluniform(fxaa_program, :frameBufSize, Vec2f0(size(iofbo)...))
    # bind(compositing_program)
    # gluniform(compositing_program, :color_texture, 0, color_attachment(iofbo.target, 1))
    bind(vao)
    draw(vao)
    return glEnable(GL_DEPTH_TEST)
end


text:
"""
# module Geometry

# Examples

```jldoctest
julia>
```
"""
module Geometry

using LinearAlgebra: det

using StaticArrays: FieldVector, SMatrix

export Point,
    Point2D, Point3D,
    Cell,
    volume,
    metric_tensor,
    inverse_metric_tensor,
    metric_tensor_velocity,
    reciprocalcell

abstract type Point{N, Float64} <: FieldVector{N, Float64} end

struct Point2D <: Point{2, Float64}
    x::Float64
    y::Float64
end

struct Point3D <: Point{3, Float64}
    x::Float64
    y::Float64
    z::Float64
end

const Cell = SMatrix{3, 3, Float64}

volume(c::Cell)::Float64 = (abs ∘ det)(c)

metric_tensor(c::Cell) = transpose(c) * c

inverse_metric_tensor(c::Cell) = (inv ∘ metric_tensor)(c)

metric_tensor_velocity(c::Cell, c_derivative::Cell) = transpose(c_derivative) * c + c * transpose(c_derivative)

reciprocalcell(c::Cell) = 2π * (transpose ∘ inv)(c)

end


text:
function discretize_dyns!(
  F,  # coeff. matrix of the sparse tableau eqs.
  nbset,  # tuple with the following elements.
    # nodes: array of Strings containing independent node names.
    # num_nodes: number of nodes.
    # branches: array of struct Branch.
    # num_branches: sum of branch phases.
    # num_dyns: number of dynamic branches.
  h  # time-step size.
  )  #=

discretizes the time derivative operators appearing in the sparse tableau matrix F given as the first argument. The information of branches, the number of nodes, the number of branches and the time step size are respectively given as the second to fifth argument, and they are used in the discretization.

Author: Taku Noda
Started on Aug. 28, 2020
Update on June 29, 2021

Copyright 2021 Central Research Institute of Electric Power Industry

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

=#

( nodes, num_nodes, branches, num_branches, num_dyns ) = nbset

noff_d = num_nodes + 2*num_branches
noff_v = num_nodes + num_branches
h2 = 0.5*h

for branch in branches
  if branch.type == "L"
    nph = branch.nph
    range = 0:( nph - 1 )
    d_range = ( noff_d + branch.p_de ) .+ range
    v_range = ( noff_v + branch.p_vars ) .+ range
    F[ d_range, v_range ] = sparse( -h2*I, nph, nph )
    F[ d_range, d_range ] = sparse( I, nph, nph )
  elseif branch.type == "C"
    nph = branch.nph
    range = 0:( nph - 1 )
    d_range = ( noff_d + branch.p_de ) .+ range
    i_range = ( num_nodes + branch.p_vars ) .+ range
    F[ d_range, i_range ] = sparse( -h2*I, nph, nph )
    F[ d_range, d_range ] = sparse( I, nph, nph )
  elseif branch.type == "NL"
    dpos = noff_d + branch.p_de
    vpos = noff_v + branch.p_vars
    F[ dpos, vpos ] = -h2
    F[ dpos, dpos ] = 1.0
  elseif branch.type == "NC"
    dpos = noff_d + branch.p_de
    ipos = num_nodes + branch.p_vars
    F[ dpos, vpos ] = -h2
    F[ dpos, dpos ] = 1.0
  end
end

end

max_stars_repo_name: takunoda/omegasim

max_issues_repo_name: takunoda/omegasim

max_forks_repo_name: takunoda/omegasim


text:
"""
**SegyToSeis**

*Convert SEGY or SU data to seis format.*

**IN**

* filename_in
* filename_out
* format="segy" (or "su")
* swap_bytes=true (flag to swap bytes)
* input_type="ibm" (or "ieee")

**OUT**

*Credits: AS, 2015*

"""
function SegyToSeis(filename_in,filename_out;format="segy",swap_bytes=true,input_type="ibm")

	if (format=="su")
		file_hsize = 0
	else
		file_hsize = 3600
		# add commands here to read text and binary headers and write them out to
		# filename_out.thead and filename_out.bhead
		stream     = open(filename_in)
		position = 3200
		seek(stream, position)
		fh = GrabFileHeader(stream)
		ntfh = swap_bytes == true ? bswap(fh.netfh) : fh.netfh
		if ntfh == -1
			error("add instructions to deal with variable extended text header")
		end
		if ntfh == 0
			file_hsize = 3600
		elseif ntfh > 0
			# file_hsize = 3200 * (ntfh+1) + 400
			file_hsize = 3200 * 1 + 400
		else
			error("unknown data format")
		end
	end
	stream = open(filename_in)
	seek(stream, segy_count["ns"] + file_hsize)
	if (swap_bytes==true)
		nt = bswap(read(stream,Int16))
	else
		nt = read(stream,Int16)
	end
	total = 60 + nt
	nx = round(Int,(filesize(stream)-file_hsize)/4/total)
	println("number of traces: ",nx)
	println("number of samples per trace: ",nt)
	h_segy = Array{SegyHeader}(undef,1)
	h_seis = Array{Header}(undef,1)
	seek(stream,file_hsize + segy_count["trace"])
	h_segy[1] = GrabSegyHeader(stream,swap_bytes,nt,file_hsize,1)
	dt = h_segy[1].dt/1000000
	extent = Extent(convert(Int32,nt),convert(Int32,nx),convert(Int32,1),convert(Int32,1),convert(Int32,1),
		   convert(Float32,0),convert(Float32,1),convert(Float32,0),convert(Float32,0),convert(Float32,0),
		   convert(Float32,dt),convert(Float32,1),convert(Float32,1),convert(Float32,1),convert(Float32,1),
		   "Time","Trace Number","","","",
		   "s","index","","","",
		   "")
	for j=1:nx
		position = file_hsize + total*(j-1)*4 + segy_count["trace"]
		seek(stream,position)
		if (input_type == "ieee")
			d = Array{Float32}(undef,nt);
			read!(stream,d)
		else
			d = Array{IBMFloat32}(undef,nt);
			read!(stream,d)
		end
		if (swap_bytes==true && input_type == "ieee")
			d = bswap_vector(d)
		end
		if (input_type != "ieee")
			d = convert(Array{Float32,1},d)
		end
		h_segy[1] = GrabSegyHeader(stream,swap_bytes,nt,file_hsize,j)
		h_seis[1] = MapHeaders(h_segy,j,"SegyToSeis")
		SeisWrite(filename_out,d,h_seis,extent,itrace=j)
	end
	close(stream)

end

function bswap_vector(a)
	for i = 1 : length(a)
		a[i] = bswap(a[i]);
	end
	return a
end