text:
not_implemented_for(x) = error("Not implemented for $(x). You might want to put:  `using Makie` into your code!")

#TODO only have one?
const Theme = Attributes

Theme(x::AbstractPlot) = x.attributes

default_theme(scene, T) = Attributes()

function default_theme(scene)
    light = Vec3f0[Vec3f0(1.0,1.0,1.0), Vec3f0(0.1,0.1,0.1), Vec3f0(0.9,0.9,0.9), Vec3f0(20,20,20)]
    Theme(
        color = theme(scene, :color),
        visible = theme(scene, :visible),
        linewidth = 1,
        light = light,
        transformation = automatic,
        model = automatic,
        alpha = 1.0,
        transparency = false,
        overdraw = false,
    )
end

#this defines which attributes in a theme should be removed if another attribute is defined by the user,
#to avoid conflicts later through the pipeline

mutual_exclusive_attributes(::Type{<:AbstractPlot}) = Dict()
"""
    `image(x, y, image)` / `image(image)`

Plots an image on range `x, y` (defaults to dimensions).
"""
@recipe(Image) do scene
    Theme(;
        default_theme(scene)...,
        colormap = [RGBAf0(0,0,0,1), RGBAf0(1,1,1,1)],
        colorrange = automatic,
        fxaa = false,
    )
end

# could be implemented via image, but might be optimized specifically by the backend
"""
    `heatmap(x, y, values)` or `heatmap(values)`

Plots a heatmap as an image on `x, y` (defaults to interpretation as dimensions).
"""
@recipe(Heatmap) do scene
    Theme(;
        default_theme(scene)...,
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        linewidth = 0.0,
        levels = 1,
        fxaa = true,
        interpolate = false
    )
end

"""
    `volume(volume_data)`

Plots a volume. Available algorithms are:
* `:iso` => IsoValue
* `:absorption` => Absorption
* `:mip` => MaximumIntensityProjection
* `:absorptionrgba` => AbsorptionRGBA
* `:indexedabsorption` => IndexedAbsorptionRGBA
"""
@recipe(Volume) do scene
    Theme(;
        default_theme(scene)...,
        fxaa = true,
        algorithm = :iso,
        absorption = 1f0,
        isovalue = 0.5f0,
        isorange = 0.05f0,
        colormap = theme(scene, :colormap),
        colorrange = (0, 1),
        color = nothing,
    )
end
mutual_exclusive_attributes(::Type{<:Volume}) =
    Dict(:colorrange => :color,
         :colormap   => :color,
         )

"""
    `surface(x, y, z)`

Plots a surface, where `(x, y, z)` are supposed to lie on a grid.
"""
@recipe(Surface) do scene
    Theme(;
        default_theme(scene)...,
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        shading = true,
        fxaa = true,
    )
end

"""
    `lines(x, y, z)` / `lines(x, y)` / or `lines(positions)`

Creates a connected line plot for each element in `(x, y, z)`, `(x, y)` or `positions`.
"""
@recipe(Lines) do scene
    Theme(;
        default_theme(scene)...,
        linewidth = 1.0,
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        linestyle = theme(scene, :linestyle),
        fxaa = false
    )
end

"""
    `linesegments(x, y, z)` / `linesegments(x, y)` / `linesegments(positions)`

Plots a line for each pair of points in `(x, y, z)`, `(x, y)`, or `positions`.

**Attributes**:
The same as for [`lines`](@ref)
"""
@recipe(LineSegments) do scene
    default_theme(scene, Lines)
end

# alternatively, mesh3d? Or having only mesh instead of poly + mesh and figure out 2d/3d via dispatch
"""
    `mesh(x, y, z)`, `mesh(mesh_object)`, `mesh(x, y, z, faces)`, or `mesh(xyz, faces)`

Plots a 3D mesh.
"""
@recipe(Mesh) do scene
    Theme(;
        default_theme(scene)...,
        fxaa = true,
        interpolate = false,
        shading = true,
        colormap = theme(scene, :colormap),
        colorrange = automatic,
    )
end

"""
    `scatter(x, y, z)` / `scatter(x, y)` / `scatter(positions)`

Plots a marker for each element in `(x, y, z)`, `(x, y)`, or `positions`.
"""
@recipe(Scatter) do scene
    Theme(;
        default_theme(scene)...,
        marker = theme(scene, :marker),
        markersize = theme(scene, :markersize),
        strokecolor = RGBA(0, 0, 0, 0),
        strokewidth = 0.0,
        glowcolor = RGBA(0, 0, 0, 0),
        glowwidth = 0.0,
        rotations = Billboard(),
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        marker_offset = automatic,
        fxaa = false,
        transform_marker = false, # Applies the plots transformation to marker
        uv_offset_width = Vec4f0(0),
        distancefield = nothing,
    )
end

"""
    `meshscatter(x, y, z)` / `meshscatter(x, y)` / `meshscatter(positions)`

Plots a mesh for each element in `(x, y, z)`, `(x, y)`, or `positions` (similar to `scatter`).
`markersize` is a scaling applied to the primitive passed as `marker`
"""
@recipe(MeshScatter) do scene
    Theme(;
        default_theme(scene)...,
        marker = Sphere(Point3f0(0), 1f0),
        markersize = theme(scene, :markersize),
        rotations = Quaternionf0(0, 0, 0, 1),
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        fxaa = true,
        shading = true
    )
end

"""
    `text(string)`

Plots a text.
"""
@recipe(Text) do scene
    Theme(;
        default_theme(scene)...,
        strokecolor = (:black, 0.0),
        strokewidth = 0,
        font = theme(scene, :font),
        align = (:left, :bottom),
        rotation = 0.0,
        textsize = 20,
        position = Point2f0(0),
    )
end

const atomic_function_symbols = (
        :text, :meshscatter, :scatter, :mesh, :linesegments,
        :lines, :surface, :volume, :heatmap, :image
)

const atomic_functions = getfield.(Ref(AbstractPlotting), atomic_function_symbols)
const Atomic{Arg} = Union{map(x-> Combined{x, Arg}, atomic_functions)...}

function color_and_colormap!(plot, intensity = plot[:color])
    if isa(intensity[], AbstractArray{<: Number})
        haskey(plot, :colormap) || error("Plot $T needs to have a colormap to allow the attribute color to be an array of numbers")
        replace_automatic!(plot, :colorrange) do
            lift(extrema_nan, intensity)
        end
        true
    else
        delete!(plot, :colorrange)
        false
    end
end

"""
    `calculated_attributes!(plot::AbstractPlot)`

Fill in values that can only be calculated when we have all other attributes filled
"""
calculated_attributes!(plot::T) where T = calculated_attributes!(T, plot)

"""
    `calculated_attributes!(trait::Type{<: AbstractPlot}, plot)`
trait version of calculated_attributes
"""
calculated_attributes!(trait, plot) = nothing

function calculated_attributes!(::Type{<: Mesh}, plot)
    need_cmap = color_and_colormap!(plot)
    need_cmap || delete!(plot, :colormap)
    return
end

function calculated_attributes!(::Type{<: Union{Heatmap, Image}}, plot)
    plot[:color] = plot[3]
    color_and_colormap!(plot)
end
function calculated_attributes!(::Type{<: Surface}, plot)
    colors = plot[3]
    if haskey(plot, :color)
        color = plot[:color][]
        if isa(color, AbstractMatrix{<: Number}) && !(color === to_value(colors))
            colors = plot[:color]
        end
    end
    color_and_colormap!(plot, colors)
end
function calculated_attributes!(::Type{<: MeshScatter}, plot)
    color_and_colormap!(plot)
end

function calculated_attributes!(::Type{<: Scatter}, plot)
    # calculate base case
    color_and_colormap!(plot)
    replace_automatic!(plot, :marker_offset) do
        # default to middle
        lift(x-> Vec2f0.((x .* (-0.5f0))), plot[:markersize])
    end
end

function calculated_attributes!(::Type{<: Union{Lines, LineSegments}}, plot)
    color_and_colormap!(plot)
    pos = plot[1][]
    # extend one color per linesegment to be one (the same) color per vertex
    # taken from @edljk  in PR #77
    if haskey(plot, :color) && isa(plot[:color][], AbstractVector) && iseven(length(pos)) && (length(pos) ÷ 2) == length(plot[:color][])
        plot[:color] = lift(plot[:color]) do cols
            map(i-> cols[(i + 1) ÷ 2], 1:(length(cols) * 2))
        end
    end
end

# # to allow one color per edge
# function calculated_attributes!(plot::LineSegments)
#     plot[:color] = lift(plot[:color], plot[1]) do c, p
#         if (length(p) ÷ 2) == length(c)
#             [c[k] for k in 1:length(c), l in 1:2]
#         else
#             c
#         end
#     end
# end

function (PT::Type{<: Combined})(parent, transformation, attributes, input_args, converted)
    PT(parent, transformation, attributes, input_args, converted, AbstractPlot[])
end
plotsym(::Type{<:AbstractPlot{F}}) where F = Symbol(typeof(F).name.mt.name)

"""
    used_attributes(args...) = ()

function used to indicate what keyword args one wants to get passed in `convert_arguments`.
Usage:
```example
    struct MyType end
    used_attributes(::MyType) = (:attribute,)
    function convert_arguments(x::MyType; attribute = 1)
        ...
    end
    # attribute will get passed to convert_arguments
    # without keyword_verload, this wouldn't happen
    plot(MyType, attribute = 2)
    #You can also use the convenience macro, to overload convert_arguments in one step:
    @keywords convert_argumetns(x::MyType; attribute = 1)
        ...
    end
```
"""
used_attributes(PlotType, args...) = ()

"""
apply for return type
    (args...,)
"""
function apply_convert!(P, attributes::Attributes, x::Tuple)
    return (plottype(P, x...), x)
end

"""
apply for return type PlotSpec
"""
function apply_convert!(P, attributes::Attributes, x::PlotSpec{S}) where S
    args, kwargs = x.args, x.kwargs
    # Note that kw_args in the plot spec that are not part of the target plot type
    # will end in the "global plot" kw_args (rest)
    for (k, v) in pairs(kwargs)
        attributes[k] = v
    end
    return (plottype(P, S), args)
end

function seperate_tuple(args::Node{<: NTuple{N, Any}}) where N
    ntuple(N) do i
        lift(args) do x
            if i <= length(x)
                x[i]
            else
                error("You changed the number of arguments. This isn't allowed!")
            end
        end
    end
end

function (PlotType::Type{<: AbstractPlot{Typ}})(scene::SceneLike, attributes::Attributes, args) where Typ
    input = to_node.(args)
    argnodes = lift(input...) do args...
        convert_arguments(PlotType, args...)
    end
    PlotType(scene, attributes, input, argnodes)
end

function (PlotType::Type{<: AbstractPlot{Typ}})(scene::SceneLike, attributes::Attributes, input, args) where Typ
    # The argument type of the final plot object is the assumened to stay constant after
    # argument conversion. This might not always hold, but it simplifies
    # things quite a bit
    ArgTyp = typeof(to_value(args))
    # construct the fully qualified plot type, from the possible incomplete (abstract)
    # PlotType
    FinalType = Combined{Typ, ArgTyp}
    plot_attributes, scene_attributes = merged_get!(()-> default_theme(scene, FinalType), plotsym(FinalType), scene, attributes)
    trans = get(plot_attributes, :transformation, automatic)
    transformation = if to_value(trans) == automatic
        Transformation(scene)
    elseif isa(to_value(trans), Transformation)
        to_value(trans)
    else
        t = Transformation(scene)
        transform!(t, to_value(trans))
        t
    end

    replace_automatic!(plot_attributes, :model) do
        transformation.model
    end
    # create the plot, with the full attributes, the input signals, and the final signal nodes.
    plot_obj = FinalType(scene, transformation, plot_attributes, input, seperate_tuple(args))
    calculated_attributes!(plot_obj)
    plot_obj, scene_attributes
end

"""
    `plot_type(plot_args...)`

The default plot type for any argument is `lines`.
Any custom argument combination that has only one meaningful way to be plotted should overload this.
e.g.:
```example
    # make plot(rand(5, 5, 5)) plot as a volume
    plottype(x::Array{<: AbstractFlot, 3}) = Volume
```
"""
plottype(plot_args...) = Combined{Any, Tuple{typeof.(to_value.(plot_args))...}} # default to dispatch to type recipes!
plottype(::AbstractVector, ::AbstractVector) = Lines
plottype(::AbstractMatrix) = Heatmap
# If the Combined has no plot func, calculate them
plottype(::Type{<: Combined{Any}}, argvalues...) = plottype(argvalues...)
# If it has something more concrete than Any, use it directly
plottype(P::Type{<: Combined{T}}, argvalues...) where T = P

"""
    plottype(P1::Type{<: Combined{T1}}, P2::Type{<: Combined{T2}})

Chooses the more concrete plot type
```example
function convert_arguments(P::PlotFunc, args...)
    ptype = plottype(P, Lines)
    ...
end
"""
plottype(P1::Type{<: Combined{Any}}, P2::Type{<: Combined{T}}) where T = P2
plottype(P1::Type{<: Combined{T}}, P2::Type{<: Combined}) where T = P1

"""
Returns the Combined type that represents the signature of `args`.
"""
function Plot(args::Vararg{Any, N}) where N
    Combined{Any, <: Tuple{args...}}
end
Base.@pure function Plot(::Type{T}) where T
    Combined{Any, <: Tuple{T}}
end
Base.@pure function Plot(::Type{T1}, ::Type{T2}) where {T1, T2}
    Combined{Any, <: Tuple{T1, T2}}
end

# all the plotting functions that get a plot type
const PlotFunc = Union{Type{Any}, Type{<: AbstractPlot}}

plot(P::PlotFunc, args...; kw_attributes...) = plot!(Scene(), P, Attributes(kw_attributes), args...)
plot!(P::PlotFunc, args...; kw_attributes...) = plot!(current_scene(), P, Attributes(kw_attributes), args...)
plot(scene::SceneLike, P::PlotFunc, args...; kw_attributes...) = plot!(Scene(scene), P, Attributes(kw_attributes), args...)
plot!(scene::SceneLike, P::PlotFunc, args...; kw_attributes...) = plot!(scene, P, Attributes(kw_attributes), args...)

plot(scene::SceneLike, P::PlotFunc, attributes::Attributes, args...; kw_attributes...) = plot!(Scene(scene), P, merge!(Attributes(kw_attributes), attributes), args...)
plot!(P::PlotFunc, attributes::Attributes, args...; kw_attributes...) = plot!(current_scene(), P, merge!(Attributes(kw_attributes), attributes), args...)
plot(P::PlotFunc, attributes::Attributes, args...; kw_attributes...) = plot!(Scene(), P, merge!(Attributes(kw_attributes), attributes), args...)

# Overload remaining functions
eval(default_plot_signatures(:plot, :plot!, :Any))

# plots to scene

plotfunc(::Combined{F}) where F = F

"""
Main plotting signatures that plot/plot! route to if no Plot Type is given
"""
function plot!(scene::SceneLike, P::PlotFunc, attributes::Attributes, args...; kw_attributes...)
    attributes = merge!(Attributes(kw_attributes), attributes)
    argvalues = to_value.(args)
    PreType = plottype(P, argvalues...)
    # plottype will lose the argument types, so we just extract the plot func
    # type and recreate the type with the argument type
    PreType = Combined{plotfunc(PreType), typeof(argvalues)}
    convert_keys = intersect(used_attributes(PreType, argvalues...), keys(attributes))
    kw_signal = if isempty(convert_keys) # lift(f) isn't supported so we need to catch the empty case
        Node(())
    else
        lift((args...)-> Pair.(convert_keys, args), getindex.(attributes, convert_keys)...) # make them one tuple to easier pass through
    end
    # call convert_arguments for a first time to get things started
    converted = convert_arguments(PreType, argvalues...; kw_signal[]...)
    # convert_arguments can return different things depending on the recipe type
    # apply_conversion deals with that!

    FinalType, argsconverted = apply_convert!(PreType, attributes, converted)
    converted_node = Node(argsconverted)
    input_nodes =  to_node.(args)
    onany(kw_signal, lift(tuple, input_nodes...)) do kwargs, args
        # do the argument conversion inside a lift
        result = convert_arguments(FinalType, args...; kwargs...)
        finaltype, argsconverted = apply_convert!(FinalType, attributes, result)
        if finaltype != FinalType
            error("Plot type changed from $FinalType to $finaltype after conversion.
                Changing the plot type based on values in convert_arguments is not allowed"
            )
        end
        converted_node[] = argsconverted
    end
    plot!(scene, FinalType, attributes, input_nodes, converted_node)
end

plot!(p::Combined) = _plot!(p)

_plot!(p::Atomic{T}) where T = p

function _plot!(p::Combined{Any, T}) where T
    args = (T.parameters...,)
    typed_args = join(string.("::", args), ", ")
    error("Plotting for the arguments ($typed_args) not defined. If you want to support those arguments, overload plot!(plot::Plot$((args...,)))")
end
function _plot!(p::Combined{X, T}) where {X, T}
    args = (T.parameters...,)
    typed_args = join(string.("::", args), ", ")
    error("Plotting for the arguments ($typed_args) not defined for $X. If you want to support those arguments, overload plot!(plot::$X{ <: $T})")
end

function plot!(scene::SceneLike, ::Type{PlotType}, attributes::Attributes, input::NTuple{N, Node}, args::Node) where {N, PlotType <: AbstractPlot}
    # create "empty" plot type - empty meaning containing no plots, just attributes + arguments
    plot_object, scene_attributes = PlotType(scene, attributes, input, args)

    nattributes, rest = merge_attributes!(scene_attributes, theme(scene))

    # TODO warn about rest - should be unused arguments!
    empty!(scene.attributes)
    # transfer the merged attributes from theme and user defined to the scene
    merge!(scene.attributes, nattributes)
    for (at1, at2) in mutual_exclusive_attributes(PlotType)
        #nothing here to get around defaults in GLVisualize
        haskey(attributes, at1) && haskey(attributes, at2) && error("$at1 conflicts with $at2, please specify only one.")
        if haskey(attributes, at1) && haskey(plot_object.attributes, at2)
            plot_object.attributes[at2] = nothing
        elseif haskey(attributes, at2) && haskey(plot_object.attributes, at1)
            plot_object.attributes[at1] = nothing
        end
    end
    # call user defined recipe overload to fill the plot type
    plot!(plot_object)

    push!(scene.plots, plot_object)

    scene[:raw][] || update_limits!(scene)
    (!scene[:raw][] || scene[:camera][] != automatic) && setup_camera!(scene)
    scene[:raw][] || add_axis!(scene, rest)
    # ! ∘ isaxis --> (x)-> !isaxis(x)
    # move axis to front, so that scene[end] gives back the last plot and not the axis!
    if !isempty(scene.plots) && isaxis(last(scene.plots))
        axis = pop!(scene.plots)
        pushfirst!(scene.plots, axis)
    end
    # call the assembly recipe, that also adds this to the scene
    # kw_args not consumed by PlotType will be passed forward to plot! as non_plot_kwargs
    #plot!(scene, plot_object, scene_attributes)
    scene
end

function plot!(scene::Combined, ::Type{PlotType}, attributes::Attributes, args...) where PlotType <: AbstractPlot
    # create "empty" plot type - empty meaning containing no plots, just attributes + arguments
    plot_object, scene_attributes = PlotType(scene, attributes, args)
    # call user defined recipe overload to fill the plot type
    plot!(plot_object)
    push!(scene.plots, plot_object)
    scene
end

function scale_scene!(scene)
    if is2d(scene)
        area = pixelarea(scene)[]
        lims = limits(scene)[]
        # not really sure how to scale 3D scenes in a reasonable way
        mini, maxi = minimum(lims), maximum(lims)
        l = ((mini[1], maxi[1]), (mini[2], maxi[2]))
        xyzfit = fit_ratio(area, l)
        s = to_ndim(Vec3f0, xyzfit, 1f0)
        scale!(scene, s)
    end
    return scene
end

function setup_camera!(scene::Scene)
    if scene[:camera][] == automatic
        cam = cameracontrols(scene)
        if cam == EmptyCamera()
            if is2d(scene)
                #@info("setting camera to 2D")
                cam2d!(scene)
            else
                #@info("setting camera to 3D")
                cam3d!(scene)
            end
        end
    elseif scene[:camera][] in (cam2d!, cam3d!, campixel!)
        scene[:camera][](scene)
    else
        error("Unrecogniced `camera` attribute type: $(typeof(scene[:camera][])). Use automatic, cam2d! or cam3d!")
    end
    scene
end

function find_in_plots(scene::Scene, key::Symbol)
    # TODO findfirst is a bit flaky... maybe merge multiple ranges + tick labels?!
    idx = findfirst(scene.plots) do plot
        !isaxis(plot) && haskey(plot, key) && plot[key][] !== automatic
    end
    if idx !== nothing
        scene.plots[idx][key]
    else
        automatic
    end
end

function add_axis!(scene::Scene, attributes = Attributes())
    show_axis = scene[:show_axis][]
    show_axis isa Bool || error("show_axis needs to be a bool")
    axistype = if scene[:axis_type][] == automatic
        is2d(scene) ? axis2d! : axis3d!
    elseif scene[:axis_type][] in (axis2d!, axis3d!)
        scene[:axis_type][]
    else
        error("Unrecogniced `axis_type` attribute type: $(typeof(scene[:axis_type][])). Use automatic, axis2d! or axis3d!")
    end

    if show_axis && !(any(isaxis, plots(scene)))
        key = axistype == axis2d! ? :axis2d : :axis3d
        axis_attributes = scene[key]
        ranges = get(attributes, :tickranges) do
            find_in_plots(scene, :tickranges)
        end
        labels = get(attributes, :ticklabels) do
            find_in_plots(scene, :ticklabels)
        end
        axistype(
            scene, axis_attributes, limits(scene),
            ticks = (ranges = ranges, labels = labels)
        )
    end
    scene
end

function add_labels!(scene::Scene)
    if plot_attributes[:show_legend][] && haskey(p.attributes, :colormap)
        legend_attributes = plot_attributes[:legend][]
        colorlegend(scene, p.attributes[:colormap], p.attributes[:colorrange], legend_attributes)
    end
    scene
end

update_limits!(scene::Scene) = update_limits!(scene, scene[:limits][], scene[:padding][])

function update_limits!(scene::Scene, limits::Automatic, padding)
    # for when scene is empty
    dlimits = data_limits(scene)
    tlims = (minimum(dlimits), maximum(dlimits))
    if !all(x-> all(isfinite, x), tlims)
        @warn "limits of scene contain non finite values: $(tlims[1]) .. $(tlims[2])"
        mini = map(x-> ifelse(isfinite(x), x, 0.0), tlims[1])
        maxi = Vec3f0(ntuple(3) do i
            x = tlims[2][i]
            ifelse(isfinite(x), x, tlims[1][i] + 1f0)
        end)
        tlims = (mini, maxi)
    end
    new_widths = Vec3f0(ntuple(3) do i
        a = tlims[1][i]; b = tlims[2][i]
        w = b - a
        # check for widths == 0.0... 3rd dimension is allowed to be 0 though.
        # TODO maybe we should allow any one dimension to be 0, and then use the other 2 as 2D
        with0 = (i != 3) && (w ≈ 0.0)
        with0 && @warn "Founds 0 width in scene limits: $(tlims[1]) .. $(tlims[2])"
        ifelse(with0, 1f0, w)
    end)
    update_limits!(scene, FRect3D(tlims[1], new_widths), padding)
end

function update_limits!(scene::Scene, new_limits::HyperRectangle, padding = Vec3f0(0))
    lims = FRect3D(new_limits)
    lim_w = widths(lims)
    # use the smallest widths for scaling, to have a consistently wide padding for all sides
    minw = if lim_w[3] ≈ 0.0
        m = min(lim_w[1], lim_w[2])
        Vec3f0(m, m, 0.0)
    else
        Vec3f0(minimum(lim_w))
    end
    padd_abs = minw .* to_ndim(Vec3f0, padding, 0.0)
    limits(scene)[] = FRect3D(minimum(lims) .- padd_abs, lim_w .+  2padd_abs)
    scene
end


text:
"""
    Bingham(M; params...)

The Bingham distribution on the manifold ``M``.

Accepted manifolds are `Sphere`, `ProjectiveSpace`, `Stiefel`, `Grassmann`, `Rotations`
and `SpecialOrthogonal`.

For manifolds with matrix points, this is also called the Matrix Bingham distribution.

# Constructors

    Bingham(M; A)

For a manifold ``M ⊂ 𝔽^{n × k}``, construct the Bingham distribution parameterized by
some positive definite matrix ``A ∈ 𝔽^{n × n}``.

The density function with respect to the normalized [`Hausdorff`](@ref) measure on ``M`` is

```math
p(x | A) = \\frac{\\exp(\\Re⟨x, Ax⟩)}{_1 F_1(\\frac{k}{2}, \\frac{n}{2}; A)},
```

where ``⟨⋅,⋅⟩`` is the Frobenius inner product, and ``_1 F_1(a, b; A)``
is a hypergeometric function with matrix argument ``A``.

Note that ``p(x | A + α I) = p(x | A)`` for all scalars ``α ∈ 𝔽``.
Hence, ``A`` can not be uniquely identified from draws from the Bingham distribution.
"""
struct Bingham{M,N,T} <: ParameterizedMeasure{N}
    manifold::M
    par::NamedTuple{N,T}
end
Bingham(M; params...) = Bingham(M, NamedTuple(params))

function Base.show(io::IO, mime::MIME"text/plain", μ::Bingham)
    return show_manifold_measure(io, mime, μ)
end

Manifolds.base_manifold(d::Bingham) = getfield(d, :manifold)

function MeasureTheory.basemeasure(μ::Bingham)
    return normalize(Hausdorff(base_manifold(μ)))
end

function MeasureTheory.logdensity(d::Bingham{M,(:A,)}, x) where {M}
    n = size(x, 1)
    k = size(x, 2)
    A = d.A
    return real(dot(x, A, x)) - logpFq((k//2,), (n//2,), A)
end

ext: jl


text:
# Integer
for op = [:add, :sub, :max, :min, :and, :or, :xor]
    fn = Symbol("scatter_$(op)!")
    atm_op = Symbol("atomic_$(op)!")
    @eval function $fn(ys::CuMatrix{T}, us::CuArray{T}, xs::CuArray{Int}) where {T<:Integer}
        function kernel!(ys, us, xs)
            li = threadIdx().y + (blockIdx().y - 1) * blockDim().y
            i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

            @inbounds if li <= length(xs) && i <= size(ys, 1)
                ind = CartesianIndices(xs)[li]
                j = Base._to_linear_index(ys, i, xs[li])
                CUDA.$atm_op(pointer(ys, j), us[i, ind])
            end

            return
        end

        thread_x = min(MAX_THREADS, size(ys, 1))
        thread_y = min(MAX_THREADS ÷ thread_x, length(xs))
        threads = (thread_x, thread_y)
        blocks = ceil.(Int, (size(ys, 1), length(xs)) ./ threads)
        @cuda blocks=blocks threads=threads kernel!(ys, us, xs)
        return ys
    end

    @eval function $fn(ys::CuArray{T}, us::CuArray{T}, xs::CuArray{<:Tuple}) where {T<:Integer}
        function kernel!(ys, us, xs)
            li = threadIdx().y + (blockIdx().y - 1) * blockDim().y
            i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

            @inbounds if li <= length(xs) && i <= size(ys, 1)
                ind = CartesianIndices(xs)[li]
                j = Base._to_linear_index(ys, i, xs[li]...)
                CUDA.$atm_op(pointer(ys, j), us[i, ind])
            end

            return
        end

        thread_x = min(MAX_THREADS, size(ys, 1))
        thread_y = min(MAX_THREADS ÷ thread_x, length(xs))
        threads = (thread_x, thread_y)
        blocks = ceil.(Int, (size(ys, 1), length(xs)) ./ threads)
        @cuda blocks=blocks threads=threads kernel!(ys, us, xs)
        return ys
    end
end

# Floating point
for op = [:add, :sub, :mul, :div, :max, :min]
    fn = Symbol("scatter_$(op)!")
    atm_op = Symbol("atomic_$(op)!")
    @eval function $fn(ys::CuMatrix{T}, us::CuArray{T}, xs::CuArray{Int}) where {T<:AbstractFloat}
        function kernel!(ys::CuDeviceArray{T}, us::CuDeviceArray{T}, xs)
            i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
            j = threadIdx().y + (blockIdx().y - 1) * blockDim().y

            @inbounds if i <= size(ys, 1) && j <= length(xs)
                ind = CartesianIndices(xs)[j]
                k = Base._to_linear_index(ys, i, xs[j])
                CUDA.$atm_op(pointer(ys, k), us[i, ind])
            end

            return
        end

        thread_i = min(MAX_THREADS, size(ys, 1))
        thread_j = min(MAX_THREADS ÷ thread_i, length(xs))
        threads = (thread_i, thread_j)
        blocks = ceil.(Int, (size(ys, 1), length(xs)) ./ threads)
        @cuda blocks=blocks threads=threads kernel!(ys, us, xs)
        return ys
    end

    @eval function $fn(ys::CuArray{T}, us::CuArray{T}, xs::CuArray{<:Tuple}) where {T<:AbstractFloat}
        function kernel!(ys::CuDeviceArray{T}, us::CuDeviceArray{T}, xs)
            i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
            j = threadIdx().y + (blockIdx().y - 1) * blockDim().y

            @inbounds if i <= size(ys, 1) && j <= length(xs)
                ind = CartesianIndices(xs)[j]
                k = Base._to_linear_index(ys, i, xs[j]...)
                CUDA.$atm_op(pointer(ys, k), us[i, ind])
            end

            return
        end

        thread_i = min(MAX_THREADS, size(ys, 1))
        thread_j = min(MAX_THREADS ÷ thread_i, length(xs))
        threads = (thread_i, thread_j)
        blocks = ceil.(Int, (size(ys, 1), length(xs)) ./ threads)
        @cuda blocks=blocks threads=threads kernel!(ys, us, xs)
        return ys
    end
end

function scatter_mean!(ys::CuMatrix{T}, us::CuArray{T}, xs::CuArray{<:IntOrTuple}) where {T<:AbstractFloat}
    yt = CUDA.zero(ys)
    ot = CUDA.zero(ys)
    os = CUDA.one.(us)
    scatter_add!(ot, os, xs)
    scatter_add!(yt, us, xs)
    ys .+= save_div.(yt, ot)
    return ys
end


text:
using JuMPContainerBenchmarks
using Test

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

max_stars_repo_name: SonntagK/JuMPContainerBenchmarks

max_issues_repo_name: SonntagK/JuMPContainerBenchmarks

max_forks_repo_name: SonntagK/JuMPContainerBenchmarks


text:
VERSION < v"0.1"

"""
Data filters:    `filter_array`, `filter_array_2`, `filter_array_time`, `getsmoothdata`, `getsmoothdata_nospace`
Data structures:  `ape_budget`, `cat_ape_budget`, `cutborders!`, `surf_quantities`, `cyclone_comp_timemean`, `Composite_Cyclone`, `Composite_Cyclone_v2`, `Composite_Cyclone_v3`
Methods:    `cyclonecompositer`, `shifter`, `smoothfilter`, `cyclonecompositer_v2`, `cyclonecompositer_v3`, `timemean_nofalseframe`, `removefalseframes`, `getapebudget`, `buoyancybudget`

"""
module AvailablePotentialEnergyFramework

using DataStructures: OrderedDict
using ImageFiltering: imfilter, imfilter!, centered, kernelfactors, mapwindow, mapwindow!, Kernel, Inner
using Images: findlocalminima, findlocalmaxima
using ImageSegmentation: SegmentedImage, segment_labels, region_adjacency_graph, seeded_region_growing, labels_map
using Interpolations: LinearInterpolation, interpolate, Gridded, Linear
using NCDatasets: Dataset, variable, defVar
using OffsetArrays: OffsetArray
using SparseArrays: SparseMatrixCSC
using Statistics: mean, median!, mean!
using Unitful: @u_str, unit, ustrip, Quantity
using LoopVectorization: @turbo, @tturbo
using Tullio: @tullio

include("apehelperfunctions.jl")
include("compositehelperfunctions.jl")
include("apebudgets.jl")
include("physicalconstants.jl")
include("physicsfunctions.jl")
include("useful_diagnostics.jl")
#include("ape_computation_from_julia_output.jl") #testing purposes
include("arrayfiltering.jl")
include("datamanagement.jl")
export
# Filters
    filter_array!,
    filter_array_2!,
    filter_array_time,
    filter_array,
    getsmoothdata!,
    getsmoothdata_nospace,
    #Data structures
    cat_ape_budget,
    cutborders!,
    #Methods
#    findlocalmaxima,
    cyclonecompositer,
    shifter,
    shifter!,
    smooth_and_mask,
    smooth_and_mask!,
    getapebudget,
    buoyancybudget,
    add_allcyclones!,
    averageallindistance,
    detect_cyclones,
    detect_cyclones!,
    get_diabatic_as_residual_buoyancy,
    run_distributed_test,
    ##Constants
    R,
    Dryair,
    Watervapor,
    Liquidwater,
    epsilon,
    g,
    #Physics functions
    compute_N2,
    compute_mse,
    get_mse_budget,
    get_vorticity,
    get_okubo_weiss,
    get_divergence,
    get_saturation_vapor_pressure,
    get_partial_vapor_pressure,
    get_mixing_ratio,
    get_virtual_temperature,
    get_lifted_condensation_level,
    get_specific_entropy,
    get_potential_temperature,
    get_virtual_temperature,
    mixing_ratio_to_specific_humidity,
    specific_humidity_to_mixing_ratio,
    get_buoyancy_of_lifted_parcel,
    surface_sensible_heat_flux_to_buoyancy,
    surface_latent_heat_flux_to_buoyancy,
    get_buoyancy,
    radiative_heating_rate_to_buoyancy,
    get_density_ideal_gas,
    #math functions
    integrate_vertically,
    average_precipitation_per_pw_bin,
    average_precipitation_per_pw_bin_dayang,
    #Files and datamanagement
    smooth_vars_and_write_to_netcdf!,
    create_APE_netcdf,
    set_netcdf_var!
end


text:
"""
    EGPpowerfit()

Estimate the first extended GP model of Naveau et al. (2016) parameters by maximum likelihood.
"""
function EGPpowerfit(data::Array{<:Real,1};
    initialvalues::Vector{<:Real}=Float64[],
    censoring::Real=0,
    covariate::Array{<:Real,1}=Float64[],  # Vecteur de la variable explicative
    scalecov::Bool=true, # variable explicative sur σ ?
    lowertailcov::Bool=false,  # variable explicative sur κ ?
    uppertailcov::Bool=false)  # variable explicative sur ξ ?)

    if !isempty(covariate)
        return EGPnonstatpowerfit(data, covariate, initialvalues=initialvalues, scalecov=scalecov, lowertailcov=lowertailcov, uppertailcov=uppertailcov, censoring=censoring)
    else
        if isempty(initialvalues)
            σ₀, ξ₀, κ₀ = [1, 0.15, 1]
            initialvalues = [σ₀, ξ₀, κ₀]
        end

        function loglike(σ::Real, ξ::Real, κ::Real)
            if κ <= 0 || σ <= 0
                return -Inf
            else
                pd = EGPpower(σ, ξ, κ)

                r_data = data[data .>= censoring]
                l_data = fill(censoring, count(data .< censoring))

                ll = sum(logcdf.(pd, l_data)) + sum(logpdf.(pd, r_data))

                return ll
            end
        end

        fobj(θ) = -loglike(θ...)

        res = optimize(fobj, initialvalues)

        if Optim.converged(res)
            σ̂, ξ̂, κ̂ = [Optim.minimizer(res)[1], Optim.minimizer(res)[2], Optim.minimizer(res)[3]]
        else
            @warn "The maximum likelihood algorithm did not find a solution. Maybe try with different initial values or with another method. The returned values are the initial values."
            σ̂, ξ̂, κ̂  = [initialvalues[1], initialvalues[2], initialvalues[3]]
        end

        fittedmodel = EGPpower(σ̂, ξ̂, κ̂)

        return fittedmodel
    end
end

function EGPpowerfit(data::TimeArray;
    initialvalues::Vector{<:Real}=Float64[],
    censoring::Real=0,
    covariate::TimeArray=TimeArray(DateTime[], Float64[]),  # Vecteur de la variable explicative
    scalecov::Bool=true,
    lowertailcov::Bool=false,  # variable explicative sur κ ?
    uppertailcov::Bool=false)  # variable explicative sur ξ ?)

    if isempty(covariate)
        return EGPpowerfit(values(data), initialvalues=initialvalues, scalecov=scalecov, lowertailcov=lowertailcov, uppertailcov=uppertailcov, censoring=censoring)
    else
        return EGPpowerfit(values(data), covariate=values(covariate), initialvalues=initialvalues,
        scalecov=scalecov, lowertailcov=lowertailcov, uppertailcov=uppertailcov, censoring=censoring)
    end
end

"""
    EGPpowermixfit()

Estimate the second extended GP model of Naveau et al. (2016) parameters by maximum likelihood.
"""
function EGPpowermixfit(data::Array{<:Real,1}; initialvalues::Vector{<:Real}=Float64[], censoring::Real=0)

    #inner_optimizer = NelderMead()
    #lower = [0, -Inf, 0, 0, 0]
    #upper = [Inf, Inf, Inf, Inf, 1]

    if isempty(initialvalues)
        σ₀, ξ₀, κ₁₀, κ₂₀, p₀ = [1, .15, 1, 2, .4]  # à remplacer par une fonction (getinitialvalue)
        initialvalues = [σ₀, ξ₀, κ₁₀, κ₂₀, p₀]
    end

    function loglike(σ::Real, ξ::Real, κ₁::Real, κ₂::Real, p::Real)
        if p < 0 || p > 1 || κ₁ <= 0 || κ₂ <= 0 || σ <= 0 || κ₂ < κ₁
            return -Inf
        else
            pd = EGPpowermix(σ, ξ, κ₁, κ₂, p)

            r_data = data[data .>= censoring]
            l_data = fill(censoring, count(data .< censoring))

            ll = sum(logcdf.(pd, l_data)) + sum(logpdf.(pd, r_data))

            return ll
        end
    end

    fobj(θ) = -loglike(θ...)

    res = optimize(fobj, initialvalues)

    if Optim.converged(res)
        σ̂, ξ̂, κ̂₁, κ̂₂, p̂ = [Optim.minimizer(res)...]
    else
        @warn "The maximum likelihood algorithm did not find a solution. Maybe try with different initial values or with another method. The returned values are the initial values."
        σ̂, ξ̂, κ̂₁, κ̂₂, p̂  = [initialvalues...]
    end

    fittedmodel = EGPpowermix(σ̂, ξ̂, κ̂₁, κ̂₂, p̂)

    return fittedmodel
end

function EGPpowermixfit(data::TimeArray; initialvalues::Vector{<:Real}=Float64[], censoring::Real=0)
    return EGPpowermixfit(values(data), initialvalues=initialvalues, censoring=censoring)
end

"""
    EGPbetafit()

Estimate the third extended GP model of Naveau et al. (2016) parameters by maximum likelihood.
"""
function EGPbetafit(data::Array{<:Real,1}; initialvalues::Vector{<:Real}=Float64[], censoring::Real=0)

    if isempty(initialvalues)
        σ₀, ξ₀, δ₀ = [5, 0.15, 2]
        initialvalues = [σ₀, ξ₀, δ₀]
    end

    function loglike(σ::Real, ξ::Real, δ::Real)
        if δ <= 0 || δ > 100 || σ <= 0 || ξ > 0.99
            return -Inf
        else
            pd = EGPbeta(σ, ξ, δ)

            r_data = data[data .>= censoring]
            l_data = fill(censoring, count(data .< censoring))

            ll = sum(logcdf.(pd, l_data)) + sum(logpdf.(pd, r_data))

            return ll
        end
    end

    fobj(θ) = -loglike(θ...)

    res = optimize(fobj, initialvalues)

    if Optim.converged(res)
        σ̂, ξ̂, δ̂ = [Optim.minimizer(res)[1], Optim.minimizer(res)[2], Optim.minimizer(res)[3]]
    else
        @warn "The maximum likelihood algorithm did not find a solution. Maybe try with different initial values or with another method. The returned values are the initial values."
        σ̂, ξ̂, δ̂  = [initialvalues[1], initialvalues[2], initialvalues[3]]
    end

    fittedmodel = EGPbeta(σ̂, ξ̂, δ̂)

    return fittedmodel
end

function EGPbetafit(data::TimeArray; initialvalues::Vector{<:Real}=Float64[], censoring::Real=0)
    return EGPbetafit(values(data), initialvalues=initialvalues, censoring=censoring)
end

"""
    EGPbetapowerfit()

Estimate the fourth extended GP model of Naveau et al. (2016) parameters by maximum likelihood.
"""
function EGPbetapowerfit(data::Array{<:Real,1}; initialvalues::Vector{<:Real}=Float64[], censoring::Real=0)

    if isempty(initialvalues)
        σ₀, ξ₀, δ₀, κ₀ = [1, 0.15, 1, 1]
        initialvalues = [σ₀, ξ₀, δ₀, κ₀]
    end

    function loglike(σ::Real, ξ::Real, δ::Real, κ::Real)
        if κ <= 0 || δ <= 0 || δ > 100 || σ <= 0
            return -Inf
        else
            pd = EGPbetapower(σ, ξ, δ, κ)

            r_data = data[data .>= censoring]
            l_data = fill(censoring, count(data .< censoring))

            ll = sum(logcdf.(pd, l_data)) + sum(logpdf.(pd, r_data))

            return ll
        end
    end

    fobj(θ) = -loglike(θ...)

    res = optimize(fobj, initialvalues)

    if Optim.converged(res)
        σ̂, ξ̂, δ̂, κ̂ = [Optim.minimizer(res)[1], Optim.minimizer(res)[2], Optim.minimizer(res)[3], Optim.minimizer(res)[4]]
    else
        @warn "The maximum likelihood algorithm did not find a solution. Maybe try with different initial values or with another method. The returned values are the initial values."
        σ̂, ξ̂, δ̂, κ̂  = [initialvalues[1], initialvalues[2], initialvalues[3], initialvalues[4]]
    end

    fittedmodel = EGPbetapower(σ̂, ξ̂, δ̂, κ̂)

    return fittedmodel
end

function EGPbetapowerfit(data::TimeArray; initialvalues::Vector{<:Real}=Float64[], censoring::Real=0)
    return EGPbetapowerfit(values(data), initialvalues=initialvalues, censoring=censoring)
end

"""
    EGPnormalfit()

Estimate the extended GP model of Gamet and Jalbert (2020) parameters by maximum likelihood.
"""
function EGPnormalfit(data::Array{<:Real,1}; initialvalues::Vector{<:Real}=Float64[])

    if isempty(initialvalues)
        logσ₀, ξ₀, logκ₀ = [0, 0.15, 0]  # à remplacer par une fonction (getinitialvalue)
        initialvalues = [logσ₀, ξ₀, logκ₀]
    end

    function loglike(logσ::Real, ξ::Real, logκ::Real)
        pd = EGPnormal(exp(logσ), ξ, exp(logκ))
        ll = sum(logpdf.(pd, data))
        return ll
    end

    fobj(θ) = -loglike(θ...)

    res = optimize(fobj, initialvalues)

    if Optim.converged(res)
        σ̂, ξ̂, κ̂ = [exp(Optim.minimizer(res)[1]), Optim.minimizer(res)[2], exp(Optim.minimizer(res)[3])]
    else
        @warn "The maximum likelihood algorithm did not find a solution. Maybe try with different initial values or with another method. The returned values are the initial values."
        σ̂, ξ̂, κ̂  = [exp(initialvalues[1]), initialvalues[2], exp(initialvalues[3])]
    end

    fittedmodel = EGPnormal(σ̂, ξ̂, κ̂)

    return fittedmodel
end

function EGPnormalfit(data::TimeArray; initialvalues::Vector{<:Real}=Float64[])
    return EGPnormalfit(values(data), initialvalues=initialvalues)
end


text:
using DataStructures, Test
import DataStructures: SparseIntSet

@testset "SparseIntSet" begin
    @testset "Construction, collect" begin
        data_in = (1,5,100)
        s = SparseIntSet(data_in)
        data_out = collect(s)
        @test all(map(d->in(d,data_out), data_in))
        @test length(data_out) == length(data_in)
    end

    @testset "eltype, empty" begin
        @test eltype(SparseIntSet()) == Int
        @test eltype(typeof(SparseIntSet())) == Int
        @test isequal(empty(SparseIntSet([1,2,3])), SparseIntSet())
    end

    @testset "Core Functionality" begin
        s = SparseIntSet([1,2,10,20,200,300,1000,10000,10002])
        @test last(s) == 10002
        @test first(s) == 1
        @test length(s) == 9
        @test pop!(s) == 10002
        @test length(s) == 8
        @test popfirst!(s) == 1
        @test length(s) == 7
        @test !in(1,s)
        @test !in(10002,s)
        @test in(10000,s)
        @test_throws ArgumentError first(SparseIntSet())
        @test_throws ArgumentError last(SparseIntSet())
        t = copy(s)
        s = SparseIntSet()
        push!(s, 1, 2, 100)
        @test 1 in s
        @test !(3 in s)
        @test 2 in s
        @test 100 in s
        @test !(101 in s)
        @test !(1000 in s)
        @test first(s) == 1
        @test last(s) == 100
        @test s == SparseIntSet([1, 2, 100])
        push!(s, 1000)
        @test [i for i in s] == [1, 2, 100, 1000]
        @test pop!(s) == 1000
        @test s == SparseIntSet([1, 2, 100])
        push!(s, 5000)
        push!(s, 2000)
        pop!(s, 5000)
        @test s.reverse[end] === DataStructures.NULL_INT_PAGE
        b = 1:1000
        s = SparseIntSet(b)
        @test collect(s) == collect(b)
        @test length(s) == length(b)
        @test pop!(s, 100) == 100
        @test_throws BoundsError pop!(s, 100)
        @test pop!(s, 100, 1) == 1
        @test pop!(s, 99, 1) == 99
        @test !in(500000, s)
        @test !in(99, s)
    end

    @testset "setdiff / symdiff" begin
        @test setdiff(SparseIntSet([1, 2, 3, 4]), SparseIntSet([2, 4, 5, 6])) == SparseIntSet([1, 3])
    end

    @testset "setdiff!" begin
        s2 = SparseIntSet([1, 2, 3, 4])
        setdiff!(s2, SparseIntSet([2, 4, 5, 6]))

        @test s2 == SparseIntSet([1, 3])
    end

    @testset "issue #7851" begin
        @test_throws DomainError SparseIntSet(-1)
        @test !(-1 in SparseIntSet(1:10))
    end
    @testset "Copy, copy!, empty" begin
        s1 = SparseIntSet([1,2,3])
        s2 = empty(s1)
        push!(s2, 10000)
        @test !in(10000, s1)
        copy!(s2, s1)
        @test !in(10000, s2)
        push!(s2, 10000)
        @test !in(10000, s1)
        s3 = copy(s2)
        push!(s3, 1000)
        @test !in(1000, s2)
        pop!(s3, 1000)
        pop!(s2, 10000)
        @test in(10000, s3)
        pop!(s3, 10000)
        @test s3 == s2 == s1
        @test collect(s3) == collect(s2) == [1,2,3]

    end

    @testset "Push, union" begin
        # Push, union
        s1 = SparseIntSet()
        @test_throws DomainError push!(s1, -1)
        push!(s1, 1, 10, 100, 1000)
        @test collect(s1) == [1, 10, 100, 1000]
        push!(s1, 606)
        @test collect(s1) == [1, 10, 100, 1000, 606]
        s2 = SparseIntSet()
        @test s2 === union!(s2, s1)
        s3 = SparseIntSet([1, 10, 100])
        union!(s3, [1, 606, 1000])
        s4 = union(SparseIntSet([1, 100, 1000]), SparseIntSet([10, 100, 606]))
        @test s1 == s2 == s3 == s4
    end

    @testset "pop!, delete!" begin
        s = SparseIntSet(1:2:10)
        @test pop!(s, 1) == 1
        @test !(1 in s)
        @test_throws BoundsError pop!(s, 1)
        @test_throws ArgumentError pop!(s, -1)
        @test_throws ArgumentError pop!(s, -1, 1)
        @test pop!(s, 1, 0) == 0
        is = copy(s.packed)
        for i in is; pop!(s, i); end
        @test isempty(s)
        push!(s, 1:2:10...)
        @test pop!(s) == 9
        @test pop!(s) == 7
        @test popfirst!(s) == 1
        @test popfirst!(s) == 5
        @test collect(s) == [3]
        empty!(s)
        @test isempty(s)
    end

    @testset "Intersect" begin
        @test isempty(intersect(SparseIntSet()))
        @test isempty(intersect(SparseIntSet(1:10), SparseIntSet()))
        @test isempty(intersect(SparseIntSet(), SparseIntSet(1:10)))

        @test intersect(SparseIntSet([1,2,3])) == SparseIntSet([1,2,3])

        @test intersect(SparseIntSet(1:7), SparseIntSet(3:10)) ==
              intersect(SparseIntSet(3:10), SparseIntSet(1:7)) == SparseIntSet(3:7)

        @test intersect!(SparseIntSet(1:10), SparseIntSet(1:4), 1:5, [1,2,10]) == SparseIntSet(1:2)
    end

    @testset "Setdiff" begin
        s1 = SparseIntSet(1:100)
        setdiff!(s1, SparseIntSet(1:2:100))
        s2 = setdiff(SparseIntSet(1:100), SparseIntSet(1:2:100))
        @test s1 == s2 == SparseIntSet(2:2:100)

        s1 = SparseIntSet(1:10)
        s2 = SparseIntSet([1:2; 6:100])
        @test setdiff(s1, s2) == setdiff(s1, [1:2; 6:100]) == SparseIntSet(3:5)
    end

    @testset "Subsets, equality" begin
        @test SparseIntSet(2:2:10) < SparseIntSet(1:10)
        @test !(SparseIntSet(2:2:10) < SparseIntSet(2:2:10))
        @test SparseIntSet(2:2:10) <= SparseIntSet(2:10)
        @test SparseIntSet(2:2:10) <= SparseIntSet(2:2:10)
    end
    @testset "zip" begin
    	a = SparseIntSet([1,2,3,5, 6, 9, 12, 24])
    	b = SparseIntSet([6, 12, 24, 1000, 2000, 3000])
    	c = SparseIntSet(2:2:100)
    	d = SparseIntSet(6:3:100)
    	e = SparseIntSet((6, 12))
    	s1 = 0
    	it = zip(a, b, c, d, e)
    	for (ia, ib, ic, id, ie) in it
	    	s1 += a.packed[ia] + b.packed[ib] + c.packed[ic] + d.packed[id] + e.packed[ie]
	    end
	    @test s1 == 5*(6+12)
    	s1 = 0
    	it = zip(a, b, c, d, exclude=(e,))
    	for (ia, ib, ic, id) in it
	    	s1 += a.packed[ia] + b.packed[ib] + c.packed[ic] + d.packed[id]
	    end
	    @test s1 == 4*24
	end
        if VERSION >= v"1.1"
            @test zip() isa Iterators.Zip  # issue 621
        end

end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@doc raw"""
    CircleGroup <: GroupManifold{Circle{ℂ},MultiplicationOperation}

The circle group is the complex circle ([`Circle(ℂ)`](@ref)) equipped with
the group operation of complex multiplication ([`MultiplicationOperation`](@ref)).
"""
const CircleGroup = GroupManifold{ℂ,Circle{ℂ},MultiplicationOperation}

CircleGroup() = GroupManifold(Circle{ℂ}(), MultiplicationOperation())

Base.show(io::IO, ::CircleGroup) = print(io, "CircleGroup()")

invariant_metric_dispatch(::CircleGroup, ::ActionDirection) = Val(true)

default_metric_dispatch(::MetricManifold{ℂ,CircleGroup,EuclideanMetric}) = Val(true)

adjoint_action(::CircleGroup, p, X) = X

adjoint_action!(::CircleGroup, Y, p, X) = copyto!(Y, X)

function _compose(G::CircleGroup, p::AbstractVector, q::AbstractVector)
    return map(compose, repeated(G), p, q)
end

_compose!(G::CircleGroup, x, p, q) = copyto!(x, compose(G, p, q))

identity_element(G::CircleGroup) = 1.0
identity_element(::CircleGroup, p::Number) = one(p)
identity_element(::CircleGroup, p::AbstractArray) = map(i -> one(eltype(p)), p)

Base.inv(G::CircleGroup, p::AbstractVector) = map(inv, repeated(G), p)

function inverse_translate(
    ::CircleGroup,
    p::AbstractVector,
    q::AbstractVector,
    ::LeftAction,
)
    return map(/, q, p)
end
function inverse_translate(
    ::CircleGroup,
    p::AbstractVector,
    q::AbstractVector,
    ::RightAction,
)
    return map(/, q, p)
end

lie_bracket(::CircleGroup, X, Y) = zero(X)

lie_bracket!(::CircleGroup, Z, X, Y) = fill!(Z, 0)

translate_diff(::GT, p, q, X, ::ActionDirection) where {GT<:CircleGroup} = map(*, p, X)
function translate_diff(
    ::CircleGroup,
    ::Identity{MultiplicationOperation},
    q,
    X,
    ::ActionDirection,
)
    return X
end

function translate_diff!(G::CircleGroup, Y, p, q, X, conv::ActionDirection)
    return copyto!(Y, translate_diff(G, p, q, X, conv))
end

function exp_lie(::CircleGroup, X)
    return map(X) do imθ
        θ = imag(imθ)
        sinθ, cosθ = sincos(θ)
        return Complex(cosθ, sinθ)
    end
end

exp_lie!(G::CircleGroup, q, X) = (q .= exp_lie(G, X))

function log_lie(::CircleGroup, q)
    return map(q) do z
        cosθ, sinθ = reim(z)
        θ = atan(sinθ, cosθ)
        return θ * im
    end
end
log_lie(::CircleGroup, e::Identity{MultiplicationOperation}) = 0.0 * im

_log_lie!(G::CircleGroup, X, q) = (X .= log_lie(G, q))

function number_of_coordinates(G::CircleGroup, B::AbstractBasis)
    return number_of_coordinates(base_manifold(G), B)
end

@doc raw"""
    RealCircleGroup <: GroupManifold{Circle{ℝ},AdditionOperation}

The real circle group is the real circle ([`Circle(ℝ)`](@ref)) equipped with
the group operation of addition ([`AdditionOperation`](@ref)).
"""
const RealCircleGroup = GroupManifold{ℝ,Circle{ℝ},AdditionOperation}

RealCircleGroup() = GroupManifold(Circle{ℝ}(), AdditionOperation())

Base.show(io::IO, ::RealCircleGroup) = print(io, "RealCircleGroup()")

invariant_metric_dispatch(::RealCircleGroup, ::ActionDirection) = Val(true)

default_metric_dispatch(::MetricManifold{ℝ,RealCircleGroup,EuclideanMetric}) = Val(true)

_compose(::RealCircleGroup, p, q) = sym_rem(p + q)
function _compose(G::RealCircleGroup, p::AbstractVector, q::AbstractVector)
    return map(compose, repeated(G), p, q)
end

function _compose!(::RealCircleGroup, x, p, q)
    x .= sym_rem.(p .+ q)
    return x
end

identity_element(G::RealCircleGroup) = 0.0
identity_element(::RealCircleGroup, p::AbstractArray) = map(i -> zero(eltype(p)), p)

Base.inv(G::RealCircleGroup, p::AbstractVector) = map(inv, repeated(G), p)

function inverse_translate(
    ::RealCircleGroup,
    p::AbstractVector,
    q::AbstractVector,
    ::LeftAction,
)
    return map((x, y) -> sym_rem(x - y), q, p)
end
function inverse_translate(
    ::RealCircleGroup,
    p::AbstractVector,
    q::AbstractVector,
    ::RightAction,
)
    return map((x, y) -> sym_rem(x - y), q, p)
end

function exp_lie(::RealCircleGroup, X)
    return sym_rem(X)
end

exp_lie!(G::RealCircleGroup, q, X) = (q .= sym_rem(X))


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
using FoodCachingPlotting, Stipple, StippleUI, StipplePlotly, Genie, DataFrames,
      Random, PlotlyJS, FoodCachingExperiments, FoodCachingModels
import FoodCachingExperiments: EXPERIMENTS, CLAYTON0103_EXPERIMENTS
const M = FoodCachingModels

Base.@kwdef mutable struct Model <: ReactiveModel
    plot::R{Bool} = false
    birdid::R{Int} = 0
    birdid_options::R{Vector{Int}} = []
    imgurl::R{String} = ""
    filter::R{String} = ""
    seed::R{String} = "random"
    models::Vector{Any} = []
    list::R{String} = ""
    layout::R{Any} = nothing
    plot_data::R{Any} = nothing
    dist_layout::R{Any} = nothing
    dist_plots::R{Any} = nothing
    summary1_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary1_plots::R{Any} = nothing
    summary2_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary2_plots::R{Any} = nothing
    summary3_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary3_plots::R{Any} = nothing
    summary4_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary4_plots::R{Any} = nothing
    summary5_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary5_plots::R{Any} = nothing
    summary6_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary6_plots::R{Any} = nothing
    summary7_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary7_plots::R{Any} = nothing
    summary8_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary8_plots::R{Any} = nothing
    summary9_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary9_plots::R{Any} = nothing
    summary10_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary10_plots::R{Any} = nothing
    summary11_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary11_plots::R{Any} = nothing
    summary12_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary12_plots::R{Any} = nothing
    summary13_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary13_plots::R{Any} = nothing
    summary14_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary14_plots::R{Any} = nothing
    summary15_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary15_plots::R{Any} = nothing
    summary16_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary16_plots::R{Any} = nothing
    summary17_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary17_plots::R{Any} = nothing
    summary18_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary18_plots::R{Any} = nothing
    summary19_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary19_plots::R{Any} = nothing
    summary20_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary20_plots::R{Any} = nothing
    summary21_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary21_plots::R{Any} = nothing
    summary22_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary22_plots::R{Any} = nothing
    summary23_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary23_plots::R{Any} = nothing
    summary24_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary24_plots::R{Any} = nothing
    summary25_layout::R{Any} = PlotlyJS.Layout(height = 10, width = 10)
    summary25_plots::R{Any} = nothing
    experiment::R{String} = ""
    datadir::R{String} = joinpath(@__DIR__, "..", "data")
    files::R{Vector{String}} = []
    opt_experiment::Vector{String} = string.(setdiff(collect(keys(EXPERIMENTS)),
                                                     CLAYTON0103_EXPERIMENTS))
    model::R{String} = ""
    opt_models::Vector{String} = ["Baseline", "MotivationalControl",
                                  "EpisodicLikeMemory", "PlasticCaching",
                                  "ReplayAndPlan"]
    commit::R{String} = ""
    opt_commit::R{Vector{String}} = [""]
    id::R{String} = ""
    opt_id::R{Vector{String}} = [""]
end

function updatefiles(model, dir)
    if ispath(dir)
        model.files[] = readdir(dir)
    end
end

function runmodel(m, experiment, trackedfields, seed)
    p = M.parameters(M.Model{:euler};
                      m.p.fixed...,
                      tracker = () -> M.Tracker(trackedfields = trackedfields))
    m = M.Population(m.m, m.s, m.l, m.u, p, m.dist, M.init(p))
    Random.seed!(seed)
    ms = rand(m, nbirds(experiment))
    data = run!(experiment, ms)
    ms, data
end

function default_trackedfields(p)
    pd = Dict(p...)
    res = [(:cage, :cacheableitems),
     (:cage, :eatableitems),
     (:cage, :ismdpresent),
     (:cage, :trays),
     :actions]
    if pd[:agent] ∈ [M.SpecSatAgent, M.PlasticCachingAgent]
        push!(res, (:agent, :hungermodel, :hunger))
    end
#     if haskey(pd, :cacheparams) && pd[:cacheparams] <: M.PlanningAgentParams
#         push!(res, (:agent, :cacheparams, :current))
#         push!(res, (:agent, :cacheparams, :memory))
#     end
    if haskey(pd, :cacheparams) && pd[:cacheparams] <: M.PlasticCachingAgentParams
        push!(res, (:agent, :cacheparams, :trayweights))
    end
    res
end

function parse_seed(seed)
    seed == "best" && return 0
    seed == "random" && return UInt(time_ns())
    parse(UInt, seed)
end

stringify(x) = x == "" ? "" : "_$x"
function setfilename()
    f = "$(model.model[])_$(model.experiment[])$(stringify(model.id[]))$(stringify(model.commit[])).bson.zstd"
    if f ∈ model.files[]
        model.list[] = f
    end
end

function plot_bird(ms, id)
    println("plotting bird $id")
    tmp = FoodCachingPlotting.plotexperiment(ms[id])
    model.layout[] = tmp.plot.layout
    model.plot_data[] = tmp.plot.data
end

function titles(m)
    vcat([fill(string(k), M.NestedStructInitialiser.free_param_length(v))
          for (k, v) in m.p.free]...)
end

function plot_dists(m; x = 0:.01:1)
    n = length(m.m)
    n1 = ceil(Int, sqrt(n))
    n2 = ceil(Int, n / n1)
    xaxis = Dict([Symbol("xaxis", i) => attr(automargin = false,
                                             title = attr(text = n, standoff = 0))
                  for (i, n) in enumerate(titles(m))]...)
    model.dist_layout[] = PlotlyJS.Layout(; grid = attr(ygap = .6,
                                                        rows = n1,
                                                        columns = n2,
                                                        pattern = "independent"),
                                           showlegend = false,
                                           height = 900,
                                           xaxis...)
    model.dist_plots[] = [PlotData(x = (u - l) * x .+ l,
                                   y = M.Distributions.pdf(m.dist(d, s), x),
                                   yaxis = "y$i", xaxis = "x$i")
                          for (i, u, l, d, s) in zip(1:length(m.u), m.u, m.l, m.m, m.s)]
end

function plot_summary(e, data)
    summary_plots = plot_compare(e, FoodCachingExperiments.summarize(e, data))
    for i in 1:25
        if i <= length(summary_plots)
            getproperty(model, Symbol("summary$(i)_plots"))[] = summary_plots[i]
            getproperty(model, Symbol("summary$(i)_layout"))[] = PlotlyJS.Layout(height = 300, width = 600)
        else
            getproperty(model, Symbol("summary$(i)_layout"))[] = PlotlyJS.Layout(height = 10, width = 10)
        end
    end
end

model = Stipple.init(Model())

updatefiles(model, model.datadir[])

on(model.birdid) do id
    isempty(id) || plot_bird(model.models, id)
end

on(model.datadir) do dir
    updatefiles(model, dir)
end

on(model.plot) do plot
    println("plot changed $plot")
    plot || return
    setfilename()
    f = model.list[]
    println("Loading $f ...")
    m = load(joinpath(model.datadir[], f))
    plot_dists(m)
    println("loaded")
    e = Symbol(model.experiment[])
    println("running")
    ms, data = runmodel(m, e,
                        default_trackedfields(m.p.fixed),
                        parse_seed(model.seed[]))
    model.models = ms
    plot_summary(e, data)
    model.birdid_options[] = 1:length(ms)
    model.birdid[] = 1
    model.plot[] = false
end

function ui()
    page(vm(model), class = "container", [input("", @bind(:datadir),
                                                  style = "width: 500px"),
         row(class = "st-br",
             [cell(class = "st-module", [
                   h6("model"),
                   quasar(:select, "", @bind(:model), options = :opt_models)]),
              cell(class = "st-module", [
                   h6("experiment"),
                   quasar(:select, "", @bind(:experiment), options = :opt_experiment)]),
              cell(class = "st-module", [
                   h6("commit"),
                   quasar(:select, "", @bind(:commit), options = :opt_commit)]),
              cell(class = "st-module", [
                   h6("id"),
                   quasar(:select, "", @bind(:id), options = :opt_id)]),
              cell(class = "st-module", [
                   h6("bird"),
                   quasar(:select, "", @bind(:birdid), options = :birdid_options)]),
              cell(class = "st-module", [
                   h6("seed"), "<br>",
                      input("", @bind(:seed), style = "width: 75px")]),
                     row(class = "justify-end",
                              btn("plot", @click("plot = true"),
                                  style = "background: #FF0000; color: white")
                ),
               ]),
         row(cell(class = "st-module", [
            StipplePlotly.plot(:plot_data, layout = :layout)
           ])),
         row(cell(class = "st-module", [
            StipplePlotly.plot(:summary1_plots, layout = :summary1_layout),
            StipplePlotly.plot(:summary2_plots, layout = :summary2_layout)
           ])),
         row(cell(class = "st-module", [
            StipplePlotly.plot(:summary3_plots, layout = :summary3_layout),
            StipplePlotly.plot(:summary4_plots, layout = :summary4_layout)
           ])),
         row(cell(class = "st-module", [
            StipplePlotly.plot(:summary5_plots, layout = :summary5_layout),
            StipplePlotly.plot(:summary6_plots, layout = :summary6_layout)
           ])),
         row(cell(class = "st-module", [
            StipplePlotly.plot(:summary7_plots, layout = :summary7_layout),
            StipplePlotly.plot(:summary8_plots, layout = :summary8_layout)
           ])),
         row(cell(class = "st-module", [
            StipplePlotly.plot(:summary9_plots, layout = :summary9_layout),
            StipplePlotly.plot(:summary10_plots, layout = :summary10_layout)
           ])),
         row(cell(class = "st-module", [
            StipplePlotly.plot(:summary11_plots, layout = :summary11_layout),
            StipplePlotly.plot(:summary12_plots, layout = :summary12_layout)
           ])),
         row(cell(class = "st-module", [
            StipplePlotly.plot(:summary13_plots, layout = :summary13_layout),
            StipplePlotly.plot(:summary14_plots, layout = :summary14_layout)
           ])),
         row(cell(class = "st-module", [
            StipplePlotly.plot(:summary15_plots, layout = :summary15_layout),
            StipplePlotly.plot(:summary16_plots, layout = :summary16_layout)
           ])),
         row(cell(class = "st-module", [
            StipplePlotly.plot(:summary17_plots, layout = :summary17_layout),
            StipplePlotly.plot(:summary18_plots, layout = :summary18_layout)
           ])),
         row(cell(class = "st-module", [
            StipplePlotly.plot(:summary19_plots, layout = :summary19_layout),
            StipplePlotly.plot(:summary20_plots, layout = :summary20_layout)
           ])),
         row(cell(class = "st-module", [
            StipplePlotly.plot(:summary21_plots, layout = :summary21_layout),
            StipplePlotly.plot(:summary22_plots, layout = :summary22_layout)
           ])),
         row(cell(class = "st-module", [
            StipplePlotly.plot(:summary23_plots, layout = :summary23_layout),
            StipplePlotly.plot(:summary24_plots, layout = :summary24_layout)
           ])),
         row(cell(class = "st-module", [
            StipplePlotly.plot(:summary25_plots, layout = :summary25_layout),
           ])),
         row(cell(class = "st-module", [
            StipplePlotly.plot(:dist_plots, layout = :dist_layout)
           ])),
         ]) |> html
end

route("/", ui)
Genie.config.server_host = "127.0.0.1"
up()


text:
const tmp_dir = joinpath(expanduser("~"), ".JITrench")
const dot_file_path = joinpath(tmp_dir, "tmp_graph.dot")
colors = Dict(
    "func" => "lightblue",
    "var" => "orange",
    "user_defined_var" => "orange"
)

"""
    PNGContainer
A structure for handling images.
You can display an image by using it as a return value of a cell or by explicitly `display(::PNGContainer)` in Jupyter.

"""
struct PNGContainer
    content
end

function Base.show(io::IO, ::MIME"image/png", c::PNGContainer)
    write(io, c.content)
end

function _dot_var(var::Variable{T}) where T

    label = ""
    value = var.values
    if value !== nothing
        var_size = size(value)
        if isempty(var_size)
            try value !== nothing
                label *= "$(value)"
            catch
                label *= "nothing"
            end
        else
            label *= "shape: $(var_size) \n type: $(T))"
        end
    end
    if var.creator === nothing
        color = colors["user_defined_var"]
        shape = "polygon, sides=8"
    else
        color = colors["var"]
        shape = "ellipse"
    end
    if var.name === nothing
        dot_var = "$(objectid(var)) [shape=$shape, label=\"$label\", color=\"$color\", style=filled]\n"
    else
        name = var.name * "
"
        dot_var = "$(objectid(var)) [shape=$shape, label=<$name$label>, color=\"$color\", style=filled]\n"
    end
    return dot_var
end

function _dot_func(f::DiffableFunction)
    f_type = typeof(f)
    txt = "$(objectid(f)) [label=\"$(f_type)\", color=\"$(colors["func"])\", style=filled, shape=box]\n"
    for x in f.grad_field.inputs
        txt *= "$(objectid(x)) -> $(objectid(f))\n"
    end
    for y in f.grad_field.outputs
        txt *= "$(objectid(f)) -> $(objectid(y))\n"
    end
    return txt
end

function get_dot_graph(var, title)
    txt = ""
    funcs = []
    seen_set = Set{DiffableFunction}()
    push!(funcs, var.creator)
    txt = _dot_var(var)
    while !(isempty(funcs))
        f = pop!(funcs)
        txt *= _dot_func(f)
        for x in f.grad_field.inputs
            txt *= _dot_var(x)
            if x.creator !== nothing && (!(x.creator in seen_set))
                push!(seen_set, x.creator)
                push!(funcs, x.creator)
            end
        end
    end
    return "digraph g {
            graph [
                labelloc=\"t\";
                label= \"$(title)\"
            ];
                $txt
            }"
end

function plot_tmp_dir()
    extension = "png"
    to_file = joinpath(tmp_dir, "graph." * extension)
    cmd = `dot $(dot_file_path) -T $(extension) -o $(to_file)`
    run(cmd)
    return to_file
end

"""
    plot_graph(var::Variable; to_file="", title="")
Draw a computational graph with y as the end. This requires the installation of graphviz.

# Arguments
- `to_file`: Specifies the file name where the image will be saved.
If it is an empty string, it will return an image of type `PNGContainer`. See also: [`PNGContainer`](@ref)
- `title`: Title of graph.

```julia-repl
julia> x = Variable(10)
name: nothing
values: 10
creator: User-Defined(nothing)

julia> y = x + 3
name: nothing
values: 13
creator: JITrench.Add

julia> JITrench.plot_graph(y, to_file="graph.png")
```
"""
function plot_graph(var::Variable; to_file="", title="")
    dot_graph = get_dot_graph(var, title)
    (!(ispath(tmp_dir))) && (mkdir(tmp_dir))
    open(dot_file_path, "w") do io
        write(io, dot_graph)
    end

    if to_file == ""
        png_file_path = plot_tmp_dir()
            c = open(png_file_path) do io
            PNGContainer(read(io))
        end
        return c
    else
        extension = split(to_file, ".")[end]
        cmd = `dot $(dot_file_path) -T $(extension) -o $(to_file)`
        run(cmd)
    end
end


text:
@testset "Write/Read Data and Checkpoints" begin
    @testset "Checkpoints D1Q3" begin
        typeset = [Float16 Float32 Float64]
        for T in typeset
            testdist = JuSwalbe.DistributionD1Q3(f0=ones(T, 10), f1=fill(T(0.1),10), f2=fill(T(-0.1),10))
            savecheckpoints(testdist)

            @test isfile("data/tmp/checkpoint_0.bson")
            readdist = loadcheckpoint("data/tmp/checkpoint_0.bson")
            @test testdist.f0 == readdist.f0
            @test testdist.f1 == readdist.f1
            @test testdist.f2 == readdist.f2

            rm("data", recursive=true)
        end
    end
    @testset "Checkpoints D2Q9" begin
        typeset = [Float16 Float32 Float64]
        for T in typeset
            testdist = simpledistD2Q9(5,5, T=T)
            @test isa(testdist, JuSwalbe.DistributionD2Q9{Matrix{T}})
            savecheckpoints(testdist)

            @test isfile("data/tmp/checkpoint_0.bson")
            readdist = loadcheckpoint("data/tmp/checkpoint_0.bson")
            @test testdist.f0 == readdist.f0
            @test testdist.f1 == readdist.f1
            @test testdist.f2 == readdist.f2
            @test testdist.f3 == readdist.f3
            @test testdist.f4 == readdist.f4
            @test testdist.f5 == readdist.f5
            @test testdist.f6 == readdist.f6
            @test testdist.f7 == readdist.f7
            @test testdist.f8 == readdist.f8

            rm("data", recursive=true)
        end
    end
    @testset "Checkpoint D1Q3" begin
        typeset = [Float16 Float32 Float64]
        for T in typeset
            testdist = JuSwalbe.DistributionD1Q3(f0=ones(T, 10), f1=fill(T(0.1),10), f2=fill(T(-0.1),10))
            savecheckpoint(testdist)

            @test isfile("data/tmp/checkpoint_0.bson")
            readdist = loadcheckpoint("data/tmp/checkpoint_0.bson")
            @test testdist.f0 == readdist.f0
            @test testdist.f1 == readdist.f1
            @test testdist.f2 == readdist.f2

            rm("data", recursive=true)
        end
    end
    @testset "Checkpoints D2Q9" begin
        typeset = [Float16 Float32 Float64]
        for T in typeset
            testdist = simpledistD2Q9(5,5, T=T)
            @test isa(testdist, JuSwalbe.DistributionD2Q9{Matrix{T}})
            savecheckpoint(testdist)

            @test isfile("data/tmp/checkpoint_0.bson")
            readdist = loadcheckpoint("data/tmp/checkpoint_0.bson")
            @test testdist.f0 == readdist.f0
            @test testdist.f1 == readdist.f1
            @test testdist.f2 == readdist.f2
            @test testdist.f3 == readdist.f3
            @test testdist.f4 == readdist.f4
            @test testdist.f5 == readdist.f5
            @test testdist.f6 == readdist.f6
            @test testdist.f7 == readdist.f7
            @test testdist.f8 == readdist.f8

            rm("data", recursive=true)
        end
    end
    @testset "Save height data one dimension" begin
        for type in [Float16 Float32 Float64]
            mom = JuSwalbe.Macroquant(height = ones(type, 10), velocity = zeros(type, 10), pressure = zeros(type, 10), energy = zeros(type, 10))
            @test isa(mom, JuSwalbe.Macroquant{Vector{type}, Vector{type}})
            height2file(mom)
            @test isfile("data/tmp/height_0.bson")
            h = BSON.load("data/tmp/height_0.bson")
            @test length(h[:height]) == 10
            @test mom.height == h[:height]
        end
        rm("data", recursive=true)
    end
    @testset "Save height data two dimensions" begin
        for type in [Float16 Float32 Float64]
            mom = JuSwalbe.Macroquant(height = ones(type, (5,5)), velocity = simpleTwovector(5,5, T=type), pressure = zeros(type, (5,5)), energy = zeros(type, (5,5)))
            @test isa(mom, JuSwalbe.Macroquant{Matrix{type}, JuSwalbe.Twovector{Matrix{type}}})
            height2file(mom)
            @test isfile("data/tmp/height_0.bson")
            h = BSON.load("data/tmp/height_0.bson")
            @test size(h[:height]) == (5,5)
            @test mom.height == h[:height]
        end
        rm("data", recursive=true)
    end
    @testset "Save velocity data one dimension" begin
        for type in [Float16 Float32 Float64]
            mom = JuSwalbe.Macroquant(height = ones(type, 10), velocity = fill(type(0.1), 10), pressure = zeros(type, 10), energy = zeros(type, 10))
            @test isa(mom, JuSwalbe.Macroquant{Vector{type}, Vector{type}})
            velocity2file(mom)
            @test isfile("data/tmp/velocity_0.bson")
            v = BSON.load("data/tmp/velocity_0.bson")
            @test length(v[:velocity]) == 10
            @test mom.velocity == v[:velocity]
        end
        rm("data", recursive=true)
    end
    @testset "Save velocity data two dimensions" begin
        for type in [Float16 Float32 Float64]
            mom = JuSwalbe.Macroquant(height = ones(type, (5,5)), velocity = simpleTwovector(5,5, T=type), pressure = zeros(type, (5,5)), energy = zeros(type, (5,5)))
            @test isa(mom, JuSwalbe.Macroquant{Matrix{type}, JuSwalbe.Twovector{Matrix{type}}})
            velocity2file(mom)
            @test isfile("data/tmp/velocity_0.bson")
            v = BSON.load("data/tmp/velocity_0.bson")
            @test size(v[:velocity_x]) == (5,5)
            @test size(v[:velocity_y]) == (5,5)
            @test mom.velocity.x == v[:velocity_x]
            @test mom.velocity.y == v[:velocity_y]
        end
        rm("data", recursive=true)
    end
    @testset "Save velocity and height data one dimension" begin
        for type in [Float16 Float32 Float64]
            mom = JuSwalbe.Macroquant(height = ones(type, 10), velocity = fill(type(0.1), 10), pressure = zeros(type, 10), energy = zeros(type, 10))
            @test isa(mom, JuSwalbe.Macroquant{Vector{type}, Vector{type}})
            velocityandheight2file(mom)
            @test isfile("data/tmp/height_velocity_0.bson")
            hv = BSON.load("data/tmp/height_velocity_0.bson")
            @test length(hv[:velocity]) == 10
            @test length(hv[:height]) == 10
            @test mom.velocity == hv[:velocity]
            @test mom.height == hv[:height]
        end
        rm("data", recursive=true)
    end
    @testset "Save velocity and height data two dimensions" begin
        for type in [Float16 Float32 Float64]
            mom = JuSwalbe.Macroquant(height = ones(type, (5,5)), velocity = simpleTwovector(5,5, T=type), pressure = zeros(type, (5,5)), energy = zeros(type, (5,5)))
            @test isa(mom, JuSwalbe.Macroquant{Matrix{type}, JuSwalbe.Twovector{Matrix{type}}})
            velocityandheight2file(mom)
            @test isfile("data/tmp/height_velocity_0.bson")
            hv = BSON.load("data/tmp/height_velocity_0.bson")
            @test size(hv[:height]) == (5,5)
            @test size(hv[:velocity_x]) == (5,5)
            @test size(hv[:velocity_y]) == (5,5)
            @test mom.height == hv[:height]
            @test mom.velocity.x == hv[:velocity_x]
            @test mom.velocity.y == hv[:velocity_y]
        end
        rm("data", recursive=true)
    end
end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
using Documenter, FactorGraph

makedocs(
    sitename = "FactorGraph",
    modules = [FactorGraph],
    clean = false,
    doctest = false,
    format = Documenter.HTML(
        assets = ["assets/style.css"],
        prettyurls = prettyurls = get(ENV, "CI", nothing) == "true",
        collapselevel = 2,
    ),
    pages = [
        "Home" => "index.md",
        "Continuous Framework" => [
        "Input Data" => "man/continuousInput.md",
        "General Factor Graph" => [
            "Graphical Model" => "man/continuousModel.md",
            "Inference" => "man/continuousInference.md",
            "Output Data" => "man/continuousOutput.md",
        ],
        "Tree Factor Graph" => [
            "Graphical Model" => "man/continuousTreeModel.md",
            "Inference" => "man/continuousTreeInference.md",
            "Output Data" => "man/continuousTreeOutput.md",
        ],
        "Utility Functions" => "man/utility.md"
        ],
        "Discrete Framework" => [
        "Input Data" => "man/discreteInput.md",
        "Tree Factor Graph" => [
            "Graphical Model" => "man/discreteTreeModel.md",
            "Inference" => "man/discreteTreeInference.md",
            "Output Data" => "man/discreteTreeOutput.md",
        ],
        ],
        "Theoretical Background" => [
            "Inference in Factor Graphs" => "man/theoreticalInference.md",
            "Continuous Variables" => "man/theoreticalBelief.md"
        ],
    ],
)

deploydocs(
    repo = "github.com/mcosovic/FactorGraph.jl.git",
    target = "build",
)


text:
using Documenter, Iex

CI_FLG = get(ENV, "CI", nothing) == "true"

makedocs(
    modules = [Iex],
    format = Documenter.HTML(
        prettyurls = CI_FLG
    ),
    sitename="IexJulia")

deploydocs(
    repo = "github.com/GodolphinCapitalManagement/Iex.git",
    push_preview = true,
)

max_stars_repo_name: GodolphinCapitalManagement/Iex

max_issues_repo_name: GodolphinCapitalManagement/Iex

max_forks_repo_name: GodolphinCapitalManagement/Iex


text:
FactSheet

ISIN = XS1391063424
Name = BURFORD CAPITAL PLC 6.125% GTD BDS 26/10/24
IssueDate = 20-Apr-2016
MaturityDate = 26-Oct-2024
RedemptionType = single / periodic / call option / put option
Redemption = 26-Oct-2024

CouponRate = +6.125%
CouponType = Fixed
IndexInflationDate = null
ReferenceRate = null
CouponDates = missing
CouponFrequency = 6 months
PreviousCouponPaymentDate = 19-Oct-2018
CouponExDate = null

DenominationCurrency = GBP
UnitOfQuotation = 100

Calendar = missing
SettlementCycle = T+3
ConventionCalculate = ACT/ACT ISMA
ConventionBusinessDays = missing