text:
using Base64
using UUIDs
using Random

const snapsvgjs = joinpath(@__DIR__, "..", "deps", "snap.svg-min.js")

# Packages can insert extra XML namespaces here to be defined in the output
# SVG.
const xmlns = Dict()

# All svg (in our use) coordinates are in millimeters. This number gives the
# largest deviation from the true position allowed in millimeters.
const eps = 0.01

# Format a floating point number into a decimal string of reasonable precision.
function svg_fmt_float(x::Fractional)
    a = @sprintf("%0.8f", round(x / eps) * eps)
    n = length(a)

    while a[n] == '0'
        n -= 1
    end

    if a[n] == '.'
        n -= 1
    end

    a[1:n]
end

# A much faster version of svg_fmt_float. This does not allocate any
# temporary buffers, because it writes directly to the output.
function svg_print_float(io::IO, x::AbstractFloat)
    ndig = 2

    if isfinite(x)
        if x < 0
            write(io, '-')
            x = abs(x)
        end
        x = round(x/eps)*eps
        xt = trunc(UInt, x)
        dx = x - convert(Float64, xt)
        0 <= dx < 1 || error("Formatting overflow")
        svg_print_uint(io, xt, 1)  # width=1 prints 0.2 instead of .2
        dxi = round(UInt, dx/eps)
        if dxi != 0
            write(io, '.')
            svg_print_uint(io, dxi, ndig, true)
        end
    elseif isnan(x)
        write(io, "NaN")
    elseif x == Inf
        write(io, "Inf")
    elseif x == -Inf
        write(io, "-Inf")
    end
end

let a = Array{UInt8}(undef, 20)
    global svg_print_uint
    function svg_print_uint(io::IO, x::Unsigned, width = 0, drop = false)
        n = length(a)
        while x > 0 && n > 0
            x, r = divrem(x, 10)
            a[n] = r
            n -= 1
        end
        n == 0 && error("Formatting overflow")
        for i = 1:width-(length(a)-n)
            write(io, '0')
        end
        last = length(a)
        if drop
            while last > n && a[last] == 0
                last -= 1
            end
        end
        for i = n+1:last
            write(io, '0'+a[i])
        end
    end
end

# Format a color for SVG.
svg_fmt_color(c::Color) = string("#", hex(c))
svg_fmt_color(c::Nothing) = "none"

# Replace newlines in a string with the appropriate SVG tspan tags.
function svg_newlines(input::AbstractString)
    output = IOBuffer()
    inputs = split(input,"\n")
    write(output, inputs[1])
    for mat in inputs[2:end]
        write(output, """<tspan x="0" dy="1.2em">""", mat, "</tspan>")
    end
    return String(take!(output))
end

# Javascript in a <script> tag in SVG needs to escape '"' and '<'.
#=function escape_script(js::AbstractString)=#
    #=return replace(replace(js, "&", "&"), "<", "<")=#
#=end=#

# Javascript contained to CDATA block needs to avoid ']]'
escape_script(js::AbstractString) = replace(js, "]]"=>"] ]")

# When subtree rooted at a context is drawn, it pushes its property children
# in the form of a property frame.
mutable struct SVGPropertyFrame
    # Vector properties in this frame.
    vector_properties::Dict{Type, Property}

    # True if this property frame has scalar properties. Scalar properties are
    # emitted as a group (<g> tag) that must be closed when the frame is popped.
    has_scalar_properties::Bool

    # True if this property frame includes a link (<a> tag) that needs
    # to be closed.
    has_link::Bool

    # True if the property frame include a mask (<mask> tag) than needs to be
    # closed.
    has_mask::Bool

    # True if the property frame has a scalar Clip property, which needs some
    # special handling.
    has_scalar_clip::Bool
end

SVGPropertyFrame() = SVGPropertyFrame(Dict{Type, Property}(), false, false, false, false)

const LOCAL_RNG = MersenneTwister()

mutable struct SVG <: Backend
    # Image size in millimeters.
    width::AbsoluteLength
    height::AbsoluteLength

    # Output stream.
    out::IO

    # Save output from IOBuffers to allow multiple calls to writemime
    cached_out::Union{AbstractString, Nothing}

    # Unique ID for the figure.
    id::AbstractString

    # Current level of indentation.
    indentation::Int

    # Stack of property frames (groups of properties) currently in effect.
    property_stack::Vector{SVGPropertyFrame}

    # SVG forbids defining the same property twice, so we have to keep track
    # of which vector property of which type is in effect. If two properties of
    # the same type are in effect, the one higher on the stack takes precedence.
    vector_properties::Dict{Type, Union{Property, Nothing}}

    # Clip-paths that need to be defined at the end of the document.
    clippaths::OrderedDict{ClipPrimitive, String}

    # Batched forms to be included within <def> tags.
    batches::Vector{Tuple{FormPrimitive, String}}

    # Embedded objects included immediately before the </svg> tag, such as extra
    # javascript or css.
    embobj::Set{AbstractString}

    # True when finish has been called and no more drawing should occur
    finished::Bool

    # Backend is responsible for opening/closing the file
    ownedfile::Bool

    # Filename when ownedfile is true
    filename::Union{AbstractString, Nothing}

    # Emit the graphic on finish when writing to a buffer.
    emit_on_finish::Bool

    # IDs of the SVG element currently being generated. `has_current_id` is
    # false if the element being drawn does not have an id.
    current_id::AbstractString
    has_current_id::Bool

    # A counter used to generate unique IDs
    id_count::Int

    # Filenames of javsacript to include before any JSCall code.
    jsheader::Set{AbstractString}

    # (Name, binding) pairs of javascript modules the embedded code depends on
    jsmodules::Set{Tuple{AbstractString, AbstractString}}

    # User javascript from JSCall attributes
    scripts::Vector{AbstractString}

    # Use javascript extensions to add interactivity, etc.
    withjs::Bool

    # What to do with javascript. One of:
    #    none: generate a static SVG without any javascript
    #    exclude: exclude external javascript libraries
    #    embed: embed external libraries
    #    linkabs: link to external libraries (absolute path)
    #    linkrel: link to external libraries (relative path)
    #
    jsmode::Symbol

    panelcoords::Tuple
end

function SVG(out::IO,
             width::AbsoluteLength,
             height::AbsoluteLength,
             emit_on_finish::Bool=true,
             jsmode::Symbol=:none;

             cached_out = nothing,
             id = string("img-", string(uuid4(LOCAL_RNG))[1:8]),
             indentation = 0,
             property_stack = Array{SVGPropertyFrame}(undef, 0),
             vector_properties = Dict{Type, Union{Property, Nothing}}(),
             clippaths = OrderedDict{ClipPrimitive, String}(),
             batches = Array{Tuple{FormPrimitive, String}}(undef, 0),
             embobj = Set{AbstractString}(),
             finished = false,
             ownedfile = false,
             filename = nothing,
             current_id = "",
             has_current_id = false,
             id_count = 0,
             jsheader = Set{AbstractString}(),
             jsmodules = Set{Tuple{AbstractString, AbstractString}}((("Snap.svg", "Snap"),)),
             scripts = AbstractString[],
             withjs = jsmode != :none,
             panelcoords = ())

    img = SVG(width,
              height,
              out,
              cached_out,
              id,
              indentation,
              property_stack,
              vector_properties,
              clippaths,
              batches,
              embobj,
              finished,
              ownedfile,
              filename,
              emit_on_finish,
              current_id,
              has_current_id,
              id_count,
              jsheader,
              jsmodules,
              scripts,
              withjs,
              jsmode,
              panelcoords)

    writeheader(img)
    return img
end

# Write to a file.
SVG(filename::AbstractString, width=default_graphic_width, height=default_graphic_height,
        jsmode::Symbol=:none) =
        SVG(open(filename, "w"), width, height, true, jsmode; ownedfile=true, filename=filename)

# Write to buffer.
SVG(width::MeasureOrNumber=default_graphic_width, height::MeasureOrNumber=default_graphic_height,
        emit_on_finish::Bool=true, jsmode::Symbol=:none) =
        SVG(IOBuffer(), width, height, emit_on_finish, jsmode)

canbatch(img::SVG) = true
iswithjs(img::SVG) = img.withjs

function (img::SVG)(x)
    draw(img, x)
end

# Return the next unique element ID. Sort of like gensym for SVG elements.
function genid(img::SVG)
    img.id_count += 1
    return @sprintf("%s-%d", img.id, img.id_count)
end

# Constructors that turn javascript extensions on
SVGJS(out::IO, width=default_graphic_width, height=default_graphic_height,
        emit_on_finish::Bool=true; jsmode::Symbol=:embed) =
        SVG(out, width, height, emit_on_finish, jsmode)
SVGJS(filename::AbstractString, width=default_graphic_width, height=default_graphic_height;
        jsmode::Symbol=:embed) = SVG(filename, width, height, jsmode)
SVGJS(width::MeasureOrNumber=default_graphic_width, height::MeasureOrNumber=default_graphic_height,
        emit_on_finish::Bool=true, jsmode::Symbol=:embed) =
        SVG(width, height, emit_on_finish, jsmode)

function writeheader(img::SVG)
    widthstr = svg_fmt_float(img.width.value)
    heightstr = svg_fmt_float(img.height.value)
    write(img.out,
          """
          <?xml version="1.0" encoding="UTF-8"?>
          <svg xmlns="http://www.w3.org/2000/svg"
               xmlns:xlink="http://www.w3.org/1999/xlink"
          """)

    for (ns, uri) in xmlns
        write(img.out, """     xmlns:$(ns)="$(uri)"\n""")
    end

    write(img.out,
          """
               version="1.2"
               width="$(widthstr)mm" height="$(heightstr)mm" viewBox="0 0 $(widthstr) $(heightstr)"
               stroke="$(svg_fmt_color(default_stroke_color))"
               fill="$(svg_fmt_color(default_fill_color))"
               stroke-width="$(svg_fmt_float(default_line_width.value))"
               font-size="$(svg_fmt_float(default_font_size.value))"
          """)
    img.withjs && write(img.out, "\n     id=\"$(img.id)\"")
    write(img.out, ">\n")
    write(img.out,
          """
          <defs>
            <marker id="arrow" markerWidth="15" markerHeight="7" refX="5" refY="3.5" orient="auto" markerUnits="strokeWidth">
              <path d="M0,0 L15,3.5 L0,7 z" stroke="context-stroke" fill="context-stroke"/>
            </marker>
          </defs>
          """)
#    write(img.out, "\n")
          return img
end

function register_coords(backend::SVG, box, units, transform, property::SVGClass)
    if length(property.primitives) == 1 && property.primitives[1].value == "plotpanel"
        backend.panelcoords = (box, units, transform)
    end
    nothing
end

function reset(img::SVG)
    if img.ownedfile
        img.out = open(img.filename, "w")
    else
        try
            seekstart(img.out)
        catch
            error("Backend can't be reused, since the output stream is not seekable.")
        end
    end
    writeheader(img)
    img.finished = false
end

function finish(img::SVG)
    img.finished && return

    while !isempty(img.property_stack)
        pop_property_frame(img)
    end

    for obj in img.embobj
        write(img.out, obj)
        write(img.out, "\n")
    end

    # defs
    if !isempty(img.clippaths) | !isempty(img.batches)
        write(img.out, "<defs>\n")
        img.indentation += 1
    end

    for (clippath, id) in img.clippaths
        indent(img)
        img.indentation += 1
        write(img.out, "<clipPath id=\"$(id)\">\n")
        indent(img)
        write(img.out, "<path d=\"")
        print_svg_path(img.out, clippath.points)
        write(img.out, "\" />\n")
        img.indentation -= 1
        indent(img)
        write(img.out, "</clipPath>\n")
    end

    for (primitive, id) in img.batches
        indent(img)
        print(img.out, "<g id=\"", id, "\">\n")
        img.indentation += 1
        draw(img, primitive, 1)
        img.indentation -= 1
        indent(img)
        print(img.out, "</g>\n")
    end

    if !isempty(img.clippaths) | !isempty(img.batches)
        img.indentation -= 1
        write(img.out, "</defs>\n")
    end

    if img.withjs
        if img.jsmode == :embed
            write(img.out,
                """
                <script> <![CDATA[
                $(escape_script(read(snapsvgjs, String)))
                ]]> </script>
                """)
        elseif img.jsmode == :linkabs
            write(img.out,
                """
                <script xlink:href="$(snapsvgjs)"></script>
                """)
        elseif img.jsmode == :linkrel
            write(img.out,
                """
                <script xlink:href="$(basename(snapsvgjs))"></script>
                """)
        end

        if !isempty(img.scripts) || !isempty(img.jsheader)
            if img.jsmode == :embed
                write(img.out, "<script> <![CDATA[\n")
                for script in img.jsheader
                    write(img.out, escape_script(read(script, String)), "\n")
                end
            elseif img.jsmode == :linkabs
                for script in img.jsheader
                    write(img.out,
                        """
                        <script xlink:href="$(script)"></script>
                        """)
                end
                write(img.out, "<script> <![CDATA[\n")
            elseif img.jsmode == :linkrel
                for script in img.jsheader
                    write(img.out,
                        """
                        <script xlink:href="$(basename(script))"></script>
                        """)
                end
                write(img.out, "<script> <![CDATA[\n")
            end

            mod_names = join([string("\"", name, "\"") for (name, binding) in img.jsmodules], ", ")
            mod_bindings = join([binding for (name, binding) in img.jsmodules], ", ")
            glob_mod_bindings = join([string("glob.", binding)
                                      for (name, binding) in img.jsmodules], ", ")

            write(img.out,
                """
                (function (glob, factory) {
                    // AMD support
                      if (typeof require === "function" && typeof define === "function" && define.amd) {
                        require([$(mod_names)], function ($(mod_bindings)) {
                            factory($(mod_bindings));
                        });
                      } else {
                          factory($(glob_mod_bindings));
                      }
                })(window, function ($(mod_bindings)) {
                    var fig = Snap(\"#$(img.id)\");
                """)

            for script in img.scripts
                write(img.out, escape_script(script), "\n")
            end

            write(img.out,
                """
                    });
                """)

            write(img.out, "]]> </script>\n")
        end

    end

    write(img.out, "</svg>\n")
    hasmethod(flush, (typeof(img.out),)) && flush(img.out)

    img.ownedfile && close(img.out)

    img.finished = true

    # If we are writing to a buffer. Collect the string and emit it.
    img.emit_on_finish && typeof(img.out) == IOBuffer && display(img)
end

isfinished(img::SVG) = img.finished

function show(io::IO, ::MIME"text/html", img::SVG)
    if img.cached_out === nothing
        img.cached_out = String(take!(img.out))
    end
    write(io, img.cached_out)
end

function show(io::IO, ::MIME"image/svg+xml", img::SVG)
    if img.cached_out === nothing
        img.cached_out = String(take!(img.out))
    end
    write(io, img.cached_out)
end

root_box(img::SVG) = BoundingBox(0mm, 0mm, img.width, img.height)

function indent(img::SVG)
    for i in 1:img.indentation
        write(img.out, "  ")
    end
end

# Draw

# Generate SVG path data from an array of points.
#
# Args:
#   out: Output stream.
#   points: points on the path
#   bridge_gaps: when true, remove non-finite values, rather than forming
#                separate lines.
#
# Returns:
#   A string containing SVG path data.
#
function print_svg_path(out, points::Vector{AbsoluteVec2},
                        bridge_gaps::Bool=false)
    isfirst = true
    for point in points
        x, y = point[1].value, point[2].value
        if !(isfinite(x) && isfinite(y))
            isfirst = true
            continue
        end

        if isfirst
            isfirst = false
            write(out, 'M')
            svg_print_float(out, x)
            write(out, ',')
            svg_print_float(out, y)
            write(out, " L")
        else
            write(out, ' ')
            svg_print_float(out, x)
            write(out, ' ')
            svg_print_float(out, y)
        end
    end
end

# TODO: This logic should be moved to Gadfly
# Return array of paths to draw with printpath
# array is formed by splitting by NaN values
#function make_paths(points::Vector{AbsoluteVec})
    #paths = Vector{AbsoluteVec}[]
    #nans = find(xy -> isnan(xy[1]) || isnan(xy[2]),
                #[(point[1].value, point[2].value) for point in points])
    #if length(nans) == 0
        #push!(paths, points)
    #else
        #nans = [0, nans..., length(points) + 1]
        #i, n = 1, length(nans)
        #while i <= n-1
            #if nans[i] + 1 < nans[i + 1]
                #push!(paths, points[(nans[i]+1):(nans[i+1] - 1)])
            #end
            #i += 1
        #end
    #end
    #paths
#end

# Property Printing
# -----------------

function print_property(img::SVG, property::StrokePrimitive)
    if property.color.alpha != 1.0
        @printf(img.out, " stroke=\"%s\" stroke-opacity=\"%0.3f\"",
                svg_fmt_color(color(property.color)), property.color.alpha)
    else
        @printf(img.out, " stroke=\"%s\"", svg_fmt_color(color(property.color)))
    end
end

function print_property(img::SVG, property::FillPrimitive)
    if property.color.alpha != 1.0
        @printf(img.out, " fill=\"%s\" fill-opacity=\"%0.3f\"",
                svg_fmt_color(color(property.color)), property.color.alpha)
    else
        @printf(img.out, " fill=\"%s\"", svg_fmt_color(color(property.color)))
    end
end

function print_property(img::SVG, property::StrokeDashPrimitive)
    if isempty(property.value)
        print(img.out, " stroke-dasharray=\"none\"")
    else
        print(img.out, " stroke-dasharray=\"")
        svg_print_float(img.out, property.value[1].value)
        for i in 2:length(property.value)
            print(img.out, ',')
            svg_print_float(img.out, property.value[i].value)
        end
        print(img.out, '"')
    end
end

# Format a line-cap specifier into the attribute string that SVG expects.
svg_fmt_linecap(::LineCapButt) = "butt"
svg_fmt_linecap(::LineCapSquare) = "square"
svg_fmt_linecap(::LineCapRound) = "round"

print_property(img::SVG, property::StrokeLineCapPrimitive) =
        @printf(img.out, " stroke-linecap=\"%s\"", svg_fmt_linecap(property.value))

# Format a line-join specifier into the attribute string that SVG expects.
svg_fmt_linejoin(::LineJoinMiter) = "miter"
svg_fmt_linejoin(::LineJoinRound) = "round"
svg_fmt_linejoin(::LineJoinBevel) = "bevel"

print_property(img::SVG, property::StrokeLineJoinPrimitive) =
        @printf(img.out, " stroke-linejoin=\"%s\"", svg_fmt_linejoin(property.value))

function print_property(img::SVG, property::LineWidthPrimitive)
    print(img.out, " stroke-width=\"")
    svg_print_float(img.out, property.value.value)
    print(img.out, '"')
end

function print_property(img::SVG, property::FillOpacityPrimitive)
    print(img.out, " fill-opacity=\"")
    svg_print_float(img.out, property.value)
    print(img.out, '"')
end

function print_property(img::SVG, property::StrokeOpacityPrimitive)
    print(img.out, " stroke-opacity=\"")
    svg_print_float(img.out, property.value)
    print(img.out, '"')
end

print_property(img::SVG, property::VisiblePrimitive) =
        @printf(img.out, " visibility=\"%s\"", property.value ? "visible" : "hidden")

# I may end up applying the same clip path to many forms separately, so I
# shouldn't make a new one for each applicaiton. Where should that happen?
function print_property(img::SVG, property::ClipPrimitive)
    url = clippathurl(img, property)
    @printf(img.out, " clip-path=\"url(#%s)\"", url)
end

print_property(img::SVG, property::FontPrimitive) =
        @printf(img.out, " font-family=\"%s\"", escape_string(property.family))

function print_property(img::SVG, property::FontSizePrimitive)
    print(img.out, " font-size=\"")
    svg_print_float(img.out, property.value.value)
    print(img.out, '"')
end

print_property(img::SVG, property::SVGIDPrimitive) =
        @printf(img.out, " id=\"%s\"", escape_string(property.value))
print_property(img::SVG, property::SVGClassPrimitive) =
        @printf(img.out, " class=\"%s\"", escape_string(property.value))
print_property(img::SVG, property::SVGAttributePrimitive) =
        @printf(img.out, " %s=\"%s\"", property.attribute, escape_string(property.value))

function print_property(img::SVG, property::JSIncludePrimitive)
    push!(img.jsheader, property.value)
    if property.jsmodule != nothing
        push!(img.jsmodules, property.jsmodule)
    end
end

function print_property(img::SVG, property::JSCallPrimitive)
    @assert img.has_current_id
    push!(img.scripts, @sprintf("fig.select(\"#%s\")\n   .%s;", img.current_id, property.code))
end

# Print the property at the given index in each vector property
function print_vector_properties(img::SVG, idx::Int, suppress_fill::Bool=false)
    if haskey(img.vector_properties, JSCall)
        if haskey(img.vector_properties, SVGID)
            img.current_id = img.vector_properties[SVGID].primitives[idx].value
        else
            img.current_id = genid(img)
            print_property(img, SVGIDPrimitive(img.current_id))
        end
        img.has_current_id = true
    end

    has_stroke_opacity = haskey(img.vector_properties, StrokeOpacity)
    has_fill_opacity = haskey(img.vector_properties, FillOpacity)

    for (propertytype, property) in img.vector_properties
        if property === nothing ||
           (propertytype == Fill && suppress_fill)
            continue
        end

        idx > length(property.primitives) &&
                error("Vector form and vector property differ in length. Can't distribute.")

        # let the opacity primitives clobber the alpha value in fill and stroke
        if propertytype == Fill && has_fill_opacity
           print_property(img, FillPrimitive(RGBA{Float64}(color(property.primitives[idx].color), 1.0)))
        elseif propertytype == Stroke && has_stroke_opacity
           print_property(img, StrokePrimitive(RGBA{Float64}(color(property.primitives[idx].color), 1.0)))
        else
            print_property(img, property.primitives[idx])
        end
    end

    img.has_current_id = false
end

function print_property(img::SVG, property::ArrowPrimitive)
    print(img.out, " marker-end=\"url(#arrow)\"")
end

# Form Drawing
# ------------

function draw(img::SVG, form::Form{T}) where T
    for i in 1:length(form.primitives)
        draw(img, form.primitives[i], i)
    end
end

function draw(img::SVG, prim::RectanglePrimitive, idx::Int)
    # SVG will hide rectangles with zero height or width. We'd prefer to have
    # zero width/height rectangles stroked, so this is a work-around.
    width = max(prim.width, 0.01mm)
    height = max(prim.height, 0.01mm)

    x0 = prim.corner[1] + width/2
    y0 = prim.corner[2] + height/2
    translated_path = [(-width/2,-height/2), ( width/2,-height/2),
                       ( width/2, height/2), (-width/2, height/2)]

    indent(img)

    img.indentation += 1
    print(img.out, "<g transform=\"translate(")
    svg_print_float(img.out, x0.value)
    print(img.out, ",")
    svg_print_float(img.out, y0.value)
    print(img.out, ")\"")
    print_vector_properties(img, idx)
    print(img.out, ">\n")
    indent(img)

    print(img.out, "<path d=\"")
    print_svg_path(img.out, translated_path)
    write(img.out, " z\"")
    print(img.out, " class=\"primitive\"")
    print(img.out, "/>\n")

    img.indentation -= 1
    indent(img)
    print(img.out, "</g>\n")
end

function draw(img::SVG, prim::PolygonPrimitive, idx::Int)
    n = length(prim.points)
    n <= 1 && return

    x0, y0 = prim.points[1][1], prim.points[1][2]
    for p in prim.points[2:end]
        x0 += p[1]
        y0 += p[2]
    end
    x0, y0 = x0/n, y0/n
    translated_path = [(p[1]-x0,p[2]-y0) for p in prim.points]

    indent(img)

    img.indentation += 1
    print(img.out, "<g transform=\"translate(")
    svg_print_float(img.out, x0.value)
    print(img.out, ",")
    svg_print_float(img.out, y0.value)
    print(img.out, ")\"")
    print_vector_properties(img, idx)
    print(img.out, ">\n")
    indent(img)

    write(img.out, "<path d=\"")
    print_svg_path(img.out, translated_path, true)
    write(img.out, " z\"")
    print(img.out, " class=\"primitive\"")
    write(img.out, "/>\n")

    img.indentation -= 1
    indent(img)
    print(img.out, "</g>\n")
end

function draw(img::SVG, prim::ComplexPolygonPrimitive, idx::Int)
    Compose.write(img.out, "<path d=\"")
    for ring in prim.rings
        indent(img)
        print_svg_path(img.out, ring, true)
        write(img.out, " ")
    end
    write(img.out, " z\"")
    print_vector_properties(img, idx)
    write(img.out, "/>\n")
end

function draw(img::SVG, prim::CirclePrimitive, idx::Int)
    indent(img)

    img.indentation += 1
    print(img.out, "<g transform=\"translate(")
    svg_print_float(img.out, prim.center[1].value)
    print(img.out, ",")
    svg_print_float(img.out, prim.center[2].value)
    print(img.out, ")\"")
    print_vector_properties(img, idx)
    print(img.out, ">\n")
    indent(img)

    print(img.out, "<circle cx=\"0\" cy=\"0\" r=\"")
    svg_print_float(img.out, prim.radius.value)
    print(img.out, "\"")
    print(img.out, " class=\"primitive\"")
    print(img.out, "/>\n")

    img.indentation -= 1
    indent(img)
    print(img.out, "</g>\n")
end

function draw(img::SVG, prim::EllipsePrimitive, idx::Int)
    cx = prim.center[1].value
    cy = prim.center[2].value
    rx = sqrt((prim.x_point[1].value - cx)^2 +
              (prim.x_point[2].value - cy)^2)
    ry = sqrt((prim.y_point[1].value - cx)^2 +
              (prim.y_point[2].value - cy)^2)
    theta = rad2deg(atan(prim.x_point[2].value - cy,
                         prim.x_point[1].value - cx))

    all(isfinite,[cx, cy, rx, ry, theta]) || return

    indent(img)

    img.indentation += 1
    print(img.out, "<g transform=\"translate(")
    svg_print_float(img.out, cx)
    print(img.out, ",")
    svg_print_float(img.out, cy)
    print(img.out, ")\"")
    print_vector_properties(img, idx)
    print(img.out, ">\n")
    indent(img)

    print(img.out, "<ellipse cx=\"0\" cy=\"0\" rx=\"")
    svg_print_float(img.out, rx)
    print(img.out, "\" ry=\"")
    svg_print_float(img.out, ry)
    print(img.out, '"')
    if abs(theta) > 1e-4
        print(img.out, " transform=\"rotate(")
        svg_print_float(img.out, theta)
        print(img.out, ")\"")

    end
    print(img.out, " class=\"primitive\"")
    print(img.out, "/>\n")

    img.indentation -= 1
    indent(img)
    print(img.out, "</g>\n")
end

function draw(img::SVG, prim::LinePrimitive, idx::Int)
    length(prim.points)<=1 && return

    x0, y0, n = prim.points[1][1], prim.points[1][2], 1
    for p in prim.points[2:end]
        if isfinite(p[1].value) && isfinite(p[2].value)
            x0 += p[1]
            y0 += p[2]
            n +=1
        end
    end
    x0, y0 = x0/n, y0/n
    translated_path = [(p[1]-x0,p[2]-y0) for p in prim.points]

    indent(img)

    img.indentation += 1
    print(img.out, "<g transform=\"translate(")
    svg_print_float(img.out, x0.value)
    print(img.out, ",")
    svg_print_float(img.out, y0.value)
    print(img.out, ")\"")
    print_vector_properties(img, idx, true)
    print(img.out, ">\n")
    indent(img)

    print(img.out, "<path fill=\"none\" d=\"")
    print_svg_path(img.out, translated_path, true)
    print(img.out, "\"")
    print(img.out, " class=\"primitive\"")
    print(img.out, "/>\n")

    img.indentation -= 1
    indent(img)
    print(img.out, "</g>\n")
end

function draw(img::SVG, prim::TextPrimitive, idx::Int)
    indent(img)

    img.indentation += 1
    print(img.out, "<g transform=\"translate(")
    svg_print_float(img.out, prim.position[1].value)
    print(img.out, ",")
    svg_print_float(img.out, prim.position[2].value)
    print(img.out, ")\"")
    print_vector_properties(img, idx, true)
    print(img.out, ">\n")
    indent(img)

    print(img.out, "<g class=\"primitive\">\n")
    img.indentation += 1
    indent(img)
    print(img.out, "<text")

    if prim.halign === hcenter
        print(img.out, " text-anchor=\"middle\"")
    elseif prim.halign === hright
        print(img.out, " text-anchor=\"end\"")
    end

    # NOTE: "dominant-baseline" is the correct way to vertically center text
    # in SVG, but implementations are pretty inconsistent (chrome in particular
    # does a really bad job). We fake it by shifting by some reasonable amount.
    nlines = length(split(prim.value,"\n"))-1
    if prim.valign === vcenter
#        print(img.out, " dy=\"0.35em\"")
        print(img.out, " dy=\"$(svg_fmt_float(0.35-0.6*nlines))em\"")
        #print(img.out, " style=\"dominant-baseline:central\"")
    elseif prim.valign === vtop
        print(img.out, " dy=\"0.6em\"")
        #print(img.out, " style=\"dominant-baseline:text-before-edge\"")
    elseif prim.valign === vbottom
        print(img.out, " dy=\"$(svg_fmt_float(-1.2*nlines))em\"")
    end

    if abs(prim.rot.theta) > 1e-4 || sum(abs.(prim.offset)) > 1e-4mm
        print(img.out, " transform=\"")
        if abs(prim.rot.theta) > 1e-4
            print(img.out, "rotate(")
            svg_print_float(img.out, rad2deg(prim.rot.theta))
            print(img.out, ",")
            svg_print_float(img.out, prim.rot.offset[1].value-prim.position[1].value)
            print(img.out, ", ")
            svg_print_float(img.out, prim.rot.offset[2].value-prim.position[2].value)
            print(img.out, ")")
        end
        if sum(abs.(prim.offset)) > 1e-4mm
            print(img.out, "translate(")
            svg_print_float(img.out, prim.offset[1].value)
            print(img.out, ",")
            svg_print_float(img.out, prim.offset[2].value)
            print(img.out, ")")
        end
        print(img.out, "\"")
    end

    @printf(img.out, ">%s</text>\n",
            svg_newlines(pango_to_svg(prim.value)))

    img.indentation -= 1
    indent(img)
    print(img.out, "</g>\n")
    img.indentation -= 1
    indent(img)
    print(img.out, "</g>\n")
end

function draw(img::SVG, prim::CurvePrimitive, idx::Int)
    x0, y0 = prim.anchor0[1], prim.anchor0[2]
    x0, y0 += prim.ctrl0[1], prim.ctrl0[2]
    x0, y0 += prim.ctrl1[1], prim.ctrl1[2]
    x0, y0 += prim.anchor1[1], prim.anchor1[2]
    x0, y0 = x0/4, y0/4
    translated_anchor0 = [prim.anchor0...] - [x0,y0]
    translated_ctrl0 = [prim.ctrl0...] - [x0,y0]
    translated_ctrl1 = [prim.ctrl1...] - [x0,y0]
    translated_anchor1 = [prim.anchor1...] - [x0,y0]

    indent(img)

    img.indentation += 1
    print(img.out, "<g transform=\"translate(")
    svg_print_float(img.out, x0.value)
    print(img.out, ",")
    svg_print_float(img.out, y0.value)
    print(img.out, ")\"")
    print_vector_properties(img, idx, true)
    print(img.out, ">\n")
    indent(img)

    print(img.out, "<path fill=\"none\" d=\"M")
    svg_print_float(img.out, translated_anchor0[1].value)
    print(img.out, ',')
    svg_print_float(img.out, translated_anchor0[2].value)
    print(img.out, " C")
    svg_print_float(img.out, translated_ctrl0[1].value)
    print(img.out, ',')
    svg_print_float(img.out, translated_ctrl0[2].value)
    print(img.out, ' ')
    svg_print_float(img.out, translated_ctrl1[1].value)
    print(img.out, ',')
    svg_print_float(img.out, translated_ctrl1[2].value)
    print(img.out, ' ')
    svg_print_float(img.out, translated_anchor1[1].value)
    print(img.out, ',')
    svg_print_float(img.out, translated_anchor1[2].value)
    print(img.out, '"')
    print(img.out, " class=\"primitive\"")
    print(img.out, "/>\n")

    img.indentation -= 1
    indent(img)
    print(img.out, "</g>\n")
end

function draw(img::SVG, prim::BitmapPrimitive, idx::Int)
    indent(img)

    img.indentation += 1
    print(img.out, "<g transform=\"translate(")
    svg_print_float(img.out, prim.corner[1].value)
    print(img.out, ",")
    svg_print_float(img.out, prim.corner[2].value)
    print(img.out, ")\"")
    print_vector_properties(img, idx, true)
    print(img.out, ">\n")
    indent(img)

    print(img.out, "<image x=\"0\" y=\"0\" width=\"")
    svg_print_float(img.out, prim.width.value)
    print(img.out, "\" height=\"")
    svg_print_float(img.out, prim.height.value)
    print(img.out, '"')
    print_vector_properties(img, idx)

    print(img.out, " xlink:href=\"data:", prim.mime, ";base64,")
    b64pipe = Base64EncodePipe(img.out)
    write(b64pipe, prim.data)
    close(b64pipe)
    print(img.out, "\"")
    print(img.out, " class=\"primitive\"")
    print(img.out, "/>\n")

    img.indentation -= 1
    indent(img)
    print(img.out, "</g>\n")
end

# FormBatch Drawing
# -----------------

function draw(img::SVG, batch::FormBatch)
    id = genid(img)
    push!(img.batches, (batch.primitive, id))
    for i in 1:length(batch.offsets)
        indent(img)
        print(img.out, "<use xlink:href=\"#", id, "\" x=\"")
        svg_print_float(img.out, batch.offsets[i][1].value)
        print(img.out, "\" y=\"")
        svg_print_float(img.out, batch.offsets[i][2].value)
        print(img.out, "\"")
        print_vector_properties(img, i)
        print(img.out, "/>\n")
    end
end

# Applying properties
# -------------------

# Return a URL corresponding to a ClipPrimitive
clippathurl(img::SVG, property::ClipPrimitive) = get!(() -> genid(img), img.clippaths, property)

function push_property_frame(img::SVG, properties::Vector{Property})
    isempty(properties) && return
    svgalphatest(properties)

    frame = SVGPropertyFrame()
    applied_properties = Set{Type}()
    scalar_properties = Array{Property}(undef, 0)
    isplotpanel = false
    for property in properties
        if !isrepeatable(property) && (typeof(property) in applied_properties)
            continue
        elseif isscalar(property) && isa(property, Clip)
            # clip-path needs to be in it's own group. Otherwise it can cause
            # problems if we apply a transform to the group.
            indent(img)
            img.indentation += 1
            write(img.out, "<g")
            print_property(img, property.primitives[1])
            write(img.out, ">\n")
            frame.has_scalar_clip = true
        elseif isscalar(property)
            push!(scalar_properties, property)
            push!(applied_properties, typeof(property))
            frame.has_scalar_properties = true
            img.vector_properties[typeof(property)] = nothing
        else
            frame.vector_properties[typeof(property)] = property
            img.vector_properties[typeof(property)] = property
        end
        isplotpanel |= typeof(property)==SVGClass &&
                length(property.primitives)==1 && property.primitives[1].value=="plotpanel"
    end
    push!(img.property_stack, frame)
    isempty(scalar_properties) && return

    id_needed = any([isa(property, JSCall) for property in scalar_properties])
    for property in scalar_properties
        if isa(property, SVGID)
            img.current_id = property.primitives[1].value
            img.has_current_id = true
        end
    end

    if !img.has_current_id
        img.current_id = genid(img)
        push!(scalar_properties, svgid(img.current_id))
        img.has_current_id = true
    end

    indent(img)
    write(img.out, "<g")
    for property in scalar_properties
        print_property(img, property.primitives[1])
    end
    write(img.out, ">\n");
    img.has_current_id = false
    img.indentation += 1

    if isplotpanel
        indent(img)
        write(img.out, "<metadata>\n")
        boundingBox = string(img.panelcoords[1].x0[1], ' ', img.panelcoords[1].x0[2], ' ',
                          img.panelcoords[1].a[1], ' ', img.panelcoords[1].a[2])
        unitBox = string(img.panelcoords[2].x0, ' ', img.panelcoords[2].y0, ' ',
                          img.panelcoords[2].width, ' ', img.panelcoords[2].height)
        indent(img)
        write(img.out, "  <boundingbox value=\"$boundingBox\"/>\n")
        indent(img)
        write(img.out, "  <unitbox value=\"$unitBox\"/>\n")
        indent(img)
        write(img.out, "</metadata>\n")
    end
end

function pop_property_frame(img::SVG)
    @assert !isempty(img.property_stack)
    frame = pop!(img.property_stack)

    if frame.has_scalar_properties
        img.indentation -= 1
        indent(img)
        write(img.out, "</g>")
        frame.has_link && write(img.out, "</a>")
        frame.has_mask && write(img.out, "</mask>")
        write(img.out, "\n")
    end

    if frame.has_scalar_clip
        img.indentation -= 1
        indent(img)
        write(img.out, "</g>\n")
    end

    for (propertytype, property) in frame.vector_properties
        img.vector_properties[propertytype] = nothing
        for i in length(img.property_stack):-1:1
            if haskey(img.property_stack[i].vector_properties, propertytype)
                img.vector_properties[propertytype] =
                    img.property_stack[i].vector_properties[propertytype]
            end
        end
    end
end

function Compose.draw(img::SVG, prim::ArcPrimitive, idx::Int)
    xc = prim.center[1].value
    yc = prim.center[2].value
    rx = ry = prim.radius.value
    x1 =  rx*cos(prim.angle1)
    y1 =  ry*sin(prim.angle1)
    x2 =  rx*cos(prim.angle2)
    y2 =  ry*sin(prim.angle2)
    dθ = prim.angle2 - prim.angle1
    dθ += 2π*(dθ<0)

    indent(img)
    img.indentation += 1
    print(img.out, "<g transform=\"translate(")
    svg_print_float(img.out, xc)
    print(img.out, ",")
    svg_print_float(img.out, yc)
    print(img.out, ")\"")
    print_vector_properties(img, idx)
    print(img.out, ">\n")
    indent(img)

    print(img.out, "<path d=\"M")
    prim.sector && print(img.out,"0,0 L")
    svg_print_float(img.out, x1)
    print(img.out, ",")
    svg_print_float(img.out, y1)
    print(img.out, " A")
    svg_print_float(img.out, rx)
    print(img.out, ",")
    svg_print_float(img.out, ry)
    print(img.out, ' ', 0, ' ', dθ>π ? 1 : 0,' ',1,' ')
    svg_print_float(img.out, x2)
    print(img.out, ",")
    svg_print_float(img.out, y2)
    prim.sector && print(img.out, "L0,0")
    print(img.out, '"')
    print(img.out, " class=\"primitive\"")
    print(img.out, "/>\n")

    img.indentation -= 1
    indent(img)
    print(img.out, "</g>\n")
end

# Currently for svg, you can't use a transparent fill(color) and fillopacity() together,
# so throw a warning if a user tries to do that.

function svgalphatest(properties::Vector{Property})
    has_fill_opacity = any(isa.(properties, Property{FillOpacityPrimitive}))
    !has_fill_opacity && return
    is_fill = isa.(properties, Property{FillPrimitive})
    !any(is_fill) && return
    fillproperties = properties[is_fill][1]
    has_alpha = any([alpha(x.color) for x in fillproperties.primitives].<1.0)
    has_alpha && @warn "For svg transparent colors, use either e.g. fill(RGBA(r,g,b,a)) or fillopacity(a), but not both."
end


text:
using Pkg
using Cascadia: parsehtml, Selector, getattr
using HTTP: HTTP, request
using CSV

data = CSV.File(joinpath("data", "FreeGermanTextbooks.csv") |> CSV.DataFrame
DOIs = SubString.(data[!,Symbol("DOI URL")], 16)
fields = data[!,Symbol("English Package Name")]
titles = data[!,Symbol("Book Title")]
URLs = replace.(data[!,Symbol("OpenURL")], "http://" => "https://")

const SELECTOR_EPUB = Selector(".c-button[title=\"Download this book in EPUB format\"]")
const SELECTOR_PDF = Selector(".c-button[title=\"Download this book in PDF format\"]")

function links(url)
    response = HTTP.get(url)
    html = parsehtml(String(response.body))
    epub = eachmatch(SELECTOR_EPUB, html.root)
    pdf = eachmatch(SELECTOR_PDF, html.root)
    (ISBN = match(r"(?<=isbn=).*", url).match,
     epub = isempty(epub) ? "" : string(HTTP.URI(scheme = "https",
                                                 host = HTTP.header(response.request, "Host"),
                                                 path = getattr(epub[1], "href"))),
     pdf = isempty(pdf) ? "" : string(HTTP.URI(scheme = "https",
                                               host = HTTP.header(response.request, "Host"),
                                               path = getattr(pdf[1], "href"))),
    )
end

output = Vector{NamedTuple{(:ISBN, :epub, :pdf), NTuple{3, String}}}()

for url in URLs
    println(url)
    push!(output, links(url))
end
isbn_links = CSV.DataFrame(output)
awesome = CSV.join(data, isbn_links, on = [Symbol("Electronic ISBN") => :ISBN])
CSV.write(joinpath(@__DIR__, "data", "FreeGermanTextbooksEnhanced.tsv"), awesome, delim = '\t')

ext: jl


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    isdefined(Luxor, Symbol("#@layer")) && precompile(Tuple{getfield(Luxor, Symbol("#@layer")), LineNumberNode, Module, Int})
    isdefined(Luxor, Symbol("#@setcolor_str")) && precompile(Tuple{getfield(Luxor, Symbol("#@setcolor_str")), LineNumberNode, Module, Int})
    isdefined(Luxor, Symbol("#box##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#box##kw")), NamedTuple{(:vertices,), Tuple{Bool}}, typeof(Luxor.box), Luxor.BoundingBox})
    isdefined(Luxor, Symbol("#box##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#box##kw")), NamedTuple{(:vertices,), Tuple{Bool}}, typeof(Luxor.box), Luxor.Point, Float64, Float64})
    isdefined(Luxor, Symbol("#dimension##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#dimension##kw")), NamedTuple{(:fromextension,), Tuple{Array{Int64, 1}}}, typeof(Luxor.dimension), Luxor.Point, Luxor.Point})
    isdefined(Luxor, Symbol("#dimension##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#dimension##kw")), NamedTuple{(:offset, :format), Tuple{Int64, typeof(Base.string)}}, typeof(Luxor.dimension), Luxor.Point, Luxor.Point})
    isdefined(Luxor, Symbol("#dimension##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#dimension##kw")), NamedTuple{(:offset, :fromextension, :toextension, :textverticaloffset, :texthorizontaloffset, :textgap), Tuple{Int64, Array{Int64, 1}, Array{Int64, 1}, Float64, Int64, Int64}}, typeof(Luxor.dimension), Luxor.Point, Luxor.Point})
    isdefined(Luxor, Symbol("#intersectionlines##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#intersectionlines##kw")), NamedTuple{(:crossingonly,), Tuple{Bool}}, typeof(Luxor.intersectionlines), Luxor.Point, Luxor.Point, Luxor.Point, Luxor.Point})
    isdefined(Luxor, Symbol("#ngon##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#ngon##kw")), NamedTuple{(:reversepath,), Tuple{Bool}}, typeof(Luxor.ngon), Luxor.Point, Float64, Int64, Int64, Symbol})
    isdefined(Luxor, Symbol("#ngon##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#ngon##kw")), NamedTuple{(:vertices,), Tuple{Bool}}, typeof(Luxor.ngon), Luxor.Point, Float64, Int64, Float64})
    isdefined(Luxor, Symbol("#ngonside##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#ngonside##kw")), NamedTuple{(:vertices,), Tuple{Bool}}, typeof(Luxor.ngonside), Luxor.Point, Int64, Int64})
    isdefined(Luxor, Symbol("#poly##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#poly##kw")), NamedTuple{(:close,), Tuple{Bool}}, typeof(Luxor.poly), Array{Luxor.Point, 1}, Symbol})
    isdefined(Luxor, Symbol("#squircle##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#squircle##kw")), NamedTuple{(:vertices,), Tuple{Bool}}, typeof(Luxor.squircle), Luxor.Point, Float64, Float64})
    isdefined(Luxor, Symbol("#star##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#star##kw")), NamedTuple{(:reversepath,), Tuple{Bool}}, typeof(Luxor.star), Float64, Float64, Float64, Int64, Float64, Float64, Symbol})
    isdefined(Luxor, Symbol("#star##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#star##kw")), NamedTuple{(:vertices,), Tuple{Bool}}, typeof(Luxor.star), Luxor.Point, Float64, Int64, Float64, Int64})
    isdefined(Luxor, Symbol("#star##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#star##kw")), NamedTuple{(:vertices,), Tuple{Bool}}, typeof(Luxor.star), Luxor.Point, Int64, Int64, Float64, Float64})
    isdefined(Luxor, Symbol("#star##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#star##kw")), NamedTuple{(:vertices,), Tuple{Bool}}, typeof(Luxor.star), Luxor.Point, Int64, Int64, Float64})
    isdefined(Luxor, Symbol("#text##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#text##kw")), NamedTuple{(:halign,), Tuple{Symbol}}, typeof(Luxor.text), String, Float64, Float64})
    isdefined(Luxor, Symbol("#text##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#text##kw")), NamedTuple{(:halign,), Tuple{Symbol}}, typeof(Luxor.text), String, Luxor.Point})
    isdefined(Luxor, Symbol("#textcurvecentered##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#textcurvecentered##kw")), NamedTuple{(:clockwise, :letter_spacing, :baselineshift), Tuple{Bool, Int64, Int64}}, typeof(Luxor.textcurvecentered), String, Float64, Int64, Luxor.Point})
    isdefined(Luxor, Symbol("#textwrap##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#textwrap##kw")), NamedTuple{(:rightgutter,), Tuple{Int64}}, typeof(Luxor.textwrap), String, Float64, Luxor.Point, typeof(identity)})
    isdefined(Luxor, Symbol("#textwrap##kw")) && precompile(Tuple{getfield(Luxor, Symbol("#textwrap##kw")), NamedTuple{(:rightgutter,), Tuple{Int64}}, typeof(Luxor.textwrap), String, Float64, Luxor.Point})
    precompile(Tuple{typeof(Luxor.arc), Luxor.Point, Float64, Float64, Float64})
    precompile(Tuple{typeof(Luxor.arc), Luxor.Point, Int64, Float64, Float64})
    precompile(Tuple{typeof(Luxor.background), String})
    precompile(Tuple{typeof(Luxor.between), Luxor.Point, Luxor.Point, Float64})
    precompile(Tuple{typeof(Luxor.bezierfrompoints), Luxor.Point, Luxor.Point, Luxor.Point, Luxor.Point})
    precompile(Tuple{typeof(Luxor.blendadjust), Cairo.CairoPattern, Luxor.Point, Float64, Float64, Float64})
    precompile(Tuple{typeof(Luxor.box), Luxor.Point, Float64, Float64, Symbol})
    precompile(Tuple{typeof(Luxor.box), Luxor.Point, Int64, Int64, Symbol})
    precompile(Tuple{typeof(Luxor.buildboxmap), Array{Float64, 1}, Float64, Float64, Float64, Float64})
    precompile(Tuple{typeof(Luxor.buildboxmap), Array{Float64, 1}, Float64, Float64, Float64, Int64})
    precompile(Tuple{typeof(Luxor.carc), Luxor.Point, Float64, Float64, Float64})
    precompile(Tuple{typeof(Luxor.carc), Luxor.Point, Int64, Float64, Float64})
    precompile(Tuple{typeof(Luxor.circle), Luxor.Point, Float64, Symbol})
    precompile(Tuple{typeof(Luxor.circle), Luxor.Point, Int64, Symbol})
    precompile(Tuple{typeof(Luxor.cropmarks), Luxor.Point, Float64, Float64})
    precompile(Tuple{typeof(Luxor.easeinoutcubic), Int64, Int64, Int64, Int64})
    precompile(Tuple{typeof(Luxor.easeinoutquad), Int64, Int64, Int64, Int64})
    precompile(Tuple{typeof(Luxor.easeinoutquint), Int64, Int64, Int64, Int64})
    precompile(Tuple{typeof(Luxor.easeinsine), Int64, Int64, Int64, Int64})
    precompile(Tuple{typeof(Luxor.ellipse), Luxor.Point, Float64, Float64, Symbol})
    precompile(Tuple{typeof(Luxor.ellipse), Luxor.Point, Int64, Int64, Symbol})
    precompile(Tuple{typeof(Luxor.ellipse), Luxor.Point, Luxor.Point, Int64, Symbol})
    precompile(Tuple{typeof(Luxor.insertvertices!), Array{Luxor.Point, 1}})
    precompile(Tuple{typeof(Luxor.intersectionlinecircle), Luxor.Point, Luxor.Point, Luxor.Point, Float64})
    precompile(Tuple{typeof(Luxor.intersectionlines), Luxor.Point, Luxor.Point, Luxor.Point, Luxor.Point})
    precompile(Tuple{typeof(Luxor.ispointinsidetriangle), Luxor.Point, Array{Luxor.Point, 1}})
    precompile(Tuple{typeof(Luxor.ispolyclockwise), Array{Luxor.Point, 1}})
    precompile(Tuple{typeof(Luxor.line), Luxor.Point, Luxor.Point, Symbol})
    precompile(Tuple{typeof(Luxor.line), Luxor.Point})
    precompile(Tuple{typeof(Luxor.midpoint), Luxor.Point, Luxor.Point})
    precompile(Tuple{typeof(Luxor.nearestindex), Array{Float64, 1}, Float64})
    precompile(Tuple{typeof(Luxor.ngon), Int64, Int64, Int64, Int64, Int64, Symbol})
    precompile(Tuple{typeof(Luxor.ngon), Luxor.Point, Float64, Int64, Int64, Symbol})
    precompile(Tuple{typeof(Luxor.polar), Int64, Float64})
    precompile(Tuple{typeof(Luxor.polyremovecollinearpoints), Array{Luxor.Point, 1}})
    precompile(Tuple{typeof(Luxor.polysample), Array{Luxor.Point, 1}, Int64})
    precompile(Tuple{typeof(Luxor.prettypoly), Array{Luxor.Point, 1}, Symbol, typeof(identity)})
    precompile(Tuple{typeof(Luxor.randomhue)})
    precompile(Tuple{typeof(Luxor.randompointarray), Float64, Float64, Float64, Float64, Int64})
    precompile(Tuple{typeof(Luxor.rect), Int64, Float64, Float64, Float64, Symbol})
    precompile(Tuple{typeof(Luxor.rect), Int64, Int64, Float64, Float64, Symbol})
    precompile(Tuple{typeof(Luxor.sector), Float64, Float64, Float64, Float64, Int64, Symbol})
    precompile(Tuple{typeof(Luxor.sector), Int64, Int64, Float64, Int64, Symbol})
    precompile(Tuple{typeof(Luxor.sector), Luxor.Point, Int64, Int64, Float64, Int64, Symbol})
    precompile(Tuple{typeof(Luxor.set_current_alpha), Float64})
    precompile(Tuple{typeof(Luxor.set_current_bluevalue), Float64})
    precompile(Tuple{typeof(Luxor.set_current_greenvalue), Float64})
    precompile(Tuple{typeof(Luxor.set_current_redvalue), Float64})
    precompile(Tuple{typeof(Luxor.setantialias), Int64})
    precompile(Tuple{typeof(Luxor.setcolor), ColorTypes.RGB{Float64}})
    precompile(Tuple{typeof(Luxor.sethue), ColorTypes.RGB{Float64}})
    precompile(Tuple{typeof(Luxor.setline), Float64})
    precompile(Tuple{typeof(Luxor.setline), Int64})
    precompile(Tuple{typeof(Luxor.setlinecap), String})
    precompile(Tuple{typeof(Luxor.setopacity), Float64})
    precompile(Tuple{typeof(Luxor.squircle), Luxor.Point, Float64, Float64, Symbol})
    precompile(Tuple{typeof(Luxor.star), Float64, Float64, Float64, Int64, Float64, Float64, Symbol})
    precompile(Tuple{typeof(Luxor.star), Float64, Float64, Int64, Int64, Float64, Int64, Symbol})
    precompile(Tuple{typeof(Luxor.text), String, Float64, Float64})
    precompile(Tuple{typeof(Luxor.text), String})
    precompile(Tuple{typeof(Luxor.textcentered), String, Int64, Int64})
    precompile(Tuple{typeof(Luxor.textcentered), String, Luxor.Point})
    precompile(Tuple{typeof(Luxor.translate), Float64, Float64})
    precompile(Tuple{typeof(Luxor.translate), Int64, Int64})
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    SerialExecution

Functor object for using serial execution in a lshaped algorithm. Create by supplying a [`Serial`](@ref) object through `execution` in the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

"""
struct SerialExecution{T <: AbstractFloat,
                       A <: AbstractVector,
                       F <: AbstractFeasibilityAlgorithm,
                       I <: AbstractIntegerAlgorithm} <: AbstractLShapedExecution
    subproblems::Vector{SubProblem{T,F,I}}
    decisions::DecisionMap
    subobjectives::A
    model_objectives::A
    metadata::MetaData

    function SerialExecution(structure::StageDecompositionStructure{2, 1, <:Tuple{ScenarioProblems}},
                             feasibility_strategy::AbstractFeasibilityStrategy,
                             integer_strategy::AbstractIntegerStrategy,
                             ::Type{T},
                             ::Type{A}) where {T <: AbstractFloat,
                                               A <: AbstractVector}
        F = worker_type(feasibility_strategy)
        I = worker_type(integer_strategy)
        execution = new{T,A,F,I}(Vector{SubProblem{T,F,I}}(), structure.decisions[2], A(), A(), MetaData())
        # Load subproblems
        for i in 1:num_subproblems(structure, 2)
            push!(execution.subproblems, SubProblem(
                subproblem(structure, 2, i),
                i,
                T(probability(scenario(structure, 2, i))),
                feasibility_strategy,
                integer_strategy))
        end
        return execution
    end
end

function num_thetas(lshaped::AbstractLShaped, ::SerialExecution)
    return num_thetas(num_subproblems(lshaped), lshaped.aggregation)
end

function initialize_subproblems!(execution::SerialExecution{T},
                                 scenarioproblems::ScenarioProblems,
                                 feasibility_strategy::AbstractFeasibilityStrategy,
                                 integer_strategy::AbstractIntegerStrategy) where T <: AbstractFloat
    for i in 1:num_subproblems(scenarioproblems)
        push!(execution.subproblems, SubProblem(
            subproblem(scenarioproblems, i),
            i,
            T(probability(scenario(scenarioproblems, i))),
            feasibility_strategy,
            integer_strategy))
    end
    return nothing
end

function mutate_subproblems!(mutator::Function, execution::SerialExecution)
    for subproblem in execution.subproblems
        mutator(subproblem)
    end
    return nothing
end

function finish_initilization!(lshaped::AbstractLShaped, execution::SerialExecution)
    append!(execution.subobjectives, fill(1e10, num_thetas(lshaped)))
    append!(execution.model_objectives, fill(-1e10, num_thetas(lshaped)))
    return nothing
end

function restore_subproblems!(lshaped::AbstractLShaped, execution::SerialExecution)
    for subproblem in execution.subproblems
        restore_subproblem!(subproblem)
    end
    return nothing
end

function resolve_subproblems!(lshaped::AbstractLShaped, execution::SerialExecution{T}) where T <: AbstractFloat
    # Update subproblems
    update_known_decisions!(execution.decisions, lshaped.x)
    # Assume no cuts are added
    added = false
    # Update and solve subproblems
    for subproblem in execution.subproblems
        update_subproblem!(subproblem)
        cut::SparseHyperPlane{T} = subproblem(lshaped.x, execution.metadata)
        added |= aggregate_cut!(lshaped, lshaped.aggregation, cut)
    end
    added |= flush!(lshaped, lshaped.aggregation)
    # Return current objective value and cut_added flag
    return current_objective_value(lshaped), added
end

# API
# ------------------------------------------------------------
function (execution::Serial)(structure::StageDecompositionStructure{2, 1, <:Tuple{ScenarioProblems}},
                             feasibility_strategy::AbstractFeasibilityStrategy,
                             integer_strategy::AbstractIntegerStrategy,
                             ::Type{T},
                             ::Type{A}) where {T <: AbstractFloat,
                                               A <: AbstractVector}
    return SerialExecution(structure,
                           feasibility_strategy,
                           integer_strategy,
                           T,
                           A)
end

function str(::Serial)
    return ""
end


text:
# Some Boiler-plate code for preconditioning
#
# Meaning of P:
#    P ≈ ∇²E, so the preconditioned gradient is P^{-1} ∇E
#    P can be an arbitrary type but at a minimum it MUST provide:
#        A_ldiv_B!(x, P, b) ->  x = P \ b
#        dot(x, P, b)  ->  x' P b
#
#    If `dot` is not provided, then Optim.jl will try to define it via
#        dot(x, P, b) = vecdot(x, A_mul_B!(similar(x), P, y))
#
#    finally the preconditioner can be updated after each x-update using
#        precondprep!
#    but this is passed as an argument at the moment!
#

#####################################################
#  [0] Defaults and aliases for easier reading of the code
#      these can also be over-written if necessary.

# an inner product w.r.t. a metric P (=preconditioner)
Base.dot(x, P, y) = vecdot(x, A_mul_B!(similar(x), P, y))

# default preconditioner update
precondprep!(P, x) = nothing

#####################################################
#  [1] Empty preconditioner = Identity
#

# out =  P^{-1} * A
Base.A_ldiv_B!(out, ::Void, A) = copy!(out, A)
# A' * P B
Base.dot(A, ::Void, B) = vecdot(A, B)

#####################################################
#  [2] Diagonal preconditioner
#      P = Diag(d)
#      unfortunately, Base does not implement
#      A_ldiv_B!(a, P, b) or A_mul_B! for this type, so we do it by hand
#      TODO: maybe implement this in Base
Base.A_ldiv_B!(out::Array, P::Diagonal, A::Array) = copy!(out, A ./ P.diag)
Base.dot(A::Array, P::Diagonal, B::Array) = vecdot(A, P.diag .* B)

#####################################################
#  [3] Inverse Diagonal preconditioner
#      here, P is stored by the entries of its inverse
#      TODO: maybe implement this in Base?
type InverseDiagonal
   diag
end
Base.A_ldiv_B!(out::Array, P::InverseDiagonal, A::Array) = copy!(out, A .* P.diag)
Base.dot(A::Array, P::InverseDiagonal, B::Vector) = vecdot(A, B ./ P.diag)

#####################################################
#  [4] Matrix Preconditioner
#  the assumption here is that P is given by its inverse, which is typical
#     > A_ldiv_B! is about to be moved to Base, so we need a temporary hack
#     > A_mul_B! is already in Base, which defines `dot`
#  nothing to do!
Base.A_ldiv_B!(x, P::AbstractMatrix, b) = copy!(x, P \ b)


text:
# """
#    fill_and_evaluate(net, batch; batch_size, fill_batches)

# Evaluate a neural network on a batch of inputs.

# If `fill_batches=true`, the batch is padded with dummy inputs until
# it has size `batch_size` before it is sent to the network.

# This function is typically called by inference servers.
# """
function fill_and_evaluate(net, batch; batch_size, fill_batches)
  n = length(batch)
  @assert n > 0
  if !fill_batches
    return evaluate_batch(net, batch)
  else
    nmissing = batch_size - n
    @assert nmissing >= 0
    if nmissing == 0
      return evaluate_batch(net, batch)
    else
      batch = vcat(batch, [batch[1] for _ in 1:nmissing])
      return evaluate_batch(net, batch)[1:n]
    end
  end
end

# function launch_inference_server(
#   net1::MuNetwork,
#   net2::MuNetwork;
#   num_workers,
#   batch_size,
#   fill_batches)

#   return Batchifier.launch_server(;num_workers, batch_size) do batch
#     n = length(batch) # batch is vector Any, cause h(O), f(S), and g(S,A) are evaluated

#     maskf1 = findall(b->b.netid == :f1, batch)
#     maskg1 = findall(b->b.netid == :g1, batch)
#     maskh1 = findall(b->b.netid == :h1, batch)

#     maskf2 = findall(b->b.netid == :f2, batch)
#     maskg2 = findall(b->b.netid == :g2, batch)
#     maskh2 = findall(b->b.netid == :h2, batch)
#     @assert length(maskf1) + length(maskg1) + length(maskh1) +
#             length(maskf2) + length(maskg2) + length(maskh2) == n

#     batchf1 = [b.query for b in batch[maskf1]]
#     batchg1 = [b.query for b in batch[maskg1]]
#     batchh1 = [b.query for b in batch[maskh1]]

#     batchf2 = [b.query for b in batch[maskf2]]
#     batchg2 = [b.query for b in batch[maskg2]]
#     batchh2 = [b.query for b in batch[maskh2]]

#     resf1 = isempty(maskf1) ? Nothing[] : fill_and_evaluate(net1.f, batchf1; batch_size=n, fill_batches)
#     resg1 = isempty(maskg1) ? Nothing[] : fill_and_evaluate(net1.g, batchg1; batch_size=n, fill_batches)
#     resh1 = isempty(maskh1) ? Nothing[] : fill_and_evaluate(net1.h, batchh1; batch_size=n, fill_batches)

#     resf2 = isempty(maskf2) ? Nothing[] : fill_and_evaluate(net2.f, batchf2; batch_size=n, fill_batches)
#     resg2 = isempty(maskg2) ? Nothing[] : fill_and_evaluate(net2.g, batchg2; batch_size=n, fill_batches)
#     resh2 = isempty(maskh2) ? Nothing[] : fill_and_evaluate(net2.h, batchh2; batch_size=n, fill_batches)

#     res = Vector{Union{eltype.((resf1,resg1,resh1,resf2,resg2,resh2,))...}}(undef, n)

#     res[maskf1] = resf1
#     res[maskg1] = resg1
#     res[maskh1] = resh1

#     res[maskf2] = resf2
#     res[maskg2] = resg2
#     res[maskh2] = resh2
#     return res
#   end
# end

# function launch_inference_server(
#   net::MuNetwork;
#   num_workers,
#   batch_size,
#   fill_batches)

#   return Batchifier.launch_server(;num_workers, batch_size) do batch
#     n = length(batch) # batch is vector Any, cause h(O), f(S), and g(S,A) are evaluated
#     maskf = findall(b->b.netid == :f, batch)
#     maskg = findall(b->b.netid == :g, batch)
#     maskh = findall(b->b.netid == :h, batch)
#     @assert length(maskf) + length(maskg) + length(maskh) == n

#     batchf = [b.query for b in batch[maskf]]
#     batchg = [b.query for b in batch[maskg]]
#     batchh = [b.query for b in batch[maskh]]

#     resf = isempty(maskf) ? Nothing[] : fill_and_evaluate(net.f, batchf; batch_size=n, fill_batches)
#     resg = isempty(maskg) ? Nothing[] : fill_and_evaluate(net.g, batchg; batch_size=n, fill_batches)
#     resh = isempty(maskh) ? Nothing[] : fill_and_evaluate(net.h, batchh; batch_size=n, fill_batches)

#     res = Vector{Union{eltype.((resf,resg,resh))...}}(undef, n)
#     res[maskf] = resf
#     res[maskg] = resg
#     res[maskh] = resh
#     return res
#   end
# end

function launch_inference_server(
  net::Tuple{InitialOracle, RecurrentOracle};
  num_workers,
  batch_size,
  fill_batches)

  return Batchifier.launch_server(;num_workers, batch_size) do batch
    n = length(batch) # batch is vector Any, cause h(O), f(S), and g(S,A) are evaluated
    mask_i = findall(b->b.netid == :i, batch)
    mask_r = findall(b->b.netid == :r, batch)
    @assert length(mask_i) + length(mask_r) == n

    batch_i = [b.query for b in batch[mask_i]]
    batch_r = [b.query for b in batch[mask_r]]

    if isempty(mask_i)
      return fill_and_evaluate(net[2], batch_r; batch_size=n, fill_batches)
    elseif isempty(mask_r)
      return fill_and_evaluate(net[1], batch_i; batch_size=n, fill_batches)
    else
      res_i = fill_and_evaluate(net[1], batch_i; batch_size=n, fill_batches)
      res_r = fill_and_evaluate(net[2], batch_r; batch_size=n, fill_batches)
      @assert typeof(res_i) == typeof(res_r)
      res = Vector{eltype(res_i)}(undef, n)
      res[mask_i] = res_i
      res[mask_r] = res_r
      return res
    end
  end
end

function launch_inference_server(
  net1::Tuple{InitialOracle, RecurrentOracle},
  net2::Tuple{InitialOracle, RecurrentOracle};
  num_workers,
  batch_size,
  fill_batches)

  return Batchifier.launch_server(;num_workers, batch_size) do batch
    n = length(batch) # batch is vector Any, cause h(O), f(S), and g(S,A) are evaluated
    mask_i  = findall(b->b.netid == :i, batch)
    mask_r  = findall(b->b.netid == :r, batch)
    mask_i2 = findall(b->b.netid == :i2, batch)
    mask_r2 = findall(b->b.netid == :r2, batch)
    @assert length(mask_i) + length(mask_r) + length(mask_i2) + length(mask_r2) == n

    batch_i = [b.query for b in batch[mask_i]]
    batch_r = [b.query for b in batch[mask_r]]
    batch_i2 = [b.query for b in batch[mask_i2]]
    batch_r2 = [b.query for b in batch[mask_r2]]

    res_r = isempty(mask_r) ? Nothing[] : fill_and_evaluate(net1[2], batch_r; batch_size=n, fill_batches)
    res_r2 = isempty(mask_r2) ? Nothing[] : fill_and_evaluate(net2[2], batch_r2; batch_size=n, fill_batches)
    res_i = isempty(mask_i) ? Nothing[] : fill_and_evaluate(net1[1], batch_i; batch_size=n, fill_batches)
    res_i2 = isempty(mask_i2) ? Nothing[] : fill_and_evaluate(net2[1], batch_i2; batch_size=n, fill_batches)
    # @info "types" typeof(res_r) == typeof(res_r2)
    # @assert typeof(res_i) == typeof(res_r) == typeof(res_i2) == typeof(res_r2)
    res = Vector{Union{eltype.((res_r,res_r2,res_i,res_i2))...}}(undef, n)
    res[mask_i] = res_i
    res[mask_r] = res_r
    res[mask_i2] = res_i2
    res[mask_r2] = res_r2

    return res
  end
end

ret_oracle(x) = () -> x
do_nothing!() = nothing
send_done!(reqcs...) = () -> foreach(Batchifier.client_done!, reqcs)
zipthunk(f1, f2) = () -> (f1(), f2())

function batchify_oracles(o::MuNetwork; kwargs...)
  reqc = launch_inference_server(o; kwargs...)
  make() = (;
  g = Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:g)),
  f = Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:f)),
  h = Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:h)),
  )
  return make, send_done!(reqc)
end

# function batchify_oracles(os::Tuple{MuNetwork,MuNetwork}; kwargs...)
#   reqc = launch_inference_server(os[1],os[2]; kwargs...)
#   make1() = (;
#   g = Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:g1)),
#   f = Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:f1)),
#   h = Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:h1)),
#   )
#   make2() = (;
#   g = Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:g2)),
#   f = Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:f2)),
#   h = Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:h2)),
#   )
#   return zipthunk(make1,make2), send_done!(reqc)
# end

# function batchify_oracles(os::Tuple{<:Any, AbstractNetwork}; kwargs...)
#   reqc = launch_inference_server(os[2]; kwargs...)
#   make2() = (;
#   g = Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:g)),
#   f = Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:f)),
#   h = Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:h)),
#   )
#   return zipthunk(ret_oracle(os[1]), make2), send_done!(reqc)
# end

# function batchify_oracles(os::Tuple{AbstractNetwork, <:Any}; kwargs...)
#   reqc = launch_inference_server(os[1]; kwargs...)
#   make1() = (;
#   g = Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:g)),
#   f = Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:f)),
#   h = Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:h)),
#   )
#   return zipthunk(make1, ret_oracle(os[2])), send_done!(reqc)
# end

function batchify_oracles(o::Any; kwargs...)
  return ret_oracle(o), do_nothing!
end

function batchify_oracles(o::Tuple{InitialOracle,RecurrentOracle}; kwargs...)
  reqc = launch_inference_server(o; kwargs...)
  make = make_make(reqc)
  return make, send_done!(reqc)
end

function make_make(reqc)
  return () -> (
    Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:i)),
    Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:r))
  )
end

function make_make2(reqc)
  return () -> (
    Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:i2)),
    Batchifier.BatchedOracle(reqc, q -> (query=q, netid=:r2))
  )
end

# Two players
function batchify_oracles(os::Tuple{Tuple{InitialOracle,RecurrentOracle}, Tuple{InitialOracle,RecurrentOracle}}; kwargs...)
  reqc = launch_inference_server(os[1], os[2]; kwargs...)
  make = make_make(reqc)
  make2 = make_make2(reqc)
  return zipthunk(make,make2), send_done!(reqc)
end

function batchify_oracles(os::Tuple{<:Any, <:Any}; kwargs...)
  return zipthunk(ret_oracle(os[1]), ret_oracle(os[2])), do_nothing!
end

function batchify_oracles(os::Tuple{<:Any, Tuple{InitialOracle,RecurrentOracle}}; kwargs...)
  reqc = launch_inference_server(os[2]; kwargs...)
  make = make_make(reqc)
  return zipthunk(ret_oracle(os[1]),make), send_done!(reqc)
end

function batchify_oracles(os::Tuple{Tuple{InitialOracle,RecurrentOracle}, <:Any}; kwargs...)
  reqc = launch_inference_server(os[1]; kwargs...)
  make = make_make(reqc)
  return zipthunk(ret_oracle(os[2]),make), send_done!(reqc)
end


text:
#####
##### NUTS tree sampler implementation.
#####

####
#### Trajectory and implementation
####

"""
Representation of a trajectory, ie a Hamiltonian with a discrete integrator that
also checks for divergence.
"""
struct TrajectoryNUTS{TH,Tf,S}
    "Hamiltonian."
    H::TH
    "Log density of z (negative log energy) at initial point."
    π₀::Tf
    "Stepsize for leapfrog."
    ϵ::Tf
    "Smallest decrease allowed in the log density."
    min_Δ::Tf
    "Turn statistic configuration."
    turn_statistic_configuration::S
end

function move(trajectory::TrajectoryNUTS, z, fwd)
    @unpack H, ϵ = trajectory
    leapfrog(H, z, fwd ? ϵ : -ϵ)
end

###
### proposals
###

"""
$(SIGNATURES)

Random boolean which is `true` with the given probability `exp(logprob)`, which can be `≥ 1`
in which case no random value is drawn.
"""
function rand_bool_logprob(rng::AbstractRNG, logprob)
    logprob ≥ 0 || (randexp(rng, Float64) > -logprob)
end

function calculate_logprob2(::TrajectoryNUTS, is_doubling, ω₁, ω₂, ω)
    biased_progressive_logprob2(is_doubling, ω₁, ω₂, ω)
end

function combine_proposals(rng, ::TrajectoryNUTS, z₁, z₂, logprob2::Real, is_forward)
    rand_bool_logprob(rng, logprob2) ? z₂ : z₁
end

###
### statistics for visited nodes
###

struct AcceptanceStatistic{T}
    """
    Logarithm of the sum of metropolis acceptances probabilities over the whole trajectory
    (including invalid parts).
    """
    log_sum_α::T
    "Total number of leapfrog steps."
    steps::Int
end

function combine_acceptance_statistics(A::AcceptanceStatistic, B::AcceptanceStatistic)
    AcceptanceStatistic(logaddexp(A.log_sum_α, B.log_sum_α), A.steps + B.steps)
end

"""
$(SIGNATURES)

Acceptance statistic for a leaf. The initial leaf is considered not to be visited.
"""
function leaf_acceptance_statistic(Δ, is_initial)
    is_initial ? AcceptanceStatistic(oftype(Δ, -Inf), 0) : AcceptanceStatistic(min(Δ, 0), 1)
end

"""
$(SIGNATURES)

Return the acceptance rate (a `Real` betwen `0` and `1`).
"""
acceptance_rate(A::AcceptanceStatistic) = min(exp(A.log_sum_α) / A.steps, 1)

combine_visited_statistics(::TrajectoryNUTS, v, w) = combine_acceptance_statistics(v, w)

###
### turn analysis
###

"Statistics for the identification of turning points. See Betancourt (2017, appendix)."
struct GeneralizedTurnStatistic{T}
    p♯₋::T
    p♯₊::T
    ρ::T
end

function leaf_turn_statistic(::Val{:generalized}, H, z)
    p♯ = calculate_p♯(H, z)
    GeneralizedTurnStatistic(p♯, p♯, z.p)
end

function combine_turn_statistics(::TrajectoryNUTS,
                                 x::GeneralizedTurnStatistic, y::GeneralizedTurnStatistic)
    GeneralizedTurnStatistic(x.p♯₋, y.p♯₊, x.ρ + y.ρ)
end

function is_turning(::TrajectoryNUTS, τ::GeneralizedTurnStatistic)
    # Uses the generalized NUTS criterion from Betancourt (2017).
    @unpack p♯₋, p♯₊, ρ = τ
    @argcheck p♯₋ ≢ p♯₊ "internal error: is_turning called on a leaf"
    dot(p♯₋, ρ) < 0 || dot(p♯₊, ρ) < 0
end

###
### leafs
###

function leaf(trajectory::TrajectoryNUTS, z, is_initial)
    @unpack H, π₀, min_Δ, turn_statistic_configuration = trajectory
    Δ = is_initial ? zero(π₀) : logdensity(H, z) - π₀
    isdiv = Δ < min_Δ
    v = leaf_acceptance_statistic(Δ, is_initial)
    if isdiv
        nothing, v
    else
        τ = leaf_turn_statistic(turn_statistic_configuration, H, z)
        (z, Δ, τ), v
    end
end

####
#### NUTS interface
####

"Default maximum depth for trees."
const DEFAULT_MAX_TREE_DEPTH = 10

"""
$(TYPEDEF)

Implementation of the “No-U-Turn Sampler” of Hoffman and Gelman (2014), including subsequent
developments, as described in Betancourt (2017).

# Fields

$(FIELDS)
"""
struct NUTS{S}
    "Maximum tree depth."
    max_depth::Int
    "Threshold for negative energy relative to starting point that indicates divergence."
    min_Δ::Float64
    """
    Turn statistic configuration. Currently only `Val(:generalized)` (the default) is
    supported.
    """
    turn_statistic_configuration::S
    function NUTS(; max_depth = DEFAULT_MAX_TREE_DEPTH, min_Δ = -1000.0,
                  turn_statistic_configuration = Val{:generalized}())
        @argcheck 0 < max_depth ≤ MAX_DIRECTIONS_DEPTH
        @argcheck min_Δ < 0
        S = typeof(turn_statistic_configuration)
        new{S}(Int(max_depth), Float64(min_Δ), turn_statistic_configuration)
    end
end

"""
$(TYPEDEF)

Diagnostic information for a single tree built with the No-U-turn sampler.

# Fields

Accessing fields directly is part of the API.

$(FIELDS)
"""
struct TreeStatisticsNUTS
    "Log density (negative energy)."
    π::Float64
    "Depth of the tree."
    depth::Int
    "Reason for termination. See [`InvalidTree`](@ref) and [`REACHED_MAX_DEPTH`](@ref)."
    termination::InvalidTree
    "Acceptance rate statistic."
    acceptance_rate::Float64
    "Number of leapfrog steps evaluated."
    steps::Int
    "Directions for tree doubling (useful for debugging)."
    directions::Directions
end

"""
$(SIGNATURES)

No-U-turn Hamiltonian Monte Carlo transition, using Hamiltonian `H`, starting at evaluated
log density position `Q`, using stepsize `ϵ`. Parameters of `algorithm` govern the details
of tree construction.

Return two values, the new evaluated log density position, and tree statistics.
"""
function sample_tree(rng, algorithm::NUTS, H::Hamiltonian, Q::EvaluatedLogDensity, ϵ;
                          p = rand_p(rng, H.κ), directions = rand(rng, Directions))
    @unpack max_depth, min_Δ, turn_statistic_configuration = algorithm
    z = PhasePoint(Q, p)
    trajectory = TrajectoryNUTS(H, logdensity(H, z), ϵ, min_Δ, turn_statistic_configuration)
    ζ, v, termination, depth = sample_trajectory(rng, trajectory, z, max_depth, directions)
    tree_statistics = TreeStatisticsNUTS(logdensity(H, ζ), depth, termination,
                                         acceptance_rate(v), v.steps, directions)
    ζ.Q, tree_statistics
end


text:
function interpolate_state(mdp::CarMDP, s::CarMDPState)
    # interpolate s in the grid
    vspace = get_car_vspace(mdp.env, mdp.vel_res)
    itp_car, itp_car_w = interpolate_state(mdp, s.car, vspace)
    itp_ego, itp_ego_w = interpolate_state(mdp, s.ego, vspace)
    itp_states = Vector{CarMDPState}(length(itp_ego)*length(itp_car))
    itp_w = Vector{Float64}(length(itp_states))
    k = 1
    for (i, c) in enumerate(itp_car)
        for (j, e) in enumerate(itp_ego)
            crash = is_colliding(Vehicle(c, mdp.car_type, 0), Vehicle(e, mdp.ego_type, 1))
            itp_states[k] = CarMDPState(crash, e, c, s.route)
            itp_w[k] = itp_car_w[i]*itp_ego_w[j]
            k += 1
        end
    end
    @assert sum(itp_w) ≈ 1.
    return itp_states, itp_w
end

function get_mdp_state(mdp::CarMDP, pomdp::UrbanPOMDP, s::Scene, car_id = 2)
    car_i = findfirst(car_id, s)
    car = Vehicle(get_off_the_grid(mdp), mdp.car_type, car_id)
    if car_i != nothing
        car = s[car_i]
    end
    ego = get_ego(s)
    sroute = SVector{0, Lane}()
    # find route
    sroute = nothing
    if haskey(pomdp.models, car_id) && car_i != nothing
        # find the exact route from the list of routes
        curr_route = [l.tag for l in pomdp.models[car_id].navigator.route]
        for route in get_car_routes(mdp.env)
            tags = intersect(Set(curr_route), Set(route))
            if length(tags) >= 2
                sroute = SVector{2, LaneTag}(route[1], route[end])
            elseif length(curr_route) == 1 && curr_route[1] ∈ route
                sroute = SVector{2, LaneTag}(route[1], route[end])
            end
        end
        if sroute == nothing
            println(curr_route)
        end
    else
        sroute = SVector{2, LaneTag}(LaneTag(0,0), LaneTag(0,0))
    end
    e_state = VehicleState(ego.state.posG, car_roadway(mdp.env), ego.state.v)
    c_state = VehicleState(car.state.posG, mdp.env.roadway, car.state.v)
    return CarMDPState(is_colliding(ego, car), e_state, c_state, sroute)
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
"""
    Masked{K<:Kernel} <: Kernel

A wrapper for kernels so that they are only applied along certain dimensions.

This is similar to the `active_dims` kernel attribute in the python GPy package and to the
`covMask` function in the matlab gpml package.

The implementation is very simple: any function of the kernel that takes an `X::Matrix`
input is delegated to the wrapped kernel along with a view of `X` that only includes the
active dimensions.
"""
struct Masked{K<:Kernel, DIM} <: Kernel
    kernel::K
    active_dims::SVector{DIM, Int}
end
num_dims(masked::Masked{K, DIM}) where {K, DIM} = DIM
function Masked(kern::Kernel, active_dims::AbstractVector{Int})
    dim = length(active_dims)
    K = typeof(kern)
    return Masked{K,dim}(kern, SVector{dim, Int}(active_dims))
end

@inline function cov_ij(masked::Masked, X::MatF64, i::Int, j::Int, dim::Int)
    return cov_ij(masked.kernel, @view(X[masked.active_dims, :]), i, j, num_dims(masked))
end
@inline function cov_ij(masked::Masked, X::MatF64, data::KernelData, i::Int, j::Int, dim::Int)
    return cov_ij(masked.kernel, @view(X[masked.active_dims, :]), data, i, j, num_dims(masked))
end
@inline function cov_ij(masked::Masked, X::MatF64, data::EmptyData, i::Int, j::Int, dim::Int)
    return cov_ij(masked, X, i, j, dim)
end

@inline function dKij_dθp(masked::Masked, X::MatF64, data::KernelData, i::Int, j::Int, p::Int, dim::Int)
    return dKij_dθp(masked.kernel, @view(X[masked.active_dims, :]), data, i, j, p, num_dims(masked))
end
@inline function dKij_dθp(masked::Masked, X::MatF64, i::Int, j::Int, p::Int, dim::Int)
    return dKij_dθp(masked.kernel, @view(X[masked.active_dims, :]), i, j, p, num_dims(masked))
end
@inline function dKij_dθp(masked::Masked, X::MatF64, data::EmptyData, i::Int, j::Int, p::Int, dim::Int)
    return dKij_dθp(masked, X, i, j, p, dim)
end
@inline @inbounds function dKij_dθ!(dK::VecF64, masked::Masked, X::MatF64, data::KernelData, i::Int, j::Int, dim::Int, npars::Int)
    Xmasked = @view(X[masked.active_dims,:])
    return dKij_dθ!(dK, masked.kernel, Xmasked, data, i, j, num_dims(masked), npars)
end

function Statistics.cov(masked::Masked, x1::MatF64, x2::MatF64)
    return cov(masked.kernel, view(x1,masked.active_dims,:), view(x2,masked.active_dims,:))
end
function Statistics.cov(masked::Masked, x1::VecF64, x2::VecF64)
    return cov(masked.kernel, view(x1,masked.active_dims), view(x2,masked.active_dims))
end
function Statistics.cov(masked::Masked, X::MatF64, data::KernelData)
    return cov(masked.kernel, view(X,masked.active_dims,:), data)
end
function cov!(
    s::MatF64, masked::Masked, X::MatF64, data::KernelData)
    return cov!(s, masked.kernel, view(X,masked.active_dims,:), data)
end
function KernelData(masked::Masked, X::MatF64)
    return KernelData(masked.kernel, view(X,masked.active_dims,:))
end
function kernel_data_key(masked::Masked, X::MatF64)
    k = kernel_data_key(masked.kernel, view(X,masked.active_dims,:))
    return @sprintf("%s_active=%s", k, masked.active_dims)
end

get_params(masked::Masked) = get_params(masked.kernel)
get_param_names(masked::Masked) = get_param_names(masked.kernel)
num_params(masked::Masked) = num_params(masked.kernel)
set_params!(masked::Masked, hyp) = set_params!(masked.kernel, hyp)
function grad_slice!(
    dK::MatF64, masked::Masked, X::MatF64, data::KernelData, iparam::Int)
    return grad_slice!(dK, masked.kernel, view(X,masked.active_dims,:), data, iparam)
end

# with EmptyData
function Statistics.cov(masked::Masked, X::MatF64, data::EmptyData)
    return cov(masked.kernel, view(X,masked.active_dims,:), data)
end
function cov!(s::MatF64, masked::Masked, X::MatF64, data::EmptyData)
    return cov!(s, masked.kernel, view(X,masked.active_dims,:), data)
end
function grad_slice!(dK::MatF64, masked::Masked, X::MatF64, data::EmptyData, iparam::Int)
    return grad_slice!(dK, masked.kernel, view(X,masked.active_dims,:), data, iparam)
end

# priors
get_priors(masked::Masked) = get_priors(masked.kernel)
set_priors!(masked::Masked, priors) = set_priors!(masked.kernel, priors)


text:
using MultivariateStats
using LinearAlgebra
using Test
using StableRNGs
using Statistics: mean, cov
using StatsBase: ConvergenceException

@testset "ICA" begin

    rng = StableRNG(15678)

    function generatetestdata(rng, n, k, m)
        t = range(0.0, step=10.0, length=n)
        s1 = sin.(t * 2)
        s2 = s2 = 1.0 .- 2.0 * Bool[isodd(floor(Int, x / 3)) for x in t]
        s3 = Float64[mod(x, 5.0) for x in t]

        s1 += 0.1 * randn(rng, n)
        s2 += 0.1 * randn(rng, n)
        s3 += 0.1 * randn(rng, n)

        S = hcat(s1, s2, s3)'
        @assert size(S) == (k, n)

        A = randn(rng, m, k)

        X = A * S
        mv = vec(mean(X, dims=2))
        @assert size(X) == (m, n)
        C = cov(X, dims=2)
        return X, mv, C
    end

    @testset "Auxiliary" begin

        f = MultivariateStats.Tanh(1.0)
        u, v = ones(1,1), zeros(1)
        MultivariateStats.update!(f, u, v)
        @test u[] ≈ 0.7615941559557649
        @test v[] ≈ 0.41997434161402614

        f = MultivariateStats.Tanh(1.0f0)
        u, v = ones(Float32,1,1), zeros(Float32,1)
        MultivariateStats.update!(f, u, v)
        @test u[] ≈ 0.7615942f0
        @test v[] ≈ 0.41997433f0

        f = MultivariateStats.Gaus()
        u, v = ones(1,1), ones(1)
        MultivariateStats.update!(f, u, v)
        @test u[] ≈ 0.6065306597126334
        @test v[] ≈ 0.0

        u, v = ones(Float32,1,1), ones(Float32,1)
        MultivariateStats.update!(f, u, v)
        @test u[] ≈ 0.60653067f0
        @test v[] ≈ 0.0f0

    end

    ## data
    @testset "Algorithm" begin

        # sources
        n = 1000
        k = 3
        m = 8
        X, μ, C = generatetestdata(rng, n, k, m)

        # FastICA

        M = fit(ICA, X, k; do_whiten=false, tol=Inf)
        @test isa(M, ICA)
        @test size(M) == (m,k)
        @test mean(M) == μ
        W = M.W
        @test predict(M, X) ≈ W' * (X .- μ)
        @test W'W ≈ Matrix(I, k, k)

        M = fit(ICA, X, k; do_whiten=true, tol=Inf)
        @test isa(M, ICA)
        @test size(M) == (m,k)
        @test mean(M) == μ
        W = M.W
        @test W'C * W ≈ Matrix(I, k, k)

        @test_throws ConvergenceException fit(ICA, X, k; do_whiten=true, tol=1e-8, maxiter=2)

        # Use data of different type
        XX = convert(Matrix{Float32}, X)

        MM = fit(ICA, XX, k; do_whiten=true, tol=Inf)
        @test eltype(mean(MM)) == Float32
        @test eltype(MM.W) == Float32

        MM = fit(ICA, XX, k; do_whiten=false, tol=Inf)
        @test isa(MM, ICA)
        @test eltype(mean(MM)) == Float32
        @test eltype(MM.W) == Float32
        W = MM.W
        @test predict(MM, X) ≈ W' * convert(Matrix{Float32}, (X .- μ))
        @test predict(M, XX) ≈ M.W' * (X .- μ) atol=1e-4
        @test W'W ≈ Matrix{Float32}(I, k, k)

        # input as view
        M = fit(ICA, view(XX, :, 1:400), k; do_whiten=true, tol=Inf)
        @test eltype(mean(MM)) == Float32
        @test eltype(MM.W) == Float32
    end

end


text:
@inline second(x) = @inbounds getindex(x, 2)
@inline third(x) = @inbounds getindex(x, 3)
@inline fourth(x) = @inbounds getindex(x, 4)
@inline fifth(x) = @inbounds getindex(x, 5)
@inline sixth(x) = @inbounds getindex(x, 6)
@inline seventh(x) = @inbounds getindex(x, 7)
@inline eighth(x) = @inbounds getindex(x, 8)
@inline ninth(x) = @inbounds getindex(x, 9)


text:
"""
    maketestdata(seed)

Return a simple example time series with the composing parts.
"""
function maketestdata(seed)
  Random.seed!(seed)
  ## simulate data of length....
  N_tim = 240
  NpY   = 24         # samples/year
  t     = 1:N_tim
  t     = t./NpY # your time vector

  # constant seasonal cycle
  A   = 2          # amplitude
  phi = 13 * pi/12   # initial phase
  S   = A .* cos.(2 .* pi .* t .+ phi)

  # generate a linear trend
  T = 0.1 .+ 0.2.*t
  @show T

  # some other oscillation
  a = 0.2
  b = 0.1
  C = 1 .+ a .*cos.(b*2*pi*t)

  # simple (coloured) noise
  φ = 0.3 # strengh of autocorrelation in noise
  E = randn(N_tim) .* 0.1

  for i = 2:N_tim
    E[i] = φ*E[i-1]+(1-φ)*E[i]
  end

  X = @. S * C + 2 * E + T
  t, X, [S, C, E, T]
end

"""
colominas2014_s()

First test data from Colominas 2014 et al.
http://dx.doi.org/10.1016/j.bspc.2014.06.009

"""
function colominas2014_s()
  n = 1:1000
  s1 = sinpi.(2 .* 0.255 .* (n .- 501))
  s1[1:500] .= 0
  s1[751:1000] .=0
  s2 = sinpi.(2 .* 0.065 .*(n .- 1))
  n, s1+s2, [s1,s2]
end

"""
    colominas2014_x()

Second test data from Colominas 2014 et al.
http://dx.doi.org/10.1016/j.bspc.2014.06.009
It is not exactly the same, because ϕ is not feasable since the arccos is only
defined for values between -1 and 1.
"""
function colominas2014_x()
  n = 1:1000
  fmax = 3//32
  fmin = 9//128
  ϕ = acos((fmax -fmin)/(3fmin+fmax))
  x1 = 3 .* exp.((.-((x .-500)./100).^2).*π) .* cospi.(5//8 .*(n.-1000))
  x2 = @. cospi((fmax+fmin)*(n - 1000)+(fmax-fmin)*500(sinpi(n/500)+ϕ-sin(ϕ)))
  x3 = exp.((.-((n.-1000)/200).^2) .*π) .* cospi.((7 .//128).*(n.-1000))
  n, x1+x2+x3, [x1,x2,x3]
end

"""
    fosso2014()

Make testdata from Fosso 2014 et al.
"""
function fosso(x=0:.001:2)
  x1 = @. 0.7 * sinpi(16 * x)
  x2 = @. 0.7 * sinpi(48 * x)
  x3 = @. 1.4 * sinpi(60 * x)
  x, x1 .+ x2 .+ x3, [x1, x2, x3]
end


text:
module MaxHelpingHandHorizontalRoomWrapController

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/HorizontalRoomWrapController" HorizontalRoomWrapController(x::Integer, y::Integer)

const placements = Ahorn.PlacementDict(
    "Horizontal Room Wrap Controller (max480's Helping Hand)" => Ahorn.EntityPlacement(
        HorizontalRoomWrapController
    )
)

function Ahorn.selection(entity::HorizontalRoomWrapController)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::HorizontalRoomWrapController, room::Maple.Room) = Ahorn.drawSprite(ctx, "ahorn/MaxHelpingHand/horizontal_room_wrap", 0, 0)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
function main(r::Robot)
    x=Ost
    y=Nord
    path = go_to_left_bottom_corner_and_get_path(r)
    while !isborder(r,x)
        mark_direction(r,y)
        if !isborder(r,x)
            move!(r,x)
        end
        side_vert=inverse(y)
    end
    mark_direction(r,y)
    go_to_left_bottom_corner(r)
    go_back(r,path)
end

function mark_direction(r::Robot,side::HorizonSide)
    putmarker!(r)
    while try_move(r,side) == true
        putmarker!(r)
    end
    while try_move(r,inverse(side))
        ;
    end
end

function go_to_left_bottom_corner_and_get_path(r::Robot) ::Array
    path=[]
    while !isborder(r,Sud) || !isborder(r,West)
        if !isborder(r,Sud)
            move!(r,Sud)
            push!(path,Nord)
        end
        if !isborder(r,West)
            move!(r,West)
            push!(path,Ost)
        end
    end
    return path
end

function mark_direction(r,side::HorizonSide)
    while !isborder(r,side)
        putmarker!(r)
        move!(r,side)
    end
    putmarker!(r)
end

inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2, 4))
next_clockwise(side::HorizonSide) = HorizonSide(mod(Int(side)-1,4))

function go_to_left_bottom_corner(r::Robot)
    while !isborder(r,West) || !isborder(r,Sud)
        if !isborder(r,West)
            move!(r,West)
        end
        if !isborder(r,Sud)
            move!(r,Sud)
        end
    end
end

function go_back(r,path::Array)
    n=length(path)
    while n>0
        move!(r,path[n])
        n=n-1
    end
end

function try_move(r::Robot, direct_side::HorizonSide)::Bool
    orthogonal_side = next_clockwise(direct_side)
    reverse_side = inverse(orthogonal_side)
    num_steps=0
    if !isborder(r,direct_side)
        move!(r,direct_side)
        result=true
    else
        while isborder(r,direct_side) == true
            if !isborder(r, orthogonal_side)
                move!(r, orthogonal_side)
                num_steps += 1
            else
                break
            end
        end
        if !isborder(r,direct_side)
            move!(r,direct_side)
            while isborder(r,reverse_side) == true
                move!(r,direct_side)
            end
            result = true
        else
            result = false
        end
        while num_steps>0
            num_steps=num_steps-1
            move!(r,reverse_side)
        end
    end
    return result
end

max_stars_repo_name: VladislavBukin/mirea-progs

max_issues_repo_name: VladislavBukin/mirea-progs

max_forks_repo_name: VladislavBukin/mirea-progs


text:
# License is MIT: LICENSE.md

using ModuleInterfaceTools

@api test StrLiterals

@static if V6_COMPAT
    eval_parse(s) = eval(parse(s))
else
    eval_parse(s) = Core.eval(@__MODULE__, Meta.parse(s))
end

const ErrorType = @static V6_COMPAT ? ParseError : LoadError

ts(io) = String(take!(io))

@testset "Unicode constants" begin
    @test f"\u{3}"     == "\x03"
    @test f"\u{f60}"   == "\u0f60"
    @test f"\u{2a7d}"  == "\u2a7d"
    @test f"\u{1f595}" == "\U0001f595"
end
@testset "Interpolation" begin
    scott = 123
    @test f"\(scott)" == "123"
end
@testset "\$ not interpolation" begin
    @test f"I have $10, $spj$" == "I have \$10, \$spj\$"
end
@testset "Valid quoted characters" begin
    @test f"\$" == "\$"
    @test f"\"" == "\""
    @test f"\'" == "\'"
    @test f"\\" == "\\"
    @test f"\0" == "\0"
    @test f"\a" == "\a"
    @test f"\b" == "\b"
    @test f"\e" == "\e"
    @test f"\f" == "\f"
    @test f"\n" == "\n"
    @test f"\r" == "\r"
    @test f"\t" == "\t"
    @test f"\v" == "\v"
end
@testset "Invalid quoted characters" begin
    for ch in "cdghijklmopqsuwxy"
        @test_throws ErrorType eval_parse("f\"\\$ch\"")
    end
    for ch in 'A':'Z'
        @test_throws ErrorType eval_parse("f\"\\$ch\"")
    end
end

@testset "Legacy mode only sequences" begin
    # Check for ones allowed in legacy mode
    for s in ("f\"\\x\"", "f\"\\x7f\"", "f\"\\u\"", "f\"\\U\"", "f\"\\U{123}\"")
        @test_throws ErrorType eval_parse(s)
    end
end

@testset "Legacy mode Hex constants" begin
    @test F"\x3"     == "\x03"
    @test F"\x7f"    == "\x7f"
    @test_throws ErrorType eval_parse("F\"\\x\"")
    @test_throws ErrorType eval_parse("F\"\\x!\"")
end

@testset "Legacy mode Unicode constants" begin
    @test F"\u3"     == "\x03"
    @test F"\uf60"   == "\u0f60"
    @test F"\u2a7d"  == "\u2a7d"
    @test F"\U1f595" == "\U0001f595"
    @test F"\u{1f595}" == "\U0001f595"
    @test_throws ErrorType eval_parse("F\"\\U\"")
    @test_throws ErrorType eval_parse("F\"\\U!\"")
    @test_throws ErrorType eval_parse("F\"\\U{123}\"")
end

@testset "Legacy mode valid quoted characters" begin
    @test f"\$" == "\$"
    @test f"\"" == "\""
    @test f"\'" == "\'"
    @test f"\\" == "\\"
    @test f"\0" == "\0"
    @test f"\a" == "\a"
    @test f"\b" == "\b"
    @test f"\e" == "\e"
    @test f"\f" == "\f"
    @test f"\n" == "\n"
    @test f"\r" == "\r"
    @test f"\t" == "\t"
    @test f"\v" == "\v"
end

@testset "Legacy mode \$ interpolation" begin
    scott = 123
    @test F"$scott" == "123"
    @test F"$(scott+1)" == "124"
end

@testset "Quoted \$ in Legacy mode" begin
    @test F"I have \$10, \$spj\$" == "I have \$10, \$spj\$"
end

@testset "Print macro support" begin
    io = IOBuffer()
    scott = 123
    pr"\(io)This is a test with \(scott)"
    @test ts(io) == "This is a test with 123"
end

@testset "escape, unescape" begin
    @test s_escape_string(f"' \" \\ \u{7f} \u{20ac} \u{1f596} \u{e0000}") ==
        "' \\\" \\\\ \\u{7f} € 🖖 \\u{e0000}"
    @test s_unescape_string(f"' \\\" \\\\ \\u{7f} € 🖖 \\u{e0000}") ==
        "' \" \\ \x7f € 🖖 \Ue0000"
    io = IOBuffer()
    s_print_escaped(io, f"' \" \\ \u{7f} \u{20ac} \u{1f596} \u{e0000}", "")
    @test ts(io) == f"' \" \\\\ \\u{7f} € 🖖 \\u{e0000}"
    io = IOBuffer()
    s_print_unescaped(io, f"' \" \\\\ \\u{7f} € 🖖 \\u{e0000}")
    @test ts(io) == "' \" \\ \x7f € 🖖 \Ue0000"
end


text:
#= Load these if running these tests in terminal
import Pkg; Pkg.activate("../")  # CytofRepFAM
using CytofRepFAM, Random, BSON, Test, Distributions, Flux
=#

printstyled("Test VB FAM on simulated data ...\n", color=:yellow)
@testset "Variational Inference" begin
  println("Test Variational Inference...")
  # Set random seed for generating data
  Random.seed!(2)

  # Generate Data
  L = Dict(0=>3, 1=>3)
  J = 20
  K = 5

  println("Simulate Data...")
  N = [8, 1, 2] * 1000

  I = length(N)
  mus=Dict(0 => -[1.0, 2.3, 3.5],
           1 => +[1.0, 2.0, 3.0])
  a_W=rand(K)*10
  a_eta=Dict(z => rand(L[z])*10 for z in 0:1)
  Z = CytofRepFAM.Model.genSimpleZ(J, K)
  @time dat = CytofRepFAM.Model.genData(J=J, N=N, K=K, L=L, Z=Z,
                                   beta=[-9.2, -2.3],
                                   sig2=[0.2, 0.1, 0.3],
                                   mus=mus,
                                   a_W=a_W,
                                   a_eta=a_eta,
                                   sortLambda=false,
                                   propMissingScale=0.7,
                                   eps=fill(.005, I))

  tau = .0001 # FIXME: why can't I do .001?
  use_stickbreak = false
  noisy_var = 10.0
  K_MCMC = 10
  L_MCMC = Dict(false=>5, true=>3)
  priors = CytofRepFAM.VB.Priors(K=K_MCMC, L=L_MCMC, use_stickbreak=use_stickbreak)
  mc = CytofRepFAM.Model.defaultConstants(CytofRepFAM.Model.Data(dat[:y]),
                                     K_MCMC, Dict{Int64,Int64}(L_MCMC),
                                     yQuantiles=[0.0, 0.25, 0.5],
                                     pBounds=[.05, .8, .05])

  beta = [mc.beta[:, i] for i in 1:I]
  c = CytofRepFAM.VB.Constants(I, N, J, K_MCMC, L_MCMC, tau, beta, use_stickbreak,
                          noisy_var, priors)
  y = dat[:y]

  println("fit model...")
  # niters = 20000
  # batchsize = 2000
  # to test compilation
  niters = 20
  batchsize = 200
  opt = ADAM(1e-2)
  out = CytofRepFAM.VB.fit(y=y, niters=niters, batchsize=batchsize, c=c, opt=opt,
                      seed=0, nsave=30)

  # Save results
  outdir = "results/vb"
  mkpath(outdir)
  @time BSON.bson("$(outdir)/vb-out.bson", out)

  #= Load results
  using CytofRepFAM, Flux, Distributions, BSON
  =#
  out = BSON.load("$(outdir)/vb-out.bson")
  samples = [CytofRepFAM.VB.rsample(out[:state])[2] for b in 1:100]

  @test true
end


text:
include("DenseSync.jl")
using FFTW
using LinearAlgebra
using Plots
gr()

function rungekutta(func, θ0, δt, tmax, rec_func)
    ans = []
    for _ in 1:tmax÷δt
        k1 = func(θ0)
        k2 = func(θ0 + 0.5*δt*k1)
        k3 = func(θ0 + 0.5*δt*k2)
        k4 = func(θ0 + δt*k3)
        θ0 += δt*(k1 + 2k2 + 2k3 + k4)/6.0
        append!(ans, rec_func(θ0))
    end
    return ans
end

function DenseSyncODE(ds::DenseSync, δt, tmax)
    x = append!([0], copy(ds.x))
    y = fft(x)
    function ODEVec(θs)
        es = exp.(1im*θs)
        return imag.(conj.(es).*(ifft(y.*fft(es))))
    end
    N = ds.N
    θ_twisted = [2π*l*ds.p/N for l in 0:N-1]
    function TwistedNorm(θs)
        vec = mod2pi.(θs - θ_twisted)
        vec = [(x≥π ? 2π-x : x) for x in vec]
        return norm(vec)
    end
    # initial value
    σ = 0.1π
    noises = σ*randn(N)/√(N)
    noises = noises .- sum(noises)/N
    θ0 = θ_twisted + noises
    return rungekutta(ODEVec, θ0, δt, tmax, TwistedNorm)
end

N, p = 1900, 100
ds = exact_sol(N, p)

δt, tmax = 10^(-3), 4
ts = [i*δt for i in 1:tmax÷δt]
norms = [DenseSyncODE(ds, δt, tmax) for _ in 1:5]
plot(
    ts, [log.(norm) for norm in norms],
    xlims=(0, 4), label=false,
    xlabel="t", ylabel="log norm"
)
savefig("ode.png")

max_stars_repo_name: yonesuke/DenseSync

max_issues_repo_name: yonesuke/DenseSync

max_forks_repo_name: yonesuke/DenseSync


text:
function check(d, M, as)
    # as is sorted
    # check if all the cows can have a distance longer than d
    N = length(as)

    last = 1
    for i in 1:(M-1)
        crt = last + 1
        while crt ≤ N && as[crt] - as[last] < d
            crt += 1
        end

        if crt ≥ N
            return false
        else
            last = crt
        end
    end

    return true
end

function main()
    N, M = parse.(Int, split(readline()))
    as = [parse(Int, readline()) for _ in 1:N]
    sort!(as)

    d̲, d̄ = 0, maximum(as)
    while d̄ - d̲ > 1
        d = (d̄ + d̲) ÷ 2
        if check(d, M, as)
            d̲ = d
        else
            d̄ = d
        end
    end

    if check(d̄, M, as)
        println(d̄)
    else
        println(d̲)
    end
end

main()