text:
@safetestset log_base_two_float = "log base two agrees with float version" begin
using BijectiveHilbert
for i in 1:100
    m = BijectiveHilbert.log_base2(i) + 1
    # println(i, " ", string(i, base = 2), " ", m, " ", ceil(Int, log(2, i)))
    @test((1<<(m - 1)) & i != 0)
    @test(1<<m & i == 0)
end
end

@safetestset bitrotate_left = "rotateleft cases" begin
using BijectiveHilbert
tries = [
    Any[0b1, 0b100, 2, 8],
    Any[0b1, 0b1000, 3, 8],
    Any[0b1, 0b1000, 3, 4],
    Any[0b10, 0b0001, 3, 4],
    Any[0b1000, 0b0001, 1, 4],
    Any[0b1000, 0b0010, 2, 4],
    Any[0b1110, 0b1011, 2, 4]
]
for (try_idx, trial) in enumerate(tries)
    x, y, k, n = trial
    y1 = BijectiveHilbert.rotateleft(x, k, n)
    if y != y1
        println("k ", k, " n ", n, " ", bitstring(y), " ", bitstring(y1))
    end
    @test y == y1
end
end

@safetestset bitrotate_right = "rotateright cases" begin
using BijectiveHilbert
tries = [
    Any[0b1, 0b10000000, 1, 8],
    Any[0b1, 0b01000000, 1, 7],
    Any[0b1, 0b01000000, 2, 8],
    Any[0b1, 0b00001000, 1, 4],
    Any[0b1, 0b00000100, 1, 3],
    Any[0b1, 0b00000010, 2, 3],
]
for (try_idx, trial) in enumerate(tries)
    x, y, k, n = trial
    y1 = BijectiveHilbert.rotateright(x, k, n)
    if y != y1
        println("k ", k, " n ", n, " ", bitstring(y), " ", bitstring(y1))
    end
    @test y == y1
end

end

@safetestset rotateright_random = "rotateright random" begin
    using BijectiveHilbert: fbvn1s, rotateright, rotateright_naive
    using Random
    rng = MersenneTwister(9871879147)
    for T in [UInt8, UInt16, UInt32, UInt64]
        for trial in 1:1000
            bit_max = 8 * sizeof(T)
            width = rand(rng, 2:bit_max)
            shift = rand(rng, 0:(width - 1))
            val = rand(rng, zero(T):fbvn1s(T, width))
            naive = rotateright_naive(val, shift, width)
            fancy = rotateright(val, shift, width)
            @test fancy == naive
        end
    end
end

@safetestset rotateleft_random = "rotateleft random" begin
    using BijectiveHilbert: fbvn1s, rotateleft, rotateleft_naive
    using Random
    rng = MersenneTwister(9871879147)
    for T in [UInt8, UInt16, UInt32, UInt64]
        for trial in 1:1000
            bit_max = 8 * sizeof(T)
            width = rand(rng, 2:bit_max)
            shift = rand(rng, 0:(width - 1))
            val = rand(rng, zero(T):fbvn1s(T, width))
            naive = rotateleft_naive(val, shift, width)
            fancy = rotateleft(val, shift, width)
            @test fancy == naive
        end
    end
end

@safetestset reverse_trials = "reversebits trials" begin
using BijectiveHilbert: reverse_bits
trials = [
    (0b11001, 5, 0b10011),
    (0b01001, 5, 0b10010),
    (0b00011, 5, 0b11000)
]
for (a, n, b) in trials
    r = reverse_bits(a, n)
    @test bitstring(r) == bitstring(b)
end
end

@safetestset setbits_trials = "setbits does some sets" begin
    using BijectiveHilbert: set_bits, set_bits_naive
    trials = [
        [0b0, 2, 0, 0b11, 0b11],
        [0b0, 2, 1, 0b11, 0b1100],
        [0b0, 2, 1, 0b01, 0b0100],
        [0b0, 2, 1, 0b10, 0b1000],
        [0b0, 3, 1, 0b111, 0b111000],
        [0b0, 3, 1, 0b101, 0b101000],
        [0b0, 3, 1, 0b100, 0b100000],
        [0b1, 3, 1, 0b100, 0b100001]
    ]
    for (h, n, i, w, r) in trials
        a = set_bits_naive(h, n, i, w)
        @test bitstring(a) == bitstring(r)
    end
end

@safetestset setbits_agrees = "setbits is same as hamilton" begin
    using BijectiveHilbert: set_bits, set_bits_naive
    using Random
    rng = MersenneTwister(82147985)
    for i in 1:10000
        a = rand(rng, UInt64)
        n = rand(rng, 1:7)
        i = rand(rng, 1:8)
        w = rand(rng, 0:(1<<n - 1))
        h1 = set_bits(a, n, i, w)
        h2 = set_bits_naive(a, n, i, w)
        @test h1 == h2
    end
end

@safetestset getbits_examples = "getbits trials" begin
using BijectiveHilbert: get_bits
    trials = [
        (0b1, 1, 0, 0b1),
        (0b1, 1, 1, 0b0),
        (0b1, 1, 2, 0b0),
        (0b1, 1, 7, 0b0),
        (0b1, 2, 0, 0b1),
        (0b1, 2, 1, 0b0),
        (0b100, 1, 0, 0b0),
        (0b100, 1, 1, 0b0),
        (0b100, 1, 2, 0b1),
        (0b100, 1, 3, 0b0),
        (0b100, 1, 4, 0b0),
        (0b10000000, 1, 7, 0b1),
        (0b10000000, 1, 6, 0b0),
        (0b10000000, 1, 5, 0b0),
        (0b111000, 1, 2, 0b0),
        (0b111000, 1, 3, 0b1),
        (0b111000, 1, 4, 0b1),
        (0b111000, 1, 5, 0b1),
        (0b111000, 1, 6, 0b0),
        (0b111000, 1, 7, 0b0),
        (0b111000, 2, 0, 0b0),
        (0b111000, 2, 1, 0b0),
        (0b111000, 2, 2, 0b10),
        (0b111000, 2, 3, 0b11),
        (0b111000, 2, 4, 0b11),
        (0b111000, 2, 5, 0b01),
        (0b111000, 2, 6, 0b0),
        (0b111000, 3, 0, 0b0),
        (0b111000, 3, 1, 0b100),
        (0b111000, 3, 2, 0b110),
        (0b111000, 3, 3, 0b111),
        (0b111000, 3, 4, 0b011),
        (0b111000, 3, 5, 0b001),
        (0b111000, 3, 6, 0b000)
    ]
    for (v, b, i, r) in trials
        a = get_bits(v, b, i)
        @test a == r
    end
end

@safetestset fsb_trials = "first_set_bit unsigned trials" begin
using BijectiveHilbert: first_set_bit
    trials = [
        (0b0, 0),
        (UInt64(0), 0),
        (UInt32(0), 0),
        (UInt16(0), 0),
        (UInt8(0), 0),
        (UInt64(1), 1),
        (UInt32(1), 1),
        (UInt16(1), 1),
        (UInt8(1), 1),
        (UInt128(0b1100), 3),
        (UInt64(0b1100), 3),
        (UInt32(0b1100), 3),
        (UInt16(0b1100), 3),
        (UInt8(0b1100), 3),
        (0x80, 8),
        (0x81, 1),
        (0x82, 2),
        (0x84, 3)
    ]
    for (v, i) in trials
        r = first_set_bit(v)
        @test r == i
    end
end

@safetestset fsb_unsigned = "first_set_bit signed matches signed" begin
using BijectiveHilbert: first_set_bit
using Random
trials = [Int8, Int16, Int32, Int64]
rng = MersenneTwister(349720)
for T in trials
    for i in 1:100
        v = rand(rng, T)
        a = first_set_bit(v)
        ui = parse(UInt64, "0b" * bitstring(v))
        b = first_set_bit(ui)
        @test a == b
    end
end
end


text:
type RendererState
    current::Dict{Symbol,Any}
    saved::Vector{Dict{Symbol,Any}}

    RendererState() = new(Dict{Symbol,Any}(),Dict{Symbol,Any}[])
end

function set(self::RendererState, name::Symbol, value)
    self.current[name] = value
end

function get(self::RendererState, name::Symbol, notfound=nothing)
    if haskey(self.current, name)
        return self.current[name]
    end
    for d = self.saved
        if haskey(d,name)
            return d[name]
        end
    end
    return notfound
end

function save(self::RendererState)
    unshift!(self.saved, self.current)
    self.current = Dict{Symbol,Any}()
end

function restore(self::RendererState)
    self.current = self.saved[1]
    splice!(self.saved, 1)
end

abstract Renderer

type CairoRenderer <: Renderer
    ctx::CairoContext
    state::RendererState

    function CairoRenderer(surface)
        ctx = CairoContext(surface)
        new(ctx, RendererState())
    end
end

width(r::CairoRenderer) = width(r.ctx.surface)
height(r::CairoRenderer) = height(r.ctx.surface)
boundingbox(c::CairoRenderer) = BoundingBox(0., width(c), 0., height(c))

# convert to postscipt pt = in/72
const xx2pt = @Dict( "in"=>72., "pt"=>1., "mm"=>2.835, "cm"=>28.35 )
function _str_size_to_pts(str)
    m = match(r"([\d.]+)([^\s]+)", str)
    num_xx = parse(Float64, m.captures[1])
    units = m.captures[2]
    num_pt = num_xx*xx2pt[units]
    return num_pt
end

## state commands

color_to_rgb(i::Integer) = convert(RGB, RGB24(unsigned(i)))
color_to_rgb(s::AbstractString) =parse(Colorant,s)
color_to_rgb(rgb::@compat(Tuple{Real,Real,Real})) = RGB(rgb...)
color_to_rgb(cv::Color) = convert(RGB, cv)
color_to_rgb(cv::TransparentColor) = cv

set_color(ctx::CairoContext, color) = set_source(ctx, color_to_rgb(color))

function set_clip_rect(ctx::CairoContext, bb::BoundingBox)
    rectangle(ctx, xmin(bb), ymin(bb), width(bb), height(bb))
    clip(ctx)
end

const __pl_style_func = @Dict(
    :color     => set_color,
    :linecolor => set_color,
    :fillcolor => set_color,
    :linestyle => set_line_type,
    :linekind  => set_line_type,
    :linewidth => set_line_width,
    :filltype  => set_fill_type,
    :cliprect  => set_clip_rect,
)

function set(self::CairoRenderer, key::Symbol, value)
    set(self.state, key, value)
    if key == :fontface
        fontsize = get(self, :fontsize, 12)
        set_font_face(self.ctx, "$value $(fontsize)px")
    elseif key == :fontsize
        fontface = get(self, :fontface, "sans-serif")
        set_font_face(self.ctx, "$fontface $(value)px")
    elseif haskey(__pl_style_func, key)
        __pl_style_func[key](self.ctx, value)
    end
end
set(self::CairoRenderer, key::AbstractString, value) = set(self, Symbol(key), value)

function get(self::CairoRenderer, parameter::Symbol, notfound=nothing)
    return get(self.state, parameter, notfound)
end
get(self::CairoRenderer, parameter::AbstractString, notfound=nothing) = get(self, Symbol(parameter), notfound)

function save_state(self::CairoRenderer)
    save(self.state)
    save(self.ctx)
end

function restore_state(self::CairoRenderer)
    restore(self.state)
    restore(self.ctx)
end

## drawing commands

stroke(cr::CairoRenderer) = stroke(cr.ctx)

move_to(self::CairoRenderer, px, py) = move_to(self.ctx, px, py)
line_to(self::CairoRenderer, px, py) = line_to(self.ctx, px, py)
rel_line_to(self::CairoRenderer, px, py) = rel_line_to(self.ctx, px, py)

function line(self::CairoRenderer, px, py, qx, qy)
    move_to(self.ctx, px, py)
    line_to(self.ctx, qx, qy)
    stroke(self.ctx)
end

const symbol_funcs = @Dict(
    "asterisk" => (c, x, y, r) -> (
        move_to(c, x, y+r);
        line_to(c, x, y-r);
        move_to(c, x+0.866r, y-0.5r);
        line_to(c, x-0.866r, y+0.5r);
        move_to(c, x+0.866r, y+0.5r);
        line_to(c, x-0.866r, y-0.5r)
    ),
    "circle" => (c, x, y, r) -> (
        new_sub_path(c);
        circle(c, x, y, r)
    ),
    "cross" => (c, x, y, r) -> (
        move_to(c, x+r, y+r);
        line_to(c, x-r, y-r);
        move_to(c, x+r, y-r);
        line_to(c, x-r, y+r)
    ),
    "diamond" => (c, x, y, r) -> (
        move_to(c, x, y+r);
        line_to(c, x+r, y);
        line_to(c, x, y-r);
        line_to(c, x-r, y);
        close_path(c)
    ),
    "dot" => (c, x, y, r) -> (
        new_sub_path(c);
        rectangle(c, x, y, 1., 1.)
    ),
    "plus" => (c, x, y, r) -> (
        move_to(c, x+r, y);
        line_to(c, x-r, y);
        move_to(c, x, y+r);
        line_to(c, x, y-r)
    ),
    "square" => (c, x, y, r) -> (
        new_sub_path(c);
        rectangle(c, x-0.866r, y-0.866r, 1.732r, 1.732r)
    ),
    "triangle" => (c, x, y, r) -> (
        move_to(c, x, y+r);
        line_to(c, x+0.866r, y-0.5r);
        line_to(c, x-0.866r, y-0.5r);
        close_path(c)
    ),
    "down-triangle" => (c, x, y, r) -> (
        move_to(c, x, y-r);
        line_to(c, x+0.866r, y+0.5r);
        line_to(c, x-0.866r, y+0.5r);
        close_path(c)
    ),
    "right-triangle" => (c, x, y, r) -> (
        move_to(c, x+r, y);
        line_to(c, x-0.5r, y+0.866r);
        line_to(c, x-0.5r, y-0.866r);
        close_path(c)
    ),
    "left-triangle" => (c, x, y, r) -> (
        move_to(c, x-r, y);
        line_to(c, x+0.5r, y+0.866r);
        line_to(c, x+0.5r, y-0.866r);
        close_path(c)
    ),
)

function symbols(self::CairoRenderer, x, y)
    fullname = get(self.state, :symbolkind, "circle")
    size = get(self.state, :symbolsize, 0.01)

    splitname = split(fullname)
    name = pop!(splitname)
    filled = "solid" in splitname || "filled" in splitname

    default_symbol_func = symbol_funcs["circle"]
    symbol_func = get(symbol_funcs, name, default_symbol_func)

    save(self.ctx)
    set_dash(self.ctx, Float64[])
    new_path(self.ctx)
    for i = 1:min(length(x),length(y))
        symbol_func(self.ctx, x[i], y[i], size)
    end
    if filled
        fill_preserve(self.ctx)
    end
    stroke(self.ctx)
    restore(self.ctx)
end

function curve(self::CairoRenderer, x::AbstractArray, y::AbstractArray)
    n = min(length(x), length(y))
    n > 0 || return
    new_path(self.ctx)

    lo = 1
    while lo < n
        while lo <= n && !(isfinite(x[lo]) && isfinite(y[lo]))
            lo += 1
        end

        hi = lo + 1
        while hi <= n &&  (isfinite(x[hi]) && isfinite(y[hi]))
            hi += 1
        end
        hi -= 1

        if lo < hi
            move_to(self, x[lo], y[lo])
            for i = (lo+1):hi
                line_to(self, x[i], y[i])
                if i < hi && (i & 127) == 0
                    stroke(self.ctx)
                    move_to(self, x[i], y[i])
                end
            end
            stroke(self.ctx)
        end

        lo = hi + 1
    end
end

image(r::CairoRenderer, src, x, y, w, h) = Cairo.image(r.ctx, src, x, y, w, h)

function polygon(self::CairoRenderer, points::Vector)
    polygon(self.ctx, points)
    fill(self.ctx)
end

function rectangle(self::CairoRenderer, bbox::BoundingBox, filled::Bool=true)
    rectangle(self.ctx, bbox)
    filled ? fill(self.ctx) : stroke(self.ctx)
end

function layout_text(self::CairoRenderer, str::AbstractString)
    set_latex(self.ctx, str, get(self,:fontsize))
end

function textdraw(self::CairoRenderer, x::Real, y::Real, str::AbstractString; kwargs...)
    return Cairo.text(self.ctx, x, y, set_latex(self.ctx, str, get(self,:fontsize)); markup=true, kwargs...)
end

function textwidth(self::CairoRenderer, str)
    layout_text(self, str)
    extents = get_layout_size(self.ctx)
    extents[1]
end

function textheight(self::CairoRenderer, str)
    get(self.state, :fontsize) ## XXX: kludge?
end

ext: jl


text:
# Activate local environment
using Pkg
Pkg.activate(".")

printstyled("  Running "; color = :green, bold = true)
print("model...\n")

# Load module
using Negotiations

# Run model
params = parameter_set_from_config("config.yaml")
db = load_database("db.sqlite")
model = setup_model(params, db)
rule1 = BoundedConfidence(bc = 1.0, inertia = 10.0)
rule2 = ContinuousHomophily(inertia = 20.0)
@time simulate(model, rule1, 3, db, batchname = "bounded_confidence")
@time simulate(model, rule2, 3, db, batchname = "continuous_homophily")


text:
# id = VERSION >= v"0.4-"get, :event, Nothing, (ArgsT...)) do ptr, evt_args..., closure
#    stuff
# end
function signal_connect(@nospecialize(cb::Function), w::GObject, sig::AbstractStringLike,
        ::Type{RT}, param_types::Tuple, after::Bool = false, user_data::CT = w) where {CT, RT}
    signal_connect_generic(cb, w, sig, RT, param_types, after, user_data)
end

function signal_connect_generic(@nospecialize(cb::Function), w::GObject, sig::AbstractStringLike,
        ::Type{RT}, param_types::Tuple, after::Bool = false, user_data::CT = w) where {CT, RT}  #TODO: assert that length(param_types) is correct
    callback = cfunction_(cb, RT, tuple(Ptr{GObject}, param_types..., Ref{CT}))
    ref, deref = gc_ref_closure(user_data)
    return ccall((:g_signal_connect_data, libgobject), Culong,
                 (Ptr{GObject}, Ptr{UInt8}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, GEnum),
                 w,
                 bytestring(sig),
                 callback,
                 ref,
                 deref,
                 after * GConnectFlags.AFTER)
end

# id = signal_connect(widget, :event) do obj, evt_args...
#    stuff
# end
function signal_connect(@nospecialize(cb::Function), w::GObject, sig::AbstractStringLike, after::Bool = false)
    _signal_connect(cb, w, sig, after, false, nothing, nothing)
end
function _signal_connect(@nospecialize(cb::Function), w::GObject, sig::AbstractStringLike, after::Bool, gtk_call_conv::Bool, param_types, user_data)
    @assert sizeof_gclosure > 0
    closuref = ccall((:g_closure_new_object, libgobject), Ptr{Nothing}, (Cuint, Ptr{GObject}), sizeof_gclosure::Int + GLib.WORD_SIZE * 2, w)
    closure_env = convert(Ptr{Ptr{Nothing}}, closuref + sizeof_gclosure)
    if gtk_call_conv
        env = Any[param_types, user_data]
        ref_env, deref_env = gc_ref_closure(env)
        unsafe_store!(closure_env, ref_env, 2)
        ccall((:g_closure_add_invalidate_notifier, libgobject), Nothing,
            (Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}), closuref, ref_env, deref_env)
    else
        unsafe_store!(convert(Ptr{Int}, closure_env), 0, 2)
    end
    ref_cb, deref_cb = invoke(gc_ref_closure, Tuple{Function}, cb)
    unsafe_store!(closure_env, ref_cb, 1)
    ccall((:g_closure_add_invalidate_notifier, libgobject), Nothing,
        (Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}), closuref, ref_cb, deref_cb)
    ccall((:g_closure_set_marshal, libgobject), Nothing,
        (Ptr{Nothing}, Ptr{Nothing}), closuref, JuliaClosureMarshal::Ptr{Nothing})
    return ccall((:g_signal_connect_closure, libgobject), Culong,
        (Ptr{GObject}, Ptr{UInt8}, Ptr{Nothing}, Cint), w, bytestring(sig), closuref, after)
end
function GClosureMarshal(closuref::Ptr{Nothing}, return_value::Ptr{GValue}, n_param_values::Cuint,
                         param_values::Ptr{GValue}, invocation_hint::Ptr{Nothing}, marshal_data::Ptr{Nothing})
    @assert sizeof_gclosure > 0
    closure_env = convert(Ptr{Any}, closuref + sizeof_gclosure)
    cb = unsafe_load(closure_env, 1)::Function
    gtk_calling_convention = (0 != unsafe_load(convert(Ptr{Int}, closure_env),  2))
    params = Vector{Any}(undef, n_param_values)
    g_siginterruptible(cb) do
        if gtk_calling_convention
            # compatibility mode, if we must
            param_types, user_data = unsafe_load(closure_env, 2)::Array{Any, 1}
            length(param_types) + 1 == n_param_values || error("GCallback called with the wrong number of parameters")
            for i = 1:n_param_values
                gv = mutable(param_values, i)
                gtyp = unsafe_load(gv).g_type
                # avoid auto-unboxing for some builtin types in gtk_calling_convention mode
                if g_isa(gtyp, g_type(GObject))
                    params[i] = ccall((:g_value_get_object, libgobject), Ptr{GObject}, (Ptr{GValue},), gv)
                elseif g_isa(gtyp, g_type(GBoxed))
                    params[i] = ccall((:g_value_get_boxed, libgobject), Ptr{Nothing}, (Ptr{GValue},), gv)
                elseif g_isa(gtyp, g_type(AbstractString))
                    params[i] = ccall((:g_value_get_string, libgobject), Ptr{Nothing}, (Ptr{GValue},), gv)
                else
                    params[i] = gv[Any]
                end
                if i > 1
                    params[i] = convert(param_types[i - 1], params[i])
                end
            end
            push!(params, user_data)
        else
            for i = 1:n_param_values
                params[i] = mutable(param_values, i)[Any]
            end
        end
        # note: make sure not to leak any of the GValue objects into this task switch, since many of them were alloca'd
        retval = cb(params...) # widget, args...
        if return_value != C_NULL && retval !== nothing
            gtyp = unsafe_load(return_value).g_type
            if gtyp != g_type(Nothing) && gtyp != 0
                try
                    return_value[] = gvalue(retval)
                catch
                    @async begin # make this async to prevent task switches from being present right here
                        blame(cb)
                        println("ERROR: failed to set return value of type $(typeof(retval)); did your callback return an unintentional value?")
                    end
                end
            end
        end
    end
    return nothing
end

function blame(@nospecialize(cb))
    @warn "Executing $cb:"
end

# Signals API for the cb pointer
# Gtk 2
#   https://developer.gnome.org/gtk2/stable/GObject.html#GObject-destroy
#   https://developer.gnome.org/gtk2/stable/GtkWidget.html#GtkWidget-accel-closures-changed
# Gtk 3
#   https://developer.gnome.org/gtk3/stable/GtkWidget.html#GtkWidget-accel-closures-changed

signal_handler_disconnect(w::GObject, handler_id::Culong) =
    ccall((:g_signal_handler_disconnect, libgobject), Nothing, (Ptr{GObject}, Culong), w, handler_id)

signal_handler_block(w::GObject, handler_id::Culong) =
    ccall((:g_signal_handler_block, libgobject), Nothing, (Ptr{GObject}, Culong), w, handler_id)

signal_handler_unblock(w::GObject, handler_id::Culong) =
    ccall((:g_signal_handler_unblock, libgobject), Nothing, (Ptr{GObject}, Culong), w, handler_id)

"""
    tf = signal_handler_is_connected(widget, id)

Return `true`/`false` depending on whether `widget` has a connected signal handler with
the given `id`.
"""
signal_handler_is_connected(w::GObject, handler_id::Culong) =
    ccall((:g_signal_handler_is_connected, libgobject), Cint, (Ptr{GObject}, Culong), w, handler_id) == 1

function signal_emit(w::GObject, sig::AbstractStringLike, ::Type{RT}, args...) where RT
    i = isa(sig, AbstractString) ? something(findfirst("::", sig), 0:-1) : (0:-1)
    if !isempty(i)
        detail = @quark_str sig[last(i) + 1:end]
        sig = sig[1:first(i)-1]
    else
        detail = UInt32(0)
    end
    signal_id = ccall((:g_signal_lookup, libgobject), Cuint, (Ptr{UInt8}, Csize_t), sig, G_OBJECT_CLASS_TYPE(w))
    return_value = RT === Nothing ? C_NULL : gvalue(RT)
    gvals = gvalues(w, args...)
    GC.@preserve gvals return_value ccall((:g_signal_emitv, libgobject), Nothing, (Ptr{GValue}, Cuint, UInt32, Ptr{GValue}), gvals, signal_id, detail, return_value)
    RT === Nothing ? nothing : return_value[RT]
end

g_stack = nothing # need to call g_loop_run from only one stack
const g_yielded = Ref(false) # when true, use the `g_doatomic` queue to run sigatom functions
const g_doatomic = [] # (work, notification) scheduler queue
const g_sigatom_flag = Ref(false) # keep track of Base sigatomic state
function g_sigatom(@nospecialize(f)) # calls f, where f never throws (but this function may throw)
    global g_sigatom_flag, g_stack, g_doatomic
    prev = g_sigatom_flag[]
    stk = g_stack
    ct = current_task()
    if g_yielded[]
        @assert g_stack !== nothing && g_stack != ct && !prev
        push!(g_doatomic, (f, ct))
        return wait()
    end

    if !prev
        sigatomic_begin()
        g_sigatom_flag[] = true
    end
    ret = nothing
    try
        if g_stack === ct
            ret = f()
        else
            @assert g_stack === nothing && !prev
            g_stack = ct
            ret = f()
        end
    catch err
        g_stack = stk
        @assert g_sigatom_flag[]
        if !prev
            g_sigatom_flag[] = false
            sigatomic_end() # may throw SIGINT
        end
        Base.println("FATAL ERROR: Gtk state corrupted by error thrown in a callback:")
        Base.display_error(err, catch_backtrace())
        println()
        rethrow(err)
    end
    g_stack = stk
    @assert g_sigatom_flag[]
    if !prev
        g_sigatom_flag[] = false
        sigatomic_end() # may throw SIGINT
    end
    return ret
end
macro sigatom(f)
    return quote
        g_sigatom() do
            $(esc(f))
        end
    end
end

function g_siginterruptible(f::Base.Callable, @nospecialize(cb)) # calls f (which may throw), but this function never throws
    global g_sigatom_flag, g_stack
    prev = g_sigatom_flag[]
    @assert xor(prev, (current_task() !== g_stack))
    try
        if prev
            # also know that current_task() === g_stack
            g_sigatom_flag[] = false
            sigatomic_end() # may throw SIGINT
        end
        f()
    catch err
        bt = catch_backtrace()
        @async begin # make this async to prevent task switches from being present right here
            blame(cb)
            Base.display_error(err, bt)
            println()
        end
    end
    @assert !g_sigatom_flag[]
    if prev
        sigatomic_begin()
        g_sigatom_flag[] = true
    end
    nothing
end

function g_yield(data)
    global g_yielded, g_doatomic
    while true
        g_yielded[] = true
        g_siginterruptible(yield, yield)
        g_yielded[] = false
        run_delayed_finalizers()

        if isempty(g_doatomic)
            return Int32(true)
        else
            f, t = pop!(g_doatomic)
            ret = nothing
            iserror = false
            try
                ret = f()
            catch err
                iserror = true
                ret = err
            end
            schedule(t, ret, error = iserror)
        end
    end
end

mutable struct _GPollFD
  @static Sys.iswindows() ? fd::Int : fd::Cint
  events::Cushort
  revents::Cushort
  _GPollFD(fd, ev) = new(fd, ev, 0)
end

mutable struct _GSourceFuncs
    prepare::Ptr{Nothing}
    check::Ptr{Nothing}
    dispatch::Ptr{Nothing}
    finalize::Ptr{Nothing}
    closure_callback::Ptr{Nothing}
    closure_marshal::Ptr{Nothing}
end
function new_gsource(source_funcs::_GSourceFuncs)
    sizeof_gsource = GLib.WORD_SIZE
    gsource = C_NULL
    while gsource == C_NULL
        sizeof_gsource += GLib.WORD_SIZE
        gsource = ccall((:g_source_new, GLib.libglib), Ptr{Nothing}, (Ptr{_GSourceFuncs}, Int), Ref(source_funcs), sizeof_gsource)
    end
    gsource
end

expiration = UInt64(0)
_isempty_workqueue() = isempty(Base.Workqueue)
uv_loop_alive(evt) = ccall(:uv_loop_alive, Cint, (Ptr{Nothing},), evt) != 0

function uv_prepare(src::Ptr{Nothing}, timeout::Ptr{Cint})
    global expiration, uv_pollfd
    local tmout_ms::Cint
    evt = Base.eventloop()
    if !_isempty_workqueue()
        tmout_ms = 0
    elseif !uv_loop_alive(evt)
        tmout_ms = -1
    elseif uv_pollfd.revents != 0
        tmout_ms = 0
    else
        ccall(:uv_update_time, Nothing, (Ptr{Nothing},), evt)
        tmout_ms = ccall(:uv_backend_timeout, Cint, (Ptr{Nothing},), evt)
        tmout_min::Cint = (uv_pollfd::_GPollFD).fd == -1 ? 10 : 5000
        if tmout_ms < 0 || tmout_ms > tmout_min
            tmout_ms = tmout_min
        end
    end
    timeout != C_NULL && unsafe_store!(timeout, tmout_ms)
    if tmout_ms < 0
        expiration = typemax(UInt64)
    elseif tmout_ms > 0
        now = ccall((:g_source_get_time, GLib.libglib), UInt64, (Ptr{Nothing},), src)
        expiration = convert(UInt64, now + tmout_ms * 1000)
    else #tmout_ms == 0
        expiration = UInt64(0)
    end
    Int32(tmout_ms == 0)
end
function uv_check(src::Ptr{Nothing})
    global expiration
    ex = expiration::UInt64
    if !_isempty_workqueue()
        return Int32(1)
    elseif !uv_loop_alive(Base.eventloop())
        return Int32(0)
    elseif ex == 0
        return Int32(1)
    elseif uv_pollfd.revents != 0
        return Int32(1)
    else
        now = ccall((:g_source_get_time, GLib.libglib), UInt64, (Ptr{Nothing},), src)
        return Int32(ex <= now)
    end
end
function uv_dispatch(src::Ptr{Nothing}, callback::Ptr{Nothing}, data)
    return ccall(callback, Cint, (UInt,), data)
end

sizeof_gclosure = 0
function __init__gtype__()
    ccall((:g_type_init, libgobject), Nothing, ())
    global jlref_quark = quark"julia_ref"
    global sizeof_gclosure = GLib.WORD_SIZE
    closure = C_NULL
    while closure == C_NULL
        sizeof_gclosure += GLib.WORD_SIZE
        closure = ccall((:g_closure_new_simple, libgobject), Ptr{Nothing}, (Int, Ptr{Nothing}), sizeof_gclosure, C_NULL)
    end
    ccall((:g_closure_sink, libgobject), Nothing, (Ptr{Nothing},), closure)
end

function __init__gmainloop__()
    global uv_sourcefuncs = _GSourceFuncs(
        @cfunction(uv_prepare, Cint, (Ptr{Nothing}, Ptr{Cint})),
        @cfunction(uv_check, Cint, (Ptr{Nothing},)),
        @cfunction(uv_dispatch, Cint, (Ptr{Nothing}, Ptr{Nothing}, Int)),
        C_NULL, C_NULL, C_NULL)
    src = new_gsource(uv_sourcefuncs)
    ccall((:g_source_set_can_recurse, GLib.libglib), Nothing, (Ptr{Nothing}, Cint), src, true)
    ccall((:g_source_set_name, GLib.libglib), Nothing, (Ptr{Nothing}, Ptr{UInt8}), src, "uv loop")
    ccall((:g_source_set_callback, GLib.libglib), Nothing, (Ptr{Nothing}, Ptr{Nothing}, UInt, Ptr{Nothing}),
        src, @cfunction(g_yield, Cint, (UInt,)), 1, C_NULL)

    uv_fd = Sys.iswindows() ? -1 : ccall(:uv_backend_fd, Cint, (Ptr{Nothing},), Base.eventloop())
    global uv_pollfd = _GPollFD(uv_fd, 0x1)
    if (uv_pollfd::_GPollFD).fd != -1
        ccall((:g_source_add_poll, GLib.libglib), Nothing, (Ptr{Nothing}, Ptr{_GPollFD}), src, Ref(uv_pollfd::_GPollFD))
    end

    ccall((:g_source_attach, GLib.libglib), Cuint, (Ptr{Nothing}, Ptr{Nothing}), src, C_NULL)
    ccall((:g_source_unref, GLib.libglib), Nothing, (Ptr{Nothing},), src)
    nothing
end

_g_callback(cb::Function) = Cint(cb())
function g_timeout_add(cb::Function, interval::Integer)
    callback = @cfunction(_g_callback, Cint, (Ref{Function},))
    ref, deref = gc_ref_closure(cb)
    return ccall((:g_timeout_add_full, libglib),Cint,
        (Cint, UInt32, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        0, UInt32(interval), callback, ref, deref)
end

@deprecate g_timeout_add(interval, cb, user_data) g_timeout_add(() -> cb(user_data), interval)

function g_idle_add(cb::Function)
    gtk_eventloop_f[](true)
    callback = @cfunction(_g_callback, Cint, (Ref{Function},))
    ref, deref = gc_ref_closure(cb)
    return ccall((:g_idle_add_full , libglib),Cint,
        (Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        0, callback, ref, deref)
end

@deprecate g_idle_add(cb, user_data)  g_idle_add(() -> cb(user_data))

# Shortcut for g_idle_add
macro idle_add(ex)
    quote
    g_idle_add() do
        $(esc(ex))
        return false
      end
    end
end

const exiting = Ref(false)
function __init__()
    if isdefined(GLib, :__init__bindeps__)
        GLib.__init__bindeps__()
    end
    global JuliaClosureMarshal = @cfunction(GClosureMarshal, Nothing,
        (Ptr{Nothing}, Ptr{GValue}, Cuint, Ptr{GValue}, Ptr{Nothing}, Ptr{Nothing}))
    exiting[] = false
    atexit(() -> (exiting[] = true))
    __init__gtype__()
    __init__gmainloop__()
    nothing
end

max_stars_repo_name: adamwillisMastery/MiniGtk

max_issues_repo_name: adamwillisMastery/MiniGtk

max_forks_repo_name: adamwillisMastery/MiniGtk


text:
using DailyTreasuryYieldCurve, Dates
using Test

testdf = getyieldcurves()

@testset "Nominal" begin

    @test Symbol.(names(testdf)) == DailyTreasuryYieldCurve.COLNAMES[2:end-1]

    @test testdf[1,:date] == Date(1990,1,2)

    @test isapprox(testdf[1,:y30],8.0)
end

@testset "Real" begin
    testdfr = getyieldcurves(;realrates=true)

    @test Symbol.(names(testdfr)) == DailyTreasuryYieldCurve.COLNAMESREAL[2:end]

    @test testdfr[1,:date] == Date(2003,1,2)

    @test isapprox(testdfr[1,:y5],1.752231)
end

@testset "Interpolation" begin
    # testdf = getyieldcurves()
    testri = createRateInterpolator(testdf)

    @test isapprox(testri(45,Date(2020,4,13)),0.23)

    @test isapprox(testri.([30;60],Date(2020,4,13)),[0.17;0.29])
end

@testset "Filters" begin
    # testdf = getyieldcurves()
    smalltestdf = getyieldcurves(;begdt=Date(2019,1,1),enddt=Date(2019,12,31))
    filter!(row->Date(2019,1,1) <= row.date <= Date(2019,12,31),testdf)

    @test testdf == smalltestdf
end


text:
facts("Single objective functions") do

context("Sphere") do
  p = BlackBoxOptim.example_problems["Sphere"]
  sphere = p.funcs[1]

  @fact sphere([0]) => 0

  @fact sphere([1]) => 1

  @fact sphere([1, 2]) => 5

  @fact sphere([1, 2, 3]) => 14

  @fact sphere([-1, 2, -3]) => 14

  @fact_throws sphere([])

  p2 = BlackBoxOptim.as_fixed_dim_problem(p, 3)
  @fact numdims(p2) => 3
  @fact ranges(search_space(p2)) => [(-100.0, 100.0), (-100.0, 100.0), (-100.0, 100.0)]
end

context("Schwefel2.22") do
  p = BlackBoxOptim.example_problems["Schwefel2.22"]
  schwefel2_22 = p.funcs[1]

  @fact schwefel2_22([0]) => 0

  @fact schwefel2_22([1]) => 2

  @fact schwefel2_22([1, 2]) => (1+2)+(1*2)

  @fact schwefel2_22([1, 2, 3]) => (1+2+3)+(1*2*3)

  @fact schwefel2_22([-1, 2, -3]) => (1+2+3)+(1*2*3)

  @fact_throws schwefel2_22([])

  p2 = BlackBoxOptim.as_fixed_dim_problem(p, 4)
  @fact numdims(p2) => 4
  @fact ranges(search_space(p2)) => [(-10.0, 10.0), (-10.0, 10.0), (-10.0, 10.0), (-10.0, 10.0)]
end

context("Schwefel1.2") do
  p = BlackBoxOptim.example_problems["Schwefel1.2"]
  schwefel1_2 = p.funcs[1]

  @fact schwefel1_2([0]) => 0

  @fact schwefel1_2([1]) => 1

  @fact schwefel1_2([1, 2]) => 1+9

  @fact schwefel1_2([1, 2, 3]) => 1+9+36

  @fact schwefel1_2([-1, 2, -3]) => 1+1+4

  @fact schwefel1_2([]) => 0
end

context("Schwefel2.21") do
  p = BlackBoxOptim.example_problems["Schwefel2.21"]
  schwefel2_21 = p.funcs[1]

  @fact schwefel2_21([0]) => 0

  @fact schwefel2_21([1]) => 1

  @fact schwefel2_21([1, 2]) => 2

  @fact schwefel2_21([1, 2, 3]) => 3

  @fact schwefel2_21([-1, 2, -3]) => 3

  @fact_throws schwefel2_21([])
end

context("Rosenbrock") do
  p = BlackBoxOptim.example_problems["Rosenbrock"]
  rosenbrock = p.funcs[1]

  @fact rosenbrock([1, 2]) => 100

  @fact rosenbrock([1, 2, 3]) => 201

  @fact rosenbrock([-1, 2, -3]) => 5005

  @fact_throws rosenbrock([])
end

end


text:
using TerminalLoggers
using Test

using Logging:
    LogLevel, BelowMinLevel, Debug, Info, Warn, Error, AboveMaxLevel,
    shouldlog, handle_message, min_enabled_level, catch_exceptions,
    with_logger

using ProgressLogging: Progress
using UUIDs: UUID

include("TerminalLogger.jl")
include("StickyMessages.jl")


text:
module TSML

using AutoMLPipeline
using AutoMLPipeline.BaseFilters
import AutoMLPipeline.AbsTypes: fit!, transform!

using AutoMLPipeline.AbsTypes
export Machine, Learner, Transformer, Workflow, Computer

using AutoMLPipeline.Utils
export holdout, kfold, score, infer_eltype, nested_dict_to_tuples,
       nested_dict_set!, nested_dict_merge, create_transformer,
       mergedict, getiris,
       skipmean,skipmedian,skipstd,
       aggregatorclskipmissing

export fit!, transform!,fit_transform!

# reexport common functions to Main
include("pkgdeps.jl")
using .PkgDeps

include("baseline.jl")
using .BaselineAlgos
export Baseline,Identity

using AutoMLPipeline.BaseFilters: Imputer, OneHotEncoder, Wrapper
export Imputer,OneHotEncoder,Wrapper

include("valdatefilters.jl")
using .ValDateFilters
export Matrifier,Dateifier,
       DateValizer,DateValgator,DateValNNer,DateValMultiNNer,
       CSVDateValReader, CSVDateValWriter, DateValLinearImputer
       #BzCSVDateValReader
export impute, impute!, chain, chain!, drop, drop!,
       interp, interp!, ImputeError, locf, nocb

include("statifier.jl")
using .Statifiers
export Statifier,tsmlfullstat

include("mlbase.jl")
using .MLBaseWrapper
export Standardize,standardize, standardize!,
       transform, estimate, transform,StandardScaler

using AutoMLPipeline.DecisionTreeLearners
export PrunedTree,RandomForest,Adaboost

include("normalizer.jl")
using .Normalizers
export Normalizer

include("monotonicer.jl")
using .Monotonicers
export Monotonicer,ismonotonic,dailyflips

include("cliwrapper.jl")
using .CLIWrappers
export tsmlrun

include("tsclassifier.jl")
using .TSClassifiers
export TSClassifier, getstats

include("outliernicer.jl")
using .Outliernicers
export Outliernicer

include("plotter.jl")
using .Plotters
export Plotter

include("timescaledb.jl")
using .TimescaleDBs
export TimescaleDB

using AutoMLPipeline.EnsembleMethods
export VoteEnsemble, StackEnsemble, BestLearner

using AutoMLPipeline.FeatureSelectors
export FeatureSelector, CatFeatureSelector, NumFeatureSelector, CatNumDiscriminator

using AutoMLPipeline.Pipelines
export @pipeline @pipelinex
export Pipeline, ComboPipeline

include("schema.jl")
using .Schemalizers
export Schemalizer, ML, table

using AutoMLPipeline.CrossValidators
export crossvalidate

include("argparse.jl")
using .ArgumentParsers
export tsmlmain

include("demo.jl")
using .TSMLDemo
export tsml_demo
end # module


text:
using DataFrames, Downloads

function list_notebooks()
    fil=Downloads.download("https://raw.githubusercontent.com/JuliaClimate/Notebooks/master/page/index.md")

    lines=readlines(fil);
    ii=findall( occursin.("[notebook url]",lines) )
    lines=lines[ii]

    sources=["gaelforget" "JuliaClimate" "JuliaOcean" "MarineEcosystemsJuliaCon2021.jl"]
    notebooks=DataFrame("folder" => [], "file" => [], "url" => [])

    for ii in lines
        tmp1=split(ii,"(")
        ii=findall( occursin.("[notebook url]",tmp1) )[1]
        tmp1=tmp1[ii+1]
        tmp1=split(tmp1,")")[1]

        tmp2=split(tmp1,"/")
        test=sum(tmp2[4].==sources)
        if test!==1
            println("Skipping : ")
            println(tmp1)
        end

        push!(notebooks, (tmp2[5], tmp2[end], string(tmp1) ))
    end

    notebooks
end

function download_notebooks(path,notebooks)
    !isdir(path) ? mkdir(path) : nothing
    for i in 1:size(notebooks,1)
        tmp1=joinpath(path,notebooks[i,:folder])
        !isdir(tmp1) ? mkdir(tmp1) : nothing
        tmp2=joinpath(tmp1,notebooks[i,:file])
        try
            Downloads.download(notebooks[i,:url],tmp2)
        catch e
            println("Skipping : ")
            println(tmp2)
        end
    end
end

max_issues_repo_name: gaelforget/MeshArrayNotebooks

max_forks_repo_name: JuliaClimate/Notebooks


text:
#=
Styling file for the figures of the velocity correlation paper.
=#
using PyPlot
PyPlot.rc("lines", lw = 2.6)

PyPlot.rc("font", size = 26) # set default fontsize
PyPlot.rc("axes", labelsize = 32)
PyPlot.rc("legend", fontsize = 26)

PyPlot.rc("xtick.major", size = 10)
PyPlot.rc("xtick.minor", size = 6)

# PyPlot.rc("font", family = "Times New Roman") # Serif main font
PyPlot.rc("font", family = "DejaVu Sans") # sans main font
# PyPlot.rc("mathtext", rm = "sanserif", fontset="dejavusans") # sans math font
PyPlot.rc("mathtext", rm = "serif", fontset="cm") # serif math font

const figx = 12
# default figsize (full span, 2-3 rows)
PyPlot.rc("figure", figsize = (2figx, figx))
PyPlot.rc("savefig", dpi = 300, transparent = true, format = "pdf")

mutable struct CyclicContainer
    c::Vector
    n::Int
end
CyclicContainer(c) = CyclicContainer(c, 0)

Base.length(c::CyclicContainer) = length(c.c)
Base.iterate(c::CyclicContainer, state=1) = Base.iterate(c.c, state)
Base.getindex(c::CyclicContainer, i) = c.c[(i-1)%length(c.c) + 1]
function Base.getindex(c::CyclicContainer)
    c.n += 1
    c[c.n]
end
COLORS = CyclicContainer(
["#233B43", "#499cbf", "#E84646",
 "#168e6c","#C29365", "#985CC9",
 "#822e01", "#7c6d8c"])
MARKERS = CyclicContainer(["o", "s", "^", "p", "P", "D", "X"])
LINES = CyclicContainer([ ":", "-", "--", "-."])
# Also set default color cycle
PyPlot.rc("axes", prop_cycle = matplotlib.cycler(color=COLORS.c))
PyPlot.rc("lines", markersize = 10)

_bbox(α) = Dict(:boxstyle => "round,pad=0.3", :facecolor=>"white", :alpha => α)

bbox = Dict(:boxstyle => "round,pad=0.3", :facecolor=>"white", :alpha => 1.0)

function add_identifiers!(fig = gcf(); loc = (0.985, 0.975))
    bbox = Dict(:boxstyle => "round,pad=0.3", :facecolor=>"white", :alpha => 1.0)
    for (i, ax) in enumerate(fig.get_axes())
        l = collect('a':'z')[i]
        try
            ax.text(loc..., "$(l)",
            transform = ax.transAxes, bbox = bbox, zorder = 99)
        catch err
            ax.text(loc..., 1, "$(l)",
            transform = ax.transAxes, bbox = bbox, zorder = 99)
        end
    end
end
function coolhist!(ax, data, bins, color, label = "", alpha = 0.25)
    h, b, = ax.hist(data, bins, density = true, color = color,
    alpha = alpha)

    b = 0.5(b[1:end-1] .+ b[2:end])
    ax.plot(b, h, color = color, lw = 1.0, label = label)
end

max_stars_repo_name: Datseris/VelocityMTD

max_issues_repo_name: Datseris/VelocityMTD

max_forks_repo_name: Datseris/VelocityMTD


text:
import PyPlot as Plt
using SparseArrays

function sparse_counts(vec1::Vector{<:Integer}, vec2::Vector{<:Integer})
    counts = Dict{Tuple{Int, Int}, Int}()
    for tp in zip(vec1, vec2)
        counts[tp] = get(counts, tp, 0) + 1
    end
    i1 = getindex.(keys(counts), 1)
    i2 = getindex.(keys(counts), 2)
    v = collect(values(counts))
    return sparse(i1, i2, v)
end

function vec_triu(M::AbstractMatrix{T}) where T
    m, n = size(M)
    m == n || throw(error("not square"))
    l = n*(n-1) ÷ 2
    v = zeros(T, l)
    k = 0
    @inbounds for i in 1:n
        for j in 1:(i-1)
            v[k + j] = M[j, i]
        end
        k += i-1
    end
    v
end

function method_palette(subs::Union{Vector{String}, Nothing}=nothing)
    d = Dict(
        "Baysor" => "#E69F00",
        "Watershed" => "#3093ca",
        "pciSeq" => "#1e9e00",
        "poly-A" => "#F0E442",
        "IF" => "#CC79A7",
        "Baysor, DAPI prior" => "#9566FF",
        "Baysor, IF prior" => "#6F4D0B"
    )

    (subs !== nothing) || return d
    return Dict(k => d[k] for k in subs)
end

function set_pyplot_defaults!()
    Plt.rc("axes", axisbelow=true, labelsize=16, grid=true);
    Plt.rc("axes.spines", right=false, top=false, left=true, bottom=true);
    Plt.rc("font", size=12)
    Plt.rc("grid", alpha=0.25)
    Plt.rc("legend", frameon=false, borderpad=0)
end

function upscale(image::T where T <: AbstractMatrix{TR}, ratio::Int64) where TR <: Real
    res_img = zeros(TR, size(image) .* ratio)
    ct = 1
    for c in 1:size(image, 2)
        for m1 in 1:ratio
            rt = 1
            for r in 1:size(image, 1)
                for m2 in 1:ratio
                    res_img[rt, ct] = image[r, c]
                    rt += 1
                end
            end
            ct += 1
        end
    end

    return res_img
end

max_stars_repo_name: hiraksarkar/BaysorAnalysis

max_issues_repo_name: hiraksarkar/BaysorAnalysis

max_forks_repo_name: hiraksarkar/BaysorAnalysis


text:
module Interface
using Luna
import Luna.PhysData: wlfreq
import Luna: Grid, Modes, Output, Fields
import Logging: @info, @debug

module Pulses

import Luna: Fields, Output

export AbstractPulse, CustomPulse, GaussPulse, SechPulse, DataPulse, LunaPulse

abstract type AbstractPulse end

struct CustomPulse{fT<:Fields.TimeField} <: AbstractPulse
    mode::Symbol
    polarisation
    field::fT
end

"""
    CustomPulse(;λ0, energy=nothing, power=nothing, ϕ=Float64[],
                mode=:lowest, polarisation=:linear, propagator=nothing)

A custom pulse defined by a function for use with `prop_capillary`, with either energy or
peak power specified.

# Keyword arguments
- `λ0::Number`: the central wavelength
- `Itshape::function`: a function `I(t)`` which defines the intensity/power envelope of the
                       pulse as a function of time `t`. Note that the normalisation of this
                       envelope is irrelevant as it will be re-scaled by `energy` or `power`.
- `energy::Number`: the pulse energy.
- `power::Number`: the pulse peak power (**after** applying any spectral phases).
- `ϕ::Vector{Number}`: spectral phases (CEP, group delay, GDD, TOD, ...).
- `mode::Symbol`: Mode in which this input should be coupled. Can be `:lowest` for the
                  lowest-order mode in the simulation, or a mode designation
                  (e.g. `:HE11`, `:HE12`, `:TM01`, etc.). Defaults to `:lowest`.
- `polarisation`: Can be `:linear`, `:circular`, or an ellipticity number -1 ≤ ε ≤ 1,
                  where ε=-1 corresponds to left-hand circular, ε=1 to right-hand circular,
                  and ε=0 to linear polarisation.
- `propagator`: A function `propagator!(Eω, grid)` which **mutates** its first argument to
                apply an arbitrary propagation to the pulse before the simulation starts.
"""
function CustomPulse(;mode=:lowest, polarisation=:linear, propagator=nothing, kwargs...)
    CustomPulse(mode, polarisation,
                Fields.PropagatedField(propagator, Fields.PulseField(;kwargs...)))
end

struct GaussPulse{fT<:Fields.TimeField} <: AbstractPulse
    mode::Symbol
    polarisation
    field::fT
end

"""
    GaussPulse(;λ0, τfwhm, energy=nothing, power=nothing, ϕ=Float64[], m=1,
               mode=:lowest, polarisation=:linear, propagator=nothing)

A (super)Gaussian pulse for use with `prop_capillary`, with either energy or peak power
specified.

# Keyword arguments
- `λ0::Number`: the central wavelength.
- `τfwhm::Number`: the pulse duration (power/intensity FWHM).
- `energy::Number`: the pulse energy.
- `power::Number`: the pulse peak power (**after** applying any spectral phases).
- `ϕ::Vector{Number}`: spectral phases (CEP, group delay, GDD, TOD, ...).
- `m::Int`: super-Gaussian parameter (the power in the Gaussian exponent is 2m).
            Defaults to 1.
- `mode::Symbol`: Mode in which this input should be coupled. Can be `:lowest` for the
                  lowest-order mode in the simulation, or a mode designation
                  (e.g. `:HE11`, `:HE12`, `:TM01`, etc.). Defaults to `:lowest`.
- `polarisation`: Can be `:linear`, `:circular`, or an ellipticity number -1 ≤ ε ≤ 1,
                  where ε=-1 corresponds to left-hand circular, ε=1 to right-hand circular,
                  and ε=0 to linear polarisation.
- `propagator`: A function `propagator!(Eω, grid)` which **mutates** its first argument to
                apply an arbitrary propagation to the pulse before the simulation starts.
"""
function GaussPulse(;mode=:lowest, polarisation=:linear, propagator=nothing, kwargs...)
    GaussPulse(mode, polarisation,
               Fields.PropagatedField(propagator, Fields.GaussField(;kwargs...)))
end

struct SechPulse{fT<:Fields.TimeField} <: AbstractPulse
    mode::Symbol
    polarisation
    field::fT
end

"""
    SechPulse(;λ0, τfwhm=nothing, τw=nothing, energy=nothing, power=nothing, ϕ=Float64[],
               mode=:lowest, polarisation=:linear, propagator=nothing)

A sech²(τ/τw) pulse for use with `prop_capillary`, with either `energy` or peak `power`
specified, and duration given either as `τfwhm` or `τw`.

# Keyword arguments
- `λ0::Number`: the central wavelength.
- `τfwhm::Number`: the pulse duration (power/intensity FWHM).
- `τw::Number`: "natural" pulse duration of a sech²(τ/τw) pulse.
- `energy::Number`: the pulse energy.
- `power::Number`: the pulse peak power (**after** applying any spectral phases).
- `ϕ::Vector{Number}`: spectral phases (CEP, group delay, GDD, TOD, ...)
- `mode::Symbol`: Mode in which this input should be coupled. Can be `:lowest` for the
                  lowest-order mode in the simulation, or a mode designation
                  (e.g. `:HE11`, `:HE12`, `:TM01`, etc.). Defaults to `:lowest`.
- `polarisation`: Can be `:linear`, `:circular`, or an ellipticity number -1 ≤ ε ≤ 1,
                  where ε=-1 corresponds to left-hand circular, ε=1 to right-hand circular,
                  and ε=0 to linear polarisation.
- `propagator`: A function `propagator!(Eω, grid)` which **mutates** its first argument to
                apply an arbitrary propagation to the pulse before the simulation starts.
"""
function SechPulse(;mode=:lowest, polarisation=:linear, propagator=nothing, kwargs...)
    SechPulse(mode, polarisation,
              Fields.PropagatedField(propagator, Fields.SechField(;kwargs...)))
end

struct DataPulse{fT<:Fields.TimeField} <: AbstractPulse
    mode::Symbol
    polarisation
    field::fT
end

#TODO add peak power to DataPulses
"""
    DataPulse(ω, Iω, ϕω; energy, λ0=NaN, mode=:lowest, polarisation=:linear, propagator=nothing)
    DataPulse(ω, Eω; energy, λ0=NaN, mode=:lowest, polarisation=:linear, propagator=nothing)
    DataPulse(fpath; energy, λ0=NaN, mode=:lowest, polarisation=:linear, propagator=nothing)

A custom pulse defined by tabulated data to be used with `prop_capillary`.

# Data input options
- `ω, Iω, ϕω`: arrays of angular frequency `ω` (units rad/s), spectral energy density `Iω`
               and spectral phase `ϕω`. `ϕω` should be unwrapped.
- `ω, Eω`: arrays of angular frequency `ω` (units rad/s) and the complex frequency-domain
           field `Eω`.
- `fpath`: a string containing the path to a file which contains 3 columns:
    Column 1: frequency (units of Hertz)
    Column 2: spectral energy density
    Column 3: spectral phase (unwrapped)

# Keyword arguments
- `energy::Number`: the pulse energy
- `λ0::Number`: the central wavelength (optional; defaults to the centre of mass of the
                given spectral energy density).
- `ϕ::Vector{Number}`: spectral phases (CEP, group delay, GDD, TOD, ...) to be applied to the
                       pulse (in addition to any phase already present in the data).
- `mode::Symbol`: Mode in which this input should be coupled. Can be `:lowest` for the
                  lowest-order mode in the simulation, or a mode designation
                  (e.g. `:HE11`, `:HE12`, `:TM01`, etc.). Defaults to `:lowest`.
- `polarisation`: Can be `:linear`, `:circular`, or an ellipticity number -1 ≤ ε ≤ 1,
                  where ε=-1 corresponds to left-hand circular, ε=1 to right-hand circular,
                  and ε=0 to linear polarisation.
- `propagator`: A function `propagator!(Eω, grid)` which **mutates** its first argument to
                apply an arbitrary propagation to the pulse before the simulation starts.
"""
function DataPulse(ω::AbstractVector, Iω, ϕω;
                   mode=:lowest, polarisation=:linear, propagator=nothing, kwargs...)
    DataPulse(mode, polarisation,
              Fields.PropagatedField(propagator, Fields.DataField(ω, Iω, ϕω; kwargs...)))
end

function DataPulse(ω, Eω;
                   mode=:lowest, polarisation=:linear, propagator=nothing, kwargs...)
    DataPulse(mode, polarisation,
              Fields.PropagatedField(propagator, Fields.DataField(ω, Eω; kwargs...)))
end

function DataPulse(fpath;
                   mode=:lowest, polarisation=:linear, propagator=nothing, kwargs...)
    DataPulse(mode, polarisation,
              Fields.PropagatedField(propagator, Fields.DataField(fpath; kwargs...)))
end

# Select fundamental mode from multi-mode sim or take just the single mode
modeslice(Eω::Array{ComplexF64, 2}) = Eω[:, end]
modeslice(Eω::Array{ComplexF64, 3}) = Eω[:, 1, end]

"""
    LunaPulse(output; energy, λ0=NaN, mode=:lowest, polarisation=:linear, propagator=nothing)

A pulse defined to be used with `prop_capillary` which comes from a previous `Luna`
propagation simulation.

For multi-mode simulations, only the lowest-order modes is transferred.

# Arguments
- `output::AbstractOutput`: output from a previous `Luna` simulation.

# Keyword arguments
- `energy::Number`: the pulse energy
- `λ0::Number`: the central wavelength (optional; defaults to the centre of mass of the
                given spectral energy density).
- `ϕ::Vector{Number}`: spectral phases (CEP, group delay, GDD, TOD, ...) to be applied to the
                       pulse (in addition to any phase already present in the data).
- `mode::Symbol`: Mode in which this input should be coupled. Can be `:lowest` for the
                  lowest-order mode in the simulation, or a mode designation
                  (e.g. `:HE11`, `:HE12`, `:TM01`, etc.). Defaults to `:lowest`.
- `polarisation`: Can be `:linear`, `:circular`, or an ellipticity number -1 ≤ ε ≤ 1,
                  where ε=-1 corresponds to left-hand circular, ε=1 to right-hand circular,
                  and ε=0 to linear polarisation.
- `propagator`: A function `propagator!(Eω, grid)` which **mutates** its first argument to
                apply an arbitrary propagation to the pulse before the simulation starts.
"""
function LunaPulse(o::Output.AbstractOutput; kwargs...)
    ω = o["grid"]["ω"]
    t = o["grid"]["t"]
    τ = length(t) * (t[2] - t[1])/2 # middle of old time window
    Eωm1 = modeslice(o["Eω"]) # either mode-averaged field or first mode
    DataPulse(ω, Eωm1 .* exp.(1im .* ω .* τ); kwargs...)
end

struct GaussBeamPulse{pT} <: AbstractPulse
    waist::Float64
    timepulse::pT
    polarisation
end

"""
    GaussBeamPulse(waist, timepulse)

A pulse whose shape in time is defined by the `timepulse::AbstractPulse`, and whose modal content is calculated by considering the overlap of an ideal Gaussian laser beam with 1/e² radius `waist` with the modes of the waveguide.
"""
function GaussBeamPulse(waist, timepulse)
    GaussBeamPulse(waist, timepulse, timepulse.polarisation)
end

end

"""
    prop_capillary(radius, flength, gas, pressure; λ0, λlims, trange, kwargs...)

Simulate pulse propagation in a hollow fibre using the capillary model.

# Mandatory arguments
- `radius`: Core radius of the fibre. Can be a `Number` for constant radius, or a function
    `a(z)` which returns the `z`-dependent radius.
- `flength::Number`: Length of the fibre.
- `gas::Symbol`: Filling gas species.
- `pressure`: Gas pressure. Can be a `Number` for constant pressure, a 2-`Tuple` of `Number`s
    for a simple pressure gradient, or a `Tuple` of `(Z, P)` where `Z` and `P`
    contain `z` positions and the pressures at those positions.
- `λ0`: (keyword argument) the reference wavelength for the simulation. For simple
    single-pulse inputs, this is also the central wavelength of the input pulse.
- `λlims::Tuple{<:Number, <:Number}`: The wavelength limits for the simulation grid.
- `trange::Number`: The total width of the time grid. To make the number of samples a
    power of 2, the actual grid used is usually bigger.

# Grid options
- `envelope::Bool`: Whether to use envelope fields for the simulation. Defaults to `false`.
    By default, envelope simulations ignore third-harmonic generation.
    Plasma has not yet been implemented for envelope fields.
- `δt::Number`: Time step on the fine grid used for the nonlinear interaction. By default,
    this is determined by the wavelength grid. If `δt` is given **and smaller** than the
    required value, it is used instead.

# Input pulse options
A single pulse in the lowest-order mode can be specified by the keyword arguments below.
More complex inputs can be defined by a single `AbstractPulse` or a `Vector{AbstractPulse}`.
In this case, all keyword arguments except for `λ0` are ignored.

- `λ0`: Central wavelength
- `τfwhm`: The pulse duration as defined by the full width at half maximum.
- `τw`: The "natural" pulse duration. Only available if pulseshape is `sech`.
- `ϕ`: Spectral phases to be applied to the transform-limited pulse. Elements are
    the usual polynomial phases ϕ₀ (CEP), ϕ₁ (group delay), ϕ₂ (GDD), ϕ₃ (TOD), etc.
- `energy`: Pulse energy.
- `power`: Peak power **after any spectral phases are added**.
- `pulseshape`: Shape of the transform-limited pulse. Can be `:gauss` for a Gaussian pulse
    or `:sech` for a sech² pulse.
- `polarisation`: Polarisation of the input pulse. Can be `:linear` (default), `:circular`,
    or an ellipticity number -1 ≤ ε ≤ 1, where ε=-1 corresponds to left-hand circular,
    ε=1 to right-hand circular, and ε=0 to linear polarisation. The major axis for
    elliptical polarisation is always the y-axis.
- `propagator`: A function `propagator!(Eω, grid)` which **mutates** its first argument to
                apply an arbitrary propagation to the pulse before the simulation starts.
- `shotnoise`:  If `true` (default), one-photon-per-mode quantum noise is included.

# Modes options
- `modes`: Defines which modes are included in the propagation. Can be any of:
    - a single mode signifier (default: :HE11), which leads to mode-averaged propagation
        (as long as all inputs are linearly polarised).
    - a list of mode signifiers, which leads to multi-mode propagation in those modes.
    - a `Number` `N` of modes, which simply creates the first `N` `HE` modes.
    Note that when elliptical or circular polarisation is included, each mode is present
    twice in the output, once for `x` and once for `y` polarisation.
- `model::Symbol`: Can be `:full`, which includes the full complex refractive index of the cladding
    in the effective index of the mode, or `:reduced`, which uses the simpler model more
    commonly seen in the literature. See `Luna.Capillary` for more details.
    Defaults to `:full`.
- `loss::Bool`: Whether to include propagation loss. Defaults to `true`.

# Nonlinear interaction options
- `kerr`: Whether to include the Kerr effect. Defaults to `true`.
- `raman`: Whether to include the Raman effect. Defaults to `false`.
- `plasma`: Can be one of
    - `:ADK` -- include plasma using the ADK ionisation rate.
    - `:PPT` -- include plasma using the PPT ionisation rate.
    - `true` (default) -- same as `:PPT`.
    - `false` -- ignore plasma.
    Note that plasma is only available for full-field simulations.
- `thg::Bool`: Whether to include third-harmonic generation. Defaults to `true` for
    full-field simulations and to `false` for envelope simulations.
If `raman` is `true`, then the following options apply:
    - `rotation::Bool = true`: whether to include the rotational Raman contribution
    - `vibration::Bool = true`: whether to include the vibrational Raman contribution

# Output options
- `saveN::Integer`: Number of points along z at which to save the field.
- `filepath`: If `nothing` (default), create a `MemoryOutput` to store the simulation results
    only in the working memory. If not `nothing`, should be a file path as a `String`,
    and the results are saved in a file at this location. If `scan` is passed, `filepath`
    determines the output **directory** for the scan instead.
- `scan`: A `Scan` instance defining a parameter scan. If `scan` is given`, a
    `Output.ScanHDF5Output` is used to automatically name and populate output files of
    the scan. `scanidx` must also be given.
- `scanidx`: Current scan index within a scan being run. Only used when `scan` is passed.
- `filename`: Can be used to to overwrite the scan name when running a parameter scan.
    The running `scanidx` will be appended to this filename. Ignored if no `scan` is given.
- `status_period::Number`: Interval (in seconds) between printed status updates.
"""
function prop_capillary(args...; status_period=5, kwargs...)
    Eω, grid, linop, transform, FT, output = prop_capillary_args(args...; kwargs...)
    Luna.run(Eω, grid, linop, transform, FT, output; status_period)
    output
end

"""
    prop_capillary_args(radius, flength, gas, pressure; λ0, λlims, trange, kwargs...)

Prepare to simulate pulse propagation in a hollow fibre using the capillary model. This
function takes the same arguments as `prop_capillary` but instead or running the
simulation and returning the output, it returns the required arguments for `Luna.run`,
which is useful for repeated simulations in an indentical fibre with different initial
conditions.
"""
function prop_capillary_args(radius, flength, gas, pressure;
                        λlims, trange, envelope=false, thg=nothing, δt=1,
                        λ0, τfwhm=nothing, τw=nothing, ϕ=Float64[],
                        power=nothing, energy=nothing,
                        pulseshape=:gauss, polarisation=:linear, propagator=nothing,
                        pulses=nothing,
                        shotnoise=true,
                        modes=:HE11, model=:full, loss=true,
                        raman=false, kerr=true, plasma=nothing,
                        rotation=true, vibration=true,
                        saveN=201, filepath=nothing,
                        scan=nothing, scanidx=nothing, filename=nothing)

    pol = needpol(polarisation, pulses) || needpol_modes(modes)
    @info "X+Y polarisation "* (pol ? "required." : "not required.")
    plasma = isnothing(plasma) ? !envelope : plasma
    thg = isnothing(thg) ? !envelope : thg

    gas = (gas == :He) ? :HeJ : gas

    grid = makegrid(flength, λ0, λlims, trange, envelope, thg, δt)
    mode_s = makemode_s(modes, flength, radius, gas, pressure, model, loss, pol)
    check_orth(mode_s)
    density = makedensity(flength, gas, pressure)
    resp = makeresponse(grid, gas, raman, kerr, plasma, thg, pol, rotation, vibration)
    inputs = makeinputs(mode_s, λ0, pulses, τfwhm, τw, ϕ,
                        power, energy, pulseshape, polarisation, propagator)
    inputs = shotnoise_maybe(inputs, mode_s, shotnoise)
    linop, Eω, transform, FT = setup(grid, mode_s, density, resp, inputs, pol,
                                     const_linop(radius, pressure))
    stats = Stats.default(grid, Eω, mode_s, linop, transform; gas=gas)
    output = makeoutput(grid, saveN, stats, filepath, scan, scanidx, filename)

    saveargs(output; radius, flength, gas, pressure, λlims, trange, envelope, thg, δt,
        λ0, τfwhm, τw, ϕ, power, energy, pulseshape, polarisation, propagator, pulses,
        shotnoise, modes, model, loss, raman, kerr, plasma, saveN, filepath, filename)

    return Eω, grid, linop, transform, FT, output
end

check_orth(mode::Modes.AbstractMode) = nothing
function check_orth(modes)
    if length(modes) > 1
        if !Modes.orthonormal(modes)
            ms = join(modes, "\n")
            error("The selected modes do not form an orthonormal set:\n$ms")
        end
    end
end

function saveargs(output; kwargs...)
    d = Dict{String, String}()
    for (k, v) in kwargs
        d[string(k)] = string(v)
    end
    output(d; group="prop_capillary_args")
end

function needpol(pol)
    if pol == :linear
        return false
    elseif pol == :circular
        return true
    else
        error("Polarisation must be :linear, :circular, or an ellipticity, not $pol")
    end
end
needpol(pol::Number) = true
needpol(pulse::Pulses.AbstractPulse) = needpol(pulse.polarisation)

needpol(pol, pulses::Nothing) = needpol(pol)
needpol(pol, pulse::Pulses.AbstractPulse) = needpol(pulse)
needpol(pol, pulses) = any(needpol, pulses)

needpol_modes(mode::Symbol) = false # mode average
needpol_modes(modes::Number) = false # only HE1m modes

function needpol_modes(modes::NTuple{N, Symbol}) where N
    any(modes) do mode
        md = parse_mode(mode)
        md[:kind] ≠ :HE || md[:n] > 1
    end
end

const_linop(radius::Number, pressure::Number) = Val(true)
const_linop(radius, pressure) = Val(false)

function makegrid(flength, λ0, λlims, trange, envelope, thg, δt)
    if envelope
        isnothing(thg) && (thg = false)
        Grid.EnvGrid(flength, λ0, λlims, trange; δt, thg)
    else
        Grid.RealGrid(flength, λ0, λlims, trange, δt)
    end
end

makegrid(flength, λ0::Tuple, args...) = makegrid(flength, λ0[1], args...)

function parse_mode(mode)
    ms = String(mode)
    Dict(:kind => Symbol(ms[1:2]), :n => parse(Int, ms[3]), :m => parse(Int, ms[4]))
end

function makemodes_pol(pol, args...; kwargs...)
    if pol
        if kwargs[:kind] == :HE && kwargs[:n] == 1
            return [Capillary.MarcatiliMode(args...; ϕ=0.0, kwargs...),
                    Capillary.MarcatiliMode(args...; ϕ=π/2, kwargs...)]
        else
            return [Capillary.MarcatiliMode(args...; ϕ=0.0, kwargs...)]
        end
    else
        Capillary.MarcatiliMode(args...; kwargs...)
    end
end

function makemode_s(mode::Symbol, flength, radius, gas, pressure::Number, model, loss, pol)
    makemodes_pol(pol, radius, gas, pressure; model, loss, parse_mode(mode)...)
end

function makemode_s(mode::Symbol, flength, radius, gas, pressure::Tuple{<:Number, <:Number},
                    model, loss, pol)
    coren, _ = Capillary.gradient(gas, flength, pressure...)
    makemodes_pol(pol, radius, coren; model, loss, parse_mode(mode)...)
end

function makemode_s(mode::Symbol, flength, radius, gas, pressure, model, loss, pol)
    Z, P = pressure
    coren, _ = Capillary.gradient(gas, Z, P)
    makemodes_pol(pol, radius, coren; model, loss, parse_mode(mode)...)
end

function makemode_s(modes::Int, args...)
    _flatten([makemode_s(Symbol("HE1$n"), args...) for n=1:modes])
end

function makemode_s(modes::NTuple{N, Symbol}, args...) where N
    _flatten([makemode_s(m, args...) for m in modes])
end

# Iterators.flatten recursively flattens arrays of arrays, but can't handle scalars
_flatten(modes::Vector{<:AbstractArray}) = collect(Iterators.flatten(modes))
_flatten(mode) = mode

function makedensity(flength, gas, pressure::Number)
    ρ0 = PhysData.density(gas, pressure)
    z -> ρ0
end

function makedensity(flength, gas, pressure::Tuple{<:Number, <:Number})
    _, density = Capillary.gradient(gas, flength, pressure...)
    density
end

function makedensity(flength, gas, pressure)
    _, density = Capillary.gradient(gas, pressure...)
    density
end

function makeresponse(grid::Grid.RealGrid, gas, raman, kerr, plasma, thg, pol,
                      rotation, vibration)
    out = Any[]
    if kerr
        if thg
            push!(out, Nonlinear.Kerr_field(PhysData.γ3_gas(gas)))
        else
            push!(out, Nonlinear.Kerr_field_nothg(PhysData.γ3_gas(gas), length(grid.to)))
        end
    end
    makeplasma!(out, grid, gas, plasma, pol)
    if raman
        rr = Raman.raman_response(grid.to, gas, rotation=rotation, vibration=vibration)
        if thg
            push!(out, Nonlinear.RamanPolarField(grid.to, rr))
        else
            push!(out, Nonlinear.RamanPolarField(grid.to, rr, thg=false))
        end
    end
    Tuple(out)
end

function makeplasma!(out, grid, gas, plasma::Bool, pol)
    # simple true/false => default to PPT
    plasma && makeplasma!(out, grid, gas, :PPT, pol)
end

function makeplasma!(out, grid, gas, plasma::Symbol, pol)
    ionpot = PhysData.ionisation_potential(gas)
    if plasma == :ADK
        ionrate = Ionisation.ionrate_fun!_ADK(gas)
    elseif plasma == :PPT
        ionrate = Ionisation.ionrate_fun!_PPTcached(gas, grid.referenceλ)
    else
        throw(DomainError(plasma, "Unknown ionisation rate $plasma."))
    end
    Et = pol ? Array{Float64}(undef, length(grid.to), 2) : grid.to
    push!(out, Nonlinear.PlasmaCumtrapz(grid.to, Et, ionrate, ionpot))
end

function makeresponse(grid::Grid.EnvGrid, gas, raman, kerr, plasma, thg, pol,
                      rotation, vibration)
    plasma && error("Plasma response for envelope fields has not been implemented yet.")
    isnothing(thg) && (thg = false)
    out = Any[]
    if kerr
        if thg
            ω0 = wlfreq(grid.referenceλ)
            r = Nonlinear.Kerr_env_thg(PhysData.γ3_gas(gas), ω0, grid.to)
            push!(out, r)
        else
            push!(out, Nonlinear.Kerr_env(PhysData.γ3_gas(gas)))
        end
    end
    if raman
        rr = Raman.raman_response(grid.to, gas, rotation=rotation, vibration=vibration)
        push!(out, Nonlinear.RamanPolarEnv(grid.to, rr))
    end
    Tuple(out)
end

getAeff(mode::Modes.AbstractMode) = Modes.Aeff(mode)
getAeff(modes) = Modes.Aeff(modes[1])

function makeinputs(mode_s, λ0, pulses::Nothing, τfwhm, τw, ϕ, power, energy,
                    pulseshape, polarisation, propagator)
    if pulseshape == :gauss
        return makeinputs(mode_s, λ0, Pulses.GaussPulse(;λ0, τfwhm, power=power, energy=energy,
                          polarisation, ϕ, propagator))
    elseif pulseshape == :sech
        return makeinputs(mode_s, λ0, Pulses.SechPulse(;λ0, τfwhm, τw, power=power, energy=energy,
                          polarisation, ϕ, propagator))
    else
        error("Valid pulse shapes are :gauss and :sech")
    end
end

function makeinputs(mode_s, λ0, pulses, args...)
    makeinputs(mode_s, λ0, pulses)
end

function findmode(mode_s, pulse)
    if pulse.mode == :lowest
        if pulse.polarisation == :linear
            return [1]
        else
            return [1, 2]
        end
    else
        md = parse_mode(pulse.mode)
        return _findmode(mode_s, md)
    end
end

function _findmode(mode_s::AbstractArray, md)
    return findall(mode_s) do m
        (m.kind == md[:kind]) && (m.n == md[:n]) && (m.m == md[:m])
    end
end

function makeinputs(mode_s, λ0, pulse::Pulses.GaussBeamPulse)
    k = 2π/λ0
    gauss = Fields.normalised_gauss_beam(k, pulse.waist)
    facs = [abs2(Modes.overlap(mi, gauss)) for mi in mode_s]
    fields = Any[]
    if pulse.polarisation == :linear
        for (modeidx, fac) in enumerate(facs)
            sf = scalefield(pulse.timepulse.field, fac)
            push!(fields, (mode=modeidx, fields=(sf,)))
        end
    else
        fy, fx = ellfields(pulse.timepulse)
        for (idx, fac) in enumerate(facs[1:2:end])
            sfy = scalefield(fy, fac)
            sfx = scalefield(fx, fac)
            push!(fields, (mode=2idx-1, fields=(sfy,)))
            push!(fields, (mode=2idx, fields=(sfx,)))
        end
    end
    Tuple(fields)
end

function scalefield(f::Fields.PulseField, fac)
    Fields.PulseField(f.λ0, nmult(f.energy, fac), nmult(f.power, fac), f.ϕ, f.Itshape)
end

function scalefield(f::Fields.DataField, fac)
    Fields.DataField(f.ω, f.Iω, f.ϕω, nmult(f.energy, fac), f.ϕ, f.λ0)
end

_findmode(mode_s, md) = _findmode([mode_s], md)

function makeinputs(mode_s, λ0, pulse::Pulses.AbstractPulse)
    idcs = findmode(mode_s, pulse)
    (length(idcs) > 0) || error("Mode $(pulse.mode) not found in mode list: $mode_s")
    if pulse.polarisation == :linear
        ((mode=idcs[1], fields=(pulse.field,)),)
    else
        (length(idcs) == 2) || error("Modes not set up for circular/elliptical polarisation")
        f1, f2 = ellfields(pulse)
        ((mode=idcs[1], fields=(f1,)), (mode=idcs[2], fields=(f2,)))
    end
end

function makeinputs(mode_s, λ0, pulses::Vector{<:Pulses.AbstractPulse})
    i = Tuple(collect(Iterators.flatten([makeinputs(mode_s, λ0, pii) for pii in pulses])))
    @debug join(string.(i), "\n")
    return i
end

ellphase(ϕ, pol::Symbol) = ellphase(ϕ, 1.0)

function ellphase(ϕ, ε)
    shift = π/2 * sign(ε)
    if length(ϕ) == 0
        return [shift]
    else
        out = copy(ϕ)
        out[1] += shift
        return out
    end
end

ellfac(pol::Symbol) = (1/2, 1/2) # circular
function ellfac(ε::Number)
    (-1 <= ε <= 1) || throw(DomainError(ε, "Ellipticity must be between -1 and 1."))
    (1-ε^2/(1+ε^2), ε^2/(1+ε^2))
end
# sqrt(px/py) = ε => px = ε^2*py; px+py = 1 => px = ε^2*(1-px) => px = ε^2/(1+ε^2)

nmult(x::Nothing, fac) = x
nmult(x, fac) = x*fac

function ellfields(pulse::Union{Pulses.CustomPulse, Pulses.GaussPulse, Pulses.SechPulse})
    f = pulse.field
    py, px = ellfac(pulse.polarisation)
    f1 = Fields.PulseField(f.λ0, nmult(f.energy, py), nmult(f.power, py), f.ϕ, f.Itshape)
    f2 = Fields.PulseField(f.λ0, nmult(f.energy, px), nmult(f.power, px),
                           ellphase(f.ϕ, pulse.polarisation), f.Itshape)
    f1, f2
end

function ellfields(pulse::Pulses.DataPulse)
    f = pulse.field
    py, px = ellfac(pulse.polarisation)
    f1 = Fields.DataField(f.ω, f.Iω, f.ϕω, nmult(f.energy, py), f.ϕ, f.λ0)
    f2 = Fields.DataField(f.ω, f.Iω, f.ϕω, nmult(f.energy, px),
                          ellphase(f.ϕ, pulse.polarisation), f.λ0)
    f1, f2
end

function shotnoise_maybe(inputs, mode::Modes.AbstractMode, shotnoise::Bool)
    shotnoise || return inputs
    (inputs..., (mode=1, fields=(Fields.ShotNoise(),)))
end

function shotnoise_maybe(inputs, modes, shotnoise::Bool)
    shotnoise || return inputs
    (inputs..., [(mode=ii, fields=(Fields.ShotNoise(),)) for ii in eachindex(modes)]...)
end

function setup(grid, mode::Modes.AbstractMode, density, responses, inputs, pol, c::Val{true})
    @info("Using mode-averaged propagation.")
    linop, βfun!, _, _ = LinearOps.make_const_linop(grid, mode, grid.referenceλ)

    Eω, transform, FT = Luna.setup(grid, density, responses, inputs,
    βfun!, z -> Modes.Aeff(mode, z=z))
    linop, Eω, transform, FT
end

function setup(grid, mode::Modes.AbstractMode, density, responses, inputs, pol, c::Val{false})
    @info("Using mode-averaged propagation.")
    linop, βfun! = LinearOps.make_linop(grid, mode, grid.referenceλ)

    Eω, transform, FT = Luna.setup(grid, density, responses, inputs,
                                   βfun!, z -> Modes.Aeff(mode, z=z))
    linop, Eω, transform, FT
end

needfull(modes) = !all(modes) do mode
    (mode.kind == :HE) && (mode.n == 1)
end

function setup(grid, modes, density, responses, inputs, pol, c::Val{true})
    nf = needfull(modes)
    @info(nf ? "Using full 2-D modal integral." : "Using radial modal integral.")
    linop = LinearOps.make_const_linop(grid, modes, grid.referenceλ)
    Eω, transform, FT = Luna.setup(grid, density, responses, inputs, modes,
                                   pol ? :xy : :y; full=nf)
    linop, Eω, transform, FT
end

function setup(grid, modes, density, responses, inputs, pol, c::Val{false})
    nf = needfull(modes)
    @info(nf ? "Using full 2-D modal integral." : "Using radial modal integral.")
    linop = LinearOps.make_linop(grid, modes, grid.referenceλ)
    Eω, transform, FT = Luna.setup(grid, density, responses, inputs, modes,
                                   pol ? :xy : :y; full=nf)
    linop, Eω, transform, FT
end

function makeoutput(grid, saveN, stats, filepath::Nothing, scan::Nothing, scanidx, filename)
    Output.MemoryOutput(0, grid.zmax, saveN, stats)
end

function makeoutput(grid, saveN, stats, filepath, scan::Nothing, scanidx, filename)
    Output.HDF5Output(filepath, 0, grid.zmax, saveN, stats)
end

function makeoutput(grid, saveN, stats, filepath, scan, scanidx, filename)
    isnothing(scanidx) && error("scanidx must be passed along with scan.")
    Output.ScanHDF5Output(scan, scanidx, 0, grid.zmax, saveN, stats;
                          fdir=filepath, fname=filename)
end

end


text:
@testset "Addition" begin
    @test A₀ + B₀ == B₀ + A₀ == (A₀.coef + B₀.coef) * v
    @test A₂ + B₃ == B₃ + A₂
    @test (A₁ + B₂) + B₃ == A₁ + (B₂ + B₃) == A₁ + B₂ + B₃
    @test A₂ - B₃ == -(B₃ - A₂)
end


text:
# transforms marginal univariate distributions

VecVec{T} = Union{Vector{Vector{Int64}}, Vector{Vector{T}}}

"""
  convertmarg!(X::Matrix, d::UnionAll, p::Vector{Vector})

Takes matrix X of realizations of size(X,2) = n dimensional random variable, with
uniform marginals numbered by i, and convert those marginals to common distribution
d with parameters p[i].
If `testunif = true` each marginal is tested for uniformity.

```jldoctest
julia> Random.seed!(43);

julia> x = rand(10,2);

julia> convertmarg!(x, Normal, [[0, 1],[0, 10]])

julia> x
10×2 Array{Float64,2}:
 -0.911655    4.17328
  0.756673  -14.4472
  1.22088   -11.4823
  1.43866   -13.1053
 -0.231978  -11.2415
  1.35696     6.43914
  0.949145  -26.0172
 -0.251957  -18.9723
 -0.177808Real4172
  1.70477    10.4192
```
"""
function convertmarg!(U::Matrix{T}, d::UnionAll, p::VecVec{T} = [fill([0,1], size(U, 2))...];
                                                testunif::Bool = true) where T <: Real
  for i = 1:size(U, 2)
    if testunif
      pvalue(ExactOneSampleKSTest(U[:,i],Uniform(0,1)))>0.0001 || throw(AssertionError("$i marg. not unif."))
    end
    @inbounds U[:,i] = quantile.(d(p[i]...), U[:,i])
  end
end

  # generates covariance matrix

"""
    cormatgen(n::Int = 20)

Returns symmetric correlation matrix Σ of size n x n.
Method:

    a = rand(n,n)
    b = a*a'
    c = b./maximum(b)

Example:

```jldoctest
julia> Random.seed!(43);

julia> cormatgen(2)
2×2 Array{Float64,2}:
 1.0       0.660768
 0.660768  1.0
```
"""
function cormatgen(n::Int = 20)
  a = rand(n,n)
  b = a*a'
  c = b./maximum(b)
  c .- Matrix(Diagonal(c)) .+ Matrix(1.0I, size(c)...)

end

"""
    cormatgen_rand(n::Int = 20)

Returns symmetric correlation matrix Σ of size n x n.
Method:

    a = rand(n,n)
    b = a*a'
    diagb = Matrix(Diagonal(1 ./sqrt.(LinearAlgebra.diag(b))))
    b = diagb*b*diagb

In general gives higher correlations than the cormatgen(). Example:

```jldoctest
julia> Random.seed!(43);

julia> cormatgen_rand(2)
2×2 Array{Float64,2}:
 1.0       0.879086
 0.879086  1.0
```
"""
function cormatgen_rand(n::Int = 20)
  a = rand(n,n)
  b = a*a'
  diagb = Matrix(Diagonal(1 ./sqrt.(LinearAlgebra.diag(b))))
  b = diagb*b*diagb
  (b+b')/2
end

"""
    cormatgen_constant(n::Int, α::Real)

Returns the constant correlation matrix with constant correlations equal to 0 <= α <= 1

```julia
julia> cormatgen_constant(2, 0.5)
2×2 Array{Real,2}:
 1.0  0.5
 0.5  1.0
```
"""
function cormatgen_constant(n::Int, α::Real)
  @assert 0 <= α <= 1 "α should satisfy 0 <= α <= 1"
  α .*ones(n, n) .+(1-α) .*Matrix(1.0I, n,n)
end

function random_unit_vector(dim::Int)
  result = rand(Normal(), dim, 1)
  result /= norm(result)
end

"""
    cormatgen_constant_noised(n::Int, α::Real; ϵ::Real = (1 .-α)/2.)

Returns the constant correlation matrix of size n x n  with constant correlations equal to 0 <= α <= 1
and additinal noise determinde by ϵ.

```julia
julia> Random.seed!(43);

julia> cormatgen_constant_noised(3, 0.5)
3×3 Array{Float64,2}:
 1.0       0.506271  0.285793
 0.506271  1.0       0.475609
 0.285793  0.475609  1.0
```
"""
function cormatgen_constant_noised(n::Int, α::Real; ϵ::Real = (1 .-α)/2.)
  @assert 0 <= ϵ <= 1-α "ϵ must satisfy 0 <= ϵ <= 1-α"
  result = cormatgen_constant(n, α)
  u = hcat([random_unit_vector(n) for i=1:n]...)
  result += ϵ .*(u'*u)
  result - ϵ .*Matrix(1.0I, size(result)...)
end

"""
    cormatgen_two_constant(n::Int, α::Real, β::Real)

Returns the correlation matrix of size n x n of correlations determined by two constants, first must be greater than the second.

```julia
julia> cormatgen_two_constant(6, 0.5, 0.1)
6×6 Array{Float64,2}:
 1.0  0.5  0.5  0.1  0.1  0.1
 0.5  1.0  0.5  0.1  0.1  0.1
 0.5  0.5  1.0  0.1  0.1  0.1
 0.1  0.1  0.1  1.0  0.1  0.1
 0.1  0.1  0.1  0.1  1.0  0.1
 0.1  0.1  0.1  0.1  0.1  1.0

 julia> cormatgen_two_constant(4, 0.5, 0.1)
 4×4 Array{Float64,2}:
  1.0  0.5  0.1  0.1
  0.5  1.0  0.1  0.1
  0.1  0.1  1.0  0.1
  0.1  0.1  0.1  1.0
```
"""
function cormatgen_two_constant(n::Int, α::Real, β::Real)
  @assert α > β "First argument must be greater"
  result = fill(β, (n,n))
  result[1:div(n,2),1:div(n,2)] = fill(α, (div(n,2),div(n,2)))
  result += Matrix(1.0I, size(result)...) - Matrix(Diagonal(result))
  result
end

"""
    cormatgen_two_constant_noised(n::Int, α::Real, β::Real; ϵ::Real= (1-α)/2)

Returns the correlation matrix of size n x n  of correlations determined by two constants, first must be greater than the second.
Additional noise is introduced by the parameter ϵ.

```julia
julia> Random.seed!(43);

julia> cormatgen_two_constant_noised(4, 0.5, 0.1)
4×4 Array{Float64,2}:
  1.0         0.314724   0.190368  -0.0530078
  0.314724    1.0       -0.085744   0.112183
  0.190368   -0.085744   1.0        0.138089
 -0.0530078   0.112183   0.138089   1.0
```
"""
function cormatgen_two_constant_noised(n::Int, α::Real, β::Real; ϵ::Real= (1-α)/2)
  @assert ϵ <= 1-α
  result = cormatgen_two_constant(n, α, β)
  u = hcat([random_unit_vector(n) for i=1:n]...)
  result += ϵ .*(u'*u)
  result - ϵ .*Matrix(1.0I, size(result)...)
end

"""
    cormatgen_toeplitz(n::Int, ρ::Real)

Returns the correlation matrix of size n x n of the Toeplitz structure with
maximal correlation equal to ρ.

```julia
julia> cormatgen_toeplitz(5, 0.9)
5×5 Array{Float64,2}:
 1.0     0.9    0.81  0.729  0.6561
 0.9     1.0    0.9   0.81   0.729
 0.81    0.9    1.0   0.9    0.81
 0.729   0.81   0.9   1.0    0.9
 0.6561  0.729  0.81  0.9    1.0

julia> cormatgen_toeplitz(5, 0.6)
5×5 Array{Float64,2}:
 1.0     0.6    0.36  0.216  0.1296
 0.6     1.0    0.6   0.36   0.216
 0.36    0.6    1.0   0.6    0.36
 0.216   0.36   0.6   1.0    0.6
 0.1296  0.216  0.36  0.6    1.0
```
"""
function cormatgen_toeplitz(n::Int, ρ::Real)
  @assert 0 <= ρ <= 1 "ρ needs to satisfy 0 <= ρ <= 1"
  [ρ^(abs(i-j)) for i=0:n-1, j=0:n-1]
end

"""
    cormatgen_toeplitz_noised(n::Int, ρ::Real; ϵ=(1-ρ)/(1+ρ)/2)

Returns the correlation matrix of size n x n of the Toeplitz structure with
maximal correlation equal to ρ. Additiona noise id added by the ϵ parameter.

```julia
julia> Random.seed!(43);

julia> cormatgen_toeplitz_noised(5, 0.9)
5×5 Array{Float64,2}:
 1.0       0.89656   0.812152  0.720547  0.64318
 0.89656   1.0       0.918136  0.832571  0.734564
 0.812152  0.918136  1.0       0.915888  0.822804
 0.720547  0.832571  0.915888  1.0       0.903819
 0.64318   0.734564  0.822804  0.903819  1.0
```
"""
function cormatgen_toeplitz_noised(n::Int, ρ::Float64; ϵ=(1-ρ)/(1+ρ)/2)
  @assert 0 <= ϵ <= (1-ρ)/(1+ρ) "ϵ must satisfy 0 <= ϵ <= (1-ρ)/(1+ρ)"

  result = cormatgen_toeplitz(n, ρ)
  u = hcat([random_unit_vector(n) for i=1:n]...)
  result += ϵ .*(u'*u)
  result - ϵ .*Matrix(1.0I, size(result)...)
end


text:
# Autogenerated wrapper script for nghttp2_jll for armv6l-linux-musleabihf
export libnghttp2

JLLWrappers.@generate_wrapper_header("nghttp2")
JLLWrappers.@declare_library_product(libnghttp2, "libnghttp2.so.14")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libnghttp2,
        "lib/libnghttp2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    function small()

A small cow, artist unknown

# Example
```jldoctest
julia> cowsay("A calf", cow=Cowsay.small)
 ________
< A calf >
 --------
       \\   ,__,
        \\  (..)____
           (__)    )\\
              ||--|| *
```
"""
function small(;eyes="oo", tongue="  ", thoughts="\\")
if eyes == "oo"
    eyes = ".."
end
the_cow =
"""
       $thoughts   ,__,
        $thoughts  ($eyes)____
           (__)    )\\
            $tongue||--|| *
"""
return the_cow
end


text:
using HTTP
using JSON

struct SingleDiff
    a_lines::UnitRange{Int}
    b_lines::UnitRange{Int}

    a_content::String
    b_content::String
end

struct FileDiff
    filename::String
    diffs::Vector{SingleDiff}
end

struct Diff
    files::Dict{String,FileDiff}
end

function parse_diff(diff::String)
    lines = filter(x -> x != "", split(diff, "\n"))

    files = Dict{String,FileDiff}()

    i = 1
    while i <= length(lines)
        if startswith(lines[i], "diff")
            file_diff, i = parse_file_diff(lines, i)

            @assert !haskey(files, file_diff.filename)
            files[file_diff.filename] = file_diff
        else
            i += 1
        end
    end

    return Diff(files)
end

# Returns a FileDiff, and the next index to look at
function parse_file_diff(lines, i::Int)
    @assert startswith(lines[i], "diff")
    filename = split(lines[i], " ")[3][3:end]

    diffs = Vector{SingleDiff}()

    i += 1
    while i <= length(lines) && !startswith(lines[i], "diff")
        if startswith(lines[i], "@@")
            diff, i = parse_single_diff(lines, i)
            push!(diffs, diff)
        else
            i += 1
        end
    end

    return FileDiff(filename, diffs), i
end

function parse_single_diff(lines, i::Int)
    @assert startswith(lines[i], "@@")

    line_num_tokens = split(lines[i], " ")
    a_line_str = line_num_tokens[2][2:end]
    b_line_str = line_num_tokens[3][2:end]

    a_lines = parse_diff_lines(a_line_str)
    b_lines = parse_diff_lines(b_line_str)

    a_content = ""
    b_content = ""

    i += 1
    while i <= length(lines) && !startswith(lines[i], "diff") && !startswith(lines[i], "@@")
        if startswith(lines[i], " ")
            a_content *= "\n" * lines[i][2:end]
            b_content *= "\n" * lines[i][2:end]
        elseif startswith(lines[i], "-")
            a_content *= "\n" * lines[i][2:end]
        elseif startswith(lines[i], "+")
            b_content *= "\n" * lines[i][2:end]
        else
            @error "Unexpected line in diff" line = lines[i]
        end

        i += 1
    end

    # Need to remove the first character from the diffs because it is an
    # extra newline character.
    return SingleDiff(a_lines, b_lines, a_content[2:end], b_content[2:end]), i
end

function parse_diff_lines(str)
    if findfirst(",", str) === nothing
        start_line = parse(Int, str)
        return start_line:start_line
    else
        start_line, len = parse.(Int, split(str, ","))
        return start_line:start_line+len-1
    end
end

function get_github_diff(url)
    r = HTTP.request("GET", url)
    @info "HTTP status" status = r.status
    diff = String(r.body)
end

function post_code_comment(file, lines, body, url, sha, gha_token)
    params =
        Dict{String,Any}("body" => body, "path" => "src/JuliaGHA.jl", "commit_id" => sha)

    if first(lines) >= last(lines)
        params["line"] = first(lines)
        params["side"] = "RIGHT"
    else
        params["start_line"] = first(lines)
        params["start_side"] = "RIGHT"
        params["line"] = last(lines)
        params["side"] = "RIGHT"
    end

    println(params)

    r = HTTP.request(
        "POST",
        url,
        [
            "Accept" => "application/vnd.github.comfort-fade-preview+json",
            "Content-Type" => "application/json",
            "Authorization" => "token $gha_token",
        ],
        JSON.json(params),
    )
    println(r.status)
    body = String(r.body)
    j = JSON.parse(body)
end

function main(; max_files = 3, max_diffs_per_file = 5)
    if !haskey(ENV, "GITHUB_API_URL")
        @error("Not running on GitHub Actions")
        return false
    end

    format_diff = parse_diff(read(`git diff -U0`, String))

    if length(format_diff.files) > max_files
        @error(
            "Aborting because there are more than $max_files files with changes",
            files_changed = length(format_diff.files)
        )
        return false
    end

    for (file, diff) in format_diff.files
        if length(diff.diffs) > max_diffs_per_file
            @error(
                "Aborting because there are more than $max_diffs_per_file in $file",
                num_diffs = length(diff.diffs)
            )
            return false
        end
    end

    pr_number = split(ENV["GITHUB_REF"], "/")[3]
    url = "$(ENV["GITHUB_SERVER_URL"])/$(ENV["GITHUB_REPOSITORY"])/pull/$pr_number.diff"
    github_diff = parse_diff(get_github_diff(url))

    # Check that each format diff can be accessed in the GitHub diff. If not,
    # then leaving suggusted changes will fail. This is a GitHub limitation.
    for (file, file_diff) in format_diff.files
        if !haskey(github_diff.files, file_diff.filename)
            @error(
                "A file has a formatting change that is not included in the PR diff",
                file = file_diff.filename
            )
            return false
        end

        github_file_diff = github_diff.files[file_diff.filename]
        for diff in file_diff.diffs
            format_lines = diff.a_lines

            lines_contained = false
            for github_single_diff in github_file_diff.diffs
                if issubset(format_lines, github_single_diff.b_lines)
                    lines_contained = true
                    break
                end
            end

            if !lines_contained
                @error(
                    "A file has changes in locations not included in the
                    original diff",
                    file = file_diff.filename
                )
            end
        end
    end

    url = "$(ENV["GITHUB_API_URL"])/repos/$(ENV["GITHUB_REPOSITORY"])/pulls/$pr_number/comments"
    for (file, file_diff) in format_diff.files
        for diff in file_diff.diffs
            @info "Diff" content = diff.b_content lines = diff.a_lines
            body = "```suggestion\n$(diff.b_content)\n```"
            post_code_comment(
                file,
                diff.a_lines,
                body,
                url,
                ENV["GITHUB_SHA"],
                ENV["GITHUB_TOKEN"],
            )
        end
    end

    return true
end
main()


text:
using Test: @testset, @test, @test_throws, @test_broken
using CBinding
using Sockets
using Random
using PLCTag
using PLCTag.libplctag

include(Sys.islinux() ? "with-server.jl" : "without-server.jl")


text:
const wet3d = trues(1, 2, 2)
using SparseArrays

# Macro to create grd in a simple `@Dict` call
# Thanks to Lyndon White (on slack)
macro Dict(vars...)
    kvs = Expr.(:call, :Pair, string.(vars), esc.(vars))
    Expr(:call, :Dict, kvs...)
end

function build_grd()
    # From Archer et al. [2000]
    vtot = 1.292e18     # m³
    area_ocean = 349e12 # m²
    depth_low = 100.0  # m
    depth_high = 250.0 # m
    fraction_area_high = 0.15 #
    # Infer other values
    depth_ocean = vtot / area_ocean

    # Build DY3d
    earth_perimeter = 40e6 # 40'000 km
    DYocean = earth_perimeter / 2 # distance from south pole to north pole
    DYhigh = fraction_area_high * DYocean
    DYlow = (1 - fraction_area_high) * DYocean
    DYT3d = zeros(size(wet3d))
    DYT3d[1, 1, :] .= DYhigh
    DYT3d[1, 2, :] .= DYlow
    # Build DZ3d
    DZT3d = zeros(size(wet3d))
    DZT3d[1, 1, 1] = depth_high
    DZT3d[1, 2, 1] = depth_low
    DZT3d[1, 1, 2] = depth_ocean - depth_high
    DZT3d[1, 2, 2] = depth_ocean - depth_low
    # Build DX3d
    DXT3d = vtot / (depth_ocean * DYocean) * ones(size(wet3d))
    # ZT3d
    ZT3d = cumsum(DZT3d, dims=3) - DZT3d / 2
    # ZW3d (depth at top of box)
    ZW3d = cumsum(DZT3d, dims=3) - DZT3d

    return @Dict DXT3d DYT3d DZT3d ZW3d ZT3d
end
const grd = build_grd()

function build_T(grd)
    # Mixing and overturning (F_HD, etc. see Figure 1 of Archer et al.)
    Mixing = zeros(4,4)
    Mixing[1, 2] = Mixing[2, 1] = 10e6 # between top boxes (1 and 2) (10 Sv = 10e6 cubic meters)
    Mixing[1, 3] = Mixing[3, 1] = 53e6 # between high-lat and deep (1 and 3)
    Mixing[2, 4] = Mixing[4, 2] = 1e6 # between low-lat and deep (2 and 4)
    Mixing[3, 4] = Mixing[4, 3] = 1e10 # (trick) fast mixing between deep boxes (3 and 4)
    overturning = 19e6
    Overturning = falses(4,4)
    Overturning[1, 3] = Overturning[3, 4] = Overturning[4, 2] = Overturning[2, 1] = true # Overturning circulation
    # Build T, as a divergence operator, i.e.,
    # T is positive where tracers are removed.
    T = spzeros(4,4)
    for orig in 1:4, dest in 1:4
        # Overturning part
        if Overturning[orig, dest]
            T[orig, orig] += overturning # add at origin
            T[dest, orig] -= overturning # remove at destination
        end

        if !iszero(Mixing[orig, dest])
            T[orig, orig] += Mixing[orig, dest] # add at origin
            T[dest, orig] -= Mixing[orig, dest] # remove at destination
        end
    end
    v3d = buildv3d(grd)
    v = vec(v3d) # Fine here because every point is wet :)
    V⁻¹ = d₀(v.^(-1))
    T = V⁻¹ * T
    return T
end
const T = build_T(grd)


text:
@testset "Node flux report" begin

fName = "report/nodeflux"

excelFluxReport( mpSim, 12, :in, "A junior", "B junior", "Reserve junior",
    "A senior", "B senior", "Master", filename = fName )
excelFluxReport( mpSim, 12, :out, "A junior", "B junior", "Reserve junior",
    "A senior", "B senior", "Master", filename = fName, overwrite = false )

end  # @testset "Node flux report"


text:
@testset "39.combination-sum.jl" begin
    @test combination_sum_39([2, 3, 6, 7], 7) == [[2, 2, 3], [7]]
    @test combination_sum_39([2, 3, 5], 8) == [[2, 2, 2, 2], [2, 3, 3], [3, 5]]
    @test combination_sum_39([2], 1) == Vector{Int}[]
end


text:
function exists(solver::String)
  return false
end

max_stars_repo_name: polislizarralde/jsos

max_issues_repo_name: polislizarralde/jsos

max_forks_repo_name: polislizarralde/jsos


text:
# Client logic deals with handling control frames, user requesting to send frames, state.
# The ClientProtocol type is defined below entirely synchronously. It takes input via the `handle()`
# function, which is defined for the different input types below. It performs internal logic and
# produces a call to its outbound interface.
#
# The outbound interface is composed of two abstract types `WebSocketHandler` and
# `AbstractWriterTaskProxy`. The concrete implementations will be `TaskProxy` objects, which will
# store the calls (function and arguments) and call it on a target object in another coroutine. This
# means that as long as the channels don't block, the logic will be performed concurrently with the
# the callbacks and writing to the network. This is important because the logic might have to respond
# to ping requests in a timely manner, which it might not be able to do if the callbacks block.
#
# For testing purposes the abstract outbound interface can be replaced with mock objects. This lets us
# test the logic of the WebSocket synchronously, without any asynchronicity or concurrency
# complicating things.

export ClientProtocol

#
# ClientProtocol
#

"Type for the logic of a client WebSocket."
mutable struct ClientProtocol <: AbstractClientProtocol
	# The object to which callbacks should be made. This proxy will make the callbacks
	# asynchronously.
	handler::WebSocketHandler
	# Writes frames to the socket, according to the framing details.
	framewriter::FrameWriter
	# Keeps track of when a pong is expected to be received from the server.
	ponger::AbstractPonger
	# Here we keep data collected when we get a message made up of multiple frames.
	buffer::Vector{UInt8}
	# This stores the type of the multiple frame message. This is the opcode of the first frame,
	# as the following frames have the OPCODE_CONTINUATION opcode.
	buffered_type::Opcode
	# This function cleans up the client when the connection is closed.
	client_cleanup::Function
	# Is the connection close behaviour if a close has been initiated
	closebehaviour::Union{ClosingBehaviour, Nothing}
end

ClientProtocol(handler::WebSocketHandler,
			framewriter::FrameWriter,
	        ponger::AbstractPonger,
			client_cleanup::Function) =
	ClientProtocol(handler, framewriter, ponger, Vector{UInt8}(), OPCODE_TEXT, client_cleanup, nothing)

"The state of the connection."
function protocolstate(p::ClientProtocol)
	if p.closebehaviour != nothing
		protocolstate(p.closebehaviour)
	else
		STATE_OPEN
	end
end

# Requirement
# @5_5_1-4 No frames after Close frame
# @7_2_1-3 Abnormal client initiated closure
#
# By design, we only send text, binary, and ping frames when state is open. When a ClosingBehaviour
# is in effect, then the state is not open. Therefore, we do not send any close frames during a
# connection close.

"Send a single text frame."
function handle(logic::ClientProtocol, req::SendTextFrame)
	if protocolstate(logic) == STATE_OPEN
		send(logic.framewriter, req.isfinal, req.opcode, req.data)
	end
end

"Send a single binary frame."
function handle(logic::ClientProtocol, req::SendBinaryFrame)
	if protocolstate(logic) == STATE_OPEN
		send(logic.framewriter, req.isfinal, req.opcode, req.data)
	end
end

function handle(logic::ClientProtocol, req::ClientPingRequest)
	if protocolstate(logic) == STATE_OPEN
		ping_sent(logic.ponger)
		send(logic.framewriter, true, OPCODE_PING, b"")
	end
end

function handle(p::ClientProtocol, ::PongMissed)
	failtheconnection(p, CLOSE_STATUS_ABNORMAL_CLOSE; reason="Missed expected pongs")
	p.client_cleanup()
end

"Handle a user request to close the WebSocket."
function handle(logic::ClientProtocol, req::CloseRequest)
	if protocolstate(logic) != STATE_OPEN
		clientprotocolinput(logic.closebehaviour, req)
		return
	end

	logic.closebehaviour = ClientInitiatedCloseBehaviour(logic.framewriter, logic.handler)
	closetheconnection(logic.closebehaviour)
end

"The underlying socket was closed. This is sent by the reader."
function handle(logic::ClientProtocol, socketclosed::SocketClosed)
	if protocolstate(logic) != STATE_OPEN
		clientprotocolinput(logic.closebehaviour, socketclosed)
		logic.client_cleanup()
		return
	end

	# Requirement
	# @7_2_1-2 Underlying connection is lost

	failtheconnection(logic, CLOSE_STATUS_ABNORMAL_CLOSE; issocketprobablyup=false)

	logic.client_cleanup()
end

iscontrolframe(frame::Frame) = frame.opcode in [OPCODE_CLOSE, OPCODE_PING, OPCODE_PONG]

"Handle a frame from the server."
function handle(logic::ClientProtocol, req::FrameFromServer)
	# Requirement
	# @6_2-3 Receiving a data frame

	if protocolstate(logic) != STATE_OPEN
		clientprotocolinput(logic.closebehaviour, req)
		return
	end

	# The client must close the connection if it receives a masked frame from the server.
	if req.frame.ismasked
		failtheconnection(logic, CLOSE_STATUS_PROTOCOL_ERROR; reason="Server sent a masked frame")
		return
	end

	# If a reserved bit is set, and no extensions are negotiated, then the client must fail the
	# connection.
	# Note: We do not support extensions yet, so by design no extension has been negotiated.
	if req.frame.rsv1 || req.frame.rsv2 || req.frame.rsv3
		failtheconnection(logic, CLOSE_STATUS_PROTOCOL_ERROR; reason="A reserved bit was set")
		return
	end

	# If a fragmented control frame is received, then the client must fail the connection.
	if iscontrolframe(req.frame) && !req.frame.fin
		failtheconnection(logic, CLOSE_STATUS_PROTOCOL_ERROR; reason="A fragmented control frame was received")
		return
	end

	# If a control frame with payload size > 125 bytes is received, then the client must fail the
	# connection.
	if iscontrolframe(req.frame) && length(req.frame.payload) > 125
		failtheconnection(logic, CLOSE_STATUS_PROTOCOL_ERROR; reason="Control frame payload too large")
		return
	end

	if req.frame.opcode == OPCODE_CLOSE
		handle_close(logic, req.frame)
	elseif req.frame.opcode == OPCODE_PING
		handle_ping(logic, req.frame.payload)
	elseif req.frame.opcode == OPCODE_PONG
		handle_pong(logic, req.frame.payload)
	elseif req.frame.opcode == OPCODE_TEXT
		handle_text(logic, req.frame)
	elseif req.frame.opcode == OPCODE_BINARY
		handle_binary(logic, req.frame)
	elseif req.frame.opcode == OPCODE_CONTINUATION
		handle_continuation(logic, req.frame)
	end
end

#
# Internal closing methods
#

function failtheconnection(p::ClientProtocol, status::CloseStatus;
						   issocketprobablyup=true,
						   reason::String = "")
	p.closebehaviour = FailTheConnectionBehaviour(p.framewriter, p.handler, status;
												  issocketprobablyup = issocketprobablyup,
												  reason = reason)
	closetheconnection(p.closebehaviour)
end

#
# Internal handle functions
#

function handle_close(p::ClientProtocol, frame::Frame)
	# TODO Read actual close status
	p.closebehaviour = ServerInitiatedCloseBehaviour(p.framewriter, p.handler, frame)
	closetheconnection(p.closebehaviour)
end

function handle_ping(logic::ClientProtocol, payload::AbstractVector{UInt8})
	send(logic.framewriter, true, OPCODE_PONG, payload)
end

function handle_pong(logic::ClientProtocol, ::Vector{UInt8})
	pong_received(logic.ponger)
end

function handle_text(logic::ClientProtocol, frame::Frame)
	if frame.fin
		text = String(copy(frame.payload))
		if !isvalid(text)
			failtheconnection(logic, CLOSE_STATUS_INCONSISTENT_DATA; reason="Invalid UTF-8")
			return
		end
		on_text(logic.handler, text)
	else
		start_buffer(logic, frame.payload, OPCODE_TEXT)
	end
end

function handle_binary(logic::ClientProtocol, frame::Frame)
	if frame.fin
		on_binary(logic.handler, frame.payload)
	else
		start_buffer(logic, frame.payload, OPCODE_BINARY)
	end
end

# TODO: What if we get a binary/text frame before we get a final continuation frame?
# Answer: Fail the connection. It's in the specification.
function handle_continuation(logic::ClientProtocol, frame::Frame)
	buffer(logic, frame.payload)
	if frame.fin
		if logic.buffered_type == OPCODE_TEXT
			text = String(logic.buffer)
			if !isvalid(text)
				failtheconnection(logic, CLOSE_STATUS_INCONSISTENT_DATA; reason="Multiframe message with invalid UTF-8")
				return
			end
			on_text(logic.handler, text)
		elseif logic.buffered_type == OPCODE_BINARY
			on_binary(logic.handler, logic.buffer)
			logic.buffer = Vector{UInt8}()
		end
	end
end

function start_buffer(logic::ClientProtocol, payload::AbstractVector{UInt8}, opcode::Opcode)
	logic.buffered_type = opcode
	logic.buffer = copy(payload)
end

function buffer(logic::ClientProtocol, payload::AbstractVector{UInt8})
	append!(logic.buffer, payload)
end

#
# Utilities
#

function masking!(input::AbstractVector{UInt8}, mask::AbstractVector{UInt8})
	m = 1
	for i in 1:length(input)
		input[i] = input[i] ⊻ mask[(m - 1) % 4 + 1]
		m += 1
	end
end


text:
using ArtefactModelling
using Test

@testset "ArtefactModelling.jl" begin
    # Write your tests here.
    @test 1 == Int(1.0)
end