text:
using GeometryBasics
using Statistics
using CategoricalArrays: categorical, levelcode

@cell "Test heatmap + image overlap" begin
    heatmap(RNG.rand(32, 32))
    image!(map(x -> RGBAf0(x, 0.5, 0.5, 0.8), RNG.rand(32, 32)))
    current_figure()
end

@cell "poly and colormap" begin
    # example by @Paulms from JuliaPlots/Makie.jl#310
    points = Point2f0[[0.0, 0.0], [0.1, 0.0], [0.1, 0.1], [0.0, 0.1]]
    colors = [0.0 ,0.0, 0.5, 0.0]
    fig, ax, polyplot = poly(points, color=colors, colorrange=(0.0, 1.0))
    points = Point2f0[[0.1, 0.1], [0.2, 0.1], [0.2, 0.2], [0.1, 0.2]]
    colors = [0.5,0.5,1.0,0.3]
    poly!(ax, points, color=colors, colorrange=(0.0, 1.0))
    fig
end

@cell "quiver" begin
    x = range(-2, stop=2, length=21)
    arrows(x, x, RNG.rand(21, 21), RNG.rand(21, 21), arrowsize=0.05)
end

@cell "Arrows on hemisphere" begin
    s = Sphere(Point3f0(0), 0.9f0)
    fig, ax, meshplot = mesh(s, transparency=true, alpha=0.05)
    pos = decompose(Point3f0, s)
    dirs = decompose_normals(s)
    arrows!(ax, pos, dirs, arrowcolor=:red, arrowsize=0.1, linecolor=:red)
    fig
end

@cell "image" begin
    fig = Figure()
    image(fig[1,1], Makie.logo(), axis = (; aspect = DataAspect()))
    image(fig[1, 2], RNG.rand(100, 500), axis = (; aspect = DataAspect()))
    fig
end

@cell "FEM polygon 2D" begin
    coordinates = [
        0.0 0.0;
        0.5 0.0;
        1.0 0.0;
        0.0 0.5;
        0.5 0.5;
        1.0 0.5;
        0.0 1.0;
        0.5 1.0;
        1.0 1.0;
    ]
    connectivity = [
        1 2 5;
        1 4 5;
        2 3 6;
        2 5 6;
        4 5 8;
        4 7 8;
        5 6 9;
        5 8 9;
    ]
    color = [0.0, 0.0, 0.0, 0.0, -0.375, 0.0, 0.0, 0.0, 0.0]
    poly(coordinates, connectivity, color=color, strokecolor=(:black, 0.6), strokewidth=4)
end

@cell "FEM mesh 2D" begin
    coordinates = [
        0.0 0.0;
        0.5 0.0;
        1.0 0.0;
        0.0 0.5;
        0.5 0.5;
        1.0 0.5;
        0.0 1.0;
        0.5 1.0;
        1.0 1.0;
    ]
    connectivity = [
        1 2 5;
        1 4 5;
        2 3 6;
        2 5 6;
        4 5 8;
        4 7 8;
        5 6 9;
        5 8 9;
    ]
    color = [0.0, 0.0, 0.0, 0.0, -0.375, 0.0, 0.0, 0.0, 0.0]
    fig, ax, meshplot = mesh(coordinates, connectivity, color=color, shading=false)
    wireframe!(ax, meshplot[1], color=(:black, 0.6), linewidth=3)
    fig
end

@cell "colored triangle" begin
    mesh(
        [(0.0, 0.0), (0.5, 1.0), (1.0, 0.0)], color=[:red, :green, :blue],
        shading=false
    )
end

@cell "colored triangle with poly" begin
    poly(
        [(0.0, 0.0), (0.5, 1.0), (1.0, 0.0)],
        color=[:red, :green, :blue],
        strokecolor=:black, strokewidth=2
    )
end

# @cell "Subscenes" begin
#     img = RNG.rand(RGBAf0, 100, 100)
#     scene = image(img, show_axis=false)
#     subscene = Scene(scene, IRect(100, 100, 300, 300))
#     scatter!(subscene, RNG.rand(100) * 200, RNG.rand(100) * 200, markersize=4)
#     scene
# end

@cell "scale_plot" begin
    t = range(0, stop=1, length=500) # time steps
    θ = (6π) .* t    # angles
    x =  # x coords of spiral
    y =  # y coords of spiral
    lines(t .* cos.(θ), t .* sin.(θ);
        color=t, colormap=:algae, linewidth=8, axis = (; aspect = DataAspect()))
end

@cell "Polygons" begin
    points = decompose(Point2f0, Circle(Point2f0(50), 50f0))
    fig, ax, pol = poly(points, color=:gray, strokewidth=10, strokecolor=:red)
    # Optimized forms
    poly!(ax, [Circle(Point2f0(50 + 300), 50f0)], color=:gray, strokewidth=10, strokecolor=:red)
    poly!(ax, [Circle(Point2f0(50 + i, 50 + i), 10f0) for i = 1:100:400], color=:red)
    poly!(ax, [FRect2D(50 + i, 50 + i, 20, 20) for i = 1:100:400], strokewidth=2, strokecolor=:green)
    linesegments!(ax,
        [Point2f0(50 + i, 50 + i) => Point2f0(i + 70, i + 70) for i = 1:100:400], linewidth=8, color=:purple
    )
    fig
end

@cell "Text Annotation" begin
    text(
        ". This is an annotation!",
        position=(300, 200),
        align=(:center,  :center),
        textsize=60,
        font="Blackchancery"
    )
end

@cell "Text rotation" begin
    fig = Figure()
    ax = fig[1, 1] = Axis(fig)
    pos = (500, 500)
    posis = Point2f0[]
    for r in range(0, stop=2pi, length=20)
        p = pos .+ (sin(r) * 100.0, cos(r) * 100)
        push!(posis, p)
        text!(ax, "test",
            position=p,
            textsize=50,
            rotation=1.5pi - r,
            align=(:center, :center)
        )
    end
    scatter!(ax, posis, markersize=10)
    fig
end

@cell "Standard deviation band" begin
    # Sample 100 Brownian motion path and plot the mean trajectory together
    # with a ±1σ band (visualizing uncertainty as marginal standard deviation).
    n, m = 100, 101
    t = range(0, 1, length=m)
    X = cumsum(RNG.randn(n, m), dims=2)
    X = X .- X[:, 1]
    μ = vec(mean(X, dims=1)) # mean
    lines(t, μ)              # plot mean line
    σ = vec(std(X, dims=1))  # stddev
    band!(t, μ + σ, μ - σ)   # plot stddev band
    current_figure()
end

@cell "Streamplot animation" begin
    v(x::Point2{T}, t) where T = Point2{T}(one(T) * x[2] * t, 4 * x[1])
    sf = Node(Base.Fix2(v, 0e0))
    title_str = Node("t = 0.00")
    sp = streamplot(sf, -2..2, -2..2;
                    linewidth=2, padding=(0, 0),
                    arrow_size=0.09, colormap=:magma, axis=(;title=title_str))
    Record(sp, LinRange(0, 20, 5)) do i
        sf[] = Base.Fix2(v, i)
        title_str[] = "t = $(round(i; sigdigits=2))"
    end
end

@cell "Line changing colour" begin
    fig, ax, lineplot = lines(RNG.rand(10); linewidth=10)
    N = 20
    Record(fig, 1:N; framerate=20) do i
        lineplot.color = RGBf0(i / N, (N - i) / N, 0) # animate scene
    end
end

let
    struct FitzhughNagumo{T}
        ϵ::T
        s::T
        γ::T
        β::T
    end
    @cell "streamplot" begin
        P = FitzhughNagumo(0.1, 0.0, 1.5, 0.8)
        f(x, P::FitzhughNagumo) = Point2f0(
            (x[1] - x[2] - x[1]^3 + P.s) / P.ϵ,
            P.γ * x[1] - x[2] + P.β
        )
        f(x) = f(x, P)
        streamplot(f, -1.5..1.5, -1.5..1.5, colormap=:magma)
    end
end

@cell "Transforming lines" begin
    N = 7 # number of colours in default palette
    fig = Figure()
    ax = Axis(fig)
    fig[1,1] = ax
    st = Stepper(fig)

    xs = 0:9        # data
    ys = zeros(10)
    colors = Makie.default_palettes.color[]
    plots = map(1:N) do i # plot lines
        lines!(ax,
            xs, ys;
            color=colors[i],
            limits=FRect((0, 0), (10, 10)),
            linewidth=5
        ) # plot lines with colors
    end
    fig
    Makie.step!(st)

    for (i, rot) in enumerate(LinRange(0, π / 2, N))
        Makie.rotate!(plots[i], rot)
        arc!(ax,
            Point2f0(0),
            (8 - i),
            pi / 2,
            (pi / 2 - rot);
            color=plots[i].color,
            linewidth=5,
            linestyle=:dash
        )
    end

    Makie.step!(st)
    fig
end

@cell "Errorbars x y low high" begin
    x = 1:10
    y = sin.(x)
    fig, ax, scatterplot = scatter(x, y)
    errorbars!(ax, x, y, RNG.rand(10) .+ 0.5, RNG.rand(10) .+ 0.5)
    errorbars!(ax, x, y, RNG.rand(10) .+ 0.5, RNG.rand(10) .+ 0.5, color = :red, direction = :x)
    fig
end

@cell "Rangebars x y low high" begin
    vals = -1:0.1:1

    lows = zeros(length(vals))
    highs = LinRange(0.1, 0.4, length(vals))

    fig, ax, rbars = rangebars(vals, lows, highs, color = :red)
    rangebars!(ax, vals, lows, highs, color = LinRange(0, 1, length(vals)),
        whiskerwidth = 3, direction = :x)
    fig
end

@cell "Simple pie chart" begin
    fig = Figure(resolution=(800, 800))
    pie(fig[1, 1], 1:5, color=collect(1:5), axis=(;aspect=DataAspect()))
    fig
end

@cell "Hollow pie chart" begin
    pie(1:5, color=collect(1.0:5), radius=2, inner_radius=1, axis=(;aspect=DataAspect()))
end

@cell "Open pie chart" begin
    pie(0.1:0.1:1.0, normalize=false, axis=(;aspect=DataAspect()))
end

@cell "intersecting polygon" begin
    x = LinRange(0, 2pi, 100)
    poly(Point2f0.(zip(sin.(x), sin.(2x))), color = :white, strokecolor = :blue, strokewidth = 10)
end

@cell "Line Function" begin
    x = range(0, stop=3pi)
    fig, ax, lineplot = lines(x, sin.(x))
    lines!(ax, x, cos.(x), color=:blue)
    fig
end

@cell "Grouped bar" begin
	x1         = ["a_right", "a_right", "a_right", "a_right"]
	y1         = [2, 3, -3, -2]
	grp_dodge1 = [2, 2,  1,  1]
	grp_stack1 = [1, 2,  1,  2]

	x2         = ["z_left", "z_left", "z_left", "z_left"]
	y2         = [2, 3, -3, -2]
	grp_dodge2 = [1, 2,  1,  2]
	grp_stack2 = [1, 1,  2,  2]

	perm = [1, 4, 2, 7, 5, 3, 8, 6]
	x = [x1; x2][perm]
	x = categorical(x, levels = ["z_left", "a_right"])
	y = [y1; y2][perm]
	grp_dodge = [grp_dodge1; grp_dodge2][perm]
	grp_stack = [grp_stack1; grp_stack2][perm]

	tbl = (; x = x, grp_dodge = grp_dodge, grp_stack = grp_stack, y = y)

	fig = Figure()
	ax = Axis(fig[1,1])

	barplot!(ax, levelcode.(tbl.x), tbl.y, dodge = tbl.grp_dodge, stack = tbl.grp_stack, color = tbl.grp_stack)

	ax.xticks = (1:2, ["z_left", "a_right"])

	fig
end


text:
"""
    ArrayStyle(::Type{A})

Used to customize the meaning of indexing arguments in the context of a given array `A`.

See also: [`argdims`](@ref), [`UnsafeIndex`](@ref)
"""
abstract type ArrayStyle end

struct DefaultArrayStyle <: ArrayStyle end

ArrayStyle(A) = ArrayStyle(typeof(A))
ArrayStyle(::Type{A}) where {A} = DefaultArrayStyle()

"""
    argdims(::ArrayStyle, ::Type{T})

What is the dimensionality of the indexing argument of type `T`?
"""
argdims(x, arg) = argdims(x, typeof(arg))
argdims(x, ::Type{T}) where {T} = argdims(ArrayStyle(x), T)
argdims(s::ArrayStyle, arg) = argdims(s, typeof(arg))
# single elements initially map to 1 dimension but that dimension is subsequently dropped.
argdims(::ArrayStyle, ::Type{T}) where {T} = 0
argdims(::ArrayStyle, ::Type{T}) where {T<:Colon} = 1
argdims(::ArrayStyle, ::Type{T}) where {T<:AbstractArray} = ndims(T)
argdims(::ArrayStyle, ::Type{T}) where {N,T<:CartesianIndex{N}} = N
argdims(::ArrayStyle, ::Type{T}) where {N,T<:AbstractArray{CartesianIndex{N}}} = N
argdims(::ArrayStyle, ::Type{T}) where {N,T<:AbstractArray{<:Any,N}} = N
argdims(::ArrayStyle, ::Type{T}) where {N,T<:LogicalIndex{<:Any,<:AbstractArray{Bool,N}}} = N
@generated function argdims(s::ArrayStyle, ::Type{T}) where {N,T<:Tuple{Vararg{<:Any,N}}}
    e = Expr(:tuple)
    for p in T.parameters
        push!(e.args, :(ArrayInterface.argdims(s, $p)))
    end
    Expr(:block, Expr(:meta, :inline), e)
end

"""
    UnsafeIndex(::ArrayStyle, ::Type{I})

`UnsafeIndex` controls how indices that have been bounds checked and converted to
native axes' indices are used to return the stored values of an array. For example,
if the indices at each dimension are single integers then `UnsafeIndex(array, inds)` returns
`UnsafeGetElement()`. Conversely, if any of the indices are vectors then `UnsafeGetCollection()`
is returned, indicating that a new array needs to be reconstructed. This method permits
customizing the terminal behavior of the indexing pipeline based on arguments passed
to `ArrayInterface.getindex`. New subtypes of `UnsafeIndex` should define `promote_rule`.
"""
abstract type UnsafeIndex end

struct UnsafeGetElement <: UnsafeIndex end

struct UnsafeGetCollection <: UnsafeIndex end

UnsafeIndex(x, i) = UnsafeIndex(x, typeof(i))
UnsafeIndex(x, ::Type{I}) where {I} = UnsafeIndex(ArrayStyle(x), I)
UnsafeIndex(s::ArrayStyle, i) = UnsafeIndex(s, typeof(i))
UnsafeIndex(::ArrayStyle, ::Type{I}) where {I} = UnsafeGetElement()
UnsafeIndex(::ArrayStyle, ::Type{I}) where {I<:AbstractArray} = UnsafeGetCollection()

Base.promote_rule(::Type{X}, ::Type{Y}) where {X<:UnsafeIndex,Y<:UnsafeGetElement} = X

@generated function UnsafeIndex(s::ArrayStyle, ::Type{T}) where {N,T<:Tuple{Vararg{<:Any,N}}}
    if N === 0
        return UnsafeGetElement()
    else
        e = Expr(:call, promote_type)
        for p in T.parameters
            push!(e.args, :(typeof(ArrayInterface.UnsafeIndex(s, $p))))
        end
        return Expr(:block, Expr(:meta, :inline), Expr(:call, e))
    end
end

# are the indexing arguments provided a linear collection into a multidim collection
is_linear_indexing(A, args::Tuple{Arg}) where {Arg} = argdims(A, Arg) < 2
is_linear_indexing(A, args::Tuple{Arg,Vararg{Any}}) where {Arg} = false

"""
    flatten_args(A, args::Tuple{Arg,Vararg{Any}}) -> Tuple

This method may be used to flatten out multidimensional arguments across several
dimensions prior to performing indexing if any of `args` can be flattened.

See also: [`can_flatten](@ref)

# Extended help

If one wishes to create a new multidimensional argument that is altered prior to most of
the indexing pipeline, then it must be supported via the `can_flatten` and a new instance
of `flatten_args`, such as the following:

```julia

function ArrayInterface.flatten_args(A, args::Tuple{Arg,Vararg{Any}}) where {Arg<:NewIndexer}
    return (some_flattening_method(first(args)), flatten_args(s, tail(args))...)
end
```

Note that `A` is _NOT_ specified here. If different methods are necessary to flatten
types of `NewIndexer` and these are known prior to implementation, then they should be
specified using methods called within `flatten_args`:

```julia

function ArrayInterface.flatten_args(A, args::Tuple{Arg,Vararg{Any}}) where {Arg<:NewIndexer}
    return flatten_new_indexer(A, args)
end
flatten_new_indexer(A::Array1, args) = ...
flatten_new_indexer(A::Array2, args) = ...

```
"""
@inline function flatten_args(A, args::Tuple{Arg,Vararg{Any}}) where {Arg}
    return (first(args), flatten_args(A, tail(args))...)
end
@inline function flatten_args(A, args::Tuple{Arg,Vararg{Any}}) where {N,Arg<:CartesianIndex{N}}
    return (first(args).I..., flatten_args(A, tail(args))...)
end
@inline function flatten_args(A, args::Tuple{Arg,Vararg{Any}}) where {N,Arg<:LinearIndices{N}}
    return (eachindex(first(args)), flatten_args(A, tail(args))...)
end
@inline function flatten_args(A, args::Tuple{Arg,Vararg{Any}}) where {N,Arg<:CartesianIndices{N}}
    return (first(args).indices..., flatten_args(A, tail(args))...)
end
# we preserve CartesianIndices{0} for dropping dimensions
@inline function flatten_args(A, args::Tuple{Arg,Vararg{Any}}) where {Arg<:CartesianIndices{0}}
    return (first(args), flatten_args(A, tail(args))...)
end
@inline function flatten_args(A, args::Tuple{Arg,Vararg{Any}}) where {N,Arg<:AbstractArray{Bool,N}}
    if length(args) === 1
        if IndexStyle(A) isa IndexLinear
            return (LogicalIndex{Int}(first(args)),)
        else
            return (LogicalIndex(first(args)),)
        end
    else
        return (LogicalIndex(first(args)), flatten_args(A, tail(args))...)
    end
end
flatten_args(A, args::Tuple{}) = ()

"""
    can_flatten(::Type{A}, ::Type{T}) -> Bool

Returns `true` if an argument passed during indexing can be flattened across multiple
dimensions. For example, `CartesianIndex{N}` can be flattened as a series of `Int`s
across `N` dimensions. This method is used to trigger `flatten_args` prior to indexing.
If a particular argument-array combination cannot cannot be flattened, then it should be
specified here. Otherwise, `A` should not be specified when supporting a new
multidimensional indexing type. For example, the following is the typical usage:

```julia
ArrayInterface.can_flatten(::Type{A}, ::Type{T}) where {A,T<:NewIndexer} = true
```

but, in rare instances, this may be necessary:

```julia
ArrayInterface.can_flatten(::Type{A}, ::Type{T}) where {A<:ForbiddenArray,T<:NewIndexer} = false
```

"""
can_flatten(A, x) = can_flatten(typeof(A), typeof(x))
can_flatten(::Type{A}, ::Type{T}) where {A,T} = false
can_flatten(::Type{A}, ::Type{T}) where {A,I<:CartesianIndex,T<:AbstractArray{I}} = false
can_flatten(::Type{A}, ::Type{T}) where {A,T<: CartesianIndices} = true
can_flatten(::Type{A}, ::Type{T}) where {A,N,T<:AbstractArray{Bool,N}} = N > 1
can_flatten(::Type{A}, ::Type{T}) where {A,N,T<:CartesianIndex{N}} = true
@generated function can_flatten(::Type{A}, ::Type{T}) where {A,T<:Tuple}
    for i in T.parameters
        can_flatten(A, i) && return true
    end
    return false
end

"""
    to_indices(A, args::Tuple) -> to_indices(A, axes(A), args)
    to_indices(A, axes::Tuple, args::Tuple)

Maps arguments `args` to the axes of `A`. This is done by iteratively passing each
axis and argument to [`to_index`](@ref). Unique behavior based on the type of `A` may be
accomplished by overloading `to_indices(A, args)`. Unique axis-argument behavior can
be accomplished using `to_index(axis, arg)`.
"""
@propagate_inbounds function to_indices(A, args::Tuple)
    if can_flatten(A, args)
        return to_indices(A, flatten_args(A, args))
    elseif is_linear_indexing(A, args)
        return (to_index(eachindex(IndexLinear(), A), first(args)),)
    else
        return to_indices(A, axes(A), args)
    end
end
@propagate_inbounds to_indices(A, args::Tuple{}) = to_indices(A, axes(A), ())
@propagate_inbounds function to_indices(A, axs::Tuple, args::Tuple{Arg,Vararg{Any}}) where {Arg}
    N = argdims(A, Arg)
    if N > 1
        axes_front, axes_tail = Base.IteratorsMD.split(axs, Val(N))
        return (to_multi_index(axes_front, first(args)), to_indices(A, axes_tail, tail(args))...)
    else
        return (to_index(first(axs), first(args)), to_indices(A, tail(axs), tail(args))...)
    end
end
@propagate_inbounds function to_indices(A, axs::Tuple, args::Tuple{})
    @boundscheck if length(first(axs)) != 1
        error("Cannot drop dimension of size $(length(first(axs))).")
    end
    return to_indices(A, tail(axs), args)
end
@propagate_inbounds function to_indices(A, ::Tuple{}, args::Tuple{Arg,Vararg{Any}}) where {Arg}
    return (to_index(OneTo(1), first(args)), to_indices(A, (), tail(args))...)
end
to_indices(A, axs::Tuple{}, args::Tuple{}) = ()

_multi_check_index(axs::Tuple, arg) = _multi_check_index(axs, axes(arg))
_multi_check_index(axs::Tuple, arg::LogicalIndex) = axs == axes(arg.mask)
function _multi_check_index(axs::Tuple, arg::AbstractArray{T}) where {T<:CartesianIndex}
    b = true
    for i in arg
        b &= Base.checkbounds_indices(Bool, axs, (i,))
    end
    return b
end
_multi_check_index(::Tuple{}, ::Tuple{}) = true
function _multi_check_index(axs::Tuple, args::Tuple)
    if checkindex(Bool, first(axs), first(args))
        return _multi_check_index(tail(axs), tail(args))
    else
        return false
    end
end
@propagate_inbounds function to_multi_index(axs::Tuple, arg)
    @boundscheck if !_multi_check_index(axs, arg)
        throw(BoundsError(axs, arg))
    end
    return arg
end

"""
    to_index([::IndexStyle, ]axis, arg) -> index

Convert the argument `arg` that was originally passed to `getindex` for the dimension
corresponding to `axis` into a form for native indexing (`Int`, Vector{Int}, etc.). New
axis types with unique behavior should use an `IndexStyle` trait:

```julia
to_index(axis::MyAxisType, arg) = to_index(IndexStyle(axis), axis, arg)
to_index(::MyIndexStyle, axis, arg) = ...
```
"""
@propagate_inbounds to_index(axis, arg) = to_index(IndexStyle(axis), axis, arg)
to_index(axis, arg::CartesianIndices{0}) = arg
# Colons get converted to slices by `indices`
to_index(::IndexStyle, axis, ::Colon) = indices(axis)
@propagate_inbounds function to_index(::IndexStyle, axis, arg::Integer)
    @boundscheck checkbounds(axis, arg)
    return Int(arg)
end
@propagate_inbounds function to_index(::IndexStyle, axis, arg::AbstractArray{Bool})
    @boundscheck checkbounds(axis, arg)
    return @inbounds(axis[arg])
end
@propagate_inbounds function to_index(::IndexStyle, axis, arg::AbstractArray{I}) where {I<:Integer}
    @boundscheck if !checkindex(Bool, axis, arg)
        throw(BoundsError(axis, arg))
    end
    return arg
end
@propagate_inbounds function to_index(::IndexStyle, axis, arg::AbstractRange{I}) where {I<:Integer}
    @boundscheck if !checkindex(Bool, axis, arg)
        throw(BoundsError(axis, arg))
    end
    return arg
end
function to_index(S::IndexStyle, axis, arg::Any)
    throw(ArgumentError("invalid index: IndexStyle $S does not support indices of type $(typeof(arg))."))
end

"""
    unsafe_reconstruct(A, data; kwargs...)

Reconstruct `A` given the values in `data`. New methods using `unsafe_reconstruct`
should only dispatch on `A`.
"""
function unsafe_reconstruct(A::OneTo, data; kwargs...)
    if can_change_size(A)
        return typeof(A)(data)
    else
        if data isa Slice || !(known_length(A) === nothing || known_length(A) !== known_length(data))
            return A
        else
            return OneTo(data)
        end
    end
end
function unsafe_reconstruct(A::UnitRange, data; kwargs...)
    if can_change_size(A)
        return typeof(A)(data)
    else
        if data isa Slice || !(known_length(A) === nothing || known_length(A) !== known_length(data))
            return A
        else
            return UnitRange(data)
        end
    end
end
function unsafe_reconstruct(A::OptionallyStaticUnitRange, data; kwargs...)
    if can_change_size(A)
        return typeof(A)(data)
    else
        if data isa Slice || !(known_length(A) === nothing || known_length(A) !== known_length(data))
            return A
        else
            return OptionallyStaticUnitRange(data)
        end
    end
end
function unsafe_reconstruct(A::AbstractUnitRange, data; kwargs...)
    return static_first(data):static_last(data)
end

"""
    to_axes(A, inds)
    to_axes(A, old_axes, inds) -> new_axes

Construct new axes given the corresponding `inds` constructed after
`to_indices(A, old_axes, args) -> inds`. This method iterates through each
pair of axes and indices calling [`to_axis`](@ref).
"""
@inline function to_axes(A, inds::Tuple)
    if ndims(A) === 1
        return (to_axis(axes(A, 1), first(inds)),)
    elseif is_linear_indexing(A, inds)
        return (to_axis(eachindex(IndexLinear(), A), first(inds)),)
    else
        return to_axes(A, axes(A), inds)
    end
end
to_axes(A, ::Tuple{Ax,Vararg{Any}}, ::Tuple{}) where {Ax} = ()
to_axes(A, ::Tuple{}, ::Tuple{}) = ()
@propagate_inbounds function to_axes(A, axs::Tuple{Ax,Vararg{Any}}, inds::Tuple{I,Vararg{Any}}) where {Ax,I}
    N = argdims(A, I)
    if N === 0
        # drop this dimension
        return to_axes(A, tail(axs), tail(inds))
    elseif N === 1
        return (to_axis(first(axs), first(inds)), to_axes(A, tail(axs), tail(inds))...)
    else
        # Only multidimensional AbstractArray{Bool} and AbstractVector{CartesianIndex{N}}
        # make it to this point. They collapse several dimensions into one.
        axes_front, axes_tail = Base.IteratorsMD.split(axs, Val(N))
        return (to_multi_axis(IndexStyle(A), axes_front, first(inds)), to_axes(A, axes_tail, tail(inds))...)
    end
end

"""
    to_axis(old_axis, index) -> new_axis

Construct an `new_axis` for a newly constructed array that corresponds to the
previously executed `to_index(old_axis, arg) -> index`. `to_axis` assumes that
`index` has already been confirmed to be in bounds. The underlying indices of
`new_axis` begins at one and extends the length of `index` (i.e., one-based indexing).
"""
@inline function to_axis(axis, inds)
    if !can_change_size(axis) && (known_length(inds) !== nothing && known_length(axis) === known_length(inds))
        return axis
    else
        return to_axis(IndexStyle(axis), axis, inds)
    end
end
@inline function to_axis(S::IndexStyle, axis, inds)
    return unsafe_reconstruct(axis, StaticInt(1):static_length(inds))
end
@inline function to_multi_axis(::IndexStyle, axs::Tuple, inds)
    return to_axis(eachindex(LinearIndices(axs)), inds)
end

"""
    ArrayInterface.getindex(A, args...)

Retrieve the value(s) stored at the given key or index within a collection. Creating
another instance of `ArrayInterface.getindex` should only be done by overloading `A`.
Changing indexing based on a given argument from `args` should be done through
[`flatten_args`](@ref), [`to_index`](@ref), or [`to_axis`](@ref).
"""
@propagate_inbounds getindex(A, args...) = unsafe_getindex(A, to_indices(A, args))

"""
    unsafe_getindex(A, inds)

Indexes into `A` given `inds`. This method assumes that `inds` have already been
bounds-checked.
"""
unsafe_getindex(A, inds) = unsafe_getindex(UnsafeIndex(A, inds), A, inds)
unsafe_getindex(::UnsafeGetElement, A, inds) = unsafe_get_element(A, inds)
unsafe_getindex(::UnsafeGetCollection, A, inds) = unsafe_get_collection(A, inds)

"""
    unsafe_get_element(A::AbstractArray{T}, inds::Tuple) -> T

Returns an element of `A` at the indices `inds`. This method assumes all `inds`
have been checked for being in bounds. Any new array type using `ArrayInterface.getindex`
must define `unsafe_get_element(::NewArrayType, inds)`.
"""
function unsafe_get_element(A, inds)
    throw(MethodError(unsafe_getindex, (A, inds)))
end
function unsafe_get_element(A::Array, inds)
    if length(inds) === 0
        return Base.arrayref(false, A, 1)
    elseif inds isa Tuple{Vararg{Int}}
        return Base.arrayref(false, A, inds...)
    else
        throw(MethodError(unsafe_get_element, (A, inds)))
    end
end
@inline function unsafe_get_element(A::LinearIndices, inds)
    return Int(Base._to_linear_index(A, inds...))
end
@inline function unsafe_get_element(A::CartesianIndices, inds)
    return CartesianIndex(Base._to_subscript_indices(A, inds...))
end

# This is based on Base._unsafe_getindex from https://github.com/JuliaLang/julia/blob/c5ede45829bf8eb09f2145bfd6f089459d77b2b1/base/multidimensional.jl#L755.
"""
    unsafe_get_collection(A, inds)

Returns a collection of `A` given `inds`. `inds` is assumed to have been bounds-checked.
"""
function unsafe_get_collection(A, inds)
    axs = to_axes(A, inds)
    dest = similar(A, axs)
    if map(Base.unsafe_length, axes(dest)) == map(Base.unsafe_length, axs)
        Base._unsafe_getindex!(dest, A, inds...) # usually a generated function, don't allow it to impact inference result
    else
        Base.throw_checksize_error(dest, axs)
    end
    return dest
end

can_preserve_indices(::Type{T}) where {T<:AbstractRange} = known_step(T) === 1
can_preserve_indices(::Type{T}) where {T<:Int} = true
can_preserve_indices(::Type{T}) where {T} = false

_ints2range(x::Integer) = x:x
_ints2range(x::AbstractRange) = x

# if linear indexing on multidim or can't reconstruct AbstractUnitRange
# then construct Array of CartesianIndex/LinearIndices
@generated function can_preserve_indices(::Type{T}) where {T<:Tuple}
    for index_type in T.parameters
        can_preserve_indices(index_type) || return false
    end
    return true
end

@inline function unsafe_get_collection(A::CartesianIndices{N}, inds) where {N}
    if (length(inds) === 1 && N > 1) || !can_preserve_indices(typeof(inds))
        return Base._getindex(IndexStyle(A), A, inds...)
    else
        return CartesianIndices(to_axes(A, _ints2range.(inds)))
    end
end
@inline function unsafe_get_collection(A::LinearIndices{N}, inds) where {N}
    if is_linear_indexing(A, inds)
        return @inbounds(eachindex(A)[first(inds)])
    elseif can_preserve_indices(typeof(inds))
        return LinearIndices(to_axes(A, _ints2range.(inds)))
    else
        return Base._getindex(IndexStyle(A), A, inds...)
    end
end

"""
    ArrayInterface.setindex!(A, args...)

Store the given values at the given key or index within a collection.
"""
@propagate_inbounds function setindex!(A, val, args...)
    if can_setindex(A)
        return unsafe_setindex!(A, val, to_indices(A, args))
    else
        error("Instance of type $(typeof(A)) are not mutable and cannot change " *
              "elements after construction.")
    end
end

"""
    unsafe_setindex!(A, val, inds::Tuple)

Sets indices (`inds`) of `A` to `val`. This method assumes that `inds` have already been
bounds-checked. This step of the processing pipeline can be customized by:
"""
unsafe_setindex!(A, val, inds::Tuple) = unsafe_setindex!(UnsafeIndex(A, inds), A, val, inds)
unsafe_setindex!(::UnsafeGetElement, A, val, inds::Tuple) = unsafe_set_element!(A, val, inds)
unsafe_setindex!(::UnsafeGetCollection, A, val, inds::Tuple) = unsafe_set_collection!(A, val, inds)

"""
    unsafe_set_element!(A, val, inds::Tuple)

Sets an element of `A` to `val` at indices `inds`. This method assumes all `inds`
have been checked for being in bounds. Any new array type using `ArrayInterface.setindex!`
must define `unsafe_set_element!(::NewArrayType, val, inds)`.
"""
function unsafe_set_element!(A, val, inds)
    throw(MethodError(unsafe_set_element!, (A, val, inds)))
end
function unsafe_set_element!(A::Array{T}, val, inds::Tuple) where {T}
    if length(inds) === 0
        return Base.arrayset(false, A, convert(T, val)::T, 1)
    elseif inds isa Tuple{Vararg{Int}}
        return Base.arrayset(false, A, convert(T, val)::T, inds...)
    else
        throw(MethodError(unsafe_set_element!, (A, inds)))
    end
end

# This is based on Base._unsafe_setindex!.
"""
    unsafe_set_collection!(A, val, inds)

Sets `inds` of `A` to `val`. `inds` is assumed to have been bounds-checked.
"""
@inline function unsafe_set_collection!(A, val, inds)
    return Base._unsafe_setindex!(IndexStyle(A), A, val, inds...)
end

ext: jl


text:
#############################################################################
# minimum.jl
# Compute the minimum value of an array.
# All expressions and atoms are subtpyes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################
import Base.minimum
export minimum

### Minimum Atom
struct MinimumAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int, Int}

    function MinimumAtom(x::AbstractExpr)
        if sign(x) == ComplexSign()
            error("Argument should be real instead it is $(sign(x))")
        else
            children = (x,)
            return new(:minimum, hash(children), children, (1, 1))
        end
    end
end

function sign(x::MinimumAtom)
    return sign(x.children[1])
end

# The monotonicity
function monotonicity(x::MinimumAtom)
    return (Nondecreasing(),)
end

# If we have h(x) = f o g(x), the chain rule says h''(x) = g'(x)^T f''(g(x))g'(x) + f'(g(x))g''(x);
# this represents the first term
function curvature(x::MinimumAtom)
    return ConcaveVexity()
end

function evaluate(x::MinimumAtom)
    return Base.minimum(evaluate(x.children[1]))
end

# x >= this if minimum(x) = this
# so, x - this will be in the :NonNeg cone
function conic_form!(x::MinimumAtom, unique_conic_forms::UniqueConicForms=UniqueConicForms())
    if !has_conic_form(unique_conic_forms, x)
        this = Variable()
        objective = conic_form!(this, unique_conic_forms)
        conic_form!(this <= x.children[1], unique_conic_forms)
        cache_conic_form!(unique_conic_forms, x, objective)
    end
    return get_conic_form(unique_conic_forms, x)
end

minimum(x::AbstractExpr) = MinimumAtom(x)


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni.
#

function print_JSON_render_animation(
          func::Function,
          vec_variables::Vector{String},
          iterable::Any,
          scenefile::String,
          algorithm::Renderer,
     	camera_type::Union{String, Nothing},
		camera_position::Union{Point, Vec, Nothing},
     	α::Float64,
     	width::Int64,
     	height::Int64,
          a::Float64,
          γ::Float64,
          lum::Union{Number, Nothing},
     	anim_output::String,
          samples_per_pixel::Int64,
          declare_float::Union{Dict{String, Float64}, Nothing},
          time_of_start::String,
          rendering_time_s::Float64,
     )

     dict_renderer = if typeof(algorithm) == OnOffRenderer
              Dict(
                    "algorithm" => "On-Off Renderer",
                    "background color" => algorithm.background_color,
                    "color" => algorithm.color,
               )

          elseif typeof(algorithm) == FlatRenderer
              Dict(
                    "algorithm" => "Flat Renderer",
                    "background color" => algorithm.background_color,
               )
          elseif typeof(algorithm) == PathTracer
              Dict(
                    "algorithm" => "Path-Tracing Renderer",
                    "background color" => algorithm.background_color,
                    "PCG" => Dict(
                         "initial state" => get_state(algorithm.pcg),
                         "initial sequence" => get_inc(algorithm.pcg),
                         ),
                    "number of rays" => algorithm.num_of_rays,
                    "max depth" => algorithm.max_depth,
                    "russian roulette limit" => algorithm.russian_roulette_limit,
               )
          elseif typeof(algorithm) == PointLightRenderer
               Dict(
                    "algorithm" => "Point-Light Renderer",
                    "background color" => algorithm.background_color,
                    "color" => algorithm.ambient_color,
                    "dark parameter" => algorithm.dark_parameter,
               )
          end

     tonemapping = Dict(
          "normalization" => a,
          "gamma"=> γ,
          "input average luminosity" => lum,
     )

     data = Dict(
          "function"=>string(func),
          "vector of variables" => vec_variables,
          "iterable" => string(iterable),
          "scene file" => scenefile,
          "time of start" => time_of_start,
          "animation output" => anim_output,
          "renderer" => dict_renderer,
          "samples per pixel (0 means no antialiasing)" => samples_per_pixel,
          "alpha" => α,
          "width" => width,
          "height" => height,
          "tonemapping" => tonemapping,
          "declared float from Command Line" => declare_float,
          "rendering time (in s)"=> @sprintf("%.3f", rendering_time_s),
     )

     open( join(map(x->x*".", split(anim_output,".")[1:end-1])) * "json","w") do f
          JSON.print(f, data, 4)
     end
end

function render_animation(x::(Pair{T1,T2} where {T1,T2})...)
	render_animation( parse_render_animation_settings(  Dict( pair for pair in [x...]) )... )
end

function render_animation(
          func::Function,
          vec_variables::Vector{String},
          iterable::Any,
          scenefile::String,
          renderer_model::Renderer = FlatRenderer(),
     	camera_type::Union{String, Nothing} = nothing,
		camera_position::Union{Point, Vec, Nothing} = nothing,
     	α::Float64 = 0.,
     	width::Int64 = 640,
     	height::Int64 = 480,
          a::Float64 = 0.18,
          γ::Float64 = 1.27,
          lum::Union{Number, Nothing} = nothing,
     	anim_output::String = "scene_animation.mp4",
          samples_per_pixel::Int64 = 0,
		bool_print::Bool = true,
          declare_float::Union{Dict{String,Float64}, Nothing} = nothing,
          ONLY_FOR_TESTS::Bool = false,
     )

     check_is_iterable(iterable) || throw(ArgumentError("the input $(iterable) is not iterable"))
     check_is_iterable(iterable, Number) || throw(ArgumentError("the input iterable $(iterable) does not contains numbers"))

     iterable_float = try
          convert.(Float64, iterable)
     catch
          throw(ArgumentError("the input iterable $(iterable) does not contains numbers convertable to Float64"))
     end

     hasmethod(func, Tuple{Float64}) || throw(ArgumentError("function $(func) does not have a method for Tuple{Float64})"))
     length(func(1.0)) == length(vec_variables) ||
          throw(ArgumentError("length of vec_variables $(length(vec_variables)) and func return $(length(func(1.0))) do not match"))

     (bool_print==true) && println("\n\nStart the reading of \"$(scenefile)\"...")

     scene_model = open(scenefile, "r") do stream
          if isnothing(declare_float)
               inputstream = InputStream(stream, scenefile)
               parse_scene(inputstream)
          else
               inputstream = InputStream(stream, scenefile)
               parse_scene(inputstream, declare_float)
          end
     end

     for name ∈ vec_variables
          (name ∈ keys(scene_model.float_variables)) || throw(ArgumentError("$(name) is not a float identifier defined in $(scenefile)"))
     end

     (bool_print==true) && println("\nReaded and parsed \"$(scenefile)\", now start the animation rendering...\n")

     (ONLY_FOR_TESTS==false) || (return nothing)
     time_of_start = Dates.format(now(), DateFormat("Y-m-d : H:M:S"))
     time_1 = time()

     run(`rm -rf .wip_animation`)
	run(`mkdir .wip_animation`)

	dict_gen = Dict(
               "scenefile"=>scenefile,
               "camera_type"=>camera_type,
			"camera_position"=>camera_position,
               "alpha"=>α,
			"width"=>width,
			"height"=>height,
               "normalization"=>a,
               "gamma"=>γ,
               "avg_lum"=>lum,
			"samples_per_pixel"=>samples_per_pixel,
			"set_pfm_name"=>".wip_animation/scene.pfm",
			"bool_print"=>false,
			"bool_savepfm"=>false,
			"ONLY_FOR_TESTS"=>ONLY_FOR_TESTS,
			)

     N = length(iterable_float)
     iter = Progress(N, "Frame generated: ")
     Threads.@threads for index in 1:N
          value = iterable_float[index]
          values = func(value)
          dict = Dict(x=>y for (x,y) in zip(vec_variables, values))
          new_declare_float = isnothing(declare_float) ? dict : merge(dict, declare_float)

		NNN = @sprintf "%03d" index
		dict_spec = Dict(
                         "renderer" => copy(renderer_model),
					"set_png_name"=>".wip_animation/image$(NNN).png",
                         "declare_float"=>new_declare_float,
					)
		render(parse_render_settings(merge(dict_gen, dict_spec))...)
          next!(iter)
	end

     time_2 = time()
     rendering_time_s = time_2 - time_1

     run(`ffmpeg -r 25 -f image2 -s $(width)x$(height) -i
	.wip_animation/image%03d.png -vcodec libx264
	-pix_fmt yuv420p $(anim_output)`)

	run(`rm -rf .wip_animation`)

     print_JSON_render_animation(
          func,
          vec_variables,
          iterable,
          scenefile,
          copy(renderer_model),
     	camera_type,
		camera_position,
     	α,
     	width,
     	height,
          a, γ, lum,
     	anim_output,
          samples_per_pixel,
          declare_float,
          time_of_start,
          rendering_time_s,
     )

     name_json = join(map(x->x*".", split(anim_output,".")[1:end-1])) * "json"
     (bool_print==true) && println("\nJSON file \"$(name_json)\" correctly created.")
     (bool_print==true) && println("\nEND OF RENDERING\n")
end

"""
	render_animation(
			func::Function,
               vec_variables::Vector{String},
               iterable::Any,
               scenefile::String,
               renderer_model::Renderer = FlatRenderer(),
               camera_type::Union{String, Nothing} = nothing,
               camera_position::Union{Point, Vec, Nothing} = nothing,
               α::Float64 = 0.,
               width::Int64 = 640,
               height::Int64 = 480,
               a::Float64 = 0.18,
               γ::Float64 = 1.27,
               lum::Union{Number, Nothing} = nothing,
               anim_output::String = "scene_animation.mp4",
               samples_per_pixel::Int64 = 0,
               bool_print::Bool = true,
               declare_float::Union{Dict{String,Float64}, Nothing} = nothing,
               ONLY_FOR_TESTS::Bool = false,
		)

	render_animation(x::(Pair{T1,T2} where {T1,T2})...) =
		render_animation( parse_render_animation_settings(  Dict( pair for pair in [x...]) )... )

Render the input `scenefile` as an animation with the specified options, and creates the following
three files:
- the animation (`scene_animation.mp4` is the default name, if none is specified from the command line)
- the JSON file (which has the same name of the animation and `.json` estention, so
  `scene_animation.json` is the default name), that saves some datas about input commands, rendering time etc.

This function works following this steps:
- creates an hidden directory, called ".wip_animation"; if it already exists,
  it will be destroyed and recreated.
- inside ".wpi_animation", creates the png images of the rendered image (using the
  [`render`](@ref) function with the specified projection, renderer and image
  dims); each image correspons to a frame of the future animation
- through the `ffmpeg` software, the 360 png images are converted into the
  animation mp4 file, and saved in the main directory
- the ".wpi_animation" directory and all the png images inside it are destroyed

## Arguments

The following input arguments refers to the first method presented in the signature;
it's obviously very uncomfortable to use that method, so it's recommended to take
advantage of the second one, which allows to write the input values in a dictionary
like syntax with arbitrary order and comfort. See the documentation of
[`parse_render_animation_settings`](@ref) to learn how to use the keys:

- `func::Function` : function that takes as input the frame number and returns the values
  of the variables for that frame; it must have a method that takes only one input number,
  and must return a tuple of length equals to the `vec_variables` one.

- `vec_variables::Vector{String}` : vector that contains the variable names (DEFINED IN SCENEFILE)
  that will be overridden from frame to frame; its length must equals the tuple length returned
  by `func`.

- `iterable::Any` : an iterable object, that defines the frame numbers and the total number
  of frames; its values will be given in input to `func`.

- `renderer::Renderer = FlatRenderer()` : renderer to be used in the rendering, with all
  the settings already setted (exception made for the `world`, that will be overridden
  and created here)

- `camera_type::String = "per"` : set the perspective projection view:
  - `camera_type=="per"` -> set [`PerspectiveCamera`](@ref)  (default value)
  - `camera_type=="ort"`  -> set [`OrthogonalCamera`](@ref)

- `camera_position::Union{Point, Vec} = Point(-1.,0.,0.)` : set the point of observation
  in (`X`,`Y,`Z`) coordinates, as a `Point`or a `Vec` input object.

- `α::Float64 = 0.` : angle of rotation _*IN RADIANTS*_, relative to the vertical
  (i.e. z) axis with a right-handed rule convention (clockwise rotation for entering (x,y,z)-axis
  corresponds to a positive input rotation angle)

- `width::Int64 = 640` and `height::Int64 = 480` : pixel dimensions of the demo image;
  they must be both even positive integers.

- `a::Float64 = 0.18` : normalization scale factor for the tone mapping.

- `γ::Float64 = 1.27` : gamma factor for the tone mapping.

- `lum::Union{Number, Nothing} = nothing ` : average luminosity of the image; iIf not specified or equal to 0,
  it's calculated through [`avg_lum`](@ref)

- `anim_output::String = "scene_animation.mp4"` : name of the output animation file.

- `samples_per_pixel::Int64 = 0` : number of rays per pixel to be used (antialiasing);
  it must be a perfect square positive integer (0, 1, 4, 9, ...) and if is set to
  0 (default value) is choosen, no anti-aliasing occurs, and only one pixel-centered
  ray is fired for each pixel.

- `bool_print::Bool = true` : specifies if the WIP messages of the demo
  function should be printed or not.

- `declare_float::Union{Dict{String,Float64}, Nothing} = nothing` : an option (for the
  command line in particularly) to manually override the values of the float variables in
  the scene file; each overriden variable name (the key) is associated with its float value
  (i.e. `declare_float = Dict("var1"=>0.1, "var2"=>2.5)`)

- `ONLY_FOR_TESTS::Bool = false` : it's a bool variable conceived only to
  test the correct behaviour of the renderer for the input arguments; if set to `true`,
  no rendering is made!

See also: [`Point`](@ref), [`Vec`](@ref), [`Renderer`](@ref) [`OnOffRenderer`](@ref),
[`FlatRenderer`](@ref), [`PathTracer`](@ref), [`PointLightRenderer`](@ref),
[`parse_render_animation_settings`](@ref), [`render`](@ref)
"""
render_animation


text:
# The following code is a modified version of
#
#    function make_preduce_body(var, body)
#    function make_pfor_body(var, body)
#    macro distributed(args...)
#
# in https://github.com/JuliaLang/julia/blob/9f3265399227fbfc4f0160ec3592a9262bd3eb5f/stdlib/Distributed/src/macros.jl
#
# Its license is MIT: https://julialang.org/license

using Distributed: preduce, pfor, make_preduce_body, make_pfor_body

function my_make_preduce_body(var, prebody, body)
    quote
        function (reducer, R, lo::Int, hi::Int)
            $(esc(prebody))
            $(esc(var)) = R[lo]
            ac = $(esc(body))
            if lo != hi
                for $(esc(var)) in R[(lo+1):hi]
                    ac = reducer(ac, $(esc(body)))
                end
            end
            ac
        end
    end
end

function my_make_pfor_body(var, prebody, body)
    quote
        function (R, lo::Int, hi::Int)
            $(esc(prebody))
            for $(esc(var)) in R[lo:hi]
                $(esc(body))
            end
        end
    end
end

"""
    @my_distributed

A distributed memory, parallel for loop of the form:

```julia
@my_distributed begin
    prebody
end [reducer] for var = range
    body
end
```
"""
macro my_distributed(args...)
    na = length(args)
    if na==2
        prebody, loop = args
    elseif na==3
        prebody, reducer, loop = args
    else
        throw(ArgumentError("wrong number of arguments to @my_distributed"))
    end
    if !isa(loop,Expr) || loop.head !== :for
        error("malformed @my_distributed loop")
    end
    var = loop.args[1].args[1]
    r = loop.args[1].args[2]
    body = loop.args[2]
    if na==2
        syncvar = esc(Base.sync_varname)
        return quote
            local ref = pfor($(my_make_pfor_body(var, prebody, body)), $(esc(r)))
            if $(Expr(:islocal, syncvar))
                put!($syncvar, ref)
            end
            ref
        end
    else
        return :(preduce($(esc(reducer)), $(my_make_preduce_body(var, prebody, body)), $(esc(r))))
    end
end


text:
if !@isdefined(TabLevel)
    TabLevel = ""
end
println(TabLevel*"Open OpSqrt.jl")
TabLevel=TabLevel*"    "

using Primes

IntegerTypes = Union{Integer,Complex{Integer},Complex{Int64},
                     Rational,Complex{Rational{Integer}},
                     Complex{Rational{UInt64}},
                     Complex{Rational{Int64}},Complex{Bool}}

struct OpSqrt ###Define Single Majorana Element
    ###Each majorana is characterized by a string giving the current type
    scale::Complex{Rational{Int64}}
    sqrt::Integer ###This should be non-negative
    OpSqrt(x,y) = y<0 ? error("Number $y under square-root should be non-negative") : new(x,y)
end

function OpSqrt(x::Complex)
    if imag(x) != 0
        throw(DomainError("OpSqrt for complex x=$x is not implemented yet. In the future this should be the principal root of $x"))
    else
        OpSqrt(real(x))
    end
end

function OpSqrt(x::OpSqrt)
    xsq = OpSqrt(x.sqrt)
    if typeof(xsq) == OpSqrt
        return OpSqrt(x.scale*xsq.scale,xsq.sqrt)
    else
        return x.scale*xsq
    end
end

### Defining a special charcter function
√(x) = OpSqrt(x) ###Looks nice but maybe not so useful

OpSqrt(x::AbstractFloat) = sqrt(x)
function OpSqrt(x::Integer)
    if x < 0
        throw(DomainError("OpSqrt for x=$x < 0 is not implemented yet. In the future this should be $im"))
    elseif x == 0
        return 0
    elseif x == 1
        return 1
    end
    ###Work out the powers
    ###Collect the factors
    Factors=collect(factor(x))
    #println("OpSqrt($x)")
    #println("Factors:" ,Factors)
    Scale = 1
    Sqrt = 1
    for (num,pow) in Factors
        #println("(num,pow) = " ,(num,pow))
        Scale *= num^div(pow,2)
        Sqrt  *= num^rem(pow,2)
        #println("(Scale,Sqrt) = " ,(Scale,Sqrt))
    end
    if Sqrt == 1
        return Scale
    else
        OpSqrt(Scale,Sqrt)
    end
end

function OpSqrt(x::Rational) ###How to handle rationals
    ###We use the convention that it is only the numerator that carry the square root. Thus we simply add an extra factor of the denominator at the  numerator and the denominator.
    return OpSqrt(x.num*x.den)//x.den
end

import Base.==
==(y::OpSqrt,x::Number) = (x == y)
function ==(x::Number,y::OpSqrt)
    ###Comparisson is made numerically to actual numbers
    x == y.scale*sqrt(y.sqrt)
end

function (==)(x::OpSqrt,y::OpSqrt)
    if x.scale == y.scale && x.sqrt == y.sqrt
        return true ###Both correct
    elseif x.scale != y.scale && x.sqrt == y.sqrt
        return false ###One Mismatch
    elseif x.scale == y.scale && x.sqrt != y.sqrt
        return false ###One Mismatch
    else  ## ( x.scale != y.scale && x.sqrt != y.sqrt )
        ###Both mismatch.. here the could be a reduction error
        return x.scale^2*x.sqrt == y.scale^2*y.sqrt
    end
end

import Base.*
function *(x::OpSqrt,y::OpSqrt)
    NewScale = (x.scale*y.scale)
    NewSqrt = OpSqrt(x.sqrt*y.sqrt)
    NewProd = NewScale*NewSqrt
    return NewProd
end

*(y::OpSqrt,x::IntegerTypes) = x*y
function *(x::IntegerTypes,y::OpSqrt)
    NewScale=x*y.scale
    OpSqrt(NewScale,y.sqrt)
end
*(x::AbstractFloat,y::OpSqrt) = x*y.scale*sqrt(y.sqrt)

import Base.//
//(x::OpSqrt,y::IntegerTypes) = x*(1//y)
//(y::Union{OpSqrt,IntegerTypes},x::OpSqrt) = y*(inv(x))

import Base.^
function ^(x::OpSqrt,y::Integer)
    if y < 0
        (1//OpSqrt)^(-y)
    elseif y==0
        return 1
    elseif y==1
        return 1*x ###To make sure it copied and reduced
    elseif y==2
        return x*x ###To make sure it copied and reduced
    else
        return x*(x^(y-1)) ###Recursively apply the product
    end
end

import Base.inv
inv(x::OpSqrt) = (1//(x.sqrt*x.scale))*OpSqrt(x.sqrt)

import Base.show
function Base.show(io::IO,x::OpSqrt)
    RatScale=x.scale*1//1
    if isreal(RatScale)
        RatScale = real(RatScale)
        if RatScale.num == 1
            Str = ""
        else
            Str = "$(RatScale.num)"
        end
        RatDenom = RatScale.den
    else ###This number if complex and has fractional structure
        RatDenom=RatScale.re.den*RatScale.im.den
        RatScale=(RatScale.re.num*RatScale.im.den+
                  im*RatScale.im.num*RatScale.re.den)
        if RatScale == 1
            Str = ""
        elseif RatScale == im
            Str = "im"
        elseif RatScale == -im
            Str = "-im"
        else
            Str="($(RatScale))"
        end
    end
    Str*="√($(x.sqrt))"
    if RatDenom == 1
        Str *= ""
    else
        Str *= "//$RatDenom"
    end
    print(io,Str)
end

import Base.isreal
isreal(x::OpSqrt) = isreal(x.scale)
import Base.real
real(x::OpSqrt) = OpSqrt(real(x.scale),x.sqrt)
import Base.imag
imag(x::OpSqrt) = OpSqrt(imag(x.scale),x.sqrt)

import Base.+
import Base.-
-(x::OpSqrt) = (-1)*x
-(x::IntegerTypes,y::OpSqrt)  = OpSqrt(x,1) - y
-(y::OpSqrt,x::IntegerTypes)  = y - OpSqrt(x,1)
-(x::OpSqrt,y::OpSqrt)  = x + (-y)
+(x::IntegerTypes,y::OpSqrt)  = OpSqrt(x,1) + y
+(y::OpSqrt,x::IntegerTypes)  = OpSqrt(x,1) + y

function +(x::OpSqrt,y::OpSqrt)
    ###Check for zeros
    if x==0
        return OpSqrt(y)
    elseif y==0
        return OpSqrt(x)
    end

    ###First make sure functions is reduced
    RedX = OpSqrt(x.sqrt)
    RedY = OpSqrt(y.sqrt)
    if typeof(RedX) == OpSqrt
        RSqX = RedX.sqrt
        RScX = RedX.scale
    else
        RSqX = 1
        RScX = RedX
    end
    if typeof(RedY) == OpSqrt
        RSqY = RedY.sqrt
        RScY = RedY.scale
    else
        RSqY = 1
        RScY = RedY
    end
    #println("RSqX:",RSqX)
    #println("RScX:",RScX)
    #println("RSqY:",RSqY)
    #println("RScY:",RScY)
    #println("x.scale:",x.scale)
    #println("y.scale:",y.scale)
    if RSqX != RSqY
        error("$x cannot be added to $y since the square-roots are different")
    else
        if RSqX == 1
            return x.scale*RScX+y.scale*RScY
        else
            return OpSqrt(x.scale*RScX+y.scale*RScY,RSqX)
        end
    end
end

import Base.conj
function conj(x::OpSqrt)
    cx=conj(x.scale)*OpSqrt(x.sqrt)
    return cx
end

import Base.isless
function isless(y::OpSqrt,x::Number)
    if isreal(y.scale)
        x > real(y.scale)*sqrt(y.sqrt)
    else
        error("$y is complex. Comparison makes no sense")
    end
end
function isless(x::Number,y::OpSqrt)
    if isreal(y.scale)
        x < real(y.scale)*sqrt(y.sqrt)
    else
        error("$y is complex. Comparison makes no sense")
    end
end
function isless(x::OpSqrt,y::OpSqrt)
    if !isreal(x.scale)
        error("$x is complex. Comparison makes no sense")
    elseif !isreal(y.scale)
        error("$y is complex. Comparison makes no sense")
    else
        real(x.scale)*sqrt(x.sqrt) < real(y.scale)*sqrt(y.sqrt)
    end
end

TabLevel=TabLevel[1:end-4]
println(TabLevel*"Close OpSqrt.jl")


text:
# nufft.jl

using MIRT: Anufft, nufft_init, nufft_errors
using MIRT: dtft_init
import MIRT: nufft_eltype

#using MIRT: dtft_init, map_many
#include("../utility/map_many.jl")
using NFFT
using Plots
using LinearAlgebra: norm
using LinearMapsAA: LinearMapAM, LinearMapAO
using Test: @test, @testset, @test_throws

"""
nufft_test1( ; M=30, N=20, n_shift=1.7, T=?, tol=?)
simple 1D tests
"""
function nufft_test1( ;
	M::Int = 30, N::Int = 20, n_shift::Real = 1.7,
	T::DataType = Float64, tol::Real = 2e-6,
)

	w = (rand(M) .- 0.5) * 2 * pi
	w = T.(w)
	x = randn(complex(T), N)
	sd = dtft_init(w, N ; n_shift=n_shift)
	sn = nufft_init(w, N ; n_shift=n_shift, operator=false)
	o0 = sd.dtft(x)
	y = Complex{T}.(o0 / norm(o0))
	a0 = sd.adjoint(y)
	o1 = sn.nufft(x)
	o2 = sn.A * x
	a1 = sn.adjoint(y)
	a2 = sn.A' * y
	@test norm(o1 - o0, Inf) / norm(o0, Inf) < tol
	@test isequal(o1, o2)
	@test norm(a1 - a0, Inf) / norm(a0, Inf) < tol
	@test isequal(a1, a2)
	@test Matrix(sn.A)' ≈ Matrix(sn.A') # 1D adjoint test

	A = sn.A
	@test A.name == "nufft1"
	@test A isa LinearMapAM

	sn = nufft_init(w, N, n_shift=n_shift, do_many=false, operator=true)
	o3 = sn.nufft(x)
	@test norm(o3 - o0, Inf) / norm(o0, Inf) < tol

	B = sn.A
	@test B isa LinearMapAO

	A = Anufft(w, N, n_shift=n_shift)
	@test A isa LinearMapAO
	@test A * x == o3

	A = Anufft(w', N, n_shift=n_shift) # w shape
	@test A._odim == (1,M)
	@test A * x == transpose(o3)

	sn.nufft(ones(Int,N)) # produce a "conversion" warning
	true
end

"""
nufft_test2( ; M=?, N=?, n_shift=?, T=?, tol=?)
simple 2D test
"""
function nufft_test2( ;
	M::Int = 35,
	N::Dims = (10,8),
	n_shift::AbstractVector{<:Real} = [4,3],
	T::DataType = Float64, tol::Real = 2e-6,
)

	w = []

#=
	# fft sampling
	M = prod(N)
	w1 = (2*pi) * (0:N[1]-1) / N[1]
    w2 = (2*pi) * (0:N[2]-1) / N[2]
    w1 = repeat(w1, 1, N[2])
    w2 = repeat(w2', N[1], 1)
    w = [vec(w1) vec(w2)]
	n_shift = [0,0]
=#

	w = (rand(M,2) .- 0.5) * 2 * pi

	w = T.(w)
	sd = dtft_init(w, N ; n_shift=n_shift)
	sn = nufft_init(w, N ; n_shift=n_shift, pi_error=false)

	x = randn(Complex{T}, N)
	o0 = sd.dtft(x)
	o1 = sn.nufft(x)
	@test norm(o1 - o0, Inf) / norm(o0, Inf) < tol

#=
	# fft test only
	o2 = fft(x)
	o0 = reshape(o0, N)
	o1 = reshape(o1, N)
	@show norm(o2 - o0, Inf) / norm(o0, Inf)
	@show norm(o1 - o0, Inf) / norm(o0, Inf)
	@show norm(o2 - o1, Inf) / norm(o0, Inf)
=#

	y = convert(Array{Complex{T}}, o0 / norm(o0))
	a0 = sd.adjoint(y)
	a1 = sn.adjoint(y)
	@test norm(a1 - a0, Inf) / norm(a0, Inf) < tol
	o2 = sn.A * x
	@test isequal(o1, o2)
	a2 = sn.A' * y
	@test isequal(a1, a2)
	@test Matrix(sn.A)' ≈ Matrix(sn.A') # 2D adjoint test

	@test isequal(sn.nufft(cat(dims=4, x, 2x)),
		cat(dims=3, sn.nufft(x), sn.nufft(2x)))
	@test isequal(sn.adjoint(cat(dims=3, y, 2y)),
			cat(dims=4, sn.adjoint(y), sn.adjoint(2y)))

	A = sn.A
	@test A.name == "nufft2"
	@test A.N == N

	Ao = nufft_init(w, N ; n_shift=n_shift, pi_error=false, operator=true).A
	Am = nufft_init(w, N ; n_shift=n_shift, pi_error=false, operator=false).A
	@test Ao * x == Am * vec(x)
	y = Ao * x
	@test Ao'*y == reshape(Am'*y, N)

	sn = nufft_init(w, N ; n_shift=n_shift, pi_error=false,
		do_many=false, operator=true)
	o3 = sn.nufft(x)
	@test norm(o3 - o0, Inf) / norm(o0, Inf) < tol

	A = Anufft(w, N, n_shift=n_shift)
	@test A isa LinearMapAO
	@test A * x == o3

	if M == 35
		odim = (7, 5) # test array output
		w = reshape(w, odim..., 2)
		A = Anufft(w, N, n_shift=n_shift)
		@test A._odim == odim
		@test A * x == reshape(o3, odim)
	end

	true
end

"""
nufft_plot1()

Plot worst-case error over all frequencies w between 0 and 2pi/N for various N.
"""
function nufft_plot1()
	Nlist = 2 .^ (4:9)
	elist = zeros(length(Nlist))
	for (ii, N) in enumerate(Nlist)
		w, errs = nufft_errors( ; N = N)
	#	plot(w*N/2/π, tmp)
		elist[ii] = maximum(errs)
	end
	scatter(Nlist, elist, xtick=Nlist, label="", xlabel="N", ylabel="error")
end

"""
nufft_plot_error_m( ; mlist=?)

Plot error vs NFFT sigma
"""
function nufft_plot_error_m(;
	mlist::AbstractArray{<:Int} = 3:7)
	worst = zeros(length(mlist))
	for (jm,nfft_m) = enumerate(mlist)
		_, errs = nufft_errors( ; nfft_m=nfft_m)
		worst[jm] = maximum(errs)
	end
	scatter(mlist, worst, xlabel="m", ylabel="error", label="")
end

"""
nufft_plot_error_s( ; slist=?)

Plot error vs NFFT sigma
"""
function nufft_plot_error_s( ;
	slist::AbstractArray{<:Real} = [1.5; 2:6])
	worst = zeros(length(slist))
	for (is,σ) in enumerate(slist)
		_, errs = nufft_errors(; nfft_sigma=σ)
		worst[is] = maximum(errs)
	end
	scatter(slist, worst, xlabel="σ", ylabel="error", label="")
end

"""
    nufft_plots()
various NUFFT error plots
"""
function nufft_plots()
	p1 = nufft_plot_error_s()
	p2 = nufft_plot_error_m()
	p3 = nufft_plot1()
	plot(p1, p2, p3)
end

@testset "basics" begin
	@test nufft_eltype(Bool) === Float32
	@test nufft_eltype(Float16) === Float32
	@test nufft_eltype(Float64) === Float64
	@test_throws String nufft_eltype(BigFloat)
	@test_throws String nufft_init([0], 2) # small
	@test_throws String nufft_init([0], 7) # odd
	@test_throws ArgumentError nufft_init([2π], 8) # π
end

@testset "1D" begin
	@test nufft_test1()
	@test nufft_test1(; T=Float32)
end

@testset "2D" begin
	@test nufft_test2()
	@test nufft_test2(; T=Float32)
end

@testset "plots" begin
	@test nufft_plots() isa Plots.Plot
end

#=
	# todo: 1d vs 2d
	M = 4
	N = (M,1)
	w = (0:(M-1))/M * 2 * pi
	w = [w zeros(M)]
	sd = dtft_init(w, N)
	Ad = Matrix(sd.A)
	sn1 = nufft_init(w[:,1], N[1]; pi_error=false)
	An1 = Matrix(sn1.A)
	@show maximum(abs.(An1 - Ad))
	sn2 = nufft_init(w, N; pi_error=false)
	An2 = Matrix(sn2.A)
#	@show o0 = sd.dtft(x)
#	@show o1 = sn.nufft(x)
#	@show o1-o0
#	@test maximum(abs.(o1 - o0)) < tol
=#

#=
	# todo MWE for 1D vs 2D
	M = 6
	x1 = collect((-Int(M/2)):(Int(M/2)-1))/M
	N1 = M
	p1 = NFFTPlan(x1, N1)
#	p1 = NFFTPlan(x1', (N1,)) # this works too
#	f1 = ones(N1)
	f1 = [1.,2,3,4,0,0]
	o1 = nfft(p1, complex(f1)) # slightly annoying to have to use complex() here

	N2 = (M,4) # works fine with "4" instead of "1" here and separable signal
	x2 = [x1 zeros(M)]
	p2 = NFFTPlan(x2', N2)
	f2 = f1 * [0,0,1,0]'
	o2 = nfft(p2, complex(f2)) # ditto

	display(round.(o1, digits=7))
	display(round.(o2, digits=7))

#	tmp = nufft_init(x1*2*pi, N1)
=#

#nufft(:test)
#nufft_test1(;N=41) # todo: inaccurate


text:
using BandedMatrices, LinearAlgebra, LazyArrays, Test
import BandedMatrices: banded_mul!, isbanded, AbstractBandedLayout, BandedStyle

struct PseudoBandedMatrix{T} <: AbstractMatrix{T}
    data::Array{T}
    l::Int
    u::Int
end

Base.size(A::PseudoBandedMatrix) = size(A.data)
function Base.getindex(A::PseudoBandedMatrix, j::Int, k::Int)
    l, u = bandwidths(A)
    if -l ≤ k-j ≤ u
        A.data[j, k]
    else
        zero(eltype(A.data))
    end
end
function Base.setindex!(A::PseudoBandedMatrix, v, j::Int, k::Int)
    l, u = bandwidths(A)
    if -l ≤ k-j ≤ u
        A.data[j, k] = v
    else
        error("out of band.")
    end
end

struct PseudoBandedLayout <: AbstractBandedLayout end
Base.BroadcastStyle(::Type{<:PseudoBandedMatrix}) = BandedStyle()
BandedMatrices.MemoryLayout(::PseudoBandedMatrix) = PseudoBandedLayout()
BandedMatrices.isbanded(::PseudoBandedMatrix) = true
BandedMatrices.bandwidths(A::PseudoBandedMatrix) = (A.l , A.u)
BandedMatrices.inbands_getindex(A::PseudoBandedMatrix, j::Int, k::Int) = A.data[j, k]
BandedMatrices.inbands_setindex!(A::PseudoBandedMatrix, v, j::Int, k::Int) = setindex!(A.data, v, j, k)

@testset "banded matrix interface" begin
    @test isbanded(Zeros(5,6))
    @test bandwidths(Zeros(5,6)) == (0,0)
    @test BandedMatrices.inbands_getindex(Zeros(5,6), 1,2) == 0

    @test isbanded(Eye(5))
    @test bandwidths(Eye(5)) == (0,0)
    @test BandedMatrices.inbands_getindex(Eye(5), 1,1) == 1

    A = Diagonal(ones(5,5))
    @test isbanded(A)
    @test bandwidths(A) == (0,0)
    @test BandedMatrices.inbands_getindex(A, 1,1) == 1
    BandedMatrices.inbands_setindex!(A, 2, 1,1)
    @test A[1,1] == 2
    @test A[1,2] == 0
    @test BandedMatrices.@inbands(A[1,2]) == 2

    A = SymTridiagonal([1,2,3],[4,5])
    @test isbanded(A)
    @test bandwidths(A) == (1,1)
    @test BandedMatrices.inbands_getindex(A, 1,1) == 1
    BandedMatrices.inbands_setindex!(A, 2, 1,1)
    @test A[1,1] == 2

    A = PseudoBandedMatrix(rand(5, 4), 2, 2)
    B = rand(5, 4)
    C = copy(B)

    @test Matrix(B .= 2.0 .* A .+ B) ≈ 2*Matrix(A) + C

    A = PseudoBandedMatrix(rand(5, 4), 1, 2)
    B = (z -> exp(z)-1).(A)
    @test B isa BandedMatrix
    @test bandwidths(B) == bandwidths(A)
    @test B == (z -> exp(z)-1).(Matrix(A))

    A = PseudoBandedMatrix(rand(5, 4), 1, 2)
    B = A .* 2
    @test B isa BandedMatrix
    @test bandwidths(B) == bandwidths(A)
    @test B == 2Matrix(A) == (2 .* A)

    A = PseudoBandedMatrix(rand(5, 4), 1, 2)
    B = PseudoBandedMatrix(rand(5, 4), 2, 1)

    @test A .+ B isa BandedMatrix
    @test bandwidths(A .+ B) == (2,2)
    @test A .+ B == Matrix(A) + Matrix(B)

    A = PseudoBandedMatrix(rand(5, 4), 1, 2)
    B = PseudoBandedMatrix(rand(5, 4), 2, 3)
    C = deepcopy(B)
    @test Matrix(C .= 2.0 .* A .+ C) ≈ 2*Matrix(A) + B ≈ 2*A + Matrix(B) ≈ 2*Matrix(A) + Matrix(B) ≈ 2*A + B

    y = rand(4)
    z = zeros(5)
    z .= Mul(A, y)
    @test z ≈ A*y ≈ Matrix(A)*y

    B = PseudoBandedMatrix(rand(4, 4), 2, 3)
    C = PseudoBandedMatrix(zeros(5, 4), 3, 4)
    D = zeros(5, 4)

    @test (C .= Mul(A, B)) ≈ (D .= Mul(A, B)) ≈ A*B

    @test bandwidths(BandedMatrix(A)) ==
            bandwidths(BandedMatrix{Float64}(A)) ==
            bandwidths(BandedMatrix{Float64,Matrix{Float64}}(A)) ==
            bandwidths(convert(BandedMatrix{Float64}, A)) ==
            bandwidths(convert(BandedMatrix{Float64,Matrix{Float64}},A)) ==
            bandwidths(A)
end


text:
function find_memory_type(physical_device::PhysicalDevice, type_flag, properties::MemoryPropertyFlag)
    mem_props = get_physical_device_memory_properties(physical_device)
    indices =
        findall(
            x -> properties in x.property_flags,
            mem_props.memory_types[1:mem_props.memory_type_count],
        ) .- 1
    if isempty(indices)
        error("Could not find memory with properties $properties")
    else
        ind = findfirst(i -> type_flag & 1 << i ≠ 0, indices)
        if isnothing(ind)
            error("Could not find memory with type $type_flag")
        else
            indices[ind]
        end
    end
end

function Vulkan.DeviceMemory(device::Device, memory_requirements::MemoryRequirements, properties)
    DeviceMemory(
        device,
        memory_requirements.size,
        find_memory_type(device.physical_device, memory_requirements.memory_type_bits, properties),
    )
end

buffer_size(data::AbstractVector{T}) where {T} = sizeof(T) * length(data)
buffer_size(data) = sizeof(data)

"""
Upload data to the specified memory.

!!! warning
    The `memory` must be host coherent and host visible, otherwise the operation will fail.
"""
function upload_data(memory::DeviceMemory, data::DenseArray{T}) where {T}
    memptr = unwrap(map_memory(memory.device, memory, 0, buffer_size(data)))
    GC.@preserve data unsafe_copyto!(Ptr{T}(memptr), pointer(data), length(data))
    unwrap(unmap_memory(memory.device, memory))
end

"""
Download data from the specified memory to an `Array`.

If `copy` if set to true, then all the data mapped from `memory` will be copied. If not, care should be taken to preserve the `memory` mapped and valid as long as the returned data is in use.
"""
function download_data(::Type{<:DenseArray{T}}, memory::DeviceMemory, dims; offset=0, copy=true, unmap=true) where {T}
    size = sizeof(T) * prod(dims)
    memptr = unwrap(map_memory(memory.device, memory, offset, size))
    data = unsafe_wrap(Array, convert(Ptr{T}, memptr), dims; own=false)
    if unmap
        unwrap(unmap_memory(memory.device, memory))
    end
    if copy
        deepcopy(data)
    else
        data
    end
end

"""
Allocate a `DeviceMemory` object with the specified properties and bind it to the `buffer` using memory requirements from `get_buffer_memory_requirements`.
"""
function Vulkan.DeviceMemory(buffer::Buffer, properties::MemoryPropertyFlag) where {T}
    device = buffer.device
    memreqs = get_buffer_memory_requirements(device, buffer)
    memory = DeviceMemory(device, memreqs, properties)
    unwrap(bind_buffer_memory(device, buffer, memory, 0))
    memory
end

"""
Allocate a host visible and coherent `DeviceMemory` object, bind it to the `buffer` using memory requirements from `get_buffer_memory_requirements` and upload `data` to it.
"""
function Vulkan.DeviceMemory(buffer::Buffer, data::DenseArray{T}) where {T}
    memory = DeviceMemory(buffer, MEMORY_PROPERTY_HOST_VISIBLE_BIT | MEMORY_PROPERTY_HOST_COHERENT_BIT)
    upload_data(memory, data)
    memory
end

"""
Allocate a `DeviceMemory` object and bind it to the `image` using memory requirements from `get_image_memory_requirements`.
"""
function Vulkan.DeviceMemory(image::Image, properties::MemoryPropertyFlag)
    memreqs = get_image_memory_requirements(image.device, image)
    memory = DeviceMemory(image.device, memreqs, properties)
    unwrap(bind_image_memory(image.device, image, memory, 0))
    memory
end


text:
@testset "Stats" begin
    @testset "zscore_stats" begin
        a = collect(1:5)
        @test mean(a) == 3.0
        @test std(a) == sqrt(sum(abs2.(a .- mean(a)) / (length(a) - 1)))
    end

    @testset "zscore_df" begin
        a = 1:5
        df = DataFrame(A = a)
        zscore!(df, :A, :newA)
        @test df[!, :newA] == zscore(a)
    end
end

@testset "zscore" begin
    @testset "zscore_df" begin
        df = DataFrame(
            A = 1:4,
            B = 5:8,
            C = 9:12)
        zscore!(df, [:A, :B], [:nA, :nB])
        @test isapprox(df[!, :nA], [ -1.161895003862225,
            -0.3872983346207417,
             0.3872983346207417,
             1.161895003862225])
        @test isapprox(df[!, :nB], [ -1.161895003862225,
             -0.3872983346207417,
              0.3872983346207417,
              1.161895003862225])
    end

    @testset "unzscore" begin
        df = DataFrame(
            A = 1:4,
            B = 5:8,
            C = 9:12)

        mean_A, std_A = mean_and_std(df[!, :A])
        mean_B, std_B = mean_and_std(df[!, :B])

        zscore!(df, [:A, :B], [:nA, :nB])
        @test isapprox(unzscore(df[!, :nA], mean_A, std_A), df[!, :A])
        @test isapprox(unzscore(df[!, :nB], mean_B, std_B), df[!, :B])
    end
end

@testset "Min-Max Scaling" begin
    @testset "Check bounds" begin
        df = DataFrame(
            A = 1:4,
            B = 5:8,
            C = 9:12)

        target_min = 0.0
        target_max = 10.0

        minmax_scaling!(df, [:A, :B], [:nA, :nB], target_min, target_max)

        @test maximum(df[!, :nA]) <= target_max
        @test maximum(df[!, :nB]) <= target_max
        @test minimum(df[!, :nA]) >= target_min
        @test minimum(df[!, :nB]) >= target_min
    end

    @testset "Check bounds" begin
        df = DataFrame(
            A = 1:4,
            B = 5:8,
            C = 9:12)

        target_min = 0.0
        target_max = 10.0

        minmax_scaling!(df, [:A, :B], [:nA, :nB], target_min, target_max)

        @test isapprox(unminmax_scaling(df[!, :nA], 1, 4, target_min, target_max), df[!, :A])
        @test isapprox(unminmax_scaling(df[!, :nB], 5, 8, target_min, target_max), df[!, :B])
    end
end


text:
# Copyright 2016, 2017, 2018 Martin Holters
# See accompanying license file.

import Base.deleteat!
import Base.isless
import Base.isempty
using Compat.Statistics: mean, var

mutable struct KDTree{Tcv<:AbstractVector,Tp<:AbstractMatrix}
    cut_dim::Vector{Int}
    cut_val::Tcv
    ps_idx::Vector{Int}
    ps::Tp
end

function KDTree(p::AbstractMatrix, Np=size(p,2))
    function calc_cut_idx(min_idx, max_idx)
        N = max_idx - min_idx + 1
        N2 = 2^floor(Int, log2(N-1))
        if 3(N2÷2) ≤ N
            return min_idx+N2-1
        else
            return min_idx + N - N2÷2 - 1
        end
    end

    if Np == 0
        return KDTree{Vector{eltype(p)},typeof(p)}([], [], [], p)
    end

    min_idx = zeros(Int, Np-1)
    max_idx = zeros(Int, Np-1)
    cut_idx = zeros(Int, Np-1)
    cut_dim = zeros(Int, Np-1)
    cut_val = zeros(eltype(p), Np-1)

    if Np == 1
        return KDTree{typeof(cut_val),typeof(p)}(cut_dim, cut_val, [1], p)
    end

    @static if VERSION ≥ v"0.7.0-DEV.4064"
        dim = argmax(vec(var(p[:,1:Np], dims=2)))
    else
        dim = argmax(vec(Base.var(p[:,1:Np], 2)))
    end
    p_idx = sortperm(vec(p[dim,:]))

    min_idx[1] = 1
    max_idx[1] = Np
    cut_idx[1] = calc_cut_idx(min_idx[1], max_idx[1])
    cut_dim[1] = dim
    cut_val[1] = mean(p[dim, p_idx[cut_idx[1]:cut_idx[1]+1]])

    for n in 2:Np-1
        parent_n = n ÷ 2
        if n % 2 == 0
            min_idx[n] = min_idx[parent_n]
            max_idx[n] = cut_idx[parent_n]
        else
            min_idx[n] = cut_idx[parent_n]+1
            max_idx[n] = max_idx[parent_n]
        end
        @static if VERSION ≥ v"0.7.0-DEV.4064"
            dim = argmax(vec(var(p[:,p_idx[min_idx[n]:max_idx[n]]], dims=2)))
        else
            dim = argmax(vec(Base.var(p[:,p_idx[min_idx[n]:max_idx[n]]], 2)))
        end
        idx = sortperm(vec(p[dim,p_idx[min_idx[n]:max_idx[n]]]))
        p_idx[min_idx[n]:max_idx[n]] = p_idx[idx .+ min_idx[n] .- 1]
        cut_idx[n] = calc_cut_idx(min_idx[n], max_idx[n])
        cut_dim[n] = dim
        cut_val[n] = mean(p[dim, p_idx[cut_idx[n]:cut_idx[n]+1]])
    end

    p_idx_final = zeros(Int, Np)
    for n in 1:Np
        parent_n = (n+Np-1) ÷ 2;
        if (n+Np) % 2 == 1
            p_idx_final[n] = p_idx[min_idx[parent_n]]
        else
            p_idx_final[n] = p_idx[max_idx[parent_n]]
        end
    end

    return KDTree{typeof(cut_val),typeof(p)}(cut_dim, cut_val, p_idx_final, p)
end

mutable struct AltEntry{T}
    idx::Int
    delta::Vector{T}
    delta_norm::T
end

isless(e1::AltEntry, e2::AltEntry) = isless(e1.delta_norm, e2.delta_norm)

mutable struct Alts{T}
    entries::Vector{AltEntry{T}}
    best_dist::T
    best_pidx::Int
    number_valid::Int
end

Alts(p::Vector{T}) where {T} =
    Alts([AltEntry(1, Vector{T}(undef, length(p)), zero(T))], typemax(T), 0, 1)

function init!(alts::Alts{T}, best_dist, best_pidx) where {T}
    alts.number_valid = 1
    alts.entries[1].idx = 1
    fill!(alts.entries[1].delta, zero(T))
    alts.entries[1].delta_norm = zero(T)
    alts.best_dist = best_dist
    alts.best_pidx = best_pidx
end

function siftup!(alts::Alts, i)
    entries = alts.entries
    if i > alts.number_valid
        throw(BoundsError())
    end
    parent = i ÷ 2
    @inbounds while i > 1 && entries[i] < entries[parent]
        entries[i], entries[parent] = entries[parent], entries[i]
        i = parent
        parent = i ÷ 2
    end
end

function siftdown!(alts::Alts, i)
    if i < 1
        throw(BoundsError())
    end
    entries = alts.entries
    N = alts.number_valid
    @inbounds while true
        min = i
        if 2i ≤ N && entries[2i] < entries[min]
            min = 2i
        end
        if 2i+1 ≤ N && entries[2i+1] < entries[min]
            min = 2i+1
        end
        if min == i
            break
        end
        entries[i], entries[min] = entries[min], entries[i]
        i = min
    end
end

isempty(alts::Alts) = alts.number_valid == 0
peek(alts::Alts) = alts.entries[1]

function deleteat!(alts::Alts, i::Integer)
    alts.entries[i], alts.entries[alts.number_valid] = alts.entries[alts.number_valid], alts.entries[i]
    alts.number_valid -= 1
    if i ≤ alts.number_valid
        if i==1 || alts.entries[i] > alts.entries[i ÷ 2]
            siftdown!(alts, i)
        else
            siftup!(alts, i)
        end
    end
end

function dequeue!(alts::Alts)
    e = alts.entries[1]
    deleteat!(alts, 1)
    return e
end

function enqueue!(alts::Alts{T}, new_idx::Int, ref_delta::Vector{T},
                  delta_update_dim::Int, delta_update_val::T,
                  new_delta_norm::T) where {T}
    if alts.number_valid == length(alts.entries)
        delta = copy(ref_delta)
        delta[delta_update_dim] = delta_update_val
        push!(alts.entries, AltEntry(new_idx, delta, new_delta_norm))
    else
        alts.entries[alts.number_valid+1].idx = new_idx
        copyto!(alts.entries[alts.number_valid+1].delta, ref_delta)
        alts.entries[alts.number_valid+1].delta[delta_update_dim] = delta_update_val
        alts.entries[alts.number_valid+1].delta_norm = new_delta_norm
    end
    if alts.entries[alts.number_valid+1].delta_norm < alts.best_dist
        alts.number_valid += 1
        siftup!(alts, alts.number_valid)
    end
end

function update_best_dist!(alts, dist, p_idx)
    if dist < alts.best_dist
        alts.best_dist = dist
        alts.best_pidx = p_idx
        for i in alts.number_valid:-1:1
            if alts.entries[i].delta_norm ≥ alts.best_dist
                deleteat!(alts, i)
            end
        end
    end
end

indnearest(tree::KDTree, p::AbstractVector, alt = Alts(p)) =
    indnearest(tree, p, typemax(Int), alt)

function indnearest(tree::KDTree, p::AbstractVector, max_leaves::Int,
                    alt = Alts(p))
    if length(p) ≠ size(tree.ps, 1)
        throw(DimensionMismatch())
    end
    l = 0
    p_idx = 0
    while l < max_leaves && !isempty(alt)
        best_alt = dequeue!(alt)
        idx = best_alt.idx
        delta = best_alt.delta
        delta_norm = best_alt.delta_norm

        while idx ≤ length(tree.cut_dim)
            dim = tree.cut_dim[idx]
            new_alt_delta_norm = delta_norm - delta[dim]^2 + (p[dim] - tree.cut_val[idx])^2
            if new_alt_delta_norm < alt.best_dist
                new_alt_idx = p[dim] ≤ tree.cut_val[idx] ? 2idx+1 : 2idx;
                enqueue!(alt, new_alt_idx, delta, dim, p[dim] - tree.cut_val[idx], new_alt_delta_norm)
            end
            if p[dim] ≤ tree.cut_val[idx]
                idx *= 2
            else
                idx = 2idx + 1
            end
        end
        idx -= length(tree.cut_dim)

        p_idx = tree.ps_idx[idx]
        if p_idx < 1 || p_idx > size(tree.ps, 2)
            throw(BoundsError)
        end
        dist = 0.
        @simd for i in 1:length(p)
            @inbounds dist += (p[i] - tree.ps[i, p_idx])^2
        end
        update_best_dist!(alt, dist, p_idx)

        l += 1
    end

    return alt.best_pidx
end


text:
# The bulk of LLVM's object model consists of values, which comprise a very rich type
# hierarchy.

export Value

abstract type Value end
reftype(::Type{T}) where {T<:Value} = API.LLVMValueRef

identify(::Type{Value}, ref::API.LLVMValueRef) =
    identify(Value, Val{API.LLVMGetValueKind(ref)}())
identify(::Type{Value}, ::Val{K}) where {K} = error("Unknown value kind $K")

@inline function check(::Type{T}, ref::API.LLVMValueRef) where T<:Value
    ref==C_NULL && throw(UndefRefError())
    if Base.JLOptions().debug_level >= 2
        T′ = identify(Value, ref)
        if T != T′
            error("invalid conversion of $T′ value reference to $T")
        end
    end
end

# Construct a concretely typed value object from an abstract value ref
function Value(ref::API.LLVMValueRef)
    ref == C_NULL && throw(UndefRefError())
    T = identify(Value, ref)
    return T(ref)
end

## general APIs

export llvmtype, name, name!, replace_uses!, isconstant, isundef, context

llvmtype(val::Value) = LLVMType(API.LLVMTypeOf(ref(val)))

name(val::Value) = unsafe_string(API.LLVMGetValueName(ref(val)))
name!(val::Value, name::String) = API.LLVMSetValueName(ref(val), name)

function Base.show(io::IO, val::Value)
    output = unsafe_string(API.LLVMPrintValueToString(ref(val)))
    print(io, output)
end

replace_uses!(old::Value, new::Value) =
    API.LLVMReplaceAllUsesWith(ref(old),
                               ref(new))

isconstant(val::Value) = convert(Core.Bool, API.LLVMIsConstant(ref(val)))

isundef(val::Value) = convert(Core.Bool, API.LLVMIsUndef(ref(val)))

context(val::Value) = Context(API.LLVMGetValueContext(ref(val)))

## user values

include("value/user.jl")

## constants

include("value/constant.jl")

## usage

export Use, user, value

@checked struct Use
    ref::API.LLVMUseRef
end
reftype(::Type{Use}) = API.LLVMUseRef

user(use::Use) =  Value(API.LLVMGetUser(     ref(use)))
value(use::Use) = Value(API.LLVMGetUsedValue(ref(use)))

# use iteration

export uses

struct ValueUseSet
    val::Value
end

uses(val::Value) = ValueUseSet(val)

Base.eltype(::ValueUseSet) = Use

function Base.iterate(iter::ValueUseSet, state=API.LLVMGetFirstUse(ref(iter.val)))
    state == C_NULL ? nothing : (Use(state), API.LLVMGetNextUse(state))
end

Base.IteratorSize(::ValueUseSet) = Base.SizeUnknown()


text:
export conversion_dict

# Julia file to contain the conversion dictionary from latin chars to emoji
conversion_dict = Dict{String,String}()

# Punctuation
conversion_dict["\n"] = "🔚\n"
conversion_dict[" "] = "⬛"
conversion_dict["!"] = "❗"
conversion_dict["\""] = "💬"
conversion_dict["#"] = "💷"
conversion_dict["\$"] = "💵"
conversion_dict["%"] = "➿"
conversion_dict["&"] = "🃏"
conversion_dict["'"] = "🗯"
conversion_dict["("] = "↪"
conversion_dict[")"] = "↩"
conversion_dict["*"] = "✖"
conversion_dict["+"] = "➕"
conversion_dict[","] = "🔸"
conversion_dict["-"] = "➖"
conversion_dict["."] = "❌"
conversion_dict["/"] = "➗"
conversion_dict[":"] = "✔"
conversion_dict[";"] = "☑"
conversion_dict["<"] = "↗"
conversion_dict["="] = "💠"
conversion_dict[">"] = "↖"
conversion_dict["?"] = "❓"
conversion_dict["@"] = "🌐"
conversion_dict["_"] = "🏁"

# Numbers
conversion_dict["0"] = "0⃣"
conversion_dict["1"] = "1⃣"
conversion_dict["2"] = "2⃣"
conversion_dict["3"] = "3⃣"
conversion_dict["4"] = "4⃣"
conversion_dict["5"] = "5⃣"
conversion_dict["6"] = "6⃣"
conversion_dict["7"] = "7⃣"
conversion_dict["8"] = "8⃣"
conversion_dict["9"] = "9⃣"

# Lower case
conversion_dict["a"] = "👽"
conversion_dict["b"] = "💤"
conversion_dict["c"] = "🤠"
conversion_dict["d"] = "☠"
conversion_dict["e"] = "😎"
conversion_dict["f"] = "📂"
conversion_dict["g"] = "👻"
conversion_dict["h"] = "🏠"
conversion_dict["i"] = "🏗"
conversion_dict["j"] = "🛩"
conversion_dict["k"] = "🍾"
conversion_dict["l"] = "🐞"
conversion_dict["m"] = "🐭"
conversion_dict["n"] = "🙅"
conversion_dict["o"] = "🎱"
conversion_dict["p"] = "🐧"
conversion_dict["q"] = "😤"
conversion_dict["r"] = "👾"
conversion_dict["s"] = "😭"
conversion_dict["t"] = "👔"
conversion_dict["u"] = "👥"
conversion_dict["v"] = "🐸"
conversion_dict["w"] = "🌕"
conversion_dict["x"] = "😵"
conversion_dict["y"] = "🤔"
conversion_dict["z"] = "😈"

# Upper case
conversion_dict["A"] = "👾"
conversion_dict["B"] = "🦇"
conversion_dict["C"] = "💀"
conversion_dict["D"] = "🤒"
conversion_dict["E"] = "😶"
conversion_dict["F"] = "🍞"
conversion_dict["G"] = "🔱"
conversion_dict["H"] = "🚜"
conversion_dict["I"] = "😼"
conversion_dict["J"] = "🐦"
conversion_dict["K"] = "🎫"
conversion_dict["L"] = "🚨"
conversion_dict["M"] = "🚇"
conversion_dict["N"] = "🎧"
conversion_dict["O"] = "🦉"
conversion_dict["P"] = "🐼"
conversion_dict["Q"] = "🤤"
conversion_dict["R"] = "💥"
conversion_dict["S"] = "🦈"
conversion_dict["T"] = "🌳"
conversion_dict["U"] = "🔥"
conversion_dict["V"] = "🍳"
conversion_dict["W"] = "🕸"
conversion_dict["X"] = "🥅"
conversion_dict["Y"] = "😱"
conversion_dict["Z"] = "🐠"

# Back convert emoji to ASCII
dict_conversion = Dict{String,String}()
for (key, value) in conversion_dict
  # println(key)
  # println(value)
  dict_conversion[value] = key
  dict_conversion[value]
  end
# println(dict_conversion)


text:
push!(LOAD_PATH,"../src/")

using Documenter
using DocumenterLaTeX
include("../src/LPBasic.jl")
#using LPBasic

makedocs(sitename="LPBasic.jl Documentation", authors = "André Mazal Krauss", modules = [LPBasic])


text:
# Autogenerated wrapper script for librealsense_jll for armv7l-linux-gnueabihf-cxx11
export librealsense2

using libusb_jll
JLLWrappers.@generate_wrapper_header("librealsense")
JLLWrappers.@declare_library_product(librealsense2, "librealsense2.so.2.41")
function __init__()
    JLLWrappers.@generate_init_header(libusb_jll)
    JLLWrappers.@init_library_product(
        librealsense2,
        "lib/librealsense2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export rm_multiedges!, rm_degree12, to_simplified_tnet
using LightGraphs: SimpleGraph, add_edge!

function rm_multiedges!(tn::TensorNetwork)
    n = length(tn.tensors)
    for i=1:n
        ti = tn.tensors[i]
        for j=i+1:n
            tj = tn.tensors[j]
            if count(x->x ∈ ti.labels, tj.labels) > 1
                # rm multi edges
                tn.tensors[i], tn.tensors[j] = rm_multiedge(ti, tj)
            end
        end
    end
    return tn
end

function rm_multiedge(ti, tj, li::Tuple, lj::Tuple)
    ti, tj, li, lj = rm_multiedge(ti, tj, collect(li), collect(lj))
    ti, tj, (li...,), (lj...,)
end

function rm_multiedge(t1, t2)
    ti, li, tj, lj = t1.array, t1.labels, t2.array, t2.labels
    common_edges = li ∩ lj
    dimsi = indexin(common_edges, li)
    remsi = setdiff(1:length(li), dimsi)
    dimsj = indexin(common_edges, lj)
    remsj = setdiff(1:length(lj), dimsj)
    orderi = Int[(remsi ∪ dimsi)...]
    orderj = Int[(remsj ∪ dimsj)...]
    # permute
    ti = permutedims(ti, orderi)
    tj = permutedims(tj, orderj)
    li = li[orderi]
    lj = lj[orderj]
    # reshape
    ti = reshape(ti, size(ti)[1:length(remsi)]..., :)
    tj = reshape(tj, size(tj)[1:length(remsj)]..., :)
    li = li[1:length(remsi)+1]
    lj = lj[1:length(remsj)+1]
    return similar(t1, ti, li), similar(t2, tj, lj)
end

function tn2graph(tn::TensorNetwork)
    n = length(tn.tensors)
    g = SimpleGraph(n)
    for i=1:n
        ti = tn.tensors[i]
        for j=i+1:n
            if any(x->x ∈ ti.labels, tn.tensors[j].labels)
                add_edge!(g, i, j)
            end
        end
    end
    return g
end

function rm_degree12(tn::TensorNetwork{T}) where T
    tn = copy(tn)
    n = length(tn)
    mask = ones(Bool, n)
    has_dangling_edges = true
    factor = one(T)
    while has_dangling_edges && length(tn) > 1
        has_dangling_edges = false
        for i=1:n
            mask[i] || continue
            ti = tn.tensors[i]
            if ndims(ti) <= 2
                has_dangling_edges = true
                mask[i] = false
                if ndims(ti) != 0
                    j = findfirst(k -> mask[k] && (ti.labels[1] ∈ tn.tensors[k].labels), 1:n)
                    # absorb i -> j
                    tn.tensors[j] = ti * tn.tensors[j]
                else
                    factor *= Array(tn.tensors[i].array)[]
                end
            end
        end
    end
    return factor, TensorNetwork(tn.tensors[mask])
end


text:
### A Pluto.jl notebook ###
# v0.19.0

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ ba502542-c9ec-4255-a334-047dca9455d6
begin
ROOT = @__DIR__

# Activate the local environment
using Pkg
Pkg.activate(ROOT, io = devnull)

# Load the packages
using Statistics
using Plots
using PlutoUI
using GraphRecipes
using StatsPlots

md"This is a Pluto notebook for showing off some `Plots.jl` goodness and preparing plots for publication."
end

# ╔═╡ 80d8f718-1ffe-4ed0-aefb-2e4cd476997c
@bind params confirm(
    PlutoUI.combine() do Child
        md"""
        Backend: $(Child("backend", Select(
            [:GR, :PGFPlotsX, :PlotlyJS, :PyPlot],
            default = :GR
        )))
         DPI: $(Child("dpi", NumberField(
            120:1:400,
            default = 300,
        )))
         Size factor: $(Child("size_factor", NumberField(
            1:0.1:10,
            default = 5,
        )))\
        Base size: $(Child("size_x", NumberField(
            50:600,
            default = 120,
        )))
        $(Child("size_y", NumberField(
            50:600,
            default = 80,
        )))\

        """
    end
)

# ╔═╡ d356cab0-f755-41f0-9d3b-4f291f822822
begin
# Unpack some the parameters
backend = params.backend
dpi = params.dpi

# Define the path to the plots directory
PLOTS = joinpath(ROOT, "plots", "$(backend)")
mkpath(PLOTS)

# Compute the plot size
size = params.size_factor * [params.size_x, params.size_y]

"Enable the backend"
function enable(backend)
    if backend == :GR
        gr()
    elseif backend == :PGFPlotsX
        pgfplotsx()
    elseif backend == :PlotlyJS
        plotlyjs()
    elseif backend == :PyPlot
        pyplot()
    end
end

"Save the plot"
function save(figure, backend)
    if backend == :PGFPlotsX
        savefig(joinpath(PLOTS, figure * ".pdf"))
    else
        savefig(joinpath(PLOTS, figure * ".png"))
    end
end

# Print the info
md"""
Plot size: $(size[1]), $(size[2]) px
  Final size: $(size[1] * dpi / 100), $(size[2] * dpi / 100) px
"""
end

# ╔═╡ 40983e0a-1340-4b66-bb84-c0da12c6da25
@bind plots confirm(
    MultiCheckBox(
        [
            "lines" => "Lines",
            "histogram" => "Histogram",
            "parametric" => "Parametric",
            "polar" => "Polar plots",
            "markers" => "Markers",
            "colors" => "Colors",
            "global" => "Global",
            "group" => "Group and Subplots",
            "contours" => "Contours",
            "heatmap" => "Heatmap",
            "histogram2d" => "Histogram2D",
            "3D" => "3D",
            "surface" => "3D surface",
            "heart" => "Animation (Heart)",
            "waves" => "Animation (Waves)",
            "graph" => "Graph",
            "AST" => "AST",
            "marginalkde" => "Marginal KDE",
        ],
        orientation=:column,
        select_all=true,
    )
)

# ╔═╡ b0b2cc0f-b7d1-4a41-84e0-4bc99a33aa3d
md"#### Lines"

# ╔═╡ 10c757b6-7719-45e1-8ecc-eb79de3f0e86
begin
if "lines" in plots
enable(backend)

# Plot some fake data
p_lines = plot(Plots.fakedata(50, 5); linewidth=3, legend=:topright, dpi, size)

# Save the plot
save("lines", backend)

p_lines
end
end

# ╔═╡ 9d0c5a15-874a-4cc3-a81b-7026c23b1f3f
md"#### Histogram"

# ╔═╡ 14d54583-ca41-4daf-83e7-4dc26eeac93f
begin
if "histogram" in plots
enable(backend)

# Plot a histogram
p_histogram = histogram(
    randn(1000);
    bins=:scott,
    weights=repeat(1:5, outer = 200),
    dpi,
    size,
)

# Save the plot
save("histogram", backend)

p_histogram
end
end

# ╔═╡ f2b8a98f-8815-4ecf-b726-8de0b5ef52b6
md"#### Parametric plots"

# ╔═╡ 6b5e1ed3-5936-44d2-8acb-ce10ed61a529
begin
if "parametric" in plots
enable(backend)

# Plot a parametric plot
p_parametric = plot(
    sin,
    x -> sin(2x),
    0,
    2π;
    lw=4,
    legend=false,
    fill=(0, :orange),
    dpi,
    size,
)

# Save the plot
save("parametric", backend)

p_parametric
end
end

# ╔═╡ b6f1a674-88cc-4e3d-9bef-b73de6126f1b
md"#### Polar plots"

# ╔═╡ bd7a3e17-054f-4bf5-b4e9-8b41ca038b35
begin
if "polar" in plots
enable(backend)

# Generate data
Θ = range(0, 1.5π, length = 100)
r = abs.(0.1 * randn(100) + sin.(3Θ))

# Plot a polar plot
p_polar = plot(Θ, r; proj=:polar, m=2, dpi, size)

# Save the plot
save("polar", backend)

p_polar
end
end

# ╔═╡ f0565d5a-5ab1-4dfd-9097-4543bb3ed621
md"#### Markers"

# ╔═╡ 39238fb7-c2c3-435c-a26f-b3eaf3f6c9ae
begin
if "markers" in plots
enable(backend)

# Get the markers
markers = filter((m -> m in Plots.supported_markers()), Plots._shape_keys)
markers = permutedims(markers)

# Generate a pseudo-grid
n_markers = length(markers)
x_markers = (range(0, stop = 10, length = n_markers + 2))[2:end - 1]
y_markers = repeat(reshape(reverse(x_markers), 1, :), n_markers, 1)

# Plot the markers
p_markers = scatter(
    x_markers,
    y_markers;
    m=markers,
    markersize=7,
    label=map(string, markers),
    bg=:linen,
    xlim=(0, 10),
    ylim=(0, 10),
    legend=:topright,
    legendfontsize=8,
    dpi,
    size,
)

# Save the plot
save("markers", backend)

p_markers
end
end

# ╔═╡ c5a83092-6494-4f76-8e43-66c0449c8e44
md"#### Colors"

# ╔═╡ dbc6f8a2-1e55-4eb4-b163-d5050b135ce5
begin
if "colors" in plots
enable(backend)

# Plot four random lines
p_colors = plot(
    0:10:100,
    rand(11, 4);
    w=3,
    palette=cgrad(:inferno),
    fill=0,
    α=0.6,
    legend=false,
    dpi,
    size,
)

# Plot a bunch of random points
y_colors = rand(100)
scatter!(
    y_colors,
    zcolor=abs.(y_colors .- 0.5),
    m=(:heat, 0.8, Plots.stroke(1, :green)),
    ms=10 * abs.(y_colors .- 0.5) .+ 4,
    colorbar=:best,
)

# Save the plot
save("colors", backend)

p_colors
end
end

# ╔═╡ ab076971-1f57-4193-a347-926d4004c78b
md"#### Global"

# ╔═╡ 4ad6f925-9512-4724-9ee1-d0a02b06c430
begin
if "global" in plots
enable(backend)

# Generate data
y_global = rand(20, 3)

# Plot the data
p_global = plot(
    y_global;
    xaxis=("XLABEL", (-5, 30), 0:2:20, :flip),
    background_color=RGB(0.2, 0.2, 0.2),
    legend=false,
    dpi,
    size,
)

# Mutate the plot
hline!(
    mean(y_global, dims = 1) + rand(1, 3);
    line=(4, :dash, 0.6, [:lightgreen :green :darkgreen]),
)
vline!([5, 10])
title!("TITLE")
yaxis!("YLABEL", :log10, minorgrid=true)

# Save the plot
save("global", backend)

p_global
end
end

# ╔═╡ 813396bf-f8a9-4f5c-8498-cdd770183273
md"#### Groups and Subplots"

# ╔═╡ fd5bd4ff-812b-47e4-b0b0-257b6636e5e9
begin
if "group" in plots
enable(backend)

# Generate random group associations
group = rand(map(i -> "group $(i)", 1:4), 100)

# Group some random data
p_group = plot(
    rand(100);
    layout=@layout([a b; c]),
    group,
    linetype=[:bar :scatter :steppre],
    linecolor=:match,
    legend=:topright,
    dpi,
    size,
)

# Save the plot
save("group", backend)

p_group
end
end

# ╔═╡ e9f22891-bd49-489d-8799-c48c20eae540
md"#### Contours"

# ╔═╡ 4036bb6c-9ced-4811-943b-0cccdf121a8b
begin
if "contours" in plots
enable(backend)

# Prepare data
x_contours = 1:0.5:20
y_contours = 1:0.5:10
f(x, y) = (3x + y ^ 2) * abs(sin(x) + cos(y))
X = repeat(reshape(x_contours, 1, :), length(y_contours), 1)
Y = repeat(y_contours, 1, length(x_contours))
Z = map(f, X, Y)
p1 = contour(x_contours, y_contours, f; fill=true)
p2 = contour(x_contours, y_contours, Z)

# Plot the contours
p_contours = plot(p1, p2; layout=(2, 1), dpi, size)

# Save the plot
save("contours", backend)

p_contours
end
end

# ╔═╡ 9e8e2f58-894d-46ba-9390-dbd0edfae4b0
md"#### Heatmap"

# ╔═╡ 5e32d1ce-2534-4a83-897a-be6369c50d1f
begin
if "heatmap" in plots
enable(backend)

# Prepare data
xs_histogram = [string("x", i) for i = 1:10]
ys_histogram = [string("y", i) for i = 1:4]
z_histogram = float((1:4) * reshape(1:10, 1, :))

# Plot a heatmap
p_heatmap = heatmap(
    xs_histogram,
    ys_histogram,
    z_histogram;
    aspect_ratio=1,
    dpi,
    size,
)

# Save the plot
save("heatmap", backend)

p_heatmap
end
end

# ╔═╡ 321184c0-415c-4c8f-b94c-3492d6f8782c
md"#### Histogram2D"

# ╔═╡ 42a45eea-cd2d-4953-9266-6d30e91efd82
begin
if "histogram2d" in plots
enable(backend)

# Compute the data
n_histogram2d = 10000
x_histogram2d = exp.(0.1 * randn(n_histogram2d) .+ randn(n_histogram2d) .* im)

# Plot the 2D-histogram
p_histogram2d = histogram2d(
    x_histogram2d;
    nbins=(20, 40),
    show_empty_bins=true,
    normed=true,
    aspect_ratio=1,
    dpi,
    size,
)

# Save the plot
save("histogram2d", backend)

p_histogram2d
end
end

# ╔═╡ d89da106-72d1-4579-b9f6-f044e1eee554
md"#### 3D"

# ╔═╡ 8c7bc68d-ffee-4ec8-b8bb-3c081ac193f3
begin
if "3D" in plots
enable(backend)

# Generate the data for the spiral
n_3D = 100
ts = range(0, stop = 8π, length = n_3D)
x_3D = ts .* map(cos, ts)
y_3D = (0.1ts) .* map(sin, ts)
z_3D = 1:n_3D

# Plot the spiral
p_3D = plot(
    x_3D,
    y_3D,
    z_3D;
    camera=(30, 30),
    zcolor=reverse(z_3D),
    marker=(10, 0.8, :blues, Plots.stroke(0)),
    legend=false,
    cbar=true,
    w=5,
    dpi,
    size,
)

# Plot the line
plot!(zeros(n_3D), zeros(n_3D), 1:n_3D, w = 10)

# Save the plot
save("3D", backend)

p_3D
end
end

# ╔═╡ 22c30b75-1106-426c-9a45-93b4f58a9c8d
md"#### 3D surface"

# ╔═╡ 2d129cdb-108a-4b7d-995f-33b5d5eb62d6
begin
if "surface" in plots
enable(backend)

# Generate data
x_surface = 0:π/20:4π
y_surface = 0:π/20:4π
z_surface = 5 .* sin.(x_surface) * cos.(y_surface)'

# Plot the surface
p_surface = plot(
    x_surface,
    y_surface,
    z_surface;
    st=:surface,
    dpi,
    size,
)

# Save the plot
save("surface", backend)

p_surface
end
end

# ╔═╡ e3bee7b1-9dcd-4568-b735-2051e9630243
md"#### Animation (Heart)"

# ╔═╡ 713db5a1-6c2f-43be-b041-fe74b95f3b0c
begin
if "heart" in plots
gr()

# Create a recipe for a custom type
@userplot CirclePlot
@recipe function f(cp::CirclePlot)
    x, y, i = cp.args
    n = length(x)
    inds = circshift(1:n, 1 - i)
    linewidth --> range(0, 10, length = n)
    seriesalpha --> range(0, 1, length = n)
    aspect_ratio --> 1
    label --> false
    x[inds], y[inds]
end

# Generate data for the heart
n_animation = 400
t_animation = range(0, 2π, length = n_animation)
x_animation = 16sin.(t_animation).^3
y_animation = 13cos.(t_animation) .- 5cos.(2t_animation) .-
              2cos.(3t_animation) .- cos.(4t_animation)

# Create an animation
animation = @animate for i ∈ 1:n_animation
    circleplot(
        x_animation,
        y_animation,
        i;
        line_z=1:n_animation,
        cbar=false,
        c=:reds,
        framestyle=:none,
        dpi,
        size,
    )
end when mod1(i, 10) == 5

# Create a GIF
gif(animation, joinpath(ROOT, "plots", "heart.gif"), fps=10, show_msg=false)
end
end

# ╔═╡ bccb93e1-5d66-4fc5-9da8-c6386b065a98
md"#### Animation (Waves)"

# ╔═╡ 2c556d8e-3c9f-4cba-b60a-45bc2dbe24c1
begin
if "waves" in plots
gr()

# Prepare data
x = y = range(-5, 5, length = 40)
zs = zeros(0, 40)
n = 100

# Create an animation
anim_waves = @animate for i in range(0, 2π; length = n)
    f(x, y) = sin(x + 10sin(i)) + cos(y)

    # Create a plot with 3 subplots and a custom layout
    l = @layout [a{0.7w} b; c{0.2h}]
    p = plot(x, y, f; st=[:surface, :contourf], layout=l, legend=false)

    # Induce a slight oscillating camera angle sweep, in degrees (azimuth, altitude)
    plot!(p[1], camera = (10 * (1 + cos(i)), 40))

    # Add a tracking line
    fixed_x = zeros(40)
    z = map(f, fixed_x, y)
    plot!(p[1], fixed_x, y, z, line = (:black, 5, 0.2))
    vline!(p[2], [0], line = (:black, 5))

    # Add to and show the tracked values over time
    global zs = vcat(zs, z')
    plot!(p[3], zs, α=0.2, palette=cgrad(:blues).colors, legend=false)
end

# Create a GIF
gif(anim_waves, joinpath(ROOT, "plots", "waves.gif"), fps=15, show_msg=false)
end
end

# ╔═╡ 1b04c4d4-8fc9-4cb2-afea-99c0d9f56c73
md"#### Graph `(GraphRecipes)`"

# ╔═╡ 07e936fd-d277-4e1d-be2c-51070918d226
begin
if "graph" in plots
enable(backend)

# Generate data
n_graph= 15
A = Float64[rand() < 0.5 ? 0 : rand() for i=1:n_graph, j=1:n_graph]
for i=1:n_graph
    A[i, 1:i-1] = A[1:i-1, i]
    A[i, i] = 0
end

# Plot the graph
p_graph = graphplot(
    A;
    markersize=0.2,
    node_weights=1:n_graph,
    markercolor=range(colorant"yellow", stop=colorant"red", length=n_graph),
    names=1:n_graph,
    fontsize=10,
    linecolor=:darkgrey,
    dpi,
    size,
)

# Save the plot
save("graph", backend)

p_graph
end
end

# ╔═╡ 24b90639-fdb8-4bc4-b316-f75c375e2c96
md"#### AST `(GraphRecipes)`"

# ╔═╡ f572e4c4-5fab-44c0-b8ff-0b2689aaea20
begin
if "AST" in plots
enable(backend)

# Prepare an expression
code = :(
function mysum(list)
    out = 0
    for value in list
        out += value
    end
    out
end
)

# Plot the AST
p_ast = plot(
    code;
    fontsize=12,
    shorten=0.01,
    axis_buffer=0.15,
    nodeshape=:rect,
    dpi,
    size,
)

# Save the plot
save("ast", backend)

p_ast
end
end

# ╔═╡ 900c5ddc-40b6-4415-8bf2-929c7490bd74
md"#### Marginal KDE `(StatsPlots)`"

# ╔═╡ 7189bccd-f128-41cc-b7c9-6b81796df1b5
begin
if "marginalkde" in plots
enable(backend)

# Generate some random data
x_marginalkde = randn(1024)
y_marginalkde = randn(1024)

# Plot a marginal kernel density estimation
p_marginalkde = marginalkde(x_marginalkde, x_marginalkde + y_marginalkde; dpi, size)

# Save the plot
save("marginalkde", backend)

p_marginalkde
end
end

# ╔═╡ Cell order:
# ╟─ba502542-c9ec-4255-a334-047dca9455d6
# ╟─80d8f718-1ffe-4ed0-aefb-2e4cd476997c
# ╟─d356cab0-f755-41f0-9d3b-4f291f822822
# ╟─40983e0a-1340-4b66-bb84-c0da12c6da25
# ╟─b0b2cc0f-b7d1-4a41-84e0-4bc99a33aa3d
# ╠═10c757b6-7719-45e1-8ecc-eb79de3f0e86
# ╟─9d0c5a15-874a-4cc3-a81b-7026c23b1f3f
# ╠═14d54583-ca41-4daf-83e7-4dc26eeac93f
# ╟─f2b8a98f-8815-4ecf-b726-8de0b5ef52b6
# ╟─6b5e1ed3-5936-44d2-8acb-ce10ed61a529
# ╟─b6f1a674-88cc-4e3d-9bef-b73de6126f1b
# ╠═bd7a3e17-054f-4bf5-b4e9-8b41ca038b35
# ╟─f0565d5a-5ab1-4dfd-9097-4543bb3ed621
# ╟─39238fb7-c2c3-435c-a26f-b3eaf3f6c9ae
# ╟─c5a83092-6494-4f76-8e43-66c0449c8e44
# ╟─dbc6f8a2-1e55-4eb4-b163-d5050b135ce5
# ╟─ab076971-1f57-4193-a347-926d4004c78b
# ╟─4ad6f925-9512-4724-9ee1-d0a02b06c430
# ╟─813396bf-f8a9-4f5c-8498-cdd770183273
# ╟─fd5bd4ff-812b-47e4-b0b0-257b6636e5e9
# ╟─e9f22891-bd49-489d-8799-c48c20eae540
# ╟─4036bb6c-9ced-4811-943b-0cccdf121a8b
# ╟─9e8e2f58-894d-46ba-9390-dbd0edfae4b0
# ╟─5e32d1ce-2534-4a83-897a-be6369c50d1f
# ╟─321184c0-415c-4c8f-b94c-3492d6f8782c
# ╟─42a45eea-cd2d-4953-9266-6d30e91efd82
# ╟─d89da106-72d1-4579-b9f6-f044e1eee554
# ╟─8c7bc68d-ffee-4ec8-b8bb-3c081ac193f3
# ╟─22c30b75-1106-426c-9a45-93b4f58a9c8d
# ╟─2d129cdb-108a-4b7d-995f-33b5d5eb62d6
# ╟─e3bee7b1-9dcd-4568-b735-2051e9630243
# ╟─713db5a1-6c2f-43be-b041-fe74b95f3b0c
# ╟─bccb93e1-5d66-4fc5-9da8-c6386b065a98
# ╟─2c556d8e-3c9f-4cba-b60a-45bc2dbe24c1
# ╟─1b04c4d4-8fc9-4cb2-afea-99c0d9f56c73
# ╟─07e936fd-d277-4e1d-be2c-51070918d226
# ╟─24b90639-fdb8-4bc4-b316-f75c375e2c96
# ╟─f572e4c4-5fab-44c0-b8ff-0b2689aaea20
# ╟─900c5ddc-40b6-4415-8bf2-929c7490bd74
# ╟─7189bccd-f128-41cc-b7c9-6b81796df1b5