text:
using Contour: Contour
using UUIDs
using Latexify
Base.@kwdef mutable struct PGFPlotsXPlot
    is_created::Bool = false
    was_shown::Bool = false
    the_plot::PGFPlotsX.TikzDocument = PGFPlotsX.TikzDocument()
    function PGFPlotsXPlot(is_created, was_shown, the_plot)
        pgfx_plot = new(is_created, was_shown, the_plot)
        # tikz libraries
        PGFPlotsX.push_preamble!(
            pgfx_plot.the_plot,
            "\\usetikzlibrary{arrows.meta}",
        )
        PGFPlotsX.push_preamble!(
            pgfx_plot.the_plot,
            "\\usetikzlibrary{backgrounds}",
        )
        # pgfplots libraries
        PGFPlotsX.push_preamble!(
            pgfx_plot.the_plot,
            "\\usepgfplotslibrary{patchplots}",
        )
        PGFPlotsX.push_preamble!(
            pgfx_plot.the_plot,
            "\\usepgfplotslibrary{fillbetween}",
        )
        # compatibility fixes
        # add background layer to standard layers
        PGFPlotsX.push_preamble!(
            pgfx_plot.the_plot,
            raw"""
            \pgfplotsset{%
            layers/standard/.define layer set={%
                background,axis background,axis grid,axis ticks,axis lines,axis tick labels,pre main,main,axis descriptions,axis foreground%
            }{grid style= {/pgfplots/on layer=axis grid},%
                tick style= {/pgfplots/on layer=axis ticks},%
                axis line style= {/pgfplots/on layer=axis lines},%
                label style= {/pgfplots/on layer=axis descriptions},%
                legend style= {/pgfplots/on layer=axis descriptions},%
                title style= {/pgfplots/on layer=axis descriptions},%
                colorbar style= {/pgfplots/on layer=axis descriptions},%
                ticklabel style= {/pgfplots/on layer=axis tick labels},%
                axis background@ style={/pgfplots/on layer=axis background},%
                3d box foreground style={/pgfplots/on layer=axis foreground},%
                },
            }
            """,
        )
        pgfx_plot
    end
end

## end user utility functions
function pgfx_axes(pgfx_plot::PGFPlotsXPlot)
    return pgfx_plot.the_plot.elements[1].elements
end

pgfx_preamble() = pgfx_preamble(current())
function pgfx_preamble(pgfx_plot::Plot{PGFPlotsXBackend})
    old_flag = pgfx_plot.attr[:tex_output_standalone]
    pgfx_plot.attr[:tex_output_standalone] = true
    fulltext = String(repr("application/x-tex", pgfx_plot))
    preamble = fulltext[1:(first(findfirst("\\begin{document}", fulltext)) - 1)]
    pgfx_plot.attr[:tex_output_standalone] = old_flag
    preamble
end
##

function surface_to_vecs(x::AVec, y::AVec, s::Union{AMat,Surface})
    a = Array(s)
    xn = Vector{eltype(x)}(undef, length(a))
    yn = Vector{eltype(y)}(undef, length(a))
    zn = Vector{eltype(s)}(undef, length(a))
    for (n, (i, j)) in enumerate(Tuple.(CartesianIndices(a)))
        if length(x) == size(s)[1]
            i, j = j, i
        end
        xn[n] = x[j]
        yn[n] = y[i]
        zn[n] = a[i, j]
    end
    return xn, yn, zn
end
surface_to_vecs(x::AVec, y::AVec, z::AVec) = x, y, z

function Base.push!(pgfx_plot::PGFPlotsXPlot, item)
    push!(pgfx_plot.the_plot, item)
end

function pgfx_split_extra_opts(extra)
    return (get(extra, :add, nothing), filter( x-> first(x) != :add, extra))
end
function (pgfx_plot::PGFPlotsXPlot)(plt::Plot{PGFPlotsXBackend})
    if !pgfx_plot.is_created || pgfx_plot.was_shown
        pgfx_sanitize_plot!(plt)
        # extract extra kwargs
        extra_plot, extra_plot_opt = pgfx_split_extra_opts(plt[:extra_plot_kwargs])
        the_plot = PGFPlotsX.TikzPicture(PGFPlotsX.Options(extra_plot_opt...))
        if extra_plot !== nothing
            extra_plot = wraptuple(extra_plot)
            push!(the_plot, extra_plot...)
        end
        bgc = plt.attr[:background_color_outside] == :match ?
            plt.attr[:background_color] : plt.attr[:background_color_outside]
        if bgc isa Colors.Colorant
            cstr = plot_color(bgc)
            a = alpha(cstr)
            push!(
                the_plot.options,
                "/tikz/background rectangle/.style" => PGFPlotsX.Options(
                    # "draw" => "black",
                    "fill" => cstr,
                    "draw opacity" => a,
                ),
                "show background rectangle" => nothing,
            )
        end

        for sp in plt.subplots
            bb2 = bbox(sp)
            dx, dy = bb2.x0
            sp_width = width(bb2)
            sp_height = height(bb2)
            if sp[:subplot_index] == plt[:plot_titleindex]
                x = dx + sp_width / 2 - 10mm # FIXME: get rid of magic constant
                y = dy + sp_height / 2
                pgfx_add_annotation!(the_plot, x, y, PlotText(plt[:plot_title], plottitlefont(plt)), pgfx_thickness_scaling(plt);
                    cs = "",
                    options = PGFPlotsX.Options("anchor" => "center")
                )
                continue
            end

            sp_width = width(bb2)
            sp_height = height(bb2)
            lpad = leftpad(sp) + sp[:left_margin]
            rpad = rightpad(sp) + sp[:right_margin]
            tpad = toppad(sp) + sp[:top_margin]
            bpad = bottompad(sp) + sp[:bottom_margin]
            dx += lpad
            dy += tpad
            axis_height = sp_height - (tpad + bpad)
            axis_width = sp_width - (rpad + lpad)
            title_cstr = plot_color(sp[:titlefontcolor])
            title_a = alpha(title_cstr)
            title_loc = sp[:titlelocation]
            bgc_inside = plot_color(sp[:background_color_inside])
            bgc_inside_a = alpha(bgc_inside)
            axis_opt = PGFPlotsX.Options(
                "point meta max" => get_clims(sp)[2],
                "point meta min" => get_clims(sp)[1],
                "title" => sp[:title],
                "title style" => PGFPlotsX.Options(
                        pgfx_get_title_pos(title_loc)...,
                        "font" => pgfx_font(
                            sp[:titlefontsize],
                            pgfx_thickness_scaling(sp),
                    ),
                    "color" => title_cstr,
                    "draw opacity" => title_a,
                    "rotate" => sp[:titlefontrotation],
                ),
                "legend style" => pgfx_get_legend_style(sp),
                "axis background/.style" => PGFPlotsX.Options(
                    "fill" => bgc_inside,
                    "opacity" => bgc_inside_a,
                ),
                # These are for layouting
                "anchor" => "north west",
                "xshift" => string(dx),
                "yshift" => string(-dy),
            )
            sp_width > 0 * mm ? push!(axis_opt, "width" => string(axis_width)) :
            nothing
            sp_height > 0 * mm ? push!(axis_opt, "height" => string(axis_height)) :
            nothing
            for letter in (:x, :y, :z)
                if letter != :z || RecipesPipeline.is3d(sp)
                    pgfx_axis!(axis_opt, sp, letter)
                end
            end
            # Search series for any gradient. In case one series uses a gradient set
            # the colorbar and colomap.
            # The reasoning behind doing this on the axis level is that pgfplots
            # colorbar seems to only works on axis level and needs the proper colormap for
            # correctly displaying it.
            # It's also possible to assign the colormap to the series itself but
            # then the colormap needs to be added twice, once for the axis and once for the
            # series.
            # As it is likely that all series within the same axis use the same
            # colormap this should not cause any problem.
            for series in series_list(sp)
                if hascolorbar(series)
                    cg = get_colorgradient(series)
                    cm = pgfx_colormap(get_colorgradient(series))
                    PGFPlotsX.push_preamble!(
                        pgfx_plot.the_plot,
                        """\\pgfplotsset{
                        colormap={plots$(sp.attr[:subplot_index])}{$cm},
                        }""",
                    )
                    push!(axis_opt, "colormap name" => "plots$(sp.attr[:subplot_index])")
                    if cg isa PlotUtils.CategoricalColorGradient
                        push!(
                            axis_opt,
                            "colormap access" => "piecewise const",
                            "colorbar sampled" => nothing,
                        )
                    end
                    # goto is needed to break out of col and series for
                    @goto colorbar_end
                end
            end
            @label colorbar_end

            if hascolorbar(sp)
                cticks = get_colorbar_ticks(sp)[2]
                colorbar_style = PGFPlotsX.Options(
                    "title" => sp[:colorbar_title],
                )
                if sp[:colorbar] === :top
                    push!(colorbar_style,
                        "at" => string((0.5, 1.05)),
                        "anchor" => "south",
                        "xtick" => string("{", join(cticks, ","), "}"),
                        "xticklabel pos" => "upper",
                        "xticklabel style" => pgfx_get_colorbar_ticklabel_style(sp),
                    )
                else
                    push!(colorbar_style,
                        "ytick" => string("{", join(cticks, ","), "}"),
                        "yticklabel style" => pgfx_get_colorbar_ticklabel_style(sp),
                    )
                end
                push!(
                    axis_opt,
                    string("colorbar", pgfx_get_colorbar_pos(sp[:colorbar])) => nothing,
                    "colorbar style" => colorbar_style,
                )
            else
                push!(axis_opt, "colorbar" => "false")
            end
            if RecipesPipeline.is3d(sp)
                azim, elev = sp[:camera]
                push!(axis_opt, "view" => (azim, elev))
            end
            axisf = if sp[:projection] == :polar
                # push!(axis_opt, "xmin" => 90)
                # push!(axis_opt, "xmax" => 450)
                PGFPlotsX.PolarAxis
            else
                PGFPlotsX.Axis
            end
            extra_sp, extra_sp_opt = pgfx_split_extra_opts(sp[:extra_kwargs])
            axis = axisf(merge(axis_opt, PGFPlotsX.Options(extra_sp_opt...)))
            if extra_sp !== nothing
                extra_sp = wraptuple(extra_sp)
                push!(axis, extra_sp...)
            end
            if sp[:legendtitle] !== nothing
                push!(axis, PGFPlotsX.Options("\\addlegendimage{empty legend}" => nothing))
                push!(
                    axis,
                    PGFPlotsX.LegendEntry(
                        string("\\hspace{-.6cm}{\\textbf{", sp[:legendtitle], "}}"),
                        false,
                    ),
                )
            end
            for (series_index, series) in enumerate(series_list(sp))
                # give each series a uuid for fillbetween
                series_id = uuid4()
                _pgfplotsx_series_ids[Symbol("$series_index")] = series_id
                opt = series.plotattributes
                st = series[:seriestype]
                sf = series[:fillrange]
                series_opt = PGFPlotsX.Options(
                    "color" => single_color(opt[:linecolor]),
                    "name path" => string(series_id),
                )
                extra_series, extra_series_opt = pgfx_split_extra_opts(series[:extra_kwargs])
                series_opt = merge(series_opt, PGFPlotsX.Options(extra_series_opt...))
                if RecipesPipeline.is3d(series) || st in (:heatmap, :contour) || (st == :quiver && opt[:z] !== nothing)
                    series_func = PGFPlotsX.Plot3
                else
                    series_func = PGFPlotsX.Plot
                end
                if sf !== nothing &&
                   !isfilledcontour(series) && series[:ribbon] === nothing
                    push!(series_opt, "area legend" => nothing)
                end
                pgfx_add_series!(Val(st), axis, series_opt, series, series_func, opt)
                if extra_series !== nothing
                    extra_series = wraptuple(extra_series)
                    push!(axis.contents[end], extra_series...)
                end
                # add series annotations
                anns = series[:series_annotations]
                for (xi, yi, str, fnt) in EachAnn(anns, series[:x], series[:y])
                    pgfx_add_annotation!(
                        axis,
                        xi,
                        yi,
                        PlotText(str, fnt),
                        pgfx_thickness_scaling(series),
                    )
                end
            end # for series
            # add subplot annotations
            for ann in sp[:annotations]
                pgfx_add_annotation!(
                    axis,
                    locate_annotation(sp, ann...)...,
                    pgfx_thickness_scaling(sp),
                )
            end
            push!(the_plot, axis)
            if length(plt.o.the_plot.elements) > 0
                plt.o.the_plot.elements[1] = the_plot
            else
                push!(plt.o, the_plot)
            end
        end # for subplots
        pgfx_plot.is_created = true
        pgfx_plot.was_shown = false
    end # if
    return pgfx_plot
end
## seriestype specifics
function pgfx_add_series!(axis, series_opt, series, series_func, opt)
    args = pgfx_series_arguments(series, opt)
    series_plot = series_func(series_opt, PGFPlotsX.Table(args...))
    push!(axis, series_plot)
    pgfx_add_legend!(axis, series, opt)
end

function pgfx_add_series!(::Val{:path}, axis, series_opt, series, series_func, opt)
    # treat segments
    segments = collect(series_segments(series, series[:seriestype]; check=true))
    sf = opt[:fillrange]
    for (k, segment) in enumerate(segments)
        i, rng = segment.attr_index, segment.range
        segment_opt = PGFPlotsX.Options()
        segment_opt = merge(segment_opt, pgfx_linestyle(opt, i))
        if opt[:markershape] != :none
            marker = _cycle(opt[:markershape], i)
            if marker isa Shape
                x = marker.x
                y = marker.y
                scale_factor = 0.00125
                mark_size = opt[:markersize] * scale_factor
                path = join(
                    [
                        "($(x[i] * mark_size), $(y[i] * mark_size))"
                        for i in eachindex(x)
                    ],
                    " -- ",
                )
                c = get_markercolor(series, i)
                a = get_markeralpha(series, i)
                PGFPlotsX.push_preamble!(
                    series[:plot_object].o.the_plot,
                    """
                    \\pgfdeclareplotmark{PlotsShape$(series[:series_plotindex])}{
                    \\filldraw
                    $path;
                    }
                    """,
                )
            end
            segment_opt = merge(segment_opt, pgfx_marker(opt, i))
        end
        # add fillrange
        if sf !== nothing &&
           !isfilledcontour(series)
           if sf isa Number || sf isa AVec
                pgfx_fillrange_series!(
                    axis,
                    series,
                    series_func,
                    i,
                    _cycle(sf, rng),
                    rng,
                )
            elseif sf isa Tuple && series[:ribbon] !== nothing
                for sfi in sf
                pgfx_fillrange_series!(
                    axis,
                    series,
                    series_func,
                    i,
                    _cycle(sfi, rng),
                    rng,
                )
                end
            end
            if k == 1 &&
               series[:subplot][:legend] != :none && pgfx_should_add_to_legend(series)
                pgfx_filllegend!(series_opt, opt)
            end
        end
        # handle arrows
        arrow = opt[:arrow]
        if arrow isa Arrow
            arrow_opt = merge(
                segment_opt,
                PGFPlotsX.Options("quiver" => PGFPlotsX.Options(
                    "u" => "\\thisrow{u}",
                    "v" => "\\thisrow{v}",
                    pgfx_arrow(arrow, :head) => nothing,
                    )
                )
            )
            if arrow.side == :head
                x_arrow = opt[:x][rng][end-1:end]
                y_arrow = opt[:y][rng][end-1:end]
                x_path  = opt[:x][rng][1:end-1]
                y_path  = opt[:y][rng][1:end-1]
            elseif arrow.side == :tail
                x_arrow = opt[:x][rng][2:-1:1]
                y_arrow = opt[:y][rng][2:-1:1]
                x_path  = opt[:x][rng][2:end]
                y_path  = opt[:y][rng][2:end]
            elseif arrow.side == :both
                x_arrow = opt[:x][rng][[2,1,end-1,end]]
                y_arrow = opt[:y][rng][[2,1,end-1,end]]
                x_path  = opt[:x][rng][2:end-1]
                y_path  = opt[:y][rng][2:end-1]
            end
            coordinates = PGFPlotsX.Table([
                    :x => x_arrow[1:2:end-1],
                    :y => y_arrow[1:2:end-1],
                    :u => [x_arrow[i] - x_arrow[i-1] for i in 2:2:lastindex(x_arrow)],
                    :v => [y_arrow[i] - y_arrow[i-1] for i in 2:2:lastindex(y_arrow)],
                ])
            arrow_plot =
                series_func(merge(series_opt, arrow_opt), coordinates)
            push!(axis, arrow_plot)
            coordinates = PGFPlotsX.Table(x_path, y_path)
            segment_plot =
                series_func(merge(series_opt, segment_opt), coordinates)
            push!(axis, segment_plot)
        else
            coordinates = PGFPlotsX.Table(pgfx_series_arguments(series, opt, rng)...)
            segment_plot =
                series_func(merge(series_opt, segment_opt), coordinates)
            push!(axis, segment_plot)
        end
        # fill between functions
        if sf isa Tuple && series[:ribbon] === nothing
            sf1, sf2 = sf
            @assert sf1 == series_index "First index of the tuple has to match the current series index."
            push!(
                axis,
                series_func(
                    merge(
                        pgfx_fillstyle(opt, series_index),
                        PGFPlotsX.Options("forget plot" => nothing),
                    ),
                    "fill between [of=$series_id and $(_pgfplotsx_series_ids[Symbol(string(sf2))])]",
                ),
            )
        end
        pgfx_add_legend!(axis, series, opt, k)
    end # for segments
    # get that last marker
    if !isnothing(opt[:y]) && !any(isnan, opt[:y]) && opt[:markershape] isa AVec
        additional_plot = PGFPlotsX.PlotInc(pgfx_marker(opt, length(segments) + 1), PGFPlotsX.Coordinates(tuple((last(opt[:x]), last(opt[:y])))))
        push!(axis, additional_plot)
    end
end

function pgfx_add_series!(::Val{:scatter}, axis, series_opt, series, series_func, opt)
    push!(series_opt, "only marks" => nothing)
    pgfx_add_series!(Val(:path), axis, series_opt, series, series_func, opt)
end

function pgfx_add_series!(::Val{:straightline}, axis, series_opt, series, series_func, opt)
    pgfx_add_series!(Val(:path), axis, series_opt, series, series_func, opt)
end

function pgfx_add_series!(::Val{:path3d}, axis, series_opt, series, series_func, opt)
    pgfx_add_series!(Val(:path), axis, series_opt, series, series_func, opt)
end

function pgfx_add_series!(::Val{:scatter3d}, axis, series_opt, series, series_func, opt)
    push!(series_opt, "only marks" => nothing)
    pgfx_add_series!(Val(:path), axis, series_opt, series, series_func, opt)
end

function pgfx_add_series!(::Val{:surface}, axis, series_opt, series, series_func, opt)
    alpha = get_fillalpha(series)
    push!(
        series_opt,
        "surf" => nothing,
        "mesh/rows" => length(unique(opt[:x])), # unique if its all vectors
        "mesh/cols" => length(unique(opt[:y])),
        "z buffer" => "sort",
        "opacity" => alpha === nothing ? 1.0 : alpha,
    )
    pgfx_add_series!(axis, series_opt, series, series_func, opt)
end

function pgfx_add_series!(::Val{:wireframe}, axis, series_opt, series, series_func, opt)
    push!(series_opt, "mesh" => nothing,
                      "mesh/rows" => length(opt[:x])
         )
    pgfx_add_series!(axis, series_opt, series, series_func, opt)
end

function pgfx_add_series!(::Val{:heatmap}, axis, series_opt, series, series_func, opt)
    push!(axis.options, "view" => "{0}{90}")
    push!(
        series_opt,
        "matrix plot*" => nothing,
        "mesh/rows" => length(opt[:x]),
        "mesh/cols" => length(opt[:y]),
        "point meta" => "\\thisrow{meta}",
    )
    args = pgfx_series_arguments(series, opt)
    meta = [any(!isfinite, r) ? NaN : r[3] for r in zip(args...)]
    for arg in args
        arg[(!isfinite).(arg)] .= 0
    end
    t = PGFPlotsX.Table(["x" => args[1], "y" => args[2], "z" => args[3], "meta" => meta])
    series_plot = series_func(series_opt, t)
    push!(axis, series_plot)
    pgfx_add_legend!(axis, series, opt)
end

function pgfx_add_series!(::Val{:mesh3d}, axis, series_opt, series, series_func, opt)
    ptable = join([ string(i, " ", j, " ", k, "\\\\") for (i, j, k) in zip(opt[:connections]...) ], "\n        ")
    push!(
        series_opt,
        "patch" => nothing,
        "table/row sep" => "\\\\",
        "patch table" => ptable
    )
    pgfx_add_series!(axis, series_opt, series, series_func, opt)
end

function pgfx_add_series!(::Val{:contour}, axis, series_opt, series, series_func, opt)
    push!(axis.options, "view" => "{0}{90}")
    if isfilledcontour(series)
        pgfx_add_series!(Val(:filledcontour), axis, series_opt, series, series_func, opt)
        return nothing
    end
    pgfx_add_series!(Val(:contour3d), axis, series_opt, series, series_func, opt)
    return nothing
end

function pgfx_add_series!(::Val{:filledcontour}, axis, series_opt, series, series_func, opt)
    push!(
        series_opt,
        "contour filled" => PGFPlotsX.Options(), # labels not supported
        "patch type" => "bilinear",
        "shader" => "flat",
    )
    if opt[:levels] isa Number
        push!(series_opt["contour filled"], "number" => opt[:levels])
    elseif opt[:levels] isa AVec
        push!(series_opt["contour filled"], "levels" => opt[:levels])
    end
    pgfx_add_series!(axis, series_opt, series, series_func, opt)
end

function pgfx_add_series!(::Val{:contour3d}, axis, series_opt, series, series_func, opt)
    push!(
        series_opt,
        "contour prepared" => PGFPlotsX.Options("labels" => opt[:contour_labels]),
    )
    series_opt = merge( series_opt, pgfx_linestyle(opt) )
    args = pgfx_series_arguments(series, opt)
    series_plot = series_func(series_opt, PGFPlotsX.Table(Contour.contours(args..., opt[:levels])))
    push!(axis, series_plot)
    pgfx_add_legend!(axis, series, opt)
end

function pgfx_add_series!(::Val{:quiver}, axis, series_opt, series, series_func, opt)
    if opt[:quiver] !== nothing
        push!(
            series_opt,
            "quiver" => PGFPlotsX.Options(
                "u" => "\\thisrow{u}",
                "v" => "\\thisrow{v}",
                pgfx_arrow(opt[:arrow]) => nothing,
            ),
        )
        x = opt[:x]
        y = opt[:y]
        z = opt[:z]
        if z !== nothing
            push!(series_opt["quiver"], "w" => "\\thisrow{w}")
            pgfx_axis!(axis.options, series[:subplot], :z)
            table = PGFPlotsX.Table([
                :x => x,
                :y => y,
                :z => z,
                :u => opt[:quiver][1],
                :v => opt[:quiver][2],
                :w => opt[:quiver][3],
            ])
        else
            table = PGFPlotsX.Table([
                :x => x,
                :y => y,
                :u => opt[:quiver][1],
                :v => opt[:quiver][2],
            ])
        end
    end
    series_plot = series_func(series_opt, table)
    push!(axis, series_plot)
    pgfx_add_legend!(axis, series, opt)
end

function pgfx_add_series!(::Val{:shape}, axis, series_opt, series, series_func, opt)
    push!(series_opt, "area legend" => nothing)
    series_opt = merge(series_opt, pgfx_fillstyle(opt))
    pgfx_add_series!(Val(:path), axis, series_opt, series, series_func, opt)
end

function pgfx_add_series!(::Val{:steppre}, axis, series_opt, series, series_func, opt)
    push!(series_opt, "const plot mark right" => nothing)
    pgfx_add_series!(Val(:path), axis, series_opt, series, series_func, opt)
end

function pgfx_add_series!(::Val{:stepmid}, axis, series_opt, series, series_func, opt)
    push!(series_opt, "const plot mark mid" => nothing)
    pgfx_add_series!(Val(:path), axis, series_opt, series, series_func, opt)
end

function pgfx_add_series!(::Val{:steppost}, axis, series_opt, series, series_func, opt)
    push!(series_opt, "const plot" => nothing)
    pgfx_add_series!(Val(:path), axis, series_opt, series, series_func, opt)
end

function pgfx_add_series!(::Val{:ysticks}, axis, series_opt, series, series_func, opt)
    push!(series_opt, "ycomb" => nothing)
    pgfx_add_series!(Val(:path), axis, series_opt, series, series_func, opt)
end

function pgfx_add_series!(::Val{:xsticks}, axis, series_opt, series, series_func, opt)
    push!(series_opt, "xcomb" => nothing)
    pgfx_add_series!(Val(:path), axis, series_opt, series, series_func, opt)
end

function pgfx_add_legend!(axis, series, opt, i = 1)
    if series[:subplot][:legend] != :none
        leg_entry = if opt[:label] isa AVec
            get(opt[:label], i, "")
        elseif opt[:label] isa AbstractString
            if i == 1
                get(opt, :label, "")
            else
                ""
            end
        else
            throw(ArgumentError("Malformed label. label = $(opt[:label])"))
        end
        if leg_entry == "" || !pgfx_should_add_to_legend(series)
            push!(axis.contents[end].options, "forget plot" => nothing)
        else
            leg_opt = PGFPlotsX.Options()
            legend = PGFPlotsX.LegendEntry(leg_opt, leg_entry, false)
            push!(axis, legend)
        end
    end
    return nothing
end

pgfx_series_arguments(series, opt, range) = (arg[range] for arg in pgfx_series_arguments(series, opt))
function pgfx_series_arguments(series, opt)
    st = series[:seriestype]
    return if st in (:contour, :contour3d)
        opt[:x], opt[:y], handle_surface(opt[:z])
    elseif st in (:heatmap, :surface, :wireframe)
        surface_to_vecs(opt[:x], opt[:y], opt[:z])
    elseif RecipesPipeline.is3d(st)
        opt[:x], opt[:y], opt[:z]
    elseif st == :straightline
        straightline_data(series)
    elseif st == :shape
        shape_data(series)
    elseif ispolar(series)
        theta, r = opt[:x], opt[:y]
        rad2deg.(theta), r
    else
        opt[:x], opt[:y]
    end
end
##
pgfx_get_linestyle(k) = get(
    (
        solid = "solid",
        dash = "dashed",
        dot = "dotted",
        dashdot = "dashdotted",
        dashdotdot = "dashdotdotted",
    ),
    Symbol(k),
    "solid",
)

pgfx_get_marker(k) = get(
    (
        none = "none",
        cross = "+",
        xcross = "x",
        + = "+",
        x = "x",
        utriangle = "triangle*",
        dtriangle = "triangle*",
        rtriangle = "triangle*",
        ltriangle = "triangle*",
        circle = "*",
        rect = "square*",
        star5 = "star",
        star6 = "asterisk",
        diamond = "diamond*",
        pentagon = "pentagon*",
        hline = "-",
        vline = "|",
    ),
    Symbol(k),
    "*",
)

pgfx_get_xguide_pos(k) = get(
    (
        top = "at={(0.5,1)},above,",
        right = "at={(ticklabel* cs:1.02)}, anchor=west,",
        left = "at={(ticklabel* cs:-0.02)}, anchor=east,",
    ),
    k,
    "at={(ticklabel cs:0.5)}, anchor=near ticklabel"
)

pgfx_get_yguide_pos(k) = get(
    (
        top = "at={(ticklabel* cs:1.02)}, anchor=south",
        right = "at={(1,0.5)},below,",
        bottom = "at={(ticklabel* cs:-0.02)}, anchor=north,",
    ),
    k,
    "at={(ticklabel cs:0.5)}, anchor=near ticklabel"
)

pgfx_get_legend_pos(k) = get(
    (
        top = ("at" => string((0.5, 0.98)), "anchor" => "north"),
        bottom = ("at" => string((0.5, 0.02)), "anchor" => "south"),
        left = ("at" => string((0.02, 0.5)), "anchor" => "west"),
        right = ("at" => string((0.98, 0.5)), "anchor" => "east"),
        bottomleft = ("at" => string((0.02, 0.02)), "anchor" => "south west"),
        bottomright = ("at" => string((0.98, 0.02)), "anchor" => "south east"),
        topright = ("at" => string((0.98, 0.98)), "anchor" => "north east"),
        topleft = ("at" => string((0.02, 0.98)), "anchor" => "north west"),
        outertop = ("at" => string((0.5, 1.02)), "anchor" => "south"),
        outerbottom = ("at" => string((0.5, -0.02)), "anchor" => "north"),
        outerleft = ("at" => string((-0.02, 0.5)), "anchor" => "east"),
        outerright = ("at" => string((1.02, 0.5)), "anchor" => "west"),
        outerbottomleft = ("at" => string((-0.02, -0.02)), "anchor" => "north east"),
        outerbottomright = ("at" => string((1.02, -0.02)), "anchor" => "north west"),
        outertopright = ("at" => string((1.02, 1)), "anchor" => "north west"),
        outertopleft = ("at" => string((-0.02, 1)), "anchor" => "north east"),
    ),
    Symbol(k),
    ("at" => string((1.02, 1)), "anchor" => "north west"),
)
pgfx_get_legend_pos(t::Tuple{S,T}) where {S<:Real,T<:Real} = ("at" => "{$(string(t))}", "anchor" => "north west")
pgfx_get_legend_pos(nt::NamedTuple) = ("at" => "{$(string(nt.at))}", "anchor" => string(nt.anchor))
pgfx_get_legend_pos(theta::Real) = pgfx_get_legend_pos((theta,:inner))
function pgfx_get_legend_pos(v::Tuple{S,Symbol}) where S <: Real
    (s,c) = sincosd(v[1])
    anchors = [
               "south west" "south" "south east";
               "west" "center" "east";
               "north west" "north" "north east";
              ]

    if v[2] === :inner
        rect = (0.07,0.5,1.0,0.07,0.52,1.0)
        anchor = anchors[legend_anchor_index(s),legend_anchor_index(c)]
    else
        rect = (-0.15,0.5,1.05,-0.15,0.52,1.1)
        anchor = anchors[4-legend_anchor_index(s),4-legend_anchor_index(c)]
    end
    return ("at"=>"$(string(legend_pos_from_angle(v[1],rect...)))", "anchor"=>anchor)
end

function pgfx_get_legend_style(sp)
    cstr = plot_color(sp[:background_color_legend])
    a = alpha(cstr)
    fg_alpha = alpha(plot_color(sp[:foreground_color_legend]))
    legfont = legendfont(sp)
    PGFPlotsX.Options(
        pgfx_linestyle(
            pgfx_thickness_scaling(sp),
            sp[:foreground_color_legend],
            fg_alpha,
            "solid",
        ) => nothing,
        "fill" => cstr,
        "fill opacity" => a,
        "text opacity" => alpha(plot_color(sp[:legendfontcolor])),
        "font" => pgfx_font(
            sp[:legendfontsize],
            pgfx_thickness_scaling(sp),
        ),
        "text" => plot_color(sp[:legendfontcolor]),
        "cells" => PGFPlotsX.Options("anchor" => get((left = "west", right = "east", hcenter = "center"), legfont.halign, "west")),
        pgfx_get_legend_pos(sp[:legend])...,
    )
end

pgfx_get_colorbar_pos(s) =
    get((left = " left", bottom = " horizontal", top = " horizontal"), s, "")
pgfx_get_colorbar_pos(b::Bool) = ""

pgfx_get_title_pos(s) =
    get((
        left = ("at" => "{(0,1)}", "anchor" => "south west"),
        right = ("at" => "{(1,1)}", "anchor" => "south east"),
        ), s,
        ("at" => "{(0.5,1)}", "anchor" => "south"))
pgfx_get_title_pos(t::Tuple) = ("at" => "{$(string(t))}", "anchor" => "south")
pgfx_get_title_pos(nt::NamedTuple) = ("at" => "{$(string(nt.at))}", "anchor" => string(nt.anchor))

function pgfx_get_ticklabel_style(sp, axis)
    cstr = plot_color(axis[:tickfontcolor])
    return PGFPlotsX.Options(
        "font" => pgfx_font(
            axis[:tickfontsize], pgfx_thickness_scaling(sp)
        ), "color" => cstr,
        "draw opacity" => alpha(cstr),
        "rotate" => axis[:tickfontrotation],
    )
end

function pgfx_get_colorbar_ticklabel_style(sp)
    cstr = plot_color(sp[:colorbar_tickfontcolor])
    return PGFPlotsX.Options(
        "font" => pgfx_font(
            sp[:colorbar_tickfontsize], pgfx_thickness_scaling(sp)
        ), "color" => cstr,
        "draw opacity" => alpha(cstr),
        "rotate" => sp[:colorbar_tickfontrotation],
    )
end
## --------------------------------------------------------------------------------------
pgfx_arrow(::Nothing) = "every arrow/.append style={-}"
function pgfx_arrow(arr::Arrow, side = arr.side)
    components = String[]
    head = String[]
    push!(head, "{stealth[length = $(arr.headlength)pt, width = $(arr.headwidth)pt")
    if arr.style == :open
        push!(head, ", open")
    end
    push!(head, "]}")
    head = join(head, "")
    if side == :both || side == :tail
        push!(components, head)
    end
    push!(components, "-")
    if side == :both || side == :head
        push!(components, head)
    end
    components = join(components, "")
    return "every arrow/.append style={$(components)}"
end

function pgfx_filllegend!(series_opt, opt)
    io = IOBuffer()
    PGFPlotsX.print_tex(io, pgfx_fillstyle(opt))
    style = strip(String(take!(io)), ['[', ']', ' '])
    push!(series_opt, "legend image code/.code" => """{
              \\draw[$style] (0cm,-0.1cm) rectangle (0.6cm,0.1cm);
              }""")
end

# Generates a colormap for pgfplots based on a ColorGradient
pgfx_colormap(cl::PlotUtils.AbstractColorList) = pgfx_colormap(color_list(cl))
function pgfx_colormap(v::Vector{<:Colorant})
    join(map(v) do c
        @sprintf("rgb=(%.8f,%.8f,%.8f)", red(c), green(c), blue(c))
    end, "\n")
end
function pgfx_colormap(cg::ColorGradient)
    join(map(1:length(cg)) do i
        @sprintf("rgb(%.8f)=(%.8f,%.8f,%.8f)", cg.values[i], red(cg.colors[i]), green(cg.colors[i]), blue(cg.colors[i]))
    end, "\n")
end

function pgfx_framestyle(style::Symbol)
    if style in (:box, :axes, :origin, :zerolines, :grid, :none)
        return style
    else
        default_style = get((semi = :box,), style, :axes)
        @warn "Framestyle :$style is not (yet) supported by the PGFPlotsX backend. :$default_style was chosen instead."
        default_style
    end
end

pgfx_thickness_scaling(plt::Plot) = plt[:thickness_scaling]
pgfx_thickness_scaling(sp::Subplot) = pgfx_thickness_scaling(sp.plt)
pgfx_thickness_scaling(series) = pgfx_thickness_scaling(series[:subplot])

function pgfx_fillstyle(plotattributes, i = 1)
    cstr = get_fillcolor(plotattributes, i)
    a = get_fillalpha(plotattributes, i)
    if a === nothing
        a = alpha(single_color(cstr))
    end
    PGFPlotsX.Options("fill" => cstr, "fill opacity" => a)
end

function pgfx_linestyle(linewidth::Real, color, α = 1, linestyle = :solid)
    cstr = plot_color(color, α)
    a = alpha(cstr)
    return PGFPlotsX.Options(
        "color" => cstr,
        "draw opacity" => a,
        "line width" => linewidth,
        pgfx_get_linestyle(linestyle) => nothing,
    )
end

function pgfx_linestyle(plotattributes, i = 1)
    lw = pgfx_thickness_scaling(plotattributes) * get_linewidth(plotattributes, i)
    lc = single_color(get_linecolor(plotattributes, i))
    la = get_linealpha(plotattributes, i)
    ls = get_linestyle(plotattributes, i)
    return pgfx_linestyle(lw, lc, la, ls)
end

function pgfx_font(fontsize, thickness_scaling = 1, font = "\\selectfont")
    fs = fontsize * thickness_scaling
    return string("{\\fontsize{", fs, " pt}{", 1.3fs, " pt}", font, "}")
end

# If a particular fontsize parameter is `nothing`, produce a figure that doesn't specify the
# font size, and therefore uses whatever fontsize is utilised by the doc in which the
# figure is located.
function pgfx_font(fontsize::Nothing, thickness_scaling = 1, font = "\\selectfont")
    return string("{", font, "}")
end

function pgfx_should_add_to_legend(series::Series)
    series.plotattributes[:primary] &&
    !(
        series.plotattributes[:seriestype] in (
            :hexbin,
            :bins2d,
            :histogram2d,
            :hline,
            :vline,
            :contour,
            :contourf,
            :contour3d,
            :heatmap,
            :image,
        )
    )
end

function pgfx_marker(plotattributes, i = 1)
    shape = _cycle(plotattributes[:markershape], i)
    cstr =
        plot_color(get_markercolor(plotattributes, i), get_markeralpha(plotattributes, i))
    a = alpha(cstr)
    cstr_stroke = plot_color(
        get_markerstrokecolor(plotattributes, i),
        get_markerstrokealpha(plotattributes, i),
    )
    a_stroke = alpha(cstr_stroke)
    mark_size =
        pgfx_thickness_scaling(plotattributes) *
        0.75 *
        _cycle(plotattributes[:markersize], i)
    mark_freq = !any(isnan, plotattributes[:y]) && plotattributes[:markershape] isa AVec ?
                length(plotattributes[:markershape]) : 1
    return PGFPlotsX.Options(
        "mark" => shape isa Shape ? "PlotsShape$i" : pgfx_get_marker(shape),
        "mark size" => "$mark_size pt",
        "mark repeat" => mark_freq,
        "mark options" => PGFPlotsX.Options(
            "color" => cstr_stroke,
            "draw opacity" => a_stroke,
            "fill" => cstr,
            "fill opacity" => a,
            "line width" =>
                pgfx_thickness_scaling(plotattributes) *
                0.75 *
                _cycle(plotattributes[:markerstrokewidth], i),
            "rotate" => if shape == :dtriangle
                180
            elseif shape == :rtriangle
                270
            elseif shape == :ltriangle
                90
            else
                0
            end,
            pgfx_get_linestyle(_cycle(plotattributes[:markerstrokestyle], i)) => nothing,
        ),
    )
end

function pgfx_add_annotation!(o, x, y, val, thickness_scaling = 1; cs = "axis cs:", options = PGFPlotsX.Options())
    # Construct the style string.
    cstr = val.font.color
    a = alpha(cstr)
    push!(
        o,
        join([
            "\\node",
            sprint(PGFPlotsX.print_tex, merge(
                PGFPlotsX.Options(
                    get((hcenter = "", left = "right", right = "left"), val.font.halign, "") =>
                        nothing,
                    "color" => cstr,
                    "draw opacity" => convert(Float16, a),
                    "rotate" => val.font.rotation,
                    "font" => pgfx_font(val.font.pointsize, thickness_scaling),
                ),
                options
            )),
            string(" at (", cs, x, ",", y, ") {", val.str, "};"),
        ]),
    )
end

function pgfx_add_ribbons!(axis, series, segment_plot, series_func, series_index)
    ribbon_y = series[:ribbon]
    opt = series.plotattributes
    if ribbon_y isa AVec
        ribbon_n = length(opt[:y]) ÷ length(ribbon_y)
        ribbon_yp = ribbon_ym = repeat(ribbon_y, outer = ribbon_n)
    elseif ribbon_y isa Tuple
        ribbon_ym, ribbon_yp = ribbon_y
        ribbon_nm = length(opt[:y]) ÷ length(ribbon_ym)
        ribbon_ym = repeat(ribbon_ym, outer = ribbon_nm)
        ribbon_np = length(opt[:y]) ÷ length(ribbon_yp)
        ribbon_yp = repeat(ribbon_yp, outer = ribbon_np)
    else
        ribbon_yp = ribbon_ym = ribbon_y
    end
    # upper ribbon
    rib_uuid = uuid4()
    ribbon_name_plus = "plots_rib_p$rib_uuid"
    ribbon_opt_plus = merge(
        segment_plot.options,
        PGFPlotsX.Options(
            "name path" => ribbon_name_plus,
            "color" => opt[:fillcolor],
            "draw opacity" => opt[:fillalpha],
            "forget plot" => nothing,
        ),
    )
    coordinates_plus = PGFPlotsX.Coordinates(opt[:x], opt[:y] .+ ribbon_yp)
    ribbon_plot_plus = series_func(ribbon_opt_plus, coordinates_plus)
    push!(axis, ribbon_plot_plus)
    # lower ribbon
    ribbon_name_minus = "plots_rib_m$rib_uuid"
    ribbon_opt_minus = merge(
        segment_plot.options,
        PGFPlotsX.Options(
            "name path" => ribbon_name_minus,
            "color" => opt[:fillcolor],
            "draw opacity" => opt[:fillalpha],
            "forget plot" => nothing,
        ),
    )
    coordinates_plus = PGFPlotsX.Coordinates(opt[:x], opt[:y] .- ribbon_ym)
    ribbon_plot_plus = series_func(ribbon_opt_minus, coordinates_plus)
    push!(axis, ribbon_plot_plus)
    # fill
    push!(
        axis,
        series_func(
            merge(
                pgfx_fillstyle(opt, series_index),
                PGFPlotsX.Options("forget plot" => nothing),
            ),
            "fill between [of=$(ribbon_name_plus) and $(ribbon_name_minus)]",
        ),
    )
    return axis
end

function pgfx_fillrange_series!(axis, series, series_func, i, fillrange, rng)
    fillrange_opt = PGFPlotsX.Options("line width" => "0", "draw opacity" => "0")
    fillrange_opt = merge(fillrange_opt, pgfx_fillstyle(series, i))
    push!(fillrange_opt, "mark" => "none") # no markers on fillranges
    push!(fillrange_opt, "forget plot" => nothing)
    opt = series.plotattributes
    args = RecipesPipeline.is3d(series) ? (opt[:x][rng], opt[:y][rng], opt[:z][rng]) :
        (opt[:x][rng], opt[:y][rng])
    push!(
        axis,
        PGFPlotsX.PlotInc(fillrange_opt, pgfx_fillrange_args(fillrange, args...)),
    )
    return axis
end

function pgfx_fillrange_args(fillrange, x, y)
    n = length(x)
    x_fill = [x; x[n:-1:1]; x[1]]
    y_fill = [y; _cycle(fillrange, n:-1:1); y[1]]
    return PGFPlotsX.Coordinates(x_fill, y_fill)
end

function pgfx_fillrange_args(fillrange, x, y, z)
    n = length(x)
    x_fill = [x; x[n:-1:1]; x[1]]
    y_fill = [y; y[n:-1:1]; x[1]]
    z_fill = [z; _cycle(fillrange, n:-1:1); z[1]]
    return PGFPlotsX.Coordinates(x_fill, y_fill, z_fill)
end

function pgfx_sanitize_string(p::PlotText)
    PlotText(pgfx_sanitize_string(p.str), p.font)
end
function pgfx_sanitize_string(s::AbstractString)
    s = replace(s, r"\\?\#" => "\\#")
    s = replace(s, r"\\?\%" => "\\%")
    s = replace(s, r"\\?\_" => "\\_")
    s = replace(s, r"\\?\&" => "\\&")
    s = replace(s, r"\\?\{" => "\\{")
    s = replace(s, r"\\?\}" => "\\}")
    s = map(split(s, "")) do s
            isascii(s) ? s : latexify(s)
            end |> join
end
@require LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" begin
    using .LaTeXStrings
    function pgfx_sanitize_string(s::LaTeXString)
        s = replace(s, r"\\?\#" => "\\#")
        s = replace(s, r"\\?\%" => "\\%")
        return LaTeXString(s)
    end
end
function pgfx_sanitize_plot!(plt)
        for (key, value) in plt.attr
            if value isa Union{AbstractString, AbstractVector{<:AbstractString}}
                plt.attr[key] = pgfx_sanitize_string.(value)
            end
        end
        for subplot in plt.subplots
            for (key, value) in subplot.attr
                if key == :annotations && subplot.attr[:annotations] !== nothing
                    old_ann = subplot.attr[key]
                    for i in eachindex(old_ann)
                        subplot.attr[key][i] = (old_ann[i][1], old_ann[i][2], pgfx_sanitize_string(old_ann[i][3]))
                    end
                elseif value isa Union{AbstractString, AbstractVector{<:AbstractString}}
                    subplot.attr[key] = pgfx_sanitize_string.(value)
                end
            end
        end
        for series in plt.series_list
            for (key, value) in series.plotattributes
                if key == :series_annotations && series.plotattributes[:series_annotations] !== nothing
                    old_ann = series.plotattributes[key].strs
                    for i in eachindex(old_ann)
                        series.plotattributes[key].strs[i] = pgfx_sanitize_string(old_ann[i])
                    end
                elseif value isa Union{AbstractString, AbstractVector{<:AbstractString}}
                    series.plotattributes[key] = pgfx_sanitize_string.(value)
                end
            end
        end
        ##
end
# --------------------------------------------------------------------------------------
function pgfx_axis!(opt::PGFPlotsX.Options, sp::Subplot, letter)
    axis = sp[Symbol(letter, :axis)]

    # turn off scaled ticks
    push!(
        opt,
        "scaled $(letter) ticks" => "false",
        string(letter, :label) => axis[:guide],
    )
    tick_color = plot_color(axis[:foreground_color_axis])
    push!(opt,
        "$(letter) tick style" => PGFPlotsX.Options(
            "color" => color(tick_color),
            "opacity" => alpha(tick_color),
        ),
    )
    tick_label_color = plot_color(axis[:tickfontcolor])
    push!(opt,
        "$(letter) tick label style" => PGFPlotsX.Options(
            "color" => color(tick_color),
            "opacity" => alpha(tick_color),
            "rotate" => axis[:rotation]
        ),
    )

    # set to supported framestyle
    framestyle = pgfx_framestyle(sp[:framestyle] == false ? :none : sp[:framestyle])

    # axis label position
    labelpos = ""
    if letter == :x
        labelpos = pgfx_get_xguide_pos(axis[:guide_position])
    elseif letter == :y
        labelpos = pgfx_get_yguide_pos(axis[:guide_position])
    end

    # Add label font
    cstr = plot_color(axis[:guidefontcolor])
    α = alpha(cstr)
    push!(
        opt,
        string(letter, "label style") => PGFPlotsX.Options(
            labelpos => nothing,
            "font" =>     pgfx_font(axis[:guidefontsize], pgfx_thickness_scaling(sp)),
            "color" => cstr,
            "draw opacity" => α,
            "rotate" => axis[:guidefontrotation],
        ),
    )

    # flip/reverse?
    axis[:flip] && push!(opt, "$letter dir" => "reverse")

    # scale
    scale = axis[:scale]
    if scale in (:log2, :ln, :log10)
        push!(opt, string(letter, :mode) => "log")
        scale == :ln ||
        push!(opt, "log basis $letter" => "$(scale == :log2 ? 2 : 10)")
    end

    # ticks on or off
    if axis[:ticks] in (nothing, false, :none) || framestyle == :none
        push!(opt, "$(letter)majorticks" => "false")
    elseif framestyle in (:grid, :zerolines)
        push!(opt, "$letter tick style" => PGFPlotsX.Options("draw" => "none"))
    end

    # grid on or off
    if axis[:grid] && framestyle != :none
        push!(opt, "$(letter)majorgrids" => "true")
    else
        push!(opt, "$(letter)majorgrids" => "false")
    end

    # limits
    lims = ispolar(sp) && letter == :x ? rad2deg.(axis_limits(sp, :x)) :
        axis_limits(sp, letter)
    push!(opt, string(letter, :min) => lims[1], string(letter, :max) => lims[2])

    if !(axis[:ticks] in (nothing, false, :none, :native)) && framestyle != :none
        # ticks
        ticks = get_ticks(sp, axis)
        #pgf plot ignores ticks with angle below 90 when xmin = 90 so shift values
        tick_values =
            ispolar(sp) && letter == :x ? [rad2deg.(ticks[1])[3:end]..., 360, 405] :
            ticks[1]
        push!(
            opt,
            string(letter, "tick") => string("{", join(tick_values, ","), "}"),
        )
        if axis[:showaxis] &&
           axis[:scale] in (:ln, :log2, :log10) && axis[:ticks] == :auto
            # wrap the power part of label with }
            tick_labels = Vector{String}(undef, length(ticks[2]))
            for (i, label) in enumerate(ticks[2])
                base, power = split(label, "^")
                power = string("{", power, "}")
                tick_labels[i] = string(base, "^", power)
            end
            push!(
                opt,
                string(letter, "ticklabels") =>
                    string("{\$", join(tick_labels, "\$,\$"), "\$}"),
            )
        elseif axis[:showaxis]
            tick_labels =
                ispolar(sp) && letter == :x ? [ticks[2][3:end]..., "0", "45"] :
                ticks[2]
            if axis[:formatter] in (:scientific, :auto)
                tick_labels = string.("\$", convert_sci_unicode.(tick_labels), "\$")
                tick_labels = replace.(tick_labels, Ref("×" => "\\times"))
            end
            push!(
                opt,
                string(letter, "ticklabels") =>
                    string("{", join(tick_labels, ","), "}"),
            )
        else
            push!(opt, string(letter, "ticklabels") => "{}")
        end
        if axis[:tick_direction] === :none
            push!(
                opt,
                string(letter, "tick style") => "draw=none",
            )
        else
            push!(
                opt,
                string(letter, "tick align") =>
                    (axis[:tick_direction] == :out ? "outside" : "inside"),
            )
        end
        push!(
            opt, string(letter, "ticklabel style") => pgfx_get_ticklabel_style(sp, axis)
        )
        push!(
            opt,
            string(letter, " grid style") => pgfx_linestyle(
                pgfx_thickness_scaling(sp) * axis[:gridlinewidth],
                axis[:foreground_color_grid],
                axis[:gridalpha],
                axis[:gridstyle],
            ),
        )

        # minor ticks
        # NOTE: PGFPlots would provide "minor x ticks num", but this only places minor ticks
        #       between major ticks and not outside first and last tick to the axis limits.
        #       Hence, we hack around with extra ticks. Unfortunately this conflicts with
        #       `:zerolines` framestyle hack. So minor ticks are not working with
        #       `:zerolines`.
        minor_ticks = get_minor_ticks(sp, axis, ticks)
        if minor_ticks !== nothing
            minor_ticks =
                ispolar(sp) && letter == :x ? [rad2deg.(minor_ticks)[3:end]..., 360, 405] :
                minor_ticks
            push!(
                opt,
                string("extra ", letter, " ticks") => string("{", join(minor_ticks, ","), "}"),
            )
            push!(opt, string("extra ", letter, " tick labels") => "")
            push!(
                opt,
                string("extra ", letter, " tick style") => PGFPlotsX.Options(
                    "grid" => axis[:minorgrid] ? "major" : "none",
                    string(letter, " grid style") => pgfx_linestyle(
                        pgfx_thickness_scaling(sp) * axis[:minorgridlinewidth],
                        axis[:foreground_color_minor_grid],
                        axis[:minorgridalpha],
                        axis[:minorgridstyle],
                    ),
                    "major tick length" => typeof(axis[:minorticks]) <: Integer && axis[:minorticks] > 1 || axis[:minorticks] ? "0.1cm" : "0"
                ),
            )
        end
    end

    # framestyle
    if framestyle in (:axes, :origin)
        axispos = axis[:mirror] ? "right" : framestyle == :axes ? "left" : "middle"

        if axis[:draw_arrow]
            push!(opt, string("axis ", letter, " line") => axispos)
        else
            # the * after line disables the arrow at the axis
            push!(opt, string("axis ", letter, " line*") => axispos)
        end
    end

    if framestyle == :zerolines
        push!(opt, string("extra ", letter, " ticks") => "0")
        push!(opt, string("extra ", letter, " tick labels") => "")
        push!(
            opt,
            string("extra ", letter, " tick style") => PGFPlotsX.Options(
                "grid" => "major",
                string(letter, " grid style") => pgfx_linestyle(
                    pgfx_thickness_scaling(sp),
                    axis[:foreground_color_border],
                    1.0,
                ),
            ),
        )
    end

    if !axis[:showaxis]
        push!(opt, "separate axis lines")
    end
    if !axis[:showaxis] || framestyle in (:zerolines, :grid, :none)
        push!(opt, string(letter, " axis line style") => "{draw opacity = 0}")
    else
        push!(
            opt,
            string(letter, " axis line style") => pgfx_linestyle(
                pgfx_thickness_scaling(sp),
                axis[:foreground_color_border],
                1.0,
            ),
        )
    end
end

# --------------------------------------------------------------------------------------
# display calls this and then _display, its called 3 times for plot(1:5)
# Set the (left, top, right, bottom) minimum padding around the plot area
# to fit ticks, tick labels, guides, colorbars, etc.
function _update_min_padding!(sp::Subplot{PGFPlotsXBackend})
    leg = sp[:legend]
    if leg in (:best, :outertopright, :outerright, :outerbottomright) || (leg isa Tuple && leg[1] >= 1)
        sp.minpad = (0mm, 0mm, 5mm, 0mm)
    else
        sp.minpad = (0mm, 0mm, 0mm, 0mm)
    end
end

function _create_backend_figure(plt::Plot{PGFPlotsXBackend})
    plt.o = PGFPlotsXPlot()
end

function _series_added(plt::Plot{PGFPlotsXBackend}, series::Series)
    plt.o.is_created = false
end

function _update_plot_object(plt::Plot{PGFPlotsXBackend})
    plt.o(plt)
end

for mime in ("application/pdf", "image/svg+xml")
    @eval function _show(
        io::IO,
        mime::MIME{Symbol($mime)},
        plt::Plot{PGFPlotsXBackend},
    )
        plt.o.was_shown = true
        show(io, mime, plt.o.the_plot)
    end
end

    function _show(
        io::IO,
        mime::MIME{Symbol("image/png")},
        plt::Plot{PGFPlotsXBackend},
    )
        plt.o.was_shown = true
        plt_file = tempname() * ".png"
        PGFPlotsX.pgfsave(plt_file, plt.o.the_plot; dpi=plt[:dpi])
        write(io, read(plt_file))
        rm(plt_file; force = true)
    end

function _show(
    io::IO,
    mime::MIME{Symbol("application/x-tex")},
    plt::Plot{PGFPlotsXBackend},
)
    plt.o.was_shown = true
    PGFPlotsX.print_tex(
        io,
        plt.o.the_plot,
        include_preamble = plt.attr[:tex_output_standalone],
    )
end

function _display(plt::Plot{PGFPlotsXBackend})
    plt.o.was_shown = true
    display(PGFPlotsX.PGFPlotsXDisplay(), plt.o.the_plot)
end


text:
neurons1=200 # sub-sampling for each digit
inds_all=Array{Int64}(undef,0)
items_per_neuron=Array{Int64}(undef,0)
@time for i=0:9
    global inds_all,items_per_neuron
    println(Dates.format(now(), "HH:MM:SS")," ",i)
    i1=yy_train_all.==i
    inds_all1,n_per_part1,items_per_neuron1=ANNBN.___clustering(neurons1,xx_train[i1,:],300)
    ia1=(1:i_train)[i1]
    inds_all=[inds_all;ia1[inds_all1]]
    items_per_neuron=[items_per_neuron;items_per_neuron1]
end
# plot(unique(sort(inds_all)))
n_per_part=[0;cumsum(items_per_neuron)]
neurons=10neurons1 # all neurons
plot(items_per_neuron)
minimum(items_per_neuron)
# xx_train=xx_train[inds_all,:]
# yy_train=yy_train[inds_all]
# yy_train_all=yy_train_all[inds_all]
# inds_all=1:i_train

ext: jl


text:
using Test
using MixedModelsSim
using Random
using DSP
include("test_utilities.jl")
τ = 1.5
fs = 18
evt,epoch_dat = simulate_lmm(MersenneTwister(1), τ,fs,  β=[0.,-1.],
    σs=[1,1,1],
    σ=1,
    n_sub = 20,
    n_item = 30,
    noise_type="AR-exponential");

#mres,res = fit(UnfoldLinearMixedModel,f,evt,dat[chIx:chIx,:,:].*1e6 ,times,contrasts=Dict(:category => EffectsCoding(), :condition => EffectsCoding()));

times = range(0,stop=τ-1/fs,step=1/fs)

f2 = @formula 0 ~ 1 + stimType + (1+stimType|subj) + (1|item)
mres,res = fit(UnfoldLinearMixedModel,f2,evt,epoch_dat,times)

p_df = cluster_permutation_test(mres,epoch_dat,times,2)

@test p_df[1,:pval] == 0.036

#--- Testing cluster_permutation call
tRange = 1:length(times)
nPerm = 10
permDat = cluster_permutation(mres,epoch_dat,tRange,2,nPerm)
@test size(permDat,2) == nPerm
@test size(permDat,1) == length(tRange)

tRange = 5:length(times)-5
permDat = cluster_permutation(mres,epoch_dat,tRange,2,nPerm)
@test size(permDat,1) == length(tRange)

#------ testing pymne_cluster
z_obs = [m.z for m in coefpvalues(mres.modelinfo) if String(m.coefname)==coefnames(mres.X.formulas)[2][coeffOfInterest]]

#-- clusterpermutation
obs_cluster = pymne_cluster(z_obs,4)
@test all(abs.(obs_cluster[2]) .> 4) # clusters have to be larger 4
@test obs_cluster[1][1][1].start == 13 # start of cluster
@test obs_cluster[1][1][1].stop == 14 # end of cluster

obs_cluster = pymne_cluster(z_obs,0.1)
@test all(abs.(obs_cluster[2]) .> 0.1) # all clusters have to be larger 0.1
@test length(obs_cluster[2]) == 12 # if nothing changed, 12 clusters should be found

#-- tfce
clusterFormingThreshold = Dict(:start=>0,:step=>0.2)
obs_cluster = pymne_cluster(z_obs,clusterFormingThreshold)
@test all(obs_cluster[2].>=0) # TFCE is always positive
@test findmax(obs_cluster[2])[2] == findmax(abs.(z_obs))[2] # tfce should not move maximum

#gen_han(τ,fs,3)
if 1 == 0
    # plot some design bits
    map = mapping(:subj,:dv,color=:stimType)#,linestyle=:channel)
    AlgebraOfGraphics.data(evt) * visual(Scatter) * map  |> draw
end


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

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 el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 8dbbaab7-c8a8-495a-b34e-2df3ee51be3f
using Plots

# ╔═╡ 44427f8f-4a65-44da-9062-cce6e75e2569
html"
Modeling the Scale Factor

The acceleration equation is given by"

# ╔═╡ ecd37f62-bc2d-4256-95e9-b29424a246d6
md"``(\frac{ä}{a})² = -\frac{4πG}{3c²} (ϵ(t) + 3P(t)) + \frac{Λ}{3}``"

# ╔═╡ dcdb469e-1ec3-4fec-8bf1-d429719dcd9b
md"
Assume the energy density is dominated by non-relativistic matter, i.e. ``ϵ ≈ ρc²`` and ``P << ϵ``."

# ╔═╡ 2890c2ad-8409-45bf-8fd4-a6cd426ab55d
md"For this case, ``ρ(t) ≈ ρ₀(a(t)=a₀)^{-3}``, where ``a₀ = 1`` is the scale factor of the universe today (when ``t = t₀``)."

# ╔═╡ 316b59fa-ead6-4a1a-ab2f-c2ded137a0f2
md"Recall that the critical density ``ρ_{crit} = \frac{3H₀²}{8πG}``."

# ╔═╡ 9cee1d18-0442-4dae-9e80-de38caa6163c
html"Thus:"

# ╔═╡ 9da1395a-6e7d-4bbe-a76d-5b35e6459b5a
md"``ä(t) = -\frac{ρ₀}{2ρ_{crit}} (\frac{H₀}{a(t)})² + \frac{a(t)Λ}{3}``,"

# ╔═╡ 8a031c70-9d29-44ce-80fb-028c560cb51a
html"or"

# ╔═╡ ba2770bd-33a1-40f8-ae37-9a25a082b4af
md"→ ``\frac{d²}{dy²}a(y) = -\frac{1}{2} Ωₘa(y)^{-2} + Ω_λa(y)``,"

# ╔═╡ ba921154-fce7-4297-8980-eee506b6d653
md"where ``Ωₘ = \frac{ρ₀}{ρ_{crit}}`` and ``Ω_λ = \frac{Λ}{3H₀²}`` and"

# ╔═╡ 0f0a3b82-c424-4c20-ae3f-824514fe01d6
md"``y = H₀t`` ⇛ ``\frac{d²}{dy²} a(y) = \frac{ä}{H₀²}``."

# ╔═╡ 50ff4ea4-8a4f-4179-9fc5-d554d4af4fcc
md"So first, let's create a model for any two given values ``Ωₘ`` and ``Ω_λ``:"

# ╔═╡ 400a2e78-afdb-11eb-2df6-7ff8a8e36528
@Base.kwdef mutable struct CosmoModel
	a  = 1.0
	da = 1.0
	Ωm = 0.3
	Ωλ = 0.7
	y  = 0.0
end

# ╔═╡ aae7f619-f21e-486f-88fd-6bb051c3b9a5
md"Recall that, today, ``y=0``, and"

# ╔═╡ eeb191c7-1ed5-4e8b-a408-9c4129166559
md"``a(y=0)=a₀≝1``. In addition, since"

# ╔═╡ ce002252-c4bc-4704-974f-73ff0b7139d7
md"``H₀≝\frac{ȧ}{a}``,"

# ╔═╡ 21e11baa-5239-4f05-83cf-bea77738e765
md"``\frac{da(y=0)}{dy}≝1``."

# ╔═╡ 844805b8-3e6e-496d-bce3-b33c6d26012a
md"Let's create a step function that advances this cosmological model:"

# ╔═╡ 0f9a7e58-73bd-42bd-b77c-3a3e670a5de8
function step!(c::CosmoModel, dy::Real=0.01)
	# Euler's method, what else?
	dda = -0.5*c.Ωm/c.a^2 + c.Ωλ*c.a
	c.da += dda * dy
	c.a += c.da * dy
	c.y += dy
end

# ╔═╡ 1ece7af1-2cd6-4937-9e4f-2a064d1dc859
md"We're going to measure ``-2≤y≤2``, so let's prepare the arrays for a and y, where values corresponding to ``y=0`` are going to fall squarely in the middle of the array:"

# ╔═╡ 539643eb-e04d-4515-b416-1e815b73fd77
begin
aplot = zeros(401, 3)
aplot[201, :] = [1.0 1.0 1.0] # Initial Conditions - let's just record immediately
yplot = -2:0.01:2
end

# ╔═╡ e01b33a4-f334-4633-a83c-93f5d67d2da2
html"Backward time integration:"

# ╔═╡ 070c5d93-09fc-42a5-b90b-73ea44cdad40
let
cosmos = [CosmoModel(Ωλ=0), CosmoModel(Ωm=0), CosmoModel()]
for (i, cosmo) ∈ enumerate(cosmos), y ∈ 200:-1:1
	step!(cosmo, -0.01) 	# dy < 0 for backwards integration
	aplot[y, i] = cosmo.a
end
end

# ╔═╡ 74ed32d4-6917-405e-b501-de1521070b53
html"Forward time integration:"

# ╔═╡ c44b31f2-d9e6-433f-a21f-a22ea78f9568
let
cosmos = [CosmoModel(Ωλ=0), CosmoModel(Ωm=0), CosmoModel()]
for (i, cosmo) ∈ enumerate(cosmos), y ∈ 202:401
	step!(cosmo)
	aplot[y, i] = cosmo.a
end
end

# ╔═╡ bff766db-188e-484d-8b53-ebd71d00e120
md"For each model, finding the ``y`` value at which ``a(y)=0`` will determine the estimated age of the universe:"

# ╔═╡ c756b29e-2a4a-4e4c-adce-fc633d3bc21e
izeros = [findfirst(≥(0), aplot[:, i]) for i ∈ 1:3]

# ╔═╡ 729c98bd-18a6-4395-abea-e21d899a9324
yzeros = [yplot[izero] for izero ∈ izeros]

# ╔═╡ d82029ce-8163-4de2-bdac-6ffa0a867580
let
start = min(izeros...)
y = yplot[start:end]
a = aplot[start:end, :]
plot(y, a[:, 1], lw=2, label="Ωλ=0", legend=:topleft)
plot!(y, a[:, 2], lw=2, label="Ωₘ=0")
plot!(y, a[:, 3], lw=2, label="Benchmark Model")
plot!(ylim=(0,6), xlabel="y", ylabel="a(y)", title="Scale Factor a(y) vs. y")
vline!(yzeros, label="Zeros: y = $(sort(yzeros))", linestyle=:dash)
annotate!(0, 2, "Today, y=0")
end

# ╔═╡ 3e5601e8-88b3-4062-902a-09518b592e4e
html"H0<math>≈70</math> km/s/Mpc. Converting units to 1/s:"

# ╔═╡ dec2b537-cc4d-4075-a0a0-7e27b25998ca
H₀ = 70 * 3.2407792700054e-20 	# 1 / s

# ╔═╡ ba3010b9-41ee-4532-a7f5-2992c89abb97
html"seconds2years will convert a time in seconds to a time in years:"

# ╔═╡ fdbd0b7e-26ad-45de-b42d-c407e120787b
seconds2years(s::Real) = s / (60*60*24*365)

# ╔═╡ cdefb2ec-7a6d-40a9-90ea-3c1a6bd44904
md"Converting each ``y`` value to **billions of years**:"

# ╔═╡ 1d655ed7-9d23-47bf-a204-84d4eea9e3ea
bya = seconds2years.([abs(y)/H₀ for y ∈ yzeros]) * 10^-9

# ╔═╡ c459cb1e-0a34-4b6e-b41d-b4da1ad02c33
("Ωₘ=0:", bya[1]), ("Ωλ=0:", bya[2]), ("Benchmark:", bya[3])

# ╔═╡ c0af49f7-f4c2-4ac3-b9b1-157dd0dd1ca6
md"What about models which do not expand forever?"

# ╔═╡ 998cf5a1-6888-4233-9310-6d4cc8c89537
md"Let's create a function that just forward-integrates our cosmological model:"

# ╔═╡ 5f785097-5fe2-436a-98b6-d8aaf4e1a607
function runsimulation(c::CosmoModel, maxsteps::Int=10_000)
	y = [c.y]; a = [c.a]
	for _ ∈ 1:maxsteps
		step!(c)
		c.a ≤ 0 && continue
		push!(y, c.y); push!(a, c.a)
	end
	y, a
end

# ╔═╡ 4b8b8af4-2940-40c9-be7b-2743c85b79ac
let
cosmo = CosmoModel(Ωm=1.5, Ωλ=0.0)
y, a = runsimulation(cosmo)
ymax_year = round(seconds2years(y[end]/H₀) * 10^-9, digits=2)
plot(y, a, lw=2, xlabel='y', ylabel="Scale Factor a(y)", label="Ωₘ=1.5, Ωλ=0")
title!("The 'Big Crunch' in $ymax_year Billion Years")
end

# ╔═╡ d7d10fdf-a202-4b1d-a401-d2606fbb2817
md"If ``Ωₘ = 1.5``, what is the minimum value of ``Ω_λ`` such that the universe does not collapse back in on itself?"

# ╔═╡ 5687f633-f671-4e90-bf73-12e5e26a3446
@bind Ωₛ html"<input type='range' min='0.00907228' max='0.00907236' step='0.00000001' value='0.00907232'>"

# ╔═╡ cbf27243-05d8-44e2-a764-97e97ddbda19
let
cosmo = CosmoModel(Ωm=1.5, Ωλ=Ωₛ)
y, a = runsimulation(cosmo)
plot(y, a, lw=2, xlabel='y', ylabel="Scale Factor a(y)", label="Ωₘ=1.5, Ωλ=$Ωₛ")
title!("a(y) vs y")
end

# ╔═╡ 65ebaec8-facf-4c8e-8a93-b5d4cd239233
Ωₛ

# ╔═╡ Cell order:
# ╟─44427f8f-4a65-44da-9062-cce6e75e2569
# ╟─ecd37f62-bc2d-4256-95e9-b29424a246d6
# ╟─dcdb469e-1ec3-4fec-8bf1-d429719dcd9b
# ╟─2890c2ad-8409-45bf-8fd4-a6cd426ab55d
# ╟─316b59fa-ead6-4a1a-ab2f-c2ded137a0f2
# ╟─9cee1d18-0442-4dae-9e80-de38caa6163c
# ╟─9da1395a-6e7d-4bbe-a76d-5b35e6459b5a
# ╟─8a031c70-9d29-44ce-80fb-028c560cb51a
# ╟─ba2770bd-33a1-40f8-ae37-9a25a082b4af
# ╟─ba921154-fce7-4297-8980-eee506b6d653
# ╟─0f0a3b82-c424-4c20-ae3f-824514fe01d6
# ╟─50ff4ea4-8a4f-4179-9fc5-d554d4af4fcc
# ╠═400a2e78-afdb-11eb-2df6-7ff8a8e36528
# ╟─aae7f619-f21e-486f-88fd-6bb051c3b9a5
# ╟─eeb191c7-1ed5-4e8b-a408-9c4129166559
# ╟─ce002252-c4bc-4704-974f-73ff0b7139d7
# ╟─21e11baa-5239-4f05-83cf-bea77738e765
# ╟─844805b8-3e6e-496d-bce3-b33c6d26012a
# ╠═0f9a7e58-73bd-42bd-b77c-3a3e670a5de8
# ╟─1ece7af1-2cd6-4937-9e4f-2a064d1dc859
# ╠═539643eb-e04d-4515-b416-1e815b73fd77
# ╟─e01b33a4-f334-4633-a83c-93f5d67d2da2
# ╠═070c5d93-09fc-42a5-b90b-73ea44cdad40
# ╟─74ed32d4-6917-405e-b501-de1521070b53
# ╠═c44b31f2-d9e6-433f-a21f-a22ea78f9568
# ╟─bff766db-188e-484d-8b53-ebd71d00e120
# ╠═c756b29e-2a4a-4e4c-adce-fc633d3bc21e
# ╠═729c98bd-18a6-4395-abea-e21d899a9324
# ╠═8dbbaab7-c8a8-495a-b34e-2df3ee51be3f
# ╠═d82029ce-8163-4de2-bdac-6ffa0a867580
# ╟─3e5601e8-88b3-4062-902a-09518b592e4e
# ╠═dec2b537-cc4d-4075-a0a0-7e27b25998ca
# ╟─ba3010b9-41ee-4532-a7f5-2992c89abb97
# ╠═fdbd0b7e-26ad-45de-b42d-c407e120787b
# ╟─cdefb2ec-7a6d-40a9-90ea-3c1a6bd44904
# ╠═1d655ed7-9d23-47bf-a204-84d4eea9e3ea
# ╠═c459cb1e-0a34-4b6e-b41d-b4da1ad02c33
# ╟─c0af49f7-f4c2-4ac3-b9b1-157dd0dd1ca6
# ╟─998cf5a1-6888-4233-9310-6d4cc8c89537
# ╠═5f785097-5fe2-436a-98b6-d8aaf4e1a607
# ╠═4b8b8af4-2940-40c9-be7b-2743c85b79ac
# ╟─d7d10fdf-a202-4b1d-a401-d2606fbb2817
# ╠═cbf27243-05d8-44e2-a764-97e97ddbda19
# ╟─5687f633-f671-4e90-bf73-12e5e26a3446
# ╠═65ebaec8-facf-4c8e-8a93-b5d4cd239233

max_stars_repo_name: benharmonics/GeneralPhysics

max_issues_repo_name: benharmonics/GeneralPhysics

max_forks_repo_name: benharmonics/GeneralPhysics


text:
const DEFAULT_SAMPLE_SIZE = 1000

export sourceParticles

import MonteCarloMeasurements
using MonteCarloMeasurements: Particles, StaticParticles, AbstractParticles

foreach([<=, >=, <, >]) do cmp
    MonteCarloMeasurements.register_primitive(cmp, eval)
end

export particles

particles(v::Vector{<:Real}) = Particles(v)
particles(d, N::Int=DEFAULT_SAMPLE_SIZE) = Particles(N, d)

# particles(n :: NUTS_result) = particles(n.samples)
using IterTools

function particles(v::Vector{Vector{T} }) where {T}
    map(eachindex(v[1])) do j particles([x[j] for x in v]) end
end

function particles(v::Vector{Array{T,N} }) where {T,N}
    map(CartesianIndices(v[1])) do j particles([x[j] for x in v]) end
end

function particles(s::Vector{NamedTuple{vs, T}}) where {vs, T}
    nt = NamedTuple()
    for k in keys(s[1])
        nt = merge(nt, [k => particles(getproperty.(s,k))])
    end
    nt
end

export parts

using MappedArrays, FillArrays

# Just a little helper function for particles
# https://github.com/baggepinnen/MonteCarloMeasurements.jl/issues/22
parts(d, N::Int=DEFAULT_SAMPLE_SIZE) = particles(d, N)
parts(x::Normal{P} where {P <: AbstractParticles}, N::Int=DEFAULT_SAMPLE_SIZE) = Particles(length(x.μ), Dists.Normal()) * x.σ + x.μ
parts(x::Dists.Sampleable{F,S}, N::Int=DEFAULT_SAMPLE_SIZE) where {F,S} = Particles(N,x)
parts(x::Integer, N::Int=DEFAULT_SAMPLE_SIZE) = parts(float(x))
parts(x::Real, N::Int=DEFAULT_SAMPLE_SIZE) = parts(repeat([x],N))
parts(x::AbstractArray, N::Int=DEFAULT_SAMPLE_SIZE) = Particles(x)
parts(p::Particles, N::Int=DEFAULT_SAMPLE_SIZE) = p

function parts(x::Bernoulli{P} where {P <: AbstractParticles}, N::Int=DEFAULT_SAMPLE_SIZE)
    k = length(x)
    us = Particles(k, Uniform()).particles
    ps = x.p.particles

    return Particles(collect(us .< ps))
end

# function Base.:<(x::AbstractParticles, y::Particles{T,N}) where {T, N}
#     return Particles(x.particles .< y.particles)
# end

# function parts(d::For, N::Int=DEFAULT_SAMPLE_SIZE)
#     ci = CartesianIndices(d.θ)
#     result = similar(Array{Any}, axes(ci))
#     for j in ci
#         result[j] = parts(d.f(Tuple(j)...), N)
#     end
# end

parts(d::ProductMeasure, N::Int=DEFAULT_SAMPLE_SIZE) = parts.(d.data, N)

# size
# dist

# MonteCarloMeasurements.Particles(n::Int, d::For) =
# -(a::Particles{Float64,1000}, b::Array{Float64,1}) = [a-bj for bj in b]

# Base.promote(a::Particles{T,N}, b)  where {T,N} = (a,parts(b))
# Base.promote(a, b::Particles{T,N})  where {T,N} = (parts(a),b)

# promote_rule(::Type{A}, ::Type{B}) where {A <: Real, B <: AbstractParticles{T,N}} where {T} = AbsractParticles{promote_type(A,T),N} where {N}
# promote_rule(::Type{B}, ::Type{A}) where {A <: Real, B <: AbstractParticles{T,N}} where {T} = AbsractParticles{promote_type(A,T),N} where {N}

@inline function particles(m::ConditionalModel, N::Int=DEFAULT_SAMPLE_SIZE)
    return _particles(getmoduletypencoding(m.model), m.model, m.argvals, Val(N))
end

sourceParticles(m::Model, N::Int) = sourceParticles()(m, Val(N))

export sourceParticles
function sourceParticles()

    function(_m::Model, ::Type{Val{_N}}) where {_N}
        proc(_m, st::Assign)  = :($(st.x) = $(st.rhs))
        proc(_m, st::Sample)  = :($(st.x) = parts($(st.rhs), $_N))
        proc(_m, st::Return)  = :(return $(st.rhs))
        proc(_m, st::LineNumber) = nothing

        vals = map(x -> Expr(:(=), x,x),variables(_m))

        wrap(kernel) = @q begin
            $kernel
            $(Expr(:tuple, vals...))
        end

        buildSource(_m, proc, wrap) |> MacroTools.flatten
    end
end

function Base.getindex(a::AbstractArray{P}, i::Particles) where P <: AbstractParticles
    return Particles([a[i.particles[n]][n] for n in eachindex(i.particles)])
end

function Base.getindex(a::AbstractArray{P}, i, j::Particles) where P <: AbstractParticles
    return Particles([a[i,j.particles[n]][n] for n in eachindex(j.particles)])
end

function Base.getindex(a::AbstractArray, i, j::Particles)
    return Particles([a[i,j.particles[n]] for n in eachindex(j.particles)])
end

function Base.getindex(a::AbstractArray{P}, i::Particles, j) where P <: AbstractParticles
    return Particles([a[i.particles[n], j][n] for n in eachindex(i.particles)])
end

function Base.getindex(a::AbstractArray{P}, i::Particles, j::Particles) where P <: AbstractParticles
    return Particles([a[i.particles[n], j.particles[n]][n] for n in eachindex(i.particles)])
end

# function Base.getindex(a::AbstractArray, i::Particles)
#     return Particles([a[ii] for ii in i.particles])
# end

# function Base.getindex(a::AbstractArray, i::Particles, j::Particles)
#     return Particles([a[ii,jj] for (ii,jj) in zip(i.particles, j.particles)])
# end

# function Base.getindex(a::AbstractArray, i::Particles, j)
#     return Particles([a[ii,j] for ii in i.particles])
# end

# Base.to_indices(A, ::Particles) = i.particles

@gg function _particles(M::Type{<:TypeLevel}, _m::Model, _args, _n::Val{_N}) where {_N}
    body = sourceParticles()(type2model(_m), _n) |> loadvals(_args, NamedTuple())
    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end

@gg function _particles(M::Type{<:TypeLevel}, _m::Model, _args::NamedTuple{()}, _n::Val{_N}) where {_N}
    body = sourceParticles()(type2model(_m), _n)
    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end


text:
@testset "Low-Level tests of the PortAudio ringbuf library" begin
    @testset "Can create PA ring buffer" begin
        buf = RingBuffers.PaUtilRingBuffer(8, 16)
        @test RingBuffers.PaUtil_GetRingBufferWriteAvailable(buf) == 16
        @test RingBuffers.PaUtil_GetRingBufferReadAvailable(buf) == 0
    end

    @testset "Throws on non power of two element count" begin
        @test_throws ErrorException RingBuffers.PaUtilRingBuffer(8, 15)
    end

    @testset "Can read/write to PA ring buffer" begin
        buf = RingBuffers.PaUtilRingBuffer(sizeof(Int), 16)
        writedata = collect(1:5)
        readdata = collect(6:10)

        @test RingBuffers.PaUtil_WriteRingBuffer(buf, writedata, 5) == 5
        @test RingBuffers.PaUtil_GetRingBufferWriteAvailable(buf) == 11
        @test RingBuffers.PaUtil_GetRingBufferReadAvailable(buf) == 5
        @test RingBuffers.PaUtil_ReadRingBuffer(buf, readdata, 5) == 5
        @test RingBuffers.PaUtil_GetRingBufferWriteAvailable(buf) == 16
        @test RingBuffers.PaUtil_GetRingBufferReadAvailable(buf) == 0
        @test readdata == writedata
    end

    @testset "Can flush PA ring buffer" begin
        buf = RingBuffers.PaUtilRingBuffer(sizeof(Int), 16)
        writedata = collect(1:5)
        readdata = collect(6:10)

        @test RingBuffers.PaUtil_WriteRingBuffer(buf, writedata, 5) == 5
        @test RingBuffers.PaUtil_GetRingBufferWriteAvailable(buf) == 11
        @test RingBuffers.PaUtil_GetRingBufferReadAvailable(buf) == 5
        RingBuffers.PaUtil_FlushRingBuffer(buf)
        @test RingBuffers.PaUtil_GetRingBufferWriteAvailable(buf) == 16
        @test RingBuffers.PaUtil_GetRingBufferReadAvailable(buf) == 0
    end

    @testset "PA ring buffer handles overflow/underflow" begin
        buf = RingBuffers.PaUtilRingBuffer(sizeof(Int), 8)
        writedata = collect(1:5)
        readdata = collect(6:10)

        @test RingBuffers.PaUtil_WriteRingBuffer(buf, writedata, 5) == 5
        @test RingBuffers.PaUtil_WriteRingBuffer(buf, writedata, 5) == 3
        @test RingBuffers.PaUtil_GetRingBufferWriteAvailable(buf) == 0
        @test RingBuffers.PaUtil_GetRingBufferReadAvailable(buf) == 8
        @test RingBuffers.PaUtil_ReadRingBuffer(buf, readdata, 5) == 5
        @test RingBuffers.PaUtil_ReadRingBuffer(buf, readdata, 5) == 3
        @test RingBuffers.PaUtil_GetRingBufferWriteAvailable(buf) == 8
        @test RingBuffers.PaUtil_GetRingBufferReadAvailable(buf) == 0
    end
end


text:
function readlabelset(io::IO, bytecomp::Bool, bias::Integer)
    n = read(io, Int32)
    dictf = sizehint!(Dict{Float64, String}(), n)
    dictui = sizehint!(Dict{UInt8, String}(), n)
    for i = 1:n
        k = read(io, Float64)
        nchar = read(io, Int8)
        label = String(read(io, nchar))
                     # advance io to next multiple of 8 bytes
        extras = rem(nchar + 1, 8)   # + 1 b/c of reading nchar
        if extras ≠ 0
            skip(io, 8 - extras)
        end
        dictf[k] = label
        if bytecomp && isinteger(k)
            kb = Int(k) + bias
            if 0 < kb < 252
                dictui[kb] = label
            end
        end
    end
    read(io, Int32) == 4 || throw(ArgumentError("Malformed stream"))
    inds = read!(io, Vector{Int32}(read(io, Int32)))
    dictf, dictui, inds
end

function readvaluelabels!(io::IO, vdict::VariableDictionary, bytecomp::Bool, bias::Integer)
    rectyp = Int32(3)
    nms = vdict.nms
    typs = vdict.typs
    vlabs = vdict.valuelabels
    while rectyp == 3
        dictf, dictui, inds = readlabelset(io, bytecomp, bias)
        if bytecomp && length(dictui) == length(dictf)  ## switch to compressed bytes
            for j in inds
                typs[j] = UInt8
                vlabs[nms[j]] = dictui
            end
#=        else
            for j in inds
                vlabs[j] = dictf
            end
=#
        end
        rectyp = read(io, Int32)
    end
    return rectyp
end


text:
using TreeSitter, Test

@testset "TreeSitter" begin
    @testset "Load & Parse" begin
        for lang in keys(TreeSitter.API.LANGUAGES)
            @testset "$lang" begin
                p = Parser(lang)
                tree = parse(p, "")
                @test !isempty(string(tree))
            end
        end
    end
    @testset "Tree Traversal" begin
        p = Parser(:julia)
        tree = parse(p, "[1, 2]")

        out = String[]
        traverse(tree) do node, enter
            enter && push!(out, TreeSitter.node_type(node))
        end
        @test out == ["source_file", "array_expression", "[", "number", ",", "number", "]"]

        out = String[]
        traverse(tree) do node, enter
            enter || push!(out, TreeSitter.node_type(node))
        end
        @test out == ["[", "number", ",", "number", "]", "array_expression", "source_file"]

        out = String[]
        traverse(tree, named_children) do node, enter
            enter && push!(out, TreeSitter.node_type(node))
        end
        @test out == ["source_file", "array_expression", "number", "number"]

        out = String[]
        traverse(tree, named_children) do node, enter
            enter || push!(out, TreeSitter.node_type(node))
        end
        @test out == ["number", "number", "array_expression", "source_file"]
    end
    @testset "Languages" begin
        @testset "bash" begin
            p = Parser(:bash)
            tree = parse(p, "echo {1..10}")
            @test string(tree) == "(program (command name: (command_name (word)) argument: (concatenation (word))))"
        end
        @testset "c" begin
            p = Parser(:c)
            tree = parse(p, "int x = 1;")
            @test string(tree) == "(translation_unit (declaration type: (primitive_type) declarator: (init_declarator declarator: (identifier) value: (number_literal))))"
        end
        @testset "cpp" begin
            p = Parser(:cpp)
            tree = parse(p, "using namespace std;")
            @test string(tree) == "(translation_unit (using_declaration (identifier)))"
        end
        @testset "go" begin
            p = Parser(:go)
            tree = parse(p, "var x := []int{1}")
            @test string(tree) == "(source_file (var_declaration (var_spec name: (identifier) (ERROR) type: (slice_type element: (type_identifier)))) (ERROR (int_literal)))"
        end
        @testset "html" begin
            p = Parser(:html)
            tree = parse(p, "
header
")
            @test string(tree) == "(fragment (element (start_tag (tag_name)) (text) (end_tag (tag_name))))"
        end
        @testset "java" begin
            p = Parser(:java)
            tree = parse(p, "public static void f() {}")
            @test string(tree) == "(program (local_variable_declaration_statement (local_variable_declaration (modifiers) type: (void_type) declarator: (variable_declarator name: (identifier))) (MISSING \";\")) (ERROR (formal_parameters)) (block))"
        end
        @testset "javascript" begin
            p = Parser(:javascript)
            tree = parse(p, "function f() {}")
            @test string(tree) == "(program (function_declaration name: (identifier) parameters: (formal_parameters) body: (statement_block)))"
        end
        @testset "json" begin
            p = Parser(:json)
            tree = parse(p, "{1: 1}")
            @test string(tree) == "(document (object (pair key: (number) value: (number))))"
        end
        @testset "julia" begin
            p = Parser(:julia)
            tree = parse(p, "f(x::Int) = x + 1")
            @test string(tree) == "(source_file (assignment_expression (call_expression (identifier) (argument_list (typed_expression (identifier) (identifier)))) (binary_expression (identifier) (number))))"
        end
        @testset "php" begin
            p = Parser(:php)
            tree = parse(p, "<?php\n\$x = 1;\n?>")
            @test string(tree) == "(program (php_tag) (expression_statement (assignment_expression left: (variable_name (name)) right: (integer))) (text_interpolation))"
        end
        @testset "python" begin
            p = Parser(:python)
            tree = parse(p, "1 < 2 and 2 < 3")
            @test string(tree) == "(module (expression_statement (boolean_operator left: (comparison_operator (integer) (integer)) right: (comparison_operator (integer) (integer)))))"
        end
        @testset "ruby" begin
            p = Parser(:ruby)
            tree = parse(p, "\"Hello\".method(:class).class")
            @test string(tree) == "(program (call receiver: (method_call method: (call receiver: (string) method: (identifier)) arguments: (argument_list (symbol))) method: (identifier)))"
        end
        @testset "rust" begin
            p = Parser(:rust)
            tree = parse(p, "let x: i32 = 13i32;")
            @test string(tree) == "(source_file (let_declaration pattern: (identifier) type: (primitive_type) value: (integer_literal)))"
        end
        @testset "typescript" begin
            p = Parser(:typescript)
            tree = parse(p, "let list: number[] = [1];")
            @test string(tree) == "(program (lexical_declaration (variable_declarator name: (identifier) type: (type_annotation (array_type (predefined_type))) value: (array (number)))))"
        end
    end
    @testset "Captures & Predicates" begin
        p = Parser(:julia)
        source =
        """
        const X = 1
        f(x) = x
        """
        tree = parse(p, source)
        q = query```
        (
            (identifier) @lowercase
            (#match? @lowercase "^[a-z]+$")
        )
        (
            (identifier) @uppercase
            (#match? @uppercase "^[A-Z]+$")
        )
        ```julia
        out = []
        for capture in TreeSitter.each_capture(tree, q, source)
            id = TreeSitter.capture_name(q, capture)
            literal = TreeSitter.slice(source, capture.node)
            push!(out, (id, literal))
        end
        @test out[1] == ("uppercase", "X")
        @test out[2] == ("lowercase", "f")
        @test out[3] == ("lowercase", "x")
        @test out[4] == ("lowercase", "x")
    end
    @testset "Queries" begin
        p = Parser(:c)
        source =
        """
        int main(void) {
            // comment
        }
        """
        tree = parse(p, source)
        q = Query(:c, ["highlights"])
        out = []
        for capture in TreeSitter.each_capture(tree, q, source)
            id = TreeSitter.capture_name(q, capture)
            literal = TreeSitter.slice(source, capture.node)
            push!(out, (id, literal))
        end
        @test out[1] == ("type", "int")
        @test out[2] == ("function", "main")
        @test out[3] == ("variable", "main")
        @test out[4] == ("type", "void")
        @test out[5] == ("comment", "// comment")
    end
end


text:
__precompile__()
module AuditoryFilters

using DSP, Compat
if VERSION >= v"0.7.0-DEV.602"
    import ..DSP: filt, filt!
else
    import Base: filt, filt!
end
import Base: convert
export erb_space, make_erb_filterbank, compute_modulation_cfs, make_modulation_filter,
       modulation_filterbank, gammatonegram, fft2gammatonemx,
       ERBFilterbank, ModulationFilterbank, filt

abstract type Filterbank end

nchannels(f::Filterbank) = length(f.filters)

function filt(fb::Filterbank, x)
    output = zeros(length(x), length(fb.filters))
    @inbounds begin
        for k = 1:length(fb.filters)
            output[:, k] = filt(fb.filters[k], x)
        end
    end
    return output
end

include("ERBFilterbank.jl")
include("ModulationFilterbank.jl")
include("Gammatonegram.jl")

end #module


text:
using ReinforcementLearning
using StableRNGs
using Flux
using IntervalSets

mutable struct KuhnMADDPGHook <: AbstractHook
    eval_freq::Int
    episode_counter::Int
    episode::Vector{Int}
    results::Vector{Float64}
end

function (hook::KuhnMADDPGHook)(::PostEpisodeStage, policy, env)
    hook.episode_counter += 1
    if hook.episode_counter % hook.eval_freq == 0
        push!(hook.episode, hook.episode_counter)
        push!(hook.results, reward(env, 1))
    end
end

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:MADDPG},
    ::Val{:KuhnPoker},
    ::Nothing;
    seed=123,
)
    rng = StableRNG(seed)
    env = KuhnPokerEnv()
    wrapped_env = ActionTransformedEnv(
        StateTransformedEnv(
            env;
            state_mapping = s -> [findfirst(==(s), state_space(env))],
            state_space_mapping = ss -> [[findfirst(==(s), state_space(env))] for s in state_space(env)]
            ),
        # drop the dummy action of the other agent.
        action_mapping = x -> length(x) == 1 ? x : Int(ceil(x[current_player(env)]) + 1),
    )
    ns, na = 1, 1 # dimension of the state and action.
    n_players = 2 # number of players

    init = glorot_uniform(rng)

    create_actor() = Chain(
            Dense(ns, 64, relu; init = init),
            Dense(64, 64, relu; init = init),
            Dense(64, na, tanh; init = init),
        )

    create_critic() = Chain(
        Dense(n_players * ns + n_players * na, 64, relu; init = init),
        Dense(64, 64, relu; init = init),
        Dense(64, 1; init = init),
        )

    policy = DDPGPolicy(
        behavior_actor = NeuralNetworkApproximator(
            model = create_actor(),
            optimizer = ADAM(),
        ),
        behavior_critic = NeuralNetworkApproximator(
            model = create_critic(),
            optimizer = ADAM(),
        ),
        target_actor = NeuralNetworkApproximator(
            model = create_actor(),
            optimizer = ADAM(),
        ),
        target_critic = NeuralNetworkApproximator(
            model = create_critic(),
            optimizer = ADAM(),
        ),
        γ = 0.95f0,
        ρ = 0.99f0,
        na = na,
        start_steps = 1000,
        start_policy = RandomPolicy(-0.99..0.99; rng = rng),
        update_after = 1000,
        act_limit = 0.99,
        act_noise = 0.,
        rng = rng,
    )
    trajectory = CircularArraySARTTrajectory(
        capacity = 100_000, # replay buffer capacity
        state = Vector{Int} => (ns, ),
        action = Float32 => (na, ),
    )

    agents = MADDPGManager(
        Dict((player, Agent(
            policy = NamedPolicy(player, deepcopy(policy)),
            trajectory = deepcopy(trajectory),
        )) for player in players(env) if player != chance_player(env)),
        SARTS, # trace's type
        512, # batch_size
        100, # update_freq
        0, # initial update_step
        rng
    )

    stop_condition = StopAfterEpisode(100_000, is_show_progress=!haskey(ENV, "CI"))
    hook = KuhnMADDPGHook(1000, 0, [], [])
    Experiment(agents, wrapped_env, stop_condition, hook, "# play KuhnPoker with MADDPG")
end

using Plots
ex = E`JuliaRL_MADDPG_KuhnPoker`
run(ex)
scatter(ex.hook.episode, ex.hook.results, xaxis=:log, xlabel="episode", ylabel="reward of player 1")

savefig("assets/JuliaRL_MADDPG_KuhnPoker.png") #hide

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
#=
 TODO find an interface to keep both collision checkers around: Minkowski and Parallel Axis
abstract type CollisionChecker end
struct MinkowskiChecker <: CollisionChecker end
struct ParallelAxisChecker <: CollisionChecker end
is_colliding(vehA, vehB, collision_checker::AbstractCollisionChecker = ParallelAxisChecker())

=#


text:
using ITensors
using ITensorNetworkAD
using ITensorsVisualization
using ITensorNetworkAD.ITensorNetworks:
  Models, inds_network, project_boundary, sqnorm, sqnorm_approx

function peps_tensor(; linkdim, sitedim)
  # left, right, top, bottom, site
  return randn(linkdim, linkdim, linkdim, linkdim, sitedim)
end

N = (3, 4)
ndims = length(N)

site_inds = siteinds("S=1/2", N...)
link_space = 2
inds_net = inds_network(site_inds; linkdims=link_space)

A = peps_tensor(; linkdim=link_space, sitedim=dim(first(site_inds)))
ψ = itensor.((A,), inds_net)

# Project periodic boundary indices
# onto the 1 state
state = 1
ψ = project_boundary(ψ, state)

row = 2
center = (row, :)
cutoff_ = 1e-15
maxdim_ = 100

sqnormψ = sqnorm(ψ)
sqnormψ_approx = sqnorm_approx(ψ; center=center, cutoff=cutoff_, maxdim=maxdim_)
@show noncommoninds(sqnormψ_approx...)
@visualize *(sqnormψ_approx...) contract = false

@disable_warn_order begin
  @show contract(sqnormψ_approx)[] / contract(sqnormψ)[]
end


text:
@testset "prox_l2s!" begin
#test prox_ls2!

m=randn(10)
x=randn(10)
rho=0.0
prox_l2s!(x,rho,m)
@test x==m

rho=1e10
y=deepcopy(x)
prox_l2s!(x,rho,m)
@test isapprox(x,y,rtol=1e-14)

x=Vector{Float64}(undef, 1); x[1]=2.0
m=Vector{Float64}(undef, 1); m[1]=1.0
rho=3.0
prox_l2s!(x,rho,m)
@test x[1]==(7/4)

end


text:
const THREAD_POOL = Base.RefValue{Channel{Int}}()  # defined in __init__(); the thread 1 is reserved for JobScheduler, when nthreads > 2

"""
    set_non_sticky!(t::Task)
    set_non_sticky!(j::Job)

Allow a job run in different threads. Not available in Julia 1.1.
"""
function set_non_sticky!(t::Task)
    @static if :sticky in fieldnames(Task)  # sticky controls wether a job is running in different threads, not found in Julia 1.1, found in julia 1.2 and so on.
        if nthreads() > 1
            t.sticky = false
        end
    end
end

set_non_sticky!(j::Job) = set_non_sticky!(j.task)

function set_thread_id!(t::Task, tid::Int)
    @static if :sticky in fieldnames(Task)  # sticky controls wether a job is running in different threads, not found in Julia 1.1, found in julia 1.2 and so on.
        if nthreads() > 1
            t.sticky = false
            ccall(:jl_set_task_tid, Cvoid, (Any, Cint), t, tid-1)
        end
    end
end

function set_thread_id!(j::Job, tid::Int)
    @static if :sticky in fieldnames(Task)  # sticky controls wether a job is running in different threads, not found in Julia 1.1, found in julia 1.2 and so on.
        if nthreads() > 1
            j.task.sticky = false
            j._thread_id = tid
            ccall(:jl_set_task_tid, Cvoid, (Any, Cint), j.task, tid-1)
        end
    end
end

function get_thread_id()

end

function schedule_thread(j::Job)
    @static if :sticky in fieldnames(Task)
        if nthreads() > 1
            j.task.sticky = false
            # take the next free thread... Will block/wait until a thread becomes free
            j._thread_id = take!(THREAD_POOL[])
            ccall(:jl_set_task_tid, Cvoid, (Any, Cint), j.task, j._thread_id-1)
        end
    end
    schedule(j.task)
end

"""
    free_thread(j::Job)

Make thread available again after work is done!
"""
function free_thread(j::Job)
    j._thread_id == 0 && return # do nothing. 0 means (1) nthreads = 1 (2) job is not started or the thread was freed.
    put!(THREAD_POOL[], j._thread_id)
    j._thread_id = 0
    return
end


text:
##### Chinese Restaurant Process arrival time utilities

type CRPinterarrival <: DiscreteUnivariateDistribution
  theta::Float64
  alpha::Float64
  n::Int
  k::Int
  # crp::Bool
end

Base.minimum(s::CRPinterarrival) = 1
Base.maximum(s::CRPinterarrival) = Inf

function partial(f,a...)
  (b...) -> f(a...,b...)
end

function CRP(a...)
  """
  Utility function for passing arguments `n` and `k` during sampling updates
  """
  partial(CRPinterarrival,a...)
end

Distributions.logpdf(s::CRPinterarrival,x::Int) = _logpdf(s,x)
Distributions.logpdf(s::CRPinterarrival,v::Vector{Int}) = _logpdf_batch(s,v)
function _logpdf(s::CRPinterarrival,x::Int)
  return crp_logpdf(s.theta,s.alpha,s.n,s.k,x)
end

function crp_logpdf(theta::Float64,alpha::Float64,n::Int,k::Int,x::Int)
  ka = k*alpha
  nt = n + theta
  logp = log(theta + ka) - log(nt)
  if x > 1
    nka = n - ka
    for i in 2:x
      logp += log(nka + i - 2) - log(nt + i - 1)
    end
  end
  return logp
end

function crp_logpdf(theta::Float64,alpha::Float64,n::Int,k::Int,v::Vector{Int})
  # for a vector of evaluation points
  vmax = maximum(v)
  vmin = minimum(v)
  nv = size(v,1)
  idx_all = 1:nv
  idx_j = trues(nv)
  ret = zeros(Float64,nv)

  ka = k*alpha
  nt = n + theta
  nka = n - ka
  logp = log(theta + ka) - log(nt)
  for j in 1:vmax
    j > 1 ? logp += log(nka + j - 2) - log(nt + j - 1) : nothing
    for i in idx_all[idx_j]
      if v[i]==j
        ret[i] = logp
        idx_j[i] = false
      end
    end
  end
  return ret
end

function _logpdf_batch(s::CRPinterarrival,v::Vector{Int})
  return crp_logpdf(s.theta,s.alpha,s.n,s.k,v)
end

Distributions.pdf(s::CRPinterarrival,x::Int) = _pdf(s,x)
function _pdf(s::CRPinterarrival,x::Int)
  return exp.(logpdf(s,x))
end

Distributions.cdf(s::CRPinterarrival,x::Int64) = _cdf(s,x)
function _cdf(s::CRPinterarrival,x::Int64)
  return crp_cdf(s.theta,s.alpha,s.n,s.k,x)
end

function crp_cdf(theta::Float64,alpha::Float64,n::Int,k::Int,x::Int)
  if x==0
    return 0.
  else
    ka = k*alpha
    nt = n + theta
    P1 = (theta + ka)/(nt)
    if x > 1
      nka = n - ka
      prev = P1
      P = 0
      for j in 2:x
        run = exp( log(prev) + log(nka + j - 2) - log(nt + j - 1) )
        P += run
        prev = run
      end
    end
    return x > 1 ? P1 + P : P1
  end
end

Base.Random.rand(s::CRPinterarrival) = _rand(s)
function _rand(s::CRPinterarrival)
  coin = 0
  ct = 0
  while coin != 1
    ct += 1
    p = (s.theta + s.alpha*s.k)/(s.theta + s.n + ct - 1)
    coin = rand(Bernoulli(p))
  end
  return ct
end

# slice sampling utilities
function crp_theta_logpdf(theta::Float64,alpha::Float64,k::Int,n::Int,log_prior::Function)
  """
  calculate unnormalized log-pdf proportional to `theta` in the CRP

  log_prior is a function that returns the (possibly unnormalized) prior log-probability
    of `theta`
  """
  if theta <= -alpha
    # println("theta: ",theta,"  //  alpha: ",alpha)
    return -Inf
  end
  logp = log_prior(theta,alpha)
  for j in 1:(k-1)
    logp += log(theta + j*alpha)
  end

  for m in 1:(n-1)
    logp += -log(theta + m)
  end
  return logp
end

function crp_theta_trans_logpdf(theta_trans::Float64,alpha::Float64,k::Int,n::Int,log_prior::Function)
  """
  computes log-pdf when theta has been transformed to the entire real line

  theta_trans = log(theta + alpha) (for fixed alpha)
  """
  theta = exp(theta_trans) - alpha
  return crp_theta_logpdf(theta,alpha,k,n,log_prior) + theta_trans
end

function crp_alpha_logpdf(alpha::Float64,theta::Float64,T::Vector{Int},n::Int,log_prior::Function)
  """
  calculate unnormalized log-pdf proportional to `alpha` in the CRP

  log_prior is a function that returns the (possibly unnormalized) prior log-probability
    of `alpha`
  """
  if theta <= -alpha || alpha > 1. || alpha < 0.
    return -Inf
  end
  k = size(T,1)
  logp = log_prior(alpha) + sum(log.(theta + alpha.*(1:(k-1))))
  logp += sum( lgamma.(T[2:end] .- 1 .- alpha.*(1:(k-1))) .- lgamma.(T[1:(k-1)] .- alpha.*(1:(k-1))) )
  return logp
end

function crp_alpha_trans_logpdf(alpha_trans::Float64,theta::Float64,T::Vector{Int},n::Int,log_prior::Function)
  """
  computes log-pdf when alpha has been transformed to the entire real line

  alpha_trans = log(alpha - max(0,-theta)) - log(1-alpha) for fixed theta
  """
  alpha = (exp(alpha_trans) + max(0,-theta))/(1 + exp(alpha_trans))
  return crp_alpha_logpdf(alpha,theta,T,n,log_prior) + alpha_trans - log(1 + exp(alpha_trans)) + log(1 - alpha)
end

function crp_alpha_trans_coupled_logpdf(alpha_trans::Float64,theta::Float64,PP::Vector{Int},T::Vector{Int},n::Int,log_prior::Function)
  """
  for the "coupled CRP" model, where the NTL alpha parameter is the same as the CRP alpha paramter
  computes log-pdf when alpha has been transformed to the entire real line

  alpha_trans = log(alpha - max(0,-theta)) - log(1-alpha) for fixed theta
  """
  alpha = (exp(alpha_trans) + max(0,-theta))/(1 + exp(alpha_trans))
  crp_arrivals_lpdf = crp_alpha_logpdf(alpha,theta,T,n,log_prior) # arrivals contribution
  trans_lpdf =  alpha_trans - log(1 + exp(alpha_trans)) + log(1 - alpha) # volume conrrection of transformation
  dummy_prior = x -> 0.
  crp_sbr_lpdf = ntl_alpha_logpdf(alpha,PP,T,dummy_prior) # size-biased reinforcement contribution

  return crp_arrivals_lpdf + crp_sbr_lpdf + trans_lpdf
end

# parameter updates
# uncoupled CRP
function update_crp_interarrival_params!(ia_params::Vector{Float64},PP::Vector{Int},T::Vector{Int},
  n::Int,log_prior_theta::Function,log_prior_alpha::Function,
  w_t::Float64,w_a::Float64,coupled::Bool)

  # update theta via slice sampling
  K = size(T,1)
  ss_gt = x -> crp_theta_trans_logpdf(x,ia_params[2],K,n,log_prior_theta)
  theta_trans = log(ia_params[1] + ia_params[2])
  theta_trans_ss = slice_sampling(ss_gt,w_t,theta_trans)
  ia_params[1] = exp(theta_trans_ss) - ia_params[2]

  if !coupled
    ss_ga = x -> crp_alpha_trans_logpdf(x,ia_params[1],T,n,log_prior_alpha)
  else
    ss_ga = x -> crp_alpha_trans_coupled_logpdf(x,ia_params[1],PP,T,n,log_prior_alpha)
  end
  alpha_trans = log(ia_params[2] - max(0.,-ia_params[1])) - log(1 - ia_params[2])
  alpha_trans_ss = slice_sampling(ss_ga,w_a,alpha_trans)
  ia_params[2] = (exp(alpha_trans_ss) + max(0,-ia_params[1]))/(1 + exp(alpha_trans_ss))

end

function crp_theta_loglik(theta::Float64,alpha::Float64,K::Int,n::Int)
  return sum( log.(theta .+ (1:(K-1)).*alpha) ) - sum( log.(theta .+ (1:(n-1))) )
end

function grad_crp_theta_loglik(theta::Float64,alpha::Float64,K::Int,n::Int)
  return sum( 1./(theta .+ (1:(K-1)).*alpha) ) - sum( 1./(theta .+ (1:(n-1))) )
end

function hess_crp_theta_loglik(theta::Float64,alpha::Float64,K::Int,n::Int)
  return -sum( 1./(theta .+ (1:(K-1)).*alpha).^2 ) + sum( 1./(theta .+ (1:(n-1))).^2 )
end

function initialize_crp_params(theta_prior::UnivariateDistribution,alpha_prior::UnivariateDistribution)
  alpha = mean(alpha_prior)
  theta = mean(theta_prior)
  return [theta; alpha]

end


text:
function week_test_functional()
    test_courier_bike_1()
    test_courier_car_1()
end

function test_courier_bike_1()
    info("STARTING COURIER BIKE 1")
    ok, weekly_schedule = schedule(test_bike_week_1_employees, test_bike_week_1_env)
    @test ok == true
    @test week_sum_coverage(test_bike_week_1_env) * TEST_THRESHOLD > week_hours_scheduled(weekly_schedule)
    info("ENDING COURIER BIKE 1")

end

function test_courier_car_1()
    info("STARTING COURIER CAR 1")
    ok, weekly_schedule = schedule(test_courier_car_1_employees, test_courier_car_1_env)
    @test ok == true
    @test week_sum_coverage(test_courier_car_1_env) * TEST_THRESHOLD > week_hours_scheduled(weekly_schedule)
    info("ENDING COURIER CAR 1")
end

function test_courier_incorrect_availability()
    info("STARTING COURIER INCORRECT AVAILABILITY LENGTH")
    ok, _ = schedule(test_bike_incorrect_availability_length_employees, test_bike_incorrect_availability_length_env)
    @test ok != true
    info("ENDING COURIER INCORRECT AVAILABILITY LENGTH")
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
module ECR

using AWS
using Base64
using Mocking

# Improper casing to avoid issues with Module name and AWS.AWSService
@service Ecr

function get_authorization_token(config::AWSConfig, params::AbstractDict=Dict{String,Any}())
    dep_msg = """
        `get_authorization_token(config::AWSConfig, params::AbstractDict=Dict{String,Any}())` is deprecated and will be removed.
        Use the AWS @service ECR.get_authorization_token() functionality instead.
        """
    Base.depwarn(dep_msg, :get_authorization_token)
    return Ecr.get_authorization_token(params; aws_config=config)
end

export get_login

"""
    get_login(registry_ids::Union{AbstractString, Integer}="") -> Cmd

Gets the AWS ECR authorization token and returns the corresponding docker login command.
The AWS `config` keyword parameter is optional (will use the default if it's not passed in).
"""
get_login

function get_login(registry_id::AbstractString=""; config::AWSConfig=global_aws_config())
    # Note: Although `get_authorization_token` can take multiple registry IDs at once it
    # will only return a "proxyEndpoint" for the first registry. Additionally, the
    # `aws ecr get-login` command returns a `docker login` command for each registry ID
    # passed in. Because of these factors we'll do our processing on a single registry.

    response = if !isempty(registry_id)
        @mock get_authorization_token(config, Dict("registryIds" => [registry_id]))
    else
        @mock get_authorization_token(config)
    end

    authorization_data = first(response["authorizationData"])
    token = String(base64decode(authorization_data["authorizationToken"]))
    username, password = split(token, ':')
    endpoint = authorization_data["proxyEndpoint"]

    return `docker login -u $username -p $password $endpoint`
end

function get_login(registry_id::Integer; config::AWSConfig=global_aws_config())
    return get_login(lpad(registry_id, 12, '0'); config=config)
end

end  # ECR