text:
function gendecomposition(manager::ExperimentManager, path::AbstractString,
                          options_src::AbstractDict, options_dst::AbstractDict,
                          nb_added_edges::Int;
                          seed=nothing, verbose=false,
                          dec_type=:cholesky, merge_type=:molzahn,
                          kwargs...)
    name = basename(path)[1:end - 2]
    graph = construct_network_graph(path)
    nb_added_edges = gen.add_edges!(graph, options_src, options_dst, nb_added_edges)
    options_src_features = nb_added_edges != 0 ? format_options(options_src) : Dict()
    options_dst_features = nb_added_edges != 0 ? format_options(options_dst) : Dict()
    cliques, cliquetree, chordal_graph, nb_lc, dec_data = gen.decompose(graph, dec_type; kwargs...)
    features = get_features_decomposition(chordal_graph, length(dec_data["added_edges"]), cliques, cliquetree)
    features["nb_lc"] = nb_lc
    data = Mongoc.BSON("cliques" => cliques,
                       "cliquetree" => cliquetree,
                       "options_src" => options_src_features,
                       "options_dst" => options_dst_features,
                       "nb_added_edges" => nb_added_edges,
                       "dec_type" => string(dec_type))
    if haskey(kwargs, "ordering")
        if typeof(kwargs["ordering"]) <: AbstractVector
            data["ordering"] = kwargs["ordering"]
        else
            data["ordering"] = string(kwargs["ordering"])
        end
    end
    for (key, val) in dec_data
        if typeof(val) <: Union{Number,AbstractString,Bool}
            features[key] = val
        elseif typeof(val) <: AbstractVector{T} where T <: Union{Number,AbstractString,Bool}
            data[key] = val
        end
    end
    instance_name, oid = db.push_decomposition!(manager.decompositions,
                                                name,
                                                data,
                                                features)
    db.add_decomposition_in_instance!(manager.instances, instance_name, oid)
    verbose && @info "($name, $oid, $nb_added_edges) inserted"
    return instance_name, oid
end

function gendecomposition_all(manager::ExperimentManager,
                              options_src::Dict, options_dst::Dict,
                              nb_added_edges::Int;
                              seed=nothing, verbose=false)
    instances_names = []
    oids = []
    paths_matpower = collect(db.getmatpowerpath_all(manager.instances))
    for path in paths_matpower
        instance_name, oid = gendecomposition(manager, path, options_src, options_dst, nb_added_edges; seed=seed, verbose=verbose)
        push!(instances_names, instance_name)
        push!(oids, oid)
    end
    return instances_names, oids
end

function gendecomposition_json(manager::ExperimentManager,
                               path::String, json_file::String;
                               verbose=false, dec_type=:cholesky,
                               kwargs...)
    config = JSON.parsefile(json_file)
    options_src = convertkeytosymbol(config["src"])
    options_dst = convertkeytosymbol(config["dst"])
    nb_added_edges = config["nb_added_edges"]
    seed = haskey(config, "seed") ? config["seed"] : 12345
    return gendecomposition(manager, path, options_src, options_dst, nb_added_edges; seed=seed, verbose=verbose, dec_type=dec_type)
end

function gendecomposition_all_json(manager::ExperimentManager,
                                   json_file::String; verbose=false,
                                   kwargs...)
    config = JSON.parsefile(json_file)
    options_src = convertkeytosymbol(config["src"])
    options_dst = convertkeytosymbol(config["dst"])
    nb_added_edges = config["nb_added_edges"]
    seed = haskey(config, "seed") ? config["seed"] : 12345
    return gendecomposition_all(manager, options_src, options_dst, nb_added_edges; seed=seed, verbose=verbose)
end

function gendecomposition_all_jsons(manager::ExperimentManager,
                                    json_file::String; verbose=false,
                                    kwargs...)
    ids = Dict()
    config = JSON.parsefile(json_file)
    for (key, val) in config
        options_src = convertkeytosymbol(val["src"])
        options_dst = convertkeytosymbol(val["dst"])
        nb_added_edges = val["nb_added_edges"]
        seed = haskey(val, "seed") ? val["seed"] : 12345
        ids[key] = gendecomposition_all(manager, options_src, options_dst, nb_added_edges; seed=seed, verbose=verbose)
    end
    return ids
end

function gendecomposition_mult_json(manager::ExperimentManager,
                                    json_file::String, paths_matpower::AbstractVector;
                                    verbose=false, kwargs...)
    config = JSON.parsefile(json_file)
    options_src = convertkeytosymbol(config["src"])
    options_dst = convertkeytosymbol(config["dst"])
    nb_added_edges = config["nb_added_edges"]
    seed = haskey(config, "seed") ? config["seed"] : 12345
    for path in paths_matpower
        gendecomposition(manager, path, options_src, options_dst, nb_added_edges; seed=seed, verbose=verbose)
    end
end

function genoriginal(manager::ExperimentManager, path; dec_type=:cholesky, verbose=false, kwargs...)
    return gendecomposition(manager, path, Dict(), Dict(), 0; dec_type=dec_type, verbose=verbose, kwargs...)
end

function genoriginal_mult(manager::ExperimentManager, path_matpower; dec_type=:cholesky, verbose=false, kwargs...)
    for path in paths_matpower
        genoriginal(manager, path; dec_type=dec_type, verbose=verbose, kwargs...)
    end
end

function genoriginal_all(manager::ExperimentManager; dec_type=:cholesky, verbose=false)
    paths_matpower = collect(db.getmatpowerpath_all(manager.instances))
    for path in paths_matpower
        genoriginal(manager, path; dec_type=dec_type, verbose=verbose, kwargs...)
    end
end

function genmerge_all(manager::ExperimentManager)
    for decomposition in manager.decompositions
        instance_name = decomposition["_id"]["instance_name"]
        oid = decomposition["_id"]["oid"]
        genmerge(manager, instance_name, oid; verbose=true)
    end
end

function genmerge_mult(manager::ExperimentManager, decompositions::AbstractVector; verbose=false)
    for decomposition in decompositions
        genmerge(manager, decomposition["instance_name"], decomposition["oid"])
    end
end

function genmerge(manager::ExperimentManager, instance_name::String,
                  oid::Union{String,UUID};
                  verbose=false,
                  merge_type=:molzahn,
                  dec_type=:cholesky,
                  merge_treshold_type=:nv,
                  merge_treshold::Float64=0.10,
                  kwargs...)
    decomposition = db.getdecomposition(manager.decompositions, instance_name, oid)
    cliques = copy(decomposition["cliques"])
    cliquetree = copy(decomposition["cliquetree"])
    paths = db.getpath_one(manager.instances, instance_name)
    alpha, beta, c = get_param_merge(paths["params"])
    array = merge_treshold_type == :nv ? Base.OneTo(decomposition["features"]["chordal_graph"]["nv"]) : cliques
    cliques, cliquetree, new_edges = gen.merge(cliques, cliquetree, merge_type;
                                               treshold=gen.treshold_percent(array, merge_treshold),
                                               alpha=alpha,
                                               beta=beta,
                                               m=get_nb_constraint(paths["OPF_ctr"]),
                                               nb_lc=decomposition["features"]["nb_lc"],
                                               kwargs...)
    chordal_graph = construct_graph_from_cliques(cliques)
    nb_added_edges = nv(chordal_graph) - db.getnv(manager.instances, instance_name)
    features = get_features_decomposition(chordal_graph, nb_added_edges, cliques, cliquetree)
    features["nb_lc"] = get_nb_lc(cliques, cliquetree)
    data = Mongoc.BSON("cliques" => cliques,
                       "cliquetree" => cliquetree,
                       "nb_added_edges" => -1,
                       "origin_oid" => oid,
                       "origin_nb_added_edges" => decomposition["nb_added_edges"],
                       "merge_treshold" => merge_treshold)
    instance_name, oid = db.push_decomposition!(manager.decompositions,
                                                instance_name,
                                                data,
                                                features)
    db.add_decomposition_in_instance!(manager.instances, instance_name, oid)
    db.add_merge_in_decomposition!(manager.decompositions,
                                   instance_name,
                                   decomposition["_id"]["oid"],
                                   string(oid),
                                   string(merge_type),
                                   string(merge_treshold_type),
                                   merge_treshold)
    verbose && @info "Merging of ($instance_name, $oid, $(decomposition["nb_added_edges"]) inserted"
    return instance_name, oid
end

function gencombine(manager::ExperimentManager, instance_name::String,
                    oid::Union{String,UUID}, k=3;
                    combine_type::Symbol=:intersect,
                    verbose=false,
                    neighbors_g=nothing)
    neighbors_indexes = db.getkneighbors_indexes(manager.decompositions, instance_name, oid, k, neighbors_g)
    neighbor = collect(manager.decompositions)[neighbors_indexes[rand(1:k)]]
    return gencombine(manager,
                      instance_name,
                      oid,
                      neighbor["_id"]["oid"];
                      combine_type=combine_type,
                      verbose=verbose)
end

function gencombine(manager::ExperimentManager, instance_name::String,
                    oid1::Union{String,UUID}, oid2::Union{String,UUID};
                    combine_type::Symbol=:intersect,
                    verbose=false)
    decomposition1 = db.getdecomposition(manager.decompositions, instance_name, oid1)
    decomposition2 = db.getdecomposition(manager.decompositions, instance_name, oid2)
    instance = db.getinstance(manager.instances, instance_name)
    g1 = construct_graph_from_cliques(decomposition1["cliques"])
    g2 = construct_graph_from_cliques(decomposition2["cliques"])
    combined_results = gen.combine_decomposition(g1, g2; combine_type=combine_type, instance_path=instance["paths"]["matpower"])
    isnothing(combined_results) && return nothing
    cliques, cliquetree, cg, nb_lc, data = combined_results
    merge!(data, Mongoc.BSON("cliques" => cliques,
                             "cliquetree" => cliquetree,
                             "nb_added_edges" => -2,
                             "oid1" => oid1,
                             "oid2" => oid2))
    nb_added_edges = nv(cg) - db.getnv(manager.instances, instance_name)
    features = get_features_decomposition(cg, nb_added_edges, cliques, cliquetree)
    instance_name, oid = db.push_decomposition!(manager.decompositions,
                                                instance_name,
                                                data,
                                                features)
    db.add_decomposition_in_instance!(manager.instances, instance_name, oid)
    verbose && @info "Combination of $(instance_name): $(oid1) and $(oid2) inserted"
    return instance_name, oid
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
@testset "NoiseWrapper" begin

using DiffEqNoiseProcess, Test, Random

_W = WienerProcess(0.0,0.0,0.0)

dt = 0.1
calculate_step!(_W,dt)

for i in 1:10
  accept_step!(_W,dt)
end

W2 = NoiseWrapper(_W)

dt = 0.1
calculate_step!(W2,dt)

for i in 1:10
  accept_step!(W2,dt)
end

_W = WienerProcess(0.0,0.0,0.0)

dt = 0.1
calculate_step!(_W,dt)

for i in 1:10
  accept_step!(_W,dt)
end

old_W = copy(_W.W)

W2 = NoiseWrapper(_W)
lspace =range(_W.t[1], stop=_W.t[end], length=1000)
dt = lspace[2]-lspace[1]
calculate_step!(W2,dt)
for t in lspace
  accept_step!(W2,dt)
end

@test W2.W[end] ≈ _W(W2.t[end])[1]
@test W2.Z[end] ≈ _W(W2.t[end])[2]

# Inplace

_W = WienerProcess!(0.0,zeros(4),zeros(4))

dt = 0.1
calculate_step!(_W,dt)

for i in 1:10
  accept_step!(_W,dt)
end

W2 = NoiseWrapper(_W)

dt = 0.1
calculate_step!(W2,dt)

for i in 1:10
  accept_step!(W2,dt)
end

_W = WienerProcess!(0.0,zeros(4),zeros(4))

dt = 0.1
calculate_step!(_W,dt)

for i in 1:10
  accept_step!(_W,dt)
end

W2 = NoiseWrapper(_W)
lspace =range(_W.t[1], stop=_W.t[end], length=20)
dt = lspace[2]-lspace[1]
calculate_step!(W2,dt)
for t in lspace
  accept_step!(W2,dt)
end

@test W2.W[end]≈ _W(W2.t[end])[1]
@test W2.Z[end]≈ _W(W2.t[end])[2]

@test W2.W[end] ≈ _W.W[end-1]

end

ext: jl


text:
using Dates: AbstractDateTime, argerror, validargs

# """
#     ZonedDateTime

# A `DateTime` that includes `TimeZone` information.
# """

struct ZonedDateTime <: AbstractDateTime
    utc_datetime::DateTime
    timezone::TimeZone
    zone::FixedTimeZone  # The current zone for the utc_datetime.

    function ZonedDateTime(utc_datetime::DateTime, timezone::TimeZone, zone::FixedTimeZone)
        return new(utc_datetime, timezone, zone)
    end

    function ZonedDateTime(utc_datetime::DateTime, timezone::VariableTimeZone, zone::FixedTimeZone)
        if timezone.cutoff !== nothing && utc_datetime >= timezone.cutoff
            throw(UnhandledTimeError(timezone))
        end

        return new(utc_datetime, timezone, zone)
    end
end

"""
    ZonedDateTime(dt::DateTime, tz::TimeZone; from_utc=false) -> ZonedDateTime

Construct a `ZonedDateTime` by applying a `TimeZone` to a `DateTime`. When the `from_utc`
keyword is true the given `DateTime` is assumed to be in UTC instead of in local time and is
converted to the specified `TimeZone`.  Note that when `from_utc` is true the given
`DateTime` will always exists and is never ambiguous.
"""
function ZonedDateTime(dt::DateTime, tz::VariableTimeZone; from_utc::Bool=false)
    possible = interpret(dt, tz, from_utc ? UTC : Local)

    num = length(possible)
    if num == 1
        return first(possible)
    elseif num == 0
        throw(NonExistentTimeError(dt, tz))
    else
        throw(AmbiguousTimeError(dt, tz))
    end
end

function ZonedDateTime(dt::DateTime, tz::FixedTimeZone; from_utc::Bool=false)
    utc_dt = from_utc ? dt : dt - tz.offset
    return ZonedDateTime(utc_dt, tz, tz)
end

"""
    ZonedDateTime(dt::DateTime, tz::VariableTimeZone, occurrence::Integer) -> ZonedDateTime

Construct a `ZonedDateTime` by applying a `TimeZone` to a `DateTime`. If the `DateTime` is
ambiguous within the given time zone you can set `occurrence` to a positive integer to
resolve the ambiguity.
"""
function ZonedDateTime(dt::DateTime, tz::VariableTimeZone, occurrence::Integer)
    possible = interpret(dt, tz, Local)

    num = length(possible)
    if num == 1
        return first(possible)
    elseif num == 0
        throw(NonExistentTimeError(dt, tz))
    elseif occurrence > 0
        return possible[occurrence]
    else
        throw(AmbiguousTimeError(dt, tz))
    end
end

"""
    ZonedDateTime(dt::DateTime, tz::VariableTimeZone, is_dst::Bool) -> ZonedDateTime

Construct a `ZonedDateTime` by applying a `TimeZone` to a `DateTime`. If the `DateTime` is
ambiguous within the given time zone you can set `is_dst` to resolve the ambiguity.
"""
function ZonedDateTime(dt::DateTime, tz::VariableTimeZone, is_dst::Bool)
    possible = interpret(dt, tz, Local)

    num = length(possible)
    if num == 1
        return first(possible)
    elseif num == 0
        throw(NonExistentTimeError(dt, tz))
    elseif num == 2
        mask = [isdst(zdt.zone.offset) for zdt in possible]

        # Mask is expected to be unambiguous.
        !xor(mask...) && throw(AmbiguousTimeError(dt, tz))

        occurrence = findfirst(d -> d == is_dst, mask)
        return possible[occurrence]
    else
        throw(AmbiguousTimeError(dt, tz))
    end
end

# Convenience constructors
@doc """
    ZonedDateTime(y, [m, d, h, mi, s, ms], tz, [amb]) -> DateTime

Construct a `ZonedDateTime` type by parts. Arguments `y, m, ..., ms` must be convertible to
`Int64` and `tz` must be a `TimeZone`. If the given `DateTime` is ambiguous in the given
`TimeZone` then `amb` can be supplied to resolve ambiguity.
""" ZonedDateTime

@optional function ZonedDateTime(y::Integer, m::Integer=1, d::Integer=1, h::Integer=0, mi::Integer=0, s::Integer=0, ms::Integer=0, tz::VariableTimeZone, amb::Union{Integer,Bool})
    ZonedDateTime(DateTime(y,m,d,h,mi,s,ms), tz, amb)
end

@optional function ZonedDateTime(y::Integer, m::Integer=1, d::Integer=1, h::Integer=0, mi::Integer=0, s::Integer=0, ms::Integer=0, tz::TimeZone)
    ZonedDateTime(DateTime(y,m,d,h,mi,s,ms), tz)
end

# Parsing constructor needed as part of the Dates parsing interface. Note we typically don't
# support passing in time zone information as a string since we cannot do not know if we
# need to support resolving ambiguity.
#
# Since we do not want users accidentially calling this function we'll use very specific
# type assertions:
# https://github.com/JuliaTime/TimeZones.jl/issues/187#issuecomment-473012078
function ZonedDateTime(y::Int64, m::Union{Int32, Int64}, d::Int64, h::Int64, mi::Int64, s::Int64, ms::Int64, tz::AbstractString)
    ZonedDateTime(DateTime(y,m,d,h,mi,s,ms), TimeZone(tz))
end

function ZonedDateTime(parts::Union{Period,TimeZone}...)
    periods = Period[]
    tz = nothing
    for part in parts
        if isa(part, Period)
            push!(periods, part)
        elseif tz === nothing
            tz = part
        else
            throw(ArgumentError("Multiple time zones found"))
        end
    end

    tz === nothing && throw(ArgumentError("Missing time zone"))
    return ZonedDateTime(DateTime(periods...), tz)
end

# Promotion

# Because of the promoting fallback definitions for TimeType, we need a special case for
# undefined promote_rule on TimeType types.
# Otherwise, typejoin(T,S) is called (returning TimeType) so no conversion happens, and
# isless(promote(x,y)...) is called again, causing a stack overflow.
function Base.promote_rule(::Type{T}, ::Type{S}) where {T<:TimeType, S<:ZonedDateTime}
    error("no promotion exists for ", T, " and ", S)
end

# Equality
Base.:(==)(a::ZonedDateTime, b::ZonedDateTime) = a.utc_datetime == b.utc_datetime
Base.isless(a::ZonedDateTime, b::ZonedDateTime) = isless(a.utc_datetime, b.utc_datetime)
Base.isequal(a::ZonedDateTime, b::ZonedDateTime) = isequal(a.utc_datetime, b.utc_datetime)

"""
    hash(::ZonedDateTime, h)

Compute an integer hash code for a ZonedDateTime by hashing the `utc_datetime` field.
`hash(:utc_instant, h)` is used to avoid collisions with `DateTime` hashes.
"""
function Base.hash(zdt::ZonedDateTime, h::UInt)
    h = hash(:utc_instant, h)
    h = hash(zdt.utc_datetime, h)
    return h
end

Base.typemin(::Type{ZonedDateTime}) = ZonedDateTime(typemin(DateTime), utc_tz; from_utc=true)
Base.typemax(::Type{ZonedDateTime}) = ZonedDateTime(typemax(DateTime), utc_tz; from_utc=true)

# Note: The `validargs` function is as part of the Dates parsing interface.
function Dates.validargs(::Type{ZonedDateTime}, y::Int64, m::Union{Int64, Int32}, d::Int64, h::Int64, mi::Int64, s::Int64, ms::Int64, tz::AbstractString)
    err = validargs(DateTime, y, Int64(m), d, h, mi, s, ms)
    err === nothing || return err
    istimezone(tz) || return argerror("TimeZone: \"$tz\" is not a recognized time zone")
    return argerror()
end


text:
const RETRYABLE_CODES = (
    "IncompleteSignature",
    "ThrottlingException",
    "RequestExpired",
)

function aws_retry_cond(s, e)
    if e isa AWSException && (500 <= e.cause.status <= 504 || e.code in RETRYABLE_CODES)
        debug(LOGGER, "CloudWatchLogs operation encountered $(e.code); retrying")
        return (s, true)
    elseif e isa MbedException
        debug(LOGGER, "CloudWatchLogs operation encountered $e; retrying")
        return (s, true)
    end

    return (s, false)
end

aws_retry(f) = retry(f, delays=GENERIC_AWS_DELAYS, check=aws_retry_cond)()

struct CloudWatchLogStream
    config::AWSConfig
    log_group_name::String
    log_stream_name::String
    token::Ref{Union{String, Nothing}}
end

"""
    CloudWatchLogStream(config::AWSConfig, log_group_name, log_stream_name)

Create a reference to a CloudWatch Log Stream on AWS with the log group name and log
stream name.
This constructor will automatically fetch the latest [sequence token](https://docs.aws.amazon.com/AmazonCloudWatchLogs/latest/APIReference/API_PutLogEvents.html#CWL-PutLogEvents-request-sequenceToken)
for the stream.
"""
function CloudWatchLogStream(
    config::AWSConfig,
    log_group_name::AbstractString,
    log_stream_name::AbstractString,
)
    stream = CloudWatchLogStream(
        config,
        log_group_name,
        log_stream_name,
        Ref{Union{String, Nothing}}(),
    )
    update_sequence_token!(stream)
    return stream
end

"""
    create_group(config::AWSConfig) -> String
    create_group(config::AWSConfig, log_group_name) -> String

Create a CloudWatch Log Group.
If the log group name is not provided, one is generated using a UUID4.

Returns the log group name.
"""
function create_group(
    config::AWSConfig,
    # this probably won't collide, most callers should add identifying information though
    log_group_name::AbstractString="julia-$(uuid4())";
    tags::AbstractDict{<:AbstractString, <:AbstractString}=Dict{String, String}(),
)
    if isempty(tags)
        aws_retry() do
            logs(config, "CreateLogGroup"; logGroupName=log_group_name)
        end
    else
        tags = Dict{String, String}(tags)
        aws_retry() do
            logs(config, "CreateLogGroup"; logGroupName=log_group_name, tags=tags)
        end
    end
    return String(log_group_name)
end

"""
    delete_stream(config::AWSConfig, log_group_name)

Delete a CloudWatch Log Group.
"""
function delete_group(
    config::AWSConfig,
    log_group_name::AbstractString,
)
    aws_retry() do
        logs(config, "DeleteLogGroup"; logGroupName=log_group_name)
    end
    return nothing
end

"""
    create_stream(config::AWSConfig, log_group_name) -> String
    create_stream(config::AWSConfig, log_group_name, log_stream_name) -> String

Create a CloudWatch Log Stream under a given Log Group.
If the log stream name is not provided, one is generated using a UUID4.

Returns the log stream name.
"""
function create_stream(
    config::AWSConfig,
    log_group_name::AbstractString,
    # this probably won't collide, most callers should add identifying information though
    log_stream_name::AbstractString="julia-$(uuid4())",
)
    aws_retry() do
        logs(
            config,
            "CreateLogStream";
            logGroupName=log_group_name,
            logStreamName=log_stream_name,
        )
    end
    return String(log_stream_name)
end

"""
    delete_stream(config::AWSConfig, log_group_name, log_stream_name)

Delete a CloudWatch Log Stream from a given Log Group.
"""
function delete_stream(
    config::AWSConfig,
    log_group_name::AbstractString,
    log_stream_name::AbstractString,
)
    aws_retry() do
        logs(
            config,
            "DeleteLogStream";
            logGroupName=log_group_name,
            logStreamName=log_stream_name,
        )
    end
    return nothing
end

"""
    sequence_token(config::AWSConfig, log_group_name, log_stream_name) -> Union{String, Nothing}

Return the current sequence token for the stream.
"""
sequence_token(stream::CloudWatchLogStream) = stream.token[]

function new_sequence_token(stream::CloudWatchLogStream)
    return new_sequence_token(stream.config, stream.log_group_name, stream.log_stream_name)
end

describe_log_streams(config; kwargs...) = logs(config, "DescribeLogStreams"; kwargs...)

"""
    new_sequence_token(stream::CloudWatchLogStream) -> Union{String, Nothing}
    new_sequence_token(config::AWSConfig, log_group_name, log_stream_name) -> Union{String, Nothing}

Fetch the current sequence token for the stream from AWS.

Returns `nothing` if the stream does not have a sequence token yet (e.g., if no events have
been logged).
"""
function new_sequence_token(
    config::AWSConfig,
    log_group::AbstractString,
    log_stream::AbstractString,
)::Union{String, Nothing}
    response = aws_retry() do
        @mock describe_log_streams(
            config;
            logGroupName=log_group,
            logStreamNamePrefix=log_stream,
            orderBy="LogStreamName",  # orderBy and limit will ensure we get just the one
            limit=1,                  # matching result
        )
    end

    streams = response["logStreams"]

    if isempty(streams) || streams[1]["logStreamName"] != log_stream
        msg = isempty(streams) ? nothing : "Did you mean $(streams[1]["logStreamName"])?"
        error(LOGGER, StreamNotFoundException(log_stream, log_group, msg))
    end

    return get(streams[1], "uploadSequenceToken") do
        debug(LOGGER) do
            string(
                "Log group '$log_group' stream '$log_stream' has no sequence token yet. ",
                "Using null as a default. ",
            )
        end

        return nothing
    end
end

"""
    update_sequence_token!(stream::CloudWatchLogStream) -> Union{String, Nothing}
    update_sequence_token!(stream::CloudWatchLogStream, token) -> Union{String, Nothing}

Fetch the current sequence token for the stream from AWS and store it.
Alternatively, set the token for the stream to `token`.

Returns the token.
"""
function update_sequence_token!(
    stream::CloudWatchLogStream,
    token::Union{String, Nothing}=new_sequence_token(stream),
)
    stream.token[] = token
end

function _put_log_events(stream::CloudWatchLogStream, events::AbstractVector{LogEvent})
    logs(
        stream.config,
        "PutLogEvents";
        logEvents=events,
        logGroupName=stream.log_group_name,
        logStreamName=stream.log_stream_name,
        sequenceToken=sequence_token(stream),
    )
end

"""
    submit_logs(stream::CloudWatchLogStream, events::AbstractVector{LogEvent}) -> Int

Submit a list of log events to AWS.

None of the log events can be more than 2 hours in the future, or older than 14 days or the
retention period of the log group.
If this occurs, those log messages will be rejected but the rest will succeed.

Submission of _all_ log events will fail if:

* the log events are more than 1 MiB of data
* the log events are not in chronological order by timestamp
* there are more than 10000 log events in `events`
* the log events span more than 24 hours

Returns the number of events successfully submitted.
"""
function submit_logs(stream::CloudWatchLogStream, events::AbstractVector{LogEvent})
    if length(events) > MAX_BATCH_LENGTH
        error(LOGGER, LogSubmissionException(
            "Log batch length exceeded 10000 events; submit fewer log events at once"
        ))
    end

    batch_size = sum(aws_size, events)

    if batch_size > MAX_BATCH_SIZE
        error(LOGGER, LogSubmissionException(
            "Log batch size exceeded 1 MiB; submit fewer log events at once"
        ))
    end

    if !issorted(events; by=timestamp)
        debug(LOGGER,
            "Log submission will be faster if log events are sorted by timestamp"
        )

        # a stable sort to avoid putting related messages out of order
        sort!(events; alg=MergeSort, by=timestamp)
    end

    min_timestamp, max_timestamp = extrema(timestamp(e) for e in events)
    if max_timestamp - min_timestamp > 24 * 3600 * 1000  # 24 hours in milliseconds
        error(LOGGER, LogSubmissionException(
            "Log events cannot span more than 24 hours; submit log events separately"
        ))
    end

    function retry_cond(s, e)
        if e isa AWSException
            if 500 <= e.cause.status <= 504 || e.code == "ThrottlingException"
                debug(LOGGER, "CloudWatchLogs PutLogEvents encountered $(e.code); retrying")
                return (s, true)
            elseif e.cause.status == 400 && e.code == "InvalidSequenceTokenException"
                debug(LOGGER) do
                    string(
                        "CloudWatchLogStream encountered InvalidSequenceTokenException. ",
                        "Are you logging to the same stream from multiple tasks?",
                    )
                end

                update_sequence_token!(stream)

                return (s, true)
            end
        elseif e isa MbedException
            debug(LOGGER, "CloudWatchLogs PutLogEvents encountered $e; retrying")
            return (s, true)
        end

        return (s, false)
    end

    f = retry(delays=PUTLOGEVENTS_DELAYS, check=retry_cond) do
        @mock _put_log_events(stream, events)
    end

    min_valid_event = 1
    max_valid_event = length(events)

    json_response = f()

    if haskey(json_response, "nextSequenceToken")
        update_sequence_token!(stream, json_response["nextSequenceToken"])
    end

    if haskey(json_response, "rejectedLogEventsInfo")
        rejected_info = json_response["rejectedLogEventsInfo"]

        if haskey(rejected_info, "expiredLogEventEndIndex")
            idx = Int(rejected_info["expiredLogEventEndIndex"]) + 1
            min_valid_event = max(min_valid_event, idx)

            warn(LOGGER) do
                string(
                    "Cannot log the following events, ",
                    "as they are older than the log retention policy allows: ",
                    events[1:idx-1],
                )
            end
        end

        if haskey(rejected_info, "tooOldLogEventEndIndex")
            idx = Int(rejected_info["tooOldLogEventEndIndex"]) + 1
            min_valid_event = max(min_valid_event, idx)

            warn(LOGGER) do
                string(
                    "Cannot log the following events, ",
                    "as they are more than 14 days old: ",
                    events[1:idx-1],
                )
            end
        end

        if haskey(rejected_info, "tooNewLogEventStartIndex")
            idx = Int(rejected_info["tooNewLogEventStartIndex"])
            max_valid_event = min(max_valid_event, idx)

            warn(LOGGER) do
                string(
                    "Cannot log the following events, ",
                    "as they are newer than 2 hours in the future: ",
                    events[idx+1:end],
                )
            end
        end
    end

    return length(min_valid_event:max_valid_event)
end

"""
    submit_log(stream::CloudWatchLogStream, event::LogEvent) -> Int

Call [`submit_logs`](@ref) with one event.
"""
submit_log(stream::CloudWatchLogStream, event::LogEvent) = submit_logs(stream, [event])


text:
#=##############################################################################
# DESCRIPTION
    Simulation interface connecting UVLMVehicle vehicles and maneuvers.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Oct 2019
  * License   : MIT
=###############################################################################

function solve(self::Simulation{V, M, R}, Vinf::Function,
                pfield::vpm.ParticleField, wake_coupled::Bool,
                dt::Real, rlx::Real, sigma::Real, rho::Real,
                speedofsound::Real, staticpfield::vpm.ParticleField;
                init_sol::Bool=false, vlm_fsgm=1
                ) where {V<:UVLMVehicle, M<:AbstractManeuver, R}

    # Bring the rotors back to their initial positions at beginning of this
    # time step before `precalculations`. This is needed to avoid double
    # accounting for the RPM when the solver is run
    rotate_rotors(self, -dt)

    vhcl = self.vehicle
    mnvr = self.maneuver
    t = self.t

    # TIME-STEPPING PROCEDURE:
    # -1) Solve one particle field time step
    # 0) Translate and rotate systems
    # 1) Recalculate horseshoes with kinematic velocity
    # 2) Paste previous Gamma solution to new system position after translation
    # 3) Shed semi-infinite wake after translation
    # 4) Calculate wake-induced velocity on VLM and Rotor system
    # 5) Solve VLM and Rotor system
    # 6) Shed unsteady-loading wake after new solution
    # 7) Save new solution as prev solution
    # Iterate
    #
    # On the first time step (pfield.nt==0), it only does steps (5) and (7),
    # meaning that the unsteady wake of the first time step is never shed.
    if self.nt==0
        # NOTE: VLMs and Rotors solutions are losely coupled

        # Set Vinf on system to get horseshoes
        vlm.setVinf(vhcl.system, Vinf)

        # Solve VLMs
        # TODO: Add Rotor-on-VLM induced velocity
        if init_sol
            # NOTE: Here I use the semi-infinite wake for the first step, which
            # may lead to some unphysical results when shedding unsteady loading wake
            vlm.solve(vhcl.vlm_system, Vinf; t=t, keep_sol=true)
        else
            vlm.solve(vhcl.vlm_system, Vinf; t=t, keep_sol=true,
                                                    vortexsheet=(X,t)->zeros(3))
        end

        # Calculate induced velocities to use in rotor solver
        ## Points where to calculate induced velocities
        Xs = _get_midXs(vhcl.rotor_systems)

        ## Particles for VLM-on-Rotor induced velocity
        static_particles_fun(pfield, args...) = _static_particles(pfield, vhcl.vlm_system, sigma)

        ## Evaluate VPM-on-Rotor induced velocity + static particles
        Vinds = Vvpm_on_Xs(pfield, Xs; static_particles_fun=static_particles_fun, dt=dt, fsgm=vlm_fsgm)

        # Solve Rotors
        for (si, rotors) in enumerate(vhcl.rotor_systems)

            # RPM of this rotor system
            RPM = self.RPMref*get_RPM(mnvr, si, t/self.ttot)

            for (ri, rotor) in enumerate(rotors)

                # Get velocities induced at every blade of this rotor
                this_Vinds = _parse_midXs(vhcl.rotor_systems,  Vinds, si, ri)

                VindVkin = _format_blades(this_Vinds, vhcl.rotor_systems, si, ri)

                if wake_coupled
                    vlm.solvefromV(rotor, VindVkin, Vinf, RPM, rho; t=t,
                                    sound_spd=speedofsound,
                                    tiploss_correction=vlm_fsgm==1)
                else
                    vlm.solvefromCCBlade(rotor, Vinf, RPM, rho;
                                              t=t, sound_spd=speedofsound)
                end
            end
        end

    else
        # ---------- 4) Calculate VPM velocity on VLM and Rotor system -----
        m = vlm.get_m(vhcl.vlm_system)

        # Control points of VLM
        Xs_cp_vlm = _get_Xs(vhcl.vlm_system, "CP"; t=t)

        # Mid-points along VLM filaments for calc_aerodynamicforce()
        ## Nodes of every horseshoe
        Ap = _get_Xs(vhcl.vlm_system, "Ap")
        A = _get_Xs(vhcl.vlm_system, "A")
        B = _get_Xs(vhcl.vlm_system, "B")
        Bp = _get_Xs(vhcl.vlm_system, "Bp")
        ## Midpoints of bound vortices
        ApA = (Ap .+ A)/2
        AB = (A .+ B)/2
        BBp = (B .+ Bp)/2
        ## Organize them
        Xs_ApA_AB_BBp_vlm = vcat(ApA, AB, BBp)

        # Points at rotors where to calculate induced velocities
        Xs_rotors = _get_midXs(vhcl.rotor_systems)

        # Calculate VPM velocity on all points (VLM and rotors)
        Vvpm = Vvpm_on_Xs(pfield, vcat(Xs_cp_vlm, Xs_ApA_AB_BBp_vlm, Xs_rotors); dt=dt, fsgm=vlm_fsgm)

        Vvpm_cp_vlm = Vvpm[1:m]
        Vvpm_ApA_AB_BBp_vlm = Vvpm[m+1:4*m]
        Vvpm_rotors = Vvpm[4*m+1:end]

        # Stitches previous rotor solutions for blade induced velocity
        prev_rotor_systems = _get_prev_rotor_systems(vhcl)
        allrotors = vlm.WingSystem()

        for (si, rotors) in enumerate(vhcl.rotor_systems)
            for (ri, rotor) in enumerate(rotors)

                vlm.addwing(allrotors, "S$(si)R$(ri)", rotor; reset=false)
                vlm.getHorseshoe(rotor, 1)          # Force HS calculation
                vlm._addsolution(rotor, "Gamma",    # Give previous solution
                                prev_rotor_systems[si][ri]._wingsystem.sol["Gamma"])
            end
        end

        # Particles for Rotor-on-VLM induced velocity
        static_particles_fun2(pfield, args...) = _static_particles(pfield, allrotors, sigma)

        # Evaluate Rotor-on-VLM induced velocity
        Vrotor_on_wing = Vvpm_on_Xs(staticpfield, vcat(Xs_cp_vlm, Xs_ApA_AB_BBp_vlm);
                                    static_particles_fun=static_particles_fun2, dt=dt, fsgm=vlm_fsgm)

        # Add and save VPM-on-VLM and Rotor-on-VLM induced velocity
        vlm._addsolution(vhcl.vlm_system, "Vvpm",
                                         Vvpm_cp_vlm + Vrotor_on_wing[1:m]; t=t)
        vlm._addsolution(vhcl.vlm_system, "Vvpm_ApA",
                        Vvpm_ApA_AB_BBp_vlm[1:m] + Vrotor_on_wing[m+1:2*m]; t=t)
        vlm._addsolution(vhcl.vlm_system, "Vvpm_AB",
                        Vvpm_ApA_AB_BBp_vlm[m+1:2*m] + Vrotor_on_wing[2*m+1:3*m]; t=t)
        vlm._addsolution(vhcl.vlm_system, "Vvpm_BBp",
                        Vvpm_ApA_AB_BBp_vlm[2*m+1:3*m] + Vrotor_on_wing[3*m+1:4*m]; t=t)

        # Calculate induced velocities to use in rotor solver
        ## Particles for VLM-on-Rotor induced velocity
        # NOTE: If I keep the rotor particles wouldn't I be double-counting
        # the blade induced velocity and make the solution unstable?
        # ANSWER: We are never double-accounting for the blade induced velocity
        function static_particles_fun3(pfield, args...)
            # Rotor static particles
            static_particles_fun2(pfield, args...)
            # VLM static particles
            _static_particles(pfield, vhcl.vlm_system, sigma)
        end

        ## Evaluate VLM-on-Rotor and Rotor-on-Rotor induced velocity
        Vvlmrotor_on_rotor = Vvpm_on_Xs(staticpfield, Xs_rotors;
                              static_particles_fun=static_particles_fun3, dt=dt, fsgm=vlm_fsgm)

        # Add VPM-on-Rotor, VLM-on-Rotor, and Rotor-on-Rotor induced velocity
        Vinds = Vvpm_rotors + Vvlmrotor_on_rotor

        # ---------- 5) Solve VLM system -----------------------------------
        # Wake-coupled solution
        if wake_coupled
            vlm.solve(vhcl.vlm_system, Vinf; t=t, extraVinf=_extraVinf2,
                                    keep_sol=true, vortexsheet=(X,t)->zeros(3))
        # Wake-decoupled solution
        else
            vlm.solve(vhcl.vlm_system, Vinf; t=t, extraVinf=_extraVinf1,
                                                                keep_sol=true)
        end

        # Relaxes (rlx->1) or stiffens (rlx->0) the VLM solution
        if rlx > 0
            rlxd_Gamma = rlx*vhcl.vlm_system.sol["Gamma"] +
                                (1-rlx)*_get_prev_vlm_system(vhcl).sol["Gamma"]
            vlm._addsolution(vhcl.vlm_system, "Gamma", rlxd_Gamma)
        end

        # ---------- 5) Solve Rotor system ---------------------------------

        # Solve Rotors
        for (si, rotors) in enumerate(vhcl.rotor_systems)

            # RPM of this rotor system
            RPM = self.RPMref*get_RPM(mnvr, si, t/self.ttot)

            for (ri, rotor) in enumerate(rotors)

                # Calculate kinematic velocities
                Vkin = _Vkinematic_rotor(vhcl.rotor_systems,
                                         prev_rotor_systems, si, ri, dt)

                # Get velocities induced at every blade of this rotor
                this_Vinds = _parse_midXs(vhcl.rotor_systems, Vinds, si, ri)

                VindVkin = _format_blades(this_Vinds + Vkin, vhcl.rotor_systems,
                                                                        si, ri)
                if wake_coupled
                    vlm.solvefromV(rotor, VindVkin, Vinf, RPM, rho; t=t,
                                    sound_spd=speedofsound,
                                    tiploss_correction=vlm_fsgm==1)
                else
                    vlm.solvefromCCBlade(rotor, Vinf, RPM, rho;
                                                t=t, sound_spd=speedofsound)
                end
            end
        end

        # Rotate rotors
        # NOTE: this means that the rotational velocity hadn't been included in
        # the kinematic velocities up to this point
        rotate_rotors(self, dt)
    end
end


text:
using LinearAlgebra: diag, norm, det
using Quantica: Hamiltonian, ParametricHamiltonian, nhoppings, nonsites, nsites, coordination, allsitepositions, axesflat

@testset "basic hamiltonians" begin
    presets = (LatticePresets.linear, LatticePresets.square, LatticePresets.triangular,
               LatticePresets.honeycomb, LatticePresets.cubic, LatticePresets.fcc,
               LatticePresets.bcc, LatticePresets.hcp)
    types = (Float16, Float32, Float64, ComplexF16, ComplexF32, ComplexF64)
    ts = (1, 2.0I, @SMatrix[1 2; 3 4], 1.0f0*I)
    orbs = (Val(1), Val(1), Val(2), (Val(1), Val(2)))
    for preset in presets, lat in (preset(), unitcell(preset())), type in types
        E, L = dims(lat)
        dn0 = ntuple(_ -> 1, Val(L))
        for (t, o) in zip(ts, orbs)
            @test hamiltonian(lat, onsite(t) + hopping(t; range = 1), orbitals = o, orbtype = type) isa Hamiltonian
            @test hamiltonian(lat, onsite(t) - hopping(t; dn = dn0), orbitals = o, orbtype = type) isa Hamiltonian
        end
    end
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = 1/√3))
    @test bloch(h) == h.harmonics[1].h
    # Inf range
    h = LatticePresets.square() |> unitcell(region = RegionPresets.square(5)) |>
        hamiltonian(hopping(1, range = Inf))
    @test nhoppings(h) == 600
    h = LatticePresets.square() |> hamiltonian(hopping(1, dn = (10,0), range = Inf))
    @test nhoppings(h) == 1
    @test isassigned(h, (10,0))
    h = LatticePresets.honeycomb() |> hamiltonian(onsite(1.0, sublats = :A), orbitals = (Val(1), Val(2)))
    @test Quantica.nonsites(h) == 1
    h = LatticePresets.square() |> unitcell(3) |> hamiltonian(hopping(1, indices = (1:8 .=> 2:9, 9=>1), range = 3, plusadjoint = true))
    @test nhoppings(h) == 48
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = (1, 1)))
    @test nhoppings(h) == 12
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = (1, 2/√3)))
    @test nhoppings(h) == 18
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = (2, 1)))
    @test Quantica.nhoppings(h) == 0
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = (30, 30)))
    @test Quantica.nhoppings(h) == 12
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = (10, 10.1)))
    @test Quantica.nhoppings(h) == 48
end

@testset "hamiltonian unitcell" begin
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1), region = r -> abs(r[2])<2)
    @test nhoppings(h) == 22
    h = LatticePresets.square() |> hamiltonian(hopping(1)) |> unitcell(3) |> unitcell((1,0), indices = not(1))
    @test nsites(h) == 8
    h = LatticePresets.square() |> hamiltonian(hopping(1, range = √2)) |> unitcell(5) |> unitcell((1,0), indices = 1:2:25)
    @test nhoppings(h) == 38
    @test nsites(h) == 13
    h = LatticePresets.square() |> hamiltonian(hopping(1, range = √2)) |> unitcell(5) |> unitcell(indices = 2:2:25)
    @test nhoppings(h) == 32
    @test nsites(h) == 12
    lat = LatticePresets.honeycomb(dim = Val(3)) |> unitcell(3) |> unitcell((1,1), indices = not(1))
    h = lat |> hamiltonian(hopping(1)) |> unitcell
    @test nsites(h) == 17
    h = unitcell(h, indices = not(1))
    @test nsites(h) == 16
    # dim-preserving unitcell reshaping should always preserve coordination
    lat = lattice(sublat((0.0, -0.1), (0,0), (0,1/√3), (0.0,1.6/√3)), bravais = SA[cos(pi/3) sin(pi/3); -cos(pi/3) sin(pi/3)]')
    h = lat |> hamiltonian(hopping(1, range = 1/√3))
    c = coordination(h)
    iter = CartesianIndices((-3:3, -3:3))
    for Ic´ in iter, Ic in iter
        sc = SMatrix{2,2,Int}(Tuple(Ic)..., Tuple(Ic´)...)
        iszero(det(sc)) && continue
        h´ = unitcell(h, sc)
        h´´ = unitcell(h´, 2)
        @test coordination(h´´) ≈ coordination(h´) ≈ c
        h´ = unitcell(h´, Tuple(Ic))
        h´´ = unitcell(h´, 2)
        @test coordination(h´´) ≈ coordination(h´)
    end
    h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell(2) |> unitcell(mincoordination = 2)
    @test nsites(h) == 6
    h = LP.cubic() |> hamiltonian(hopping(1)) |> unitcell(4) |> unitcell(mincoordination = 4)
    @test nsites(h) == 0
    h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell(region = RP.circle(5) & !RP.circle(2)) |> unitcell(mincoordination = 2)
    @test nsites(h) == 144
    h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell(10, region = !RP.circle(2, (0,8)))
    h´ = h |> unitcell(1, mincoordination = 2)
    @test nsites(h´) == nsites(h) - 1
    h = LP.square() |> hamiltonian(hopping(0)) |> unitcell(4, mincoordination = 2)
    @test nsites(h) == 0
    # check dn=0 invariant
    h = LP.linear() |> hamiltonian(hopping(1)) |> unitcell((1,))
    @test length(h.harmonics) == 3 && iszero(first(h.harmonics).dn)
end

@testset "hamiltonian wrap" begin
    h = LatticePresets.bcc() |> hamiltonian(hopping((r, dr) -> 1/norm(dr), range = 10))
    wh = wrap(h, phases = (1,2,3))
    @test bloch(wh) ≈ bloch(h, (1,2,3))
    h = LatticePresets.bcc() |> hamiltonian(hopping((r, dr) -> 1/norm(dr), range = 10)) |> unitcell(3)
    wh = wrap(h, phases = (1,2,3))
    @test bloch(wh) ≈ bloch(h, (1,2,3))
end

@testset "similarmatrix" begin
    types = (ComplexF16, ComplexF32, ComplexF64)
    lat = LatticePresets.honeycomb()
    for T in types
        h0 = hamiltonian(lat, onsite(I) + hopping(2I; range = 1), orbitals = (Val(1), Val(2)), orbtype = T)
        hf = flatten(h0)
        hm = Matrix(h0)
        hs = (h0, hf, hm)
        As = (SparseMatrixCSC, SparseMatrixCSC, Matrix)
        Es = (SMatrix{2,2,T,4}, T, SMatrix{2,2,T,4})
        for (h, A, E) in zip(hs, As, Es)
            sh = similarmatrix(h)
            @test sh isa A{E}
            b1 = bloch!(similarmatrix(flatten(h)), flatten(h), (1,1))
            b2 = bloch!(similarmatrix(h, flatten), h, (1,1))
            @test isapprox(b1, b2)
            for T´ in types
                E´s = E <: SMatrix ? (SMatrix{2,2,T´,4}, T´) : (T´,)
                for E´ in E´s
                    s1 = similarmatrix(h, Matrix{E´})
                    s2 = similarmatrix(h, Matrix)
                    @test s1 isa Matrix{E´}
                    @test s2 isa Matrix{E}
                    if A != Matrix
                        s1 = similarmatrix(h, SparseMatrixCSC{E´})
                        s2 = similarmatrix(h, SparseMatrixCSC)
                        @test s1 isa SparseMatrixCSC{E´}
                        @test s2 isa SparseMatrixCSC{E}
                    end
                end
            end
        end
    end

    h = LatticePresets.honeycomb() |> hamiltonian(hopping(I), orbitals = (Val(1), Val(2)))
    s = similarmatrix(h)
    @test size(s) == (2,2) && s isa SparseMatrixCSC{<:SMatrix{2,2}}
    s = similarmatrix(h, flatten)
    @test size(s) == (3,3) && s isa SparseMatrixCSC{ComplexF64}
    s = similarmatrix(h, SparseMatrixCSC)
    @test size(s) == (2,2) && s isa SparseMatrixCSC{<:SMatrix{2,2}}
    s = similarmatrix(h, SparseMatrixCSC{ComplexF16})
    @test size(s) == (3,3) && s isa SparseMatrixCSC{ComplexF16}
    s = similarmatrix(h, Matrix)
    @test size(s) == (2,2) && s isa Matrix{<:SMatrix{2,2}}
    s = similarmatrix(h, Matrix{Float64})
    @test size(s) == (3,3) && s isa Matrix{Float64}

    h = Matrix(h)
    s = similarmatrix(h)
    @test size(s) == (2,2) && s isa Matrix{<:SMatrix{2,2}}
    s = similarmatrix(h, flatten)
    @test size(s) == (3,3) && s isa Matrix{ComplexF64}
    @test_throws ArgumentError similarmatrix(h, SparseMatrixCSC)
    @test_throws ArgumentError similarmatrix(h, SparseMatrixCSC{ComplexF16})
    s = similarmatrix(h, Matrix)
    @test size(s) == (2,2) && s isa Matrix{<:SMatrix{2,2}}
    s = similarmatrix(h, Matrix{Float64})
    @test size(s) == (3,3) && s isa Matrix{Float64}

    h = LatticePresets.honeycomb() |> hamiltonian(hopping(I))
    s = similarmatrix(h)
    @test size(s) == (2,2) && s isa SparseMatrixCSC{ComplexF64}
    s = similarmatrix(h, flatten)
    @test size(s) == (2,2) && s isa SparseMatrixCSC{ComplexF64}
    s = similarmatrix(h, SparseMatrixCSC)
    @test size(s) == (2,2) && s isa SparseMatrixCSC{ComplexF64}
    s = similarmatrix(h, SparseMatrixCSC{ComplexF16})
    @test size(s) == (2,2) && s isa SparseMatrixCSC{ComplexF16}
    s = similarmatrix(h, Matrix)
    @test size(s) == (2,2) && s isa Matrix{ComplexF64}
    s = similarmatrix(h, Matrix{Float64})
    @test size(s) == (2,2) && s isa Matrix{Float64}

    h = Matrix(h)
    s = similarmatrix(h)
    @test size(s) == (2,2) && s isa Matrix{ComplexF64}
    s = similarmatrix(h, flatten)
    @test size(s) == (2,2) && s isa Matrix{ComplexF64}
    @test_throws ArgumentError similarmatrix(h, SparseMatrixCSC)
    @test_throws ArgumentError similarmatrix(h, SparseMatrixCSC{ComplexF16})
    s = similarmatrix(h, Matrix)
    @test size(s) == (2,2) && s isa Matrix{ComplexF64}
    s = similarmatrix(h, Matrix{Float64})
    @test size(s) == (2,2) && s isa Matrix{Float64}
end

@testset "orbitals and sublats" begin
    orbs = (:a, (:a,), (:a, :b, 3), ((:a, :b), :c), ((:a, :b), (:c,)), (Val(2), Val(1)),
            (:A => (:a, :b), :D => :c), :D => Val(4))
    lat = LatticePresets.honeycomb()
    for orb in orbs
        @test hamiltonian(lat, onsite(I), orbitals = orb) isa Hamiltonian
    end
    @test hamiltonian(lat, onsite(I) + hopping(@SMatrix[1 2], sublats = :A =>:B),
                      orbitals = :A => Val(2)) isa Hamiltonian
    h1 = hamiltonian(lat, onsite(I) + hopping(@SMatrix[1 2], sublats = :A =>:B),
                      orbitals = :A => Val(2))
    h2 = hamiltonian(lat, onsite(I) + hopping(@SMatrix[1 2], sublats = (:A =>:B,)),
                      orbitals = :A => Val(2))
    @test bloch(h1, (1, 2)) == bloch(h2, (1, 2))
end

@testset "onsite dimensions" begin
    lat = LatticePresets.honeycomb()
    ts = (@SMatrix[1 2; 3 4], r -> @SMatrix[r[1] 2; 3 4], 2, r -> 2r[1],
          @SMatrix[1 2; 3 4], r -> @SMatrix[r[1] 2; 3 4])
    os = (Val(1), Val(1), Val(2), Val(2), Val(3), Val(3))
    for (t, o) in zip(ts, os)
        model = onsite(t)
        @test_throws DimensionMismatch hamiltonian(lat, model, orbitals = o)
    end

    ts = (@SMatrix[1 2; 3 4], r -> @SMatrix[r[1] 2; 3 4], 2, r -> 2r[1], 3I, r -> I, 3I)
    os = (Val(2), Val(2), Val(1), Val(1), Val(3), Val(3), (Val(1), Val(3)))
    for (t, o) in zip(ts, os)
        model = onsite(t)
        @test hamiltonian(lat, model, orbitals = o) isa Hamiltonian
    end
    @test bloch(hamiltonian(lat, onsite(3I), orbitals = (Val(1), Val(3))))[1,1] ==
        @SMatrix[3 0 0; 0 0 0; 0 0 0]
    @test bloch(hamiltonian(lat, onsite(3I), orbitals = (Val(1), Val(3))))[2,2] ==
        SMatrix{3,3}(3I)
end

@testset "hopping dimensions" begin
    lat = LatticePresets.honeycomb()
    ts = (@SMatrix[1 2; 2 3], (r,dr) -> @SMatrix[r[1] 2; 3 4], 2, (r,dr) -> 2r[1],
          @SMatrix[1 2], @SMatrix[1 ;2], @SMatrix[1 2], @SMatrix[1 2; 2 3])
    os = (Val(1), Val(1), Val(2), Val(2), (Val(2), Val(1)), (Val(2), Val(1)),
         (Val(2), Val(1)), (Val(2), Val(1)))
    ss = (missing, missing, missing, missing, :B => :A, :A => :B, missing, missing)
    for (t, o, s) in zip(ts, os, ss)
        model = hopping(t, sublats = s)
        @test_throws DimensionMismatch hamiltonian(lat, model, orbitals = o)
    end
    ts = (@SMatrix[1 2], @SMatrix[1 ;2])
    os = ((Val(2), Val(1)), (Val(2), Val(1)))
    ss = (:A => :B, :B => :A)
    for (t, o, s) in zip(ts, os, ss)
        model = hopping(t, sublats = s)
        @test hamiltonian(lat, model, orbitals = o) isa Hamiltonian
    end
    @test bloch(hamiltonian(lat, hopping(3I, range = 1/√3), orbitals = (Val(1), Val(2))))[2,1] ==
        @SMatrix[3 0; 0 0]
    @test bloch(hamiltonian(lat, hopping(3I, range = 1/√3), orbitals = (Val(1), Val(2))))[1,2] ==
        @SMatrix[3 0; 0 0]
end

@testset "hermiticity" begin
    lat = LatticePresets.honeycomb()
    @test !ishermitian(hamiltonian(lat, hopping(im, sublats = :A=>:B)))
    @test !ishermitian(hamiltonian(lat, hopping(1, sublats = :A=>:B)))
    @test !ishermitian(hamiltonian(lat, hopping(1, sublats = :A=>:B, dn = (-1,0))))
    @test !ishermitian(hamiltonian(lat, hopping(im)))
    @test ishermitian(hamiltonian(lat, hopping(1)))

    @test ishermitian(hamiltonian(lat, hopping(im, sublats = :A=>:B, plusadjoint = true)))
    @test ishermitian(hamiltonian(lat, hopping(1, sublats = :A=>:B, plusadjoint = true)))
    @test ishermitian(hamiltonian(lat, hopping(1, sublats = :A=>:B, dn = (1,0), plusadjoint = true)))
    @test ishermitian(hamiltonian(lat, hopping(im, plusadjoint = true)))
    @test ishermitian(hamiltonian(lat, hopping(1, plusadjoint = true)))
end

@testset "unitcell modifiers" begin
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1) + onsite(0)) |> unitcell(2, modifiers = (@onsite!((o, r) -> 1), @hopping!(h -> 1)))
    @test diag(bloch(h)) == ComplexF64[1, 1, 1, 1, 1, 1, 1, 1]
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1) + onsite(0)) |> unitcell(2) |> unitcell(2, modifiers = @onsite!(o -> 1; indices = 3), indices = not(2))
    @test nonsites(h) == 4
    @test nsites(h) == 28
end

@testset "@onsite!" begin
    el = @SMatrix[1 2; 2 1]

    @test @onsite!(o -> 2o)(el) == 2el
    @test @onsite!(o -> 2o)(el, p = 2) == 2el
    @test @onsite!((o;) -> 2o)(el) == 2el
    @test @onsite!((o;) -> 2o)(el, p = 2) == 2el
    @test @onsite!((o; z) -> 2o)(el, z = 2) == 2el
    @test @onsite!((o; z) -> 2o)(el, z = 2, p = 2) == 2el
    @test @onsite!((o; z = 2) -> 2o)(el) == 2el
    @test @onsite!((o; z = 2) -> 2o)(el, p = 2) == 2el
    @test @onsite!((o; z = 2) -> 2o)(el, z = 1, p = 2) == 2el
    @test @onsite!((o; kw...) -> 2o)(el) == 2el
    @test @onsite!((o; kw...) -> 2o)(el, p = 2) == 2el
    @test @onsite!((o; z, kw...) -> 2o)(el, z = 2) == 2el
    @test @onsite!((o; z, kw...) -> 2o)(el, z = 2, p = 2) == 2el
    @test @onsite!((o; z, y = 2, kw...) -> 2o)(el, z = 2, p = 2) == 2el
    @test @onsite!((o; z, y = 2, kw...) -> 2o)(el, z = 2, y = 3, p = 2) == 2el

    r = SVector(0,0)

    @test @onsite!((o, r;) -> 2o)(el, r) == 2el
    @test @onsite!((o, r;) -> 2o)(el, r, p = 2) == 2el
    @test @onsite!((o, r; z) -> 2o)(el, r, z = 2) == 2el
    @test @onsite!((o, r; z) -> 2o)(el, r, z = 2, p = 2) == 2el
    @test @onsite!((o; z = 2) -> 2o)(el, r) == 2el
    @test @onsite!((o; z = 2) -> 2o)(el, r, p = 2) == 2el
    @test @onsite!((o; z = 2) -> 2o)(el, r, z = 1, p = 2) == 2el
    @test @onsite!((o, r; kw...) -> 2o)(el, r) == 2el
    @test @onsite!((o, r; kw...) -> 2o)(el, r, p = 2) == 2el
    @test @onsite!((o, r; z, kw...) -> 2o)(el, r, z = 2) == 2el
    @test @onsite!((o, r; z, kw...) -> 2o)(el, r, z = 2, p = 2) == 2el
    @test @onsite!((o, r; z, y = 2, kw...) -> 2o)(el, r, z = 2, p = 2) == 2el
    @test @onsite!((o, r; z, y = 2, kw...) -> 2o)(el, r, z = 2, y = 3, p = 2) == 2el

    @test @onsite!((o; z, y = 2, kw...) -> 2o) isa Quantica.UniformOnsiteModifier
    @test @onsite!((o, r; z, y = 2, kw...) -> 2o) isa Quantica.OnsiteModifier{2}

    @test parameters(@onsite!((o, r; z, y = 2, kw...) -> 2o)) == (:z, :y)
end

@testset "@hopping!" begin
    el = @SMatrix[1 2; 2 1]

    @test @hopping!(t -> 2t)(el) == 2el
    @test @hopping!(t -> 2t)(el, p = 2) == 2el
    @test @hopping!((t;) -> 2t)(el) == 2el
    @test @hopping!((t;) -> 2t)(el, p = 2) == 2el
    @test @hopping!((t; z) -> 2t)(el, z = 2) == 2el
    @test @hopping!((t; z) -> 2t)(el, z = 2, p = 2) == 2el
    @test @hopping!((t; z = 2) -> 2t)(el) == 2el
    @test @hopping!((t; z = 2) -> 2t)(el, p = 2) == 2el
    @test @hopping!((t; z = 2) -> 2t)(el, z = 1, p = 2) == 2el
    @test @hopping!((t; kw...) -> 2t)(el) == 2el
    @test @hopping!((t; kw...) -> 2t)(el, p = 2) == 2el
    @test @hopping!((t; z, kw...) -> 2t)(el, z = 2) == 2el
    @test @hopping!((t; z, kw...) -> 2t)(el, z = 2, p = 2) == 2el
    @test @hopping!((t; z, y = 2, kw...) -> 2t)(el, z = 2, p = 2) == 2el
    @test @hopping!((t; z, y = 2, kw...) -> 2t)(el, z = 2, y = 3, p = 2) == 2el

    r = dr = SVector(0,0)

    @test @hopping!((t, r, dr;) -> 2t)(el, r, dr) == 2el
    @test @hopping!((t, r, dr;) -> 2t)(el, r, dr, p = 2) == 2el
    @test @hopping!((t, r, dr; z) -> 2t)(el, r, dr, z = 2) == 2el
    @test @hopping!((t, r, dr; z) -> 2t)(el, r, dr, z = 2, p = 2) == 2el
    @test @hopping!((t; z = 2) -> 2t)(el, r, dr) == 2el
    @test @hopping!((t; z = 2) -> 2t)(el, r, dr, p = 2) == 2el
    @test @hopping!((t; z = 2) -> 2t)(el, r, dr, z = 1, p = 2) == 2el
    @test @hopping!((t, r, dr; kw...) -> 2t)(el, r, dr) == 2el
    @test @hopping!((t, r, dr; kw...) -> 2t)(el, r, dr, p = 2) == 2el
    @test @hopping!((t, r, dr; z, kw...) -> 2t)(el, r, dr, z = 2) == 2el
    @test @hopping!((t, r, dr; z, kw...) -> 2t)(el, r, dr, z = 2, p = 2) == 2el
    @test @hopping!((t, r, dr; z, y = 2, kw...) -> 2t)(el, r, dr, z = 2, p = 2) == 2el
    @test @hopping!((t, r, dr; z, y = 2, kw...) -> 2t)(el, r, dr, z = 2, y = 3, p = 2) == 2el

    @test @hopping!((t; z, y = 2, kw...) -> 2t) isa Quantica.UniformHoppingModifier
    @test @hopping!((t, r, dr; z, y = 2, kw...) -> 2t) isa Quantica.HoppingModifier{3}

    @test parameters(@hopping!((o, r, dr; z, y = 2, kw...) -> 2o)) == (:z, :y)
end

@testset "parametric" begin
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1) + onsite(2)) |> unitcell(10)
    T = typeof(h)
    @test parametric(h, @onsite!(o -> 2o))() isa T
    @test parametric(h, @onsite!((o, r) -> 2o))() isa T
    @test parametric(h, @onsite!((o, r; a = 2) -> a*o))() isa T
    @test parametric(h, @onsite!((o, r; a = 2) -> a*o))(a=1) isa T
    @test parametric(h, @onsite!((o, r; a = 2) -> a*o), @hopping!(t -> 2t))(a=1) isa T
    @test parametric(h, @onsite!((o, r) -> o), @hopping!((t, r, dr) -> r[1]*t))() isa T
    @test parametric(h, @onsite!((o, r) -> o), @hopping!((t, r, dr; a = 2) -> r[1]*t))() isa T
    @test parametric(h, @onsite!((o, r) -> o), @hopping!((t, r, dr; a = 2) -> r[1]*t))(a=1) isa T
    @test parametric(h, @onsite!((o, r; b) -> o), @hopping!((t, r, dr; a = 2) -> r[1]*t))(b=1) isa T
    @test parametric(h, @onsite!((o, r; b) -> o*b), @hopping!((t, r, dr; a = 2) -> r[1]*t))(a=1, b=2) isa T

    ph = LatticePresets.linear() |> hamiltonian(hopping(1)) |> parametric(@onsite!((o; k) -> o + k*I))
    # No onsites, no need to specify k
    @test ph() isa Hamiltonian
    ph = LatticePresets.linear() |> hamiltonian(onsite(0)) |> parametric(@onsite!((o; k) -> o + k*I))
    ph´ = LatticePresets.linear() |> hamiltonian(onsite(0)) |> parametric(@onsite!((o; k = 1) -> o + k*I))
    @test_throws UndefKeywordError ph()
    @test bloch(ph, (;k = 1)) == bloch(ph(k = 1)) == bloch(ph´)

    # Issue #35
    for orb in (Val(1), Val(2))
        h = LatticePresets.triangular() |> hamiltonian(hopping(I) + onsite(I), orbitals = orb) |> unitcell(10)
        ph = parametric(h, @onsite!((o, r; b) -> o+b*I), @hopping!((t, r, dr; a = 2) -> t+r[1]*I),
                       @onsite!((o, r; b) -> o-b*I), @hopping!((t, r, dr; a = 2) -> t-r[1]*I))
        @test isapprox(bloch(ph(a=1, b=2), (1, 2)), bloch(h, (1, 2)))
    end
    # Issue #37
    for orb in (Val(1), Val(2))
        h = LatticePresets.triangular() |> hamiltonian(hopping(I) + onsite(I), orbitals = orb) |> unitcell(10)
        ph = parametric(h, @onsite!(o -> o*cis(1)))
        @test ph()[][1,1] ≈ h[][1,1]*cis(1)
    end
    # Issue #54. Parametric Haldane model
    sK(dr::SVector) = sK(atan(dr[2],dr[1]))
    sK(ϕ) = 2*mod(round(Int, 6*ϕ/(2π)), 2) - 1
    ph = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = 1)) |>
         parametric(@hopping!((t, r, dr; λ) ->  λ*im*sK(dr); sublats = :A=>:A),
                    @hopping!((t, r, dr; λ) -> -λ*im*sK(dr); sublats = :B=>:B))
    @test bloch(ph(λ=1), (π/2, -π/2)) == bloch(ph, (π/2, -π/2, (;λ=1))) ≈ [4 1; 1 -4]
    # Non-numeric parameters
    ph = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = 1)) |>
         parametric(@hopping!((t, r, dr; λ, k) ->  λ*im*sK(dr+k); sublats = :A=>:A),
                    @hopping!((t, r, dr; λ, k) -> -λ*im*sK(dr+k); sublats = :B=>:B))
    @test bloch(ph(λ=1, k=SA[1,0]), (π/2, -π/2)) == bloch(ph, (π/2, -π/2, (;λ = 1, k = SA[1,0]))) ≈ [-4 1; 1 4]
    # Issue 61, type stability
    h = LatticePresets.honeycomb() |> hamiltonian(onsite(0))
    @inferred parametric(h, @onsite!((o;μ) -> o- μ))
    @inferred parametric(h, @onsite!(o->2o), @hopping!((t)->2t), @onsite!((o, r)->o+r[1]))
    @inferred parametric(h, @onsite!((o, r)->o*r[1]), @hopping!((t; p)->p*t), @onsite!((o; μ)->o-μ))

    h = LatticePresets.honeycomb() |> hamiltonian(hopping(2I), orbitals = (Val(2), Val(1)))
    ph = parametric(h, @hopping!((t; α, β = 0) -> α * t .+ β))
    b = bloch!(similarmatrix(ph, flatten), ph, (0, 0, (; α = 2)))
    @test b == [0 0 12; 0 0 0; 12 0 0]

    # @block! modifiers
    h = LatticePresets.honeycomb() |> hamiltonian(onsite(0I) + hopping(I), orbitals = (Val(1), Val(3))) |> unitcell(4)
    sites = siteindices(h, region = r->r[1]<0);
    ph = h |> parametric(@onsite!((o; λ=0) -> o - λ*I), @block!((b; λ=0) -> b + λ*I, sites))
    @test_throws ArgumentError ph(λ = 1)
    ph = h |> parametric(@onsite!((o; λ=0) -> o - λ*I), @block!((b; λ=0) -> b + λ*I, sites); check = false)
    @test ph(λ = 1) isa Hamiltonian
    sites = siteindices(h, region = r->r[1]<0, sublats = :B);
    ph = h |> parametric(@onsite!((o; λ=0) -> o - λ*I; sublats = :B), @block!((b; λ=0) -> b + λ*I, sites))
    h3 = ph(λ = 3)
    for i in Quantica.siterange(h.lattice, 1)
        @test iszero(h3[][i, i])
    end
    for i in Quantica.siterange(h.lattice, 2)
        @test i in sites ? iszero(h3[][i,i]) : h3[][i,i] == -3I
    end

    h = LatticePresets.square() |> hamiltonian(onsite(0I) + hopping(I)) |> unitcell(4)
    rows = siteindices(h, region = r->r[1]==0);
    cols = siteindices(h, region = r->r[1]==3);
    ph = h |> parametric(
        @block!((b; λ=0) -> b + λ*I, rows, cols; dn = ((0,0),(1,0))),
        @block!((b; λ=0) -> b + λ*I, cols, rows; dn = (-1,0)))
    @test ph(λ = 10)[(1,0)][collect(rows), collect(cols)] == 11I
    @test ph(λ = 10)[(-1,0)][collect(cols), collect(rows)] == 11I
    # check that reset works by calling twice
    @test ph(λ = 10)[(1,0)][collect(rows), collect(cols)] == 11I
    @test ph(λ = 10)[(-1,0)][collect(cols), collect(rows)] == 11I
end

@testset "boolean masks" begin
    for b in ((), (1,1), 4)
        h1 = LatticePresets.honeycomb() |> hamiltonian(hopping(1) + onsite(2)) |>
             supercell(b, region = RegionPresets.circle(10))
        h2 = LatticePresets.honeycomb() |> hamiltonian(hopping(1) + onsite(2)) |>
             supercell(b, region = RegionPresets.circle(20))

        @test isequal(h1 & h2, h1)
        @test isequal(h1, h2) || !isequal(h1 & h2, h2)
        @test isequal(h1, h2) || !isequal(h1 | h2, h1)
        @test  isequal(h1 | h2, h2)

        @test isequal(unitcell(h1 & h2), unitcell(h1))
        @test isequal(h1, h2) || !isequal(unitcell(h1 & h2), unitcell(h2))
        @test isequal(h1, h2) || !isequal(unitcell(h1 | h2), unitcell(h1))
        @test isequal(unitcell(h1 | h2), unitcell(h2))

        h1 = h1.lattice
        h2 = h2.lattice

        @test isequal(h1 & h2, h1)
        @test isequal(h1, h2) || !isequal(h1 & h2, h2)
        @test isequal(h1, h2) || !isequal(h1 | h2, h1)
        @test  isequal(h1 | h2, h2)

        @test isequal(unitcell(h1 & h2), unitcell(h1))
        @test isequal(h1, h2) || !isequal(unitcell(h1 & h2), unitcell(h2))
        @test isequal(h1, h2) || !isequal(unitcell(h1 | h2), unitcell(h1))
        @test isequal(unitcell(h1 | h2), unitcell(h2))
    end
end

@testset "unitcell seeds" begin
    p1 = SA[100,0]
    p2 = SA[0,20]
    lat = LatticePresets.honeycomb()
    model = hopping(1, range = 1/√3) + onsite(2)
    h1 = lat |> hamiltonian(model) |> supercell(region = r -> norm(r-p1)<3, seed = p1)
    h2 = lat |> hamiltonian(model) |> supercell(region = r -> norm(r-p2)<3, seed = p2)
    h = unitcell(h1 | h2)
    h3 = lat |> hamiltonian(model) |> unitcell(region = r -> norm(r-p1)<3 || norm(r-p2)<3, seed = p2)

    @test Quantica.nsites(h) == 130
    @test Quantica.nsites(h3) == 64
end

@testset "nrange" begin
    lat = LatticePresets.honeycomb(a0 = 2)
    @test nrange(1, lat) ≈ 2/√3
    @test nrange(2, lat) ≈ 2
    @test nrange(3, lat) ≈ 4/√3
    @test hamiltonian(lat, hopping(1)) |> nhoppings == 6
    @test hamiltonian(lat, hopping(1, range = nrange(2))) |> nhoppings == 18
    @test hamiltonian(lat, hopping(1, range = nrange(3))) |> nhoppings == 24
    @test hamiltonian(lat, hopping(1, range = (nrange(2), nrange(3)))) |> nhoppings == 18
    @test hamiltonian(lat, hopping(1, range = (nrange(20), nrange(20)))) |> nhoppings == 18
    @test hamiltonian(lat, hopping(1, range = (nrange(20), nrange(21)))) |> nhoppings == 30
end

@testset "hamiltonian algebra" begin
    h1 = LatticePresets.honeycomb() |> hamiltonian(onsite(-I) + hopping(I), orbitals = (Val(2), Val(1))) |> unitcell(4)
    h2 = LatticePresets.honeycomb() |> hamiltonian(hopping(3I, range = 1), orbitals = (Val(2), Val(1))) |> unitcell(4)
    phis = (0.3, 0.8)
    @test bloch(3.3*h1, phis) ≈ 3.3*bloch(h1, phis)
    @test bloch(3.3*h1 + h2, phis) ≈ 3.3*bloch(h1, phis) + bloch(h2, phis)
    @test bloch(h1^2, phis) ≈ bloch(h1, phis)^2
    h´ = 2.3h1^3 - 2h2*h1 - 3I
    @test Quantica.check_orbital_consistency(h´) === nothing
    @test bloch(flatten(h´), phis) ≈ 2.3bloch(flatten(h1), phis)^3 - 2bloch(flatten(h2), phis)*bloch(flatten(h1), phis) - 3I
end

@testset "transform! hamiltonians" begin
    h = LP.honeycomb(dim = 3) |> hamiltonian(hopping(1))
    h1 = copy(h)
    h2 = transform!(r -> SA[1 2 0; 2 3 0; 0 0 1] * r + SA[0,0,1], h1)
    h3 = h1 |> transform!(r -> SA[1 2 0; 2 3 0; 0 0 1] * r + SA[0,0,1])
    @test h1 === h2 === h3
    @test all(r->r[3] == 2.0, allsitepositions(h3.lattice))
    @test bloch(h, (1,2)) == bloch(h3, (1,2))
end

@testset "combine hamiltonians" begin
    h1 = LP.square(dim = Val(3)) |> hamiltonian(hopping(1))
    h2 = transform!(r -> r + SA[0,0,1], copy(h1))
    h = combine(h1, h2; coupling = hopping((r,dr) -> exp(-norm(dr)), range = √2))
    @test coordination(h) == 9
    h0 = LP.honeycomb(dim = Val(3), names = (:A, :B)) |> hamiltonian(hopping(1))
    ht = LP.honeycomb(dim = Val(3), names = (:At, :Bt)) |> hamiltonian(hopping(1)) |> transform!(r -> r + SA[0,1/√3,1])
    hb = LP.honeycomb(dim = Val(3), names = (:Ab, :Bb)) |> hamiltonian(hopping(1)) |> transform!(r -> r + SA[0,1/√3,-1])
    h = combine(hb, h0, ht; coupling = hopping((r,dr) -> exp(-norm(dr)), range = 2,
        sublats = ((:A,:B) => (:At, :Bt), (:A,:B) => (:Ab, :Bb)),  plusadjoint = true))
    @test iszero(bloch(h)[1:2, 5:6])
    h = combine(hb, h0, ht; coupling = hopping((r,dr) -> exp(-norm(dr)), range = 2))
    @test !iszero(bloch(h)[1:2, 5:6])
end

@testset "hamiltonian slices" begin
    function testslices(h::Hamiltonian, inds, φs)
        hs = h[inds...]
        @test bloch(hs, φs) ≈ bloch(h, φs)[axes(hs)...]
        @test bloch!(similarmatrix(hs, Matrix), hs, φs) ≈ bloch!(similarmatrix(h, Matrix), h, φs)[axes(hs)...]
        @test bloch!(similarmatrix(hs, Matrix{ComplexF64}), hs, φs) ≈ bloch!(similarmatrix(h, Matrix{ComplexF64}), h, φs)[axesflat(hs)...]
        @test bloch!(similarmatrix(hs, flatten), hs, φs) ≈ bloch!(similarmatrix(h, flatten), h, φs)[axesflat(hs)...]
        inds´ = notscalar.(inds, Ref(h))
        @test h[][inds´...] ≈ h[inds...][]
    end

    function testslices(h::ParametricHamiltonian, inds, φs)
        hs = h[inds...]
        @test bloch(hs, φs) ≈ bloch(h, φs)[axes(hs)...] ≈
            bloch(h(), φs)[axes(hs)...]
        @test bloch!(similarmatrix(hs, Matrix), hs, φs) ≈ bloch!(similarmatrix(h, Matrix), h, φs)[axes(hs)...] ≈
            bloch!(similarmatrix(h, Matrix), h(), φs)[axes(hs)...]
        @test bloch!(similarmatrix(hs, Matrix{ComplexF64}), hs, φs) ≈ bloch!(similarmatrix(h, Matrix{ComplexF64}), h, φs)[axesflat(hs)...] ≈
            bloch!(similarmatrix(h, Matrix{ComplexF64}), h(), φs)[axesflat(hs)...]
        @test bloch!(similarmatrix(hs, flatten), hs, φs) ≈ bloch!(similarmatrix(h, flatten), h, φs)[axesflat(hs)...] ≈
            bloch!(similarmatrix(h, flatten), h(), φs)[axesflat(hs)...]
    end

    notscalar(s::Integer, h) = s:s
    notscalar(s::Quantica.SiteSelector, h) = collect(siteindices(h, s))
    notscalar(s, h) = s

    function alltests(h, φs)
        testslices(h, (70, 2), φs)
        testslices(h, (:, 2), φs)
        testslices(h, (:, 70), φs)
        testslices(h, (2, :), φs)
        testslices(h, (70, :), φs)
        testslices(h, (1:4, 10:70), φs)
        testslices(h, (1:4, :), φs)
        testslices(h, (:, 10:70), φs)
        s1 = siteselector(; region = r -> 0<r[2]<7)
        s2 = siteselector(; sublats = :B, region = r -> 0<r[2]<7)
        s3 = siteselector(; indices = 2:70)
        testslices(h, (:, s1), φs)
        testslices(h, (s1, :), φs)
        testslices(h, (s1, s2), φs)
        testslices(h, (s3, s2), φs)
        hs = h[s3, s3]
        @test axes(hs) == (2:70, 2:70)
    end

    h = LatticePresets.honeycomb() |> hamiltonian(onsite(1I) + hopping(I), orbitals = (Val(3),Val(2))) |> unitcell((3,-4), region = r -> 0<r[2]<10)
    alltests(h, .2)
    h = LatticePresets.honeycomb() |> hamiltonian(onsite(1I) + hopping(I)) |> unitcell(20)
    alltests(h, (.2, .3))
    ph = LatticePresets.honeycomb() |> hamiltonian(onsite(1I) + hopping(I), orbitals = (Val(3),Val(2))) |> unitcell((3,-4), region = r -> 0<r[2]<10)
        parametric(@onsite!((o; k = 2) -> k*o))
    alltests(ph, .2)
    ph = LatticePresets.honeycomb() |> hamiltonian(onsite(1I) + hopping(I)) |> unitcell(20) |>
        parametric(@onsite!((o; k = 2) -> k*o))
    alltests(ph, (.2, .3))
end

max_stars_repo_name: pablosanjose/Quantica

max_issues_repo_name: pablosanjose/Quantica

max_forks_repo_name: pablosanjose/Quantica


text:
"""
    gram_schmidt(A::Array[, ret_coef_mat::Bool])

Return orthogonalized basis.
Optionally returns the Gram-Schmidt coefficient matrix using the ret_coef_mat parameter.

# Example

```julia
julia> x = [
        [ 3  -1   5]
        [-5   2  -1]
        [-3   9   2]
    ]

julia> gram_schmidt(x)
3×3 Array{Float64,2}:
  3.0      -1.0       5.0
 -3.11429   1.37143   2.14286
  2.68728   6.5689   -0.298587

```
"""
function gram_schmidt{T<:AbstractFloat}(A::Array{T, 2}; ret_coef_mat::Bool = false)
    # Assume the vectors are given as row vectors of X
    # Maybe redo later since vectors in Julia are stored in column-major order.

    X = transpose(copy(A))
    X_star = zeros(X)

    # Coefficient matrix
    μ = zeros(X)

    n_cols = size(X)[2]
    for i in 1:n_cols
        X_star[:, i] = X[:, i]
        for j in 1:i-1
            t = dot(X[:, i], X_star[:, j]) / dot(X_star[:, j], X_star[:, j])
            μ[j, i] = t
            X_star[:, i] -= (t * X_star[:, j])
        end
        μ[i,i] = 1
    end

    if ret_coef_mat
        transpose(X_star), transpose(μ)
    else
        transpose(X_star)
    end
end

function gram_schmidt{T<:Integer}(A::Array{T, 2}; ret_coef_mat::Bool = false)
    gram_schmidt(float(A), ret_coef_mat=ret_coef_mat)
end


text:
#
# Interface to DifferentialEquations.jl
#

#
# Evaluate function and Jacobian at u if they have not
# been evaluated before at u.
#
# See https://github.com/SciML/DifferentialEquations.jl/issues/521
# for discussion of another way to do this
#
function _eval_res_jac!(sys,u)
    uhash=hash(u)
    if uhash!=sys.uhash
        ur=reshape(u,sys)
        eval_and_assemble(sys,ur,ur,sys.residual,Inf)
        sys.uhash=uhash
    else
        global nd
        nd=nd+1
    end
end

"""
$(SIGNATURES)

Interpret the  discrete problem as an ODE/DAE problem. Provide the
rhs function for DifferentialEquations.jl.
"""
function eval_rhs!(du, u, sys,t)
    global nf
    nf=nf+1
    _eval_res_jac!(sys,u)
    du.=-vec(sys.residual)
    nothing
end

njac=0
nf=0
nd=0
"""
$(SIGNATURES)

Interpret the  discrete problem as an ODE/DAE problem. Provide the
jacobi matrix calculation function for DifferentialEquations.jl.
"""
function eval_jacobian!(J, u, sys,t)
    global njac
    njac=njac+1
    _eval_res_jac!(sys,u)
    # Need to implement broadcast for ExtendableSparse.
    J.=-sys.matrix.cscmatrix
    nothing
end

"""
$(SIGNATURES)

Calculate the mass matrix for use with DifferentialEquations.jl.
"""
function mass_matrix(system::AbstractSystem{Tv, Ti}) where {Tv, Ti}
    physics=system.physics
    grid=system.grid
    data=physics.data
    node=Node{Tv,Ti}(system)
    bnode=BNode{Tv,Ti}(system)
    nspecies=num_species(system)
    cellnodefactors=system.cellnodefactors
    bfacenodefactors=system.bfacenodefactors
    bgeom=grid[BFaceGeometries][1]
    bfacenodes=grid[BFaceNodes]
    bfaceregions=grid[BFaceRegions]
    nbfaces=num_bfaces(grid)

    if isdata(data)
        isbstorage=(physics.bstorage!=nofunc)

        storagewrap= function(y, u)
            y.=0
            physics.storage(y,u,node,data)
            nothing
        end

        bstoragewrap=function(y, u)
            y.=0
            physics.bstorage(y,u,bnode,data)
            nothing
        end

    else
        isbstorage=(physics.bstorage!=nofunc2)

        storagewrap= function(y, u)
            y.=0
            physics.storage(y,u,node)
            nothing
        end
        bstoragewrap=function(y, u)
            y.=0
            physics.bstorage(y,u,bnode)
            nothing
        end
    end

    result_s=DiffResults.DiffResult(Vector{Tv}(undef,nspecies),Matrix{Tv}(undef,nspecies,nspecies))
    Y=Array{Tv,1}(undef,nspecies)
    UK=Array{Tv,1}(undef,nspecies)
    cfg_s=ForwardDiff.JacobianConfig(storagewrap, Y, UK)
    cfg_bs=ForwardDiff.JacobianConfig(bstoragewrap, Y, UK)

    U=unknowns(system,inival=1)
    M=unknowns(system,inival=0)

    geom=grid[CellGeometries][1]
    cellnodes=grid[CellNodes]
    cellregions=grid[CellRegions]
    bfacenodes=grid[BFaceNodes]

    for icell=1:num_cells(grid)
        for inode=1:num_nodes(geom)
            _fill!(node,cellnodes,cellregions,inode,icell)
            for ispec=1:nspecies
                UK[ispec]=U[ispec,node.index]
            end
            ForwardDiff.jacobian!(result_s,storagewrap,Y,UK,cfg_s)
            res_stor=DiffResults.value(result_s)
            jac_stor=DiffResults.jacobian(result_s)
            K=node.index
            for idof=_firstnodedof(M,K):_lastnodedof(M,K)
                ispec=_spec(M,idof,K)
                M[ispec,K]+=cellnodefactors[inode,icell]*jac_stor[ispec,ispec]
            end
        end
    end

    if isbstorage
        for ibface=1:nbfaces
            for ibnode=1:num_nodes(bgeom)
                _fill!(bnode,bfacenodes,bfaceregions,ibnode,ibface)
                for ispec=1:nspecies
                    UK[ispec]=U[ispec,bnode.index]
                end
                ForwardDiff.jacobian!(result_s,bstoragewrap,Y,UK,cfg_bs)
                jac_bstor=DiffResults.jacobian(result_s)
                for idof=_firstnodedof(M,K):_lastnodedof(M,K)
                    ispec=_spec(M,idof,K)
                    M[ispec,K]+=bfacenodefactors[ibnode,ibface]*jac_bstor[ispec,ispec]
                end

            end
        end
    end

    return Diagonal(vec(M))
end

"""
$(SIGNATURES)

Complete the system and provide the jacobi matrix as prototype
for the jacobian for use with `DifferentialEquations.jl`.
"""
function jac_prototype(sys::AbstractSystem)
    _complete!(sys, create_newtonvectors=true)
    sys.matrix.cscmatrix
end

"""
        solve mit Union{Module, Nothing}
    package=DifferentialEquations
    solve(inival,sys,tspan,pkg=..., solver=...)

    saveat,tstops,tspan vs tsample
    => tstops=tsample
       tspan=[tstops[1],tstops[end]]
       save=:all, :sample,:nothing
    FunctionCallingCallback for collecting data
    SavingCallback
    TerminateSteadyState
    keep evolve ?
"""

"""
````
solve(DifferentialEquations, inival, system, tspan;
      solver=nothing,
      kwargs...)
````

Use a  timestepping scheme  from [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl) to perform
transient solution  of the  system.

The system must have a constant diagonal mass matrix.

Any stiff solver  capable to work with mass matrices  is possible.  If
the system  contains elliptic equations, it  needs to be a  DAE solver
able to work with mass matrices.  See a list
[here](https://diffeq.sciml.ai/stable/solvers/dae_solve/#dae_solve_full).

The  constructed  solver should  be  passed  as `solver=`  kwarg.   By
default (if `nothing` is passed), the `Rosenbrock23` solver is used.

All other  kwargs are  passed to  `DifferentialEquations.solve()`, see
[here](https://diffeq.sciml.ai/stable/basics/common_solver_opts/)  for
an overview.

In order to avoid a package dependency on DifferentialEquations,  the module
corresponding to that package is passed as the first argument.

Returns a transient solution object `sol` containing stored solutions,
see [`TransientSolution`](@ref).

"""
function solve(DiffEq::Module,
               inival::AbstractArray,
               sys::AbstractSystem,
               tspan; solver=nothing, kwargs...)

    if isnothing(solver)
        solver=DiffEq.Rosenbrock23()
    end
    global njac
    global nf
    global nd

    njac=0
    nf=0
    nd=0
    f = DiffEq.ODEFunction(eval_rhs!,
                           jac=eval_jacobian!,
                           jac_prototype=jac_prototype(sys),
                           mass_matrix=mass_matrix(sys))

    prob = DiffEq.ODEProblem(f,vec(inival),tspan,sys)
    sol = DiffEq.solve(prob,solver; kwargs...)
    # Return solution as TransientSolution such that sol[i] adheres to the
    # different storage schemes for multispecies cases.
    TransientSolution([reshape(sol.u[i],sys) for i=1:length(sol.u)] ,sol.t)
end


text:
using Test, TopicModels, Random
using TopicModels: updateSufficientStatistics, joint_log_p #non-exported fns we need

# use the equality of likelihood ratio to test that the conditional distribution is consistent with the joint distribution
@testset "LDA docs" begin
    # generate some data from LDA where the doclength is Poisson
    k = 7
    lexLength = 10
    corpLambda = 10 # poisson parameter for random doc length
    corpLength = 10
    scaleK = 0.1
    scaleL = 0.1
    Random.seed!(123)

    corpus = LdaCorpus(k, lexLength, corpLambda, corpLength, scaleK, scaleL)

    model = Model(corpus.alpha, corpus.beta, corpus)
    state = State(model, corpus)
    trainModel(model, state, 10) # update all the state variables

    # pick a random doc/word to iterate the sampler
    doc_ind = rand(1:corpLength)
    word_ind = rand(1:length(corpus.docs[doc_ind]))
    word = corpus.docs[doc_ind].terms[word_ind]

    conditional = state.conditionals[doc_ind][word_ind,:]
    oldTopic = copy(state.assignments[doc_ind][word_ind])  # the original word token

    newTopic = rand(collect(1:k)[1:end .!= oldTopic],1) # a different word token
    newTopic = Int64(newTopic[1])

    #get the original state probs
    joint_Lw = copy(joint_log_p(model,state)) # log prob of the full joint under original topic for <doc><word>
    cond_Lw = log(state.conditionals[doc_ind][word_ind,oldTopic]/sum(state.conditionals[doc_ind][word_ind,:])) # log conditional p(z=k|...)
    cond_Lw_new = log(state.conditionals[doc_ind][word_ind,newTopic]/sum(state.conditionals[doc_ind][word_ind,:])) # log conditional p(z=k|...)

    updateSufficientStatistics(word, oldTopic, doc_ind, -model.corpus.weights[doc_ind][word_ind], state) #remove counts for the old topic
    updateSufficientStatistics(word, newTopic, doc_ind, model.corpus.weights[doc_ind][word_ind], state) #update stats for new topic
    joint_Lw_new = copy(joint_log_p(model,state)) # log prob of the full joint under original topic for <doc><word>

    print("joint_Lw: ", joint_Lw, "\n")
    print("cond_Lw: ", cond_Lw, "\n")

    print("joint_Lw_new: ", joint_Lw_new, "\n")
    print("cond_Lw_new: ", cond_Lw_new, "\n")

    print("joint_LR: ", joint_Lw_new-joint_Lw, "\n")
    print("cond_LR: ", cond_Lw_new-cond_Lw, "\n")
    print("old Topic: ", oldTopic, "\n")
    print("new Topic: ", newTopic, "\n")

    @test isless(abs(joint_Lw_new-joint_Lw - cond_Lw_new+cond_Lw),1e-5)
end


text:
using RayTracingInOneJulia
using Test

@testset "RayTracingInOneJulia.jl" begin
    # Write your tests here.
end


text:
using Documenter, StorageGraphs

makedocs(;
    modules=[StorageGraphs],
    format=Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages=[
        "Home" => "index.md",
        "Adding data" => "add.md",
        "Querying the graph" => "query.md",
        "Internals" => "internals.md",
    ],
    repo="https://github.com/SebastianM-C/StorageGraphs.jl/blob/{commit}{path}#L{line}",
    sitename="StorageGraphs.jl",
    authors="Sebastian Micluța-Câmpeanu",
    assets=[],
)

deploydocs(;
    repo="github.com/SebastianM-C/StorageGraphs.jl",
)


text:
using DiffEqJump, DiffEqBase
using Test

function regular_rate(out,u,p,t)
    out[1] = (0.1/1000.0)*u[1]*u[2]
    out[2] = 0.01u[2]
end

function regular_c(dc,u,p,t,mark)
    dc[1,1] = -1
    dc[2,1] = 1
    dc[2,2] = -1
    dc[3,2] = 1
end

dc = zeros(3,2)

rj = RegularJump(regular_rate,regular_c,dc;constant_c=true)
jumps = JumpSet(rj)

prob = DiscreteProblem([999.0,1.0,0.0],(0.0,250.0))
jump_prob = JumpProblem(prob,Direct(),rj)
sol = solve(jump_prob,SimpleTauLeaping();dt=1.0)
sol = solve(jump_prob,RegularSSA())


text:
# PELT algorithm

@doc """
# Description
Runs the PELT algorithm using a specified cost function for a given penalty value to find the position and number of changepoints

# Arguments
* `segment_cost::Function`: Calculates cost between two specified indices
* `n::Int`: Length of time series
* `pen::Float64`: This is the penalty value used to avoid over/under -fitting the model.

# Returns
* `(CP::Vector{Int}, cost::Float64)`:
  * `CP::Vector{Int}`: Vector of indices of detected changepoints
  * `cost::Float64`: Cost of optimal segmentation

# Example
```
n, λ = 1000, 100
μ, σ = Normal(0.0, 10.0), 1.0
sample, cps = @changepoint_sampler n λ Normal(μ, σ)  # Sample changepoints
seg_cost = NormalMeanChange(sample, σ)               # Create segment cost function
pelt_cps, pelt_cost = PELT(seg_cost, n)              # Run PELT
```

# See also
CROPS, @PELT

# References
Killick, R., Fearnhead, P. and Eckley, I.A. (2012) Optimal detection of changepoints with a linear computational cost, JASA 107(500), 1590-1598
""" ->
function PELT( segment_cost::Function , n::Int; pen::Float64 = log(n) )

    # F[t] is optimal cost of segmentation upto time t
    F = Array(Float64, n+1)
    F[1] = - pen
    F[2] = 0

    # last chpt prior to time t
    chpts =  Array(Int64, n)
    chpts[1] = 0
    # vector of candidate chpts at t
    R = Array(Int64, 1)
    R = [0]

    for t in 2:n
        cpt_cands = R
        seg_costs = Array(Float64, length(cpt_cands))
        for i in 1:length(cpt_cands)
            seg_costs[i] = segment_cost(cpt_cands[i], t)
        end

        F[t+1] , tau = findmin( F[cpt_cands+1] + seg_costs + pen )
        chpts[t] = cpt_cands[tau]

        # pruning step
        ineq_prune = F[cpt_cands+1] + seg_costs .< F[t+1]
        R =  push!( cpt_cands[ineq_prune] , t-1 )

    end

    # get changepoints
    CP = Array(Int64,0)
    last = chpts[n]
    push!(CP,last)
    while last > 0
      last = chpts[last]
      push!(CP,last)
    end
    sort!(CP)

    return CP , F[n+1]

end


text:
using Test, SegmentIntersections

@testset "Test queue" begin
    Q = EventQueue()
    s1 = Segment(1, 2, 3, 4)
    s2 = Segment(5, 6, 3, 4)
    s3 = Segment(4, 3, 2, 7)
    s4 = Segment(4, 3, 2, 4)
    s5 = Segment(0, 1, 2, 7)
    insert!(Q, s1)
    insert!(Q, s2)
    insert!(Q, s3)
    insert!(Q, s4)
    insert!(Q, s5)
    e = fetch_event!(Q)
    @test e.point == Point(2,7)
    @test e.segments == Set([s3, s5])
    e = fetch_event!(Q)
    @test e.point == Point(5,6)
    @test e.segments == Set([s2])
    e = fetch_event!(Q)
    @test e.point == Point(2,4)
    @test e.segments == Set([s4])
    e = fetch_event!(Q)
    @test e.point == Point(3,4)
    @test e.segments == Set([s1])
end


text:
#! /usr/bin/julia

# Rosetta Code, Variable-length quantity

type VLQ
    q::Array{Uint8,1}
end

function VLQ{T<:Integer}(n::T)
    q = uint8(digits(n, 128))
    for i in 2:length(q)
        q[i] |= 0x80
    end
    VLQ(reverse(q))
end

function Base.uint(vlq::VLQ)
    q = reverse(vlq.q)
    n = shift!(q)
    p = one(Uint64)
    for i in q
        p *= 0x80
        n += p*(i&0x7f)
    end
    return n
end

test = [0x00200000, 0x001fffff, 0x00000000, 0x0000007f,
        0x00000080, 0x00002000, 0x00003fff, 0x00004000,
        0x08000000, 0x0fffffff]

for i in test
    q = VLQ(i)
    j = uint(q)
    print(@sprintf "    0x%x => " i)
    print(@sprintf "[%s]" join(["0x"*hex(r, 2) for r in q.q], ", "))
    println(@sprintf " => 0x%x" j)
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using Statistics
using JLD2
using Printf
using GLMakie

# using Oceananigans.Grids
# using Oceananigans.Utils: prettytime, hours, day, days, years

function geographic2cartesian(λ, φ, r=1)
    λ_azimuthal = λ .+ 180  # Convert to λ ∈ [0°, 360°]
    φ_azimuthal = 90 .- φ   # Convert to φ ∈ [0°, 180°] (0° at north pole)

    x = @. r * cosd(λ_azimuthal) * sind(φ_azimuthal)
    y = @. r * sind(λ_azimuthal) * sind(φ_azimuthal)
    z = @. r * cosd(φ_azimuthal)

    return x, y, z
end

ds_cs = jldopen("cubed_sphere_face_waves.jld2")

λᶜᶜᵃ = ds_cs["grid/λᶜᶜᵃ"][2:33, 2:33]
φᶜᶜᵃ = ds_cs["grid/φᶜᶜᵃ"][2:33, 2:33]

xc, yc, zc = geographic2cartesian(λᶜᶜᵃ, φᶜᶜᵃ)

iterations = parse.(Int, keys(ds_cs["timeseries/t"]))

iter = Node(0)

plot_title = @lift @sprintf("η′ on a cubed sphere face: time = %s", prettytime(ds_cs["timeseries/t/" * string($iter)]))

η = @lift ds_cs["timeseries/η/" * string($iter)][:, :, 1]

fig = Figure(resolution = (1920, 1080))

ax = fig[1, 1] = LScene(fig) # make plot area wider
wireframe!(ax, Sphere(Point3f0(0), 0.99f0), show_axis=false)
sf = surface!(ax, xc, yc, zc, color=η, colormap=:balance, colorrange=(-0.01, 0.01))
rotate_cam!(ax.scene, (3π/4, π/6, 0))
zoom!(ax.scene, (0, 0, 0), 5, false)
# fig[2, 2 + 3*(n-1)] = Label(fig, statenames[n], textsize = 50) # put names in center

cb1 = fig[1, 2] = Colorbar(fig, sf, label="η′", width=30)

supertitle = fig[0, :] = Label(fig, plot_title, textsize=50)

record(fig, "cubed_sphere_waves.mp4", iterations, framerate=60) do i
    @info "Animating iteration $i/$(iterations[end])..."
    iter[] = i
end

close(ds_cs)


text:
using ReTest
using Banyan, BanyanArrays
using FilePathsBase, AWSS3, HDF5
using LibGit2

global jobs_for_testing = Dict()

function destroy_all_jobs_for_testing()
    global jobs_for_testing
    for (job_config_hash, job_id) in jobs_for_testing
        destroy_job(job_id)
        delete!(jobs_for_testing, job_config_hash)
    end
end

function get_branch_name()
    prepo = LibGit2.GitRepo(realpath(joinpath(@__DIR__, "../..")))
    phead = LibGit2.head(prepo)
    branchname = LibGit2.shortname(phead)
    branchname
end

function use_job_for_testing(
    f::Function;
    sample_rate = 2,
    max_exact_sample_length = 50,
    with_s3fs = nothing,
    scheduling_config_name = "default scheduling",
)
    haskey(ENV, "BANYAN_CLUSTER_NAME") || error(
        "Please specify the Banyan cluster to use for testing with the BANYAN_CLUSTER_NAME environment variable",
    )

    # This will be a more complex hash if there are more possible ways of
    # configuring a job for testing. Different sample rates are typically used
    # to test different data sizes. Stress tests may need a much greater sample
    # rate.
    job_config_hash = sample_rate

    # Set the job and create a new one if needed
    global jobs_for_testing
    set_job(
        if haskey(jobs_for_testing, job_config_hash)
            jobs_for_testing[job_config_hash]
        else
            create_job(
                cluster_name = ENV["BANYAN_CLUSTER_NAME"],
                nworkers = 2,
                sample_rate = sample_rate,
                print_logs = true,
                url = "https://github.com/banyan-team/banyan-julia.git",
                branch = get(ENV, "BANYAN_JULIA_BRANCH", get_branch_name()),
                directory = "banyan-julia/BanyanArrays/test",
                dev_paths = [
                    "banyan-julia/Banyan",
                    "banyan-julia/BanyanArrays"
                ],
                force_pull = get(ENV, "BANYAN_FORCE_CLONE", "0") == "0",
                force_clone = get(ENV, "BANYAN_FORCE_CLONE", "0") == "1",
                force_install = get(ENV, "BANYAN_FORCE_INSTALL", "0") == "1",
                store_logs_on_cluster=get(ENV, "BANYAN_STORE_LOGS_ON_CLUSTER", "0") == "1"
            )
        end
    )
    # If selected job has already failed, this will throw an error.
    jobs_for_testing[job_config_hash] = get_job_id()

    # Set the maximum exact sample length
    ENV["BANYAN_MAX_EXACT_SAMPLE_LENGTH"] = string(max_exact_sample_length)

    # Force usage of S3FS if so desired
    if !isnothing(with_s3fs)
        ENV["BANYAN_USE_S3FS"] = with_s3fs ? "1" : "0"
    end

    configure_scheduling(name = scheduling_config_name)

    try
        f()
    catch
        # We will destroy the job if any error occurs. This is because we can't
        # properly intercept errors that happen in tests. If an error occurs,
        # the whole test suite exits and we don't have an opportunity to delete
        # stray jobs. This ensures that jobs are destroyed. In later tests sets,
        # `get_job()` is called which ensures that the job hasn't yet been
        # destroyed or failed.
        destroy_all_jobs_for_testing()
        rethrow()
        # If no errors occur, we will destroy all jobs in the `finally...` block.
    end
end

global data_for_testing = false

function use_data(data_src = "S3")
    global data_for_testing

    if !data_for_testing && data_src == "S3"
        original = h5open(
            download(
                "https://support.hdfgroup.org/ftp/HDF5/examples/files/exbyapi/h5ex_d_fillval.h5",
            ),
        )
        with_downloaded_path_for_reading(
            joinpath(
                S3Path("s3://$(get_cluster_s3_bucket_name())", config = Banyan.get_aws_config()),
                "fillval.h5",
            ),
            for_writing = true,
        ) do f
            new = h5open(string(f), "w")
            new["DS1"] = repeat(original["DS1"][:, :], 100, 100)
            close(new)
        end
        close(original)

        # rm(get_s3fs_path(joinpath(get_cluster_s3_bucket_name(), "fillval_copy.h5")), force=true)
        rm(
            joinpath(
                S3Path("s3://$(get_cluster_s3_bucket_name())", config = Banyan.get_aws_config()),
                "fillval_copy.h5",
            ),
            force = true,
        )
        # TODO: Maybe fsync here so that the directory gets properly updated
    end

    data_for_testing = true
end

include("sample_computation.jl")
include("mapreduce.jl")
include("hdf5.jl")
include("black_scholes.jl")

try
    runtests(Regex.(ARGS)...)
finally
    # Destroy jobs to clean up
    destroy_all_jobs_for_testing()
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
"""
A superpackage for handling common astrodynamics
problems. See the __Extended Help__ section
for more information!

# Extended Help

## License
$(LICENSE)

## Imports
$(IMPORTS)

## Exports
$(EXPORTS)
"""
module GeneralAstrodynamics

using Reexport, Requires
using DocStringExtensions

@template (FUNCTIONS, METHODS, MACROS) =
    """
    $(SIGNATURES)
    $(DOCSTRING)
    $(METHODLIST)
    """

@template (TYPES, CONSTANTS) =
    """
    $(TYPEDEF)
    $(DOCSTRING)
    """

include(joinpath("Calculations", "Calculations.jl"))
include(joinpath("CoordinateFrames", "CoordinateFrames.jl"))
include(joinpath("States", "States.jl"))

@reexport using .Calculations
@reexport using .CoordinateFrames
@reexport using .States

function __init__()
    @require DifferentialEquations="0c46a032-eb83-5123-abaf-570d42b7fbaa" begin
        include(joinpath("Propagation", "Propagation.jl"))
        @reexport using .Propagation
    end
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
        include(joinpath("Visualizations", "Visualizations.jl"))
        @reexport using .Visualizations
    end
end

end # module


text:
using LBLRTM
using Unitful
using Test

@testset "Planck function and its inverse" begin
    ν = 1000u"cm^-1"
    T_0 = 300u"K"
    L = planck(ν, T_0)
    T_1 = radiance_to_bt(ν, L)
    T_2 = radiance_to_bt(ν)(L)

    @test T_0 == T_1
    @test T_0 == T_2
end


text:
function generateTransFluxPlot( fluxReport::DataFrame, timeFactor::Float64, showPlot::Bool, savePlot::Bool, filename::String )

    # Generate the plot.
    plotTitle = string( "Flux plot of transitions" )
    plotData = Matrix( fluxReport[:, 3:end] )
    labels = hcat( string.( names( fluxReport )[3:end] )... )
    ymax = max( maximum( plotData ), 20 )

    plt = plot( fluxReport[:, :timePoint] / timeFactor, plotData, show=false,
        size=(960, 540), title=plotTitle, label=labels, lw=2, ylim=[0, ymax] )

    # Show plot if needed.
    if showPlot
        gui( plt )
    end  # if showPlot

    # Save plot if needed.
    if savePlot
        # If the user doesn't specify a filename, a random one is generated.
        if filename == ""
            dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
            filename = string( "Transition flux plot  ", dateStr )
        end  # if filename == ""

        if !ispath( dirname( filename ) )
            mkpath( dirname( filename ) )
        end  # if !ispath( dirname( filename ) )

        extension = split( filename, "." )
        extension = length( extension ) == 1 ? "" : extension[end]
        filename = string( filename, extension ∈ extensions ? "" : ".svg" )
        savefig( plt, filename )
    end  # if savePlot

end  # generateTransFluxPlot( fluxReport, timeFactor, showPlot, savePlot,
     #   filename )


text:
"""
Suffixes in csv file names created by StanSample.jl.

$(SIGNATURES)

### Required arguments
```julia
* `model`                    : SampleModel
```

Returns a vector with available chain suffixes.
"""
function available_chains(m::SampleModel)
  suffix_array = AbstractString[]
  for i in 1:m.num_julia_chains   # Number of exec processes
      for k in 1:m.num_cpp_chains   # Number of cpp chains handled in cmdstan

        if (m.use_cpp_chains && m.check_num_chains) ||
          !m.use_cpp_chains || m.num_cpp_chains == 1
          if m.use_cpp_chains && m.num_cpp_chains > 1
            csvfile_suffix = "$(k)"
          else
            #if m.use_cpp_chains && m.num_cpp_chains == 1
              csvfile_suffix = "$(i)"
            #else
            #  csvfile_suffix = "$(k)"
            #end
          end
        else
          if i == 1
            csvfile_suffix = "$(i)_$(k)"
          else
            csvfile_suffix = "$(i)_$(k + i - 1)"
          end
        end
        append!(suffix_array, [csvfile_suffix])
      end
  end
  Dict(:chain => collect(1:length(suffix_array)), :suffix => suffix_array)
end

export
  available_chains

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
using Test
using JLD2
using Symbolics
using LinearAlgebra
using Random
using BenchmarkTools
using FiniteDiff
using RoboDojo

# Utils
include("utils.jl")

# Solver
include("solver/lu.jl")
include("solver/random_qp.jl")
include("solver/soc.jl")

# Robots
include("robots/hopper.jl")
include("robots/biped.jl")
include("robots/quadruped.jl")

# Simulator
include("simulator/utils.jl")
include("simulator/simulate.jl")
include("simulator/codegen.jl")
include("simulator/dynamics.jl")


text:
# compare output from KernelDensity.jl with sinc kernel
# to what Comte-Butucea yields. (should be the same)
# also look at how DeLaValleePoussin performs

using MinimaxCalibratedEBayes
using Test
using Distributions
using Random

prior_dist = MixtureModel(Normal, [(-0.3,0.5), (1.05,0.5)])
marginal_dist = marginalize(prior_dist, StandardNormalSample(NaN))

n = 100_000
Random.seed!(1)
Xs = rand(marginal_dist, n)
extrema(Xs)

certainty_kde_res_dlv = certainty_banded_KDE(Xs, -4, 4; kernel=DeLaValleePoussinKernel)
certainty_kde_res_sinc = certainty_banded_KDE(Xs, -4, 4; kernel=SincKernel)

fitted_kde = certainty_kde_res[:fitted_kde]

f_true = pdf.(marginal_dist, fitted_kde.x)
C∞ = certainty_kde_res[:C∞]*1.01

true_C = maximum(abs.(f_true .- fitted_kde.density))

using MosekTools

xgrid = -3:0.01:3
σ = 0.5
gmix_tmp = GaussianMixturePriorClass(σ, xgrid, Mosek.Optimizer)

Z1 = DiscretizedStandardNormalSample(-6:0.01:6)
LinCalibSinc = DiscretizedAffineEstimator(Z1, certainty_kde_res_sinc[:kernel])
LinCalibDLV = DiscretizedAffineEstimator(Z1, certainty_kde_res_sinc[:kernel])

max_bias
plot(LinCalib)

max_bias_res = worst_case_bias(LinCalib,
                               Z1,
                               gmix_tmp,
                               MarginalDensityTarget(StandardNormalSample(0.0));
                               maximization=true)
max_bias_res[:max_bias]

using Plots
pgfplotsx()

y = fitted_kde.density
infty_bound = C∞
lower_lims = y .- max.(y .- infty_bound, 0)
upper_lims = fill(infty_bound, length(lower_lims))

plot(fitted_kde.x, fitted_kde.density,
       linestyle=:dash,
       ribbon = (lower_lims, upper_lims),
       fillcolor="red", fillalpha=0.3)
plot!(fitted_kde.x, f_true, color="black", linestyle=:solid)
#----------------------------------------
#---- Test with actual Sinc Kernel-------
#----------------------------------------
f_sinc = sinc_kde(Xs, marginal_grid, SincKernel)

x_tst_1 = f_sinc.x[700]
f_tst_1 = f_sinc.density[700]

# Let us compare to the result from comte_butucea
f_cbt_1 = estimate(Xs, ComteButucea,
  MarginalDensityTarget(x_tst_1), marginal_grid)

@test f_cbt_1 ≈ f_tst_1 atol=0.001

#----- Check default dispatch ---------------------------------
dv_kde = sinc_kde(Xs, marginal_grid, DeLaValleePoussinKernel)
dv_kde_auto_dispatch = sinc_kde(Xs, marginal_grid)
@test dv_kde.density == dv_kde_auto_dispatch.density

#--------------------------------------------------------------
#---- Tests with both Sinc and DeLaValleePoussin Kernels-------
#--- Are our bands longer than true bands?
#--------------------------------------------------------------

f_true = pdf.(Ref(d_true), marginal_grid)

for T in [SincKernel, DeLaValleePoussinKernel]
  f_kde= sinc_kde(Xs, marginal_grid, T)
  true_C = maximum(abs.(f_true .- f_kde.density))

  f_nb = fit(BinnedMarginalDensityNeighborhood, Xs, marginal_grid, T)
  f_nb.C_std

  f_nb_ds = fit(BinnedMarginalDensityNeighborhood, Xs, ds, T)
  @show true_C, f_nb_ds.C_std, f_nb_ds.C_bias, f_nb_ds.C_std + f_nb_ds.C_bias
  @test true_C <= f_nb_ds.C_std + f_nb_ds.C_bias
end