text:
module EqualizerFilters

import DSP.Biquad
import DSP.SecondOrderSections
using DSP

include("equalizerAPOstrings.jl")
include("DSPjlBiquads.jl")

export LP, HP, BP, NO, AP, PK, LS, HS, Biquad, eqAPOstring, SecondOrderSections

#default sampling rate
fs = 48e3

"""
    EqualizerFilters.sampling_rate(custom_sampling_rate)

Set the sampling rate for calculation of the biquads in Hz.

Defaults to 48 kHz.
"""
function sampling_rate(custom_sampling_rate)
    global fs = custom_sampling_rate
end

"""
    EqualizerFilters.sampling_rate()

Return the current sampling rate for biquad calculation in Hz.

Defaults to 48 kHz.
"""
function sampling_rate()
    fs
end

"""
    biquadPrecalculations(dbGain, freq, srate, Q)

Returns `(A, sn, cs, alpha, beta)` that have to be calculated for all
filter types.

The parameter A is not needed for gain-less filters.
The calling functions simply pass 0 dB Gain and do not use the returned parameter
`A`.
"""
function biquadPrecalculations(dbGain, freq, srate, Q)
    A = 10^(dbGain/40)
    omega = 2pi * freq / srate
    sn = sin(omega)
    cs = cos(omega)
    alpha = sn / (2 * Q)
    beta = 2 * sqrt(A) * alpha
    return A, sn, cs, alpha, beta
end

"""
    LP(f, Q, fs=sampling_rate)

Return lowpass biquad at frequency `f` in Hz with quality factor `Q`.

This second order lowpass falls of with 12 dB/oct.

If not given in the arguments, the sampling rate defaults to 48 kHz or can be
set with [`EqualizerFilters.sampling_rate`](@ref).
"""
function LP(f, Q, fs=fs)
    A, sn, cs, alpha, beta = biquadPrecalculations(0, f, fs, Q)
    b0 = (1 - cs) / 2
    b1 = 1 - cs
    b2 = (1 - cs) / 2
    a0 = 1 + alpha
    a1 = -2 * cs
    a2 = 1 - alpha
    biquadCoefs = [b0, b1, b2, a1, a2]./a0
    return Biquad(biquadCoefs...)
end

"""
    HP(f, Q, fs=sampling_rate)

Return highpass biquad at frequency `f` in Hz with quality factor `Q`.

This second order highpass falls of with 12 dB/oct.

If not given in the arguments, the sampling rate defaults to 48 kHz or can be
set with [`EqualizerFilters.sampling_rate`](@ref).
"""
function HP(f, Q, fs=fs)
    A, sn, cs, alpha, beta = biquadPrecalculations(0, f, fs, Q)
    b0 = (1 + cs) / 2;
    b1 = -(1 + cs);
    b2 = (1 + cs) / 2;
    a0 = 1 + alpha;
    a1 = -2 * cs;
    a2 = 1 - alpha;
    biquadCoefs = [b0, b1, b2, a1, a2]./a0
    return Biquad(biquadCoefs...)
end

"""
    BP(f, Q, fs=sampling_rate)

Return bandpass biquad at frequency `f` in Hz with quality factor `Q`.

If not given in the arguments, the sampling rate defaults to 48 kHz or can be
set with [`EqualizerFilters.sampling_rate`](@ref).
"""
function BP(f, Q, fs=fs)
    A, sn, cs, alpha, beta = biquadPrecalculations(0, f, fs, Q)
    b0 = alpha;
    b1 = 0;
    b2 = -alpha;
    a0 = 1 + alpha;
    a1 = -2 * cs;
    a2 = 1 - alpha;
    biquadCoefs = [b0, b1, b2, a1, a2]./a0
    return Biquad(biquadCoefs...)
end

"""
    NO(f, Q, fs=sampling_rate)

Return notch biquad at frequency `f` in Hz with quality factor `Q`.

If not given in the arguments, the sampling rate defaults to 48 kHz or can be
set with [`EqualizerFilters.sampling_rate`](@ref).
"""
function NO(f, Q, fs=fs)
    A, sn, cs, alpha, beta = biquadPrecalculations(0, f, fs, Q)
    b0 = 1;
    b1 = -2 * cs;
    b2 = 1;
    a0 = 1 + alpha;
    a1 = -2 * cs;
    a2 = 1 - alpha;
    biquadCoefs = [b0, b1, b2, a1, a2]./a0
    return Biquad(biquadCoefs...)
end

"""
    AP(f, Q, fs=sampling_rate)

Return allpass biquad at frequency `f` in Hz with quality factor `Q`.

This second order allpass rotates the phase by a total of ``2π``.

If not given in the arguments, the sampling rate defaults to 48 kHz or can be
set with [`EqualizerFilters.sampling_rate`](@ref).
"""
function AP(f, Q, fs=fs)
    A, sn, cs, alpha, beta = biquadPrecalculations(0, f, fs, Q)
    b0 = 1 - alpha;
    b1 = -2 * cs;
    b2 = 1 + alpha;
    a0 = 1 + alpha;
    a1 = -2 * cs;
    a2 = 1 - alpha;
    biquadCoefs = [b0, b1, b2, a1, a2]./a0
    return Biquad(biquadCoefs...)
end

"""
    PK(f, dbGain, Q, fs=sampling_rate)

Return peaking biquad with Gain `dbGain` in dB at frequency `f` in Hz
with quality factor `Q`.

In this implementation the Gain at 0 Hz and at the Nyquist frequency will
always be 0 dB. For an Orfanidis-like peak EQ with analog-inspired
Nyquist-gain see the package LakeBiquads.jl.

If not given in the arguments, the sampling rate defaults to 48 kHz or can be
set with [`EqualizerFilters.sampling_rate`](@ref).
"""
function PK(f, dbGain, Q, fs=fs)
    A, sn, cs, alpha, beta = biquadPrecalculations(dbGain, f, fs, Q)
    b0 = 1 + (alpha * A);
    b1 = -2 * cs;
    b2 = 1 - (alpha * A);
    a0 = 1 + (alpha / A);
    a1 = -2 * cs;
    a2 = 1 - (alpha / A);
    biquadCoefs = [b0, b1, b2, a1, a2]./a0
    return Biquad(biquadCoefs...)
end

"""
    LS(f, dbGain, Q, fs=sampling_rate)

Return low shelf biquad with Gain `dbGain` in dB, center frequency `f` in Hz and
quality factor `Q`.

If not given in the arguments, the sampling rate defaults to 48 kHz or can be
set with [`EqualizerFilters.sampling_rate`](@ref).
"""
function LS(f, dbGain, Q, fs=fs)
    A, sn, cs, alpha, beta = biquadPrecalculations(dbGain, f, fs, Q)
    b0 = A * ((A + 1) - (A - 1) * cs + beta);
    b1 = 2 * A * ((A - 1) - (A + 1) * cs);
    b2 = A * ((A + 1) - (A - 1) * cs - beta);
    a0 = (A + 1) + (A - 1) * cs + beta;
    a1 = -2 * ((A - 1) + (A + 1) * cs);
    a2 = (A + 1) + (A - 1) * cs - beta;
    biquadCoefs = [b0, b1, b2, a1, a2]./a0
    return Biquad(biquadCoefs...)
end

"""
    HS(f, dbGain, Q, fs=sampling_rate)

Return high shelf biquad with Gain `dbGain` in dB, center frequency `f` in Hz and
quality factor `Q`.

If not given in the arguments, the sampling rate defaults to 48 kHz or can be
set with [`EqualizerFilters.sampling_rate`](@ref).
"""
function HS(f, dbGain, Q, fs=fs)
    A, sn, cs, alpha, beta = biquadPrecalculations(dbGain, f, fs, Q)
    b0 = A * ((A + 1) + (A - 1) * cs + beta);
    b1 = -2 * A * ((A - 1) + (A + 1) * cs);
    b2 = A * ((A + 1) + (A - 1) * cs - beta);
    a0 = (A + 1) - (A - 1) * cs + beta;
    a1 = 2 * ((A - 1) - (A + 1) * cs);
    a2 = (A + 1) - (A - 1) * cs - beta;
    biquadCoefs = [b0, b1, b2, a1, a2]./a0
    return Biquad(biquadCoefs...)
end

end # module


text:
using ALASMPCCmod

ALASMPCCmod_success=true
println("On teste le module ALASMPCCmod")

#Bilan : donner une sortie bilan
if ALASMPCCmod_success==true
 println("ALASMPCCmod.jl passes the test !")
else
 println("ALASMPCCmod.jl contains some error")
end

ext: jl


text:
""" prepare_line_list( linelist_fn, all_spectra, pipeline; recalc, convert_air_to_vacuum, orders_to_use, Δv_to_avoid_tellurics, v_center_to_avoid_tellurics )
linelist_fn is the full path to the line list file
linelist_fn is loaded as a VALD mask if it contains substring "VALD", and linelist_fn is loaded as an ESPRESSO mask if it contains the substring "espresso.mas".
convert_air_to_vacuum is a boolean (defualt value true) that determines whether to convert to vaccum wavelengths when loading VALD and ESPRESSO masks.
"""
function prepare_line_list( linelist_fn::String, all_spectra::AbstractVector{SpecT}, pipeline::PipelinePlan; recalc::Bool = false, convert_air_to_vacuum::Bool = true,
         orders_to_use = RvSpectML.orders_to_use_default(first(all_spectra).inst),
         Δv_to_avoid_tellurics::Real = EchelleInstruments.default_Δv_to_avoid_tellurics, v_center_to_avoid_tellurics::Real = 0.0, verbose::Bool = false ) where { SpecT <: AbstractSpectra }
   @assert length(linelist_fn) >= 1
   @assert length(all_spectra) >= 1
   local line_list_df
   if need_to(pipeline,:read_line_list) || recalc
      lambda_range_with_good_data = get_λ_range(all_spectra)
      if verbose println("# Reading line list for CCF: ", linelist_fn, ".")  end
      if occursin(r"espresso\.mas$", linelist_fn)
         line_list_df = EchelleCCFs.read_linelist_espresso(linelist_fn, convert_air_to_vacuum=convert_air_to_vacuum)
      elseif occursin(r"VALD", linelist_fn)
         line_list_df = EchelleCCFs.read_linelist_vald(linelist_fn, convert_air_to_vacuum=convert_air_to_vacuum)
      else
         line_list_df = EchelleCCFs.read_linelist_rvspectml(linelist_fn)
      end
      inst = first(all_spectra).inst
      inst_module = get_inst_module(first(all_spectra).inst)
      if verbose   println("# extrema(λ) = ", extrema(line_list_df.lambda), "   pre-filter" )   end
      λmin_data = maximum(map(obsid->all_spectra[obsid].λ[min_col_default(inst,first(orders_to_use)),first(orders_to_use)],1:length(all_spectra)))
      λmax_data = minimum(map(obsid->all_spectra[obsid].λ[max_col_default(inst,last(orders_to_use)),last(orders_to_use)],1:length(all_spectra)))
      if verbose   println("# λmin_data = ",λmin_data, "  λmax_data = ",λmax_data)   end
      line_list_df = inst_module.filter_line_list(line_list_df,first(all_spectra).inst, λmin=λmin_data, λmax=λmax_data)
      if verbose   println("# extrema(λ) = ", extrema(line_list_df.lambda), "   post-filter" )   end

      #=
      println("# Adding col-based λ boundaries...")
      line_list_df = add_line_boundaries_to_line_list(line_list_df,  Δv_to_avoid_tellurics=Δv_to_avoid_tellurics, v_center_to_avoid_tellurics=v_center_to_avoid_tellurics )

      order_info = get_order_info(all_spectra, orders_to_use=orders_to_use)
      #println("# order_info contains: ", names(order_info))
      #println("# line_list_df contains: ", names(line_list_df))
      line_list_df = RvSpectML.assign_lines_to_orders(line_list_df, order_info)
      =#

      if eltype(all_spectra) <: EchelleInstruments.AnyEXPRES
         RvSpectMLBase.discard_pixel_mask(all_spectra)
         RvSpectMLBase.discard_excalibur_mask(all_spectra)
      end
      set_cache!(pipeline,:read_line_list,deepcopy(line_list_df))
      dont_need_to!(pipeline,:read_line_list);
    end

    if need_to(pipeline,:clean_line_list_tellurics) || recalc
      if verbose println("# Removing lines with telluric contamination.")  end
      @assert !need_to(pipeline,:read_line_list)
      @assert !need_to(pipeline,:read_spectra)
      line_list_df = make_clean_line_list_from_tellurics(line_list_df, all_spectra,
               Δv_to_avoid_tellurics = Δv_to_avoid_tellurics, v_center_to_avoid_tellurics=v_center_to_avoid_tellurics)
      #if verbose   println("# extrema(λ) = ", extrema(line_list_df.lambda), "   post-clean-tellurics" )   end
      if typeof(first(all_spectra).inst) <: EchelleInstruments.AnyEXPRES
         # RvSpectML.discard_tellurics(all_spectra)  # Keep, since individual line fits use the tellurics info data later
      elseif typeof(first(all_spectra).inst) <: EchelleInstruments.AnyNEID
         # Nothing to do for NEID
      else
         @warn("Removing lines with telluric contamination currently only works with EXPRES data.")
      end

      set_cache!(pipeline,:clean_line_list_tellurics,deepcopy(line_list_df))
      dont_need_to!(pipeline,:clean_line_list_tellurics);
    end

    if need_to(pipeline,:calc_line_list_snr_weights) || recalc
      #order_info = get_order_info(all_spectra, orders_to_use=orders_to_use)
      #println("# order_info contains: ", names(order_info))
      #println("# line_list_df contains: ", names(line_list_df))
      #line_list_df = RvSpectML.assign_lines_to_orders(line_list_df, order_info)
      if verbose   println("# line_list extrema(λ) = ", extrema(line_list_df.lambda)) end #, "   post-assign-to-orders" )   end
      RvSpectML.calc_snr_weights_for_lines!(line_list_df, all_spectra)
      set_cache!(pipeline,:calc_line_list_snr_weights,line_list_df)
      dont_need_to!(pipeline,:calc_line_list_snr_weights);
    end

  if     has_cache(pipeline,:calc_line_list_snr_weights) return read_cache(pipeline,:calc_line_list_snr_weights)
  elseif has_cache(pipeline,:clean_line_list_tellurics)  return read_cache(pipeline,:clean_line_list_tellurics)
  elseif has_cache(pipeline,:read_line_list)             return read_cache(pipeline,:read_line_list)
  else   @error("Invalid pipeline state.")               end
end


text:
# this is a convenience class to connect Transformations in a chain.
# One could represent an ANN like:
#   affine(n1,n2) --> relu(n2) --> affine(n2,n3) --> logistic(n3)
# The output of a chain could be fed into a loss model for backprop
type Chain{T,P<:Params,PREP<:PreprocessStep} <: Learnable
    nin::Int
    nout::Int
    input::SumNode{T,1}
    output::OutputNode{T,1}
    ts::Vector{Transformation}
    params::P
    prep::PREP
    grad_calc::Symbol
    Bs::Vector{Nullable{Matrix{T}}}  # these are the fixed/random matrices Bᵢ for the DFA method
end

Chain(ts::Transformation...; kw...) = Chain(Float64, ts...; kw...)
Chain{T}(::Type{T}, ts::Transformation...; kw...) = Chain(T, convert(Array{Transformation}, collect(ts)); kw...)

function Chain{T,TR<:Transformation}(::Type{T}, ts::AbstractVector{TR};
                                     prep::PreprocessStep = NoPreprocessing(),
                                     grad_calc::Symbol = :backprop)
    link_nodes!(ts)

    # initialize DFA method
    Bs = Nullable{Matrix{T}}[]
    if grad_calc == :dfa
        for (i,t) in enumerate(ts)
            nΘ = params_length(t)
            Bi = if nΘ == 0 || i >= length(ts)-1
                # don't do anything for non-learnable transformations,
                # or the final transformation (we'll just backprop that one)
                Nullable{Matrix{T}}()
            else
                ny = output_length(t)
                # B = zeros(T, input_length(t), output_length(ts[end]))
                B = zeros(T, output_length(t), output_length(ts[end]))
                initialize_weights!(B)
                fill!(params(t), zero(T))  # reset params to 0... should we keep this?
                Nullable{Matrix{T}}(B)
            end
            push!(Bs, Bi)
        end
    end

    Chain(
        input_length(ts[1]),
        output_length(ts[end]),
        input_node(ts[1]),
        output_node(ts[end]),
        ts,
        consolidate_params(T, ts),
        prep,
        grad_calc,
        Bs
    )
end

Base.getindex(chain::Chain, i) = chain.ts[i]
Base.length(chain::Chain) = length(chain.ts)
Base.size(chain::Chain) = size(chain.ts)
Base.endof(chain::Chain) = length(chain)

function Base.show{T}(io::IO, chain::Chain{T})
    println(io, "Chain{$T}(")
    for t in chain.ts
        print(io, "   ")
        show(io, t)
        println(io)
    end
    print(io, ") ")
end

# # add a transformation to the end of the pipeline
# function Base.push!(chain::Chain, t::Transformation)
#     # the new input should be linked to:
#     #   - the chain input if it's first, otherwise
#     #   - the last transformation's output
#     link_nodes!(t.input, isempty(chain.ts) ? chain.input : chain.ts[end].output)
#
#     # we always update the chain output to match the output of the last transformation
#     link_nodes!(t.output, chain.output)
#
#     push!(chain.ts, t)
#     chain
# end

# now that it's set up, just go forward to transform, or backwards to backprop

function transform!(chain::Chain)
    if !isa(chain.prep, NoPreprocessing)
        # update the whitener and transform the input data
        copy!(input_value(chain.prep), input_value(chain))
        learn!(chain.prep)
        copy!(input_value(chain), transform!(chain.prep))
    end
    for (i,t) in enumerate(chain.ts)
        i > 1 && transform!(input_node(t))
        transform!(t)
    end
    output_value(chain)
end

function grad!(chain::Chain)
    if chain.grad_calc == :backprop
        for (i,t) in enumerate(reverse(chain.ts))
            i > 1 && grad!(output_node(t))
            grad!(t)
        end
        input_grad(chain)
    elseif chain.grad_calc == :dfa
        #= going through the list:
            - skip non-learnables
            - final transform uses normal backprop
            - other learnables use DFA: ∇yᵢ = Bᵢ∇yₘ
        =#
        ∇out = output_grad(chain)
        n = length(chain.ts)
        # @show n ∇out
        for (j,t) in enumerate(reverse(chain.ts))
            i = n - j + 1
            # @show i, t
            if isnull(chain.Bs[i])
                # non-learnables get backprop
                i < n && grad!(output_node(t))
                # @show "null" grad!(t)
                grad!(t)
            else
                if i == n
                    # final node gets backprop
                    grad!(t)
                else
                    # other learnables get DFA
                    Bᵢ = get(chain.Bs[i])
                    # @show Bᵢ
                    grad!(t, Bᵢ * ∇out)
                    # ∇y = if isa(t, Affine) || isa(t, Linear)
                    #     w = t.params.views[1]
                    #     w * (Bᵢ * ∇out)
                    # else
                    #     error("Only Affine and Linear Learnables are supported with DFA")
                    # end
                    # grad!(t, ∇y)
                end
            end
            # @show output_grad(t) grad(t) input_grad(t)
        end
        input_grad(chain)
    else
        error("Unsupported grad_calc: ", chain.grad_calc)
    end
end

# transform!(chain::Chain) = (foreach(transform!, chain.ts); chain.output.val)
# grad!(chain::Chain) = foreach(grad!, reverse(chain.ts))

function reset_params!{T}(chain::Chain{T}, θ::AbstractVector, ∇::AbstractVector)
    chain.params = consolidate_params(T, chain.ts, θ=θ, ∇=∇)
end

link_nodes!(chain::Chain, i_from::Int, i_to::Int) = link_nodes!(chain.ts[i_from], chain.ts[i_to])


text:
mutable struct Scene <: AbstractScene
    events::Events

    px_area::Node{IRect2D}
    camera::Camera
    camera_controls::RefValue
    limits::Node{FRect3D}

    transformation::Transformation

    plots::Vector{AbstractPlot}
    theme::Attributes
    children::Vector{Scene}
    current_screens::Vector{AbstractScreen}

    function Scene(
            events::Events,
            px_area::Node{IRect2D},
            camera::Camera,
            camera_controls::RefValue,
            limits::Node,
            transformation::Transformation,
            plots::Vector{AbstractPlot},
            theme::Attributes,
            children::Vector{Scene},
            current_screens::Vector{AbstractScreen},
        )
        obj = new(events, px_area, camera, camera_controls, limits, transformation, plots, theme, children, current_screens)
        finalizer(obj) do obj
            # save_print("Freeing scene")
            close_all_nodes(obj.events)
            close_all_nodes(obj.transformation)
            for field in (:px_area, :limits)
                close(getfield(obj, field), true)
            end
            disconnect!(obj.camera)
            empty!(obj.theme)
            empty!(obj.children)
            empty!(obj.current_screens)
            return
        end
        obj
    end
end

# Just indexing into a scene gets you plot 1, plot 2 etc
Base.iterate(scene::Scene, idx = 1) = idx <= length(scene) ? (scene[idx], idx + 1) : nothing
Base.length(scene::Scene) = length(scene.plots)
Base.lastindex(scene::Scene) = length(scene.plots)
getindex(scene::Scene, idx::Integer) = scene.plots[idx]
GeometryTypes.widths(scene::Scene) = widths(to_value(pixelarea(scene)))
struct Axis end

child(scene::Scene) = Scene(scene, value(pixelarea(scene)))

"""
Creates a subscene with a pixel camera
"""
function cam2d(scene::Scene)
    sub = child(scene)
    cam2d!(sub)
    sub
end
function campixel(scene::Scene)
    sub = child(scene)
    campixel!(sub)
    sub
end

function getindex(scene::Scene, ::Type{Axis})
    for plot in scene
        isaxis(plot) && return plot
    end
    nothing
end

"""
Each argument can be named for a certain plot type `P`. Falls back to `arg1`, `arg2`, etc.
"""
function argument_names(plot::P) where P <: AbstractPlot
    argument_names(P, length(plot.converted))
end

function argument_names(::Type{<: AbstractPlot}, num_args::Integer)
    # this is called in the indexing function, so let's be a bit efficient
    ntuple(i-> Symbol("arg$i"), num_args)
end

function Base.empty!(scene::Scene)
    empty!(scene.plots)
    disconnect!(scene.camera)
    scene.limits[] = FRect3D(Vec3f0(0), Vec3f0(0))
    scene.camera_controls[] = EmptyCamera()
    empty!(scene.theme)
    merge!(scene.theme, _current_default_theme)
    empty!(scene.children)
    empty!(scene.current_screens)
end

limits(scene::Scene) = scene.limits
limits(scene::SceneLike) = scene.parent.limits

# Since we can use Combined like a scene in some circumstances, we define this alias
theme(x::SceneLike, args...) = theme(x.parent, args...)
theme(x::Scene) = x.theme
theme(x::Scene, key) = x.theme[key]
theme(x::AbstractPlot, key) = x.attributes[key]
theme(::Nothing, key::Symbol) = current_default_theme()[key]

Base.push!(scene::Combined, subscene) = nothing # Combined plots add themselves uppon creation
function Base.push!(scene::Scene, plot::AbstractPlot)
    push!(scene.plots, plot)
    plot.parent[] = scene
    for screen in scene.current_screens
        insert!(screen, scene, plot)
    end
end
function Base.push!(scene::Scene, plot::Combined)
    push!(scene.plots, plot)
    for screen in scene.current_screens
        insert!(screen, scene, plot)
    end
end

function connect!(scene::Scene, child::Scene)

end

function Base.push!(scene::Scene, child::Scene)
    push!(scene.children, child)
    disconnect!(child.camera)
    nodes = map([:view, :projection, :projectionview, :resolution, :eyeposition]) do field
        lift(getfield(scene.camera, field)) do val
            push!(getfield(child.camera, field), val)
        end
    end
    cameracontrols!(child, nodes)
end

events(scene::Scene) = scene.events
events(scene::SceneLike) = events(scene.parent)

camera(scene::Scene) = scene.camera
camera(scene::SceneLike) = camera(scene.parent)

cameracontrols(scene::Scene) = scene.camera_controls[]
cameracontrols(scene::SceneLike) = cameracontrols(scene.parent)

cameracontrols!(scene::Scene, cam) = (scene.camera_controls[] = cam)
cameracontrols!(scene::SceneLike, cam) = cameracontrols!(parent(scene), cam)

pixelarea(scene::Scene) = scene.px_area
pixelarea(scene::SceneLike) = pixelarea(scene.parent)

plots(scene::SceneLike) = scene.plots

const _forced_update_scheduled = Ref(false)
function must_update()
    val = _forced_update_scheduled[]
    _forced_update_scheduled[] = false
    val
end
function force_update!()
    _forced_update_scheduled[] = true
end

const current_global_scene = Ref{Any}()

if Sys.iswindows()
    function _primary_resolution()
        # ccall((:GetSystemMetricsForDpi, :user32), Cint, (Cint, Cuint), 0, ccall((:GetDpiForSystem, :user32), Cuint, ()))
        # ccall((:GetSystemMetrics, :user32), Cint, (Cint,), 17)
        dc = ccall((:GetDC, :user32), Ptr{Cvoid}, (Ptr{Cvoid},), C_NULL)
        ntuple(2) do i
            Int(ccall((:GetDeviceCaps, :gdi32), Cint, (Ptr{Cvoid}, Cint), dc, (2 - i) + 117))
        end
    end
else
    # TODO implement osx + linux
    _primary_resolution() = (1920, 1080) # everyone should have at least a hd monitor :D
end

function primary_resolution()
    # Since this is pretty low level and os specific + we can't test on all possible
    # computers, I assume we'll have bugs here. Let's not sweat about it too much,
    # we just use primary_resolution to have a good estimate for a default window resolution
    # if this fails, only thing happening will be a too small/big window when the user doesn't give any resolution.
    try
        _primary_resolution()
    catch e
        warn("Could not retrieve primary monitor resolution. A default resolution of (1920, 1080) is assumed!
        Error: $(sprint(io->showerror(io, e))).")
        (1920, 1080)
    end
end
reasonable_resolution() = primary_resolution() .÷ 2

function current_scene()
    if isassigned(current_global_scene)
        current_global_scene[]
    else
        Scene()
    end
end

Scene(::Nothing) = Scene()

const minimal_default = Attributes(
    font = "Dejavu Sans",
    backgroundcolor = RGBAf0(1,1,1,1),
    color = :black,
    colormap = :viridis,
    resolution = reasonable_resolution(),
    visible = true
)

const _current_default_theme = copy(minimal_default)

function current_default_theme(; kw_args...)
    copy = Attributes(_current_default_theme...)
    merge(copy, Attributes(;kw_args...))
end

function set_theme!(new_theme::Attributes)
    empty!(_current_default_theme)
    merge!(_current_default_theme, minimal_default, new_theme)
    return
end
function set_theme!(;kw_args...)
    set_theme!(Attributes(; kw_args...))
end

function Scene(;
        kw_args...
    )
    events = Events()
    theme = current_default_theme(; kw_args...)
    resolution = theme[:resolution][]
    px_area = foldp(IRect(0, 0, resolution), events.window_area) do v0, w_area
        wh = widths(w_area)
        wh = (wh == Vec(0, 0)) ? widths(v0) : wh
        IRect(0, 0, wh)
    end
    scene = Scene(
        events,
        px_area,
        Camera(px_area),
        RefValue{Any}(EmptyCamera()),
        node(:scene_limits, FRect3D(Vec3f0(0), Vec3f0(1))),
        Transformation(),
        AbstractPlot[],
        theme,
        Scene[],
        AbstractScreen[]
    )
    current_global_scene[] = scene
    scene
end

function Scene(
        scene::Scene;
        events = scene.events,
        px_area = scene.px_area,
        cam = scene.camera,
        camera_controls = scene.camera_controls,
        boundingbox = Node(AABB(Vec3f0(0), Vec3f0(1))),
        transformation = Transformation(scene),
        theme = Theme(),
        current_screens = scene.current_screens
    )
    child = Scene(
        events,
        px_area,
        cam,
        camera_controls,
        boundingbox,
        transformation,
        AbstractPlot[],
        merge(current_default_theme(), theme),
        Scene[],
        current_screens
    )
    push!(scene.children, child)
    child
end

function Scene(scene::Scene, area)
    events = scene.events
    px_area = signal_convert(Signal{IRect2D}, area)
    child = Scene(
        events,
        px_area,
        Camera(px_area),
        RefValue{Any}(EmptyCamera()),
        node(:scene_limits, FRect3D(Vec3f0(0), Vec3f0(1))),
        Transformation(),
        AbstractPlot[],
        copy(current_default_theme()),
        Scene[],
        scene.current_screens
    )
    push!(scene.children, child)
    child
end

"""
Fetches all plots sharing the same camera
"""
plots_from_camera(scene::Scene) = plots_from_camera(scene, scene.camera)
function plots_from_camera(scene::Scene, camera::Camera, list = AbstractPlot[])
    append!(list, scene.plots)
    for child in scene.children
        child.camera == camera && plots_from_camera(child, camera, list)
    end
    list
end

"""
Flattens all the combined plots and returns a Vector of Atomic plots
"""
function flatten_combined(plots::Vector, flat = AbstractPlot[])
    for elem in plots
        if (elem isa Combined)
            flatten_combined(elem.plots, flat)
        else
            push!(flat, elem)
        end
    end
    flat
end

function insertplots!(screen::AbstractDisplay, scene::Scene)
    for elem in scene.plots
        insert!(screen, scene, elem)
    end
    foreach(s-> insertplots!(screen, s), scene.children)
end
update_cam!(scene::Scene, bb::AbstractCamera, rect) = nothing

function center!(scene::Scene, padding = 0.01)
    bb = boundingbox(scene)
    w = widths(bb)
    padd = w .* padding
    bb = FRect3D(minimum(bb) .- padd, w .+ 2padd)
    update_cam!(scene, bb)
    force_update!()
    scene
end
parent_scene(x::Combined) = parent_scene(parent(x))
parent_scene(x::Scene) = x

Base.isopen(x::SceneLike) = events(x).window_open[]


text:
@testset "haskey" begin
    p = Properties()
    @test haskey(p, "key1") == false
    setprop(p, "key1", "value1")
    @test haskey(p, "key1") == true
end

@testset "keys&values" begin
    p = Properties()
    @test length(keys(p)) == 0
    @test length(values(p)) == 0
    setprop(p, "key1", "val1")
    setprop(p, "key2", "val2")
    pkeys = keys(p)
    pvals = values(p)
    @test getprop.(Ref(p), pkeys) == collect(pvals)
    @test length(pkeys) == 2
end

@testset "pairs" begin
    p = Properties()
    @test pairs(p) == Dict{String,String}()
    @test length(pairs(p)) == 0

    setprop(p, "key1", "val1")
    setprop(p, "key2", "val2")
    @test length(pairs(p)) == 2
    @test get.(Ref(pairs(p)), keys(p), "def") == collect(values(p))

    for (k, v) in pairs(p)
        @test haskey(p, k)
    end
end

@testset "others" begin
    # this test only makes sure, that the function can run,
    # functionality is not tested
    p = Properties()
    @test p == sizehint!(p, 15)

    @test keytype(p) == String
    @test valtype(p) == String
end

@testset "delete!" begin
    p = Properties()
    setprop(p, "key1", "val1")
    setprop(p, "key2", "val2")
    @test getprop(p, "key1") == "val1"
    @test getprop(p, "key2") == "val2"

    @test p == delete!(p, "key1")
    @test getprop(p, "key1") == nothing

    @test getprop(p, "key2") == "val2"
    @test length(keys(p)) == 1
end


text:
module TestCustomSequence

using PSCOPF

using Test
using Dates
using DataStructures

@testset verbose=true "test_custom_sequence" begin

    @testset "test_custom_sequence" begin
        network = PSCOPF.Networks.Network()
        TS = PSCOPF.create_target_timepoints(DateTime("2015-01-01T11:00:00"))
        ECH = [DateTime("2015-01-01T07:00:00"),
                DateTime("2015-01-01T10:00:00"),
                DateTime("2015-01-01T10:30:00")]
        schedule_history = Vector{PSCOPF.Schedule}()
        uncertainties = PSCOPF.Uncertainties(
            DateTime("2015-01-01T07:00:00") => PSCOPF.UncertaintiesAtEch(),
            DateTime("2015-01-01T10:00:00") => PSCOPF.UncertaintiesAtEch(),
            DateTime("2015-01-01T10:30:00") => PSCOPF.UncertaintiesAtEch()
            )

        struct MockMarket <: PSCOPF.AbstractMarket
        end
        function PSCOPF.run(runnable::MockMarket, ech, firmness, TS, context::PSCOPF.AbstractContext)
            return nothing
        end
        #affects_market_schedule default to true cause <:AbstractMarket
        function PSCOPF.update_market_schedule!(context::PSCOPF.AbstractContext, ech, result, firmness, runnable::MockMarket)
            market_schedule = PSCOPF.get_market_schedule(context)
            market_schedule.decider_type = PSCOPF.DeciderType(runnable)
            market_schedule.decision_time = ech
            push!(schedule_history, deepcopy(market_schedule))
        end

        struct MockTSO <: PSCOPF.AbstractTSO
        end
        function PSCOPF.run(runnable::MockTSO, ech, firmness, TS, context::PSCOPF.AbstractContext)
            return nothing
        end
        #affects_tso_schedule default to true cause <:AbstractTSO
        function PSCOPF.update_tso_schedule!(context::PSCOPF.AbstractContext, ech, result, firmness,  runnable::MockTSO)
            tso_schedule = PSCOPF.get_tso_schedule(context)
            tso_schedule.decider_type = PSCOPF.DeciderType(runnable)
            tso_schedule.decision_time = ech
            push!(schedule_history, deepcopy(tso_schedule))
        end
        #affects_tso_actions_schedule default to true cause <:AbstractTSO
        function PSCOPF.update_tso_actions!(context::PSCOPF.AbstractContext, ech, result, firmness, runnable::MockTSO)
            nothing
        end

        sequence = PSCOPF.Sequence(SortedDict(
            ECH[1]     => [MockMarket(), MockTSO()],
            ECH[2]     => [MockTSO()],
            ECH[3]     => [MockTSO()]
        ))

        mode = PSCOPF.ManagementMode("test_sequencing", Dates.Minute(0))
        exec_context = PSCOPF.PSCOPFContext(network, TS, mode,
                                            SortedDict{String,PSCOPF.GeneratorState}(), uncertainties)

        PSCOPF.run!(exec_context, sequence, check_context=false)

        @test length(schedule_history) == 4 # one for each executed step

        @test PSCOPF.is_market(schedule_history[1].decider_type)
        @test schedule_history[1].decision_time == ECH[1]
        @test PSCOPF.is_tso(schedule_history[2].decider_type)
        @test schedule_history[2].decision_time == ECH[1]
        @test PSCOPF.is_tso(schedule_history[3].decider_type)
        @test schedule_history[3].decision_time == ECH[2]
        @test PSCOPF.is_tso(schedule_history[4].decider_type)
        @test schedule_history[4].decision_time == ECH[3]
    end

end

end #module

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
# Calculate time derivative
function rhs!(dg::Dg2D, t_stage, uses_mpi::Val{true})
  # Start to receive MPI data
  @timeit timer() "start MPI receive" start_mpi_receive!(dg)

  # Prolong solution to MPI interfaces
  @timeit timer() "prolong2mpiinterfaces" prolong2mpiinterfaces!(dg)

  # Start to send MPI data
  @timeit timer() "start MPI send" start_mpi_send!(dg)

  # Reset u_t
  @timeit timer() "reset ∂u/∂t" dg.elements.u_t .= 0

  # Calculate volume integral
  @timeit timer() "volume integral" calc_volume_integral!(dg)

  # Prolong solution to interfaces
  @timeit timer() "prolong2interfaces" prolong2interfaces!(dg)

  # Calculate interface fluxes
  @timeit timer() "interface flux" calc_interface_flux!(dg)

  # Prolong solution to boundaries
  @timeit timer() "prolong2boundaries" prolong2boundaries!(dg)

  # Calculate boundary fluxes
  @timeit timer() "boundary flux" calc_boundary_flux!(dg, t_stage)

  # Prolong solution to mortars
  @timeit timer() "prolong2mortars" prolong2mortars!(dg)

  # Calculate mortar fluxes
  @timeit timer() "mortar flux" calc_mortar_flux!(dg)

  # Finish to receive MPI data
  @timeit timer() "finish MPI receive" finish_mpi_receive!(dg)

  # Calculate MPI interface fluxes
  @timeit timer() "MPI interface flux" calc_mpi_interface_flux!(dg)

  # Calculate surface integrals
  @timeit timer() "surface integral" calc_surface_integral!(dg)

  # Apply Jacobian from mapping to reference element
  @timeit timer() "Jacobian" apply_jacobian!(dg)

  # Calculate source terms
  @timeit timer() "source terms" calc_sources!(dg, dg.source_terms, t_stage)

  # Finish to send MPI data
  @timeit timer() "finish MPI send" finish_mpi_send!(dg)
end

# Count the number of MPI interfaces that need to be created
function count_required_mpi_interfaces(mesh::TreeMesh2D, cell_ids)
  count = 0

  # Iterate over all cells
  for cell_id in cell_ids
    for direction in 1:n_directions(mesh.tree)
      # If no neighbor exists, current cell is small or at boundary and thus we need a mortar
      if !has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # Skip if neighbor has children
      neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id]
      if has_children(mesh.tree, neighbor_cell_id)
        continue
      end

      # Skip if neighbor is on this rank -> create regular interface instead
      if mpi_isparallel() && is_own_cell(mesh.tree, neighbor_cell_id)
        continue
      end

      count += 1
    end
  end

  return count
end

# Create MPI interface container, initialize interface data, and return interface container for further use
function init_mpi_interfaces(cell_ids, mesh::TreeMesh2D, ::Val{NVARS}, ::Val{POLYDEG}, elements) where {NVARS, POLYDEG}
  # Initialize container
  n_mpi_interfaces = count_required_mpi_interfaces(mesh, cell_ids)
  mpi_interfaces = MpiInterfaceContainer2D{NVARS, POLYDEG}(n_mpi_interfaces)

  # Connect elements with interfaces
  init_mpi_interface_connectivity!(elements, mpi_interfaces, mesh)

  return mpi_interfaces
end

function start_mpi_receive!(dg::Dg2D)
  for (index, d) in enumerate(dg.mpi_neighbor_ranks)
    dg.mpi_recv_requests[index] = MPI.Irecv!(dg.mpi_recv_buffers[index], d, d, mpi_comm())
  end
end

# Initialize connectivity between elements and interfaces
function init_mpi_interface_connectivity!(elements, mpi_interfaces, mesh::TreeMesh2D)
  # Reset interface count
  count = 0

  # Iterate over all elements to find neighbors and to connect via mpi_interfaces
  for element_id in 1:nelements(elements)
    # Get cell id
    cell_id = elements.cell_ids[element_id]

    # Loop over directions
    for direction in 1:n_directions(mesh.tree)
      # If no neighbor exists, current cell is small and thus we need a mortar
      if !has_neighbor(mesh.tree, cell_id, direction)
        continue
      end

      # Skip if neighbor has children
      neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id]
      if has_children(mesh.tree, neighbor_cell_id)
        continue
      end

      # Skip if neighbor is on this MPI rank -> create regular interface instead
      if mpi_isparallel() && is_own_cell(mesh.tree, neighbor_cell_id)
        continue
      end

      # Create interface between elements
      count += 1
      mpi_interfaces.local_element_ids[count] = element_id

      if direction in (2, 4) # element is "left" of interface, remote cell is "right" of interface
        mpi_interfaces.remote_sides[count] = 2
      else
        mpi_interfaces.remote_sides[count] = 1
      end

      # Set orientation (x -> 1, y -> 2)
      if direction in (1, 2) # x-direction
        mpi_interfaces.orientations[count] = 1
      else # y-direction
        mpi_interfaces.orientations[count] = 2
      end
    end
  end

  @assert count == nmpiinterfaces(mpi_interfaces) ("Actual interface count ($count) does not match "
                                                   * "expectations $(nmpiinterfaces(mpi_interfaces))")
end

# Initialize connectivity between MPI neighbor ranks
function init_mpi_neighbor_connectivity(elements, mpi_interfaces, mesh::TreeMesh2D)
  tree = mesh.tree

  # Determine neighbor ranks and sides for MPI interfaces
  neighbor_ranks = fill(-1, nmpiinterfaces(mpi_interfaces))
  # The global interface id is the smaller of the (globally unique) neighbor cell ids, multiplied by
  # number of directions (2 * ndims) plus direction minus one
  global_interface_ids = fill(-1, nmpiinterfaces(mpi_interfaces))
  for interface_id in 1:nmpiinterfaces(mpi_interfaces)
    orientation = mpi_interfaces.orientations[interface_id]
    remote_side = mpi_interfaces.remote_sides[interface_id]
    # Direction is from local cell to remote cell
    if orientation == 1 # MPI interface in x-direction
      if remote_side == 1 # remote cell on the "left" of MPI interface
        direction = 1
      else # remote cell on the "right" of MPI interface
        direction = 2
      end
    else # MPI interface in y-direction
      if remote_side == 1 # remote cell on the "left" of MPI interface
        direction = 3
      else # remote cell on the "right" of MPI interface
        direction = 4
      end
    end
    local_element_id = mpi_interfaces.local_element_ids[interface_id]
    local_cell_id = elements.cell_ids[local_element_id]
    remote_cell_id = tree.neighbor_ids[direction, local_cell_id]
    neighbor_ranks[interface_id] = tree.mpi_ranks[remote_cell_id]
    if local_cell_id < remote_cell_id
      global_interface_ids[interface_id] = 2 * ndims(tree) * local_cell_id + direction - 1
    else
      global_interface_ids[interface_id] = (2 * ndims(tree) * remote_cell_id +
                                            opposite_direction(direction) - 1)
    end
  end

  # Get sorted, unique neighbor ranks
  mpi_neighbor_ranks = unique(sort(neighbor_ranks))

  # Sort interfaces by global interface id
  p = sortperm(global_interface_ids)
  neighbor_ranks .= neighbor_ranks[p]
  interface_ids = collect(1:nmpiinterfaces(mpi_interfaces))[p]

  # For each neighbor rank, init connectivity data structures
  mpi_neighbor_interfaces = Vector{Vector{Int}}(undef, length(mpi_neighbor_ranks))
  for (index, d) in enumerate(mpi_neighbor_ranks)
    mpi_neighbor_interfaces[index] = interface_ids[findall(x->(x == d), neighbor_ranks)]
  end

  # Sanity check that we counted all interfaces exactly once
  @assert sum(length(v) for v in mpi_neighbor_interfaces) == nmpiinterfaces(mpi_interfaces)

  return mpi_neighbor_ranks, mpi_neighbor_interfaces
end

# Initialize MPI data structures
function init_mpi_data_structures(mpi_neighbor_interfaces, ::Val{NDIMS}, ::Val{NVARS},
                                  ::Val{POLYDEG}) where {NDIMS, NVARS, POLYDEG}
  data_size = NVARS * (POLYDEG + 1)^(NDIMS - 1)
  mpi_send_buffers = Vector{Vector{Float64}}(undef, length(mpi_neighbor_interfaces))
  mpi_recv_buffers = Vector{Vector{Float64}}(undef, length(mpi_neighbor_interfaces))
  for index in 1:length(mpi_neighbor_interfaces)
    mpi_send_buffers[index] = Vector{Float64}(undef, length(mpi_neighbor_interfaces[index]) * data_size)
    mpi_recv_buffers[index] = Vector{Float64}(undef, length(mpi_neighbor_interfaces[index]) * data_size)
  end

  mpi_send_requests = Vector{MPI.Request}(undef, length(mpi_neighbor_interfaces))
  mpi_recv_requests = Vector{MPI.Request}(undef, length(mpi_neighbor_interfaces))

  return mpi_send_buffers, mpi_recv_buffers, mpi_send_requests, mpi_recv_requests
end

function prolong2mpiinterfaces!(dg::Dg2D)
  equation = equations(dg)

  Threads.@threads for s in 1:dg.n_mpi_interfaces
    local_element_id = dg.mpi_interfaces.local_element_ids[s]
    if dg.mpi_interfaces.orientations[s] == 1 # interface in x-direction
      if dg.mpi_interfaces.remote_sides[s] == 1 # local element in positive direction
        for j in 1:nnodes(dg), v in 1:nvariables(dg)
          dg.mpi_interfaces.u[2, v, j, s] = dg.elements.u[v,          1, j, local_element_id]
        end
      else # local element in negative direction
        for j in 1:nnodes(dg), v in 1:nvariables(dg)
          dg.mpi_interfaces.u[1, v, j, s] = dg.elements.u[v, nnodes(dg), j, local_element_id]
        end
      end
    else # interface in y-direction
      if dg.mpi_interfaces.remote_sides[s] == 1 # local element in positive direction
        for i in 1:nnodes(dg), v in 1:nvariables(dg)
          dg.mpi_interfaces.u[2, v, i, s] = dg.elements.u[v, i,          1, local_element_id]
        end
      else # local element in negative direction
        for i in 1:nnodes(dg), v in 1:nvariables(dg)
          dg.mpi_interfaces.u[1, v, i, s] = dg.elements.u[v, i, nnodes(dg), local_element_id]
        end
      end
    end
  end
end

function start_mpi_send!(dg::Dg2D)
  data_size = nvariables(dg) * nnodes(dg)^(ndims(dg) - 1)

  for d in 1:length(dg.mpi_neighbor_ranks)
    send_buffer = dg.mpi_send_buffers[d]

    for (index, s) in enumerate(dg.mpi_neighbor_interfaces[d])
      first = (index - 1) * data_size + 1
      last =  (index - 1) * data_size + data_size

      if dg.mpi_interfaces.remote_sides[s] == 1 # local element in positive direction
        @views send_buffer[first:last] .= vec(dg.mpi_interfaces.u[2, :, :, s])
      else # local element in negative direction
        @views send_buffer[first:last] .= vec(dg.mpi_interfaces.u[1, :, :, s])
      end
    end
  end

  # Start sending
  for (index, d) in enumerate(dg.mpi_neighbor_ranks)
    dg.mpi_send_requests[index] = MPI.Isend(dg.mpi_send_buffers[index], d, mpi_rank(), mpi_comm())
  end
end

function finish_mpi_receive!(dg::Dg2D)
  data_size = nvariables(dg) * nnodes(dg)^(ndims(dg) - 1)

  # Start receiving and unpack received data until all communication is finished
  d, _ = MPI.Waitany!(dg.mpi_recv_requests)
  while d != 0
    recv_buffer = dg.mpi_recv_buffers[d]

    for (index, s) in enumerate(dg.mpi_neighbor_interfaces[d])
      first = (index - 1) * data_size + 1
      last =  (index - 1) * data_size + data_size

      if dg.mpi_interfaces.remote_sides[s] == 1 # local element in positive direction
        @views vec(dg.mpi_interfaces.u[1, :, :, s]) .= recv_buffer[first:last]
      else # local element in negative direction
        @views vec(dg.mpi_interfaces.u[2, :, :, s]) .= recv_buffer[first:last]
      end
    end

    d, _ = MPI.Waitany!(dg.mpi_recv_requests)
  end
end

# Calculate and store the surface fluxes (standard Riemann and nonconservative parts) at an MPI interface
# OBS! Regarding the nonconservative terms: 1) currently only needed for the MHD equations
#                                           2) not implemented for MPI
calc_mpi_interface_flux!(dg::Dg2D) = calc_mpi_interface_flux!(dg.elements.surface_flux_values,
                                                              have_nonconservative_terms(dg.equations),
                                                              dg)

function calc_mpi_interface_flux!(surface_flux_values, nonconservative_terms::Val{false}, dg::Dg2D)
  @unpack surface_flux_function = dg
  @unpack u, local_element_ids, orientations, remote_sides = dg.mpi_interfaces

  Threads.@threads for s in 1:dg.n_mpi_interfaces
    # Get local neighboring element
    element_id = local_element_ids[s]

    # Determine interface direction with respect to element:
    if orientations[s] == 1 # interface in x-direction
      if remote_sides[s] == 1 # local element in positive direction
        direction = 1
      else # local element in negative direction
        direction = 2
      end
    else # interface in y-direction
      if remote_sides[s] == 1 # local element in positive direction
        direction = 3
      else # local element in negative direction
        direction = 4
      end
    end

    for i in 1:nnodes(dg)
      # Call pointwise Riemann solver
      u_ll, u_rr = get_surface_node_vars(u, dg, i, s)
      flux = surface_flux_function(u_ll, u_rr, orientations[s], equations(dg))

      # Copy flux to local element storage
      for v in 1:nvariables(dg)
        surface_flux_values[v, i, direction, element_id] = flux[v]
      end
    end
  end
end

function finish_mpi_send!(dg::Dg2D)
  MPI.Waitall!(dg.mpi_send_requests)
end

function analyze_solution(dg::Dg2D, mesh::TreeMesh, time, dt, step, runtime_absolute,
                          runtime_relative, uses_mpi::Val{true}; solver_gravity=nothing)
  equation = equations(dg)

  # General information
  mpi_println()
  mpi_println("-"^80)
  mpi_println(" Simulation running '", get_name(equation), "' with POLYDEG = ", polydeg(dg))
  mpi_println("-"^80)
  mpi_println(" #timesteps:     " * @sprintf("% 14d", step) *
              "               " *
              " run time:       " * @sprintf("%10.8e s", runtime_absolute))
  mpi_println(" dt:             " * @sprintf("%10.8e", dt) *
              "               " *
              " PID:            " * @sprintf("%10.8e s", runtime_relative))
  mpi_println(" sim. time:      " * @sprintf("%10.8e", time) *
              "               " *
              " PID × #ranks:   " * @sprintf("%10.8e s", runtime_relative * mpi_nranks()))

  # Level information (only show for AMR) #TODO MPI add when AMR is enabled
  # if parameter("amr_interval", 0)::Int > 0 && mpi_isroot()
  #   levels = Vector{Int}(undef, dg.n_elements)
  #   for element_id in 1:dg.n_elements
  #     levels[element_id] = mesh.tree.levels[dg.elements.cell_ids[element_id]]
  #   end
  #   min_level = minimum(levels)
  #   max_level = maximum(levels)

  #   mpi_println(" #elements:      " * @sprintf("% 14d", dg.n_elements))
  #   for level = max_level:-1:min_level+1
  #     mpi_println(" ├── level $level:    " * @sprintf("% 14d", count(x->x==level, levels)))
  #   end
  #   mpi_println(" └── level $min_level:    " * @sprintf("% 14d", count(x->x==min_level, levels)))
  # end
  # mpi_println()

  # Open file for appending and store time step and time information
  if dg.save_analysis && mpi_isroot()
    f = open(dg.analysis_filename, "a")
    @printf(f, "% 9d", step)
    @printf(f, "  %10.8e", time)
    @printf(f, "  %10.8e", dt)
  end

  # Calculate and print derived quantities (error norms, entropy etc.)
  # Variable names required for L2 error, Linf error, and conservation error
  if mpi_isroot()
    if any(q in dg.analysis_quantities for q in
          (:l2_error, :linf_error, :conservation_error, :residual))
      print(" Variable:    ")
      for v in 1:nvariables(equation)
        @printf("   %-14s", varnames_cons(equation)[v])
      end
      println()
    end
  end

  # Calculate L2/Linf errors, which are also returned by analyze_solution
  l2_error, linf_error = calc_error_norms(dg, time, Val(true))

  if mpi_isroot()
    # L2 error
    if :l2_error in dg.analysis_quantities
      print(" L2 error:    ")
      for v in 1:nvariables(equation)
        @printf("  % 10.8e", l2_error[v])
        dg.save_analysis && @printf(f, "  % 10.8e", l2_error[v])
      end
      println()
    end

    # Linf error
    if :linf_error in dg.analysis_quantities
      print(" Linf error:  ")
      for v in 1:nvariables(equation)
        @printf("  % 10.8e", linf_error[v])
        dg.save_analysis && @printf(f, "  % 10.8e", linf_error[v])
      end
      println()
    end
  end

  # Conservation errror
  if :conservation_error in dg.analysis_quantities
    # Calculate state integrals
    state_integrals = integrate(dg.elements.u, dg)

    # Store initial state integrals at first invocation
    if isempty(dg.initial_state_integrals)
      dg.initial_state_integrals = zeros(nvariables(equation))
      if mpi_isroot()
        # Only set on MPI root; all other ranks do not get any value from `integrate`
        dg.initial_state_integrals .= state_integrals
      end
    end

    if mpi_isroot()
      print(" |∑U - ∑U₀|:  ")
      for v in 1:nvariables(equation)
        err = abs(state_integrals[v] - dg.initial_state_integrals[v])
        @printf("  % 10.8e", err)
        dg.save_analysis && @printf(f, "  % 10.8e", err)
      end
      println()
    end
  end

  # Residual (defined here as the vector maximum of the absolute values of the time derivatives)
  if :residual in dg.analysis_quantities
    mpi_print(" max(|Uₜ|):   ")
    for v in 1:nvariables(equation)
      # Calculate maximum absolute value of Uₜ
      res = maximum(abs, view(dg.elements.u_t, v, :, :, :))
      res = MPI.Reduce!(Ref(res), max, mpi_root(), mpi_comm())
      mpi_isroot() && @printf("  % 10.8e", res[])
      mpi_isroot() && dg.save_analysis && @printf(f, "  % 10.8e", res[])
    end
    mpi_println()
  end

  # L2/L∞ errors of the primitive variables
  if :l2_error_primitive in dg.analysis_quantities || :linf_error_primitive in dg.analysis_quantities
    l2_error_prim, linf_error_prim = calc_error_norms(cons2prim, dg, time, Val(true))

    if mpi_isroot()
      print(" Variable:    ")
      for v in 1:nvariables(equation)
        @printf("   %-14s", varnames_prim(equation)[v])
      end
      println()

      # L2 error
      if :l2_error_primitive in dg.analysis_quantities
        print(" L2 error prim.: ")
        for v in 1:nvariables(equation)
          @printf("%10.8e   ", l2_error_prim[v])
          dg.save_analysis && @printf(f, "  % 10.8e", l2_error_prim[v])
        end
        println()
      end

      # L∞ error
      if :linf_error_primitive in dg.analysis_quantities
        print(" Linf error pri.:")
        for v in 1:nvariables(equation)
          @printf("%10.8e   ", linf_error_prim[v])
          dg.save_analysis && @printf(f, "  % 10.8e", linf_error_prim[v])
        end
        println()
      end
    end
  end

  # Entropy time derivative
  if :dsdu_ut in dg.analysis_quantities
    dsdu_ut = calc_entropy_timederivative(dg, time)
    if mpi_isroot()
      print(" ∑∂S/∂U ⋅ Uₜ: ")
      @printf("  % 10.8e", dsdu_ut)
      dg.save_analysis && @printf(f, "  % 10.8e", dsdu_ut)
      println()
    end
  end

  # Entropy
  if :entropy in dg.analysis_quantities
    s = integrate(dg, dg.elements.u) do i, j, element_id, dg, u
      cons = get_node_vars(u, dg, i, j, element_id)
      return entropy(cons, equations(dg))
    end
    if mpi_isroot()
      print(" ∑S:          ")
      @printf("  % 10.8e", s)
      dg.save_analysis && @printf(f, "  % 10.8e", s)
      println()
    end
  end

  # Total energy
  if :energy_total in dg.analysis_quantities
    e_total = integrate(dg, dg.elements.u) do i, j, element_id, dg, u
      cons = get_node_vars(u, dg, i, j, element_id)
      return energy_total(cons, equations(dg))
    end
    if mpi_isroot()
      print(" ∑e_total:    ")
      @printf("  % 10.8e", e_total)
      dg.save_analysis && @printf(f, "  % 10.8e", e_total)
      println()
    end
  end

  # Kinetic energy
  if :energy_kinetic in dg.analysis_quantities
    e_kinetic = integrate(dg, dg.elements.u) do i, j, element_id, dg, u
      cons = get_node_vars(u, dg, i, j, element_id)
      return energy_kinetic(cons, equations(dg))
    end
    if mpi_isroot()
      print(" ∑e_kinetic:  ")
      @printf("  % 10.8e", e_kinetic)
      dg.save_analysis && @printf(f, "  % 10.8e", e_kinetic)
      println()
    end
  end

  # Internal energy
  if :energy_internal in dg.analysis_quantities
    e_internal = integrate(dg, dg.elements.u) do i, j, element_id, dg, u
      cons = get_node_vars(u, dg, i, j, element_id)
      return energy_internal(cons, equations(dg))
    end
    if mpi_isroot()
      print(" ∑e_internal: ")
      @printf("  % 10.8e", e_internal)
      dg.save_analysis && @printf(f, "  % 10.8e", e_internal)
      println()
    end
  end

  # Magnetic energy #TODO MPI add when MHD is enabled
  # if :energy_magnetic in dg.analysis_quantities
  #   e_magnetic = integrate(dg, dg.elements.u) do i, j, element_id, dg, u
  #     cons = get_node_vars(u, dg, i, j, element_id)
  #     return energy_magnetic(cons, equations(dg))
  #   end
  #   if mpi_isroot()
  #     print(" ∑e_magnetic: ")
  #     @printf("  % 10.8e", e_magnetic)
  #     dg.save_analysis && @printf(f, "  % 10.8e", e_magnetic)
  #     println()
  #   end
  # end

  # Potential energy #TODO MPI add when Euler-gravity is enabled
  # if :energy_potential in dg.analysis_quantities
  #   # FIXME: This should be implemented properly for multiple coupled solvers
  #   @assert !isnothing(solver_gravity) "Only works if gravity solver is supplied"
  #   @assert dg.initial_conditions == initial_conditions_jeans_instability "Only works with Jeans instability setup"

  #   e_potential = integrate(dg, dg.elements.u, solver_gravity.elements.u) do i, j, element_id, dg, u_euler, u_gravity
  #     cons_euler = get_node_vars(u_euler, dg, i, j, element_id)
  #     cons_gravity = get_node_vars(u_gravity, solver_gravity, i, j, element_id)
  #     # OBS! subtraction is specific to Jeans instability test where rho_0 = 1.5e7
  #     return (cons_euler[1] - 1.5e7) * cons_gravity[1]
  #   end
  #   if mpi_isroot()
  #     print(" ∑e_pot:      ")
  #     @printf("  % 10.8e", e_potential)
  #     dg.save_analysis && @printf(f, "  % 10.8e", e_potential)
  #     println()
  #   end
  # end

  # Solenoidal condition ∇ ⋅ B = 0 #TODO MPI add when MHD is enabled
  # if :l2_divb in dg.analysis_quantities || :linf_divb in dg.analysis_quantities
  #   l2_divb, linf_divb = calc_mhd_solenoid_condition(dg, time)
  # end
  # if mpi_isroot()
  #   # L2 norm of ∇ ⋅ B
  #   if :l2_divb in dg.analysis_quantities
  #     print(" L2 ∇ ⋅B:     ")
  #     @printf("  % 10.8e", l2_divb)
  #     dg.save_analysis && @printf(f, "  % 10.8e", l2_divb)
  #     println()
  #   end
  #   # Linf norm of ∇ ⋅ B
  #   if :linf_divb in dg.analysis_quantities
  #     print(" Linf ∇ ⋅B:   ")
  #     @printf("  % 10.8e", linf_divb)
  #     dg.save_analysis && @printf(f, "  % 10.8e", linf_divb)
  #     println()
  #   end
  # end

  # Cross helicity #TODO MPI add when MHD is enabled
  # if :cross_helicity in dg.analysis_quantities
  #   h_c = integrate(dg, dg.elements.u) do i, j, element_id, dg, u
  #     cons = get_node_vars(u, dg, i, j, element_id)
  #     return cross_helicity(cons, equations(dg))
  #   end
  #   if mpi_isroot()
  #     print(" ∑H_c:        ")
  #     @printf("  % 10.8e", h_c)
  #     dg.save_analysis && @printf(f, "  % 10.8e", h_c)
  #     println()
  #   end
  # end

  if mpi_isroot()
    println("-"^80)
    println()

    # Add line break and close analysis file if it was opened
    if dg.save_analysis
      println(f)
      close(f)
    end
  end

  # Return errors for EOC analysis
  return l2_error, linf_error
end

# OBS! Global results are only calculated on MPI root
@inline calc_error_norms(dg::Dg2D, t, uses_mpi) = calc_error_norms(cons2cons, dg, t, uses_mpi)
function calc_error_norms(func, dg::Dg2D, t, uses_mpi::Val{true})
  l2_error, linf_error = calc_error_norms(func, dg, t, Val(false))

  # Since the local L2 norm is already normalized and square-rooted, we need to undo this first
  global_l2_error = Vector(l2_error.^2 .* dg.analysis_total_volume)
  global_linf_error = Vector(linf_error)
  MPI.Reduce!(global_l2_error, +, mpi_root(), mpi_comm())
  MPI.Reduce!(global_linf_error, max, mpi_root(), mpi_comm())
  l2_error = convert(typeof(l2_error), global_l2_error)
  linf_error = convert(typeof(linf_error), global_linf_error)

  l2_error = @. sqrt(l2_error / dg.analysis_total_volume)

  return l2_error, linf_error
end

#TODO MPI add when MHD is enabled
# function calc_mhd_solenoid_condition(dg::Dg2D, t, mpi_parallel::Val{true})
#   l2_divb, linf_divb = calc_mhd_solenoid_condition(func, dg, t, Val(false))
#
#   # Since the local L2 norm is already normalized and square-rooted, we need to undo this first
#   global_l2_divb = Vector(l2_divb.^2 .* dg.analysis_total_volume)
#   global_linf_divb = Vector(linf_divb)
#   MPI.Reduce!(global_l2_divb, +, mpi_root(), mpi_comm())
#   MPI.Reduce!(global_linf_divb, max, mpi_root(), mpi_comm())
#   l2_divb = convert(typeof(l2_divb), global_l2_divb)
#   linf_divb = convert(typeof(linf_divb), global_linf_divb)
#
#   l2_divb = @. sqrt(l2_divb / dg.analysis_total_volume)
#
#   return l2_divb, linf_divb
# end

# OBS! Global results are only calculated on MPI root, all other domains receive `nothing`
function integrate(func, dg::Dg2D, uses_mpi::Val{true}, args...; normalize=true)
  integral = integrate(func, dg, Val(false), args...; normalize=normalize)
  integral = MPI.Reduce!(Ref(integral), +, mpi_root(), mpi_comm())

  return mpi_isroot() ? integral[] : integral
end

# Calculate stable time step size
function calc_dt(dg::Dg2D, cfl, uses_mpi::Val{true})
  min_dt = calc_dt(dg, cfl, Val(false))
  min_dt = MPI.Allreduce!(Ref(min_dt), min, mpi_comm())[]

  return min_dt
end


text:
export take_until

import Base: show

"""
    take_until(notifier::S)

Creates a take operator, which returns an Observable
that emits the values emitted by the source Observable until a `notifier` Observable emits a value or a completion event.

# Arguments
- `notifier::S`: The Observable whose first emitted value will cause the output Observable of `take_until` to stop emitting values from the source Observable.

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples

```julia
using Rocket

source = interval(100) |> take_until(timer(1000))

subscribe!(source, logger())
;

# output

[LogActor] Data: 0
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Data: 4
[LogActor] Data: 5
[LogActor] Data: 6
[LogActor] Data: 7
[LogActor] Data: 8
[LogActor] Completed
```

```julia
using Rocket

source = interval(100)

subscribe!(source |> take_until(source |> filter(i -> i == 3)), logger())
;

# output

[LogActor] Data: 0
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Completed
```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
take_until(notifier::N) where N = TakeUntilOperator{N}(notifier)

struct TakeUntilOperator{N} <: InferableOperator
    notifier :: N
end

function on_call!(::Type{L}, ::Type{L}, operator::TakeUntilOperator{N}, source) where { L, N }
    return proxy(L, source, TakeUntilProxy{N}(operator.notifier))
end

operator_right(operator::TakeUntilOperator, ::Type{L}) where L = L

struct TakeUntilProxy{N} <: SourceProxy
    notifier :: N
end

source_proxy!(::Type{L}, proxy::TakeUntilProxy{N}, source::S) where { L, N, S } = TakeUntilSource{L, N, S}(proxy.notifier, source)

mutable struct TakeUntilInnerActor{L, A} <: Actor{L}
    actor         :: A
    isdisposed    :: Bool
    ssubscription :: Teardown
    nsubscription :: Teardown

    TakeUntilInnerActor{L, A}(actor::A) where { L, A } = new(actor, false, voidTeardown, voidTeardown)
end

function on_next!(actor::TakeUntilInnerActor{L}, data::L) where L
    if !actor.isdisposed
        next!(actor.actor, data)
    end
end

function on_error!(actor::TakeUntilInnerActor, err)
    if !actor.isdisposed
        __dispose(actor)
        error!(actor.actor, err)
    end
end

function on_complete!(actor::TakeUntilInnerActor)
    if !actor.isdisposed
        __dispose(actor)
        complete!(actor.actor)
    end
end

function __dispose(actor::TakeUntilInnerActor)
    actor.isdisposed = true
    unsubscribe!(actor.ssubscription)
    unsubscribe!(actor.nsubscription)
end

struct TakeUntilInnerActorGuard{I} <: Actor{Any}
    inner :: I
end

on_next!(guard::TakeUntilInnerActorGuard, _)    = complete!(guard.inner)
on_error!(guard::TakeUntilInnerActorGuard, err) = error!(guard.inner, err)
on_complete!(guard::TakeUntilInnerActorGuard)   = complete!(guard.inner)

@subscribable struct TakeUntilSource{L, N, S} <: Subscribable{L}
    notifier :: N
    source   :: S
end

function on_subscribe!(observable::TakeUntilSource{L}, actor::A) where { L, A }
    inner = TakeUntilInnerActor{L, A}(actor)
    guard = TakeUntilInnerActorGuard(inner)

    inner.nsubscription = subscribe!(observable.notifier, guard)

    if !inner.isdisposed
        inner.ssubscription = subscribe!(observable.source, inner)
        # Extra check for cases when __dispose has been called during main subscription
        if inner.isdisposed
            unsubscribe!(inner.ssubscription)
        end
    end

    return TakeUntilSubscription(inner)
end

struct TakeUntilSubscription{A} <: Teardown
    inner :: A
end

as_teardown(::Type{<:TakeUntilSubscription}) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::TakeUntilSubscription)
    __dispose(subscription.inner)
    return nothing
end

Base.show(io::IO, ::TakeUntilOperator)              = print(io, "TakeUntilOperator()")
Base.show(io::IO, ::TakeUntilProxy)                 = print(io, "TakeUntilProxy()")
Base.show(io::IO, ::TakeUntilInnerActor{L}) where L = print(io, "TakeUntilInnerActor($L)")
Base.show(io::IO, ::TakeUntilSource{L})     where L = print(io, "TakeUntilSource($L)")
Base.show(io::IO, ::TakeUntilSubscription)          = print(io, "TakeUntilSubscription()")


text:
const AMD_DBGAPI_SYMBOL_VERSION = "AMD_DBGAPI_0.19"
const GlobalAddress = UInt64
const Notifier = Cint
const ClientProcessID = Ptr{Cvoid}
const QueuePacketID = UInt64
const WatchpointID = UInt32
const LaneID = UInt32
const SegmentAddress = UInt64
const ClientThreadID = Ptr{Cvoid}

# Perhaps useful to know macro definitions for object handle values...
# AMD_DBGAPI_ARCHITECTURE_NONE ( amd_dbgapi_architecture_id_t { 0 } )
# AMD_DBGAPI_PROCESS_NONE ( amd_dbgapi_process_id_t { 0 } )
# AMD_DBGAPI_CODE_OBJECT_NONE ( amd_dbgapi_code_object_id_t { 0 } )
# AMD_DBGAPI_AGENT_NONE ( amd_dbgapi_agent_id_t { 0 } )
# AMD_DBGAPI_QUEUE_NONE ( amd_dbgapi_queue_id_t { 0 } )
# AMD_DBGAPI_DISPATCH_NONE ( amd_dbgapi_dispatch_id_t { 0 } )
# AMD_DBGAPI_WAVE_NONE ( amd_dbgapi_wave_id_t { 0 } )
# AMD_DBGAPI_DISPLACED_STEPPING_NONE ( amd_dbgapi_displaced_stepping_id_t { 0 } )
# AMD_DBGAPI_WATCHPOINT_NONE ( amd_dbgapi_watchpoint_id_t ( - 1 ) )
# AMD_DBGAPI_REGISTER_CLASS_NONE ( amd_dbgapi_register_class_id_t { 0 } )
# AMD_DBGAPI_REGISTER_NONE ( amd_dbgapi_register_id_t { 0 } )
# AMD_DBGAPI_LANE_NONE ( amd_dbgapi_lane_id_t ( - 1 ) )
# AMD_DBGAPI_ADDRESS_CLASS_NONE ( amd_dbgapi_address_class_id_t { 0 } )
# AMD_DBGAPI_EVENT_NONE ( amd_dbgapi_event_id_t { 0 } )
# AMD_DBGAPI_SHARED_LIBRARY_NONE ( amd_dbgapi_shared_library_id_t { 0 } )
# AMD_DBGAPI_BREAKPOINT_NONE ( amd_dbgapi_breakpoint_id_t { 0 } )


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

@testset "PureBasis" begin

##
@info("--------- Testing PureBasis ----------")
using LinearAlgebra, Test
using ACE
using ACE.OrthPolys: OrthPolyBasis, OneOrthogonal
using JuLIP: evaluate, evaluate_d, evaluate!, evaluate_d!

##

@info("Generate many-body 1-orth basis")

deg = 10  # degree
r0, rcut = 0.5, 2.5
trans = PolyTransform(2, r0)
t0, tcut = trans(r0), trans(rcut)
tdf = range(t0, tcut, length=1000)
ww = ones(length(tdf))
J = OrthPolyBasis(deg,  0, trans(r0), 2, trans(rcut), tdf, ww)
P = ACE.TransformedPolys(J, trans, r0, rcut)

spec = SparseSHIP(2, deg-1; wL = 1.5)
aceB = SHIPBasis(spec, P)

pureB = ACE.PureBasis(aceB)
pureB1 = ACE.PureBasis(aceB, 1)
pureB2 = ACE.PureBasis(aceB, 2)

##

Rs = ACE.rand_vec(aceB.J, 2)
Zs = zeros(Int16, 2)
b = evaluate(aceB, Rs, Zs, 0)
bp = pureB2(Rs)
norm(bp)
# @assert length(b) == length(aceB) == length(bp)

##

function condtest(basis, nneigs;
                  nsamples = 10_000,
                  randfun = ()->ACE.rand_vec(aceB.J, nneigs) )
   G = zeros(length(basis), length(basis))
   for _ = 1:nsamples
      Rs = randfun()
      b = basis(Rs)
      G += b * b' / nsamples
   end
   D = Diagonal( diag(G).^(-0.5) )
   Gscal = D * G * D
   return cond(Gscal)
end

##

@info("Conditioning of a pure 1N basis, 1 neig")
@show condtest(pureB1, 1)

@info("Conditioning of a pure 1N basis, 5 neighbours")
@show condtest(pureB1, 5)

@info("Conditioning of a pure 2N basis, 2 neigs")
@show condtest(pureB2, 2; nsamples=100_000)

@info("Conditioning of a pure 2N basis, 5 neigs")
@show condtest(pureB2, 5; nsamples=100_000)

@info("Conditioning of the full basis, 2 neigs")
@show condtest(pureB, 2; nsamples=100_000)

@info("Conditioning of the full basis, 5 neigs")
@show condtest(pureB, 5; nsamples=100_000)

##

# now we try the same with a 1-orthogonal basis
@info("Construct 1-orthogonal basis")
J1 = OneOrthogonal(J)
P1 = ACE.TransformedPolys(J1, trans, r0, rcut)
spec = SparseSHIP(2, deg-1; wL = 1.5,
                  filterfcn = ν -> ACE.OrthPolys.filter_oneorth(ν, J1))
aceBoo = SHIPBasis(spec, P1)

pureB1oo = ACE.PureBasis(aceBoo, 1)
pureB2oo = ACE.PureBasis(aceBoo, 2)
pureBoo = ACE.PureBasis(aceBoo)

##

@info("Conditioning of One-Orth pure 1N basis, 1 neig")
@show condtest(pureB1oo, 1)

@info("Conditioning of  One-Orth pure 1N basis, 5 neigs")
@show condtest(pureB1oo, 5)

@info("Conditioning of  One-Orth pure 2N basis, 2 neigs")
@show condtest(pureB2oo, 2; nsamples=100_000)

@info("Conditioning of  One-Orth pure 2N basis, 5 neigs")
@show condtest(pureB2oo, 5; nsamples=100_000)

@info("Conditioning of the full One-Orth basis, 2 neigs")
@show condtest(pureBoo, 2; nsamples=100_000)

@info("Conditioning of the full One-Orth basis, 5 neigs")
@show condtest(pureBoo, 5; nsamples=100_000)

##

end


text:
function selectbyproperty(pkproperties, pmin, pmax)
    keep = trues(length(pkproperties))
    if pmin !== nothing
        keep .&= (pmin .<= pkproperties)
    end
    if pmax !== nothing
        keep .&= (pkproperties .<= pmax)
    end
    keep
end

function selectbypeakdistance(pkindices, priority, distance)
    npkindices = length(pkindices)
    keep = trues(npkindices)

    prioritytoposition = fsortperm(priority)
    for i in npkindices:-1:1
        j = prioritytoposition[i]
        (keep[j] == 0) && continue

        k = j-1
        while (1 <= k) && ((pkindices[j]-pkindices[k]) < distance)
            keep[k] = 0
            k -= 1
        end

        k = j+1
        while (k<=npkindices) && ((pkindices[k]-pkindices[j]) < distance)
            keep[k] = 0
            k += 1
        end
    end
    keep
end

function argwlenasexpected(value)
    if value === nothing
        value = -1
    elseif 1 < value
        value = ceil(Int, value)
    else
        throw(ArgumentError("`wlen` must be larger than 1, was $(value)"))
    end
    value
end


text:
using Makie

function main()
    x = range(0, stop=2π, length=40)
    scene = lines(x,sin,color=:blue)
    axis = scene[Axis]
    axis.grid.linecolor=((:green, 1), (:blue, 0.5))
    axis[:names].textcolor = ((:red, 1.0), (:blue, 1.0))
    axis[:names][:axisnames] = ("x", "y=cos(x)")
    scene
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# GLFW Gamepad functions
# https://www.glfw.org/docs/latest/input_guide.html#gamepad

const AXIS_LEFT_X       = 1
const AXIS_LEFT_Y       = 2
const AXIS_RIGHT_X      = 3
const AXIS_RIGHT_Y      = 4

const BUTTON_DPAD_UP    = 12
const BUTTON_DPAD_RIGHT = 13
const BUTTON_DPAD_DOWN  = 14
const BUTTON_DPAD_LEFT  = 15

struct GLFWgamepadstate
	buttons::SVector{15, UInt8}  # unsigned char buttons[15];    // GLFW_PRESS or GLFW_RELEASE
	axes::SVector{6, Float32}    # float axes[6];                // -1.0 to 1.0 inclusive
end

#= GetGamepadState
This function retrieves the state of the specified joystick remapped to an Xbox-like gamepad.
If the specified joystick is not present or does not have a gamepad mapping this function will return GLFW_FALSE but will not generate an error. Call glfwJoystickPresent to check whether it is present regardless of whether it has a mapping.
Not all devices have all the buttons or axes provided by GLFWgamepadstate. Unavailable buttons and axes will always report GLFW_RELEASE and 0.0 respectively.
Possible errors include GLFW_NOT_INITIALIZED and GLFW_INVALID_ENUM.
Added in version 3.3.
=#

GetGamepadState(joy::GLFW.Joystick, state) = Bool(ccall((:glfwGetGamepadState, GLFW.libglfw), Cint, (Cint, Ref{GLFWgamepadstate}), joy, state))

function handle_scalar_settings(s::mjSim, buttons, axes)
    config  = s.robot.controller.config
    command = s.robot.command

    # Velocity
    command.horizontal_velocity[1] = axes[AXIS_LEFT_Y] * -config.max_x_velocity
    command.horizontal_velocity[2] = axes[AXIS_LEFT_X] * -config.max_y_velocity

    # Yaw
    command.yaw_rate = axes[AXIS_RIGHT_X] * -config.max_yaw_rate

    # Pitch
    set_pitch!(s.robot, axes[AXIS_RIGHT_Y])

    # Roll
    dpadx = Int(buttons[BUTTON_DPAD_RIGHT]) - Int(buttons[BUTTON_DPAD_LEFT])
    dpadx != 0 && adjust_roll!(s.robot, dpadx)

    # Height
    dpady = Int(buttons[BUTTON_DPAD_UP]) - Int(buttons[BUTTON_DPAD_DOWN])
    dpady != 0 && adjust_height!(s.robot, dpady)
end

prev_buttons = @SVector fill(0x00, 15)
button_commands = SVector{15, RobotCmd}(
    CYCLE_HOP,          # GLFW_GAMEPAD_BUTTON_CROSS
    NO_COMMAND, NO_COMMAND, NO_COMMAND,
    TOGGLE_ACTIVATION,  # GLFW_GAMEPAD_BUTTON_LEFT_BUMPER
    TOGGLE_TROT,        # GLFW_GAMEPAD_BUTTON_RIGHT_BUMPER
    NO_COMMAND, NO_COMMAND, NO_COMMAND, NO_COMMAND,
    NO_COMMAND, NO_COMMAND, NO_COMMAND, NO_COMMAND, NO_COMMAND)

function handle_behavior_state_change(s::mjSim, buttons)
    triggered = buttons .& (prev_buttons .⊻ buttons)
    global prev_buttons = deepcopy(buttons)

    for (doit, robotcmd) in zip(triggered, button_commands)
        if Bool(doit)
            # println("Executing Robot Command: $(repr(robotcmd))")
            execute_robotcmd(s, robotcmd)
        end
    end
end

function gamepad(s::mjSim, joy::GLFW.Joystick)
    state = Ref{GLFWgamepadstate}()

    if GetGamepadState(joy, state)
        buttons = state[].buttons
        axes = state[].axes

        # println("Gamepad state: $(repr(state))")
        #message = "Gamepad state: $buttons $(map((x) -> round(x; digits=2), axes))"
        #send(s.socket, message)

        handle_behavior_state_change(s, buttons)
        handle_scalar_settings(s, buttons, axes)
    end
end


text:
"""
    sequencingRRAB(setting::RRBSetting, Union{ExactArc,SimpleArc}, packet)

Call the correct rrb method based on the settings

# Arguments
- `setting::RRBSetting`: Which function to call
- `arc::Union{ExactArc,SimpleArc}`: The current arc
- `packet`: Data with the necesary info
"""
function sequencingRRAB(setting::RRBSetting, arc::Union{ExactArc,SimpleArc}, packet)
    if setting == sop
        return sopRRAB(arc,packet[1],packet[2])
    end
end


text:
using Test
using Statistics
using RoME

@testset "Testing tree init prior usage" begin

fg = generateCanonicalFG_Circle(4;graphinit=true)
# deleteVariable!.(fg, [:x3, :x4, :l1])

# fg.solverParams.showtree = true
# fg.solverParams.drawtree = true
# fg.solverParams.dbg = true
fg.solverParams.useMsgLikelihoods = true

# now move the graph prior somewhere else (ie all init is wrong)
deleteFactor!(fg, :x0f1)
prpo = PriorPose2(MvNormal([5.,0.,0.0], 0.01*Matrix{Float64}(LinearAlgebra.I,3,3)))
addFactor!(fg, [:x0], prpo)

solveTree!(fg);

X4 = getBelief(fg, :x4) |> getPoints

@error("Must first fix IIF #913")
# @test_broken 2.0 < Statistics.mean(X4[1,:])
@test -1.5 < Statistics.mean(X4[2,:]) < 1.5

end


text:
module Sample

using ..Space, ..Util
import Random

export defrandalg, condomegasample, randsample

"Default sampling algorithm"
function defrandalg end

"`n` samples from `ω::ΩT` such that `yₓ(ω)` is true where `yₓ` are conditions of `x`"
omegasample(rng, ΩT, x, n; kwargs) =
  randsample(rng, ΩT, conditions(x), n; kwargs...)

"`n` samples from `ω::ΩT` such that `y(ω)` is true"
condomegasample(rng, ΩT, y, n; alg = defrandalg(rng, x, n), kwargs...) =
  map

"""
`n` random samples from `x` such that `yₓ(ω)` is true where `yₓ` are conditions of `x
using `alg` algorithm

```
using Distributions
x = Normal(0, 1)
randsample(x |ᶜ x >ₛ 2.0, 3; alg = RejectionSample)
'''
"""
randsample(rng::Random.AbstractRNG, x, n; alg = defrandalg(x), ΩT = defΩ(), kwargs...) =
  randsample(rng, ΩT, x, n, alg; kwargs...)

# Convenience methods

"`randsample` using GLOBAL_RNG"
randsample(x, n::Integer; kwargs...) =
  randsample(Random.GLOBAL_RNG, x, n; kwargs...)

randsample(x; kwargs...) =
  first(randsample(x, 1, kwargs...))

randsample(rng::Random.AbstractRNG, xs::Tuple, args...; kwargs...) =
  randsample(rng, ω -> mapf(ω, xs), args...; kwargs...)

# There are two different notions of ΩT
# The first is the return type of the omegas
# Second is Omegas used

end


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

mutable struct IoK8sApiAuditregistrationV1alpha1AuditSinkSpec <: SwaggerModel
    policy::Any # spec type: Union{ Nothing, IoK8sApiAuditregistrationV1alpha1Policy } # spec name: policy
    webhook::Any # spec type: Union{ Nothing, IoK8sApiAuditregistrationV1alpha1Webhook } # spec name: webhook

    function IoK8sApiAuditregistrationV1alpha1AuditSinkSpec(;policy=nothing, webhook=nothing)
        o = new()
        validate_property(IoK8sApiAuditregistrationV1alpha1AuditSinkSpec, Symbol("policy"), policy)
        setfield!(o, Symbol("policy"), policy)
        validate_property(IoK8sApiAuditregistrationV1alpha1AuditSinkSpec, Symbol("webhook"), webhook)
        setfield!(o, Symbol("webhook"), webhook)
        o
    end
end # type IoK8sApiAuditregistrationV1alpha1AuditSinkSpec

const _property_map_IoK8sApiAuditregistrationV1alpha1AuditSinkSpec = Dict{Symbol,Symbol}(Symbol("policy")=>Symbol("policy"), Symbol("webhook")=>Symbol("webhook"))
const _property_types_IoK8sApiAuditregistrationV1alpha1AuditSinkSpec = Dict{Symbol,String}(Symbol("policy")=>"IoK8sApiAuditregistrationV1alpha1Policy", Symbol("webhook")=>"IoK8sApiAuditregistrationV1alpha1Webhook")
Base.propertynames(::Type{ IoK8sApiAuditregistrationV1alpha1AuditSinkSpec }) = collect(keys(_property_map_IoK8sApiAuditregistrationV1alpha1AuditSinkSpec))
Swagger.property_type(::Type{ IoK8sApiAuditregistrationV1alpha1AuditSinkSpec }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAuditregistrationV1alpha1AuditSinkSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiAuditregistrationV1alpha1AuditSinkSpec }, property_name::Symbol) =  _property_map_IoK8sApiAuditregistrationV1alpha1AuditSinkSpec[property_name]

function check_required(o::IoK8sApiAuditregistrationV1alpha1AuditSinkSpec)
    (getproperty(o, Symbol("policy")) === nothing) && (return false)
    (getproperty(o, Symbol("webhook")) === nothing) && (return false)
    true
end

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


text:
# This script actually solves the TextbookNK model with a risk-adjusted linearization
# and times the methods, if desired
using RiskAdjustedLinearizations, JLD2, LinearAlgebra, Test
include("textbook_nk.jl")
out = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "..", "test", "reference", "textbook_nk_ss_output.jld2"), "r")

# Settings
autodiff              = false
algorithm             = :relaxation
euler_equation_errors = false
test_price_dispersion = false # check if price dispersion in steady state is always bounded below by 1

# Set up
m_nk = TextbookNK() # create parameters
m = textbook_nk(m_nk) # instantiate risk-adjusted linearization
autodiff_method = autodiff ? :forward : :central

# Solve!
solve!(m; algorithm = :deterministic, autodiff = autodiff_method)
@test m.z ≈ out["z_det"]
@test m.y ≈ out["y_det"]
@test m.Ψ ≈ out["Psi_det"]
z_det = copy(m.z)
y_det = copy(m.y)
Ψ_det = copy(m.Ψ)

solve!(m; algorithm = algorithm, autodiff = autodiff_method)
@test m.z ≈ out["z"]
@test m.y ≈ out["y"]
@test m.Ψ ≈ out["Ψ"]

if test_price_dispersion
    π̃_ss_vec = log.(range(1 - .005, stop = 1 + .005, length = 10)) # On annualized basis, range from -2% to 2% target inflation
    det_soln = Dict()
    sss_soln = Vector{RiskAdjustedLinearization}(undef, length(π̃_ss_vec))

    for (i, π̃_ss) in enumerate(π̃_ss_vec)
        local m_nk = TextbookNK(; π̃_ss = π̃_ss)
        local m = textbook_nk(m_nk)
        solve!(m; algorithm = :deterministic, verbose = :none)
        det_soln[i] = Dict()
        det_soln[i][:z] = copy(m.z)
        det_soln[i][:y] = copy(m.y)
        det_soln[i][:Ψ] = copy(m.Ψ)
        solve!(m; algorithm = algorithm, verbose = :none)
        sss_soln[i] = m
    end

    det_v = exp.([det_soln[i][:z][3] for i in 1:length(det_soln)])
    sss_v = exp.([sss_soln[i].z[3] for i in 1:length(sss_soln)])
    @test all(det_v .> 1.)
    @test all(sss_v .> 1.)
end

if euler_equation_errors
    # Load shocks. Using CRW ones b/c that model also has 2 standard normal random variables
    shocks = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "..", "test", "reference", "crw_shocks.jld2"), "r")["shocks"]

    # With this simple model, the Euler equation for bonds holds exactly
    @test abs(euler_equation_error(m, nk_cₜ, (a, b, c, d) -> nk_logSDFxR(a, b, c, d; β = m_nk.β, σ = m_nk.σ),
                                   nk_𝔼_quadrature, shocks, summary_statistic = x -> norm(x, Inf))) ≈ 0.
end

nothing


text:
# ------------------------------------------------------------------------------------------
# ### Notebook 1
# #### 1.1
# Look up docs for the `convert` function.
# ------------------------------------------------------------------------------------------

?convert

# ------------------------------------------------------------------------------------------
# #### 1.2
# Assign `365` to a variable named `days`. Convert `days` to a float.
# ------------------------------------------------------------------------------------------

days = 365
convert(Float64, days)

# ------------------------------------------------------------------------------------------
# #### 1.3
# See what happens when you execute
#
# ```julia
# convert(Int64, '1')
# ```
# and
#
# ```julia
# parse(Int64, '1')
# ```
#
# What's the difference?
# ------------------------------------------------------------------------------------------

convert(Int64, '1')

parse(Int64, '1')

# ------------------------------------------------------------------------------------------
# ### Notebook 2
# #### 2.1
# Create a string that says "hi" 1000 times.
# ------------------------------------------------------------------------------------------

"hi"^1000

# ------------------------------------------------------------------------------------------
# #### 2.2
# Declare two variables
#
# ```julia
# a = 3
# b = 4
# ```
# and use them to create two strings:
# ```julia
# "3 + 4"
# "7"
# ```
# ------------------------------------------------------------------------------------------

a = 3; b = 4

"$a + $b"

"$(a + b)"

# ------------------------------------------------------------------------------------------
# ### Notebook 3
#
# #### 3.1
# Create an array, `a_ray`, with the following code:
#
# ```julia
# a_ray = [1, 2, 3]
# ```
#
# Add the number `4` to the end of this array and then remove it.
# ------------------------------------------------------------------------------------------

a_ray = [1, 2, 3]
push!(a_ray, 4)

# ------------------------------------------------------------------------------------------
# #### 3.2
# Try to add "Emergency" as key to `myphonebook` with the value `string(911)` with the
# following code
# ```julia
# myphonebook["Emergency"] = 911
# ```
#
# Why doesn't this work?
# ------------------------------------------------------------------------------------------

myphonebook = Dict("Jenny" => "867-5309", "Ghostbusters" => "555-2368")

myphonebook["Emergency"] = 911
#=

Julia constructed `myphonebook` to match the arguments
it was passed; both of which had `String`s for both their
keys and values. We see that myphonebook is a
`Dict{String,String} with 2 entries`. This means Julia
will not accept integers as values in myphonebook.

=#

# This will work:
myphonebook["Emergency"] = "911"

# ------------------------------------------------------------------------------------------
# #### 3.3
# Create a new dictionary called `flexible_phonebook` that has Jenny's number stored as a
# string and Ghostbusters' number stored as an integer.
# ------------------------------------------------------------------------------------------

flexible_phonebook = Dict("Jenny" => "867-5309", "Ghostbusters" => 5552368)

# ------------------------------------------------------------------------------------------
# #### 3.4
# Add the key "Emergency" with the value `911` (an integer) to `flexible_phonebook`.
# ------------------------------------------------------------------------------------------

flexible_phonebook["Emergency"] = 911

# ------------------------------------------------------------------------------------------
# #### 3.5
# Why can we add an integer as a value to `flexible_phonebook` but not `myphonebook`? How
# could we have initialized `myphonebook` so that it would accept integers or strings as
# values?
# ------------------------------------------------------------------------------------------

#=

We constructed `flexible_phonebook` with two key-value
pairs whose values were different types, so Julia needed
to construct a dictionary that could hold `Any` value.
Unlike myphonebook, flexible_phonebook is a
`Dict{String,Any} with 2 entries`.

To avoid this, we could have initialized myphonebook a
an empty dictionary and added entries later. Or we
could have explicitly told Julia that we wanted a
dictionary that accepted objects of type `Any` as
values. See examples!
=#

myphonebook = Dict()

# Alternatively we could use a parametric type constructor to tell Julia explicitly
# what types we want to populate our dictionary with.
myphonebook = Dict{String, Any}("Jenny" => "867-5309", "Ghostbusters" => "555-2368")

# ------------------------------------------------------------------------------------------
# ### Notebook 4
#
# #### 4.1
# Write a loop that prints the squares of integers between 1 and 100.
# ------------------------------------------------------------------------------------------

for i in 1:100
   println(i^2)
end

# ------------------------------------------------------------------------------------------
# #### 4.2
#
# Add to the code above a bit to create a dictionary, `squares` that holds integers and
# their squares as key, value pairs such that
#
# ```julia
# squares[10] == 100
# ```
# ------------------------------------------------------------------------------------------

squares = Dict()
for i in 1:100
   squares[i] = i^2
end
squares

# ------------------------------------------------------------------------------------------
# #### 4.3
# Use an array comprehension to create an an array that stores the squares of all integers
# between 1 and 100.
# ------------------------------------------------------------------------------------------

[i^2 for i in 1:100]

# ------------------------------------------------------------------------------------------
# ### Notebook 5
#
# #### 5.1
# Write a conditional statement that prints a number if the number is even and the string
# "odd" if the number is odd.
# ------------------------------------------------------------------------------------------

n = 3
if (n % 2) == 0
    println(n)
else
    println("odd")
end

# ------------------------------------------------------------------------------------------
# #### 5.2
# Rewrite the code from 5.1 using a ternary operator.
# ------------------------------------------------------------------------------------------

(n % 2) == 0 ? println(n) : println("odd")

# ------------------------------------------------------------------------------------------
# ### Notebook 6
#
# #### 6.1
# Write a function that adds 1 to its input.
# ------------------------------------------------------------------------------------------

x -> x + 1

add1(x) = x + 1

# ------------------------------------------------------------------------------------------
# #### 6.2
# Use `map` or `broadcast` to increment every element of matrix `A` by `1`.
# ------------------------------------------------------------------------------------------

A = rand(5,5)
map(x -> x + 1, A)

broadcast(add1, A)

# ------------------------------------------------------------------------------------------
# #### 6.3
# Use the broadcast dot syntax to increment every element of matrix `A` by `1`.
# ------------------------------------------------------------------------------------------

add1.(A)

A .+ 1

# ------------------------------------------------------------------------------------------
# ### Notebook 7
#
# #### 7.1
# Load the Primes package (source code at https://github.com/JuliaMath/Primes.jl).
# ------------------------------------------------------------------------------------------

# using Pkg; Pkg.add("Primes")
using Primes

# ------------------------------------------------------------------------------------------
# #### 7.2
# Verify that you can now use the function `primes` to grab all prime numbers under
# 1,000,000.
# ------------------------------------------------------------------------------------------

primes(1000000)

# ------------------------------------------------------------------------------------------
# ### Notebook 8
#
# #### 8.1
# Given
# ```julia
# x = -10:10
# ```
# plot y vs. x for $y = x^2$.
# ------------------------------------------------------------------------------------------

using Plots; gr()
x = -10:10
y = x .^ 2
plot(x, y)

plot(x -> x ^ 2, x)

# ------------------------------------------------------------------------------------------
# #### 8.2
# Execute the following code
# ------------------------------------------------------------------------------------------

p1 = plot(x, x)
p2 = plot(x, x.^2)
p3 = plot(x, x.^3)
p4 = plot(x, x.^4)
plot(p1,p2,p3,p4,layout=(2,2),legend=false)

# ------------------------------------------------------------------------------------------
# and then create a $4x1$ plot that uses `p1`, `p2`, `p3`, and `p4` as subplots.
# ------------------------------------------------------------------------------------------

p1 = plot(x, x)
p2 = plot(x, x.^2)
p3 = plot(x, x.^3)
p4 = plot(x, x.^4)
plot(p1,p2,p3,p4,layout=(4,1),legend=false)

# ------------------------------------------------------------------------------------------
# ### Notebook 9
#
# #### 9.1
#
# Extend the function `foo`, adding a method that takes only one input argument, which is of
# type `Bool`, and prints "foo with one boolean!"
# ------------------------------------------------------------------------------------------

foo(x::Bool) = println("foo with one boolean!")

# ------------------------------------------------------------------------------------------
# #### 9.2
#
# Check that the method being dispatched when you execute
# ```julia
# foo(true)
# ```
# is the one you wrote.
# ------------------------------------------------------------------------------------------

foo(true)

@which foo(true)

# ------------------------------------------------------------------------------------------
# ### Notebook 10
#
# #### 10.1
# Take the inner product (or "dot" product) of a vector `v` with itself.
# ------------------------------------------------------------------------------------------

v = [1, 2, 3]
v' * v

using LinearAlgebra
dot(v, v)

# ------------------------------------------------------------------------------------------
# #### 10.2
# Take the outer product of a vector v with itself.
# ------------------------------------------------------------------------------------------

v * v'

# ------------------------------------------------------------------------------------------
# ### Notebook 11
#
# #### 11.1
#
# What are the eigenvalues of matrix A?
#
# ```
# A =
# [
#  140   97   74  168  131
#   97  106   89  131   36
#   74   89  152  144   71
#  168  131  144   54  142
#  131   36   71  142   36
# ]
# ```
# ------------------------------------------------------------------------------------------

A =
[
 140   97   74  168  131
  97  106   89  131   36
  74   89  152  144   71
 168  131  144   54  142
 131   36   71  142   36
]

eigdec = eigen(A)
eigdec.values

# ------------------------------------------------------------------------------------------
# #### 11.2
# Create a `Diagonal` matrix from the eigenvalues of `A`.
# ------------------------------------------------------------------------------------------

Diagonal(eigdec.values)

# ------------------------------------------------------------------------------------------
# #### 11.3
# Create a `LowerTriangular` matrix from `A`.
# ------------------------------------------------------------------------------------------

LowerTriangular(A)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
function conjugatefrankwolfe(network::AbstractNetwork,
                             trips::Matrix,
                             costfn::Function;
                             basedon=:link,
                             errtol=1e-4,
                             δ=1e-3)
    # initialize
    flows = allornothing(network, trips, costfn; basedon=basedon)

    # start iteration
    err = 1.
    oldflows = copy(flows)
    while err > errtol
        ## find target solution
        tmp = costfn(flows; returnitems=[:costs, :derivs])
        linkcosts, H = tmp[:costs], Diagonal(tmp[:derivs])
        shortflows = allornothing(network, trips, linkcosts; basedon=basedon)

        newdir = reduce(+, shortflows - flows, dims=2:length(flows))
        olddir = reduce(+, oldflows - flows, dims=2:length(flows))
        α = (olddir' * H * newdir)/(olddir' * H * (newdir - olddir))
        α = isnan(α) ? 0. : (α < 0.) ? 0. : (α > (1. - δ)) ? (1. - δ) : α
        conjflows = α * oldflows + (1-α) * shortflows

        ## find stepsize
        conjdirfull = conjflows - flows
        μ = bisection((μ) -> sum(costfn(flows + μ * conjdirfull).costs .* conjdirfull))

        ## calculate new link travel times
        flows += μ * conjdirfull

        ## calculate error
        err = sum(flows .* linkcosts)/sum(shortflows .* linkcosts) - 1.
        oldflows = conjflows
    end
    return (flows, err)
end


text:
using PlotMesh
using Test
using VisualRegressionTests

istravis = "TRAVIS" ∈ keys(ENV)

function testfun(fname)
    plotmesh("test_geom.h5";outputfn=fname)
end

@testset "PlotMesh.jl" begin
     #@test test_images(VisualTest(testfun, "test_geom.png"), popup=!istravis) |> save_comparison |> success
     @visualtest testfun "test_geom.png" false #!istravis
end


text:
using SimpleBibTeX
using Test

@testset "SimpleBibTeX.jl" begin
    bibstr = """
        @article{Rudin:1997,
            author    = {Rudin, S. and Dyakonov, M.},
            title     = {Edge and strip plasmons in a two-dimensional electron fluid},
            journal   = {Phys. Rev. B},
            year      = {1997},
            volume    = {55},
            pages     = {4684},
            doi       = {10.1103/PhysRevB.55.4684},
            timestamp = {2015.07.03}
        }
        @article{Abajo:1997,
            author    = {Garc\\'\\i{}a de Abajo, {\\relax F.J} and Aizpurua, J.},
            title     = {Numerical simulation of electron energy loss near inhomogeneous dielectrics},
            journal   = {Phys. Rev.~B},
            year      = {1997},
            volume    = {56},
            pages     = {15873-15884},
            doi       = {10.1103/PhysRevB.56.15873},
            timestamp = {2015.08.19}
        }
        """

    B, preamble = parsebibtex(bibstr, isfilename=false)
    C = B["Abajo:1997"]
    # test basic functionality
    @test length(B) == 2
    @test B["Abajo:1997"]["volume"] == "56"
    @test B["Rudin:1997"]["pages"] == "4684"
    @test sort(collect(keys(B))) == ["Abajo:1997", "Rudin:1997"]
    @test preamble === ""

    # test show methods
    display_string(x) = (io=IOBuffer(); show(io, MIME"text/plain"(), x); String(take!(io)))

    strB = display_string(B)
    strC =  display_string(C)
    ioC  = IOBuffer(strC)

    @test (readline(IOBuffer(strB)) ==                    # Bibliography summary
           "Bibliography with 2 Citations:" ==
           summary(B)*':')
    @test readline(ioC) == "  Abajo:1997 (Article)"       # Citation "header"
    @test readline(ioC)[1:3] == "  ┌"                     # Citation opening bracket

    # test preamble
    B′, preamble′ = parsebibtex(bibstr*"\n@preamble{A nice preamble}", isfilename=false)
    @test all(values(B) .== values(B′))
    @test keys(B) == keys(B′)
    @test preamble′ == "A nice preamble"

    # test citation keys
    println(Set(keys(C)))
    ks = ["volume","author","year","journal","pages","doi","title","timestamp"]
    @test Set(keys(C)) == Set(ks)

    # Corner-case: an balanced brackets + unmatched/unbalanced \{ + empty in a field entry
    bibstr = """
        @article{test:1,
            author = {Cuppa, {J.O}},
            title  = {A tale of an unmatched \$\\}\$ parenthesis},
            journal = {},
        }
        """

    B, _ = parsebibtex(bibstr, isfilename=false)
    C = B["test:1"]
    @test C["author"] == "Cuppa, {J.O}" # benign, matched brackets
    @test C["title"] == "A tale of an unmatched \$\\}\$ parenthesis" # unmatched bracket
    @test C["journal"] == "" # empty field
end


text:
using Test
using Dictionaries
using Indexing

include("Indices.jl")
include("Dictionary.jl")
include("ArrayIndices.jl")
include("ArrayDictionary.jl")
include("PairDictionary.jl")
include("FillDictionary.jl")
include("indexing.jl")
include("foreach.jl")
include("map.jl")
include("broadcast.jl")
include("filter.jl")
include("find.jl")
include("reverse.jl")
include("show.jl")