text:
function LineAxis(parent::Scene; kwargs...)

    attrs = merge!(Attributes(kwargs), default_attributes(LineAxis))

    decorations = Dict{Symbol, Any}()

    @extract attrs (endpoints, limits, flipped, ticksize, tickwidth,
        tickcolor, tickalign, ticks, tickformat, ticklabelalign, ticklabelrotation, ticksvisible,
        ticklabelspace, ticklabelpad, labelpadding,
        ticklabelsize, ticklabelsvisible, spinewidth, spinecolor, label, labelsize, labelcolor,
        labelfont, ticklabelfont, ticklabelcolor,
        labelvisible, spinevisible, trimspine, flip_vertical_label, reversed,
        minorticksvisible, minortickalign, minorticksize, minortickwidth, minortickcolor, minorticks)

    pos_extents_horizontal = lift(endpoints) do endpoints
        if endpoints[1][2] == endpoints[2][2]
            horizontal = true
            extents = (endpoints[1][1], endpoints[2][1])
            position = endpoints[1][2]
            return (position, extents, horizontal)
        elseif endpoints[1][1] == endpoints[2][1]
            horizontal = false
            extents = (endpoints[1][2], endpoints[2][2])
            position = endpoints[1][1]
            return (position, extents, horizontal)
        else
            error("OldAxis endpoints $(endpoints[1]) and $(endpoints[2]) are neither on a horizontal nor vertical line")
        end
    end

    ticksnode = Observable(Point2f[])
    ticklines = linesegments!(
        parent, ticksnode, linewidth = tickwidth, color = tickcolor, linestyle = nothing,
        visible = ticksvisible, inspectable = false
    )
    decorations[:ticklines] = ticklines
    translate!(ticklines, 0, 0, 10)

    minorticksnode = Observable(Point2f[])
    minorticklines = linesegments!(
        parent, minorticksnode, linewidth = minortickwidth, color = minortickcolor,
        linestyle = nothing, visible = minorticksvisible, inspectable = false
    )
    decorations[:minorticklines] = minorticklines
    realticklabelalign = Observable{Tuple{Symbol, Symbol}}()

    map!(realticklabelalign, ticklabelalign, pos_extents_horizontal, flipped, ticklabelrotation) do al, (pos, ex, hor), fl, rot
        if al !== automatic
            return al
        end
        if rot == 0 || !(rot isa Real)
            if hor
                (:center, fl ? :bottom : :top)
            else
                (fl ? :left : :right, :center)
            end
        elseif rot ≈ pi/2
            if hor
                (fl ? :left : :right, :center)
            else
                (:center, fl ? :top : :bottom)
            end
        elseif rot ≈ -pi/2
            if hor
                (fl ? :right : :left, :center)
            else
                (:center, fl ? :bottom : :top)
            end
        elseif rot > 0
            if hor
                (fl ? :left : :right, fl ? :bottom : :top)
            else
                (fl ? :left : :right, :center)
            end
        elseif rot < 0
            if hor
                (fl ? :right : :left, fl ? :bottom : :top)
            else
                (fl ? :left : :right, :center)
            end
        end
    end

    ticklabelannosnode = Observable(Tuple{AbstractString, Point2f}[])
    ticklabels = nothing
    ticklabel_ideal_space = Observable{Float32}()
    map!(ticklabel_ideal_space, ticklabelannosnode, ticklabelalign, ticklabelrotation, ticklabelfont, ticklabelsvisible) do args...
        maxwidth = if pos_extents_horizontal[][3]
                # height
                ticklabelsvisible[] ? (ticklabels === nothing ? 0f0 : height(Rect2f(boundingbox(ticklabels)))) : 0f0
            else
                # width
                ticklabelsvisible[] ? (ticklabels === nothing ? 0f0 : width(Rect2f(boundingbox(ticklabels)))) : 0f0
        end
        # in case there is no string in the annotations and the boundingbox comes back all NaN
        if !isfinite(maxwidth)
            maxwidth = zero(maxwidth)
        end
        maxwidth
    end

    attrs[:actual_ticklabelspace] = 0f0
    actual_ticklabelspace = attrs[:actual_ticklabelspace]

    onany(ticklabel_ideal_space, ticklabelspace) do idealspace, space
        s = if space == automatic
            idealspace
        else
            space
        end
        if s != actual_ticklabelspace[]
            actual_ticklabelspace[] = s
        end
    end

    tickspace = lift(ticksvisible, ticksize, tickalign) do ticksvisible,
            ticksize, tickalign

        ticksvisible ? max(0f0, ticksize * (1f0 - tickalign)) : 0f0
    end

    labelgap = lift(spinewidth, tickspace, ticklabelsvisible, actual_ticklabelspace,
        ticklabelpad, labelpadding) do spinewidth, tickspace, ticklabelsvisible,
            actual_ticklabelspace, ticklabelpad, labelpadding

        spinewidth + tickspace +
            (ticklabelsvisible ? actual_ticklabelspace + ticklabelpad : 0f0) +
            labelpadding
    end

    labelpos = lift(pos_extents_horizontal, flipped, labelgap) do (position, extents, horizontal), flipped, labelgap

        # fullgap = tickspace[] + labelgap

        middle = extents[1] + 0.5f0 * (extents[2] - extents[1])

        x_or_y = flipped ? position + labelgap : position - labelgap

        if horizontal
            Point2(middle, x_or_y)
        else
            Point2(x_or_y, middle)
        end
    end

    labelalign = lift(pos_extents_horizontal, flipped, flip_vertical_label) do (position, extents, horizontal), flipped, flip_vertical_label
        if horizontal
            (:center, flipped ? :bottom : :top)
        else
            (:center, if flipped
                    flip_vertical_label ? :bottom : :top
                else
                    flip_vertical_label ? :top : :bottom
                end
            )
        end
    end

    labelrotation = lift(pos_extents_horizontal, flip_vertical_label) do (position, extents, horizontal), flip_vertical_label
        if horizontal
            0f0
        else
            if flip_vertical_label
                Float32(-0.5pi)
            else
                Float32(0.5pi)
            end
        end
    end

    labeltext = text!(
        parent, label, textsize = labelsize, color = labelcolor,
        position = labelpos, visible = labelvisible,
        align = labelalign, rotation = labelrotation, font = labelfont,
        markerspace = :data, inspectable = false
    )

    decorations[:labeltext] = labeltext

    tickvalues = Observable(Float32[])

    tickvalues_labels_unfiltered = lift(pos_extents_horizontal, limits, ticks, tickformat, attrs.scale) do (position, extents, horizontal),
            limits, ticks, tickformat, scale
        get_ticks(ticks, scale, tickformat, limits...)
    end

    tickpositions = Observable(Point2f[])
    tickstrings = Observable(AbstractString[])

    onany(tickvalues_labels_unfiltered, reversed) do tickvalues_labels_unfiltered, reversed

        tickvalues_unfiltered, tickstrings_unfiltered = tickvalues_labels_unfiltered

        position, extents_uncorrected, horizontal = pos_extents_horizontal[]

        extents = reversed ? reverse(extents_uncorrected) : extents_uncorrected

        px_o = extents[1]
        px_width = extents[2] - extents[1]

        lim_o = limits[][1]
        lim_w = limits[][2] - limits[][1]

        # if labels are given manually, it's possible that some of them are outside the displayed limits
        # we only check approximately because otherwise because of floating point errors, ticks can be dismissed sometimes
        i_values_within_limits = findall(tickvalues_unfiltered) do tv
            (limits[][1] <= tv || limits[][1] ≈ tv) &&
             (tv <= limits[][2] || tv ≈ limits[][2])
        end

        tickvalues[] = tickvalues_unfiltered[i_values_within_limits]

        scale = attrs.scale[]
        tickvalues_scaled = scale.(tickvalues[])

        tick_fractions = (tickvalues_scaled .- scale(limits[][1])) ./ (scale(limits[][2]) - scale(limits[][1]))

        tick_scenecoords = px_o .+ px_width .* tick_fractions

        tickpos = if horizontal
            [Point(x, position) for x in tick_scenecoords]
        else
            [Point(position, y) for y in tick_scenecoords]
        end

        # now trigger updates
        tickpositions[] = tickpos

        tickstrings[] = tickstrings_unfiltered[i_values_within_limits]
    end

    minortickvalues = Observable(Float32[])
    minortickpositions = Observable(Point2f[])

    onany(tickvalues, minorticks) do tickvalues, minorticks
        minortickvalues[] = get_minor_tickvalues(minorticks, attrs.scale[], tickvalues, limits[]...)
    end

    onany(minortickvalues) do minortickvalues
        position, extents_uncorrected, horizontal = pos_extents_horizontal[]

        extents = reversed[] ? reverse(extents_uncorrected) : extents_uncorrected

        px_o = extents[1]
        px_width = extents[2] - extents[1]

        lim_o = limits[][1]
        lim_w = limits[][2] - limits[][1]

        scale = attrs.scale[]
        tickvalues_scaled = scale.(minortickvalues)

        tick_fractions = (tickvalues_scaled .- scale(limits[][1])) ./ (scale(limits[][2]) - scale(limits[][1]))

        tick_scenecoords = px_o .+ px_width .* tick_fractions

        minortickpositions[] = if horizontal
            [Point(x, position) for x in tick_scenecoords]
        else
            [Point(position, y) for y in tick_scenecoords]
        end
    end

    onany(minortickpositions, minortickalign, minorticksize, spinewidth) do tickpositions,
        tickalign, ticksize, spinewidth

        position, extents, horizontal = pos_extents_horizontal[]

        if horizontal
            tickstarts = [tp + (flipped[] ? -1f0 : 1f0) * Point2f(0f0, tickalign * ticksize - 0.5f0 * spinewidth) for tp in tickpositions]
            tickends = [t + (flipped[] ? -1f0 : 1f0) * Point2f(0f0, -ticksize) for t in tickstarts]
            minorticksnode[] = interleave_vectors(tickstarts, tickends)
        else
            tickstarts = [tp + (flipped[] ? -1f0 : 1f0) * Point2f(tickalign * ticksize - 0.5f0 * spinewidth, 0f0) for tp in tickpositions]
            tickends = [t + (flipped[] ? -1f0 : 1f0) * Point2f(-ticksize, 0f0) for t in tickstarts]
            minorticksnode[] = interleave_vectors(tickstarts, tickends)
        end
    end

    onany(tickstrings, labelgap, flipped) do tickstrings, labelgap, flipped
        # tickspace is always updated before labelgap
        # tickpositions are always updated before tickstrings
        # so we don't need to lift those

        position, extents, horizontal = pos_extents_horizontal[]

        nticks = length(tickvalues[])

        ticklabelgap = spinewidth[] + tickspace[] + ticklabelpad[]

        shift = if horizontal
            Point2f(0f0, flipped ? ticklabelgap : -ticklabelgap)
        else
            Point2f(flipped ? ticklabelgap : -ticklabelgap, 0f0)
        end

        ticklabelpositions = tickpositions[] .+ Ref(shift)
        ticklabelannosnode[] = collect(zip(tickstrings, ticklabelpositions))
    end

    onany(tickpositions, tickalign, ticksize, spinewidth) do tickpositions,
            tickalign, ticksize, spinewidth

        position, extents, horizontal = pos_extents_horizontal[]

        if horizontal
            tickstarts = [tp + (flipped[] ? -1f0 : 1f0) * Point2f(0f0, tickalign * ticksize - 0.5f0 * spinewidth) for tp in tickpositions]
            tickends = [t + (flipped[] ? -1f0 : 1f0) * Point2f(0f0, -ticksize) for t in tickstarts]
            ticksnode[] = interleave_vectors(tickstarts, tickends)
        else
            tickstarts = [tp + (flipped[] ? -1f0 : 1f0) * Point2f(tickalign * ticksize - 0.5f0 * spinewidth, 0f0) for tp in tickpositions]
            tickends = [t + (flipped[] ? -1f0 : 1f0) * Point2f(-ticksize, 0f0) for t in tickstarts]
            ticksnode[] = interleave_vectors(tickstarts, tickends)
        end
    end

    linepoints = lift(pos_extents_horizontal, flipped, spinewidth, trimspine, tickpositions, tickwidth) do (position, extents, horizontal),
            flipped, sw, trimspine, tickpositions, tickwidth

        if !trimspine || length(tickpositions) < 2
            if horizontal
                y = position
                p1 = Point2f(extents[1] - 0.5sw, y)
                p2 = Point2f(extents[2] + 0.5sw, y)
                [p1, p2]
            else
                x = position
                p1 = Point2f(x, extents[1] - 0.5sw)
                p2 = Point2f(x, extents[2] + 0.5sw)
                [p1, p2]
            end
        else
            [tickpositions[1], tickpositions[end]] .+ [
                (horizontal ? Point2f(-0.5f0 * tickwidth, 0) : Point2f(0, -0.5f0 * tickwidth)),
                (horizontal ? Point2f(0.5f0 * tickwidth, 0) : Point2f(0, 0.5f0 * tickwidth)),
            ]
        end
    end

    decorations[:axisline] = lines!(parent, linepoints, linewidth = spinewidth, visible = spinevisible,
        color = spinecolor, inspectable = false, linestyle = nothing)
    translate!(decorations[:axisline], 0, 0, 20)

    protrusion = lift(ticksvisible, label, labelvisible, labelpadding, labelsize, tickalign, tickspace, ticklabelsvisible, actual_ticklabelspace, ticklabelpad, labelfont, ticklabelfont) do ticksvisible,
            label, labelvisible, labelpadding, labelsize, tickalign, tickspace, ticklabelsvisible,
            actual_ticklabelspace, ticklabelpad, labelfont, ticklabelfont
        position, extents, horizontal = pos_extents_horizontal[]

        label_is_empty = iswhitespace(label) || isempty(label)

        real_labelsize = if label_is_empty
            0f0
        else
            horizontal ? boundingbox(labeltext).widths[2] : boundingbox(labeltext).widths[1]
        end

        labelspace = (labelvisible && !label_is_empty) ? real_labelsize + labelpadding : 0f0
        # tickspace = ticksvisible ? max(0f0, xticksize * (1f0 - xtickalign)) : 0f0
        tickspace = (ticksvisible && !isempty(ticklabelannosnode[])) ? tickspace : 0f0

        ticklabelgap = (ticklabelsvisible && actual_ticklabelspace > 0) ? actual_ticklabelspace + ticklabelpad : 0f0

        together = tickspace + ticklabelgap + labelspace
    end

    # trigger whole pipeline once to fill tickpositions and tickstrings
    # etc to avoid empty ticks bug #69
    limits[] = limits[]

    # in order to dispatch to the correct text recipe later (normal text, latex, etc.)
    # we need to have the ticklabelannosnode populated once before adding the annotations
    ticklabels = text!(
        parent,
        ticklabelannosnode,
        align = realticklabelalign,
        rotation = ticklabelrotation,
        textsize = ticklabelsize,
        font = ticklabelfont,
        color = ticklabelcolor,
        visible = ticklabelsvisible,
        markerspace = :data,
        inspectable = false)
    decorations[:ticklabels] = ticklabels

    # HACKY: the ticklabels in the string need to be updated
    # before other stuff is triggered by them, which accesses the
    # ticklabel boundingbox (which needs to be updated already)
    # so we move the new listener from text! to the front
    pushfirst!(
        ticklabelannosnode.listeners,
        pop!(ticklabelannosnode.listeners))

    # trigger calculation of ticklabel width once, now that it's not nothing anymore
    notify(ticklabelsvisible)

    LineAxis(parent, protrusion, attrs, decorations, tickpositions, tickvalues, tickstrings, minortickpositions, minortickvalues)
end

function tight_ticklabel_spacing!(la::LineAxis)

    horizontal = if la.attributes.endpoints[][1][2] == la.attributes.endpoints[][2][2]
        true
    elseif la.attributes.endpoints[][1][1] == la.attributes.endpoints[][2][1]
        false
    else
        error("endpoints not on a horizontal or vertical line")
    end

    tls = la.elements[:ticklabels]
    maxwidth = if horizontal
            # height
            tls.visible[] ? height(Rect2f(boundingbox(tls))) : 0f0
        else
            # width
            tls.visible[] ? width(Rect2f(boundingbox(tls))) : 0f0
    end
    la.attributes.ticklabelspace = maxwidth
end

function iswhitespace(str)
    match(r"^\s*$", str) !== nothing
end

function Base.delete!(la::LineAxis)
    for (_, d) in la.elements
        if d isa AbstractPlot
            delete!(d.parent, d)
        else
            delete!(d)
        end
    end
end

"""
    get_ticks(ticks, scale, formatter, vmin, vmax)

Base function that calls `get_tickvalues(ticks, vmin, max)` and
`get_ticklabels(formatter, tickvalues)` and returns a tuple
`(tickvalues, ticklabels)`.
For custom ticks / formatter combinations, this method can be overloaded
directly, or both `get_tickvalues` and `get_ticklabels` separately.
"""
function get_ticks(ticks, scale, formatter, vmin, vmax)
    tickvalues = get_tickvalues(ticks, scale, vmin, vmax)
    ticklabels = get_ticklabels(formatter, tickvalues)
    return tickvalues, ticklabels
end

# automatic with identity scaling uses WilkinsonTicks by default
get_tickvalues(::Automatic, ::typeof(identity), vmin, vmax) = get_tickvalues(WilkinsonTicks(5, k_min = 3), vmin, vmax)

# fall back to identity if not overloaded scale function is used with automatic
get_tickvalues(::Automatic, F, vmin, vmax) = get_tickvalues(automatic, identity, vmin, vmax)

# fall back to non-scale aware behavior if no special version is overloaded
get_tickvalues(ticks, scale, vmin, vmax) = get_tickvalues(ticks, vmin, vmax)

function get_ticks(ticks_and_labels::Tuple{Any, Any}, any_scale, ::Automatic, vmin, vmax)
    n1 = length(ticks_and_labels[1])
    n2 = length(ticks_and_labels[2])
    if n1 != n2
        error("There are $n1 tick values in $(ticks_and_labels[1]) but $n2 tick labels in $(ticks_and_labels[2]).")
    end
    ticks_and_labels
end

function get_ticks(tickfunction::Function, any_scale, formatter, vmin, vmax)
    result = tickfunction(vmin, vmax)
    if result isa Tuple{Any, Any}
        tickvalues, ticklabels = result
    else
        tickvalues = result
        ticklabels = get_ticklabels(formatter, tickvalues)
    end
    return tickvalues, ticklabels
end

_logbase(::typeof(log10)) = "10"
_logbase(::typeof(log2)) = "2"
_logbase(::typeof(log)) = "e"

function get_ticks(::Automatic, scale::Union{typeof(log10), typeof(log2), typeof(log)},
        any_formatter, vmin, vmax)
    get_ticks(LogTicks(WilkinsonTicks(5, k_min = 3)), scale, any_formatter, vmin, vmax)
end

# log ticks just use the normal pipeline but with log'd limits, then transform the labels
function get_ticks(l::LogTicks, scale::Union{typeof(log10), typeof(log2), typeof(log)}, ::Automatic, vmin, vmax)
    ticks_scaled = get_tickvalues(l.linear_ticks, identity, scale(vmin), scale(vmax))

    ticks = Makie.inverse_transform(scale).(ticks_scaled)

    labels_scaled = get_ticklabels(
        # avoid unicode superscripts in ticks, as the ticks are converted
        # to superscripts in the next step
        xs -> Showoff.showoff(xs, :plain),
        ticks_scaled
    )
    labels = _logbase(scale) .* Makie.UnicodeFun.to_superscript.(labels_scaled)

    (ticks, labels)
end

# function get_ticks(::Automatic, scale::typeof(Makie.logit), any_formatter, vmin, vmax)
#     get_ticks(LogitTicks(WilkinsonTicks(5, k_min = 3)), scale, any_formatter, vmin, vmax)
# end

logit_10(x) = Makie.logit(x) / log(10)
expit_10(x) = Makie.logistic(log(10) * x)

# function get_ticks(l::LogitTicks, scale::typeof(Makie.logit), ::Automatic, vmin, vmax)

#     ticks_scaled = get_tickvalues(l.linear_ticks, identity, logit_10(vmin), logit_10(vmax))

#     ticks = expit_10.(ticks_scaled)

#     base_labels = get_ticklabels(automatic, ticks_scaled)

#     labels = map(ticks_scaled, base_labels) do t, bl
#         if t == 0
#             "¹/₂"
#         elseif t < 0
#             "10" * Makie.UnicodeFun.to_superscript(bl)
#         else
#             "1-10" * Makie.UnicodeFun.to_superscript("-" * bl)
#         end
#     end

#     (ticks, labels)
# end

"""
    get_tickvalues(lt::LinearTicks, vmin, vmax)

Runs a common tick finding algorithm to as many ticks as requested by the
`LinearTicks` instance.
"""
get_tickvalues(lt::LinearTicks, vmin, vmax) = locateticks(vmin, vmax, lt.n_ideal)

"""
    get_tickvalues(tickvalues, vmin, vmax)

Convert tickvalues to a float array by default.
"""
get_tickvalues(tickvalues, vmin, vmax) = convert(Vector{Float64}, tickvalues)

# function get_tickvalues(l::LogitTicks, vmin, vmax)
#     ticks_scaled = get_tickvalues(l.linear_ticks, identity, logit_10(vmin), logit_10(vmax))
#     expit_10.(ticks_scaled)
# end

function get_tickvalues(l::LogTicks, scale, vmin, vmax)
    ticks_scaled = get_tickvalues(l.linear_ticks, scale(vmin), scale(vmax))
    Makie.inverse_transform(scale).(ticks_scaled)
end

"""
    get_ticklabels(::Automatic, values)

Gets tick labels by applying `Showoff.showoff` to `values`.
"""
get_ticklabels(::Automatic, values) = Showoff.showoff(values)

"""
    get_ticklabels(formatfunction::Function, values)

Gets tick labels by applying `formatfunction` to `values`.
"""
get_ticklabels(formatfunction::Function, values) = formatfunction(values)

"""
    get_ticklabels(formatstring::AbstractString, values)

Gets tick labels by formatting each value in `values` according to a `Formatting.format` format string.
"""
get_ticklabels(formatstring::AbstractString, values) = [Formatting.format(formatstring, v) for v in values]

function get_ticks(m::MultiplesTicks, any_scale, ::Automatic, vmin, vmax)
    dvmin = vmin / m.multiple
    dvmax = vmax / m.multiple
    multiples = MakieLayout.get_tickvalues(LinearTicks(m.n_ideal), dvmin, dvmax)

    multiples .* m.multiple, Showoff.showoff(multiples) .* m.suffix
end

function get_minor_tickvalues(i::IntervalsBetween, scale, tickvalues, vmin, vmax)
    vals = Float64[]
    length(tickvalues) < 2 && return vals
    n = i.n

    if i.mirror
        firstinterval = tickvalues[2] - tickvalues[1]
        stepsize = firstinterval / n
        v = tickvalues[1] - stepsize
        prepend!(vals, v:-stepsize:vmin)
    end

    for (lo, hi) in zip(@view(tickvalues[1:end-1]), @view(tickvalues[2:end]))
        interval = hi - lo
        stepsize = interval / n
        v = lo
        for i in 1:n-1
            v += stepsize
            push!(vals, v)
        end
    end

    if i.mirror
        lastinterval = tickvalues[end] - tickvalues[end-1]
        stepsize = lastinterval / n
        v = tickvalues[end] + stepsize
        append!(vals, v:stepsize:vmax)
    end

    vals
end

# for log scales, we need to step in log steps at the edges
function get_minor_tickvalues(i::IntervalsBetween, scale::Union{typeof(log), typeof(log2), typeof(log10)}, tickvalues, vmin, vmax)

    vals = Float64[]
    length(tickvalues) < 2 && return vals
    n = i.n

    invscale = Makie.inverse_transform(scale)

    if i.mirror
        firstinterval_scaled = scale(tickvalues[2]) - scale(tickvalues[1])
        stepsize = firstinterval_scaled / n
        prevtick = invscale(scale(tickvalues[1]) - firstinterval_scaled)
        stepsize = (tickvalues[1] - prevtick) / n
        v = tickvalues[1] - stepsize
        prepend!(vals, v:-stepsize:vmin)
    end

    for (lo, hi) in zip(@view(tickvalues[1:end-1]), @view(tickvalues[2:end]))
        interval = hi - lo
        stepsize = interval / n
        v = lo
        for i in 1:n-1
            v += stepsize
            push!(vals, v)
        end
    end

    if i.mirror
        lastinterval_scaled = scale(tickvalues[end]) - scale(tickvalues[end-1])
        nexttick = invscale(scale(tickvalues[end]) + lastinterval_scaled)
        stepsize = (nexttick - tickvalues[end]) / n
        v = tickvalues[end] + stepsize
        append!(vals, v:stepsize:vmax)
    end

    vals
end

function get_minor_tickvalues(v::AbstractVector{<:Real}, _, _, _, _)
    Float32.(v)
end


text:
using FixedPoint, Test
@testset "Vector Equation" begin
    Ts = LinRange(0.01, 2.0, 500)
    βs = 1 ./ Ts
    mag = afps(
        x -> tanh.(βs .* x),
        zero(βs) .+ 1,
        grad_norm = x -> maximum(abs, x),
        iters = 5000,
    )
    @test mag.error < 1e-4
end

@testset "Scalar Equation" begin
    s = afps(x -> 2 - x^2 + x, 1.3)
    @test s.x ≈ √2
end

@testset "Inplace Vector Equation" begin
    Ts = LinRange(0.01, 2.0, 500)
    βs = 1 ./ Ts
    function f!(out,x)
        @. out = tanh(βs * x)
    end
    x  = zero(βs) .+ 1
    mag = afps!(
        f!,
        x,
        grad_norm = x -> maximum(abs, x),
        iters = 5000,
    )
    @test mag.error < 1e-4
    @test maximum(abs, x .- tanh.(βs .* x)) < 1e-4
end

@testset "Inplace Equation" begin
    x = [1.3]
    function f!(out, x)
        out[1] = 2 - x[1]^2 + x[1]
    end
    s = afps!(f!, x)
    @test s.x[1] ≈ √2
end

ext: jl


text: include("myfile.jl")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

##############
# test_hybsv #
##############

function test_hybsv(elty)
    A = rand(elty,m,m)
    A = triu(A)
    x = rand(elty,m)
    alpha = rand(elty)
    beta = rand(elty)
    d_x = CudaArray(x)
    d_A = CudaSparseMatrixCSR(sparse(A))
    d_A = CUSPARSE.switch2hyb(d_A)
    d_y = CUSPARSE.sv('N','T','U',alpha,d_A,d_x,'O')
    h_y = to_host(d_y)
    y = A\(alpha * x)
    @test_approx_eq(y,h_y)

    d_y = UpperTriangular(d_A) \ d_x
    h_y = to_host(d_y)
    @test h_y ≈ A\x

    d_x = CudaArray(rand(elty,n))
    info = CUSPARSE.sv_analysis('N','T','U',d_A,'O')
    @test_throws(DimensionMismatch, CUSPARSE.sv_solve('N','U',alpha,d_A,d_x,info,'O'))
    A = sparse(rand(elty,m,n))
    d_A = CudaSparseMatrixCSR(A)
    d_A = CUSPARSE.switch2hyb(d_A)
    @test_throws(DimensionMismatch, CUSPARSE.sv_analysis('T','T','U',d_A,'O'))
    CUSPARSE.cusparseDestroySolveAnalysisInfo(info)
end

types = [Float32,Float64,Complex64,Complex128]
for elty in types
    tic()
    test_hybsv(elty)
    println("hybsv took ", toq(), " for ", elty)
end


text:
# Defining a function
life_universe_everything(args...) = 42

max_stars_repo_path: 2 Tips and tricks/3_share-and-instantiate-projects/my-project/somedefinitions.jl

max_stars_repo_name: crstnbr/julia-workshop


text:
# Test objective functions
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, Test, LinearAlgebra

###################################################################################################
# MSE loss

Y = randn(Float32, 28, 28, 3, 1)
X = randn(Float32, 28, 28, 3, 1)
X0 = randn(Float32, 28, 28, 3, 1)
dX = X - X0

function loss_mse(X, Y)
    f = mse(X, Y)
    g = ∇mse(X, Y)
    return f, g
end

f0, g = loss_mse(X0, Y)
h = .1f0
maxiter = 6

err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test mse loss\n")
for j=1:maxiter
    f = loss_mse(X0 + h*dX, Y)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dX, g))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

###################################################################################################
# (Negative) Log-likelihood

X = randn(Float32, 28, 28, 3, 1)
X0 = randn(Float32, 28, 28, 3, 1)
dX = X - X0

function loss_mle(X)
    f = log_likelihood(X; μ=1f0, σ=2f0)
    g = ∇log_likelihood(X; μ=1f0, σ=2f0)
    return f, g
end

f0, g = loss_mle(X0)
h = .1f0
maxiter = 6

err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)

print("\nGradient test log likelihood\n")
for j=1:maxiter
    f = loss_mle(X0 + h*dX)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(dX, g))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)


text:
module SphericalBesselFunctions

using LinearAlgebra
using Formatting

include("continued_fractions.jl")
include("bessels.jl")
include("coulombs.jl")

end # module


text:
#=
Scope information for `EXRR`.

NOTE:
- StaticLint.jl adds various scope-related information for `EXPR.meta.scope` field,
    but we currently just need to know whether `EXPR` introduces a scope or not.
    So let's just make a function to check that, and don't fill `EXPR` with scope information.
- Adapted from https://github.com/julia-vscode/StaticLint.jl/blob/cd8935a138caf18385c46db977b52c5ac9e90809/src/scope.jl
=#

function hasscope(x::EXPR)
    t = typof(x)

    # NOTE: added conditions below when adapted
    if t === CSTParser.TupleH && (p = parentof(x)) !== nothing && !hasscope(p)
        return true
    elseif iswhereclause(x)
        return true
    elseif t === CSTParser.MacroCall
        return true
    elseif t === CSTParser.Quote
        return true
    # NOTE: end

    elseif t === CSTParser.BinaryOpCall
        k = kindof(x.args[2])
        if k === Tokens.EQ && CSTParser.is_func_call(x.args[1])
            return true
        elseif k === Tokens.EQ && typof(x.args[1]) === CSTParser.Curly
            return true
        elseif k === Tokens.ANON_FUNC
            return true
        else
            return false
        end
    # # NOTE: commented out when adapted
    # elseif t === CSTParser.WhereOpCall
    #     # unless in func def signature
    #     return !_in_func_def(x)
    elseif t === CSTParser.FunctionDef ||
           t === CSTParser.Macro ||
           t === CSTParser.For ||
           t === CSTParser.While ||
           t === CSTParser.Let ||
           t === CSTParser.Generator || # and Flatten?
           t === CSTParser.Try ||
           t === CSTParser.Do ||
           t === CSTParser.ModuleH ||
           t === CSTParser.BareModule ||
           t === CSTParser.Abstract ||
           t === CSTParser.Primitive ||
           t === CSTParser.Mutable ||
           t === CSTParser.Struct
        return true
    end

    return false
end

# # NOTE: commented out when adapted
# # only called in WhereOpCall
# function _in_func_def(x::EXPR)
#     # check 1st arg contains a call (or op call)
#     ex = x.args[1]
#     while true
#         if typof(ex) === CSTParser.WhereOpCall ||
#            (
#             typof(ex) === CSTParser.BinaryOpCall &&
#             kindof(ex.args[2]) === CSTParser.Tokens.DECLARATION
#            )
#             ex = ex.args[1]
#         elseif typof(ex) === CSTParser.Call ||
#                (
#                 typof(ex) === CSTParser.BinaryOpCall &&
#                 !(kindof(ex.args[2]) === CSTParser.Tokens.DOT)
#                ) ||
#                typof(ex) == CSTParser.UnaryOpCall #&& kindof(ex.args[1]) == CSTParser.Tokens.MINUS
#             break
#         else
#             return false
#         end
#     end
#     # check parent is func def
#     ex = x
#     while true
#         if !(parentof(ex) isa EXPR)
#             return false
#         elseif typof(parentof(ex)) === CSTParser.WhereOpCall ||
#                typof(parentof(ex)) === CSTParser.InvisBrackets
#             ex = parentof(ex)
#         elseif typof(parentof(ex)) === CSTParser.FunctionDef ||
#                (
#                 typof(parentof(ex)) === CSTParser.BinaryOpCall &&
#                 kindof(parentof(ex).args[2]) === CSTParser.Tokens.EQ
#                )
#             return true
#         else
#             return false
#         end
#     end
#     return false
# end


text:
ranges, my_ticket, all_tickets = split(read("day16/input.txt", String), "\n\n")

ranges = split(ranges, "\n")
my_ticket = split(my_ticket, "\n")[2]
all_tickets = split(all_tickets, "\n")[2:end]

valid_nums = Set{Int}()

for r in ranges
    s1, e1, s2, e2 = parse.(Int, match(r"(\d+)-(\d+) or (\d+)-(\d+)", r).captures)
    push!(valid_nums, s1:e1..., s2:e2...)
end

out = 0
for t in all_tickets
    t_nums = parse.(Int, split(t, ","))
    global out += sum(filter(n -> !(n in valid_nums), t_nums))
end

println(out)


text:
using IJuliaTimeMachine
using Test

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

max_stars_repo_name: danspielman/IJuliaTimeMachine

max_issues_repo_name: danspielman/IJuliaTimeMachine

max_forks_repo_name: danspielman/IJuliaTimeMachine


text:
"""
    test_bqgql_rossbywave(; kwargs...)

Evolvesa a Rossby wave and compares with the analytic solution.
"""
function test_bqgql_rossbywave(stepper, dt, nsteps, dev::Device=CPU())
  nx = 64
   β = 2.0
  Lx = 2π
   μ = 0.0
   ν = 0.0
   T = Float64

  # the following if statement is called so that all the cases of
  # Problem() function are tested
  if stepper=="ForwardEuler"
    eta = zeros(dev, T, (nx, nx))
  else
    eta(x, y) = 0*x
  end

  prob = BarotropicQGQL.Problem(dev; nx=nx, Lx=Lx, eta=eta, β=β, μ=μ, ν=ν, stepper=stepper, dt=dt)
  sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

  x, y = gridpoints(grid)

  # the Rossby wave initial condition
   ampl = 1e-2
  kwave, lwave = 3 * 2π / grid.Lx, 2 * 2π / grid.Ly
      ω = - params.β * kwave / (kwave^2 + lwave^2)
     ζ0 = @. ampl * cos(kwave * x) * cos(lwave * y)
    ζ0h = rfft(ζ0)

  BarotropicQGQL.set_zeta!(prob, ζ0)

  stepforward!(prob, nsteps)
  dealias!(sol, grid)
  BarotropicQGQL.updatevars!(prob)

  ζ_theory = @. ampl * cos(kwave * (x - ω / kwave * clock.t)) * cos(lwave * y)

  return isapprox(ζ_theory, vars.zeta, rtol = grid.nx * grid.ny * nsteps * 1e-12)
end

"""
    test_stochasticforcingbudgets(dev; kwargs...)

Tests if the energy budget is closed for BarotropicQG problem with stochastic forcing.
"""
function test_bqgql_stochasticforcingbudgets(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1, T=Float64)
  n, L  = 256, 2π
  ν, nν = 1e-7, 2
  μ = 1e-1
  dt, tf = 0.005, 0.1/μ
  nt = round(Int, tf/dt)

  gr = TwoDGrid(dev, n, L)
  x, y = gridpoints(gr)

  # Forcing
  kf, dkf = 12.0, 2.0
  ε = 0.1

  CUDA.@allowscalar Kr = ArrayType(dev)([ gr.kr[i] for i=1:gr.nkr, j=1:gr.nl])

  forcingcovariancespectrum = zeros(dev, T, (gr.nkr, gr.nl))
  @. forcingcovariancespectrum = exp(-(sqrt(gr.Krsq) - kf)^2 / (2 * dkf^2))
  CUDA.@allowscalar @. forcingcovariancespectrum[gr.Krsq .< 2^2] = 0
  CUDA.@allowscalar @. forcingcovariancespectrum[gr.Krsq .> 20^2] = 0
  CUDA.@allowscalar @. forcingcovariancespectrum[Kr .< 2π/L] = 0
  ε0 = parsevalsum(forcingcovariancespectrum .* gr.invKrsq / 2, gr) / (gr.Lx * gr.Ly)
  forcingcovariancespectrum .= ε / ε0 * forcingcovariancespectrum

  Random.seed!(1234)

  function calcF!(F, sol, t, clock, vars, params, grid)
    eta = ArrayType(dev)(exp.(2π * im * rand(T, size(sol))) / sqrt(clock.dt))
    CUDA.@allowscalar eta[1, 1] = 0
    @. F = eta * sqrt(forcingcovariancespectrum)
    return nothing
  end

  prob = BarotropicQGQL.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt,
   stepper="RK4", calcF=calcF!, stochastic=true)

  BarotropicQGQL.set_zeta!(prob, 0*x)
  E = Diagnostic(BarotropicQGQL.energy,      prob, nsteps=nt)
  D = Diagnostic(BarotropicQGQL.dissipation, prob, nsteps=nt)
  R = Diagnostic(BarotropicQGQL.drag,        prob, nsteps=nt)
  W = Diagnostic(BarotropicQGQL.work,        prob, nsteps=nt)
  diags = [E, D, W, R]

  stepforward!(prob, diags, round(Int, nt))

  BarotropicQGQL.updatevars!(prob)

  dEdt_numerical = (E[2:E.i] - E[1:E.i-1]) / prob.clock.dt

  # If the Ito interpretation was used for the work
  # then we need to add the drift term
  # dEdt_computed = W[2:E.i] + ε - D[1:E.i-1] - R[1:E.i-1]      # Ito
  dEdt_computed = W[2:E.i] - D[1:E.i-1] - R[1:E.i-1]        # Stratonovich

  return isapprox(dEdt_numerical, dEdt_computed, rtol=1e-3)
end

"""
    test_bqgql_deterministicforcingbudgets(dev ; kwargs...)

Tests if the energy budget is closed for BarotropicQGQL problem with deterministic forcing.
"""
function test_bqgql_deterministicforcingbudgets(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1)
  n, L  = 256, 2π
  ν, nν = 1e-7, 2
  μ = 1e-1
  dt, tf = 0.005, 0.1/μ
  nt = round(Int, tf/dt)

  gr = TwoDGrid(dev, n, L)
  x, y = gridpoints(gr)
  k₀, l₀ = 2π/gr.Lx, 2π/gr.Ly

  # Forcing = 0.01cos(4x)cos(5y)cos(2t)
  f = @. 0.01 * cos(4k₀*x) * cos(5l₀*y)
  fh = rfft(f)

  function calcF!(Fh, sol, t, clock, vars, params, grid)
    @. Fh = fh*cos(2t)
    return nothing
  end

  prob = BarotropicQGQL.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt,
   stepper="RK4", calcF=calcF!, stochastic=false)

  BarotropicQGQL.set_zeta!(prob, 0*x)

  E = Diagnostic(BarotropicQGQL.energy,      prob, nsteps=nt)
  D = Diagnostic(BarotropicQGQL.dissipation, prob, nsteps=nt)
  R = Diagnostic(BarotropicQGQL.drag,        prob, nsteps=nt)
  W = Diagnostic(BarotropicQGQL.work,        prob, nsteps=nt)
  diags = [E, D, W, R]

  stepforward!(prob, diags, nt)

  BarotropicQGQL.updatevars!(prob)

  dEdt_numerical = (E[3:E.i] - E[1:E.i-2]) / (2 * prob.clock.dt)
  dEdt_computed  = W[2:E.i-1] - D[2:E.i-1] - R[2:E.i-1]

  residual = dEdt_numerical - dEdt_computed

  return isapprox(dEdt_numerical, dEdt_computed, atol=1e-10)
end

"""
    test_bqgql_nonlinearadvection(dt, stepper, dev; kwargs...)

Tests the advection term by timestepping a test problem with timestep `dt` and timestepper
`stepper`. The test problem is derived by picking a solution ζf (with associated
streamfunction ψf) for which the advection term J(ψf, ζf) is non-zero. Next, a
forcing Ff is derived according to Ff = ∂ζf/∂t + J(ψf, ζf) - νΔζf. One solution
to the vorticity equation forced by this Ff is then ζf. (This solution may not
be realized, at least at long times, if it is unstable.)
"""
function test_bqgql_advection(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=1e-2, nν=1, μ=0.0)
  n, L  = 128, 2π
  ν, nν = 1e-2, 1
   μ = 0.0
  tf = 1.0
  nt = round(Int, tf/dt)

  gr = TwoDGrid(dev, n, L)
  x, y = gridpoints(gr)

  ψf = @.    cos(3y) +  sin(2x)*cos(2y) +  2sin(x)*cos(3y)
  qf = @. - 9cos(3y) - 8sin(2x)*cos(2y) - 20sin(x)*cos(3y)

  Ff = @. ν*( 81cos(3y) + 200cos(3y)*sin(x) + 64cos(2y)*sin(2x) ) -
    3sin(3y)*(-16cos(2x)*cos(2y) - 20cos(x)*cos(3y)) -
 27sin(3y)*(2cos(2x)*cos(2y) + 2cos(x)*cos(3y)) + 0*(-8cos(x)*cos(3y)*sin(2x)*sin(2y) +
 24*cos(2x)*cos(2y)*sin(x)*sin(3y))

  Ffh = -rfft(Ff)

  # Forcing
  function calcF!(Fh, sol, t, cl, v, p, g)
    Fh .= Ffh
    return nothing
  end

  prob = BarotropicQGQL.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt, stepper=stepper, calcF=calcF!)

  BarotropicQGQL.set_zeta!(prob, qf)

  stepforward!(prob, round(Int, nt))

  BarotropicQGQL.updatevars!(prob)

  return isapprox(prob.vars.zeta + prob.vars.Zeta, qf, rtol = 1e-13)
end

"""
    test_bqgql_energyenstrophy(dev)

Tests the energy and enstrophy function for a BarotropicQGQL problem.
"""
function test_bqgql_energyenstrophy(dev::Device=CPU())
  nx, Lx  = 64, 2π
  ny, Ly  = 64, 3π
  gr = TwoDGrid(dev, nx, Lx, ny, Ly)
  k₀, l₀ = 2π/Lx, 2π/Ly # fundamental wavenumbers
  x, y = gridpoints(gr)

  energy_calc = 29/9
  enstrophy_calc = 2701/162

    eta = @. cos(10k₀*x)*cos(10l₀*y)
  ψ₀ = @. sin(2k₀*x)*cos(2l₀*y) + 2sin(k₀*x)*cos(3l₀*y)
  ζ₀ = @. -((2k₀)^2+(2l₀)^2)*sin(2k₀*x)*cos(2l₀*y) - (k₀^2+(3l₀)^2)*2sin(k₀*x)*cos(3l₀*y)

  prob = BarotropicQGQL.Problem(dev; nx=nx, Lx=Lx, ny=ny, Ly=Ly, eta=eta, stepper="ForwardEuler")
  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

  BarotropicQGQL.set_zeta!(prob, ζ₀)
  BarotropicQGQL.updatevars!(prob)

  energyζ₀ = BarotropicQGQL.energy(prob)
  enstrophyζ₀ = BarotropicQGQL.enstrophy(prob)

  return isapprox(energyζ₀, energy_calc, rtol=1e-13) && isapprox(enstrophyζ₀, enstrophy_calc, rtol=1e-13) && BarotropicQGQL.addforcing!(prob.timestepper.N, sol, cl.t, cl, v, p, g)==nothing
end

function test_bqgql_problemtype(dev, T)
  prob = BarotropicQGQL.Problem(dev; T=T)

  A = ArrayType(dev)

  return (typeof(prob.sol)<:A{Complex{T}, 2} && typeof(prob.grid.Lx)==T && eltype(prob.grid.x)==T && typeof(prob.vars.u)<:A{T, 2})
end


text:
@testset "GNNGraph" begin

    @testset "Constructor: adjacency matrix" begin
        A = sprand(10, 10, 0.5)
        sA, tA, vA = findnz(A)

        g = GNNGraph(A, graph_type=GRAPH_T)
        s, t = edge_index(g)
        v = get_edge_weight(g)
        @test s == sA
        @test t == tA
        @test v == vA

        g = GNNGraph(Matrix(A), graph_type=GRAPH_T)
        s, t = edge_index(g)
        v = get_edge_weight(g)
        @test s == sA
        @test t == tA
        @test v == vA

        g = GNNGraph([0 0 0
                      0 0 1
                      0 1 0], graph_type=GRAPH_T)
        @test g.num_nodes == 3
        @test g.num_edges == 2

        g = GNNGraph([0 1 0
                      1 0 0
                      0 0 0], graph_type=GRAPH_T)
        @test g.num_nodes == 3
        @test g.num_edges == 2
    end

    @testset "Constructor: integer" begin
        g = GNNGraph(10, graph_type=GRAPH_T)
        @test g.num_nodes == 10
        @test g.num_edges == 0

        g2 = rand_graph(10, 30, graph_type=GRAPH_T)
        G = typeof(g2)
        g = G(10)
        @test g.num_nodes == 10
        @test g.num_edges == 0
    end

    @testset "symmetric graph" begin
        s = [1, 1, 2, 2, 3, 3, 4, 4]
        t = [2, 4, 1, 3, 2, 4, 1, 3]
        adj_mat =  [0  1  0  1
                    1  0  1  0
                    0  1  0  1
                    1  0  1  0]
        adj_list_out =  [[2,4], [1,3], [2,4], [1,3]]
        adj_list_in =  [[2,4], [1,3], [2,4], [1,3]]

        # core functionality
        g = GNNGraph(s, t; graph_type=GRAPH_T)
        if TEST_GPU
            g_gpu = g |> gpu
        end

        @test g.num_edges == 8
        @test g.num_nodes == 4
        @test nv(g) == g.num_nodes
        @test ne(g) == g.num_edges
        @test collect(edges(g)) |> sort == collect(zip(s, t)) |> sort
        @test sort(outneighbors(g, 1)) == [2, 4]
        @test sort(inneighbors(g, 1)) == [2, 4]
        @test is_directed(g) == true
        s1, t1 = sort_edge_index(edge_index(g))
        @test s1 == s
        @test t1 == t
        @test vertices(g) == 1:g.num_nodes

        @test sort.(adjacency_list(g; dir=:in)) == adj_list_in
        @test sort.(adjacency_list(g; dir=:out)) == adj_list_out

        @testset "adjacency_matrix" begin
            @test adjacency_matrix(g) == adj_mat
            @test adjacency_matrix(g; dir=:in) == adj_mat
            @test adjacency_matrix(g; dir=:out) == adj_mat

            if TEST_GPU
                # See https://github.com/JuliaGPU/CUDA.jl/pull/1093
                mat_gpu = adjacency_matrix(g_gpu)
                @test mat_gpu isa ACUMatrix{Int}
                @test Array(mat_gpu) == adj_mat
            end
        end

        @testset "normalized_laplacian" begin
            mat = normalized_laplacian(g)
            if TEST_GPU
                mat_gpu = normalized_laplacian(g_gpu)
                @test mat_gpu isa ACUMatrix{Float32}
                @test Array(mat_gpu) == mat
            end
        end

        @testset "scaled_laplacian" begin
            if TEST_GPU
                @test_broken begin
                    mat = scaled_laplacian(g)
                    mat_gpu = scaled_laplacian(g_gpu)
                    @test mat_gpu isa ACUMatrix{Float32}
                    @test Array(mat_gpu) == mat
                end
            end
        end

        @testset "constructors" begin
            adjacency_matrix(g; dir=:out) == adj_mat
            adjacency_matrix(g; dir=:in) == adj_mat
        end

        if TEST_GPU
            @testset "functor" begin
                s_cpu, t_cpu = edge_index(g)
                s_gpu, t_gpu = edge_index(g_gpu)
                @test s_gpu isa CuVector{Int}
                @test Array(s_gpu) == s_cpu
                @test t_gpu isa CuVector{Int}
                @test Array(t_gpu) == t_cpu
            end
        end
    end

    @testset "asymmetric graph" begin
        s = [1, 2, 3, 4]
        t = [2, 3, 4, 1]
        adj_mat_out =  [0  1  0  0
                        0  0  1  0
                        0  0  0  1
                        1  0  0  0]
        adj_list_out =  [[2], [3], [4], [1]]

        adj_mat_in =   [0  0  0  1
                        1  0  0  0
                        0  1  0  0
                        0  0  1  0]
        adj_list_in =  [[4], [1], [2], [3]]

        # core functionality
        g = GNNGraph(s, t; graph_type=GRAPH_T)
        if TEST_GPU
            g_gpu = g |> gpu
        end

        @test g.num_edges == 4
        @test g.num_nodes == 4
        @test collect(edges(g)) |> sort == collect(zip(s, t)) |> sort
        @test sort(outneighbors(g, 1)) == [2]
        @test sort(inneighbors(g, 1)) == [4]
        @test is_directed(g) == true
        @test is_directed(typeof(g)) == true
        s1, t1 = sort_edge_index(edge_index(g))
        @test s1 == s
        @test t1 == t

        # adjacency
        @test adjacency_matrix(g) ==  adj_mat_out
        @test adjacency_list(g) ==  adj_list_out
        @test adjacency_matrix(g, dir=:out) ==  adj_mat_out
        @test adjacency_list(g, dir=:out) ==  adj_list_out
        @test adjacency_matrix(g, dir=:in) ==  adj_mat_in
        @test adjacency_list(g, dir=:in) ==  adj_list_in
    end

    @testset "Graphs.jl constructor" begin
        lg = random_regular_graph(10, 4)
        @test !Graphs.is_directed(lg)
        g = GNNGraph(lg)
        @test g.num_edges == 2*ne(lg) # g in undirected
        @test Graphs.is_directed(g)
        for e in Graphs.edges(lg)
            i, j = src(e), dst(e)
            @test has_edge(g, i, j)
            @test has_edge(g, j, i)
        end

        @testset "SimpleGraph{Int32}" begin
            g = GNNGraph(SimpleGraph{Int32}(6), graph_type=GRAPH_T)
            @test g.num_nodes == 6
        end
    end

    @testset "Features" begin
        g = GNNGraph(sprand(10, 10, 0.3), graph_type=GRAPH_T)

        # default names
        X = rand(10, g.num_nodes)
        E = rand(10, g.num_edges)
        U = rand(10, g.num_graphs)

        g = GNNGraph(g, ndata=X, edata=E, gdata=U)
        @test g.ndata.x === X
        @test g.edata.e === E
        @test g.gdata.u === U

        # Check no args
        g = GNNGraph(g)
        @test g.ndata.x === X
        @test g.edata.e === E
        @test g.gdata.u === U

        # multiple features names
        g = GNNGraph(g, ndata=(x2=2X, g.ndata...), edata=(e2=2E, g.edata...), gdata=(u2=2U, g.gdata...))
        @test g.ndata.x === X
        @test g.edata.e === E
        @test g.gdata.u === U
        @test g.ndata.x2 ≈ 2X
        @test g.edata.e2 ≈ 2E
        @test g.gdata.u2 ≈ 2U

        # Dimension checks
        @test_throws AssertionError GNNGraph(erdos_renyi(10, 30), edata=rand(29), graph_type=GRAPH_T)
        @test_throws AssertionError GNNGraph(erdos_renyi(10, 30), edata=rand(2, 29), graph_type=GRAPH_T)
        @test_throws AssertionError GNNGraph(erdos_renyi(10, 30), edata=(; x=rand(30), y=rand(29)), graph_type=GRAPH_T)

        # Copy features on reverse edge
        e = rand(30)
        g = GNNGraph(erdos_renyi(10,  30), edata=e, graph_type=GRAPH_T)
        @test g.edata.e == [e; e]

        # non-array global
        g = rand_graph(10,  30, gdata="ciao", graph_type=GRAPH_T)
        @test g.gdata.u == "ciao"

        # vectors stays vectors
        g = rand_graph(10,  30, ndata=rand(10),
                                edata=rand(30),
                                gdata=(u=rand(2), z=rand(1), q=1),
                                graph_type=GRAPH_T)
        @test size(g.ndata.x) == (10,)
        @test size(g.edata.e) == (30,)
        @test size(g.gdata.u) == (2, 1)
        @test size(g.gdata.z) == (1,)
        @test g.gdata.q === 1

        # Error for non-array ndata
        @test_throws AssertionError rand_graph(10,  30, ndata="ciao", graph_type=GRAPH_T)
        @test_throws AssertionError rand_graph(10,  30, ndata=1, graph_type=GRAPH_T)
    end

    @testset "LearnBase and DataLoader compat" begin
        n, m, num_graphs = 10, 30, 50
        X = rand(10, n)
        E = rand(10, m)
        U = rand(10, 1)
        data = [rand_graph(n, m, ndata=X, edata=E, gdata=U, graph_type=GRAPH_T)
                for _ in 1:num_graphs]
        g = Flux.batch(data)

        @testset "batch then pass to dataloader" begin
            @test LearnBase.getobs(g, 3) == getgraph(g, 3)
            @test LearnBase.getobs(g, 3:5) == getgraph(g, 3:5)
            @test StatsBase.nobs(g) == g.num_graphs

            d = Flux.Data.DataLoader(g, batchsize=2, shuffle=false)
            @test first(d) == getgraph(g, 1:2)
        end

        @testset "pass to dataloader and collate" begin
            @test LearnBase.getobs(data, 3) == getgraph(g, 3)
            @test LearnBase.getobs(data, 3:5) == getgraph(g, 3:5)
            @test StatsBase.nobs(data) == g.num_graphs

            d = Flux.Data.DataLoader(data, batchsize=2, shuffle=false)
            @test first(d) == getgraph(g, 1:2)
        end
    end

    @testset "Graphs.jl integration" begin
        g = GNNGraph(erdos_renyi(10, 20), graph_type=GRAPH_T)
        @test g isa Graphs.AbstractGraph
    end

    @testset "==" begin
        g1 = rand_graph(5, 6, ndata=rand(5), edata=rand(6), graph_type=GRAPH_T)
        @test g1 == g1
        @test g1 == deepcopy(g1)
        @test g1 !== deepcopy(g1)

        g2 = GNNGraph(g1, graph_type=GRAPH_T)
        @test g1 == g2
        @test g1 === g2 # this is true since GNNGraph is immutable

        g2 = GNNGraph(g1, ndata=rand(5), graph_type=GRAPH_T)
        @test g1 != g2
        @test g1 !== g2

        g2 = GNNGraph(g1, edata=rand(6), graph_type=GRAPH_T)
        @test g1 != g2
        @test g1 !== g2
    end

    @testset "hash" begin
        g1 = rand_graph(5, 6, ndata=rand(5), edata=rand(6), graph_type=GRAPH_T)
        @test hash(g1) == hash(g1)
        @test hash(g1) == hash(deepcopy(g1))
        @test hash(g1) == hash(GNNGraph(g1, ndata=g1.ndata, graph_type=GRAPH_T))
        @test hash(g1) == hash(GNNGraph(g1, ndata=g1.ndata, graph_type=GRAPH_T))
        @test hash(g1) != hash(GNNGraph(g1, ndata=rand(5), graph_type=GRAPH_T))
        @test hash(g1) != hash(GNNGraph(g1, edata=rand(6), graph_type=GRAPH_T))
    end
end


text:
using JET, InteractiveUtils

const CC = JET.CC

import .CC:
    widenconst, ⊑

import JET:
    JETInterpreter,
    AbstractGlobal,
    analyze_file,
    analyze_text,
    get_result,
    ToplevelConfig,
    virtual_process,
    gen_virtual_module,
    ToplevelErrorReport,
    InferenceErrorReport,
    ExceptionReport,
    print_reports

for sym in Symbol.(last.(Base.Fix2(split, '.').(string.(vcat(subtypes(JET, ToplevelErrorReport),
                                                             subtypes(JET, InferenceErrorReport),
                                                             subtypes(JET, ExceptionReport),
                                                             )))))
    Core.eval(@__MODULE__, :(import JET: $(sym)))
end

import Base:
    Fix1,
    Fix2

import Base.Meta:
    isexpr

"""
    @fixturedef

Creates a virtual module and defines fixtures from toplevel expression `ex`
"""
macro fixturedef(ex)
    @assert isexpr(ex, :block)
    return quote let
        vmod = $gen_virtual_module()
        for x in $(ex.args)
            Core.eval(vmod, x)
        end
        vmod # return virtual module
    end end
end

"""
    @analyze_toplevel [jetconfigs...] ex
    @analyze_toplevel [jetconfigs...] vmod ex

Enters JET analysis from toplevel expression `ex`, and returns the analysis result.
"""
macro analyze_toplevel(xs...)
    jetconfigs = filter(iskwarg, xs)
    xs′ = filter(!iskwarg, xs)
    n = length(xs′)
    @assert 1 ≤ n ≤ 2
    actualmod = __module__
    virtualmod, ex = n == 1 ? (gen_virtual_module(actualmod), first(xs′)) : (xs′...,)
    return _analyze_toplevel(ex, __source__, actualmod, virtualmod, jetconfigs)
end

iskwarg(@nospecialize(x)) = isexpr(x, :(=))

function _analyze_toplevel(ex, lnn, actualmod, virtualmod, jetconfigs)
    toplevelex = (isexpr(ex, :block) ?
                  Expr(:toplevel, lnn, ex.args...) : # flatten here
                  Expr(:toplevel, lnn, ex)
                  ) |> QuoteNode
    return :(let
        actualmod = $(esc(actualmod))
        virtualmod = $(esc(virtualmod))
        interp = JETInterpreter(; $(map(esc, jetconfigs)...))
        config = ToplevelConfig(; $(map(esc, jetconfigs)...))
        $virtual_process($toplevelex,
                         $(string(lnn.file)),
                         actualmod,
                         interp,
                         config,
                         virtualmod,
                         )
    end)
end

is_concrete(mod, sym) = isdefined(mod, sym) && !isa(getfield(mod, sym), AbstractGlobal)
is_abstract(mod, sym) = isdefined(mod, sym) && isa(getfield(mod, sym), AbstractGlobal)
isa_abstract(x, @nospecialize(typ)) = isa(x, AbstractGlobal) && x.t ⊑ typ

# JET will try to concretize global variable when its type is a constant at analysis time,
# but the starategy is a bit complicated right now and may change in the future
# these utilities allow robust testing to check if a object is successfully analyzed by JET whichever it's concretized or abstracted
is_analyzed(mod, sym) = isdefined(mod, sym) # essentially, `is_concrete(mod, sym) || is_abstract(mod, sym)`
isa_analyzed(x, @nospecialize(typ)) = isa_abstract(x, typ) || isa(x, typ)

# fresh execution/benchmark tools
include(normpath(@__DIR__, "..", "benchmark", "JETBenchmarkUtils.jl"))
using .JETBenchmarkUtils


text:
##

GtkIDE.reload()

##

Gtk.GLib.nth

#eval( Core.atdoc(:(Gtk.A.rand)) )

##

#    sc = Gtk.G_.style_context(t.view)
#    push!(sc, provider, 600)

sp = GtkIDE.main_window.style_and_language_manager.main_style
style = GtkSourceWidget.style(sp,"text")

mc = GtkIDE.MarkdownTextViews.MarkdownColors(
    Gtk.get_gtk_property(style,:foreground,String),
    Gtk.get_gtk_property(style,:background,String),
    Gtk.get_gtk_property(GtkSourceWidget.style(sp,"def:note"),:foreground,String),
    Gtk.get_gtk_property(style,:background,String),
)

v = GtkIDE.MarkdownTextViews.MarkdownTextView(string(@doc collect),mc)

#GtkSourceWidget.GtkSourceStyle(

w = Gtk.GtkWindow("asd",200,300) |> v
showall(w)

##
s = """
# fill *asd*

```
fill(x, dims)
```

Create an array filled with the value `x`. For example, `fill(1.0, (5,5))` returns a 5×5 array of floats, with each element initialized to `1.0`.

Test *italic*

"""

m = string(@doc collect)

m = Base.Markdown.parse(m)

el = m.content[1]

##

using Gtk
import Gtk.GtkTextIter

import Markdown

function tag(buffer,what,i,j)
    Gtk.apply_tag(buffer,what,
        GtkTextIter(buffer,i) , GtkTextIter(buffer,j)
    )
end

function insert_MD!(buffer,m::Markdown.Header,i)
    ip = i

    insert!(buffer,"    ")
    i += 4
    for el in m.text
        i = insert_MD!(buffer,el,i)
    end
    tag(buffer, "h1", ip, i)
    i
end

function insert_MD!(buffer,m::String,i)
    insert!(buffer,m)
    i += length(m)
end

function insert_MD!(buffer,m::Markdown.LaTeX,i)
    i = insert_MD!(buffer,m.formula,i)
end

function insert_MD!(buffer,m::Markdown.Paragraph,i)
#    insert!(buffer,"\n\n")
#    i += 2
    for el in m.content
        i = insert_MD!(buffer,el,i)
    end
    i
end

function insert_MD!(buffer,m::Markdown.Code,i)
    insert!(buffer,m.code)
    tag(buffer, "code", i, i+sizeof(m.code))
    i += length(m.code)
end

function insert_MD!(buffer,m::Markdown.List,i)
    for it in m.items
        insert!(buffer,"    - ")
        i += 6
        for el in it
            i = insert_MD!(buffer,el,i)
        end
        insert!(buffer,"\n")
        i += 1
    end
    i
end

function insert_MD!(buffer,m::Markdown.Italic,i)

    ip = i
    for el in m.text
        i = insert_MD!(buffer,el,i)
    end
    tag(buffer, "italic", ip, i)
    i
end

function insert_MD!(buffer,m,i)
    if isdefined(m,:text)
        for el in m.text
            i = insert_MD!(buffer,el,i)
        end
    end
    if isdefined(m,:content)
        for el in m.content
            i = insert_MD!(buffer,el,i)
        end
    end
    i
end

function insert_MD!(buffer,m::Markdown.MD)
    i = 1
    for el in m.content
        i = insert_MD!(buffer,el,i)
        insert!(buffer,"\n\n")
        i += 2
    end
end

mutable struct MarkdownTextView <: GtkTextView

    handle::Ptr{Gtk.GObject}
    view::GtkTextView
    buffer::GtkTextBuffer

    function MarkdownTextView(txt::Markdown.MD)

        buffer = GtkTextBuffer()

        view = GtkTextView(buffer)

        set_gtk_property!(view,:margin_left,1)
        set_gtk_property!(view,:monospace,true)
        set_gtk_property!(view,:wrap_mode,true)

        Gtk.create_tag(buffer, "normal", font="13")
        Gtk.create_tag(buffer, "h1", font="Bold 15")
        Gtk.create_tag(buffer, "h2", font="bold 14")
        Gtk.create_tag(buffer, "bold", font="bold")
        Gtk.create_tag(buffer, "italic", font="italic")
        Gtk.create_tag(buffer, "code", font="bold", background="#eee")

        insert_MD!(buffer,m)
        tag(buffer,"normal",1,length(buffer))

        n = new(view.handle,view,buffer)
        Gtk.gobject_move_ref(n, view)
    end

end

##

#buffer = GtkTextBuffer()
#set_gtk_property!(buffer,:text,"wesh")

view = MarkdownTextView(m)

w = GtkWindow("",600,400)
push!(w,view)
showall(w)


text:
#   Unit examples of the FourierAnalysis Package for julia language
#
#   MIT License
#   Copyright (c) 2019-2022,
#   Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This example shows how to compute and plot time-frequency analytic signal,
#   amplitude and phase objects and how to extract inforation from them.

#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #
#                                                                             #
#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #

using FourierAnalysis, FFTW, LinearAlgebra, Statistics, Plots, Plots.Measures

# add module for reading the two EEG text files to be used ater
push!(LOAD_PATH, @__DIR__)
using IOtxt

# Get EEG file names with complete path (they have extension .txt)
S=getFilesInDir(@__DIR__; ext=(".txt",))

# read the two EEG data files and put them in a Matrix object
X1=readEEG(S[1])
X2=readEEG(S[2])

####################
# test time-frequency

using Plots, FourierAnalysis
sr, t, bandwidth=128, 512, 2
h=taper(harris4, t)
x1=sinusoidal(10, 8, sr, t, 0)
x2=sinusoidal(10, 19, sr, t, 0)
x=Vector((x1+x2).*h.y+randn(t))
y1=sinusoidal(10, 6, sr, t, 0)
y2=sinusoidal(10, 16, sr, t, 0)
y=Vector((y1+y2).*h.y+randn(t))

# Vector of Time-Frequency Object for x and y
𝒀 = TFanalyticsignal([x, y], sr, t)

# compute the mean in a TF region (8:12Hz and samples 1:128) for the first object
y=mean(𝒀[1], (8, 12), (1, 128))

# extract a TF region (8:12Hz and samples 1:128) for the first object
Z=extract(𝒀[1], (8, 12), (1, 128))

# compute the mean in a TF region (8:12Hz and samples 1:128) for the two objects
Y=mean(𝒀, (8, 12), (1, 128))
# do the same computation without checking homogeneity of the two objects in 𝒀
Y=mean(𝒀, (8, 12), (1, 128); check=false)

# extract the data in a TF region (8:12Hz and samples 1:128) for the two objects
E=extract(𝒀, (8, 12), (8, 12))
# do the same computation without checking homogeneity of the two objects in 𝒀
E=extract(𝒀, (8, 12), (8, 12); check=false)

# compute time-frequency object for vector x
Y=TFanalyticsignal(x, sr, t, bandwidth; fmax=32)

# gather useful attributes to obtain nice heatmpap plots
tfArgs=(right_margin = 2mm,
        top_margin = 2mm,
        xtickfont = font(10, "Times"),
        ytickfont = font(10, "Times"))

# plot the real part of the AS
heatmap(tfAxes(Y)..., real(Y.y); c=:fire, tfArgs...)

# ...the imaginary part of the AS
heatmap(tfAxes(Y)..., imag(Y.y); c=:bluesreds, tfArgs...)

# ...the amplitude of the AS
heatmap(tfAxes(Y)..., amplitude(Y.y); c=:amp, tfArgs...)

# or generate a TFAmplitude object
A=TFamplitude(Y)
# and plot it (with different colors)
heatmap(tfAxes(A)..., A.y; c=:fire, tfArgs...)

# ...the amplitude of the AS smoothed in the freq. dim.
heatmap(tfAxes(Y)...,
        amplitude(smooth(hannSmoother, noSmoother, Y).y);
        c=:amp, tfArgs...)

# ...the amplitude of the AS smoothed in freq. and time
heatmap(tfAxes(Y)...,
        amplitude(smooth(blackmanSmoother, blackmanSmoother, Y).y);
        c=:amp, tfArgs...)

# ...the phase of the AS
heatmap(tfAxes(Y)..., phase(Y.y);
        c=:bluesreds, tfArgs...)

#or, generate a TFPhase object
ϴ=TFphase(Y)
# and plot it
heatmap(tfAxes(ϴ)..., ϴ.y;
        c=:fire, tfArgs...)

# compute and plot phase in [0, 2π]
heatmap(tfAxes(Y)..., TFphase(Y; func=x->x+π).y;
        c=:amp, tfArgs...)

# compute and plot unwrapped phase
heatmap(tfAxes(Y)..., TFphase(Y; unwrapped=true).y;
        c=:amp, tfArgs...)

# smooth time-frequency analytic signal: smooth frequency
Z=smooth(blackmanSmoother, noSmoother, Y)
# plot amplitude of smoothed analytic signal
heatmap(tfAxes(Z)..., amplitude(Z.y); c=:amp, tfArgs...)

# not equivalently (!), you can create an amplitude object and smooth them:
# in this case the amplitude is smoothed, not the analytic signal
A=smooth(blackmanSmoother, noSmoother, TFamplitude(Y))
heatmap(tfAxes(A)..., A.y; c=:amp, tfArgs...)

# Smoothing raw phase estimates is unappropriate
# since the phase is a discontinous function, however it makes sense to smooth
# phase if the phase is unwrapped.
heatmap(tfAxes(Y)...,
        smooth(blackmanSmoother, noSmoother, TFphase(Y; unwrapped=true)).y;
        c=:amp, tfArgs...)

# smooth AS: smooth both frequency and time
E=smooth(blackmanSmoother, blackmanSmoother, Y)

# plot amplitude of smoothed analytic signal
heatmap(tfAxes(E)..., amplitude(E.y);
        c=:fire, tfArgs...)

# plot phase of smoothed analytic signal
heatmap(tfAxes(E)..., phase(E.y);
        c=:bluesreds, tfArgs...)

# not equivalently (!), create amplitude and phase objects and smooth them
A=smooth(blackmanSmoother, blackmanSmoother, TFamplitude(Y))
heatmap(tfAxes(A)..., A.y;
        c=:fire, tfArgs...)

ϴ=smooth(blackmanSmoother, blackmanSmoother, TFphase(Y, unwrapped=true))
heatmap(tfAxes(ϴ)..., ϴ.y;
        c=:fire, tfArgs...)

# smooth again
ϴ=smooth(blackmanSmoother, blackmanSmoother, ϴ)
heatmap(tfAxes(ϴ)..., ϴ.y;
        c=:fire, tfArgs...)
# and again ...

# you may also create all these objects already smoothed, for example
# create directly smoothed AS
Y=TFanalyticsignal(x, sr, t, bandwidth;
                   fmax=32,
                   fsmoothing=hannSmoother,
                   tsmoothing=hannSmoother)

# plot amplitude of smoothed analytic signal
heatmap(tfAxes(Y)..., amplitude(Y.y);
   c=:amp, tfArgs...)

# create directly smoothed Amplitude
A=TFamplitude(x, sr, t, bandwidth;
         fmax=32,
         fsmoothing=hannSmoother,
         tsmoothing=hannSmoother)

# plot smoothed amplitude
heatmap(tfAxes(A)..., A.y;
   c=:amp, tfArgs...)

# compute a TFAnalyticSignal object with non-linear AS
Y=TFanalyticsignal(x, sr, t, bandwidth; fmax=32, nonlinear=true)

# check that it is non-linear
Y.nonlinear

# check that the amplitude is now 1.0 everywhere
norm(amplitude(Y.y)-ones(eltype(Y.y), size(Y.y))) # must be zero

# plot non-linear phase
heatmap(tfAxes(Y)..., phase(Y.y);
   c=:bkr, tfArgs...)

# get the center frequencies of TFAmplitude object A
A.flabels

# extract the amplitude in a time-frequency region
extract(A, (2, 10), (1, 256))
# extract the amplitude in a time-frequency region at only one frequency
extract(A, 10, (1, 256))
# extract the amplitude at one temporal sample at one frequency
extract(A, 10, 12)
# or
extract(A, 10.0, 12)
# extract amplitude at one temporal sample in a frequency band-pass
extract(A, (10, 12), 12)
#or
extract(A, (10.0, 12.0), 12)
# extract amplitude at one temporal sample and all frequencies
extract(A, :, 12)
# NB: the extract function work in the same way
#     for objects TFAnayticSignal and TFPhase

# compute the mean in a time-frequency region:
mean(A, (2, 10), (1, 256))
# is equivalent to (but may be less efficient than)
mean(extract(A, (2, 10), (1, 256)))

# using column sign for extracting all time samples
extract(A, (2, 10), :)
# This:
extract(A, :, :)
# is equivalent to this:
A.y
# but if you don't need to extract all frequencies,
# use the extract function to control what frequencies will be extracted:
# This
extract(A, (4, 5), 10)
# is not equivalent to this
A.y[4:5, 10]
# since the `extract` function finds the correct rows corresponding
# to the sought frequencies (in Hz), while A.y[4:5, 10]
# just returns the element [4:5, 10] in the TF amplitude object

mean(A, 10, (1, 256))
mean(A, 10, 12)
# this is not equivalent to the following
A.y[10, 12]
# since the `mean` function also finds the correct rows corresponding
# to the sought frequencies (in Hz), while A.y[10, 12]
# just returns the element [10, 12] in the TF amplitude object

# Although the first center frequency in A is 2Hz, its
# band-pass region is 1-3Hz, therefore the frequency range 1:10 is accepted
mean(A, (1, 10), (1, 256))
# but this result in an error (see the REPL) since 0.5 is out of range
mean(A, (0.5, 10), (1, 256))

# using a colon sign for time range
a=mean(A, (1, 10), :)
# using an integer for time range (indicates one specific sample)
a=mean(A, (1, 10), 16)

# using a colon sign for frequency range
a=mean(A, :, (1, 16))
# using a real number for frequency range
a=mean(A, 8.5, (1, 16))

# This
a=mean(A, :, :)
# is equivalent to this
a=mean(A.y)

###############################################################

Pz=15
𝐱=[X1[:, Pz], X2[:, Pz]] # get the two times-series at electrode Pz
𝐘=TFanalyticsignal(𝐱, sr, t, bandwidth; fmax=32, nonlinear=false)
𝐀=TFamplitude(𝐘)
heatmap(tfAxes(𝐀[1])..., 𝐀[1].y; c=:amp, tfArgs...)
heatmap(tfAxes(𝐀[2])..., 𝐀[2].y; c=:fire, tfArgs...)

# plot the power over time from instantaneous amplitude
plot([sum(𝐀[2].y[:, t]) for t=1:size(𝐀[2].y, 2)])
# plot the smoothed power over time from the instantaneous amplitude
A=smooth(noSmoother, blackmanSmoother, 𝐀[2])
plot([sum(A.y[:, t]) for t=1:size(A.y, 2)])

# 𝐀 holds 3 TFAmplitude objects;
# The following will create a 3-Vector of matrices
# holding the amplitude in the TF region defined by frequency 8:12 and samples 1:128:
extract(𝐀, (8, 12), (1, 128))

# use of a column is allowed for the frequency range as all the amplitudes have the same number of frequencies:
extract(𝐀, :, (1, 128))

# however, the use of a column is not allowed for the time range as the number
# of samples is different in the three TFAmplitude objects
extract(𝐀, (8, 12), :)

# the same goes for the mean function
mean(𝐀, (8, 12), :)

# to obtain the grand-average across all objects:
a=mean(mean(𝐀, (8, 12), (1, 128)))


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# ––––––––––
# Paragraphs
# ––––––––––

type Paragraph
    content
end

Paragraph() = Paragraph([])

function paragraph(stream::IO, md::MD)
    buffer = IOBuffer()
    p = Paragraph()
    push!(md, p)
    skipwhitespace(stream)
    prev_char = '\n'
    while !eof(stream)
        char = read(stream, Char)
        if char == '\n' || char == '\r'
            char == '\r' && !eof(stream) && Char(peek(stream)) == '\n' && read(stream, Char)
            if prev_char == '\\'
                write(buffer, '\n')
            elseif blankline(stream) || parse(stream, md, breaking = true)
                break
            else
                write(buffer, ' ')
            end
        else
            write(buffer, char)
        end
        prev_char = char
    end
    p.content = parseinline(seek(buffer, 0), md)
    return true
end

# –––––––
# Headers
# –––––––

type Header{level}
    text
end

Header(s, level::Int) = Header{level}(s)
Header(s) = Header(s, 1)

@breaking true ->
function hashheader(stream::IO, md::MD)
    withstream(stream) do
        eatindent(stream) || return false
        level = 0
        while startswith(stream, '#') level += 1 end
        level < 1 || level > 6 && return false

        c = ' '
        # Allow empty headers, but require a space
        !eof(stream) && (c = read(stream, Char); !(c in " \n")) &&
            return false

        if c != '\n' # Empty header
            h = strip(readline(stream))
            h = match(r"(.*?)( +#+)?$", h).captures[1]
            buffer = IOBuffer()
            print(buffer, h)
            push!(md.content, Header(parseinline(seek(buffer, 0), md), level))
        else
            push!(md.content, Header("", level))
        end
        return true
    end
end

function setextheader(stream::IO, md::MD)
    withstream(stream) do
        eatindent(stream) || return false
        header = strip(readline(stream))
        header == "" && return false

        eatindent(stream) || return false
        underline = strip(readline(stream))
        length(underline) < 3 && return false
        u = underline[1]
        u in "-=" || return false
        all(c -> c == u, underline) || return false
        level = (u == '=') ? 1 : 2

        push!(md.content, Header(parseinline(header, md), level))
        return true
    end
end

# ––––
# Code
# ––––

type Code
    language::String
    code::String
end

Code(code) = Code("", code)

function indentcode(stream::IO, block::MD)
    withstream(stream) do
        buffer = IOBuffer()
        while !eof(stream)
            if startswith(stream, "    ") || startswith(stream, "\t")
                write(buffer, readline(stream, chomp=false))
            elseif blankline(stream)
                write(buffer, '\n')
            else
                break
            end
        end
        code = String(take!(buffer))
        !isempty(code) && (push!(block, Code(rstrip(code))); return true)
        return false
    end
end

# --------
# Footnote
# --------

type Footnote
    id::String
    text
end

function footnote(stream::IO, block::MD)
    withstream(stream) do
        regex = r"^\[\^(\w+)\]:"
        str = startswith(stream, regex)
        if isempty(str)
            return false
        else
            ref = match(regex, str).captures[1]
            buffer = IOBuffer()
            write(buffer, readline(stream, chomp=false))
            while !eof(stream)
                if startswith(stream, "    ")
                    write(buffer, readline(stream, chomp=false))
                elseif blankline(stream)
                    write(buffer, '\n')
                else
                    break
                end
            end
            content = parse(seekstart(buffer)).content
            push!(block, Footnote(ref, content))
            return true
        end
    end
end

# ––––––
# Quotes
# ––––––

type BlockQuote
    content
end

BlockQuote() = BlockQuote([])

# TODO: Laziness
@breaking true ->
function blockquote(stream::IO, block::MD)
    withstream(stream) do
        buffer = IOBuffer()
        empty = true
        while eatindent(stream) && startswith(stream, '>')
            startswith(stream, " ")
            write(buffer, readline(stream, chomp=false))
            empty = false
        end
        empty && return false

        md = String(take!(buffer))
        push!(block, BlockQuote(parse(md, flavor = config(block)).content))
        return true
    end
end

# -----------
# Admonitions
# -----------

type Admonition
    category::String
    title::String
    content::Vector
end

@breaking true ->
function admonition(stream::IO, block::MD)
    withstream(stream) do
        # Admonition syntax:
        #
        # !!! category "optional explicit title within double quotes"
        #     Any number of other indented markdown elements.
        #
        #     This is the second paragraph.
        #
        startswith(stream, "!!! ") || return false
        # Extract the category of admonition and its title:
        category, title =
            let untitled = r"^([a-z]+)$",          # !!! <CATEGORY_NAME>
                titled   = r"^([a-z]+) \"(.*)\"$", # !!! <CATEGORY_NAME> "<TITLE>"
                line     = strip(readline(stream))
                if ismatch(untitled, line)
                    m = match(untitled, line)
                    # When no title is provided we use CATEGORY_NAME, capitalising it.
                    m.captures[1], ucfirst(m.captures[1])
                elseif ismatch(titled, line)
                    m = match(titled, line)
                    # To have a blank TITLE provide an explicit empty string as TITLE.
                    m.captures[1], m.captures[2]
                else
                    # Admonition header is invalid so we give up parsing here and move
                    # on to the next parser.
                    return false
                end
            end
        # Consume the following indented (4 spaces) block.
        buffer = IOBuffer()
        while !eof(stream)
            if startswith(stream, "    ")
                write(buffer, readline(stream, chomp=false))
            elseif blankline(stream)
                write(buffer, '\n')
            else
                break
            end
        end
        # Parse the nested block as markdown and create a new Admonition block.
        nested = parse(String(take!(buffer)), flavor = config(block))
        push!(block, Admonition(category, title, nested.content))
        return true
    end
end

# –––––
# Lists
# –––––

type List
    items::Vector{Any}
    ordered::Int # `-1` is unordered, `>= 0` is ordered.

    List(x::AbstractVector, b::Integer) = new(x, b)
    List(x::AbstractVector) = new(x, -1)
    List(b::Integer) = new(Any[], b)
end

List(xs...) = List(vcat(xs...))

isordered(list::List) = list.ordered >= 0

const BULLETS = r"^ {0,3}(\*|\+|-)( |$)"
const NUM_OR_BULLETS = r"^ {0,3}(\*|\+|-|\d+(\.|\)))( |$)"

@breaking true ->
function list(stream::IO, block::MD)
    withstream(stream) do
        bullet = startswith(stream, NUM_OR_BULLETS; eat = false)
        indent = isempty(bullet) ? (return false) : length(bullet)
        # Calculate the starting number and regex to use for bullet matching.
        initial, regex =
            if ismatch(BULLETS, bullet)
                # An unordered list. Use `-1` to flag the list as unordered.
                -1, BULLETS
            elseif ismatch(r"^ {0,3}\d+(\.|\))( |$)", bullet)
                # An ordered list. Either with `1. ` or `1) ` style numbering.
                r = contains(bullet, ".") ? r"^ {0,3}(\d+)\.( |$)" : r"^ {0,3}(\d+)\)( |$)"
                Base.parse(Int, match(r, bullet).captures[1]), r
            else
                # Failed to match any bullets. This branch shouldn't actually be needed
                # since the `NUM_OR_BULLETS` regex should cover this, but we include it
                # simply for thoroughness.
                return false
            end

        # Initialise the empty list object: either ordered or unordered.
        list = List(initial)

        buffer = IOBuffer() # For capturing nested text for recursive parsing.
        newline = false     # For checking if we have two consecutive newlines: end of list.
        count = 0           # Count of list items. Used to check if we need to push remaining
                            # content in `buffer` after leaving the `while` loop.
        while !eof(stream)
            if startswith(stream, "\n")
                if newline
                    # Double newline ends the current list.
                    pushitem!(list, buffer)
                    break
                else
                    newline = true
                    println(buffer)
                end
            else
                newline = false
                if startswith(stream, " "^indent)
                    # Indented text that is part of the current list item.
                    print(buffer, readline(stream, chomp=false))
                else
                    matched = startswith(stream, regex)
                    if isempty(matched)
                        # Unindented text meaning we have left the current list.
                        pushitem!(list, buffer)
                        break
                    else
                        # Start of a new list item.
                        count += 1
                        count > 1 && pushitem!(list, buffer)
                        print(buffer, readline(stream, chomp=false))
                    end
                end
            end
        end
        count == length(list.items) || pushitem!(list, buffer)
        push!(block, list)
        return true
    end
end
pushitem!(list, buffer) = push!(list.items, parse(String(take!(buffer))).content)

# ––––––––––––––
# HorizontalRule
# ––––––––––––––

type HorizontalRule
end

function horizontalrule(stream::IO, block::MD)
   withstream(stream) do
       n, rule = 0, ' '
       while !eof(stream)
           char = read(stream, Char)
           char == '\n' && break
           isspace(char) && continue
           if n==0 || char==rule
               rule = char
               n += 1
           else
               return false
           end
       end
       is_hr = (n ≥ 3 && rule in "*-")
       is_hr && push!(block, HorizontalRule())
       return is_hr
   end
end


text:
using BenchmarkTools, Plots, SNOPT7

T = Float64

# options
max_con_viol = 1.0e-8
verbose=false

opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,
    live_plotting=:off)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
    opts_uncon=opts_ilqr,
    cost_tolerance=1.0e-5,
    cost_tolerance_intermediate=1.0e-2,
    constraint_tolerance=max_con_viol,
    penalty_scaling=100.,
    penalty_initial=1.)

opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
    feasibility_tolerance=max_con_viol)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,
    opts_al=opts_al,
    opts_pn=opts_pn,
    projected_newton=true,
    projected_newton_tolerance=1.0e-4);

opts_ipopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=Ipopt.Optimizer(),
    feasibility_tolerance=max_con_viol)

opts_snopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=SNOPT7.Optimizer(),
    feasibility_tolerance=max_con_viol)

# ALTRO w/ Newton
prob_altro = copy(Problems.acrobot)
@time p1, s1 = solve(prob_altro, opts_altro)
@benchmark p1, s1 = solve($prob_altro, $opts_altro)
plot(p1.X,title="Acrobot state (ALTRO)")
plot(p1.U,title="Acrobot control (ALTRO)")
max_violation_direct(p1)

# DIRCOL w/ Ipopt
prob_ipopt = copy(Problems.acrobot)
rollout!(prob_ipopt)
prob_ipopt = update_problem(prob_ipopt,model=Dynamics.acrobot_model) # get continuous time model
@time p2, s2 = solve(prob_ipopt, opts_ipopt)
@benchmark p2, s2 = solve($prob_ipopt, $opts_ipopt)
plot(p2.X,title="Acrobot state (Ipopt)")
plot(p2.U,title="Acrobot control (Ipopt)")
max_violation_direct(p2)

# DIRCOL w/ SNOPT
prob_snopt = copy(Problems.acrobot)
rollout!(prob_snopt)
prob_snopt = update_problem(prob_snopt,model=Dynamics.acrobot_model) # get continuous time model
@time p3, s3 = solve(prob_snopt, opts_snopt)
@benchmark p3, s3 = solve($prob_snopt, $opts_snopt)
plot(p3.X,title="Acrobot state (SNOPT)")
plot(p3.U,title="Acrobot control (SNOPT)")
max_violation_direct(p3)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
strcmp(s1, s2) = s1 == s2
strcat(s...) = prod(s)
num2str(n) = "$n"
lower(s) = lowercase(s)
repmat(v::AbstractArray,s...) = repeat(v, s...)
repmat(v::Real,s...) = repeat([v], s...)

function atan2(y, x)
    return atan(y, x)
end

function permute(args...)
    permutedims(args...)
end

# spdiags([reshape(-D,nJ,1) reshape(D,nJ,1)],[-1 0] ,nJ,nJ)
function spdiags(mat, opt::Real, m, n)
    return spdiagm(m, n, 0 => mat[:,1])
end
function spdiags(mat, opt::AbstractArray, m, n)
    @assert length(opt) == 2
    # @show size(mat)
    # @show opt
    # error("Debugging spdiags")
    if opt[1] == -1 && opt[2] == 0
        ev = mat[:,1][1:m-1] # TODO: verify
        dv = mat[:,2][1:m]
#       return Bidiagonal(dv, ev, :L)
        return spdiagm(m, n, 0 => dv, -1 => ev)
    elseif opt[1] == 0 && opt[2] == 1
        dv = mat[:,1][1:m]
        ev = mat[:,2][2:m] # TODO: verify
#       return Bidiagonal(dv, ev, :U)
        return spdiagm(m, n, 0 => dv, 1 => ev)
    else
        error("Uncaught case")
    end
end
function spdiags1(mat, opt::AbstractArray, m, n)
    @assert length(opt) == 2
    # @show size(mat)
    # @show opt
    # error("Debugging spdiags")
    if opt[1] == -1 && opt[2] == 0
        ev = mat[:,1][1:m-1] # TODO: verify
        dv = mat[:,2][1:m]
        return Bidiagonal(dv, ev, :L)
#       return spdiagm(m, n, 0 => dv, -1 => ev)
    elseif opt[1] == 0 && opt[2] == 1
        dv = mat[:,1][1:m]
        ev = mat[:,2][2:m] # TODO: verify
        return Bidiagonal(dv, ev, :U)
#       return spdiagm(m, n, 0 => dv, 1 => ev)
    else
        error("Uncaught case")
    end
end
# function spdiags(dv, v, m, n)
#     Bidiagonal(dv, ev::V, uplo::Symbol)
# end


text:
#=win.jl - replay the game fromstart to finish, at every step checking if the game
is won. If the game is won, print ”B” for Black winning, and ”W” for white winning. If black resigned,
print ”R”. If white resigned, print ”r”. If the game is on, print ”?”. If the game is a draw, print
”D”.
Accepts 1 command line argument,<filename> => database
=#
include("board.jl")
include("dParse.jl")
module win
using BM
  using SQLite

function winner()
  database = ARGS[1] #/path/to/database/file {string}
  db = SQLite.DB(database)  #Opens the database gamefile
  board = BM.startGame("standard")
  state = board.state
  res = "?"
  maxMove = parse(get(SQLite.query(db, """SELECT max(move_number) from moves;""")[1,1]))
  for x in 1:maxMove  #iterates through each row of the database
    dataMove = SQLite.query(db, """SELECT move_number, move_type, targetx, targety,sourcex,sourcey FROM moves WHERE "move_number" = '$x'""")
    move_type = get(dataMove[1,2])
    if !isnull(dataMove[1,3])
    targetx = parse(get(dataMove[1,3]))
    targety = parse(get(dataMove[1,4]))
    end

    #Case 1: Resigned
    if move_type == "resign"
      if iseven(x) == true
        res="r"
      else #x is odd
        res="R"
      end

    #Case 2: Win
    #Game is won on the move that captures a king
    #I implemented win where, if the target X Y location of the current move is a king, it means the king is captpured so the other team wins
  elseif ( get(state,(targetx,targety),('x','x'))[1] == 'k') #uses function from spuare.jl
      if (get(state,(targetx,targety),('x','x'))[2] == 'b')
        res = "W"
      else #board[targetx][targety].team == "w"
        res = "B"
      end
    end
    if !(isnull(dataMove[1,5])) #type:move
	    sourcex=parse(get(dataMove[1,5]))
	    sourcey=parse(get(dataMove[1,6]))
      state[(targetx,targety)] = deepcopy(state[(sourcex,sourcey)]) #Updates the board before next move
      delete!(state,(sourcex,sourcey))
    end
  end
  return res
end #End function
res = winner()
println(res)
end

max_stars_repo_name: kinoroy/Shogi

max_issues_repo_name: kinoroy/Shulia

max_forks_repo_name: kinoroy/Shulia


text:
using Pardiso
using SparseArrays
ps = PardisoSolver()

A = sparse(rand(2, 2))
B = rand(2, 2)
X = zeros(2, 2)
solve!(ps, X, A, B)


text:
using FixedEffects, CategoricalArrays, Test

v1 = ["ok", missing, "ok", "first"]
v2 = [1, 2, 3, 1]
v3 = [missing, 1, 1, missing]

@test levels(group(v1)) == [1, 2]
@test group(v1).refs == [1, 0, 1, 2]

@test group(v2) == categorical(v2)

@test levels(group(v3)) ==  [1]
@test group(v3).refs ==  [0, 1, 1, 0]

@test levels(group(v1, v2, v3)) == [1]
@test group(group(v1, v2, v3)).refs == [0, 0, 1, 0]

show(FixedEffect(v1, v2))
show(FixedEffect(v1; interaction = v2))
show(FixedEffect(v1, v3 ;interaction = v2))

# test different syntaxes
v1 = categorical(collect(1:1000))
v2 = categorical(fill(1, 1000))
@test group(v1, v2) == collect(1:1000)
@test group(v1, v2) == collect(1:1000)
@test group(v1, v2) == collect(1:1000)


text:
"""
    felix_header(; npoints, bandwidth, frequency)

Create Felix NMR header format for 2D NMR experiment with 2D data matrix size
`npoints = (t2_npoints, t1_npoints)`, spectral width `bandwidth = (t2_bw,
t1_bw)` and given spectrometer `frequency` (unit-compatible with Hz).

# Notes about Felix .dat formats.

The Felix `.dat` format exists to convey raw spectrometer data, spectrometer
configuration and processing parameters to Felix. Felix variables can be
listed within the Felix interface by typing the `lis` ("list symbols")
command.  See <www.felixnmr.com/doc/felix2004L/cref/Output/A_CommandRef_1_all.html>.

Unfortunately `.dat` is rather poorly documented.  Here's what I (Chris
Foster) have been able to deduce about it by looking at old source code and
writing exploratory data files followed by dumping the Felix symbols with the
`lis` command.

`.dat` comes in two (or perhaps three?) flavours:

1. Felix "old .dat format": A format based on FORTRAN sequential, unformatted
   access which is system dependent, with the exact detail depending on the
   system FORTRAN compiler.  We don't deal with this format here, but you can
   look in the GAMMA NMR simulator for a certian amount of example code. It's
   unknown whether Felix2007 can still deal with this format.

2. Felix "new .dat format": A format which is machine and compiler
   independent, and which comes in "frames" (each defining a dimension?),
   rather poorly documented in fxrw.c and wn.f files which come in the felix
   distribution at `gifts/felix/newformat`, and slightly better documented by
   `gifts/felix/BRUKER/mips/ux2flx.f`.  Felix appears to use this format if
   you use "Process1D->Open and Process".

3. Multidimensional `.dat` format for Felix2007: A machine/compiler
   independent format defining parameters for multidimensional NMR
   experiments up to 6 dimensions. This appears to have the same pre-header
   as format 2 above, and it's unclear whether this is distinct or part of
   the same thing. Felix seems to use this for "ProcessND->Open and Process 2D".

We partially implement format 3 here.
"""
function felix_header(; npoints, bandwidth, frequency)
    # Initial part of header
    iwords = zeros(UInt32, 256)

    iwords[1] = 1    # number of frames
    iwords[2] = 0    # data format (0 = IEEE floating point,  1 = 32-bit integer)
    iwords[3] = 1    # ?
    iwords[4] = 32   # frame size
    iwords[5] = 210  # Felix Version number

    fwords = reinterpret(Float32, iwords)

    # ---------------------- 1D data section ------------------------
    # Section for 1D data ??
    # Still doesn't seem quite right...
    iwords[96] = npoints[1]  # number of points
    iwords[97] = 1          # data type (0=real, 1=complex, 2=integer)
    #iwords[98] = domain      # domain (0=time, 1=frequency)
    #iwords[99] = axtype      # axis type (0=none, 1=points, 3=ppm)
    # For 1D?
    fwords[111] = ustrip(u"Hz", bandwidth[1]) # swidth
    fwords[112] = ustrip(u"MHz", frequency)   # sfreq
    #fwords[113] = refpt       # refpt
    #fwords[114] = refHz       # ref
    #fwords[116] = phase0[1]
    #fwords[117] = phase1[1]

    # Old bruker dsp acqusition "decimation" params
    #iwords[118] = dspfvs
    #iwords[119] = decim
    #iwords[120] = grpdly

    # ---------------------- 2D data section ------------------------
    # Need to set the following value to 1, or felix won't read the
    # parameters further down. Is this a flag for ND data?
    iwords[121] = 1

    # The names of parameters on the right hand side here are the names of the
    # associated variables inside Felix.

    # Data size
    iwords[122] = npoints[1]  # nacqs1
    iwords[123] = npoints[2]  # nacqs2
    #iwords[124] = nacqs3
    #iwords[125] = nacqs4
    #iwords[126] = nacqs5
    #iwords[127] = nacqs6

    # Complex vs real acquisition flag "quadf == Quadrature Flag?"
    iwords[128:133] .= 1
    #iwords[128] = quadf1
    #iwords[129] = quadf2
    #iwords[130] = quadf3
    #iwords[131] = quadf4
    #iwords[132] = quadf5
    #iwords[133] = quadf6

    # Spectrometer Frequency (in MHz)  a1rsf1 ... a1rsf6
    fwords[134] = ustrip(u"MHz", frequency)
    fwords[135] = ustrip(u"MHz", frequency)
    #fwords[136] = a1rsf3
    #fwords[137] = a1rsf4
    #fwords[138] = a1rsf5
    #fwords[139] = a1rsf6

    # Sweep width (Hz)   # a1rsw1 .. a1rsw6
    fwords[140] = ustrip(u"Hz", bandwidth[1])
    fwords[141] = ustrip(u"Hz", bandwidth[2])
    #fwords[142] = a1rsw3
    #fwords[143] = a1rsw4
    #fwords[144] = a1rsw5
    #fwords[145] = a1rsw6

    # Spectral referencing (refpt / refHz
    #fwords[146] = a1rpv1
    #fwords[147] = a1rpv2
    #fwords[148] = a1rpv3
    #fwords[149] = a1rpv4
    #fwords[150] = a1rpv5
    #fwords[151] = a1rpv6
    #fwords[152] = a1rpn1
    #fwords[153] = a1rpn2
    #fwords[154] = a1rpn3
    #fwords[155] = a1rpn4
    #fwords[156] = a1rpn5
    #fwords[157] = a1rpn6

    # Phase correction params
    #fwords[158] = d1ph0
    #fwords[159] = d1ph1
    #fwords[160] = d2ph0
    #fwords[161] = d2ph1
    #fwords[162] = d3ph0
    #fwords[163] = d3ph1
    #fwords[164] = d4ph0
    #fwords[165] = d4ph1
    #fwords[166] = d5ph0
    #fwords[167] = d5ph1
    #fwords[168] = d6ph0
    #fwords[169] = d6ph1

    # ?
    #iwords[170] = splmod
    #iwords[171] = qrtord

    # Temperature
    #fwords[172] = rtemp

    # Date and acqusition session
    #iwords[173] = crmon
    #iwords[174] = crday
    #iwords[175] = cryear
    #iwords[176] = acqsch
    #iwords[177] = aqmod1
    #iwords[178] = aqmod2
    #iwords[179] = aqmod3
    #iwords[180] = aqmod4
    #iwords[181] = aqmod5
    #iwords[182] = aqmod6

    # ???
    #iwords[183] = fflag1
    #iwords[184] = fflag2
    #iwords[185] = fflag3
    #iwords[186] = fflag4
    #iwords[187] = fflag5
    #iwords[188] = fflag6

    iwords
end

function write_felix_data(io, data)
    for i=1:size(data,2)
        # We add a `conj` here because Felix seems to expect positive
        # frequencies to be negative offsets from the spectrometer reference.
        # (TODO: Maybe this is the standard convention for NMR data? It seems
        # confusing.)
        fid = conj.(ComplexF32.(data[:,i]))
        # Felix format counts data in 32-bit words
        nwords = 2*length(fid)
        write(io, convert(Int32, nwords))
        write(io, fid)
    end
end

"""
    save_felix(fname, data; bandwidth, frequency)

Create Felix NMR file for 2D NMR experiment with 2D data matrix size
`npoints = size(data)`, spectral width `bandwidth = (t2_bw, t1_bw)` and given
spectrometer `frequency`, which should be unit-compatible with Hz.

`data[i,:]` is assumed to contain the real time FIDs as acquired in the
standard 2D COSY experiment (ie, the "t2" dimension).
"""
save_felix(fname::String, data; kwargs...) = open(fname, "w") do io
    save_felix(io, data; kwargs...)
end

function save_felix(io::IO, data; bandwidth, frequency)
    header = felix_header(npoints=size(data), bandwidth=bandwidth, frequency=frequency)

    # Endian marker
    write(io, UInt32(0x04030201))
    # Header size
    write(io, UInt32(length(header)))
    # Header data
    write(io, header)

    write_felix_data(io, data)
end

#-------------------------------------------------------------------------------
#=
"""
Read a spectrum file in Felix ".mat" format

(Felix is software for 2D spectroscopy.)
"""
function read_felix_spectrum(fname)
    @warn "read_felix_spectrum is mostly untested!" max_log=1
    header,rawdata = open(fname) do io
        read(io, 2^14), read(io, 2^22)
    end

    # Felix spectral data appears to be stored in big endian Float32 in the last
    # section of the file.
    #
    # Furthermore, the data layout is not a simple strided matrix, but looks like a
    # set of 128x32 image tiles.
    #
    # It also looks like this data might be magnitude only (??)
    tilesize = (128,32)
    ntiles = (16,16)
    spectrum_size = tilesize .* ntiles
    data_tiles = reshape(ntoh.(reinterpret(Float32, rawdata)), tilesize..., ntiles...)

    # Reshape data from tiled format into a normal strided format
    reshape(permutedims(data_tiles, (1,3,2,4)), spectrum_size...)
end
=#

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module AnnotationRDF
export annotation_to_rdf, annotation_rdf_node

using Serd
using Catlab, Catlab.WiringDiagrams

using ...Doctrine, ...Ontology
using ..RDFUtils
using ..ConceptRDF: generator_rdf_node
using ..WiringRDF

const R = RDF.Resource

# Constants
###########

const language_properties = Dict(
  :class => "annotatedClass",
  :function => "annotatedFunction",
  :method => "annotatedMethod",
)

# RDF
#####

""" Convert annotation into triples for RDF/OWL ontology.
"""
function annotation_to_rdf(annotation::ObAnnotation, prefix::RDF.Prefix; kw...)
  node = annotation_rdf_node(annotation, prefix)
  stmts = RDF.Statement[
    [ RDF.Triple(node, R("rdf","type"), R("monocl","TypeAnnotation")) ];
    annotation_name_to_rdf(annotation, prefix);
    annotation_language_to_rdf(annotation, prefix);
  ]

  # Definition as expression, assuming it's a basic object.
  gen_node = generator_rdf_node(annotation.definition, prefix)
  push!(stmts, RDF.Triple(node, R("monocl","codeDefinition"), gen_node))

  # Slot annotations.
  for (i, hom) in enumerate(annotation.slots)
    slot = annotation.language[:slots][i]["slot"]
    slot_node = R(prefix.name, "$(node.name):slot$i")
    append!(stmts, [
      RDF.Triple(node, R("monocl","annotatedSlot"), slot_node),
      RDF.Triple(slot_node, R("rdf","type"), R("monocl","SlotAnnotation")),
      RDF.Triple(slot_node, R("monocl","codeSlot"), RDF.Literal(slot)),
    ])
    if head(hom) == :generator
      gen_node = generator_rdf_node(hom, prefix)
      push!(stmts, RDF.Triple(slot_node, R("monocl","codeDefinition"), gen_node))
    end
  end

  stmts
end

function annotation_to_rdf(annotation::HomAnnotation, prefix::RDF.Prefix;
                           include_provenance::Bool=true,
                           include_wiring_diagrams::Bool=true)
  node = annotation_rdf_node(annotation, prefix)
  stmts = RDF.Statement[
    [ RDF.Triple(node, R("rdf","type"), R("monocl","FunctionAnnotation")) ];
    annotation_name_to_rdf(annotation, prefix);
    annotation_language_to_rdf(annotation, prefix);
    annotation_domain_to_rdf(annotation, prefix);
  ]

  # Definition as expression, if it's a basic morphism.
  if head(annotation.definition) == :generator
    gen_node = generator_rdf_node(annotation.definition, prefix)
    push!(stmts, RDF.Triple(node, R("monocl","codeDefinition"), gen_node))
  end

  # Definition as wiring diagram.
  if include_wiring_diagrams
    diagram = to_wiring_diagram(annotation.definition)
    diagram_name = "$(node.name):diagram"
    root_node, diagram_stmts = semantic_graph_to_rdf(diagram,
      expr -> generator_rdf_node(expr, prefix);
      box_rdf_node = box -> R(prefix.name, "$diagram_name:$box"),
      port_rdf_node = (box, port) -> R(prefix.name, "$diagram_name:$box:$port"),
      wire_rdf_node = wire -> R(prefix.name, "$diagram_name:$wire"),
      include_provenance = include_provenance)
    append!(stmts, [
      [ RDF.Triple(node, R("monocl","codeDefinition"), root_node) ];
      diagram_stmts;
    ])
  end

  stmts
end

""" Convert annotation's name to RDF.
"""
function annotation_name_to_rdf(annotation::Annotation, prefix::RDF.Prefix)
  node = annotation_rdf_node(annotation, prefix)
  name = annotation.name
  RDF.Statement[
    RDF.Triple(node, R("monocl","annotatedLanguage"), RDF.Literal(name.language)),
    RDF.Triple(node, R("monocl","annotatedPackage"), RDF.Literal(name.package)),
    RDF.Triple(node, R("monocl","id"), RDF.Literal(name.id)),
  ]
end

""" Convert annotation's language-specific data to RDF.
"""
function annotation_language_to_rdf(annotation::Annotation, prefix::RDF.Prefix)
  node = annotation_rdf_node(annotation, prefix)
  stmts = RDF.Statement[]
  for key in intersect(keys(language_properties), keys(annotation.language))
    value = annotation.language[key]
    values = value isa AbstractArray ? value : [ value ]
    append!(stmts, [
      RDF.Triple(node, R("monocl", language_properties[key]), RDF.Literal(v))
      for v in values
    ])
  end
  stmts
end

""" Convert annotation's language-specific domain and codomain data to RDF.
"""
function annotation_domain_to_rdf(annotation::HomAnnotation, prefix::RDF.Prefix)
  node = annotation_rdf_node(annotation, prefix)
  cell_node = name -> R(prefix.name, "$(node.name):$name")
  inputs, outputs = annotation.language[:inputs], annotation.language[:outputs]
  nin, nout = length(inputs), length(outputs)
  owl_inputs_outputs(node, cell_node, nin, nout, index=true) do cell, is_input, i
    data = (is_input ? inputs : outputs)[i]
    slot = data["slot"]
    [ RDF.Triple(cell, R("monocl","codeSlot"), RDF.Literal(slot)) ]
  end
end

""" Create RDF node for annotation.
"""
function annotation_rdf_node(annotation::Annotation, prefix::RDF.Prefix)::RDF.Node
  name = annotation.name
  node_name = join([name.language, name.package, name.id], ":")
  R(prefix.name, node_name)
end

end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph