text:
"""
    layoutable(Axis, fig_or_scene; bbox = nothing, kwargs...)

Creates an `Axis` object in the parent `fig_or_scene` which consists of a child scene
with orthographic projection for 2D plots and axis decorations that live in the
parent.
"""
function layoutable(::Type{<:Axis}, fig_or_scene::Union{Figure, Scene}; bbox = nothing, kwargs...)

    topscene = get_topscene(fig_or_scene)

    default_attrs = default_attributes(Axis, topscene).attributes
    theme_attrs = subtheme(topscene, :Axis)
    attrs = merge!(merge!(Attributes(kwargs), theme_attrs), default_attrs)

    @extract attrs (
        title, titlefont, titlesize, titlegap, titlevisible, titlealign, titlecolor,
        xlabel, ylabel, xlabelcolor, ylabelcolor, xlabelsize, ylabelsize,
        xlabelvisible, ylabelvisible, xlabelpadding, ylabelpadding,
        xticklabelsize, xticklabelcolor, yticklabelsize, xticklabelsvisible, yticklabelsvisible,
        xticksize, yticksize, xticksvisible, yticksvisible,
        xticklabelspace, yticklabelspace, yticklabelcolor, xticklabelpad, yticklabelpad,
        xticklabelrotation, yticklabelrotation, xticklabelalign, yticklabelalign,
        xtickalign, ytickalign, xtickwidth, ytickwidth, xtickcolor, ytickcolor,
        xpanlock, ypanlock, xzoomlock, yzoomlock,
        spinewidth, xtrimspine, ytrimspine,
        xgridvisible, ygridvisible, xgridwidth, ygridwidth, xgridcolor, ygridcolor,
        xgridstyle, ygridstyle,
        aspect, halign, valign, xticks, yticks, xtickformat, ytickformat, panbutton,
        xpankey, ypankey, xzoomkey, yzoomkey,
        xaxisposition, yaxisposition,
        bottomspinevisible, leftspinevisible, topspinevisible, rightspinevisible,
        bottomspinecolor, leftspinecolor, topspinecolor, rightspinecolor,
        backgroundcolor,
        xlabelfont, ylabelfont, xticklabelfont, yticklabelfont,
        flip_ylabel, xreversed, yreversed,
        xminorticksvisible, xminortickalign, xminorticksize, xminortickwidth, xminortickcolor, xminorticks,
        yminorticksvisible, yminortickalign, yminorticksize, yminortickwidth, yminortickcolor, yminorticks,
        xminorgridvisible, yminorgridvisible, xminorgridwidth, yminorgridwidth,
        xminorgridcolor, yminorgridcolor, xminorgridstyle, yminorgridstyle,
        limits
    )

    decorations = Dict{Symbol, Any}()

    protrusions = Observable(GridLayoutBase.RectSides{Float32}(0,0,0,0))
    layoutobservables = LayoutObservables{Axis}(attrs.width, attrs.height, attrs.tellwidth, attrs.tellheight, halign, valign, attrs.alignmode;
        suggestedbbox = bbox, protrusions = protrusions)

    # initialize either with user limits, or pick defaults based on scales
    # so that we don't immediately error
    targetlimits = Observable{Rect2f}(defaultlimits(limits[], attrs.xscale[], attrs.yscale[]))
    finallimits = Observable{Rect2f}(targetlimits[])

    # the first thing to do when setting a new scale is
    # resetting the limits because simply through expanding they might be invalid for log
    # but we don't have the axis here yet, so we make this nice and ugly ref for it
    this_axis = Ref{Union{Nothing, Axis}}(nothing)
    onany(attrs.xscale, attrs.yscale) do _, _
        isnothing(this_axis[]) || reset_limits!(this_axis[])
    end

    on(targetlimits) do lims
        # this should validate the targetlimits before anything else happens with them
        # so there should be nothing before this lifting `targetlimits`
        # we don't use finallimits because that's one step later and you
        # already shouldn't set invalid targetlimits (even if they could
        # theoretically be adjusted to fit somehow later?)
        # and this way we can error pretty early
        validate_limits_for_scales(lims, attrs.xscale[], attrs.yscale[])
    end

    scenearea = sceneareanode!(layoutobservables.computedbbox, finallimits, aspect)

    scene = Scene(topscene, px_area=scenearea)

    background = poly!(topscene, scenearea, color = backgroundcolor, strokewidth = 0, inspectable = false)
    translate!(background, 0, 0, -100)
    decorations[:background] = background

    block_limit_linking = Observable(false)

    xaxislinks = Axis[]
    yaxislinks = Axis[]

    xgridnode = Observable(Point2f[])
    xgridlines = linesegments!(
        topscene, xgridnode, linewidth = xgridwidth, visible = xgridvisible,
        color = xgridcolor, linestyle = xgridstyle, inspectable = false
    )
    # put gridlines behind the zero plane so they don't overlay plots
    translate!(xgridlines, 0, 0, -10)
    decorations[:xgridlines] = xgridlines

    xminorgridnode = Observable(Point2f[])
    xminorgridlines = linesegments!(
        topscene, xminorgridnode, linewidth = xminorgridwidth, visible = xminorgridvisible,
        color = xminorgridcolor, linestyle = xminorgridstyle, inspectable = false
    )
    # put gridlines behind the zero plane so they don't overlay plots
    translate!(xminorgridlines, 0, 0, -10)
    decorations[:xminorgridlines] = xminorgridlines

    ygridnode = Observable(Point2f[])
    ygridlines = linesegments!(
        topscene, ygridnode, linewidth = ygridwidth, visible = ygridvisible,
        color = ygridcolor, linestyle = ygridstyle, inspectable = false
    )
    # put gridlines behind the zero plane so they don't overlay plots
    translate!(ygridlines, 0, 0, -10)
    decorations[:ygridlines] = ygridlines

    yminorgridnode = Observable(Point2f[])
    yminorgridlines = linesegments!(
        topscene, yminorgridnode, linewidth = yminorgridwidth, visible = yminorgridvisible,
        color = yminorgridcolor, linestyle = yminorgridstyle, inspectable = false
    )
    # put gridlines behind the zero plane so they don't overlay plots
    translate!(yminorgridlines, 0, 0, -10)
    decorations[:yminorgridlines] = yminorgridlines

    onany(attrs.xscale, attrs.yscale) do xsc, ysc
        scene.transformation.transform_func[] = (xsc, ysc)
    end
    notify(attrs.xscale)

    onany(finallimits, xreversed, yreversed, scene.transformation.transform_func) do lims, xrev, yrev, t
        nearclip = -10_000f0
        farclip = 10_000f0

        left, bottom = Makie.apply_transform(t, Point(minimum(lims)))
        right, top = Makie.apply_transform(t, Point(maximum(lims)))

        leftright = xrev ? (right, left) : (left, right)
        bottomtop = yrev ? (top, bottom) : (bottom, top)

        projection = Makie.orthographicprojection(
            leftright...,
            bottomtop..., nearclip, farclip)

        Makie.set_proj_view!(camera(scene), projection, Makie.Mat4f(Makie.I))
    end

    xaxis_endpoints = lift(xaxisposition, scene.px_area) do xaxisposition, area
        if xaxisposition == :bottom
            bottomline(Rect2f(area))
        elseif xaxisposition == :top
            topline(Rect2f(area))
        else
            error("Invalid xaxisposition $xaxisposition")
        end
    end

    yaxis_endpoints = lift(yaxisposition, scene.px_area) do yaxisposition, area
        if yaxisposition == :left
            leftline(Rect2f(area))
        elseif yaxisposition == :right
            rightline(Rect2f(area))
        else
            error("Invalid yaxisposition $yaxisposition")
        end
    end

    xaxis_flipped = lift(x->x == :top, xaxisposition)
    yaxis_flipped = lift(x->x == :right, yaxisposition)

    xspinevisible = lift(xaxis_flipped, bottomspinevisible, topspinevisible) do xflip, bv, tv
        xflip ? tv : bv
    end
    xoppositespinevisible = lift(xaxis_flipped, bottomspinevisible, topspinevisible) do xflip, bv, tv
        xflip ? bv : tv
    end
    yspinevisible = lift(yaxis_flipped, leftspinevisible, rightspinevisible) do yflip, lv, rv
        yflip ? rv : lv
    end
    yoppositespinevisible = lift(yaxis_flipped, leftspinevisible, rightspinevisible) do yflip, lv, rv
        yflip ? lv : rv
    end
    xspinecolor = lift(xaxis_flipped, bottomspinecolor, topspinecolor) do xflip, bc, tc
        xflip ? tc : bc
    end
    xoppositespinecolor = lift(xaxis_flipped, bottomspinecolor, topspinecolor) do xflip, bc, tc
        xflip ? bc : tc
    end
    yspinecolor = lift(yaxis_flipped, leftspinecolor, rightspinecolor) do yflip, lc, rc
        yflip ? rc : lc
    end
    yoppositespinecolor = lift(yaxis_flipped, leftspinecolor, rightspinecolor) do yflip, lc, rc
        yflip ? lc : rc
    end

    xlims = Observable(xlimits(finallimits[]))
    ylims = Observable(ylimits(finallimits[]))

    on(finallimits) do lims
        nxl = xlimits(lims)
        nyl = ylimits(lims)

        if nxl != xlims[]
            xlims[] = nxl
        end
        if nyl != ylims[]
            ylims[] = nyl
        end
    end

    xaxis = LineAxis(topscene, endpoints = xaxis_endpoints, limits = xlims,
        flipped = xaxis_flipped, ticklabelrotation = xticklabelrotation,
        ticklabelalign = xticklabelalign, labelsize = xlabelsize,
        labelpadding = xlabelpadding, ticklabelpad = xticklabelpad, labelvisible = xlabelvisible,
        label = xlabel, labelfont = xlabelfont, ticklabelfont = xticklabelfont, ticklabelcolor = xticklabelcolor, labelcolor = xlabelcolor, tickalign = xtickalign,
        ticklabelspace = xticklabelspace, ticks = xticks, tickformat = xtickformat, ticklabelsvisible = xticklabelsvisible,
        ticksvisible = xticksvisible, spinevisible = xspinevisible, spinecolor = xspinecolor, spinewidth = spinewidth,
        ticklabelsize = xticklabelsize, trimspine = xtrimspine, ticksize = xticksize,
        reversed = xreversed, tickwidth = xtickwidth, tickcolor = xtickcolor,
        minorticksvisible = xminorticksvisible, minortickalign = xminortickalign, minorticksize = xminorticksize, minortickwidth = xminortickwidth, minortickcolor = xminortickcolor, minorticks = xminorticks, scale = attrs.xscale,
        )
    decorations[:xaxis] = xaxis

    yaxis = LineAxis(topscene, endpoints = yaxis_endpoints, limits = ylims,
        flipped = yaxis_flipped, ticklabelrotation = yticklabelrotation,
        ticklabelalign = yticklabelalign, labelsize = ylabelsize,
        labelpadding = ylabelpadding, ticklabelpad = yticklabelpad, labelvisible = ylabelvisible,
        label = ylabel, labelfont = ylabelfont, ticklabelfont = yticklabelfont, ticklabelcolor = yticklabelcolor, labelcolor = ylabelcolor, tickalign = ytickalign,
        ticklabelspace = yticklabelspace, ticks = yticks, tickformat = ytickformat, ticklabelsvisible = yticklabelsvisible,
        ticksvisible = yticksvisible, spinevisible = yspinevisible, spinecolor = yspinecolor, spinewidth = spinewidth,
        trimspine = ytrimspine, ticklabelsize = yticklabelsize, ticksize = yticksize, flip_vertical_label = flip_ylabel, reversed = yreversed, tickwidth = ytickwidth,
            tickcolor = ytickcolor,
        minorticksvisible = yminorticksvisible, minortickalign = yminortickalign, minorticksize = yminorticksize, minortickwidth = yminortickwidth, minortickcolor = yminortickcolor, minorticks = yminorticks, scale = attrs.yscale,
        )

    decorations[:yaxis] = yaxis

    xoppositelinepoints = lift(scene.px_area, spinewidth, xaxisposition) do r, sw, xaxpos
        if xaxpos == :top
            y = bottom(r)
            p1 = Point2(left(r) - 0.5sw, y)
            p2 = Point2(right(r) + 0.5sw, y)
            [p1, p2]
        else
            y = top(r)
            p1 = Point2(left(r) - 0.5sw, y)
            p2 = Point2(right(r) + 0.5sw, y)
            [p1, p2]
        end
    end

    yoppositelinepoints = lift(scene.px_area, spinewidth, yaxisposition) do r, sw, yaxpos
        if yaxpos == :right
            x = left(r)
            p1 = Point2(x, bottom(r) - 0.5sw)
            p2 = Point2(x, top(r) + 0.5sw)
            [p1, p2]
        else
            x = right(r)
            p1 = Point2(x, bottom(r) - 0.5sw)
            p2 = Point2(x, top(r) + 0.5sw)
            [p1, p2]
        end
    end

    xoppositeline = lines!(topscene, xoppositelinepoints, linewidth = spinewidth,
        visible = xoppositespinevisible, color = xoppositespinecolor, inspectable = false,
        linestyle = nothing)
    decorations[:xoppositeline] = xoppositeline
    translate!(xoppositeline, 0, 0, 20)
    yoppositeline = lines!(topscene, yoppositelinepoints, linewidth = spinewidth,
        visible = yoppositespinevisible, color = yoppositespinecolor, inspectable = false,
        linestyle = nothing)
    decorations[:yoppositeline] = yoppositeline
    translate!(yoppositeline, 0, 0, 20)

    on(xaxis.tickpositions) do tickpos
        pxheight = height(scene.px_area[])
        offset = xaxisposition[] == :bottom ? pxheight : -pxheight
        opposite_tickpos = tickpos .+ Ref(Point2f(0, offset))
        xgridnode[] = interleave_vectors(tickpos, opposite_tickpos)
    end

    on(yaxis.tickpositions) do tickpos
        pxwidth = width(scene.px_area[])
        offset = yaxisposition[] == :left ? pxwidth : -pxwidth
        opposite_tickpos = tickpos .+ Ref(Point2f(offset, 0))
        ygridnode[] = interleave_vectors(tickpos, opposite_tickpos)
    end

    on(xaxis.minortickpositions) do tickpos
        pxheight = height(scene.px_area[])
        offset = xaxisposition[] == :bottom ? pxheight : -pxheight
        opposite_tickpos = tickpos .+ Ref(Point2f(0, offset))
        xminorgridnode[] = interleave_vectors(tickpos, opposite_tickpos)
    end

    on(yaxis.minortickpositions) do tickpos
        pxwidth = width(scene.px_area[])
        offset = yaxisposition[] == :left ? pxwidth : -pxwidth
        opposite_tickpos = tickpos .+ Ref(Point2f(offset, 0))
        yminorgridnode[] = interleave_vectors(tickpos, opposite_tickpos)
    end

    titlepos = lift(scene.px_area, titlegap, titlealign, xaxisposition, xaxis.protrusion) do a,
            titlegap, align, xaxisposition, xaxisprotrusion

        x = if align == :center
            a.origin[1] + a.widths[1] / 2
        elseif align == :left
            a.origin[1]
        elseif align == :right
            a.origin[1] + a.widths[1]
        else
            error("Title align $align not supported.")
        end

        yoffset = top(a) + titlegap + (xaxisposition == :top ? xaxisprotrusion : 0f0)

        Point2(x, yoffset)
    end

    titlealignnode = lift(titlealign) do align
        (align, :bottom)
    end

    titlet = text!(
        topscene, title,
        position = titlepos,
        visible = titlevisible,
        textsize = titlesize,
        align = titlealignnode,
        font = titlefont,
        color = titlecolor,
        space = :data,
        inspectable = false)
    decorations[:title] = titlet

    function compute_protrusions(title, titlesize, titlegap, titlevisible, spinewidth,
                topspinevisible, bottomspinevisible, leftspinevisible, rightspinevisible,
                xaxisprotrusion, yaxisprotrusion, xaxisposition, yaxisposition)

        left, right, bottom, top = 0f0, 0f0, 0f0, 0f0

        if xaxisposition == :bottom
            bottom = xaxisprotrusion
        else
            top = xaxisprotrusion
        end

        titlespace = if !titlevisible || iswhitespace(title)
            0f0
        else
            boundingbox(titlet).widths[2] + titlegap
        end
        top += titlespace

        if yaxisposition == :left
            left = yaxisprotrusion
        else
            right = yaxisprotrusion
        end

        GridLayoutBase.RectSides{Float32}(left, right, bottom, top)
    end

    onany(title, titlesize, titlegap, titlevisible, spinewidth,
            topspinevisible, bottomspinevisible, leftspinevisible, rightspinevisible,
            xaxis.protrusion, yaxis.protrusion, xaxisposition, yaxisposition) do args...
        protrusions[] = compute_protrusions(args...)
    end

    # trigger first protrusions with one of the observables
    title[] = title[]

    # trigger bboxnode so the axis layouts itself even if not connected to a
    # layout
    layoutobservables.suggestedbbox[] = layoutobservables.suggestedbbox[]

    mouseeventhandle = addmouseevents!(scene)
    scrollevents = Observable(ScrollEvent(0, 0))
    keysevents = Observable(KeysEvent(Set()))

    on(scene.events.scroll) do s
        if is_mouseinside(scene)
            scrollevents[] = ScrollEvent(s[1], s[2])
            return Consume(true)
        end
        return Consume(false)
    end

    # TODO this should probably just forward KeyEvent from Makie
    on(scene.events.keyboardbutton) do e
        keysevents[] = KeysEvent(scene.events.keyboardstate)
        return Consume(false)
    end

    interactions = Dict{Symbol, Tuple{Bool, Any}}()

    ax = Axis(fig_or_scene, layoutobservables, attrs, decorations, scene,
        xaxislinks, yaxislinks, targetlimits, finallimits, block_limit_linking,
        mouseeventhandle, scrollevents, keysevents, interactions, Cycler())
    this_axis[] = ax

    function process_event(event)
        for (active, interaction) in values(ax.interactions)
            if active
                maybe_consume = process_interaction(interaction, event, ax)
                maybe_consume == Consume(true) && return Consume(true)
            end
        end
        return Consume(false)
    end

    on(process_event, mouseeventhandle.obs)
    on(process_event, scrollevents)
    on(process_event, keysevents)

    register_interaction!(ax, :rectanglezoom, RectangleZoom(ax))

    register_interaction!(ax, :limitreset, LimitReset())

    register_interaction!(ax,
        :scrollzoom,
        ScrollZoom(0.1, Ref{Any}(nothing), Ref{Any}(0), Ref{Any}(0), 0.2))

    register_interaction!(ax,
        :dragpan,
        DragPan(Ref{Any}(nothing), Ref{Any}(0), Ref{Any}(0), 0.2))

    # these are the user defined limits
    on(limits) do mlims
        reset_limits!(ax)
    end

    # these are the limits that we try to target, but they can be changed for correct aspects
    on(targetlimits) do tlims
        update_linked_limits!(block_limit_linking, xaxislinks, yaxislinks, tlims)
    end

    # compute limits that adhere to the limit aspect ratio whenever the targeted
    # limits or the scene size change, because both influence the displayed ratio
    onany(scene.px_area, targetlimits) do pxa, lims
        adjustlimits!(ax)
    end

    # trigger limit pipeline once, with manual finallimits if they haven't changed from
    # their initial value as they need to be triggered at least once to correctly set up
    # projection matrices etc.
    fl = finallimits[]
    notify(limits)
    if fl == finallimits[]
        notify(finallimits)
    end

    ax
end

"""
    reset_limits!(ax; xauto = true, yauto = true)

Resets the axis limits depending on the value of `ax.limits`.
If one of the two components of limits is nothing,
that value is either copied from the targetlimits if `xauto` or `yauto` is false,
respectively, or it is determined automatically from the plots in the axis.
If one of the components is a tuple of two numbers, those are used directly.
"""
function reset_limits!(ax; xauto = true, yauto = true, zauto = true)
    mlims = convert_limit_attribute(ax.limits[])

    if ax isa Axis
        mxlims, mylims = mlims::Tuple{Any, Any}
    elseif ax isa Axis3
        mxlims, mylims, mzlims = mlims::Tuple{Any, Any, Any}
    else
        error()
    end

    xlims = if isnothing(mxlims) || mxlims[1] === nothing || mxlims[2] === nothing
        l = if xauto
            xautolimits(ax)
        else
            minimum(ax.targetlimits[])[1], maximum(ax.targetlimits[])[1]
        end
        if mxlims === nothing
            l
        else
            lo = mxlims[1] === nothing ? l[1] : mxlims[1]
            hi = mxlims[2] === nothing ? l[2] : mxlims[2]
            (lo, hi)
        end
    else
        convert(Tuple{Float32, Float32}, tuple(mxlims...))
    end
    ylims = if isnothing(mylims) || mylims[1] === nothing || mylims[2] === nothing
        l = if yauto
            yautolimits(ax)
        else
            minimum(ax.targetlimits[])[2], maximum(ax.targetlimits[])[2]
        end
        if mylims === nothing
            l
        else
            lo = mylims[1] === nothing ? l[1] : mylims[1]
            hi = mylims[2] === nothing ? l[2] : mylims[2]
            (lo, hi)
        end
    else
        convert(Tuple{Float32, Float32}, mylims)
    end

    if ax isa Axis3
        zlims = if isnothing(mzlims) || mzlims[1] === nothing || mzlims[2] === nothing
            l = if zauto
                zautolimits(ax)
            else
                minimum(ax.targetlimits[])[3], maximum(ax.targetlimits[])[3]
            end
            if mzlims === nothing
                l
            else
                lo = mzlims[1] === nothing ? l[1] : mzlims[1]
                hi = mzlims[2] === nothing ? l[2] : mzlims[2]
                (lo, hi)
            end
        else
            convert(Tuple{Float32, Float32}, mzlims)
        end
    end

    if !(xlims[1] <= xlims[2])
        error("Invalid x-limits as xlims[1] <= xlims[2] is not met for $xlims.")
    end
    if !(ylims[1] <= ylims[2])
        error("Invalid y-limits as ylims[1] <= ylims[2] is not met for $ylims.")
    end
    if ax isa Axis3
        if !(zlims[1] <= zlims[2])
            error("Invalid z-limits as zlims[1] <= zlims[2] is not met for $zlims.")
        end
    end

    if ax isa Axis
        ax.targetlimits[] = BBox(xlims..., ylims...)
    elseif ax isa Axis3
        ax.targetlimits[] = Rect3f(
            Vec3f(xlims[1], ylims[1], zlims[1]),
            Vec3f(xlims[2] - xlims[1], ylims[2] - ylims[1], zlims[2] - zlims[1]),
        )
    end
    nothing
end

# this is so users can do limits = (left, right, bottom, top)
function convert_limit_attribute(lims::Tuple{Any, Any, Any, Any})
    (lims[1:2], lims[3:4])
end

function convert_limit_attribute(lims::Tuple{Any, Any})
    lims
end
can_be_current_axis(ax::Axis) = true

function validate_limits_for_scales(lims::Rect, xsc, ysc)
    mi = minimum(lims)
    ma = maximum(lims)
    xlims = (mi[1], ma[1])
    ylims = (mi[2], ma[2])

    if !validate_limits_for_scale(xlims, xsc)
        error("Invalid x-limits $xlims for scale $xsc which is defined on the interval $(defined_interval(xsc))")
    end
    if !validate_limits_for_scale(ylims, ysc)
        error("Invalid y-limits $ylims for scale $ysc which is defined on the interval $(defined_interval(ysc))")
    end
    nothing
end

validate_limits_for_scale(lims, scale) = all(x -> x in defined_interval(scale), lims)

palettesyms(cycle::Cycle) = [c[2] for c in cycle.cycle]
attrsyms(cycle::Cycle) = [c[1] for c in cycle.cycle]

function get_cycler_index!(c::Cycler, P::Type)
    if !haskey(c.counters, P)
        c.counters[P] = 1
    else
        c.counters[P] += 1
    end
end

function get_cycle_for_plottype(allattrs, P)::Cycle
    psym = MakieCore.plotsym(P)

    plottheme = Makie.default_theme(nothing, P)

    cdt = Makie.current_default_theme()
    cycle_raw = if haskey(allattrs, :cycle)
        allattrs.cycle[]
    elseif haskey(cdt, psym) && haskey(cdt[psym], :cycle)
        cdt[psym].cycle[]
    else
        haskey(plottheme, :cycle) ? plottheme.cycle[] : nothing
    end

    if isnothing(cycle_raw)
        Cycle([])
    elseif cycle_raw isa Cycle
        cycle_raw
    else
        Cycle(cycle_raw)
    end
end

function add_cycle_attributes!(allattrs, P, cycle::Cycle, cycler::Cycler, palette::Attributes)
    # check if none of the cycled attributes of this plot
    # were passed manually, because we don't use the cycler
    # if any of the cycled attributes were specified manually
    no_cycle_attribute_passed = !any(keys(allattrs)) do key
        any(syms -> key in syms, attrsyms(cycle))
    end

    # check if any attributes were passed as `Cycled` entries
    # because if there were any, these are looked up directly
    # in the cycler without advancing the counter etc.
    manually_cycled_attributes = filter(keys(allattrs)) do key
        allattrs[key][] isa Cycled
    end

    # if there are any manually cycled attributes, we don't do the normal
    # cycling but only look up exactly the passed attributes
    if !isempty(manually_cycled_attributes)
        # an attribute given as Cycled needs to be present in the cycler,
        # otherwise there's no cycle in which to look up a value
        for k in manually_cycled_attributes
            if k ∉ palettesyms(cycle)
                error("Attribute `$k` was passed with an explicit `Cycled` value, but $k is not specified in the cycler for this plot type $P.")
            end
        end

        palettes = [palette[sym][] for sym in palettesyms(cycle)]

        for sym in manually_cycled_attributes
            isym = findfirst(syms -> sym in syms, attrsyms(cycle))
            index = allattrs[sym][].i
            # replace the Cycled values with values from the correct palettes
            # at the index inside the Cycled object
            allattrs[sym] = if cycle.covary
                palettes[isym][mod1(index, length(palettes[isym]))]
            else
                cis = CartesianIndices(Tuple(length(p) for p in palettes))
                n = length(cis)
                k = mod1(index, n)
                idx = Tuple(cis[k])
                isym
                palettes[isym][idx[isym]]
            end
        end

    elseif no_cycle_attribute_passed
        index = get_cycler_index!(cycler, P)

        palettes = [palette[sym][] for sym in palettesyms(cycle)]

        for (isym, syms) in enumerate(attrsyms(cycle))
            for sym in syms
                allattrs[sym] = if cycle.covary
                    palettes[isym][mod1(index, length(palettes[isym]))]
                else
                    cis = CartesianIndices(Tuple(length(p) for p in palettes))
                    n = length(cis)
                    k = mod1(index, n)
                    idx = Tuple(cis[k])
                    palettes[isym][idx[isym]]
                end
            end
        end
    end
end

function Makie.plot!(
        la::Axis, P::Makie.PlotFunc,
        attributes::Makie.Attributes, args...;
        kw_attributes...)

    allattrs = merge(attributes, Attributes(kw_attributes))

    cycle = get_cycle_for_plottype(allattrs, P)
    add_cycle_attributes!(allattrs, P, cycle, la.cycler, la.palette)

    plot = Makie.plot!(la.scene, P, allattrs, args...)

    # some area-like plots basically always look better if they cover the whole plot area.
    # adjust the limit margins in those cases automatically.
    needs_tight_limits(plot) && tightlimits!(la)

    reset_limits!(la)
    plot
end

function Makie.plot!(P::Makie.PlotFunc, ax::Axis, args...; kw_attributes...)
    attributes = Makie.Attributes(kw_attributes)
    Makie.plot!(ax, P, attributes, args...)
end

needs_tight_limits(@nospecialize any) = false
needs_tight_limits(::Union{Heatmap, Image}) = true
function needs_tight_limits(c::Contourf)
    # we know that all values are included and the contourf is rectangular
    # otherwise here it could be in an arbitrary shape
    return c.levels[] isa Int
end

function expandbboxwithfractionalmargins(bb, margins)
    newwidths = bb.widths .* (1f0 .+ margins)
    diffs = newwidths .- bb.widths
    neworigin = bb.origin .- (0.5f0 .* diffs)
    return Rect2f(neworigin, newwidths)
end

limitunion(lims1, lims2) = (min(lims1..., lims2...), max(lims1..., lims2...))

function expandlimits(lims, margin_low, margin_high, scale)
    # expand limits so that the margins are applied at the current axis scale
    limsordered = (min(lims[1], lims[2]), max(lims[1], lims[2]))
    lims_scaled = scale.(limsordered)

    w_scaled = lims_scaled[2] - lims_scaled[1]
    d_low_scaled = w_scaled * margin_low
    d_high_scaled = w_scaled * margin_high
    inverse = Makie.inverse_transform(scale)
    lims = inverse.((lims_scaled[1] - d_low_scaled, lims_scaled[2] + d_high_scaled))

    # guard against singular limits from something like a vline or hline
    if lims[2] - lims[1] ≈ 0
        # this works for log as well
        # we look at the distance to zero in scaled space
        # then try to center the value between that zero and the value
        # that is the same scaled distance away on the other side
        # which centers the singular value optically
        zerodist = abs(scale(lims[1]))

        # for 0 in linear space this doesn't work so here we just expand to -1, 1
        if zerodist ≈ 0 && scale === identity
            lims = (-one(lims[1]), one(lims[1]))
        else
            lims = inverse.(scale.(lims) .+ (-zerodist, zerodist))
        end
    end
    lims
end

function getlimits(la::Axis, dim)
    # find all plots that don't have exclusion attributes set
    # for this dimension
    if !(dim in (1, 2))
        error("Dimension $dim not allowed. Only 1 or 2.")
    end

    function exclude(plot)
        # only use plots with autolimits = true
        to_value(get(plot, dim == 1 ? :xautolimits : :yautolimits, true)) || return true
        # only use visible plots for limits
        return !to_value(get(plot, :visible, true))
    end
    # get all data limits, minus the excluded plots
    boundingbox = Makie.data_limits(la.scene, exclude)
    # if there are no bboxes remaining, `nothing` signals that no limits could be determined
    Makie.isfinite_rect(boundingbox) || return nothing

    # otherwise start with the first box
    mini, maxi = minimum(boundingbox), maximum(boundingbox)
    return (mini[dim], maxi[dim])
end

getxlimits(la::Axis) = getlimits(la, 1)
getylimits(la::Axis) = getlimits(la, 2)

function update_linked_limits!(block_limit_linking, xaxislinks, yaxislinks, tlims)

    thisxlims = xlimits(tlims)
    thisylims = ylimits(tlims)

    # only change linked axis if not prohibited from doing so because
    # we're currently being updated by another axis' link
    if !block_limit_linking[]

        bothlinks = intersect(xaxislinks, yaxislinks)
        xlinks = setdiff(xaxislinks, yaxislinks)
        ylinks = setdiff(yaxislinks, xaxislinks)

        for link in bothlinks
            otherlims = link.targetlimits[]
            if tlims != otherlims
                link.block_limit_linking[] = true
                link.targetlimits[] = tlims
                link.block_limit_linking[] = false
            end
        end

        for xlink in xlinks
            otherlims = xlink.targetlimits[]
            otherxlims = limits(otherlims, 1)
            otherylims = limits(otherlims, 2)
            if thisxlims != otherxlims
                xlink.block_limit_linking[] = true
                xlink.targetlimits[] = BBox(thisxlims[1], thisxlims[2], otherylims[1], otherylims[2])
                xlink.block_limit_linking[] = false
            end
        end

        for ylink in ylinks
            otherlims = ylink.targetlimits[]
            otherxlims = limits(otherlims, 1)
            otherylims = limits(otherlims, 2)
            if thisylims != otherylims
                ylink.block_limit_linking[] = true
                ylink.targetlimits[] = BBox(otherxlims[1], otherxlims[2], thisylims[1], thisylims[2])
                ylink.block_limit_linking[] = false
            end
        end
    end
end

"""
    autolimits!(la::Axis)

Reset manually specified limits of `la` to an automatically determined rectangle, that depends on the data limits of all plot objects in the axis, as well as the autolimit margins for x and y axis.
"""
function autolimits!(ax::Axis)
    ax.limits[] = (nothing, nothing)
    return
end

function autolimits(ax::Axis, dim::Integer)
    # try getting x limits for the axis and then union them with linked axes
    lims = getlimits(ax, dim)

    for link in ax.xaxislinks
        if isnothing(lims)
            lims = getlimits(link, dim)
        else
            newlims = getlimits(link, dim)
            if !isnothing(newlims)
                lims = limitunion(lims, newlims)
            end
        end
    end

    dimsym = dim == 1 ? :x : :y
    scale = getproperty(ax, Symbol(dimsym, :scale))[]
    margin = getproperty(ax.attributes, Symbol(dimsym, :autolimitmargin))[]
    if !isnothing(lims)
        if !validate_limits_for_scale(lims, scale)
            error("Found invalid x-limits $lims for scale $(scale) which is defined on the interval $(defined_interval(scale))")
        end
        lims = expandlimits(lims, margin[1], margin[2], scale)
    end

    # if no limits have been found, use the targetlimits directly
    if isnothing(lims)
        lims = limits(ax.targetlimits[], dim)
    end
    return lims
end

xautolimits(ax::Axis) = autolimits(ax, 1)
yautolimits(ax::Axis) = autolimits(ax, 2)

"""
    linkaxes!(a::Axis, others...)

Link both x and y axes of all given `Axis` so that they stay synchronized.
"""
function linkaxes!(a::Axis, others...)
    linkxaxes!(a, others...)
    linkyaxes!(a, others...)
end

function adjustlimits!(la)
    asp = la.autolimitaspect[]
    target = la.targetlimits[]

    # in the simplest case, just update the final limits with the target limits
    if isnothing(asp)
        la.finallimits[] = target
        return
    end

    area = la.scene.px_area[]
    xlims = (left(target), right(target))
    ylims = (bottom(target), top(target))

    size_aspect = width(area) / height(area)
    data_aspect = (xlims[2] - xlims[1]) / (ylims[2] - ylims[1])

    aspect_ratio = data_aspect / size_aspect

    correction_factor = asp / aspect_ratio

    if correction_factor > 1
        # need to go wider

        marginsum = sum(la.xautolimitmargin[])
        ratios = if marginsum == 0
            (0.5, 0.5)
        else
            (la.xautolimitmargin[] ./ marginsum)
        end

        xlims = expandlimits(xlims, ((correction_factor - 1) .* ratios)..., identity) # don't use scale here?
    elseif correction_factor < 1
        # need to go taller

        marginsum = sum(la.yautolimitmargin[])
        ratios = if marginsum == 0
            (0.5, 0.5)
        else
            (la.yautolimitmargin[] ./ marginsum)
        end
        ylims = expandlimits(ylims, (((1 / correction_factor) - 1) .* ratios)..., identity) # don't use scale here?
    end

    bbox = BBox(xlims[1], xlims[2], ylims[1], ylims[2])
    la.finallimits[] = bbox
    return
end

"""
    linkxaxes!(a::Axis, others...)

Link the x axes of all given `Axis` so that they stay synchronized.
"""
function linkxaxes!(a::Axis, others...)
    axes = Axis[a; others...]

    for i in 1:length(axes)-1
        for j in i+1:length(axes)
            axa = axes[i]
            axb = axes[j]

            if axa ∉ axb.xaxislinks
                push!(axb.xaxislinks, axa)
            end
            if axb ∉ axa.xaxislinks
                push!(axa.xaxislinks, axb)
            end
        end
    end
    # update limits because users will expect to see the effect
    reset_limits!(a)
end

"""
    linkyaxes!(a::Axis, others...)

Link the y axes of all given `Axis` so that they stay synchronized.
"""
function linkyaxes!(a::Axis, others...)
    axes = Axis[a; others...]

    for i in 1:length(axes)-1
        for j in i+1:length(axes)
            axa = axes[i]
            axb = axes[j]

            if axa ∉ axb.yaxislinks
                push!(axb.yaxislinks, axa)
            end
            if axb ∉ axa.yaxislinks
                push!(axa.yaxislinks, axb)
            end
        end
    end
    # update limits because users will expect to see the effect
    reset_limits!(a)
end

"""
Keeps the ticklabelspace static for a short duration and then resets it to its previous
value. If that value is Makie.automatic, the reset will trigger new
protrusions for the axis and the layout will adjust. This is so the layout doesn't
immediately readjust during interaction, which would let the whole layout jitter around.
"""
function timed_ticklabelspace_reset(ax::Axis, reset_timer::Ref,
        prev_xticklabelspace::Ref, prev_yticklabelspace::Ref, threshold_sec::Real)

    if !isnothing(reset_timer[])
        close(reset_timer[])
    else
        prev_xticklabelspace[] = ax.xticklabelspace[]
        prev_yticklabelspace[] = ax.yticklabelspace[]

        ax.xticklabelspace = ax.elements[:xaxis].attributes.actual_ticklabelspace[]
        ax.yticklabelspace = ax.elements[:yaxis].attributes.actual_ticklabelspace[]
    end

    reset_timer[] = Timer(threshold_sec) do t
        reset_timer[] = nothing

        ax.xticklabelspace = prev_xticklabelspace[]
        ax.yticklabelspace = prev_yticklabelspace[]
    end

end

"""
    hidexdecorations!(la::Axis; label = true, ticklabels = true, ticks = true, grid = true,
        minorgrid = true, minorticks = true)

Hide decorations of the x-axis: label, ticklabels, ticks and grid.
"""
function hidexdecorations!(la::Axis; label = true, ticklabels = true, ticks = true, grid = true,
        minorgrid = true, minorticks = true)
    if label
        la.xlabelvisible = false
    end
    if ticklabels
        la.xticklabelsvisible = false
    end
    if ticks
        la.xticksvisible = false
    end
    if grid
        la.xgridvisible = false
    end
    if minorgrid
        la.xminorgridvisible = false
    end
    if minorticks
        la.xminorticksvisible = false
    end
end

"""
    hideydecorations!(la::Axis; label = true, ticklabels = true, ticks = true, grid = true,
        minorgrid = true, minorticks = true)

Hide decorations of the y-axis: label, ticklabels, ticks and grid.
"""
function hideydecorations!(la::Axis; label = true, ticklabels = true, ticks = true, grid = true,
        minorgrid = true, minorticks = true)
    if label
        la.ylabelvisible = false
    end
    if ticklabels
        la.yticklabelsvisible = false
    end
    if ticks
        la.yticksvisible = false
    end
    if grid
        la.ygridvisible = false
    end
    if minorgrid
        la.yminorgridvisible = false
    end
    if minorticks
        la.yminorticksvisible = false
    end
end

"""
    hidedecorations!(la::Axis)

Hide decorations of both x and y-axis: label, ticklabels, ticks and grid.
"""
function hidedecorations!(la::Axis; label = true, ticklabels = true, ticks = true, grid = true,
        minorgrid = true, minorticks = true)
    hidexdecorations!(la; label = label, ticklabels = ticklabels, ticks = ticks, grid = grid,
        minorgrid = minorgrid, minorticks = minorticks)
    hideydecorations!(la; label = label, ticklabels = ticklabels, ticks = ticks, grid = grid,
        minorgrid = minorgrid, minorticks = minorticks)
end

"""
    hidespines!(la::Axis, spines::Symbol... = (:l, :r, :b, :t)...)

Hide all specified axis spines. Hides all spines by default, otherwise choose
with the symbols :l, :r, :b and :t.
"""
function hidespines!(la::Axis, spines::Symbol... = (:l, :r, :b, :t)...)
    for s in spines
        @match s begin
            :l => (la.leftspinevisible = false)
            :r => (la.rightspinevisible = false)
            :b => (la.bottomspinevisible = false)
            :t => (la.topspinevisible = false)
            x => error("Invalid spine identifier $x. Valid options are :l, :r, :b and :t.")
        end
    end
end

function tight_yticklabel_spacing!(la::Axis)
    tight_ticklabel_spacing!(la.elements[:yaxis])
end

function tight_xticklabel_spacing!(la::Axis)
    tight_ticklabel_spacing!(la.elements[:xaxis])
end

function tight_ticklabel_spacing!(la::Axis)
    tight_xticklabel_spacing!(la)
    tight_yticklabel_spacing!(la)
end

function Base.show(io::IO, ::MIME"text/plain", ax::Axis)
    nplots = length(ax.scene.plots)
    println(io, "Axis with $nplots plots:")

    for (i, p) in enumerate(ax.scene.plots)
        println(io, (i == nplots ? " ┗━ " : " ┣━ ") * string(typeof(p)))
    end
end

function Base.show(io::IO, ax::Axis)
    nplots = length(ax.scene.plots)
    print(io, "Axis ($nplots plots)")
end

function Makie.xlims!(ax::Axis, xlims)
    if length(xlims) != 2
        error("Invalid xlims length of $(length(xlims)), must be 2.")
    elseif xlims[1] == xlims[2]
        error("Can't set x limits to the same value $(xlims[1]).")
    elseif all(x -> x isa Real, xlims) && xlims[1] > xlims[2]
        xlims = reverse(xlims)
        ax.xreversed[] = true
    else
        ax.xreversed[] = false
    end

    ax.limits.val = (xlims, ax.limits[][2])
    reset_limits!(ax, yauto = false)
    nothing
end

function Makie.ylims!(ax::Axis, ylims)
    if length(ylims) != 2
        error("Invalid ylims length of $(length(ylims)), must be 2.")
    elseif ylims[1] == ylims[2]
        error("Can't set y limits to the same value $(ylims[1]).")
    elseif all(x -> x isa Real, ylims) && ylims[1] > ylims[2]
        ylims = reverse(ylims)
        ax.yreversed[] = true
    else
        ax.yreversed[] = false
    end

    ax.limits.val = (ax.limits[][1], ylims)
    reset_limits!(ax, xauto = false)
    nothing
end

Makie.xlims!(ax, low, high) = Makie.xlims!(ax, (low, high))
Makie.ylims!(ax, low, high) = Makie.ylims!(ax, (low, high))
Makie.zlims!(ax, low, high) = Makie.zlims!(ax, (low, high))

Makie.xlims!(low::Optional{<:Real}, high::Optional{<:Real}) = Makie.xlims!(current_axis(), low, high)
Makie.ylims!(low::Optional{<:Real}, high::Optional{<:Real}) = Makie.ylims!(current_axis(), low, high)
Makie.zlims!(low::Optional{<:Real}, high::Optional{<:Real}) = Makie.zlims!(current_axis(), low, high)

Makie.xlims!(ax = current_axis(); low = nothing, high = nothing) = Makie.xlims!(ax, low, high)
Makie.ylims!(ax = current_axis(); low = nothing, high = nothing) = Makie.ylims!(ax, low, high)
Makie.zlims!(ax = current_axis(); low = nothing, high = nothing) = Makie.zlims!(ax, low, high)

"""
    limits!(ax::Axis, xlims, ylims)

Set the axis limits to `xlims` and `ylims`.
If limits are ordered high-low, this reverses the axis orientation.
"""
function limits!(ax::Axis, xlims, ylims)
    Makie.xlims!(ax, xlims)
    Makie.ylims!(ax, ylims)
end

"""
    limits!(ax::Axis, x1, x2, y1, y2)

Set the axis x-limits to `x1` and `x2` and the y-limits to `y1` and `y2`.
If limits are ordered high-low, this reverses the axis orientation.
"""
function limits!(ax::Axis, x1, x2, y1, y2)
    Makie.xlims!(ax, x1, x2)
    Makie.ylims!(ax, y1, y2)
end

"""
    limits!(ax::Axis, rect::Rect2)

Set the axis limits to `rect`.
If limits are ordered high-low, this reverses the axis orientation.
"""
function limits!(ax::Axis, rect::Rect2)
    xmin, ymin = minimum(rect)
    xmax, ymax = maximum(rect)
    Makie.xlims!(ax, xmin, xmax)
    Makie.ylims!(ax, ymin, ymax)
end

function limits!(args...)
    limits!(current_axis(), args...)
end

function Base.delete!(ax::Axis, plot::AbstractPlot)
    delete!(ax.scene, plot)
    ax
end

function Base.empty!(ax::Axis)
    for plot in copy(ax.scene.plots)
        delete!(ax, plot)
    end
    ax
end

Makie.transform_func(ax::Axis) = Makie.transform_func(ax.scene)

# these functions pick limits for different x and y scales, so that
# we don't pick values that are invalid, such as 0 for log etc.
function defaultlimits(userlimits::Tuple{Real, Real, Real, Real}, xscale, yscale)
    BBox(userlimits...)
end

defaultlimits(l::Tuple{Any, Any, Any, Any}, xscale, yscale) = defaultlimits(((l[1], l[2]), (l[3], l[4])), xscale, yscale)

function defaultlimits(userlimits::Tuple{Any, Any}, xscale, yscale)
    xl = defaultlimits(userlimits[1], xscale)
    yl = defaultlimits(userlimits[2], yscale)
    BBox(xl..., yl...)
end

defaultlimits(limits::Nothing, scale) = defaultlimits(scale)
defaultlimits(limits::Tuple{Real, Real}, scale) = limits
defaultlimits(limits::Tuple{Real, Nothing}, scale) = (limits[1], defaultlimits(scale)[2])
defaultlimits(limits::Tuple{Nothing, Real}, scale) = (defaultlimits(scale)[1], limits[2])
defaultlimits(limits::Tuple{Nothing, Nothing}, scale) = defaultlimits(scale)

defaultlimits(::typeof(log10)) = (1.0, 1000.0)
defaultlimits(::typeof(log2)) = (1.0, 8.0)
defaultlimits(::typeof(log)) = (1.0, exp(3.0))
defaultlimits(::typeof(identity)) = (0.0, 10.0)
defaultlimits(::typeof(sqrt)) = (0.0, 100.0)
defaultlimits(::typeof(Makie.logit)) = (0.01, 0.99)
defaultlimits(::typeof(Makie.pseudolog10)) = (0.0, 100.0)
defaultlimits(::Makie.Symlog10) = (0.0, 100.0)

defined_interval(::typeof(identity)) = OpenInterval(-Inf, Inf)
defined_interval(::Union{typeof(log2), typeof(log10), typeof(log)}) = OpenInterval(0.0, Inf)
defined_interval(::typeof(sqrt)) = Interval{:closed,:open}(0, Inf)
defined_interval(::typeof(Makie.logit)) = OpenInterval(0.0, 1.0)
defined_interval(::typeof(Makie.pseudolog10)) = OpenInterval(-Inf, Inf)
defined_interval(::Makie.Symlog10) = OpenInterval(-Inf, Inf)


text:
# place here configuration options that will be set for all environments
config.app_is_api = false

ext: jl

max_stars_repo_name: essenciary/genie-todo-mvc

max_issues_repo_name: essenciary/genie-todo-mvc

max_forks_repo_name: essenciary/genie-todo-mvc


text:
using DifferentialEquations
using ModelingToolkit
# using DataDrivenDiffEq
using LinearAlgebra, DiffEqSensitivity, Optim
using DiffEqFlux, Flux
using Plots, Revise
gr()

function triNN!(du,u,p,t,dens,cons)
    # unpack rates and constants
    nᵣ,nₓ,n₃ = u
    k₁,k₋₁,k₂,k₋₂ = cons
    mᵣ,mₗ,mₓ,A = dens
    z = L(u,p)
    # model
    du[1] = dnᵣ = A*k₁*mᵣ*mₗ - k₋₁*nᵣ - z[1] + z[2]
    du[2] = dnₓ = A*k₂*mₓ*mₗ - k₋₂*nₓ - z[3] + z[4]
    du[3] = dn₃ = z[1] + z[3] - z[2] - z[4]

end
function triNN!(du,u,p,t,dens,cons)
    # unpack rates and constants
    nᵣ,nₓ,n₃ = u
    k₁,k₋₁,k₂,k₋₂ = cons
    mᵣ,mₗ,mₓ,A = dens
    z =  L(u,p)
    # model
    du[1] = dnᵣ = A*k₁*mᵣ*mₗ - k₋₁*nᵣ - z[1] + z[2]
    du[2] = dnₓ = A*k₂*mₓ*mₗ - k₋₂*nₓ - z[3] + z[4]
    du[3] = dn₃ = z[1] + z[3] - z[2] - z[4]

end
function triNN!(du,u,p,t,dens,cons)
    # unpack rates and constants
    nᵣ,nₓ,n₃ = u
    k₁,k₋₁,k₂,k₋₂ = cons
    mᵣ,mₗ,mₓ,A = dens
    z = L(u,p)
    # model
    du[1] = dnᵣ = A*k₁*mᵣ*mₗ - k₋₁*nᵣ + z[1]
    du[2] = dnₓ = A*k₂*mₓ*mₗ - k₋₂*nₓ + z[2]
    du[3] = dn₃ = z[1] + z[2]

end

L = FastChain(FastDense(3, 32, tanh),FastDense(32, 32, tanh), FastDense(32, 2))
L = FastChain(FastDense(3, 32, tanh),FastDense(32, 50, tanh),FastDense(50, 32, tanh), FastDense(32, 3))

# Define the experimental parameter
tspan = (0.0,16.1)
u0 = Float32[0.0,0.0,0.0]

# Xₙ = n_tc1 #load data
Xₙ = [0.02020270731751947
 0.0
 0.06187540371808745
 0.10536051565782635
 0.083381608939051
 0.19845093872383823
 0.3011050927839216
 0.3011050927839216
 0.3566749439387324
 0.38566248081198473
 0.5108256237659907
 0.6931471805599453
 0.8209805520698303
 0.7339691750802005
 0.8209805520698303
 0.7339691750802005]
 t = [0.25
  0.25
  0.5
  0.5
  1.0
  1.0
  2.0
  2.0
  4.0
  4.0
  6.0
  6.0
  8.0
  8.0
 16.0
 16.0]
densities_tc1 = [15.0,38.0,10.0,1.0]
cons_tc = [ 0.0006692541890287495,0.8662361534770547,1.169965568192585e-6,0.4]

function predict(θ)

    tmp_prob = remake(prob_nn,u0=u0,p=θ)
    tmp_sol =  solve(tmp_prob, VCABM(), saveat = t,
                  abstol=1e-5, reltol=1e-5)
                  # sensealg=BacksolveAdjoint())
                  # backsolve=true)

    Σ_sol = sum(Array(tmp_sol),dims=1)
end

function loss(θ)
    pred = predict(θ)
    # L = sum(abs2, (Xₙ .- pred').*t), pred
    L = sum(abs2, Xₙ .- pred'), pred
    # prediction = pred
    # return L, prediction
end

callback(θ,l,pred) = begin
    push!(losses, l)
    @show l
    if length(losses)%10==0
        println("Current loss after $(length(losses)) iterations: $(losses[end])")
        fig = plot(t,pred')
        scatter!(fig,t,Xₙ)
        display(fig)
    end
    false
end

# function run_model(learning_rate)
p = initial_params(L)
f = (du,u,p,t) -> triNN!(du,u, p,t,densities_tc1,cons_tc)
prob_nn = ODEProblem(f,u0, tspan, p)
sol_nn = solve(prob_nn, Rodas4(), saveat = t)

# plot(solution)
plot(sol_nn)
summ = reduce(vcat,sum(sol_nn,dims=1))
h = plot!(sol_nn.t,summ,linecolor=:black)
# Test
@time loss(p)
# L,predition = loss(p)

# p = initial_params(L)
losses = []
# First train with ADAM for better convergence
res1 = DiffEqFlux.sciml_train(loss, p, ADAM(0.03), cb=callback, maxiters = 300)

# Train with BFGS
res2 = DiffEqFlux.sciml_train(loss, res1.minimizer, BFGS(initial_stepnorm=0.001),
                                cb=callback, maxiters = 100)

println("Final training loss after $(length(losses)) iterations: $(losses[end])")

# Plot the losses
plot(losses, yaxis = :log, xaxis = :log, xlabel = "Iterations", ylabel = "Loss")
end

# run_model(0.005)
# Plot the data and the approximation
NNsolution = predict(res2.minimizer)
# Trained on noisy data vs real solution
plot(t, NNsolution')
scatter!(t, Xₙ)

tmp_prob = remake(prob_nn,u0=u0,p=res2.minimizer)
tmp_sol =  solve(tmp_prob, Tsit5(), saveat = 0.1,
              abstol=1e-5, reltol=1e-5,
              sensealg=BacksolveAdjoint())
              # backsolve=true)

Σ_sol = sum(Array(tmp_sol),dims=1)

plot(tmp_sol.t,Σ_sol')

max_stars_repo_path: neuralODE modeling/NNode_try1.jl

max_issues_repo_path: neuralODE modeling/NNode_try1.jl

max_forks_repo_path: neuralODE modeling/NNode_try1.jl


text:
using Test
using DandelionWebSockets.Proxy
import DandelionWebSockets.Proxy: write, close

struct MockWriter <: IO
    channel::Channel{Frame}
    closechannel::Channel{Bool}
    exception::Union{Exception, Nothing}

    MockWriter() = new(Channel{Frame}(0), Channel{Bool}(0), nothing)
    MockWriter(ex::Exception) = new(Channel{Frame}(0), Channel{Bool}(0), ex)
end

function write(m::MockWriter, frame::Frame)
    if m.exception == nothing
        put!(m.channel, frame)
    else
        throw(get(m.exception))
    end
end

takeframe!(m::MockWriter) = take!(m.channel)
close(m::MockWriter) = put!(m.closechannel, true)
isclosed!(m::MockWriter) = take!(m.closechannel)

@testset "WriterProxy            " begin
    @testset "WriterProxy works in a separate task" begin
        mockwriter = MockWriter()
        proxywriter = WriterProxy(mockwriter)

        frame = Frame(true, 0, 0, 0, OPCODE_TEXT, true, 1, 0, b"\x01\x02\x03\x04", b"1")
        write(proxywriter, frame)

        written_frame = takeframe!(mockwriter)
        @test written_frame == frame
    end

    @testset "Stop WriterProxy" begin
        mockwriter = MockWriter()
        proxywriter = WriterProxy(mockwriter)

        stopproxy(proxywriter)

        @test isopen(proxywriter.channel) == false
    end

    @testset "Close the socket via WriterProxy" begin
        mockwriter = MockWriter()
        proxywriter = WriterProxy(mockwriter)

        close(proxywriter)

        @test isclosed!(mockwriter)
    end

    @testset "Streams throws exception on write; Socket is closed" begin
        writer = MockWriter(EOFError())
        proxywriter = WriterProxy(writer)

        frame = Frame(true, 0, 0, 0, OPCODE_TEXT, true, 1, 0, b"\x01\x02\x03\x04", b"1")
        write(proxywriter, frame)

        @test isclosed!(writer)
    end
end


text:
module TestSingularLRRL

using Modia
@usingModiaPlot

SingularLRRL = Model(
    v0   = 10,
    L1_L = 0.1,
    L2_L = 0.2,
    R1_R = 10,
    R2_R = 20,

    L2_p_i = Var(init = 0.1),
    R2_p_i = Var(start = 0.0),
    V_n_v = Var(start = 0.0),

    equations =:[
        # Inductor 1
        L1_v = L1_p_v - L1_n_v
        0 = L1_p_i + L1_n_i
        L1_L*der(L1_p_i) = L1_v

        # Inductor 2
        L2_v = L2_p_v - L2_n_v
        0 = L2_p_i + L2_n_i
        L2_L*der(L2_p_i) = L2_v

        # Resistor 1
        R1_v = R1_p_v - R1_n_v
        0 = R1_p_i + R1_n_i
        R1_v = R1_R*R1_p_i

        # Resistor 2
        R2_v = R2_p_v - R2_n_v
        0 = R2_p_i + R2_n_i
        R2_v = R2_R*R2_p_i

        # Voltage source
        V_v = V_p_v - V_n_v
        0   = V_p_i + V_n_i

        # Connect equations
        #   connect(V.p, L1.p)
        #   connect(L1.n, R1.p)
        #   connect(L1.n, R2.p)
        #   connect(R1.n, L2.p)
        #   connect(R2.n, L2.p)
        #   connect(L2.n, V.n)
        V_p_v = L1_p_v
        0 = V_p_i + L1_p_i
        L1_n_v = R1_p_v
        L1_n_v = R2_p_v
        0 = L1_n_i + R1_p_i + R2_p_i
        R1_n_v = L2_p_v
        R2_n_v = L2_p_v
        0 = R1_n_i + R2_n_i + L2_p_i
        L2_n_v = V_n_v
        0 = L2_n_i + V_n_i #+ ground_v_i
        #V_n_v = 0
        V_v = v0 ]
)

singularLRRL = @instantiateModel(SingularLRRL, unitless=true)

simulate!(singularLRRL, Tsit5(), stopTime = 1.0, requiredFinalStates = [1.4999999937318995])

plot(singularLRRL, [("L1_p_i", "L2_p_i", "R1_p_i", "R2_p_i"), ("L1_v", "L2_v", "R1_v", "R2_v")])

include("../models/Electric.jl")

setLogMerge(false)

SingularLRRL2 = Model(
    R1 = Resistor | Map(R=10u"Ω"),
    R2 = Resistor | Map(R=20u"Ω", start = Map(v=0.0u"V")),
    L1 = Inductor | Map(L=0.1u"H", i = Var(init=nothing)),
    L2 = Inductor | Map(L=0.2u"H", start = Map(v=0.0u"V")),
    V = ConstantVoltage | Map(V=10u"V"),
    connections = :[
        (V.p, L1.p)
        (L1.n, R1.p, R2.p)
        (L2.p, R1.n, R2.n)
        (L2.n, V.n) ]
)

#=
singularLRRL2 = @instantiateModel(SingularLRRL2, unitless=true, log=true, logCode=true, logExecution=true)

simulate!(singularLRRL2, Tsit5(), stopTime = 1.0)

plot(singularLRRL2, [("L1.p.i", "L2.p.i", "R1.p.i", "R2.p.i"), ("L1.v", "L2.v", "R1.v", "R2.v")])
=#

SingularLRRL3 = Model(
    R1 = Resistor | Map(R=1),
    R2 = Resistor | Map(R=2),
    R3 = Resistor | Map(R=3),
    R4 = Resistor | Map(R=4),
    L1 = Inductor | Map(L=0.1),
    L2 = Inductor | Map(L=0.2, i = Var(init=nothing)),
    V = ConstantVoltage | Map(V=10),
#    ground = Ground,
    connections = :[
        (V.p, L1.p)
        (L1.n, R1.p, R2.p)
        (R1.n, R2.n, R3.p, R4.p)
        (L2.p, R3.n, R4.n)
        (L2.n, V.n) ] #, ground.p) ]
)

singularLRRL3 = @instantiateModel(SingularLRRL3, unitless=true, log=false, logDetails=false)

simulate!(singularLRRL3, Tsit5(), stopTime = 1.0, requiredFinalStates = [4.198498632779839])

plot(singularLRRL3, [("L1.i", "L2.i", "R1.i", "R2.i", "R3.i", "R4.i"), ("L1.v", "L2.v", "R1.v", "R2.v", "R3.v", "R4.v")])

end


text:
function save_mesh(::Val{:FLUENT}, io::IO, m::AirfoilMesh)
    NI, NJ = size(m.grid)
    NV = NI * NJ - m.ncon
    NF = NI * (NJ - 1) + NJ * (NI - 1) - max(m.ncon - 1, 0)
    NC = (NI - 1) * (NJ - 1)

    NLI = Matrix(LinearIndices((NI, NJ)))
    NCI = Dict{Int,CartesianIndex{2}}()
    # try do remove duplicated
    if m.ncon > 0
        if m.topo == :O
            @assert m.ncon == NJ
            NLI[:, NJ] .= NLI[:, 1]
        elseif m.topo == :C
            NLI[NI:-1:(NI - m.ncon + 1), 1] .= NLI[1:(m.ncon), 1]
        end
        old_indicies = sort(collect(Set(NLI)))
        @assert length(old_indicies) == NV
        o2n = Dict{Int,Int}()
        for i in eachindex(old_indicies)
            o2n[old_indicies[i]] = i
        end
        for i in eachindex(NLI)
            NLI[i] = o2n[NLI[i]]
        end
    end
    for j in 1:NJ
        for i in 1:NI
            li = NLI[i, j]
            if !haskey(NCI, li)
                NCI[li] = CartesianIndex(i, j)
            end
        end
    end

    CLI = LinearIndices((NI - 1, NJ - 1))
    f_int = Vector{NTuple{4,Int}}()
    f_wall = Vector{NTuple{4,Int}}()
    f_far = Vector{NTuple{4,Int}}()
    f_min = Vector{NTuple{4,Int}}()
    f_max = Vector{NTuple{4,Int}}()
    # INTERIOR
    for j in 1:(NJ - 1)
        for i in 2:(NI - 1)
            p0 = NLI[i, j]
            p1 = NLI[i, j + 1]
            c0 = CLI[i - 1, j]
            c1 = CLI[i, j]
            push!(f_int, (p0, p1, c0, c1))
        end
    end
    for j in 2:(NJ - 1)
        for i in 1:(NI - 1)
            p0 = NLI[i, j]
            p1 = NLI[i + 1, j]
            c0 = CLI[i, j]
            c1 = CLI[i, j - 1]
            push!(f_int, (p0, p1, c0, c1))
        end
    end
    # J = 1
    if m.topo == :C
        for i in 1:(m.ncon - 1)
            p0 = NLI[i, 1]
            p1 = NLI[i + 1, 1]
            c0 = CLI[i, 1]
            c1 = CLI[NI - i, 1]
            push!(f_int, (p0, p1, c0, c1))
        end
        for i in (m.ncon):(NI - m.ncon)
            p0 = NLI[i, 1]
            p1 = NLI[i + 1, 1]
            c0 = CLI[i, 1]
            push!(f_wall, (p0, p1, c0, 0))
        end
    else
        for i in 1:(NI - 1)
            p0 = NLI[i, 1]
            p1 = NLI[i + 1, 1]
            c0 = CLI[i, 1]
            push!(f_wall, (p0, p1, c0, 0))
        end
    end
    # J = NJ
    for i in 1:(NI - 1)
        p0 = NLI[i + 1, NJ] # reversed at J=NJ
        p1 = NLI[i, NJ]
        c0 = CLI[i, NJ - 1]
        c1 = 0
        push!(f_far, (p0, p1, c0, c1))
    end
    # I = 1
    for j in 1:(NJ - 1)
        p0 = NLI[1, j + 1] # reversed
        p1 = NLI[1, j]
        c0 = CLI[1, j]
        if m.topo == :O
            c1 = CLI[NI - 1, j]
            push!(f_int, (p0, p1, c0, c1))
        else
            push!(f_min, (p0, p1, c0, 0))
        end
    end
    # I = NI
    if m.topo != :O
        for j in 1:(NJ - 1)
            p0 = NLI[NI, j]
            p1 = NLI[NI, j + 1]
            c0 = CLI[NI - 1, j]
            push!(f_max, (p0, p1, c0, 0))
        end
    end
    NFINT = length(f_int)
    I_INT = (1, NFINT)
    NFWAL = length(f_wall)
    I_WAL = (I_INT[2] + 1, I_INT[2] + NFWAL)
    NFFAR = length(f_far)
    I_FAR = (I_WAL[2] + 1, I_WAL[2] + NFFAR)
    NFMIN = length(f_min)
    I_MIN = (I_FAR[2] + 1, I_FAR[2] + NFMIN)
    NFAMX = length(f_max)
    I_MAX = (I_MIN[2] + 1, I_MIN[2] + NFAMX)
    @assert I_MAX[2] == NF
    # informations
    _fluent_comment(io, "Exported from AirfoilMeshers.jl")
    _fluent_comment(io, "Dim = 2")
    println(io, "(2 2)")
    _fluent_comment(io, "# of nodes: $(NV)")
    println(io, "(10 (0 1 $(string(NV, base=16)) 0 2))")
    _fluent_comment(io, "# of cells: $(NC)")
    println(io, "(12 (0 1 $(string(NC, base=16)) 0))")
    _fluent_comment(io, "# of faces: $(NF)")
    println(io, "(13 (0 1 $(string(NF, base=16)) 0))")
    # Zone 1: Nodes
    izone = 1
    _fluent_comment(io, "Zone $izone: nodes [$NV]")
    println(io, "(10 ($izone 1 $(string(NV, base=16)) 1 2)(")
    for li in 1:NV
        p = m.grid[CartesianIndex(NCI[li])]
        @printf io "%24.16E %24.16E\n" p[1] p[2]
    end
    println(io, "))")
    # Zone 2: Cells
    izone += 1
    _fluent_comment(io, "Zone $izone: cells [$NC]")
    println(io, "(12 ($izone 1 $(string(NC, base=16)) 1 3))")
    println(io, "(45 ($izone fluid flowfield)())")
    # Zone 3: Interior
    izone += 1
    _fluent_comment(io, "Zone $izone: interior faces $(I_INT[1]):$(I_INT[2])")
    println(
        io, "(13 ($izone $(string(I_INT[1], base=16)) $(string(I_INT[2], base=16)) 2 2)("
    )
    for f in f_int
        println(io, join(string.(f, base=16), ' '))
    end
    println(io, "))")
    println(io, "(45 ($izone interior interior-flow)())")
    # Zone 4: Wall
    izone += 1
    _fluent_comment(io, "Zone $izone: wall faces $(I_WAL[1]):$(I_WAL[2])")
    println(
        io, "(13 ($izone $(string(I_WAL[1], base=16)) $(string(I_WAL[2], base=16)) 3 2)("
    )
    for f in f_wall
        println(io, join(string.(f, base=16), ' '))
    end
    println(io, "))")
    println(io, "(45 ($izone wall wall)())")
    # Zone 5: Far
    izone += 1
    _fluent_comment(io, "Zone $izone: far faces $(I_FAR[1]):$(I_FAR[2])")
    println(
        io, "(13 ($izone $(string(I_FAR[1], base=16)) $(string(I_FAR[2], base=16)) a 2)("
    )
    for f in f_far
        println(io, join(string.(f, base=16), ' '))
    end
    println(io, "))")
    println(io, "(45 ($izone velocity-inlet far)())")
    if NFMIN > 0
        # Zone 6: IMin
        izone += 1
        _fluent_comment(io, "Zone $izone: imin faces $(I_MIN[1]):$(I_MIN[2])")
        println(
            io,
            "(13 ($izone $(string(I_MIN[1], base=16)) $(string(I_MIN[2], base=16)) 5 2)(",
        )
        for f in f_min
            println(io, join(string.(f, base=16), ' '))
        end
        println(io, "))")
        println(io, "(45 ($izone pressure-outlet imin)())")
        # Zone 7: IMax
        izone += 1
        _fluent_comment(io, "Zone $izone: imax faces $(I_MAX[1]):$(I_MAX[2])")
        println(
            io,
            "(13 ($izone $(string(I_MAX[1], base=16)) $(string(I_MAX[2], base=16)) 5 2)(",
        )
        for f in f_max
            println(io, join(string.(f, base=16), ' '))
        end
        println(io, "))")
        println(io, "(45 ($izone pressure-outlet imax)())")
    end

    return nothing
end

function _fluent_comment(io::IO, msg::AbstractString, tag::Int=0)
    return println(io, "($tag \"$msg\")")
end


text:
### Remove the long-offset receiver traces, only use +-3km for imaging
## author: Ziyi Yin (ziyi.yin@gatech.edu)

using DrWatson
@quickactivate "CompassTimeLapseCCS"

using PyPlot, JLD2, JUDI

JLD2.@load "../data/d_lin.jld2" d_lin q
JLD2.@load "../model/OBN.jld2"

d_lin_cut = Array{judiVector{Float32,Array{Float32,2}},1}(undef, 5)

for i = 1:5
    xloc_cut = Array{Array{Float32,1}}(undef, d_lin[i].nsrc)
    zloc_cut = Array{Array{Float32,1}}(undef, d_lin[i].nsrc)
    data_cut = Array{Array{Float32,2}}(undef, d_lin[i].nsrc)
    for j = 1:d_lin[i].nsrc
        obn = OBN_loc[j]
        xrec = d_lin[i][j].geometry.xloc[1]
        inds = findall(abs.(xrec.-obn).<=3000f0)
        data_cut[j] = d_lin[i].data[j][:,inds]
        xloc_cut[j] = d_lin[i].geometry.xloc[j][inds]
        zloc_cut[j] = d_lin[i].geometry.zloc[j][inds]
    end
    geometry_cut = Geometry(xloc_cut, d_lin[i].geometry.yloc, zloc_cut; dt=d_lin[i].geometry.dt[1], t=d_lin[i].geometry.t[1])
    d_lin_cut[i] = judiVector(geometry_cut,data_cut)
end

JLD2.@save "../data/d_lin_cut.jld2" d_lin_cut q


text:
abstract type AbstractNode end
abstract type AbstractBranch end

push!(branches::Array{AbstractBranch,1}, b::AbstractBranch) = Base.pushfirst!(branches, b)
push!(branches::Array{AbstractBranch,1}, b::Nothing) = branches

macro abstract_node_fields()
    return esc(quote
        id::Int
        parent::Union{Nothing,AbstractNode}
        depth::Int
        branch::Union{Nothing,T}
        solution_status::MOI.TerminationStatusCode
        bound::Real
        solution::Dict{Any,Real}
    end)
end

# empty branching function
branch(node::AbstractNode)::Vector{AbstractNode} = []

# empty bounding function
function bound!(node::AbstractNode) end

# empty heuristics function
function heuristics!(node::AbstractNode) end

# apply changes (branched information) from ancestor to node
function apply_changes!(node::AbstractNode, ancestor::AbstractNode) end

apply_changes!(node::AbstractNode) = apply_changes!(node, node)


text:
#using FemAdjoint
using Test, JLD2, UnicodePlots
using LinearAlgebra, ForwardDiff, Arpack
using SparseArrays, SparseDiffTools, Symbolics
using FemAdjoint

# load mesh and data for tests
meshfile = "$(@__DIR__)/data/diskmesh.jld2"
JLD2.@load(meshfile, p, t)
np, nt = size(p, 1), size(t, 1)
ε = 1e-5
u, v = 2 * rand(np) .- 1, 2 * rand(np) .- 1
x, dp = p, 2 * rand(size(p)...) .- 1
xp, xm = p .+ ε * dp, p .- ε * dp
SK, SM = FemAdjoint.assembKM_P12D(p, t)
IK, JK = FemAdjoint.indKM_sparse(t)
IM, JM = FemAdjoint.indKM_sparse(t)
K = sparse(IK, JK, SK, np, np)
M = sparse(IM, JM, SM, np, np)

@testset "mesh data" begin
    @test size(p, 1) == 560
    @test size(t, 1) == 1045
end

@testset "first assembly" begin
    @test length(IK) == length(JK)
    @test length(IK) == length(SK)
    @test length(IM) == length(JM)
    @test length(IM) == length(SM)
end

@testset "test eval function" begin
    u, v = 2 * rand(np) .- 1, 2 * rand(np) .- 1
    @test abs(dot(K * u, v) - FemAdjoint.costKproduv(p, t, u, v)) < 1e-8
end

f = x -> FemAdjoint.costKproduv(x, t, u, v)
fp, fm = f(xp), f(xm)
@time g = ForwardDiff.gradient(f, x)
@time g = ForwardDiff.gradient(f, x)

@testset "test eval gradient dot function" begin
    @test abs(dot(g, dp) - (fp - fm) / (2 * ε)) < 1e-4
end
display(K)
print("\n")
println("directional derivative $(abs(dot(g, dp)))")

#-------------------------------------------------------------------------------
# square mesh / eigenvalue test
meshfile = "$(@__DIR__)/data/squaremesh.jld2"
JLD2.@load(meshfile, p, t, jacK, jacM, ps, ts, jacKs, jacMs)
np, nt = size(p, 1), size(t, 1)
SK, SM = FemAdjoint.assembKM_P12D(p, t)
IK, JK = FemAdjoint.indKM_sparse(t)
IM, JM = FemAdjoint.indKM_sparse(t)
K = sparse(IK, JK, SK, np, np)
M = sparse(IM, JM, SM, np, np)
Ib = unique(FemAdjoint.btri(t)[:])
K += sparse(Ib, Ib, fill(1e8, length(Ib)), np, np) # Dirichlet conditions

λ, U = eigs(K, M, nev = 10, which = :SM, tol = 0., maxiter = 30_000)
@testset "test eigenvalue" begin
    @test abs(λ[1] - pi ^ 2 / 2) < 1e-2
end
# plot
FemAdjoint.plotunicode(U[:, 4])

#-------------------------------------------------------------------------------
#p, t, jac = ps, ts, jacs
y = similar(FemAdjoint.assembKM_P12D(p, t)[1])
SK, SM = FemAdjoint.assembKM_P12D(p, t)
SKi, SMi = similar(SK), similar(SM)
FemAdjoint.assembK_P12D_inplace!(p, t, SKi)
FemAdjoint.assembM_P12D_inplace!(p, t, SMi)
@testset "test in place assemble" begin
    @test maximum(abs.(SK .- SKi)) < 1e-8
    @test maximum(abs.(SM .- SMi)) < 1e-8
end
fK(y, x) = FemAdjoint.assembK_P12D_inplace!(x, t, y)
# to speedup test use a precomputed jac matrix
#@time sparsity_pattern = Symbolics.jacobian_sparsity(fs, y, p[:])
#jacK = Float64.(sparse(sparsity_pattern))
colors = matrix_colors(jacK)
@time J = forwarddiff_color_jacobian!(jacK, fK, p[:], colorvec = colors)
@time J = forwarddiff_color_jacobian!(jacK, fK, p[:], colorvec = colors)
display(J)

#-------------------------------------------------------------------------------
# test adjoint approach for costKproduv cost function
include("test_costKproduv.jl")
#-------------------------------------------------------------------------------
# test adjoint approach for costnormU cost function
include("test_costnormU.jl")
#-------------------------------------------------------------------------------
# test direct adjoint approach for costnormU cost function (definitely simpler!!)
include("test_costnormU_direct.jl")


text:
using Distributions
using Base.Test

# RealInterval
r = RealInterval(1.5, 4.0)
@test minimum(r) == 1.5
@test maximum(r) == 4.0


text:
# this should exist?
"""
    filename(path)

Sring -> String

Extract core name of the file from path

E.g. filename("/some/path/to/file_name.ext") -> "file_name"
"""
filename(path) = splitext(basename(path))[1]

# this should be in PropDicts?
"""
    propdict(json_file)

AbstractString -> PropDict

Construct a PropDict based on given <json_file>.

1) Why can't I name it function PropDict()?
2) Why doesn't this already exist in PropDicts?
I find the PropDicts.read(PropDict, String) format kinda cumbersome
"""
function propdict(json_file::AbstractString)
    PropDicts.read(PropDict, json_file)
end

"""
    remove_negative(value)

Real -> Real

Replace negative values by zeros.

Used to remove a glitch in SSD pulses as a quick fix
    while the glitch is being debugged.
"""
function remove_negative(value::Real)
    value < 0 ? 0 : value
end

"""
    preamp_gain(preamp, noise_model)

GenericPreAmp, NoiseFromData -> Real

Calculate gain of <preamp> based on its max energy and
    the offset observed in data baselines contained in <noise_model>

In principle we should not do this in this simulation, and user has to provide
    the precise parameters of the electronics chain used in producing data
    the baselines from which are contained in <noise_model>.
"""
function preamp_gain(preamp::PreAmp, noise_model::NoiseFromData)
    # average offset in baselines
    offset = mean(mean.(noise_model.baseline_catalog.waveform.value))
    # calculate gain based on offset in data baselines
    (typemax(UInt16) - offset) * germanium_ionization_energy / uconvert(u"eV", preamp.max_e)
end

"""
    preamp_gain(preamp, ::NoiseFromSim)

GenericPreAmp -> Float64

Do nothing and return original value of gain parameter in <preamp>,
    since when NoiseFromSim model is used, gain is provided by the user
    (or calculated based on user given max energy)
"""
function preamp_gain(preamp::GenericPreAmp, ::NoiseFromSim)
    preamp.gain
end

# """
#     trigger_threshold(trigger, noise_model, preamp)

# Trigger, NoiseFromSim, GenericPreAmp -> Real

# Calculate trigger threshold based on whether
#     the user provided it in the simulation settings
# """
# function trigger_threshold(trigger::Trigger, noise_model::NoiseFromSim, preamp::PreAmp)
#     # if trigger threshold not given, take noise level as reference
#     threshold = trigger.threshold_keV == 0u"keV" ? noise_model.noise_σ * 3 : trigger.threshold_keV
#     uconvert(u"eV", threshold) / germanium_ionization_energy * preamp.gain
# end

"""
    trigger_threshold(trigger, preamp, ::NoiseFromSim)

Trigger, GenericPreAmp -> Real

In NoiseFromSim setting, non-zero trigger threshold in keV MUST be given by the user.
Calculate final threshold based the value contained in <trigger>, and <preamp> gain.

"""
function trigger_threshold(trigger::Trigger, preamp::GenericPreAmp, ::NoiseFromSim)
    uconvert(u"eV", trigger.threshold_keV) / germanium_ionization_energy * preamp.gain
    # threshold = trigger.threshold_keV == 0u"keV" ? noise_model.noise_σ * 3 : trigger.threshold_keV
    # uconvert(u"eV", threshold) / germanium_ionization_energy * preamp.gain
end

"""
    trigger_threshold(trigger, preamp, noise_model)

Trigger, GenericPreAmp, NoiseFromData -> Real

In NoiseFromData setting, if trigger threshold in keV is not provided,
    it is calculated based on noise levels in the baselines contained in
    <noise_model>.
Otherwise the final threshold is calculated based on <preamp> gain.
"""
function trigger_threshold(trigger::Trigger, preamp::GenericPreAmp, noise_model::NoiseFromData)
    trigger.threshold_keV == 0u"keV" ? std(noise_model.baseline_catalog.waveform[1].value) * 3 : uconvert(u"eV", trigger.threshold_keV) / germanium_ionization_energy * preamp.gain
end