text:
# # Stratified Kelvin-Helmholtz instability
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, JLD2, CairoMakie"
# ```

# ## The physical domain
#
# We simulate Kelvin-Helmholtz instability in two-dimensions in ``x, z``
# and therefore assign `Flat` to the `y` direction,

using Oceananigans

grid = RectilinearGrid(size=(64, 64), x=(-5, 5), z=(-5, 5),
                       topology=(Periodic, Flat, Bounded))

# # The basic state
#
# We're simulating the instability of a sheared and stably-stratified basic state
# ``U(z)`` and ``B(z)``. Two parameters define our basic state: the Richardson number
# defined as,
#
# ```math
# Ri = \frac{∂_z B}{(∂_z U)^2} ,
# ```
#
# and the width of the stratification layer, ``h``.

shear_flow(x, y, z, t) = tanh(z)

stratification(x, y, z, t, p) = p.h * p.Ri * tanh(z / p.h)

U = BackgroundField(shear_flow)

B = BackgroundField(stratification, parameters=(Ri=0.1, h=1/4))

# Our basic state thus has a thin layer of stratification in the center of
# the channel, embedded within a thicker shear layer surrounded by unstratified fluid.

using CairoMakie

zF = znodes(Face, grid)
zC = znodes(Center, grid)

Ri, h = B.parameters

fig = Figure(resolution = (850, 450))

ax = Axis(fig[1, 1], xlabel = "U(z)", ylabel = "z")
lines!(ax, shear_flow.(0, 0, zC, 0), zC; linewidth = 3)

ax = Axis(fig[1, 2], xlabel = "B(z)")
lines!(ax, [stratification(0, 0, z, 0, (Ri=Ri, h=h)) for z in zC], zC; linewidth = 3, color = :red)

ax = Axis(fig[1, 3], xlabel = "Ri(z)")
lines!(ax, [Ri * sech(z / h)^2 / sech(z)^2 for z in zF], zF; linewidth = 3, color = :black) # Ri(z)= ∂_z B / (∂_z U)²; derivatives computed by hand

# In unstable flows it is often useful to determine the dominant spatial structure of the
# instability and the growth rate at which the instability grows.
# If the simulation idealizes a physical flow, this can be used to make
# predictions as to what should develop and how quickly.
# Since these instabilities are often attributed to a linear instability,
# we can determine information about the structure and the growth rate of the instability by analyzing the linear operator
# that governs small perturbations about a base state, or by solving for the linear dynamics.

# Here, we discuss first briefly linear instabilities and how one can obtain growth rates and structures of most unstable
# modes via eigenanalysis. Then we present an alternative method for approximating the eigenanalysis results when
# one does not have access to the linear dynamics or the linear operator about the base state.

# # Linear Instabilities
#
# The base state ``U(z)``, ``B(z)`` is a solution of the inviscid equations of motion. Whether the base state is
# stable or not is determined by whether small perturbations about this base state grow or decay. To formalize this,
# we study the linearized dynamics satisfied by perturbations about the base state:
# ```math
# \partial_t \Phi = L \Phi \, .
# ```
# where ``\Phi = (u, v, w, b)`` is a vector of the perturbation velocities ``u, v, w`` and perturbation buoyancy ``b``
# and ``L`` a linear operator that depends on the base state, ``L = L(U(z), B(z))`` (the `background_fields`).
# Eigenanalysis of the linear operator ``L`` determines the stability of the base state, such as the Kelvin-Helmholtz
# instability. That is, by using the ansantz
# ```math
# \Phi(x, y, z, t) = \phi(x, y, z) \, \exp(\lambda t) \, ,
# ```
# then ``\lambda`` and ``\phi`` are respectively eigenvalues and eigenmodes of ``L``, i.e., they obey
# ```math
# L \, \phi_j = \lambda_j \, \phi_j \quad j=1,2,\dots \, .
# ```
# From hereafter we'll use the convention that the eigenvalues are ordered according to their real part,
# ``\mathrm{Re}(\lambda_1) \ge \mathrm{Re}(\lambda_2) \ge \dotsb``.
#
# Remarks:
#
# As we touched upon briefly above, Oceananigans.jl, does not include a linearized version of the equations.
# Furthermore, Oceananigans.jl does not give us access to the linear operator ``L`` so that we can perform eigenanalysis.
# Below we discuss an alternative way of approximating the eigenanalysis results.
# The method boils down to solving the nonlinear equations while continually renormalize
# the magnitude of the perturbations to ensure that nonlinear terms
# (terms that are quadratic or higher in perturbations) remain negligibly small,
# i.e.,much smaller than the background flow.

# # The power method algorithm
#
# Successive application of ``L`` to a random initial state will eventually render it parallel
# with eigenmode ``\phi_1``:
# ```math
# \lim_{n \to \infty} L^n \Phi \propto \phi_1 \, .
# ```
# Of course, if ``\phi_1`` is an unstable mode (i.e., ``\sigma_1 = \mathrm{Re}(\lambda_1) > 0``), then successive application
# of ``L`` will lead to exponential amplification. (Similarly, if ``\sigma_1 < 0``, successive application of ``L`` will
# lead to exponential decay of ``\Phi`` down to machine precision.) Therefore, after each
# application of the linear operator ``L``, we rescale the output ``L \Phi`` back to a pre-selected amplitude.
#
# So, we initialize a `simulation` with random initial conditions with amplitude much less than those of
# the base state (which are ``O(1)``). Instead of "applying" ``L`` on our initial state, we evolve the
# (approximately) linear dynamics for interval ``\Delta \tau``. We measure how much the energy has grown
# during that interval, rescale the perturbations back to original energy amplitude and repeat.
# After some iterations the state will converge to the most unstable eigenmode.
#
# In summary, each iteration of the power method includes:
# - compute the perturbation energy, ``E_0``,
# - evolve the system for a time-interval ``\Delta \tau``,
# - compute the perturbation energy, ``E_1``,
# - determine the exponential growth of the most unstable mode during the interval
#   ``\Delta \tau`` as  ``\log(E_1 / E_0) / (2 \Delta \tau)``,
# - repeat the above until growth rate converges.
#
# By fiddling a bit with ``\Delta t`` we can get convergence after only a few iterations.
#
# Let's apply all these to our example.

# # The model

model = NonhydrostaticModel(timestepper = :RungeKutta3,
                              advection = UpwindBiasedFifthOrder(),
                                   grid = grid,
                               coriolis = nothing,
                      background_fields = (u=U, b=B),
                                closure = ScalarDiffusivity(ν=2e-4, κ=2e-4),
                               buoyancy = BuoyancyTracer(),
                                tracers = :b)

# We have included a "pinch" of viscosity and diffusivity in anticipation of what will follow furtherdown:
# viscosity and diffusivity will ensure numerical stability when we evolve the unstable mode to the point it becomes nonlinear.

# For this example, we take ``\Delta \tau = 15``.

simulation = Simulation(model, Δt=0.1, stop_iteration=150)

# Now some helper functions that will be used during for the power method algorithm.
#
# First a function that evolves the state for ``\Delta \tau`` and measure the energy growth
# over that period.

"""
    grow_instability!(simulation, energy)

Grow an instability by running `simulation`.

Estimates the growth rate ``σ`` of the instability
using the fractional change in volume-mean kinetic energy,
over the course of the `simulation`

``
energy(t₀ + Δτ) / energy(t₀) ≈ exp(2 σ Δτ)
``

where ``t₀`` is the starting time of the simulation and ``t₀ + Δτ``
the ending time of the simulation. We thus find that the growth rate
is measured by

``
σ = log(energy(t₀ + Δτ) / energy(t₀)) / (2 Δτ) .
``
"""
function grow_instability!(simulation, energy)
    ## Initialize
    simulation.model.clock.iteration = 0
    t₀ = simulation.model.clock.time = 0
    compute!(energy)
    energy₀ = energy[1, 1, 1]

    ## Grow
    run!(simulation)

    ## Analyze
    compute!(energy)
    energy₁ = energy[1, 1, 1]
    Δτ = simulation.model.clock.time - t₀

    ## ½(u² + v²) ~ exp(2 σ Δτ)
    σ = growth_rate = log(energy₁ / energy₀) / 2Δτ

    return growth_rate
end
nothing # hide

# Finally, we write a function that rescales the state. The rescaling is done via computing the
# kinetic energy and then rescaling all flow fields so that the kinetic energy assumes a targetted value.
#
# (Measuring the perturbation growth via the kinetic energy works fine _unless_ an unstable mode _only_ has
# buoyancy structure. In that case, the total perturbation energy is more adequate.)

"""
    rescale!(model, energy; target_kinetic_energy = 1e-3)

Rescales all model fields so that `energy = target_kinetic_energy`.
"""
function rescale!(model, energy; target_kinetic_energy = 1e-6)
    compute!(energy)
    rescale_factor = √(target_kinetic_energy / energy[1, 1, 1])

    for f in merge(model.velocities, model.tracers)
        f .*= rescale_factor
    end

    return nothing
end

using Printf

# Another helper function for the power method,

"""
    convergence(σ)

Check if the growth rate has converged. If the array `σ` has at least 2 elements then returns the
relative difference between ``σ[end]`` and ``σ[end-1]``.
"""
convergence(σ) = length(σ) > 1 ? abs((σ[end] - σ[end-1]) / σ[end]) : 9.1e18 # pretty big (not Inf tho)
nothing # hide

# and the main function that performs the power method iteration.

"""
    estimate_growth_rate(simulation, energy, ω; convergence_criterion=1e-3)

Estimates the growth rate iteratively until the relative change
in the estimated growth rate ``σ`` falls below `convergence_criterion`.

Returns ``σ``.
"""
function estimate_growth_rate(simulation, energy, ω, b; convergence_criterion=1e-3)
    σ = []
    power_method_data = []
    compute!(ω)
    push!(power_method_data, (ω=collect(interior(ω)[:, 1, :]), b=collect(interior(b)[:, 1, :]), σ=deepcopy(σ)))

    while convergence(σ) > convergence_criterion
        compute!(energy)

        @info @sprintf("About to start power method iteration %d; kinetic energy: %.2e", length(σ)+1, energy[1, 1, 1])
        push!(σ, grow_instability!(simulation, energy))
        compute!(energy)

        @info @sprintf("Power method iteration %d, kinetic energy: %.2e, σⁿ: %.2e, relative Δσ: %.2e",
                       length(σ), energy[1, 1, 1], σ[end], convergence(σ))

        compute!(ω)
        rescale!(simulation.model, energy)
        push!(power_method_data, (ω=collect(interior(ω)[:, 1, :]), b=collect(interior(b)[:, 1, :]), σ=deepcopy(σ)))
    end

    return σ, power_method_data
end
nothing # hide

# # Eigenplotting
#
# A good algorithm wouldn't be complete without a good visualization,

u, v, w = model.velocities
b = model.tracers.b

perturbation_vorticity = Field(∂z(u) - ∂x(w))

xω, yω, zω = nodes(perturbation_vorticity)
xb, yb, zb = nodes(b)

# # Rev your engines...
#
# We initialize the power iteration with random noise and rescale to have a `target_kinetic_energy`

using Random, Statistics

mean_perturbation_kinetic_energy = Field(Average(1/2 * (u^2 + w^2)))
noise(x, y, z) = randn()
set!(model, u=noise, w=noise, b=noise)
rescale!(simulation.model, mean_perturbation_kinetic_energy, target_kinetic_energy=1e-6)
growth_rates, power_method_data = estimate_growth_rate(simulation, mean_perturbation_kinetic_energy, perturbation_vorticity, b)

@info "Power iterations converged! Estimated growth rate: $(growth_rates[end])"

# # Powerful convergence
#
# We animate the power method steps. A scatter plot illustrates how the growth rate converges
# as the power method iterates.

n = Observable(1)

fig = Figure(resolution=(800, 600))

kwargs = (xlabel="x", ylabel="z", limits = ((xω[1], xω[end]), (zω[1], zω[end])), aspect=1,)

ω_title(t) = t === nothing ? @sprintf("vorticity") : @sprintf("vorticity at t = %.2f", t)
b_title(t) = t === nothing ? @sprintf("buoyancy")  : @sprintf("buoyancy at t = %.2f", t)

ax_ω = Axis(fig[2, 1]; title = ω_title(nothing), kwargs...)

ax_b = Axis(fig[2, 3]; title = b_title(nothing), kwargs...)

ωₙ = @lift power_method_data[$n].ω
bₙ = @lift power_method_data[$n].b

σₙ = @lift [(i-1, i==1 ? NaN : growth_rates[i-1]) for i in 1:$n]

ω_lims = @lift (-maximum(abs, power_method_data[$n].ω) - 1e-16, maximum(abs, power_method_data[$n].ω) + 1e-16)
b_lims = @lift (-maximum(abs, power_method_data[$n].b) - 1e-16, maximum(abs, power_method_data[$n].b) + 1e-16)

hm_ω = heatmap!(ax_ω, xω, zω, ωₙ; colorrange = ω_lims, colormap = :balance)
Colorbar(fig[2, 2], hm_ω)

hm_b = heatmap!(ax_b, xb, zb, bₙ; colorrange = b_lims, colormap = :balance)
Colorbar(fig[2, 4], hm_b)

eigentitle(σ, t) = length(σ) > 0 ? @sprintf("Iteration #%i; growth rate %.2e", length(σ), σ[end]) : @sprintf("Initial perturbation fields")
σ_title = @lift eigentitle(power_method_data[$n].σ, nothing)

ax_σ = Axis(fig[1, :];
            xlabel = "Power iteration",
            ylabel = "Growth rate",
            title = σ_title,
            xticks = 1:length(power_method_data)-1,
            limits = ((0.5, length(power_method_data)-0.5), (-0.25, 0.25)))

scatter!(ax_σ, σₙ; color = :blue)

frames = 1:length(power_method_data)

record(fig, "powermethod.mp4", frames, framerate=1) do i
       @info "Plotting frame $i of $(frames[end])..."
       n[] = i
end

nothing #hide

# ![](powermethod.mp4)

# # Now for the fun part
#
# Now we simulate the nonlinear evolution of the eigenmode
# we've isolated for a few e-folding times ``1/\sigma``,

## Reset the clock
model.clock.iteration = 0
model.clock.time = 0

estimated_growth_rate = growth_rates[end]

simulation.stop_time = 5 / estimated_growth_rate
simulation.stop_iteration = 9.1e18 # pretty big (not Inf tho)

## Rescale the eigenmode
initial_eigenmode_energy = 5e-5
rescale!(simulation.model, mean_perturbation_kinetic_energy, target_kinetic_energy=initial_eigenmode_energy)

# Let's save and plot the perturbation vorticity and buoyancy and also the total vorticity and
# buoyancy (perturbation + basic state). It'll be also neat to plot the kinetic energy time-series
# and confirm it grows with the estimated growth rate.

total_vorticity = Field(∂z(u) + ∂z(model.background_fields.velocities.u) - ∂x(w))

total_b = Field(b + model.background_fields.tracers.b)

simulation.output_writers[:vorticity] =
    JLD2OutputWriter(model, (ω=perturbation_vorticity, Ω=total_vorticity, b=b, B=total_b, KE=mean_perturbation_kinetic_energy),
                     schedule = TimeInterval(0.10 / estimated_growth_rate),
                     filename = "kelvin_helmholtz_instability.jld2",
                     overwrite_existing = true)

# And now we...

@info "*** Running a simulation of Kelvin-Helmholtz instability..."
run!(simulation)

# ## Pretty things
#
# First we plot the nonlinear equilibration of the perturbation fields together
# with the evolution of the kinetic energy,

@info "Making a neat movie of stratified shear flow..."

filepath = simulation.output_writers[:vorticity].filepath

ω_timeseries = FieldTimeSeries(filepath, "ω")
b_timeseries = FieldTimeSeries(filepath, "b")
Ω_timeseries = FieldTimeSeries(filepath, "Ω")
B_timeseries = FieldTimeSeries(filepath, "B")
KE_timeseries = FieldTimeSeries(filepath, "KE")

times = ω_timeseries.times

t_final = times[end]

n = Observable(1)

ωₙ = @lift interior(ω_timeseries, :, 1, :, $n)
bₙ = @lift interior(b_timeseries, :, 1, :, $n)

fig = Figure(resolution=(800, 600))

kwargs = (xlabel="x", ylabel="z", limits = ((xω[1], xω[end]), (zω[1], zω[end])), aspect=1,)

title = @lift @sprintf("t = %.2f", times[$n])

ax_ω = Axis(fig[2, 1]; title = "perturbation vorticity", kwargs...)

ax_b = Axis(fig[2, 3]; title = "perturbation buoyancy", kwargs...)

ax_KE = Axis(fig[3, :];
             yscale = log10,
             limits = ((0, t_final), (initial_eigenmode_energy, 1e-1)),
             xlabel = "time")

fig[1, :] = Label(fig, title, textsize=24, tellwidth=false)

ω_lims = @lift (-maximum(abs, interior(ω_timeseries, :, 1, :, $n)) - 1e-16, maximum(abs, interior(ω_timeseries, :, 1, :, $n)) + 1e-16)
b_lims = @lift (-maximum(abs, interior(b_timeseries, :, 1, :, $n)) - 1e-16, maximum(abs, interior(b_timeseries, :, 1, :, $n)) + 1e-16)

hm_ω = heatmap!(ax_ω, xω, zω, ωₙ; colorrange = ω_lims, colormap = :balance)
Colorbar(fig[2, 2], hm_ω)

hm_b = heatmap!(ax_b, xb, zb, bₙ; colorrange = b_lims, colormap = :balance)
Colorbar(fig[2, 4], hm_b)

tₙ = @lift times[1:$n]
KEₙ = @lift KE_timeseries[1:$n]

lines!(ax_KE, [0, t_final], @. initial_eigenmode_energy * exp(2 * estimated_growth_rate * [0, t_final]);
       label = "~ exp(2 σ t)",
       linewidth = 2,
       color = :black)

lines!(ax_KE, times, KE_timeseries[:];
       label = "perturbation kinetic energy",
       linewidth = 4, color = :blue, alpha = 0.4)

KE_point = @lift Point2f[(times[$n], KE_timeseries[$n][1, 1, 1])]

scatter!(ax_KE, KE_point;
         marker = :circle, markersize = 16, color = :blue)

frames = 1:length(times)

record(fig, "kelvin_helmholtz_instability_perturbations.mp4", frames, framerate=8) do i
    @info "Plotting frame $i of $(frames[end])..."
    n[] = i
end
nothing #hide

# ![](kelvin_helmholtz_instability_perturbations.mp4)

# And then the same for total vorticity & buoyancy of the fluid.
n = Observable(1)

Ωₙ = @lift interior(Ω_timeseries, :, 1, :, $n)
Bₙ = @lift interior(B_timeseries, :, 1, :, $n)

fig = Figure(resolution=(800, 600))

kwargs = (xlabel="x", ylabel="z", limits = ((xω[1], xω[end]), (zω[1], zω[end])), aspect=1,)

title = @lift @sprintf("t = %.2f", times[$n])

ax_Ω = Axis(fig[2, 1]; title = "total vorticity", kwargs...)

ax_B = Axis(fig[2, 3]; title = "total buoyancy", kwargs...)

ax_KE = Axis(fig[3, :];
             yscale = log10,
             limits = ((0, t_final), (initial_eigenmode_energy, 1e-1)),
             xlabel = "time")

fig[1, :] = Label(fig, title, textsize=24, tellwidth=false)

hm_Ω = heatmap!(ax_Ω, xω, zω, Ωₙ; colorrange = (-1, 1), colormap = :balance)
Colorbar(fig[2, 2], hm_Ω)

hm_B = heatmap!(ax_B, xb, zb, Bₙ; colorrange = (-0.05, 0.05), colormap = :balance)
Colorbar(fig[2, 4], hm_B)

tₙ = @lift times[1:$n]
KEₙ = @lift KE_timeseries[1, 1, 1, 1:$n]

lines!(ax_KE, [0, t_final], @. initial_eigenmode_energy * exp(2 * estimated_growth_rate * [0, t_final]);
       label = "~ exp(2 σ t)",
       linewidth = 2,
       color = :black)

lines!(ax_KE, times, KE_timeseries[:];
       label = "perturbation kinetic energy",
       linewidth = 4, color = :blue, alpha = 0.4)

KE_point = @lift Point2f[(times[$n], KE_timeseries[$n][1, 1, 1])]

scatter!(ax_KE, KE_point;
         marker = :circle, markersize = 16, color = :blue)

axislegend(ax_KE; position = :rb)

record(fig, "kelvin_helmholtz_instability_total.mp4", frames, framerate=8) do i
       @info "Plotting frame $i of $(frames[end])..."
       n[] = i
end
nothing #hide

# ![](kelvin_helmholtz_instability_total.mp4)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus <: SwaggerModel
    evaluationError::Any # spec type: Union{ Nothing, String } # spec name: evaluationError
    incomplete::Any # spec type: Union{ Nothing, Bool } # spec name: incomplete
    nonResourceRules::Any # spec type: Union{ Nothing, Vector{IoK8sApiAuthorizationV1beta1NonResourceRule} } # spec name: nonResourceRules
    resourceRules::Any # spec type: Union{ Nothing, Vector{IoK8sApiAuthorizationV1beta1ResourceRule} } # spec name: resourceRules

    function IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus(;evaluationError=nothing, incomplete=nothing, nonResourceRules=nothing, resourceRules=nothing)
        o = new()
        validate_property(IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus, Symbol("evaluationError"), evaluationError)
        setfield!(o, Symbol("evaluationError"), evaluationError)
        validate_property(IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus, Symbol("incomplete"), incomplete)
        setfield!(o, Symbol("incomplete"), incomplete)
        validate_property(IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus, Symbol("nonResourceRules"), nonResourceRules)
        setfield!(o, Symbol("nonResourceRules"), nonResourceRules)
        validate_property(IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus, Symbol("resourceRules"), resourceRules)
        setfield!(o, Symbol("resourceRules"), resourceRules)
        o
    end
end # type IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus

const _property_map_IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus = Dict{Symbol,Symbol}(Symbol("evaluationError")=>Symbol("evaluationError"), Symbol("incomplete")=>Symbol("incomplete"), Symbol("nonResourceRules")=>Symbol("nonResourceRules"), Symbol("resourceRules")=>Symbol("resourceRules"))
const _property_types_IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus = Dict{Symbol,String}(Symbol("evaluationError")=>"String", Symbol("incomplete")=>"Bool", Symbol("nonResourceRules")=>"Vector{IoK8sApiAuthorizationV1beta1NonResourceRule}", Symbol("resourceRules")=>"Vector{IoK8sApiAuthorizationV1beta1ResourceRule}")
Base.propertynames(::Type{ IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus }) = collect(keys(_property_map_IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus))
Swagger.property_type(::Type{ IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus }, property_name::Symbol) =  _property_map_IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus[property_name]

function check_required(o::IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus)
    (getproperty(o, Symbol("incomplete")) === nothing) && (return false)
    (getproperty(o, Symbol("nonResourceRules")) === nothing) && (return false)
    (getproperty(o, Symbol("resourceRules")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAuthorizationV1beta1SubjectRulesReviewStatus }, name::Symbol, val)
end

ext: jl


text:
using Statistics
using FFTW
using Plots
using BenchmarkTools
using Profile
using LinearAlgebra
using Measures
using ProgressMeter
push!(LOAD_PATH, pwd()*"/main")
using DHC_2DUtils

function DHC_compute_temp(image::Array{Float64,2}, filter_hash::Dict, filter_hash2::Dict=filter_hash;
    doS2::Bool=true, doS20::Bool=false, norm=true, iso=false, FFTthreads=2, normS1::Bool=false, normS1iso::Bool=false)
    # image        - input for WST
    # filter_hash  - filter hash from fink_filter_hash
    # filter_hash2 - filters for second order.  Default to same as first order.
    # doS2         - compute S2 coeffs
    # doS20        - compute S2 coeffs
    # norm         - scale to mean zero, unit variance
    # iso          - sum over angles to obtain isotropic coeffs

    # Use 2 threads for FFT
    FFTW.set_num_threads(FFTthreads)

    # array sizes
    (Nx, Ny)  = size(image)
    if Nx != Ny error("Input image must be square") end
    (Nf, )    = size(filter_hash["filt_index"])
    if Nf == 0  error("filter hash corrupted") end
    @assert Nx==filter_hash["npix"] "Filter size should match npix"
    @assert Nx==filter_hash2["npix"] "Filter2 size should match npix"
    @assert (normS1 && normS1iso) != 1 "normS1 and normS1iso are mutually exclusive"

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, Nf)
    if doS2  S2  = zeros(Float64, Nf, Nf) end  # traditional 2nd order
    if doS20 S20 = zeros(Float64, Nf, Nf) end  # real space correlation
    anyM2 = doS2 | doS20
    anyrd = doS2 | doS20             # compute real domain with iFFT

    # allocate image arrays for internal use
    if anyrd im_rd_0_1  = Array{Float64, 3}(undef, Nx, Ny, Nf) end

    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(norm_im .* norm_im)/(Nx*Ny)
    if norm
        norm_im ./= sqrt(Nx*Ny*S0[2])
    else
        norm_im = copy(image)
    end

    append!(out_coeff,S0[:])

    ## 1st Order
    im_fd_0 = fft(norm_im)  # total power=1.0

    # unpack filter_hash
    f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
    f_val   = filter_hash["filt_value"]

    zarr = zeros(ComplexF64, Nx, Ny)  # temporary array to fill with zvals

    # make a FFTW "plan" for an array of the given size and type
    if anyrd
        P = plan_ifft(im_fd_0) end  # P is an operator, P*im is ifft(im)

    ## Main 1st Order and Precompute 2nd Order
    for f = 1:Nf
        S1tot = 0.0
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i
        # for (ind, val) in zip(f_i, f_v)   # this is slower!
        if length(f_i) > 0
            for i = 1:length(f_i)
                ind       = f_i[i]
                zval      = f_v[i] * im_fd_0[ind]
                S1tot    += abs2(zval)
                zarr[ind] = zval        # filter*image in Fourier domain
            end
            S1[f] = S1tot/(Nx*Ny)  # image power
            if anyrd
                im_rd_0_1[:,:,f] .= abs2.(P*zarr) end
            zarr[f_i] .= 0
        end
    end

    append!(out_coeff, iso ? filter_hash["S1_iso_mat"]*S1 : S1)

    if normS1iso
        S1iso = vec(reshape(filter_hash["S1_iso_mat"]*S1,(1,:))*filter_hash["S1_iso_mat"])
    end

    # we stored the abs()^2, so take sqrt (this is faster to do all at once)
    if anyrd im_rd_0_1 .= sqrt.(im_rd_0_1) end

    Mat2 = filter_hash["S2_iso_mat"]
    if doS2
        f_ind2   = filter_hash2["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val2   = filter_hash2["filt_value"]

        ## Traditional second order
        for f1 = 1:Nf
            thisim = fft(im_rd_0_1[:,:,f1])  # Could try rfft here
            # Loop over f2 and do second-order convolution
            if normS1
                normS1pwr = S1[f1]
            elseif normS1iso
                normS1pwr = S1iso[f1]
            else
                normS1pwr = 1
            end
            for f2 = 1:Nf
                f_i = f_ind2[f2]  # CartesianIndex list for filter
                f_v = f_val2[f2]  # Values for f_i
                # sum im^2 = sum(|fft|^2/npix)
                S2[f1,f2] = sum(abs2.(f_v .* thisim[f_i]))/(Nx*Ny)/normS1pwr
            end
        end
        append!(out_coeff, iso ? Mat2*S2[:] : S2[:])
    end

    # Real domain 2nd order
    if doS20
        Amat = reshape(im_rd_0_1, Nx*Ny, Nf)
        S20  = Amat' * Amat
        append!(out_coeff, iso ? Mat2*S20[:] : S20[:])
    end

    return out_coeff
end

filter_hash = fink_filter_hash(1,8,nx=256,wd=2)

DHC_compute_temp(rand(256,256), filter_hash, iso=false, normS1=true)

DHC_compute_temp(rand(256,256), filter_hash, iso=false, normS1iso=true)

DHC_compute_temp(rand(256,256), filter_hash, iso=false, normS1iso=true, normS1=true)

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
# Collection of basic types to use during routing for clarity.

"""
Type to access channels in the routing taskgraph. Essentially, it is just a
wrapper for an integer, but typed to allow safer and clearer usage.
"""
struct ChannelIndex
    idx::Int
end

Base.getindex(A::AbstractArray, i::ChannelIndex) = getindex(A, i.idx)
Base.setindex!(A::AbstractArray, X, i::ChannelIndex) = setindex!(A, X, i.idx)

"""
Indices of a [`RoutingGraph`](@ref) that can serve as either start or stop
vertices (depending on the context) of one branch of a [`RoutingChannel`](@ref).
"""
struct PortVertices
    indices::Vector{Int64}
end

Base.getindex(A::PortVertices, i) = getindex(A.indices, i)
Base.setindex!(A::PortVertices, X, i) = setindex!(A.incides, X, i)
Base.iterate(A::PortVertices, args...) = iterate(A.indices, args...)
Base.length(A::PortVertices) = length(A.indices)


text:
using Routino
using Test

@testset "Routino.jl" begin
    if isfile(Routino.BINARY)
        wayp1 = (lo=-1.7587022270945216, la=53.67050050773988) # HD2
        wayp2 = (lo=-1.6091688317085435, la=53.74083763716689) # pcode_centre(idx, "LS27 0AL")
        @test distance(wayp1, wayp2, Routino.BINARY) == (km = 19.3, mins = 13)
    end
    if isfile(Routino.LIB)
        router = Router("foot")
        add_waypoint!(-1.7587022270945216, 53.67050050773988, router)
        add_waypoint!(-1.6091688317085435, 53.74083763716689, router)
        out = quickest_route(router)
        @test out != C_NULL
        @test typeof(out) == Ptr{Routino.COutput}
        @test Routino.walk_to_distance(out) == (km = 15, mins = 229)
    end

end


text:
using DashBootstrapComponents

text_input = html_div([
    dbc_label("Text"),
    dbc_input(placeholder = "Input goes here...", type = "text"),
    dbc_formtext("Type something in the box above"),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Miscellaneous functions used by all back-ends.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export compact_type_str

"""
    compact_type_str(T)

Return a string with a compact representation of type `T`.
"""
compact_type_str(T) = string(T)

function compact_type_str(T::Union)
    if VERSION < v"1.3.0"
        str = if T >: Missing
            string(Core.Compiler.typesubtract(T, Missing)) * "?"
        else
            string(T)
        end
    else
        str = T >: Missing ? string(nonmissingtype(T)) * "?" : string(T)
    end

    str = replace(str, "Union" => "U")
    return str
end

################################################################################
#                              Private functions
################################################################################

"""
    _aprint(buf::IO, [v,] indentation = 0, nspace = 2)

Print the variable `v` to the buffer `buf` at the indentation level
`indentation`. Each level has `nspaces` spaces.

If `v` is not present, then only the indentation spaces will be printed.
"""
function _aprint(buf::IO, v::String, indentation::Int = 0, nspaces::Int = 2)
    tokens  = split(v, '\n')
    padding = " "^(indentation * nspaces)
    ntokens = length(tokens)

    @inbounds for i = 1:ntokens
        # If the token is empty, then we do nothing to avoid unnecessary white
        # spaces.
        if length(tokens[i]) != 0
            print(buf, padding)
            print(buf, tokens[i])
        end
        i != ntokens && println(buf)
    end

    return nothing
end

function _aprint(buf::IO, indentation::Int = 0, nspaces::Int = 2)
    padding = " "^(indentation*nspaces)
    print(buf, padding)

    return nothing
end

"""
    _aprintln(buf::IO, [v,] indentation = 0, nspaces = 2)

Same as `_aprint`, but a new line will be added at the end.
"""
function _aprintln(buf::IO, v::String, indentation::Int = 0, nspaces::Int = 2)
    tokens  = split(v, '\n')
    padding = " "^(indentation * nspaces)
    ntokens = length(tokens)

    @inbounds for i = 1:ntokens
        # If the token is empty, then we do nothing to avoid unnecessary white
        # spaces.
        if length(tokens[i]) != 0
            print(buf, padding)
            println(buf, tokens[i])
        else
            println(buf)
        end
    end

    return nothing
end

function _aprintln(buf::IO, indentation::Int = 0, nspaces::Int = 2)
    padding = " "^(indentation*nspaces)
    println(buf, padding)

    return nothing
end

"""
    _process_hlines(hlines::Union{Symbol,AbstractVector}, body_hlines::AbstractVector, num_printed_rows::Int, noheader::Bool)

Process the horizontal lines in `hlines` and `body_hlines` considering the
number of printed rows `num_printed_rows` and if the header is present
(`noheader`).

It returns a vector of `Int` stating where the horizontal lines must be drawn.
"""
@inline function _process_hlines(
    hlines::Symbol,
    body_hlines::AbstractVector,
    num_printed_rows::Int,
    noheader::Bool
)
    if hlines == :all
        vhlines = collect(0:1:num_printed_rows + !noheader)
    elseif hlines == :none
        vhlines = Int[]
    else
        error("`hlines` must be `:all`, `:none`, or a vector of integers.")
    end

    return _process_hlines(vhlines, body_hlines, num_printed_rows, noheader)
end

function _process_hlines(
    hlines::AbstractVector,
    body_hlines::AbstractVector,
    num_printed_rows::Int,
    noheader::Bool
)
    # The symbol `:begin` is replaced by 0, the symbol `:header` by the line
    # after the header, and the symbol `:end` is replaced by the last row.
    hlines = replace(
        hlines,
        :begin  => 0,
        :header => noheader ? -1 : 1,
        :end    => num_printed_rows + !noheader
    )

    # All numbers less than 1 and higher or equal the number of printed rows
    # must be # removed from `body_hlines`.
    body_hlines = filter(x -> (x ≥ 1) && (x < num_printed_rows), body_hlines)

    # Merge `hlines` with `body_hlines`.
    hlines = unique(vcat(hlines, body_hlines .+ !noheader))
    #                                               ^
    #                                               |
    # If we have header, then the index in `body_hlines` must be incremented.

    # Make sure that the compiler knows that this function always returns a
    # vector of `Int`s.
    ret::Vector{Int} = Vector{Int}(hlines)

    return ret
end

"""
    _process_vlines(vlines::AbstractVector, num_printed_cols::Int)

Process the vertical lines `vlines` considerering the number of printed columns
`num_printed_cols`.

It returns a vector of `Int` stating where the vertical lines must be drawn.
"""
@inline function _process_vlines(vlines::Symbol, num_printed_cols::Int)
    # Process `vlines`.
    if vlines == :all
        vvlines = collect(0:1:num_printed_cols)
    elseif vlines == :none
        vvlines = Int[]
    else
        error("`vlines` must be `:all`, `:none`, or a vector of integers.")
    end

    return _process_vlines(vvlines, num_printed_cols)
end

function _process_vlines(vlines::AbstractVector, num_printed_cols::Int)
    # Process `vlines`.
    if vlines == :all
        vlines = collect(0:1:num_printed_cols)
    elseif vlines == :none
        vlines = Int[]
    elseif !(typeof(vlines) <: AbstractVector)
        error("`vlines` must be `:all`, `:none`, or a vector of integers.")
    end

    # The symbol `:begin` is replaced by 0 and the symbol `:end` is replaced by
    # the last column.
    vlines = replace(
        vlines,
        :begin => 0,
        :end   => num_printed_cols
    )

    return Vector{Int}(vlines)
end

"""
    _str_escaped(str::AbstractString)

Return the escaped string representation of `str`.

"""
function _str_escaped(str::AbstractString)
    # NOTE: Here we cannot use `escape_string(str)` because it also adds the
    # character `"` to the list of characters to be escaped.
    return sprint(escape_string, str, "", sizehint = lastindex(str))
end


text:
# flips bits between ith and jth bit in n
flipbits(n, i, j) = n ⊻ (2^(j - i + 1) - 1) << (i - 1)

# checks if given pattern has same number of home and away games
isPattern(nslots, pat) =
    count_ones(pat) == count_zeros(pat) - sizeof(pat) * 8 + nslots

ishard(c::T) where {T<:Constraint} = c.pen == 0

@inline teams(c::Constraint) = c.teams
@inline teams(c::CAcon) = IdxSet(c.teams1.binrepr | c.teams2.binrepr)
@inline teams(c::GAcon) = unique(Iterators.flatten(c.meetings))

@inline patconteam(c::CAcon) = only(c.teams1) # get team for a patetrn con
@inline patconteam(c::BRcon) = only(c.teams)

@inline patsetconteams(c::CAcon) = c.teams1
@inline patsetconteams(c::BRcon) = c.teams

@inline slots(c::Constraint) = c.slots

@inline constraints(i::Instance) =
    vcat(i.constrCA, i.constrBR, i.constrGA, i.constrFA, i.constrSE)

# base case
@inline is_pattern_constraint(rr2::RR2, con::Constraint) = false

@inline is_pattern_constraint(rr2::RR2, con::BRcon) =
    is_pattern_constraint(con::BRcon)
@inline is_pattern_constraint(con::BRcon) = length(con.teams) == 1
@inline is_pattern_constraint(rr2::RR2, con::CAcon) =
    length(con.teams1) == 1 && length(con.teams2) == rr2.nteams

@inline is_pattern_set_constraint(rr2::RR2, con::Constraint) = false
@inline is_pattern_set_constraint(rr2::RR2, con::CAcon) =
    length(con.teams1) > 1 && length(con.teams2) == rr2.nteams
@inline is_pattern_set_constraint(rr2::RR2, con::BRcon) = length(con.teams) > 1

function patternset(sched::Matrix{Int})
    return [
        sum(2^(s - 1) for s in 1:size(sched, 2) if sign(sched[i, s]) == 1) for
        i = 1:size(sched, 1)
    ]
end

"""
    is_schedule(rr2::RR2, S::Array{Int, 2})

Checks if a schedule satisfies all base constraints.

# Arguments

- `rr2::RR2`: Format data of the problem instance
- `S::Array{Int,2}`: The potential schedule to be verified. Assuming teams `i`
    and `j` meet in slot `s` we have `S[i,s] == j`, if team `i` plays at home
    and `S[i,s] == -j` if team `i` plays away.

# Notes

This function is purely intended for debugging purposes.
"""
function is_schedule(rr2::RR2, S::Array{Int,2})
    @assert size(S, 1) == rr2.nteams && size(S, 2) == rr2.nslots

    for s = 1:rr2.nslots
        # exactly rr2.nteams÷2 teams play at home in each slot
        if count(x -> sign(x) == 1, S[:, s]) != rr2.nteams ÷ 2
            return false
        end
        # each team participates in each slot exactly once
        if sort(abs.(S[:, s])) != 1:rr2.nteams
            return false
        end
    end

    # each team plays against each other team once at home and once away
    for r = 1:rr2.nteams
        if sort(S[r, :]) !=
           vcat(-rr2.nteams:-r-1, -r+1:-1, 1:r-1, r+1:rr2.nteams)
            return false
        end
    end
    # check if first half of tournament is a single round robin
    if rr2.isphased
        for i = 1:rr2.nteams
            if sort(abs.(S[i, 1:rr2.nslots÷2])) != vcat(1:i-1, i+1:rr2.nteams)
                return false
            end
        end
    end

    return true
end


text:
struct API
    mod::Module
    var::Symbol
end

Base.fullname(api::API) = (fullname(api.mod)..., api.var)
Base.Module(api::API) = api.mod
Base.nameof(api::API) = api.var

"""
    PublicAPI.of(provider::Module; [recursive = true]) -> apis::Vector

List public API from the `provider` module.

Each element `api` of `apis` supports the following accessor functions:

* `Module(api) :: Module`: module in which the API is defined
* `nameof(api) :: Symbol`: the name of the API in the module
* `fullname(api) :: Tuple{Vararg{Symbol}}`: the components of the
  fully-qualified name; i.e., `(:Package, :SubModule, :function)` for
  `Package.SubModule.function`.

The `provider` module itself is not included in the `apis`.

# Keyword Arguments
- `recursive::Bool = true`: Include public APIs from public sub-modules.
"""
function PublicAPI.of(provider::Module; recursive::Bool = true)
    apis = Vector{API}()
    function sweep(m::Module)
        function mayberecurse(n::Symbol)
            if recursive
                y = try
                    getfield(m, n)
                catch
                    nothing
                end
                if y isa Module && parentmodule(y) === m && y !== m
                    sweep(y)
                end
            end
        end
        registry = nothing
        if isdefined(m, API_REGISTRY_NAME)
            registry = getfield(m, API_REGISTRY_NAME)
            if registry isa Vector{Symbol}
                for n in registry
                    push!(apis, API(m, n))
                    mayberecurse(n)
                end
            else
                @error "Malformed API registry found in module `$m`" registry
            end
        end
        for n in names(m)
            if registry !== nothing
                n in registry && continue
            end
            n === nameof(m) && continue  # avoid the module to be listed twice
            if isdefined(m, n)
                push!(apis, API(m, n))
                mayberecurse(n)
            else
                # Should error?
                @error "Exported but undefined: `$m.$n`"
            end
        end
    end
    sweep(provider)
    return apis
end


text:
# module for defining world coordinate system
# puddle can be put

using Plots
pyplot()

include(joinpath(split(@__FILE__, "src")[1], "src/model/puddle/puddle.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/model/goal/goal.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/model/robot/differential_wheeled_robot/differential_wheeled_robot.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/model/robot/warp_robot/warp_robot.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/model/agent/agent.jl"))

mutable struct PuddleWorld
  x_min
  x_max
  y_min
  y_max
  objects
  delta_time
  end_time
  is_test
  save_path
  puddles
  robots
  goals

  # init
  function PuddleWorld(x_min::Float64, x_max::Float64,
                       y_min::Float64, y_max::Float64;
                       delta_time=0.1, end_time=30,
                       is_test=false, save_path=nothing)
    self = new()
    self.x_min = x_min
    self.x_max = x_max
    self.y_min = y_min
    self.y_max = y_max
    self.objects = []
    self.delta_time = delta_time
    self.end_time = end_time
    self.is_test = is_test
    self.save_path = save_path
    self.puddles = []
    self.robots = []
    self.goals = []
    return self
  end
end

function append(self::PuddleWorld, obj)
  push!(self.objects, obj)

  if typeof(obj) == Puddle
    push!(self.puddles, obj)
  end

  if typeof(obj) == DifferentialWheeledRobot || typeof(obj) == WarpRobot
    push!(self.robots, obj)
  end

  if typeof(obj) == Goal
    push!(self.goals, obj)
  end
end

function puddle_depth(self::PuddleWorld, pose)
  if length(self.puddles) > 0
    return sum([p.depth*inside(p, pose)  for p in self.puddles])
  else
    return 0.0
  end
end

function one_step(self::PuddleWorld, delta_time)
  plot([], [], aspect_ratio=true, xlabel="X", ylabel="Y",
       xlims=(self.x_min, self.x_max), ylims=(self.y_min, self.y_max),
       legend=false)

  for r in self.robots
    if typeof(r.agent) != Agent
      r.agent.puddle_depth = puddle_depth(self, r.pose)
      for g in self.goals
        if inside(g, r.pose)
          r.agent.in_goal = true
          r.agent.final_value = g.value
        end
      end
    end
  end

  for obj in self.objects
    draw!(obj)
  end

  annotate!(-3.5, 4.5, "t = $(delta_time)", "black")
end

function draw(self::PuddleWorld)
  if self.is_test
    for t in 0:self.delta_time:10
      one_step(self, t)
    end
  else
    anime = @animate for t in 0:self.delta_time:self.end_time
      one_step(self, t)
    end
    gif(anime, fps=15, joinpath(split(@__FILE__, "src")[1], self.save_path))
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
module TestLookup

include("preamble.jl")

using SymSpellChecker: SuggestItem, Verbosity,
    VerbosityALL, VerbosityTOP, VerbosityCLOSEST

@testset "basic lookup" begin
    d = SymSpell()

    push!(d, "a", 5)
    push!(d, "bc", 2)

    @test_throws ArgumentError lookup(d, "xyz", max_edit_distance = 3)

    @test isempty(lookup(d, "qwerty"))
    @test lookup(d, "qwerty", include_unknown = true) == [SuggestItem("qwerty", 3, 0)]
    @test lookup(d, "bc") == [SuggestItem("bc", 0, 2)]
    @test lookup(d, "qwe", max_edit_distance = 0, include_unknown = true) == [SuggestItem("qwe", 1, 0)]
    @test lookup(d, "abc", max_edit_distance = 1) == [SuggestItem("bc", 1, 2)]

    @test first(d["abc"]) == "bc"
end

@testset "basic options" begin
    d = SymSpell()

    push!(d, "a", 5)
    push!(d, "bc", 2)

    set_options!(d, include_unknown = true)
    @test first(d["qwerty"]) == "qwerty"

    set_options!(d, include_unknown = false, transfer_casing = true, max_edit_distance = 2)
    @test first(d["Bde"]) == "Bc"

    @test length(lookup(d, "ac", verbosity = "all")) == 2
end

@testset "words with shared prefix should retain counts" begin
    d = SymSpell()
    push!(d, "pipe", 5)
    push!(d, "pips", 10)

    result = lookup(d, "pipe", verbosity = VerbosityALL, max_edit_distance = 1)
    @test length(result) == 2
    @test result[1].term == "pipe"
    @test result[1].count == 5
    @test result[2].term == "pips"
    @test result[2].count == 10

    result = lookup(d, "pips", verbosity = VerbosityALL, max_edit_distance = 1)
    @test length(result) == 2
    @test result[1].term == "pips"
    @test result[1].count == 10
    @test result[2].term == "pipe"
    @test result[2].count == 5

    result = lookup(d, "pip", verbosity = VerbosityALL, max_edit_distance = 1)
    @test length(result) == 2
    @test result[1].term == "pips"
    @test result[1].count == 10
    @test result[2].term == "pipe"
    @test result[2].count == 5
end

@testset "verbosity should control lookup results" begin
    d = SymSpell()
    push!(d, "steam", 1)
    push!(d, "steams", 2)
    push!(d, "steem", 3)

    result = lookup(d, "steems", verbosity = VerbosityTOP, max_edit_distance = 2)
    @test length(result) == 1

    result = lookup(d, "steems", verbosity = VerbosityCLOSEST, max_edit_distance = 2)
    @test length(result) == 2

    result = lookup(d, "steems", verbosity = VerbosityALL, max_edit_distance = 2)
    @test length(result) == 3
end

@testset "lookup should return most frequent" begin
    d = SymSpell()
    push!(d, "steama", 4)
    push!(d, "steamb", 6)
    push!(d, "steamc", 2)

    result = lookup(d, "stream", verbosity = VerbosityTOP, max_edit_distance = 2)
    @test length(result) == 1
    @test result[1].term == "steamb"
    @test result[1].count == 6
end

@testset "lookup should find exact match" begin
    d = SymSpell()
    push!(d, "steama", 4)
    push!(d, "steamb", 6)
    push!(d, "steamc", 2)

    result = lookup(d, "streama", verbosity = VerbosityTOP, max_edit_distance = 2)
    @test length(result) == 1
    @test result[1].term == "steama"
end

@testset "lookup should not return non word delete" begin
    d = SymSpell(count_threshold = 10)
    push!(d, "pawn", 10)
    result = lookup(d, "paw", verbosity = VerbosityTOP, max_edit_distance = 0)
    @test length(result) == 0

    result = lookup(d, "awn", verbosity = VerbosityTOP, max_edit_distance = 0)
    @test length(result) == 0
end

@testset "lookup should not return low count word" begin
    d = SymSpell(count_threshold = 10)
    push!(d, "pawn", 1)
    result = lookup(d, "pawn", verbosity = VerbosityTOP, max_edit_distance = 0)
    @test length(result) == 0
end

@testset "lookup should not return low count word that are also delete word" begin
    d = SymSpell(count_threshold = 10)
    push!(d, "flame", 20)
    push!(d, "flam", 1)
    result = lookup(d, "flam", verbosity = VerbosityTOP, max_edit_distance = 0)
    @test length(result) == 0
end

@testset "lookup max edit distance too large" begin
    d = SymSpell(count_threshold = 10)
    push!(d, "flame", 20)
    push!(d, "flam", 1)
    @test_throws ArgumentError lookup(d, "flam", verbosity = VerbosityTOP, max_edit_distance = 3)
end

@testset "lookup include unknown" begin
    d = SymSpell(count_threshold = 10)
    push!(d, "flame", 20)
    push!(d, "flam", 1)
    result = lookup(d, "flam", verbosity = VerbosityTOP, max_edit_distance = 0, include_unknown = true)
    @test length(result) == 1
    @test result[1].term == "flam"
end

@testset "lookup avoid exact match early exit" begin
    max_edit_distance = 2
    d = SymSpell(count_threshold = 10, max_dictionary_edit_distance = max_edit_distance)
    push!(d, "flame", 20)
    push!(d, "flam", 1)
    result = lookup(d, "24th", verbosity = VerbosityALL, max_edit_distance = max_edit_distance,
        ignore_token = r"\d{2}\w*\b")
    @test length(result) == 1
    @test result[1].term == "24th"
end

@testset "no duplicates in lookup" begin
    d = SymSpell()
    push!(d, "bank", 10)
    result = lookup(d, "xbank", verbosity = VerbosityALL)

    @test length(result) == 1
end

@testset "shouldn't switch to far away words" begin
    d = SymSpell()
    push!(d, "border", 200)
    push!(d, "bored", 40)
    push!(d, "board", 50)
    push!(d, "for", 100)
    result = lookup(d, "bord")

    @test result[1].term == "board"
end

@testset "should correctly check words at the prefix edge" begin
    d = SymSpell(max_dictionary_edit_distance = 1, prefix_length = 5)
    push!(d, "called", 10)

    result = lookup(d, "callvd", verbosity = VerbosityTOP)
    @test length(result) > 0
    @test result[1].term == "called"

    d = SymSpell(max_dictionary_edit_distance = 1, prefix_length = 5)
    push!(d, "calleds", 10)
    result = lookup(d, "callvds", verbosity = VerbosityTOP)
    @test length(result) > 0
    @test result[1].term == "calleds"

    d = SymSpell(max_dictionary_edit_distance = 1, prefix_length = 5)
    push!(d, "calleds", 10)
    result = lookup(d, "callvdx", verbosity = VerbosityTOP)
    @test length(result) == 0
end

@testset "should correctly search utf keywords" begin
    d = SymSpell(max_dictionary_edit_distance = 1)

    push!(d, "привет", 10)
    @test length(d["прифет"]) > 0
    @test length(d["пфетик"]) == 0
end

end # module


text:
"""
Functions to solve the 1d1d scrap - off-gas model
Here we use an explicit euler (EE), which is not very elegant or fast, but allows
easy modifications atm - and should deliver accurate (physical) solution if stable ...

DISCLAIMER:
The newer solver_1d1d_diffeq solver should be used, which uses DifferentialEquations.jl
with the Tsit5() solver and should be generally faster and more robust.
This solver is only left here for development and comparision.

Be careful when using it, there may be instabillities and some tweaking of timestep
sizes may be necessary to get correct results ...
"""

function solve1d1d_ee(m::M1d1d)

    print_log(m1d1d_logfile, "Using explicit euler solver!")
    # Init Results according to used meshes
    results = Vector{Results1d1d}(undef,0)

    mesh_tmp = nothing

    T_H = collect(0:0.05:Tmax)

    scrap_conduction = false
    scrap_conduction_temp = false
    phase_start_time = 0.0

    mesh = M1d1dMesh(m, 1)
    T = M1d1dVals(mesh, init_val = m.bulk[1].T_charged)
    Tnew = M1d1dVals(mesh, init_val = m.bulk[1].T_charged)

    for p=1:1:length(m.sim_phase)
        t_save = m.sim_phase[p].dt_save

        mesh = M1d1dMesh(m, p)
        printInfo(mesh; tabs = 1)

        # just make sure that meshes do not change between the phases
        # because we than would have to interpolate between the meshes if they
        # change:
        if p > 1
            if mesh_tmp != mesh
                error("Mesh change between different phases is not supported atm")
                print_log(m1d1d_logfile, "Error: Unsupported mesh change")
            else
                mesh_tmp = mesh
            end
        else
            mesh_tmp = mesh
        end

        no_cells = mesh.no_cells
        dx_f = mesh.solid_cells.len

        A_f = mesh.gas_cells.dA
        A_s = mesh.solid_cells.dA

        V_f = mesh.gas_cells.vol
        V_s = mesh.solid_cells.vol
        V_fs = V_f + V_s

        if !isnothing(mesh.solid_1d_submesh)
            no_sub_cells = mesh.solid_1d_submesh.no_cells
            dx_s = mesh.solid_1d_submesh.sub_cells.len
            A_s_s =  mesh.solid_1d_submesh.sub_cells.dA
            V_s_sub_c = mesh.solid_1d_submesh.sub_cells.vol
        end

        heat_flow_exchng_sf = zeros(Float64, m.mesh_settings.no_cells_bulk)

        println("Starting Simulation of Phase $p")

        mf_f = m.offgas[p].mass_flow
        mf_s = m.bulk[p].mass_flow
        # TODO Check solid and fluid are equal

        println("Mass flow fluid: ", mf_f)
        println("Mass flow solid: ", mf_s)

        @assert isapprox(V_fs, dx * m.bulk[p].cross_sect_surf_area)
        @assert  m.bulk[p].scrap_pieces_contact_factor <= 1.0 +
                eps(typeof(m.bulk[p].scrap_pieces_contact_factor))

        # Fluid to Solid Surface (and vice versa)
        no_scrap_pieces_per_vol::Float64 = ( m.bulk[p].density /
                (m.bulk[p].scrap_props.vol * pX_bounded(
                m.bulk[p].scrap_props.thermal_props.rho, m.bulk[p].T_charged)) )
        surface_of_scrap_per_vol::Float64 = (
                m.bulk[p].scrap_pieces_contact_factor * no_scrap_pieces_per_vol
                * m.bulk[p].scrap_props.surf_area )

        A_fs =  surface_of_scrap_per_vol * V_fs

        lam_f =  m.offgas[p].thermal_props.lambda
        rho_f = m.offgas[p].thermal_props.rho
        cp_f = m.offgas[p].thermal_props.cp_true

        lam_s = m.bulk[p].scrap_props.thermal_props.lambda
        rho_s = m.bulk[p].scrap_props.thermal_props.rho
        cp_s = m.bulk[p].scrap_props.thermal_props.cp_true

        # working but actually slows down computation ...
        # println("Tabulation fluid material properties")
        # lam_f = TabulatedMaterialPropertyOfX(m.offgas[p].thermal_props.lambda, Trange)
        # rho_f = TabulatedMaterialPropertyOfX(m.offgas[p].thermal_props.rho, Trange)
        # cp_f = TabulatedMaterialPropertyOfX(m.offgas[p].thermal_props.cp_true, Trange)

        # lam_s = TabulatedMaterialPropertyOfX(m.bulk[p].scrap_props.thermal_props.lambda, Trange)
        # rho_s = TabulatedMaterialPropertyOfX(m.bulk[p].scrap_props.thermal_props.rho, Trange)
        # cp_s = TabulatedMaterialPropertyOfX(m.bulk[p].scrap_props.thermal_props.cp_true, Trange)

        alpha_fs = m.heat_transfer[p]

        println("Integrating H of T curve")
        H_T_f = H_of_T.(cp_f, T_H)
        H_T_s = H_of_T.(cp_s, T_H)

        # just to be sure ...
        check_monolithic_growing(H_T_f)
        check_monolithic_growing(H_T_s)

        println("Check Radiation Settings ...")

        if m.sim_phase[p].submodels.scrap_discretization_model_enabled
            println("Scrap conduction enabled")

            if isnothing(mesh.solid_1d_submesh)
                error("No submesh generated for scrap part heating!")
            end

            scrap_conduction = true
            # 2.0 due to symmetry
            mf_sub_s::Float64  = mf_s / (mesh.solid_1d_submesh.no_cells * 2.0)
        else
            scrap_conduction = false
        end

        if p > 1 && scrap_conduction && !scrap_conduction_temp
            println("Interpolating scrap temperatures to sub cells")
            # set prev temperature of T.s[i] to all T.s_sub[i][j]
            @inbounds for i = 1:no_cells
                for j = 1:no_sub_cells
                    T.s_sub[i][j] = T.s[i]
                end
            end
        end

        if p > 1
            scrap_conduction_temp = scrap_conduction
            phase_start_time += m.sim_phase[p-1].duration
        end

        cross_radiation = m.sim_phase[p].submodels.cross_radiation_model_enabled

        if cross_radiation
            println("Cross radiation enabled")
        end

        front_radiation = !isnothing(
                m.sim_phase[p].submodels.furnace_side_rad_model)

        if front_radiation
            println("Front radiation enabled")
        end

        if cross_radiation || front_radiation

            c = m.sim_phase[p].submodels.cross_radiation_constant_c

            @assert c <= (1.0 + eps(typeof(c)))  && c >= (0.0 - eps(typeof(c)))

            A_fs_rad = c * 0.5 * A_fs

            A_cross = m.bulk[p].cross_sect_surf_area

            eps_s = m.bulk[p].scrap_props.thermal_props.epsilon

            max_b_fac_layers = floor(Int64, no_cells/2)
            b_fac_over_layers = zeros(max_b_fac_layers)

            if A_fs_rad > A_cross/pX_bounded(eps_s, Tmax)
                # TODO: we asume here pX_bounded(eps_s, Tmax) == epsilon_max
                # me may want to enhance the way this limit here is set ...

                @warn   "A_fs_rad bigger than  A_cross/epsilon_max, this would lead "
                        "to unphysical results, consider choosing a smaller " *
                        "cross radiation constant c. Limiting A_fs_rad to A_cross/epsilon_max"
                A_fs_rad =  A_cross/pX_bounded(eps_s, Tmax)
                k = 1

                b_fac_over_layers[1] = 1
            elseif A_fs_rad > 0
                mean_overlap = calc_mean_overlap(A_cross, A_fs_rad, 1, 10_000)

                for j = 1:length(b_fac_over_layers)
                    if j == 1
                        b_fac_over_layers[j] = 1 # with itself
                    elseif j == 2
                        b_fac_over_layers[j] = mean_overlap  # with direct neighbor
                    else
                        b_fac_over_layers[j] = (1-mean_overlap).^(j-2) *
                                mean_overlap # thrugh cells bfore
                    end
                end

                println("b_fac_over_layers: ", b_fac_over_layers[b_fac_over_layers .> 1E-3])

                k = findfirst( x -> x < 1E-3, b_fac_over_layers)
                println("k: ",  k)
            else
                error("Wrong settings of cross_radiation_constant_c?")
            end

            if isnothing(k)
                @warn "b_fac_over_layers[" * string(max_b_fac_layers) *
                        "] is very high, maybe there is " *
                        "something wrong? k will be set to maximum of " *
                        string(max_b_fac_layers)
                k = max_b_fac_layers
            end
        end

        progress = nothing
        if !isapprox(m.sim_phase[p].duration, 0)
            progress = Progress(ceil(Int64, m.sim_phase[p].duration), 1)
        end
        t::Float64 = 0.0

        dt = set_timestepsize(T, m , p, mesh)
        n_time_iters = 0
        println("Initial timestep size: ", dt)
        ran_seconds = 0

        while t <= m.sim_phase[p].duration

            if !isnothing(progress) && ceil(Int64,t) > ran_seconds
                update!(progress, ceil(Int64,t))
                ran_seconds = ceil(Int64,t)
            end

            dt = set_timestepsize(T, m , p, mesh)

            @inbounds for i = 1:1:no_cells
                # heat exchange fluid - solid (gas - scrap)
                heat_flow_exchng_fs = 0.0
                if !scrap_conduction
                    heat_flow_exchng_fs = (A_fs * htc_fs(alpha_fs) *
                            (T.s[i] - T.f[i]))
                else
                    heat_flow_exchng_fs = (A_fs * htc_fs(alpha_fs) *
                            (T.s_sub[i][1] - T.f[i]))
                end

                lambda_dT(TI, T2)::Float64 = harm_mean(pX_bounded(lam_f,T2),
                        pX_bounded(lam_f, TI)) * (T2 - TI)

                # heat conduction through fluid
                heat_flow_cond_f = 0.0
                if i != 1 && i != no_cells
                    heat_flow_cond_f = A_f/dx_f * ( lambda_dT(T.f[i], T.f[i-1]) +
                        lambda_dT(T.f[i], T.f[i+1]) )
                elseif i == no_cells
                    heat_flow_cond_f = A_f/dx_f * lambda_dT(T.f[i], T.f[i-1])
                elseif i == 1
                    heat_flow_cond_f = A_f/dx_f * lambda_dT(T.f[i], T.f[i+1])
                end

                # movement of fluid (gas)
                #integration only for polynomial cp atm, see matprop_functons.jl
                f_dHdt_f_mvmnt(T2)::Float64 = ( mf_f *
                        integrate_pX(cp_f, T.f[i], T2) )

                dHdt_f_mvmnt = 0.0

                if i != 1
                    dHdt_f_mvmnt = f_dHdt_f_mvmnt(T.f[i-1])
                else
                    dHdt_f_mvmnt = f_dHdt_f_mvmnt(m.offgas[p].T_in)
                end

                dH_f::Float64 = dt * (dHdt_f_mvmnt + heat_flow_exchng_fs +
                        heat_flow_cond_f)

                # get the new T from the changed enthalpy of the cell
                Tnew.f[i] = T_of_H( H_T_f, T_H, ( H_of_T(cp_f, T.f[i]) +
                        dH_f/(V_f * pX_bounded(rho_f, T.f[i])) ) )

                heat_flow_exchng_sf[i] = - heat_flow_exchng_fs

                if Tnew.s[i] > Twarn
                    @warn ("Temperatures getting to high possible numerical error")
                end
            end

            @inbounds for i = no_cells:-1:1
                    # radiation models:
                    local heat_flow_rad_front::Float64
                    local heat_flow_rad_in_cross::Float64
                    # front radiation model
                    heat_flow_rad_front = 0.0
                    if front_radiation
                        if i <= k
                            if scrap_conduction
                                T_front = T.s_sub[i][1]
                            else
                                T_front = T.s[i]
                            end

                            if i == 1
                                heat_flow_rad_front= heat_flow_front_rm(
                                        m.sim_phase[p].submodels.furnace_side_rad_model,
                                        pX_bounded(eps_s, T_front), b_fac_over_layers[i],
                                        A_fs_rad, T_front)
                            else
                                heat_flow_rad_front = heat_flow_front_rm(
                                        m.sim_phase[p].submodels.furnace_side_rad_model,
                                        pX_bounded(eps_s, T_front), b_fac_over_layers[i+1],
                                        A_fs_rad, T_front)
                            end
                        end
                    end

                    # cross radiation model
                    heat_flow_rad_in_cross = 0.0
                    if cross_radiation
                        range1 = vcat(collect(i-k:1:i-1), collect(i+1:1:i+k))

                        for r in range1
                            if r > 0 && r <= no_cells
                                if scrap_conduction
                                    T_r = T.s_sub[r][1]
                                    T_i = T.s_sub[i][1]
                                else
                                    T_r = T.s[r]
                                    T_i = T.s[i]
                                end

                                r_dist = abs(i-r)

                                heat_flow_rad_in_cross +=
                                        hf_radiation_par_surf_b_weighted(
                                        pX_bounded(eps_s, T_r), pX_bounded(eps_s, T_i),
                                        b_fac_over_layers[r_dist+1], T_r, T_i )
                            end
                        end
                    end

                if scrap_conduction
                    lambda_dT(TI, T2)::Float64 = ( harm_mean(pX_bounded(lam_s,T2),
                            pX_bounded(lam_s, TI)) * (T2 - TI) )

                    # scrap conduction
                    for j = 1:1:no_sub_cells
                        heat_flow_sub_cond_s = 0.0
                        if j != 1 && j != no_sub_cells
                            heat_flow_sub_cond_s = A_s_s/dx_s *
                                    ( lambda_dT(T.s_sub[i][j], T.s_sub[i][j-1]) +
                                      lambda_dT(T.s_sub[i][j], T.s_sub[i][j+1]) )
                        elseif j == no_sub_cells
                            # symmetry
                                heat_flow_sub_cond_s = A_s_s/dx_s *
                                        lambda_dT(T.s_sub[i][j], T.s_sub[i][j-1])
                        elseif j == 1
                            # convection
                            heat_flow_sub_cond_s = A_s_s/dx_s *
                                    lambda_dT(T.s_sub[i][j],
                                    T.s_sub[i][j+1])
                        end

                        local dHdt_s_mvmnt::Float64
                        dHdt_s_mvmnt = 0.0
                        #scrap movement
                        if i != no_cells
                            dHdt_s_mvmnt = (mf_sub_s * integrate_pX(cp_s,
                                    T.s_sub[i][j], T.s_sub[i+1][j]))
                        else
                            dHdt_s_mvmnt = (mf_sub_s * integrate_pX(cp_s,
                                    T.s_sub[i][j], m.bulk[p].T_charged))
                        end

                        local dH_s::Float64
                        dH_s = 0.0
                        if j == 1
                            # 0.5 because symmetry
                            dH_s = dt * ( 0.5 * (heat_flow_exchng_sf[i] +
                                    heat_flow_rad_front + heat_flow_rad_in_cross) +
                                    heat_flow_sub_cond_s + dHdt_s_mvmnt)
                        else
                            dH_s = dt * (dHdt_s_mvmnt + heat_flow_sub_cond_s)
                        end

                        Tnew.s_sub[i][j] = T_of_H(H_T_s, T_H,
                                ( H_of_T(cp_s, T.s_sub[i][j]) +
                                dH_s/(V_s_sub_c * pX_bounded(rho_s, T.s_sub[i][j])) ) )

                        if Tnew.s_sub[i][j] > Twarn
                            @warn ("Temperatures getting to high possible numerical error")
                        end
                    end
                else
                    #integration only for polynomial cp atm, see matprop_functons.jl
                    f_dHdt_s_mvmnt(T2)::Float64 = (mf_s * integrate_pX(cp_s, T.s[i], T2))

                    local dHdt_s_mvmnt::Float64
                    dHdt_s_mvmnt = 0.0
                    #scrap movement
                    if i != mesh.no_cells
                        dHdt_s_mvmnt = f_dHdt_s_mvmnt(T.s[i+1])
                    else
                        dHdt_s_mvmnt = f_dHdt_s_mvmnt(m.bulk[p].T_charged)
                    end

                    # final new temperature
                    dH_s = dt * ( dHdt_s_mvmnt + heat_flow_exchng_sf[i]
                            + heat_flow_rad_front + heat_flow_rad_in_cross )

                    Tnew.s[i] = T_of_H(H_T_s, T_H, (H_of_T(cp_s, T.s[i]) +
                            dH_s/(V_s * pX_bounded(rho_s, T.s[i]))) )

                    if Tnew.s[i] > Twarn
                        @warn ("Temperatures getting to high possible numerical error")
                    end
                end
            end

            # T.f[:] = copy(Tnew.f)
            unsafe_update_v1_to_v2!(T.f, Tnew.f, no_cells)

            if scrap_conduction
                # T.s_sub[:] = deepcopy(Tnew.s_sub)
                unsafe_update_vv1_to_vv2!(T.s_sub, Tnew.s_sub, no_cells, no_sub_cells )

                # get average scrap cell temperature for post prossing ...
                @inbounds for i = 1:no_cells
                    # usually factor 2 because of symmetry but its cuts out
                    m_sum = 0.0
                    H_T_sub_sum = 0.0
                    @simd for j=1:no_sub_cells
                        m_sum += pX_bounded(rho_s,T.s_sub[i][j])  * V_s_sub_c
                        H_T_sub_sum += H_of_T(cp_s, T.s_sub[i][j]) *
                                pX_bounded(rho_s, T.s_sub[i][j]) * V_s_sub_c
                    end

                    T.s[i] = T_of_H(H_T_s, T_H, H_T_sub_sum/(m_sum))
                end
            else
                # T.s[:] = copy(Tnew.s)
                unsafe_update_v1_to_v2!(T.s, Tnew.s, no_cells)
            end

            t += dt

            if t >= t_save || t > m.sim_phase[p].duration
            push!(results, Results1d1d(phase_start_time + t, deepcopy(T)))
                t_save = t + m.sim_phase[p].dt_save
            end

            n_time_iters += 1
        end

        println("Finsihed p ", p, " after ", n_time_iters, " timesteps!")
    end

    return results
end

"""
set_timestepsize:
set (hoefully) stable according to CFL and Fo Criteria
TODO: Radiation may lead to instabilities, we could to iterative time step Criteria
but  that takes also it time and as EE is not the best after all, we should implement
better method in the future, by now hardcode stable timestep factor ...
"""
@inline function dt_cfl_min(rho::MaterialPropertyOfX, c_vol::Real, m_flow::Real, temp_vec::Vector{<:Real})

    dt_min = 10_000_000.0 # big number
    if !isapprox(m_flow, 0)
        for temp in temp_vec
            local dt
            dt = c_vol * pX_bounded(rho, temp)/m_flow
            if (dt_min > dt)
                dt_min = dt
            end
        end
    end

    return dt_min
end

@inline function dt_fo_min(rho::T1, lam::T2, cp::T3, dx::Real, temp_vec::Vector{<:Real}) where  {T1,T2,T3 <: MaterialPropertyOfX}
    dt_min = 10_000_000.0 # big number
    for temp in temp_vec
        local dt
        dt = (0.5 * pX_bounded(rho, temp) * pX_bounded(cp, temp) * dx^2)/
                pX_bounded(lam, temp)
        if (dt_min > dt)
            dt_min = dt
        end
    end
    return dt_min
end

@inline function set_timestepsize(T::M1d1dVals, m::M1d1d, p::Int64, mesh::M1d1dMesh)::Float64

    lam_f =  m.offgas[p].thermal_props.lambda
    rho_f = m.offgas[p].thermal_props.rho
    cp_f = m.offgas[p].thermal_props.cp_true

    lam_s = m.bulk[p].scrap_props.thermal_props.lambda
    rho_s = m.bulk[p].scrap_props.thermal_props.rho
    cp_s = m.bulk[p].scrap_props.thermal_props.cp_true

    # CFL & FO like timestep size ...
    dtmin_1 = dt_cfl_min(rho_f, mesh.gas_cells.vol, m.offgas[p].mass_flow, T.f)
    dtmin_2 = dt_cfl_min(rho_s, mesh.solid_cells.vol, m.bulk[p].mass_flow, T.s)

    dtmin_3 = dt_fo_min(rho_f, lam_f, cp_f, mesh.gas_cells.len, T.s)

    dtmin_4 = 10_000_000.0 # big number
    if m.sim_phase[p].submodels.scrap_discretization_model_enabled
        T_s_min = minimum(T.s_sub[end])
        T_s_max = maximum(T.s_sub[1])

        if T_s_min < T_s_max
            steps = 5.0
            if (T_s_max - T_s_min)/steps > 1.0
                T_s = collect(T_s_min:steps:T_s_max)
            else
                T_s = [(T_s_max + T_s_min)/2.0]
            end
        else
            # equal?
            T_s = [T_s_max]
        end

        dtmin_4 = dt_fo_min(rho_s, lam_s, cp_s, mesh.solid_1d_submesh.sub_cells.len, T_s)
    end

    dt = minimum((dtmin_1, dtmin_2, dtmin_3, dtmin_4))

    if !isnothing(m.sim_phase[p].submodels.furnace_side_rad_model)
        # hardcoded due to EE instabilities with radiation ...
        dt = 0.05 * dt
        if m.sim_phase[p].submodels.scrap_discretization_model_enabled
            dt = dt/(2.0*mesh.solid_1d_submesh.no_cells)
        end
    else
        dt = 0.48 * dt
    end

    if dt < 1.0E-6
        dt = 1.0E-6
        println("dt1: ", dts[1], "| dt2: ", dts[2], "| dt3: ", dts[3], "| dt4: "
                , dts[4])
        @warn "dt smaller dt_min, maybe something wrong here ..."

        if dt < 0
            error("negative timestep calculated, check bounding of pX functions")
        end
    end

    return dt
end


text:
export DynamicObstacles
using Random

mutable struct DynamicObstacles <: AbstractGridWorld
    world::GridWorldBase{Tuple{Empty,Wall,Obstacle,Goal}}
    agent_pos::CartesianIndex{2}
    agent::Agent
    num_obstacle::Int
    obs_pos_array::Array{CartesianIndex{2},1}
    n::Int
    rng::AbstractRNG
end

function DynamicObstacles(;n=8, agent_start_pos=CartesianIndex(2,2), agent_start_dir=RIGHT, num_obstacle=nothing, rng=Random.GLOBAL_RNG)
    objects = (EMPTY, WALL, OBSTACLE, GOAL)
    w = GridWorldBase(objects, n, n)

    w[EMPTY, 2:n-1, 2:n-1] .= true
    w[WALL, [1,n], 1:n] .= true
    w[WALL, 1:n, [1,n]] .= true
    w[GOAL, n-1, n-1] = true
    w[EMPTY, n-1, n-1] = false

    #w[OBSTACLE,1:n,1:n] .= false
    #num obstacles is n-3 only for n=5 or n=6, it need not be a function of n otherwise
    if num_obstacle === nothing
        num_obstacle = n-3
    end

    obs_pos_array = Array{CartesianIndex{2},1}(undef,num_obstacle)

    obstacles_placed = 0
    while obstacles_placed < num_obstacle
        obs_pos = CartesianIndex(rand(rng, 2:n-1), rand(rng, 2:n-1))
        if (obs_pos == agent_start_pos) || (w[OBSTACLE, obs_pos] == true) || (obs_pos == CartesianIndex(n-1,n-1))
            continue
        else

            w[OBSTACLE, obs_pos] = true
            w[EMPTY, obs_pos] = false
            obstacles_placed = obstacles_placed + 1
            obs_pos_array[obstacles_placed] = obs_pos

        end
    end

    DynamicObstacles(w, agent_start_pos, Agent(dir=agent_start_dir), num_obstacle, obs_pos_array, n, rng)
end

function (w::DynamicObstacles)(::MoveForward)
    dir = get_dir(w.agent)
    dest = dir(w.agent_pos)
    if !w.world[WALL, dest]

        obstacles_replaced = 0

        while obstacles_replaced < w.num_obstacle
            old_pos = w.obs_pos_array[obstacles_replaced+1]
            next_pos = CartesianIndex(old_pos[1]+rand(w.rng, [-1,0,1]),old_pos[2]+rand(w.rng, [-1,0,1]))
            if (!w.world[WALL, next_pos]) && (next_pos != CartesianIndex(w.n-1,w.n-1))
                flag = 0
                #flag indicates whether
                #1)the new-position of the kth obstacle overlaps the new-positions of [1,k-1] obstacles

                for i = 0:obstacles_replaced-1
                    if (w.obs_pos_array[i+1]==next_pos)

                        flag = 1
                        break
                    end
                end

                if flag == 0

                    w.obs_pos_array[obstacles_replaced+1] = next_pos
                    w.world[OBSTACLE, old_pos] = false
                    w.world[EMPTY, old_pos] = true
                    obstacles_replaced += 1
                end
            end
        end
        #println(w.obs_pos_array)
        for obs in w.obs_pos_array
            w.world[OBSTACLE, obs] = true
            w.world[EMPTY, obs] = false
        end
        if w.world[OBSTACLE, dest]
            #end the game
            println("You crashed")
        end

        w.agent_pos = dest
    end
    w
end


text:
# Autogenerated wrapper script for Perl_jll for x86_64-apple-darwin14
export libperl, perl

using Readline_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "DYLD_FALLBACK_LIBRARY_PATH"
LIBPATH_default = "~/lib:/usr/local/lib:/lib:/usr/lib"

# Relative path to `libperl`
const libperl_splitpath = ["lib", "libperl.dylib"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libperl_path = ""

# libperl-specific global declaration
# This will be filled out by __init__()
libperl_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libperl = "@rpath/libperl.dylib"

# Relative path to `perl`
const perl_splitpath = ["bin", "perl"]

# This will be filled out by __init__() for all products, as it must be done at runtime
perl_path = ""

# perl-specific global declaration
function perl(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(perl_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"Perl")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (Readline_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (Readline_jll.LIBPATH_list,))

    global libperl_path = normpath(joinpath(artifact_dir, libperl_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libperl_handle = dlopen(libperl_path)
    push!(LIBPATH_list, dirname(libperl_path))

    global perl_path = normpath(joinpath(artifact_dir, perl_splitpath...))

    push!(PATH_list, dirname(perl_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')
end  # __init__()


text:
export Cos

"""
`Cos([domainType=Float64::Type,] dim_in::Tuple)`

Creates a cosine non-linear operator with input dimensions `dim_in`:
```math
\\cos (\\mathbf{x} ).
```

"""
struct Cos{T,N} <: NonLinearOperator
	dim::NTuple{N,Int}
end

function Cos(DomainType::Type, DomainDim::NTuple{N,Int}) where {N}
	Cos{DomainType,N}(DomainDim)
end

Cos(DomainDim::NTuple{N,Int}) where {N} = Cos{Float64,N}(DomainDim)
Cos(DomainDim::Vararg{Int}) = Cos{Float64,length(DomainDim)}(DomainDim)

function mul!(y::AbstractArray{T,N}, L::Cos{T,N}, x::AbstractArray{T,N}) where {T,N}
	y .= cos.(x)
end

function mul!(y::AbstractArray,
              J::AdjointOperator{Jacobian{A,TT}},
              b::AbstractArray) where {T,N, A<: Cos{T,N}, TT <: AbstractArray{T,N}}
    L = J.A
    y .= -conj.(sin.(L.x)).*b
end

fun_name(L::Cos) = "cos"

size(L::Cos) = (L.dim, L.dim)

domainType(L::Cos{T,N}) where {T,N} = T
codomainType(L::Cos{T,N}) where {T,N} = T


text:
# This file exemplifies the workflow from data input to optimization result generation
using ClustForOpt
using Clp
## LOAD DATA ##
state="GER_1" # or "GER_18" or "CA_1" or "TX_1"
years=[2015, 2016] #2016 works for GER_1 and CA_1, GER_1 can also be used with 2006 to 2016 and, GER_18 is 2015 TX_1 is 2008
# laod ts-data
ts_input_data = load_timeseries_data("CEP", state; T=24, years=years) #CEP
# load cep-data
cep_data = load_cep_data_provided(state)

## CLUSTERING ##
# run aggregation with kmeans
ts_clust_data = run_clust(ts_input_data;method="kmeans",representation="centroid",n_init=100,n_clust=5) # default k-means make sure that n_init is high enough otherwise the results could be crap and drive you crazy

# run aggregation with kmeans and have periods segmented
ts_seg_data = run_clust(ts_input_data;method="kmeans",representation="centroid",n_init=1000,n_clust=5, n_seg=4) # default k-means make sure that n_init is high enough otherwise the results could be crap and drive you crazy

# run no aggregation just get ts_full_data
ts_full_data = run_clust(ts_input_data;method="kmeans",representation="centroid",n_init=1,n_clust=365) # default k-means


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

export BrickletTilt
"""
Detects inclination of Bricklet (tilt switch open/closed)
"""
mutable struct BrickletTilt <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export get_tilt_state
"""
    $(SIGNATURES)

Returns the current tilt state. The state can either be

* 0 = Closed: The ball in the tilt switch closes the circuit.
* 1 = Open: The ball in the tilt switch does not close the circuit.
* 2 = Closed Vibrating: The tilt switch is in motion (rapid change between open and close).

.. image:: /Images/Bricklets/bricklet_tilt_mechanics.jpg
   :scale: 100 %
   :alt: Tilt states
   :align: center
   :target: ../../_images/Bricklets/bricklet_tilt_mechanics.jpg
"""
function get_tilt_state(device::BrickletTilt)
    return device.deviceInternal.get_tilt_state()
end

export enable_tilt_state_callback
"""
    $(SIGNATURES)

Enables the :cb:`Tilt State` callback.
"""
function enable_tilt_state_callback(device::BrickletTilt)
    device.deviceInternal.enable_tilt_state_callback()
end

export disable_tilt_state_callback
"""
    $(SIGNATURES)

Disables the :cb:`Tilt State` callback.
"""
function disable_tilt_state_callback(device::BrickletTilt)
    device.deviceInternal.disable_tilt_state_callback()
end

export is_tilt_state_callback_enabled
"""
    $(SIGNATURES)

Returns *true* if the :cb:`Tilt State` callback is enabled.
"""
function is_tilt_state_callback_enabled(device::BrickletTilt)
    return device.deviceInternal.is_tilt_state_callback_enabled()
end

export get_identity
"""
    $(SIGNATURES)

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

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

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

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


text:
module MCluster
using DataFrames

    mutable struct MC
        label
        n::Int64
        ls::Array{Float64}
        ss::Array{Float64}
        centroid::Array{Float64}

        function MC()
            x = new()
        end

        function MC(sample::Array{T, 1}, label) where {T<:Number}
            v_size = length(sample)
            x = new()
            x.n = 1
            x.label = label
            x.ls = zeros(T, v_size)
            x.ss = zeros(T, v_size)
            for i = 1:v_size
                x.ls[i] = x.ls[i] + sample[i]
                x.ss[i] = x.ss[i] + sample[i] ^ 2
            end
            x.centroid = copy(sample)
            return x
        end
    end

    function mc_set_centroid(mc::MC)
        mc.centroid = mc.ls / mc.n
    end

    function predict_r(mc::MC, sample::Array{T}) where {T<:Number}
        ss = mc.ss + sample .^2
        ls = mc.ls + sample
        n = mc.n + 1
        return abs(sum((ss / n) - (ls / n).^ 2)) ^ (1/2)
    end

    function append!(mc::MC, sample::Array{T}) where {T<:Number}
        mc.n = mc.n + 1
        mc.ls = mc.ls + sample
        mc.ss = mc.ss + sample .^ 2
        mc_set_centroid(mc)
    end

    function merge(mc_a::MC, mc_b::MC)
        new_mc = MC()
        new_mc.n = mc_a.n + mc_b.n
        new_mc.ls = mc_a.ls + mc_b.ls
        new_mc.ss = mc_a.ss + mc_b.ss
        mc_set_centroid(new_mc)
        return new_mc

    end
end #module


text:
export getreducedmargin, updatereducedmargin

"""
$(TYPEDSIGNATURES)

Computes the reduced margin of a set of multi-indices.
See Baptista, R., Zahm, O., & Marzouk, Y. (2020).
An adaptive transport framework for joint and conditional density estimation. arXiv preprint arXiv:2009.10303.
"""
function getreducedmargin(midx::Array{Int64,2})

    if isempty(midx)
        return zeros(Int64,0,0)
    end

    n, d = size(midx)

    # add [0 0 0] to compute the reduced margin
    if !any(zeros(Int64,d) in eachslice(midx; dims = 1))
            idx = vcat(zeros(Int64,1,d), midx)
    else
            idx = midx
    end

    n, d = size(idx)

    neighbours = zeros(Int64,(n,d,d))
    eye = Matrix(1*I,d,d)

    @inbounds for i=1:d
        neighbours[:,i,:]  .= idx
    end

    @inbounds for i=1:n
        neighbours[i,:,:]  .+= eye
    end

    neighbours = reshape(neighbours,(n*d,d))

    margintmp = setdiff(eachslice(neighbours; dims = 1), eachslice(idx;dims = 1))
    cardinal = size(margintmp,1)
    margin = zeros(Int64,cardinal,d)
    @inbounds for i=1:cardinal
        margin[i,:] = margintmp[i]
    end

    # margin = margin[sortperm(margin[:, 1]), :] # sorted by the 1st column
    margin = sortslices(margin; dims = 1)

    neighbours = zeros(Int64, (cardinal,d,d))

    @inbounds for i=1:d
        neighbours[:,i,:]  .= margin
    end

    @inbounds for i=1:cardinal
        neighbours[i,:,:]  .-= eye
    end

    neighbours = reshape(neighbours, cardinal*d,d)

    ok = Bool[any(x in eachslice(idx; dims = 1)) for x in eachslice(neighbours; dims = 1)]

    isout = Bool[any(x.<0) for x in eachslice(neighbours; dims = 1)]

    ok .|=  isout

    ok = reshape(ok, (cardinal, d))

    keep = Bool[all(x) for x in eachslice(ok; dims = 1)]

    return margin[keep,:]
end

"""
$(TYPEDSIGNATURES)

Updates the reduced margin.

See Baptista, R., Zahm, O., & Marzouk, Y. (2020).
An adaptive transport framework for joint and conditional density estimation. arXiv preprint arXiv:2009.10303.
"""
function updatereducedmargin(lowerset::Array{Int64,2}, reduced_margin::Array{Int64,2}, idx::Int64)
    d = size(lowerset, 2)
    nr = size(reduced_margin, 1)
    @assert idx <= nr "idx is larger than dimension of the reduced_margin"
    newidx = reduced_margin[idx:idx,:]
#     @show lowerset
#     @show reduced_margin
#     @show newidx
    lowerset = vcat(lowerset, newidx)
#     @show lowerset
    # Remove the idx-th line of reducedmargin
    toremove = Bool[!(i in idx) for i=1:nr]
    # reduced_margin = reduced_margin[1:end .!= idx,:]
    reduced_margin = reduced_margin[toremove,:]

#     @show toremove
#     @show reduced_margin

    eye =  Matrix(1*I,d,d)
    # Update the reduced margin
    candidate  = repeat(newidx, outer = (d,1)) + eye

#     @show candidate
    ok = falses(d)
#     @show ok

    @inbounds for i = 1:d
#         @show i
#         @show repeat(candidate[i,:]', outer = (d,1))
#         @show size(repeat(candidate[i,:]', outer = (d,1)))
        parents_of_candidate = repeat(candidate[i,:]', outer = (d,1)) - eye
#         @show parents_of_candidate
        tokeep = Bool[!any(x .< 0) for x in eachslice(parents_of_candidate; dims = 1)]
#         @show tokeep
        parents_of_candidate = parents_of_candidate[tokeep,:]
#         @show parents_of_candidate

#         @show [any(x in eachslice(lowerset; dims = 1)) for x in eachslice(parents_of_candidate; dims = 1)]
        ok[i] = all([any(x in eachslice(lowerset; dims = 1)) for x in eachslice(parents_of_candidate; dims = 1)])
#        @show ok
    end
    #
    candidate = candidate[ok,:]
    #
    # # Add the candidates to the reduced margin
    reduced_margin = vcat(reduced_margin, candidate)

    # reduced_margin  = sortslices(reduced_margin; dims = 1)
    # lowerset = sortslices(lowerset; dims = 1)
    return lowerset, reduced_margin
end


text:
module SlurmAddAllocatedProcs

using ClusterManagers
export addprocs_slurm_allocated
export parse_slurm_tasks_per_node

function parse_slurm_tasks_per_node(SLURM_TASKS_PER_NODE)
    v = split(SLURM_TASKS_PER_NODE, ",")
    pattern = r"([0-9]+)\(x([0-9]+)\)"
    np = 0
    n = 0
    for el in v
        el_Int = tryparse(Int, el)
        if el_Int isa Int
            np += el_Int
            n += 1
        else
            rm = match(pattern, el)
            if rm !== nothing
                ni = parse(Int, rm[2])
                n += ni
                p_ni = parse(Int, rm[1])
                np += p_ni*ni
            end
        end
    end
    if n == 0
        error("Could not infer the number of tasks")
    end
    return np
end

"""
    addprocs_slurm_allocated()

Add workers when run in a batch Slurm job. The number of tasks to use is inferred from the environment
variables set by the resource allocation step (eg. using `sbatch` flags).
"""
function addprocs_slurm_allocated()
    SLURM_TASKS_PER_NODE = get(ENV, "SLURM_TASKS_PER_NODE", "")
    np = parse_slurm_tasks_per_node(SLURM_TASKS_PER_NODE)
    addprocs_slurm(np)
end

end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

struct DotDict
    inner::Dict
end

DotDict() = DotDict(Dict())

function Base.setproperty!(d::DotDict, key::Symbol, value)
    setindex!(getfield(d, :inner), value, key)
end

function Base.getproperty(d::DotDict, key::Symbol)
    (key == :inner ? getfield(d, :inner) : d.inner[key])
end

function Base.getindex(d::DotDict, key::Int64)
    d.inner[Symbol(key)]
end

function Base.getindex(d::DotDict, key::Symbol)
    d.inner[key]
end

function Base.keys(d::DotDict)
    keys(d.inner)
end

function Base.values(d::DotDict)
    values(d.inner)
end

function Base.iterate(d::DotDict)
    iterate(values(d.inner))
end

function Base.iterate(d::DotDict, v::Int64)
    iterate(values(d.inner), v)
end

function Base.length(d::DotDict)
    length(values(d.inner))
end

function Base.show(io::IO, d::DotDict)
    print(io, "DotDict with $(length(keys(d.inner))) entries:\n")
    count = 0
    for k in keys(d.inner)
        count += 1
        if count > 10
            print(io, "  ...\n")
            break
        end
        print(io, "  :$(k) => $(d.inner[k])\n")
    end
end

function recursive_to_dot_dict(el)
    if typeof(el) == Dict{String, Any}
        return DotDict(Dict(Symbol(k) => recursive_to_dot_dict(el[k]) for k in keys(el)))
    else
        return el
    end
end

export recursive_to_dot_dict


text:
APERTURES = [
    CircularAperture,
    CircularAnnulus
]

PARAMS = [
    (3),
    (3, 5)
]

###########################
# Some helpers for testing
area(c::CircularAperture) = π * c.r^2
area(c::CircularAnnulus) = π * (c.r_out^2 - c.r_in^2)

@testset "outside - $AP" for (AP, params) in zip(APERTURES, PARAMS)
    data = ones(10, 10)
    aperture = AP(-60, 60, params...)
    @test aperture_photometry(aperture, data).aperture_sum ≈ 0
end

@testset "inside zeros - $AP" for (AP, params) in zip(APERTURES, PARAMS)
    data = zeros(40, 40)
    aperture = AP(20.0, 20.0, params...)

    table_cent = aperture_photometry(aperture, data, method = :center)
    table_sub = aperture_photometry(aperture, data, method = (:subpixel, 10))
    table_ex = aperture_photometry(aperture, data, method = :exact)

    @test table_ex.aperture_sum ≈ 0
    @test table_sub.aperture_sum ≈ 0
    @test table_cent.aperture_sum ≈ 0

end

@testset "inside ones - $AP" for (AP, params) in zip(APERTURES, PARAMS)
    data = ones(40, 40)
    aperture = AP(20.0, 20.0, params...)

    table_cent = aperture_photometry(aperture, data, method = :center)
    table_sub = aperture_photometry(aperture, data, method = (:subpixel, 10))
    table_ex = aperture_photometry(aperture, data, method = :exact)

    true_flux = area(aperture)

    @test table_ex.aperture_sum ≈ true_flux
    @test table_sub.aperture_sum ≈ table_ex.aperture_sum atol = 0.1
    @test table_cent.aperture_sum < table_ex.aperture_sum

end

function test_aperture(data, aperture)
    error = ones(size(data))

    table_cent = aperture_photometry(aperture, data, error, method = :center)
    table_sub = aperture_photometry(aperture, data, error, method = (:subpixel, 12))
    table_ex = aperture_photometry(aperture, data, error, method = :exact)

    true_flux = area(aperture)
    true_err = sqrt(true_flux)

    @test table_ex.aperture_sum ≈ true_flux
    @test table_sub.aperture_sum ≈ table_ex.aperture_sum atol = 0.1
    @test table_cent.aperture_sum < table_ex.aperture_sum

    @test table_ex.aperture_sum_err ≈ true_err
    @test table_sub.aperture_sum_err ≈ table_ex.aperture_sum_err atol = 0.1
    @test table_cent.aperture_sum_err < table_ex.aperture_sum_err
end

@testset "errors - CircularAperture" begin
    data = ones(40, 40)
    aperture = CircularAperture(20, 20, 10)
    test_aperture(data, aperture)

    # apertures = [aperture, CircularAperture(25, 25, 10)]
    # test_aperture(data, apertures)
end

@testset "errors - CircularAnnulus" begin
    data = ones(40, 40)
    aperture = CircularAnnulus(20, 20, 8, 10)
    test_aperture(data, aperture)

    # apertures = [aperture, CircularAnnulus(25, 25, 8, 10)]
    # test_aperture(data, apertures)
end

@testset "partial overlap" begin
    data = ones(20, 20)
    error = ones(size(data))
    positions = [10.5 10.5; 1 1; 1 20; 20 1; 20 20]
    apertures = [CircularAperture(positions[i, :], 5) for i in axes(positions, 1)]

    table = aperture_photometry(apertures, data, error)
    @test table.aperture_sum[1] ≈ 25π
    @test all(table.aperture_sum[2:end] .< 25π)
end