text:
# ----------- #
# Blink setup #
# ----------- #

mutable struct ElectronDisplay <: AbstractPlotlyDisplay
    divid::Base.Random.UUID
    w::Nullable{Any}
end

const ElectronPlot = SyncPlot{ElectronDisplay}

ElectronDisplay(divid::Base.Random.UUID) = ElectronDisplay(divid, Nullable())
ElectronDisplay(p::Plot) = ElectronDisplay(p.divid)
ElectronPlot(p::Plot) = ElectronPlot(p, ElectronDisplay(p.divid))

fork(jp::ElectronPlot) = ElectronPlot(fork(jp.plot), ElectronDisplay())

# define some core methods on the display

"""
Return true if the display has been initialized and is still open, else false
"""
isactive(ed::ElectronDisplay) = isnull(ed.w) ? false : Blink.active(get(ed.w))

"""
If the display is active, close the window
"""
Base.close(ed::ElectronDisplay) = isactive(ed) && close(get(ed.w))

"Variable we can query to see if the plot has finished rendering"
function done_var(ed::ElectronDisplay)
    Symbol("done_$(split(string(ed.divid), "-")[1])")
end

"variable to store the svg data in"
function svg_var(ed::ElectronDisplay)
    Symbol("svg_$(split(string(ed.divid), "-")[1])")
end

"""
Return true if the display is active and the plot has finished rendering, else
false
"""
function Base.isready(ed::ElectronDisplay)
    !isactive(ed) && return false
    Blink.js(ed, done_var(ed))
end

# helper method to obtain a Blink window, given some options
get_window(opts::Dict) = Juno.isactive() ? Juno.Atom.blinkplot() : Window(opts)

"""
Initialize a Blink window for displaying the a plot. This will return an
existing window if it has been initialized and is still active, otherwise it
creates one.

Part of the creation process here is loading the plotly javascript.
"""
function get_window(ed::ElectronDisplay; kwargs...)
    if !isnull(ed.w) && active(get(ed.w))
        w = get(ed.w)
    else
        w = get_window(Dict(kwargs))
        ed.w = w
        # load the js here
        Blink.load!(w, _js_path)
        Blink.loadjs!(w, _mathjax_cdn_path)
        Blink.js(w, :($(done_var(ed)) = false))
        Blink.js(w, :($(svg_var(ed)) = undefined))

        # make sure plotly has been loaded
        done = false
        while !done
            try
                @js w Plotly
                done = true
            catch e
                continue
            end
        end
    end
    w
end

# the corresponding methods for ElectronPlot can (for the most part)
# just forward on to the view
Base.isready(p::ElectronPlot) = isready(p.view)

done_var(p::ElectronPlot) = done_var(p.view)
svg_var(p::ElectronPlot) = svg_var(p.view)

"""
Close the plot window (if active) and reset the fields and reset the `w` field
on the display to be `Nullable()`.
"""
function Base.close(p::ElectronPlot)
    close(p.view)
    p.view.w = Nullable()
    nothing
end

"""
Use the plot width and height to open a window of the correct size
"""
function get_window(p::ElectronPlot; kwargs...)
    w, h = size(p.plot)
    get_window(p.view; width=w, height=h, kwargs...)
end

Base.display(p::ElectronPlot; show = true, resize::Bool=_autoresize[1]) =
  display_blink(p; show = show, resize = resize)

function display_blink(p::ElectronPlot; show=true, resize::Bool=_autoresize[1],)
    w = get_window(p, show=show)

    done = done_var(p)
    svg = svg_var(p)

    if !resize
        code = """
        <script>
        (function(){
        var gd = Plotly.d3.select("body")
            .append("div")
            .attr("id", "$(p.plot.divid)")
            .node();
        var plot_json = $(json(p.plot));
        var data = plot_json.data;
        var layout = plot_json.layout;
        Plotly.newPlot(gd, data, layout).then(function(gd) {
                 $(done) = true;
                 return Plotly.toImage(gd, {"format": "svg"});
            }).then(function(data) {
                var svg_data = data.replace(/^data:image\\/svg\\+xml,/, "");
                $(svg) = decodeURIComponent(svg_data);
             });
         })();
        </script>
        """
        @js(w, Blink.fill("body", $code))
    else
        #= this is the same as above, with a few differences:

        - when we create the div, set some style attributes so the div itself
        grows with the display window
        - After obtaining the image data, we delete the width and height
        properties on the layout so that if the user happened to have set them
        on the Julia side, they will be removed -- allowing the plot to always
        fill the div
        =#
        magic = """
        <script>
        (function() {
            var WIDTH_IN_PERCENT_OF_PARENT = 100
            var HEIGHT_IN_PERCENT_OF_PARENT = 100;
            var gd = Plotly.d3.select('body')
                .append('div').attr("id", "$(p.plot.divid)")
                .style({
                    width: WIDTH_IN_PERCENT_OF_PARENT + '%',
                    'margin-left': (100 - WIDTH_IN_PERCENT_OF_PARENT) / 2 + '%',
                    height: HEIGHT_IN_PERCENT_OF_PARENT + 'vh',
                    'margin-top': (100 - HEIGHT_IN_PERCENT_OF_PARENT) / 2 + 'vh'
                })
                .node();
            var plot_json = $(json(p.plot));
            var data = plot_json.data;
            var layout = plot_json.layout;
            Plotly.newPlot(gd, data, layout).then(function(gd) {
                     $(done) = true;
                     var img_data = Plotly.toImage(gd, {"format": "svg"});
                     delete gd.layout.width
                     delete gd.layout.height
                     return img_data
                }).then(function(img_data) {
                    var svg_data = img_data.replace(/^data:image\\/svg\\+xml,/, "");
                    $(svg) = decodeURIComponent(svg_data);
                 });
            window.onresize = function() {
                Plotly.Plots.resize(gd);
                };
            }
        )();
        </script>
        """
        @js(w, Blink.fill("body", $magic))
    end
    p.plot
end

## API Methods for ElectronDisplay
function _img_data(p::ElectronPlot, format::String; show::Bool=false)
    _formats = ["png", "jpeg", "webp", "svg"]
    if !(format in _formats)
        error("Unsupported format $format, must be one of $_formats")
    end

    opened_here = !isactive(p.view)
    opened_here && display(p; show=show, resize=false)

    out = @js p.view begin
        ev = Plotly.Snapshot.toImage(this, d("format"=>$format))
        @new Promise(resolve -> ev.once("success", resolve))
    end
    opened_here && close(p)
    out
end

function svg_data(p::ElectronPlot, format="png", robust=false)
    opened_here = !isactive(p.view)
    opened_here && (display(p; show=false, resize=false); sleep(0.1))

    out = nothing
    while out === nothing
        out = @js p $(svg_var(p))
        # wait for plot to render
        sleep(0.2)
    end

    opened_here && close(p)

    return out
end

function Blink.js(p::ElectronDisplay, code::JSString; callback=true)
    if !isactive(p)
        return
    end

    Blink.js(get_window(p),
             :(Blink.evalwith(document.getElementById($(string(p.divid))),
                              $(Blink.jsstring(code)))),
             callback=callback)
end

Blink.js(p::ElectronPlot, code::JSString; callback=true) =
    Blink.js(p.view, code; callback=callback)

# Methods from javascript API (docstrings found in api.jl)
relayout!(p::ElectronDisplay, update::Associative=Dict(); kwargs...) =
    @js_ p Plotly.relayout(this, $(merge(update, prep_kwargs(kwargs))))

restyle!(p::ElectronDisplay, ind::Int, update::Associative=Dict(); kwargs...) =
    @js_ p Plotly.restyle(this, $(merge(update, prep_kwargs(kwargs))), $(ind-1))

restyle!(p::ElectronDisplay, inds::AbstractVector{Int}, update::Associative=Dict(); kwargs...) =
    @js_ p Plotly.restyle(this, $(merge(update, prep_kwargs(kwargs))), $(inds-1))

restyle!(p::ElectronDisplay, update=Dict(); kwargs...) =
    @js_ p Plotly.restyle(this, $(merge(update, prep_kwargs(kwargs))))

addtraces!(p::ElectronDisplay, traces::AbstractTrace...) =
    @js_ p Plotly.addTraces(this, $traces)

addtraces!(p::ElectronDisplay, where::Int, traces::AbstractTrace...) =
    @js_ p Plotly.addTraces(this, $traces, $(where-1))

deletetraces!(p::ElectronDisplay, traces::Int...) =
    @js_ p Plotly.deleteTraces(this, $(collect(traces)-1))

movetraces!(p::ElectronDisplay, to_end::Int...) =
    @js_ p Plotly.moveTraces(this, $(collect(to_end)-1))

movetraces!(p::ElectronDisplay, src::AbstractVector{Int}, dest::AbstractVector{Int}) =
    @js_ p Plotly.moveTraces(this, $(src-1), $(dest-1))

redraw!(p::ElectronDisplay) =
    @js_ p Plotly.redraw(this)

purge!(p::ElectronDisplay) =
    @js_ p Plotly.purge(this)

to_image(p::ElectronDisplay; kwargs...) =
    @js p Plotly.toImage(this, $(Dict(kwargs)))

download_image(p::ElectronDisplay; kwargs...) =
    @js p Plotly.downloadImage(this, $(Dict(kwargs)))

# unexported (by plotly.js) api methods
extendtraces!(ed::ElectronDisplay, update::Associative=Dict(),
              indices::Vector{Int}=[1], maxpoints=-1;) =
    @js_ ed Plotly.extendTraces(this, $(prep_kwargs(update)), $(indices-1), $maxpoints)

prependtraces!(ed::ElectronDisplay, update::Associative=Dict(),
               indices::Vector{Int}=[1], maxpoints=-1;) =
    @js_ ed Plotly.prependTraces(this, $(prep_kwargs(update)), $(indices-1), $maxpoints)


text:
module RandomSequence

using SPPL

ns = @sppl (debug) begin
    X = array(3)
    W = array(3)
    X[0] ~ SPPL.Normal()
    for i in 1 : 2
        if X[i - 1] > 0
            X[i] ~ SPPL.Normal(0, 1)
            W[i] ~ SPPL.Atomic(0)
        else
            W[i] ~ SPPL.Bernoulli(0.5)
            W[i] == 0 ? X[i] .> SPPL.Fraction(1, 2) * X[i-1]^2 + X[i-1] : X[i] ~ SPPL.Normal(0, 1)
        end
    end
end

random_sequence = ns.model

# Observe X[1] > 0.
random_sequence_given_X1 = condition(model, ns.X[1] > 0)
println(probability(random_sequence, ns.X[1] > 0))
println(probability(random_sequence_given_X1, ns.X[1] > 0))

# Observe x[2] > 0.
random_sequence_given_X2 = random_sequence.condition(ns.X[2] > 0)
println(probability(random_sequence, ns.X[1] > 0))
println(probability(random_sequence_given_X2, ns.X[1] > 0))

# Compute mutual information.
mutual_information(random_sequence, ns.X[0]>0, ns.X[1]>0) |> println
mutual_information(random_sequence, ns.X[0]>0, ns.X[2]>0) |> println

end # module

ext: jl


text:
module Parameters

abstract type AbstractRRTMGPParameters end
const ARP = AbstractRRTMGPParameters

Base.@kwdef struct RRTMGPParameters{FT} <: ARP
    grav::FT
    molmass_dryair::FT
    molmass_water::FT
    gas_constant::FT
    kappa_d::FT
    Stefan::FT
    avogad::FT
end

# Method wrappers
for var in fieldnames(RRTMGPParameters)
    @eval $var(ps::ARP) = ps.$var
end

# Derived parameters
R_d(ps::ARP) = gas_constant(ps) / molmass_dryair(ps)
cp_d(ps::ARP) = R_d(ps) / kappa_d(ps)

end # module

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
"""
    generate_output_from_url(choice::AbstractPolyChoice, url, params, sink)

Internal function to dynamically generate output from a user specified AbstractPolyChoice.
"""
function generate_output_from_url(choice::AbstractPolyChoice, url, params, sink)
    request_json = HTTP.get(url, query=params).body |> JSON3.read

    return apply_choice(choice, request_json, sink)
end

"""
    apply_choice(::NoSinkYesResults, x, sink)

Internal function to apply a choice for no sink and yes results scenarios.
"""
function apply_choice(::NoSinkYesResults, x, sink)
    return x.results
end

"""
    apply_choice(::NoSinkNoResults, x, sink)

Internal function to apply a choice for no sink and no results scenarios.
"""
function apply_choice(::NoSinkNoResults, x, sink)
    return x
end

"""
    apply_choice(::YesSinkYesResults, x, sink)

Internal function to apply a choice for yes sink and yes results scenarios.
"""
function apply_choice(::YesSinkYesResults, x, sink)
    try
        return x.results |> jsontable |> sink
    catch
        return x.results |> x -> sink([x])
    end
end

"""
    apply_choice(::YesSinkNoResults, x, sink)

Internal function to apply a choice for no sink and yes results scenarios.
"""
function apply_choice(::YesSinkYesResults, x, sink::Nothing)
    return apply_choice(NoSinkYesResults(), x, sink)
end

"""
    apply_choice(::YesSinkNoResults, x, sink)

Internal function to apply a choice for yes sink and no results scenarios.
"""
function apply_choice(::YesSinkNoResults, x, sink)
    try
        return x |> jsontable |> sink
    catch
        return x |> x -> sink([x])
    end
end

"""
    apply_choice(::YesSinkNoResults, x, sink)

Internal function to apply a choice for no sink and no results scenarios.
"""
function apply_choice(::YesSinkNoResults, x, sink::Nothing)
    return apply_choice(NoSinkNoResults(), x, sink)
end

"""
    apply_choice(::NoSinkYesTickers, x, sink)

Internal function to apply a choice for no sink and yes tickers scenarios.
"""
function apply_choice(::NoSinkYesTickers, x, sink)
    return x.tickers
end

"""
    apply_choice(::YesSinkYesTickers, x, sink)

Internal function to apply a choice for yes sink and yes tickers scenarios.
"""
function apply_choice(::YesSinkYesTickers, x, sink)
    return x.tickers |> jsontable |> sink
end

"""
    apply_choice(::YesSinkYesTickers, x, sink::Nothing)

Internal function to apply a choice for no sink and yes tickers scenarios.
"""
function apply_choice(::YesSinkYesTickers, x, sink::Nothing)
    return apply_choice(NoSinkYesTickers(), x, sink)
end

"""
    apply_choice(::NoSinkYesTicker, x, sink)

Internal function to apply a choice for no sink and yes tickers scenarios.
"""
function apply_choice(::NoSinkYesTicker, x, sink)
    return x.ticker
end

"""
    apply_choice(::YesSinkYesLast, x, sink)

Internal function to apply a choice for yes sink and yes last scenarios.
"""
function apply_choice(::YesSinkYesLast, x, sink)
    return x.last |> x -> sink([x])
end

"""
    apply_choice(::NoSinkYesLast, x, sink)

Internal function to apply a choice for no sink and yes last scenarios.
"""
function apply_choice(::NoSinkYesLast, x, sink)
    return x.last
end

"""
    apply_choice(::YesSinkYesLast, x, sink::Nothing)

Internal function to apply a choice for no sink and yes last scenarios.
"""
function apply_choice(::YesSinkYesLast, x, sink::Nothing)
    return apply_choice(NoSinkYesLast(), x, sink)
end

"""
    apply_choice(::YesSinkYesTicker, x, sink)

Internal function to apply a choice for yes sink and yes ticker scenarios.
"""
function apply_choice(::YesSinkYesTicker, x, sink)
    try
        return x.ticker |> jsontable |> sink
    catch
        return x.ticker |> x -> sink([x])
    end
end

"""
    apply_choice(::YesSinkYesTicker, x, sink::Nothing)

Internal function to apply a choice for yes sink and yes ticker scenarios.
"""
function apply_choice(::YesSinkYesTicker, x, sink::Nothing)
    return x.ticker
end

"""
    apply_choice(::NoSinkYesData, x, sink)

Internal function to apply a choice for no sink and yes data scenarios.
"""
function apply_choice(::NoSinkYesData, x, sink)
    return x.data
end

"""
    apply_choice(::YesSinkYesData, x, sink)

Internal function to apply a choice for yes sink and yes data scenarios.
"""
function apply_choice(::YesSinkYesData, x, sink)
    try
        return x.data |> jsontable |> sink
    catch
        return x.data |> x -> sink([x])
    end
end

"""
    apply_choice(::YesSinkYesData, x, sink::Nothing)

Internal function to apply a choice for no sink and yes data scenarios.
"""
function apply_choice(::YesSinkYesData, x, sink::Nothing)
    return apply_choice(NoSinkYesData(), x, sink)
end

"""
    apply_choice(::YesSinkYesTicks, x, sink)

Internal function to apply a choice for yes sink and yes ticks scenarios.
"""
function apply_choice(::YesSinkYesTicks, x, sink)
    return x.ticks |> jsontable |> sink
end

"""
    apply_choice(::NoSinkYesTicks, x, sink)

Internal function to apply a choice for no sink and yes ticks scenarios.
"""
function apply_choice(::NoSinkYesTicks, x, sink)
    return x.ticks
end

"""
    apply_choice(::YesSinkYesTicks, x, sink::Nothing)

Internal function to apply a choice for no sink and yes ticks scenarios.
"""
function apply_choice(::YesSinkYesTicks, x, sink::Nothing)
    return apply_choice(NoSinkYesTicks(), x, sink)
end


text:
datadir = "datasets/tr_to_en"

if !isdir(datadir)
    download("http://www.phontron.com/data/qi18naacl-dataset.tar.gz", "qi18naacl-dataset.tar.gz")
    run(`tar xzf qi18naacl-dataset.tar.gz`)
end

if !isdefined(Main, :tr_vocab)
    tr_vocab = Vocab("$datadir/tr.train", mincount=5)
    en_vocab = Vocab("$datadir/en.train", mincount=5)
    tr_train = TextReader("$datadir/tr.train", tr_vocab)
    en_train = TextReader("$datadir/en.train", en_vocab)
    tr_dev = TextReader("$datadir/tr.dev", tr_vocab)
    en_dev = TextReader("$datadir/en.dev", en_vocab)
    tr_test = TextReader("$datadir/tr.test", tr_vocab)
    en_test = TextReader("$datadir/en.test", en_vocab)
    dtrn = MTData(tr_train, en_train)
    ddev = MTData(tr_dev, en_dev)
    dtst = MTData(tr_test, en_test)
end

max_stars_repo_name: denizyuret/research

max_issues_repo_name: denizyuret/research

max_forks_repo_name: denizyuret/research


text:
using Test
using DynamicalBilliards
using DynamicalBilliards.Testing

function test_finite_collisions(p, bd, N = 1e6)
    n = 0
    i, t = next_collision(p, bd)
    s = Set(i)
    @testset "$(tag(p)) 0 < t < Inf" begin
        while n < N
            i, t = bounce!(p, bd)
            if 0 < t < Inf
                @test true
            else
                error("$(tag(p, bd)) collision time $(t) with obst.index $(i)")
            end
            n += 1
            push!(s, i)
        end
    end
    if typeof(p) <: Particle && !isperiodic(bd)
        iprev = 0
        n = 0
        @testset "collide with each obstacle ONCE!" begin
            while n < 1000
                i, t = bounce!(p, bd)
                if i != 4
                    if i != iprev
                        @test true
                    else
                        error("$(tag(p)) consecutive collisions with $(i)")
                    end
                end
                iprev = i
                n += 1
            end
        end
    end
    @testset "$(tag(p, bd)) collided with all obstacles" begin
        @test s == Set(1:length(bd))
    end
end

function test_maximum_distance(p, bd, thres = 1e-12, N = 1e6)
    acc = acclevel(p, bd)
    n = 0
    t, i = next_collision(p, bd)
    maxds = fill(-Inf, length(bd))
    minds = fill(Inf, length(bd))
    thresholds = fill(0, length(bd))
    isp = isperiodic(bd)
    @testset "$(tag(p, bd)) max/min distance" begin
        while n < N

            i, t, cp = next_collision(p, bd)
            d = distance(cp, bd[i])

            d < minds[i] && (minds[i] = d)
            d > maxds[i] && (maxds[i] = d)

            i, t = bounce!(p, bd)

            @test abs(d) < acc
            abs(d) > thres && (thresholds[i] += 1)

            n += 1
        end
        println(tag(p, bd))
        println("minds: ", minds)
        println("maxds: ", maxds)
        println("more than $(thres): ", thresholds, " (acc = $(acc))")
        println()
    end
end

billiards_testset("Finite Collision time", test_finite_collisions)
billiards_testset("min/max distance", test_maximum_distance, 1e-12)


text:
using Distributions, Plots

function visualize_gamma()
    a = 0.01:0.01:4
    b = 0.01:0.01:4
    A = repeat(a', length(b), 1)
    B = repeat(b, 1, length(a))
    Z = Gamma.(A, B)

    xx = range(0.01, 4, length=100)
    p = contour(a, b, pdf.(Z,0.01), title="Gamma PDF for x=0.01")

    anim = @animate for x in xx
        contour(a, b, pdf.(Z,x), title="Gamma PDF for x=$(round(x,digits=4))")
    end
    gif(anim, "tmp.gif", fps=10)
    try
        if "BROWSER" in keys(ENV)
            browser = ENV["BROWSER"]
            cmd = `open tmp.gif -a $browser`
            run(cmd)
        else
            browser = "$(ENV["HOME"])/Applications/Google Chrome.app/"
            cmd = `open tmp.gif -a $browser`
            run(cmd)
        end
    finally
        return anim
    end
end
visualize_gamma()

max_stars_repo_name: dysonance/Preferences

max_issues_repo_name: dysonance/Preferences

max_forks_repo_name: dysonance/Preferences


text:
module MockAWS

using Base64
using GitHub
using JSON
using Mocking
using OrderedCollections: OrderedDict

export @patches, generate_patches

const REPO_NAME = "aws/aws-sdk-js"

"""
Attach all patches to a high-level AWS.jl `service`.
"""
macro patches(service::Symbol)
    service_name = joinpath(@__DIR__, "patches", lowercase(string(service)) * ".jl")

    # Main is the toplevel module
    # This expects to mirror where AWS.@service places its modules
    return Expr(:toplevel, :(Main.$(service).include($service_name)))
end

"""
Get all aws/aws-sdk-js files
"""
function _get_aws_sdk_js_files(auth::GitHub.OAuth2)
    master_tree = tree(REPO_NAME, "master"; auth=auth)
    apis_sha = [t for t in master_tree.tree if t["path"]=="apis"][1]["sha"]
    files = tree(REPO_NAME, apis_sha)

    return files.tree
end

"""
Filter only examples.json files
"""
function _get_example_files(files::Vector, auth::GitHub.OAuth2)
    example_files = filter(f -> endswith(f["path"], ".examples.json"), files)

    seen_services = Set{String}()
    latest_versions = OrderedDict[]

    for service in reverse(example_files)
        service_name, _ = @mock _get_service_and_version(service["path"], files, auth)

        if !(service_name in seen_services)
            push!(seen_services, service_name)
            push!(latest_versions, service)
        end
    end

    return latest_versions
end

"""
Find the SHA hash for a give filename
"""
function _get_sha(filename::String, files::Vector)
    for file in files
        if file["path"] == filename
            return file["sha"]
        end
    end
end

"""
Get service name and version from normal.json file
"""
function _get_service_and_version(filename::String, files::Vector, auth::GitHub.OAuth2)
    try
        filename = replace(filename, "examples"=>"normal")
        sha = _get_sha(filename, files)
        contents = blob(REPO_NAME, sha; auth=auth)
        contents = JSON.parse(String(base64decode(contents.content)))

        metadata = contents["metadata"]
        service = lowercase(metadata["serviceId"])
        service = replace(service, " " => '_')
        version = metadata["apiVersion"]

        return (service, version)
    catch e
        if e isa BoundsError
            throw(InvalidFileName("$filename is an invalid AWS JSON filename."))
        else
            rethrow()
        end
    end
end

"""
Format function names to match Julia styling
"""
function _format_name(function_name::AbstractString)
    # Replace a string of uppercase characters with themselves prefaced by an underscore
    # [A-Z](?![A-Z]) => Match a single uppercase character that is not followed by another uppercase character
    # |(A-Z]{1,})    => Match 1-Infinite amounts of uppercase characters
    function_name = replace(function_name, r"[A-Z](?![A-Z])|([A-Z]{1,})" => s"_\g<0>")

    # Lowercase the entire string
    function_name = lowercase(function_name)

    # Chop off the leading underscore
    return if startswith(function_name, "_")
        chop(function_name, head=1, tail=0)
    else
        function_name
    end
end

"""
Generate patches for endpoints
"""
function _generate_endpoint_patches(endpoint_name::String, endpoint_examples::Array)
    patches = String[]

    for (i, example) in enumerate(endpoint_examples)
        if haskey(example, "output")
            endpoint_name = _format_name(endpoint_name)

            patch = """
            $(endpoint_name)_$(i)_patch = @patch function $(endpoint_name)(a...; b...)
                return $(example["output"])
            end
            """

            push!(patches, patch)
        end
    end

    return patches
end

"""
Generate the patch file
"""
function _create_service_patches(service_name::String, service::Dict{String, <:Any})
    filename = joinpath(@__DIR__, "patches", "$(service_name).jl")
    service_examples = service["examples"]

    endpoints = collect(keys(service_examples))
    patches = String[]

    for endpoint in endpoints
        endpoint_patches = _generate_endpoint_patches(endpoint, service_examples[endpoint])

        patches = vcat(patches, endpoint_patches)
        sort!(patches)
    end

    if !isempty(patches)
        template = """
        using Mocking

        $(join(patches, "\n"))
        """

        open(filename, "w") do f
            print(f, template)
        end
    end
end

"""
Generate patches for various AWS services based off of the `aws/aws-sdk-js` repository.
"""
function generate_patches()
    auth = GitHub.authenticate(ENV["GITHUB_AUTH"])

    files = _get_aws_sdk_js_files(auth)
    example_files = _get_example_files(files, auth)

    for file in example_files
        service_name, _ = _get_service_and_version(file["path"], files, auth)
        contents = blob(REPO_NAME, file["sha"]; auth=auth)
        contents = JSON.parse(String(base64decode(contents.content)))

        @info "Generating patches for $(service_name)"

        _create_service_patches(service_name, contents)
    end
end

end


text:
function trinary_to_decimal(str::AbstractString)

end


text:
import EarlyStopping: done, update, update_training

# dummy criterion: stop when loss = training loss
struct Dummy <: StoppingCriterion end

@test update_training(Dummy(), 1.0) === nothing
@test update_training(Dummy(), 1.0, 42.0) == 42.0

update_training(c::Dummy, loss) = (training=loss, loss=nothing)
update_training(c::Dummy, loss, state) = (training=loss, loss=state.loss)
update(c::Dummy, loss, state) = (training=state.training, loss=loss)
done(c::Dummy, state) = state.training == state.loss

stopper  = EarlyStopper(Dummy())
@test !done!(stopper, 1.0, training=true)
@test done!(stopper, 1.0)

true


text:
"""
    AbstractClusteringMethod{P}

Abstract supertype for all clustering types, with partition of type `P`.

### Notes

A clustering method defines a function which maps reach-sets to one or several
reach-sets. The mapping can possibly be over-approximative, i.e. such that the
set union of the input reach-sets is included in the set union of the output reach-sets.

By taking the convex hull of the input reach-sets one can reduce the number of
outputs sets to a single one, overapproximately. This is the method that corresponds
to the `LazyClustering` type. However, in some cases it is conveinent to do other
types of transformations, such as:

- Return several reach-sets that are obtained by grouping the input in a way
  defined by a partition. For example, a vector of ten input sets can be clustered
  in two groups of five sets each, or five groups of two sets each, etc.
- Use different set representations such as boxes or zonotopes.
- Further post-process the output of the convexification by splitting into smaller
  sets, eg. box splitting.

Each concrete subtype of `AbstractClusteringMethod` has a parameter `P` that defines
what type of clustering strategy is applied. The method should be accessed with the
`partition` getter function.

The following strategies are implemented
at the interface level:

- If `P` is of type `Missing`: no partition is applied
- If `P` is of type integer: the partition corresponds to gruping the into the given integer
                             number of sets (or as close as possible)
- If `P` is of type vector of vectors: the given partition is applied

### Examples

`LazyClustering([1:2, 3:10])` groups the first two reach-sets in one cluster and
the third to tenth reach-sets in another cluster.
"""
abstract type AbstractClusteringMethod{P} end

_partition(method::AbstractClusteringMethod{Missing}, idx) = nfolds(idx, 1)
_partition(method::AbstractClusteringMethod{<:Integer}, idx) = nfolds(idx, partition(method))
_partition(method::AbstractClusteringMethod{VT}, idx) where {D<:Integer, VTi<:AbstractVector{D}, VT<:AbstractVector{VTi}} = [idx[vi] for vi in partition(method)]

# partition an array x into n parts of equal size (if possible); if n is too large
# we split into length(x) parts
function nfolds(x::AbstractVector, n::Int)
    p = length(x)
    m = min(n, p)
    s = p / m
    [x[round(Int64, (i-1)*s)+1:min(p,round(Int64, i*s))] for i in 1:m]
end

# =====================================
# Identity (no clustering)
# =====================================

"""
    NoClustering{P} <: AbstractClusteringMethod{P}

No-op clustering method.
"""
struct NoClustering{P} <: AbstractClusteringMethod{P}
#
end

partition(::NoClustering{Missing}) = missing

NoClustering() = NoClustering{Missing}()

function cluster(F, idx, ::NoClustering)
    return view(F, idx)
end

# =====================================
# Lazy convexification clustering
# =====================================

"""
    LazyClustering{P} <: AbstractClusteringMethod{P}

Cluster according to the given partition by applying a lazy representation of the
convex hull.
"""
struct LazyClustering{P, T} <: AbstractClusteringMethod{P}
    partition::P
    convex::T
end

partition(method::LazyClustering) = method.partition

LazyClustering(; convex::Bool=true) = LazyClustering(missing, Val(convex))
LazyClustering(nchunks::D; convex::Bool=true) where {D<:Integer} = LazyClustering{D, typeof(Val(convex))}(nchunks, Val(convex))
LazyClustering(partition::VT) where {D<:Integer, VTi<:AbstractVector{D}, VT<:AbstractVector{VTi}} = LazyClustering{VT, typeof(Val(convex))}(partition, Val(convex))

function cluster(F, idx, ::LazyClustering{Missing, Val{true}})
    return [convexify(view(F, idx))]
end

# FIXME return a UnionSetArray of reach-sets
function cluster(F, idx, ::LazyClustering{Missing, Val{false}})
    return [view(F, idx)]
end

function cluster(F, idx, method::LazyClustering{P, Val{true}}) where {P}
    p = _partition(method, idx)
    convF = [convexify(view(F, cj)) for cj in p]
end

function cluster(F, idx, method::LazyClustering{P, Val{false}}) where {P}
    p = _partition(method, idx)
    convF = [view(F, cj) for cj in p]
end

# for Taylor model flowpipes we preprocess it with a zonotopic overapproximation
function cluster(F::Flowpipe{N, TaylorModelReachSet{N}}, idx, method::LazyClustering{P, Val{true}}) where {N, P}
    Fz = overapproximate(Flowpipe(view(F, idx)), Zonotope)
    return cluster(Fz, 1:length(idx), method) # Fx is now indexed from 1 ... length(idx)
end

# ambiguity fix
function cluster(F::Flowpipe{N, TaylorModelReachSet{N}}, idx, method::LazyClustering{P, Val{false}}) where {N, P}
    Fz = overapproximate(Flowpipe(view(F, idx)), Zonotope)
    return cluster(Fz, 1:length(idx), method)
end

# =====================================
# Box clustering
# =====================================

"""
    BoxClustering{PI, PO} <: AbstractClusteringMethod{P}

### Notes

This method first takes a lazy convex hull for the given partition, then computes
a tight hyperrectangular approximation for each element in the partition.
"""
struct BoxClustering{PI, PO} <: AbstractClusteringMethod{PI}
    inpartition::PI
    outpartition::PO
end

partition(method::BoxClustering) = method.inpartition

BoxClustering() = BoxClustering(missing, missing)
BoxClustering(nchunks::D) where {D<:Integer} = BoxClustering{D, Missing}(nchunks, missing)
BoxClustering(partition::VT) where {D<:Integer, VTi<:AbstractVector{D}, VT<:AbstractVector{VTi}} = BoxClustering{VT, Missing}(partition, missing)

_out_partition(method::BoxClustering{PI, Missing}, idx, n) where {PI} = fill(1, n)
_out_partition(method::BoxClustering{PI, <:Integer}, idx, n) where {PI} = fill(method.outpartition, n)
_out_partition(method::BoxClustering{PI, VT}, idx, n) where {PI, D<:Integer, VT<:AbstractVector{D}} = method.outpartition

# do not partition neither the input nor the output
function cluster(F::Flowpipe{N, <:AbstractLazyReachSet}, idx, ::BoxClustering{Missing, Missing}) where {N}
    convF = convexify(view(F, idx))
    return [overapproximate(convF, Hyperrectangle)]
end

# do not partition neither the input but split the output
function cluster(F::Flowpipe{N, <:AbstractLazyReachSet}, idx, method::BoxClustering{Missing, PO}) where {N, PO}
    convF = convexify(view(F, idx))
    Y = overapproximate(convF, Hyperrectangle)
    p = _out_partition(method, idx, dim(F))
    return split(Y, p)
end

# partition the input array and do not partition the output
#function cluster(F::Flowpipe{N, Union{<:ReachSet, <:SparseReachSet}}, idx, method::BoxClustering{PI, Missing}) where {N, PI}
function cluster(F::Flowpipe{N, <:AbstractLazyReachSet}, idx, method::BoxClustering{PI, Missing}) where {N, PI}
    p = _partition(method, idx)
    convF = [convexify(view(F, cj)) for cj in p]
    return [overapproximate(Xc, Hyperrectangle) for Xc in convF]
end

# partition the input array and do not partition the output
function cluster(F::Flowpipe{N, <:AbstractLazyReachSet}, idx, method::BoxClustering{PI, PO}) where {N, PI, PO}
    p = _partition(method, idx)
    convF = [convexify(view(F, cj)) for cj in p]
    Y = [overapproximate(Xc, Hyperrectangle) for Xc in convF]
    pout = _out_partition(method, idx, dim(F))
    Ys = reduce(vcat, [split(y, pout) for y in Y])
    return Ys
end

# for Taylor model flowpipes we preprocess it with a zonotopic overapproximation
function cluster(F::Flowpipe{N, TaylorModelReachSet{N}}, idx, method::BoxClustering) where {N}
    Fz = overapproximate(Flowpipe(view(F, idx)), Zonotope)

    # Fx is now indexed from 1 ... length(idx)
    return cluster(Fz, 1:length(idx), method)
end

# =====================================
# Box clustering with a vector output
# =====================================

struct BoxVecClustering{P} <: AbstractClusteringMethod{P}
    part::P
end

BoxVecClustering() = BoxVecClustering{Missing}(missing)

function cluster(F, idx, cl::BoxVecClustering)
    out = cluster(F, idx, BoxClustering(cl.part))
    return [out]
end

# =====================================
# Zonotope clustering
# =====================================

"""
    ZonotopeClustering{P} <: AbstractClusteringMethod{P}

### Notes

This method first takes a lazy convex hull for the given partition, then computes
a zonotope overapproximation of the convex hull.
"""
struct ZonotopeClustering{P} <: AbstractClusteringMethod{P}
    partition::P
end

ZonotopeClustering() = ZonotopeClustering(missing)

# for the generalization to > 2 ses, we iteratively apply the overapprox of the
# CH of two zonotopes, cf LazySets #2154
function cluster(F::Flowpipe{N, RT, VRT}, idx, ::ZonotopeClustering{Missing}) where {N, ZT<:Zonotope, RT<:ReachSet{N, ZT}, VRT<:AbstractVector{RT}}
    if length(idx) == 1
        return F[idx]
    elseif length(idx) == 2
        X = ConvexHull(set(F[idx[1]]), set(F[idx[2]]))
        Y = overapproximate(X, Zonotope) # TODO pass algorithm
        return [ReachSet(Y, tspan(F[idx]))]
    else
        Zaux = overapproximate(ConvexHull(set(F[idx[1]]), set(F[idx[2]])), Zonotope)
        for k in 3:length(idx)
            Zaux = overapproximate(ConvexHull(Zaux, set(F[idx[k]])), Zonotope)
        end
        return [ReachSet(Zaux, tspan(F[idx]))]
    end
end

# convexify and convert to vrep
#C = ReachabilityAnalysis.convexify(sol[end-aux+1:end])
#Cvertex = convex_hull(vcat([vertices_list(Z) for Z in LazySets.array(set(C))]...)) |> VPolygon

function _cartesian_product_standard(Z, blocks)
    @assert length.(blocks) == [2, 2, 1]

    # shuffle columns to bring to standard form
    M = Z.generators[:, [1, 2, 5, 6, 3, 4, 7, 8, 9]]

    # add extra column of zeros
    Mst = hcat(M[:, 1:4], zeros(5), M[:, 5:9])
    Zst = Zonotope(Z.center, Mst)
end

# TODO wrap into clustering strategy and pass outer partition
function _cluster_decompose_zono(fp::AbstractVector{<:TaylorModelReachSet}, blocks, dirs)
    # zonotopic enclosure
    fpz = [overapproximate(R, Zonotope) for R in fp]
    _cluster_decompose_zono(fpz, blocks, dirs)
end

function _cluster_decompose_zono(fpz::AbstractVector{RT}, blocks, dirs) where {N, ZT<:Zonotope, RT<:ReachSet{N, ZT}}

    # concrete projection onto each block
    πfpz = [[project(set(R), bi) for R in fpz] for bi in blocks]

    # reconstruct zonotope enclosure
    ZTD = Zonotope{Float64, Vector{Float64}, Matrix{Float64}}
    zono_blk = Vector{ZTD}(undef, length(blocks))

    for (i, V) in enumerate(πfpz)

        Vch = ConvexHullArray(V)
        dim_blk = dim(Vch)
        @assert dim_blk == 1 || dim_blk == 2

        if dim_blk == 2
            # epsilon-close approximation
            Yi = overapproximate(Vch, 1e-3)

            # zonotope overapproximation
            Yi = overapproximate(Yi, Zonotope, dirs)

            # order reduction
            Yi = reduce_order(Yi, 2)

        elseif dim_blk == 1
            Yi = overapproximate(Vch, LazySets.Interval)
            Yi = convert(Zonotope, Yi)

        end

        Yi = Zonotope(Vector(Yi.center), Matrix(Yi.generators))
        zono_blk[i] = Yi
    end

    # recomposition into a higher-dimensional zonotope
    out = concretize(CartesianProductArray(zono_blk))

    # bring to standard form
    out_st = _cartesian_product_standard(out, blocks)

    return out_st
end

# =====================================
# Lazy union set array clustering
# =====================================

"""
    UnionClustering{P} <: AbstractClusteringMethod{P}

Cluster according to the given partition by applying a lazy representation of the
set union.
"""
struct UnionClustering{P} <: AbstractClusteringMethod{P}
    partition::P
end

partition(method::UnionClustering) = method.partition

UnionClustering() = UnionClustering(missing)
UnionClustering(nchunks::D) where {D<:Integer} = UnionClustering{D}(nchunks)
UnionClustering(partition::VT) where {D<:Integer, VTi<:AbstractVector{D}, VT<:AbstractVector{VTi}} = UnionClustering{VT}(partition)

function cluster(F::Flowpipe{N, <:AbstractLazyReachSet}, idx, method::UnionClustering{Missing}) where {N}
    Fidx = view(F, idx)
    Δt = tspan(Fidx)
    Uidx = UnionSetArray([set(R) for R in Fidx])
    return [ReachSet(Uidx, Δt)]
end

function cluster(F::Flowpipe, idx, method::UnionClustering{P}) where {P}
    p = _partition(method, idx)
    return [cluster(F, cj, UnionClustering()) for cj in p]
end

# for Taylor model flowpipes we preprocess it with a zonotopic overapproximation
function cluster(F::Flowpipe{N, <:TaylorModelReachSet{N}}, idx, method::UnionClustering{P}) where {N, P}
    Fz = overapproximate(Flowpipe(view(F, idx)), Zonotope)

    # Fx is now indexed from 1 ... length(idx)
    return cluster(Fz, 1:length(idx), method)
end


text:
const message_severities = Dict(
    "debug" => VK_DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT,
    "info" => VK_DEBUG_UTILS_MESSAGE_SEVERITY_INFO_BIT_EXT,
    "warn" => VK_DEBUG_UTILS_MESSAGE_SEVERITY_WARNING_BIT_EXT,
    "error" => VK_DEBUG_UTILS_MESSAGE_SEVERITY_ERROR_BIT_EXT,
)

const message_types = Dict(
    "general" => VK_DEBUG_UTILS_MESSAGE_TYPE_GENERAL_BIT_EXT,
    "validation" => VK_DEBUG_UTILS_MESSAGE_TYPE_VALIDATION_BIT_EXT,
    "performance" => VK_DEBUG_UTILS_MESSAGE_TYPE_PERFORMANCE_BIT_EXT,
)

const message_types_r = reverse(message_types)
const message_severities_r = reverse(message_severities)

"""
    @vk_log VK_DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT str

Converts a Vulkan log to a Julia log
"""
macro vk_log(log_type_symbol, str)
    quote
        log_type = getfield(LibVulkan, log_type_symbol)
        if log_type == VK_DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT
            @debug($(esc(str)))
        elseif log_type == VK_DEBUG_UTILS_MESSAGE_SEVERITY_INFO_BIT_EXT
            @info($(esc(str)))
        elseif log_type == VK_DEBUG_UTILS_MESSAGE_SEVERITY_WARNING_BIT_EXT
            @warn($(esc(str)))
        elseif log_type == VK_DEBUG_UTILS_MESSAGE_SEVERITY_ERROR_BIT_EXT
            @error($(esc(str)))
        else
            throw(ArgumentError("Unknown logging type $log_type"))
        end
    end
end

function debug_callback(message_severity::VkDebugUtilsMessageSeverityFlagBitsEXT,
    message_type::VkDebugUtilsMessageTypeFlagsEXT,
    message,
    )
    @warn "$message_severity: $message_type: $message"
    @info "YYYOU LOGG"
    # @vk_log message_severity "$message_type: $message"
    return 0
end

function generate_debug_callback_interface(callback_f::Function)
    function debug_callback_interface(message_severity::VkDebugUtilsMessageSeverityFlagBitsEXT,
        message_type::VkDebugUtilsMessageTypeFlagsEXT,
        pCallbackData::Ptr{VkDebugUtilsMessengerCallbackDataEXT},
        pUserData::Ptr{Cvoid})
        callback_f(message_severity, message_type, unsafe_load(pCallbackData).pMessage)
    end
end

callback_func = generate_debug_callback_interface(debug_callback)
callback_func_c = @cfunction(callback_func, Cint, (UInt32, UInt32, Ptr{Cvoid}, Ptr{Cvoid}))

function debug_utils_messenger(inst; severity = "debug", types = ["general", "validation", "performance"])
    severity_hierarchy = ["debug", "info", "warn", "error"]
    index = findfirst(severity_hierarchy .== severity)
    severity_bits = [message_severities[key] for key in severity_hierarchy[index:end]]
    type_bits = [message_types[key] for key in types]
    messenger_info = Ref(VkDebugUtilsMessengerCreateInfoEXT(VK_STRUCTURE_TYPE_DEBUG_UTILS_MESSENGER_CREATE_INFO_EXT, C_NULL, 0, |(severity_bits...), |(type_bits...), callback_func_c, C_NULL))
    fname = "vkCreateDebugUtilsMessengerEXT"
    cfunc = Container(vkGetInstanceProcAddr(inst, pointer(fname)), Any[Ref(fname)])
    messenger = Ref{VkDebugUtilsMessengerEXT}()
    func = cfunc[]
    @vkcheck ccall(func, VkResult, (VkInstance, Ptr{VkDebugUtilsMessengerCreateInfoEXT}, Ptr{Nothing}, Ptr{VkDebugUtilsMessengerEXT}), inst, unsafe_pointer(messenger_info), C_NULL, messenger)
    # @vkcheck unsafe_load(func)(inst, unsafe_pointer(messenger_info), C_NULL, messenger)
    Container(messenger[], Any[messenger_info, cfunc])
end

# using Libdl
# find_library("libVkLayer_khronos_validation.so", ["/home/belmant/.julia/dev/Builders/Vulkan-ValidationLayers/build/layers/", "/home/belmant/.julia/dev/Builders/Vulkan-ValidationLayers/build/"])


text:
struct UniformFunctionSpace{vdim,sdim}
    vbasis::TensorProductBasis{vdim}
    sbasis::TensorProductBasis{sdim}
    vtpq::QuadratureRule{vdim}
    ftpq::Vector{QuadratureRule{sdim}}
    vquads
    fquads
    fnormals
    icoeffs
    iquad
    imap
    inormals
    isactiveface
    dofsperelement
    function UniformFunctionSpace(vbasis::TensorProductBasis{vdim},
        sbasis::TensorProductBasis{sdim},vtpq,ftpq,vquads,fquads,fnormals,
        icoeffs,iquad,imap,inormals) where {vdim,sdim}

            @assert sdim == vdim-1
            nphase,ncells = size(vquads)

            nfaces = number_of_faces(vdim)

            @assert length(ftpq) == nfaces
            @assert length(icoeffs) == ncells
            @assert size(fquads) == (nphase,ncells)
            @assert length(fnormals) == nfaces
            @assert all([length(n) == vdim for n in fnormals])
            nfaces = number_of_faces(vdim)
            @assert all(length.(fquads) .== nfaces)

            NHF = number_of_basis_functions(sbasis)
            dofsperelement = NHF*vdim

            isactiveface = [active_faces(fquads[p,c]) for p = 1:nphase,c = 1:ncells]

            return new{vdim,sdim}(vbasis,sbasis,vtpq,ftpq,vquads,
                fquads,fnormals,icoeffs,iquad,imap,inormals,isactiveface,
                dofsperelement)

        end
end

function UniformFunctionSpace(dgmesh,polyorder,quadorder,coeffs,levelset)

    vdim = dimension(dgmesh)
    sdim = vdim-1
    nfaces = number_of_faces(vdim)

    facemaps = dgmesh.facemaps
    fnormals = reference_normals()

    vbasis = TensorProductBasis(vdim,polyorder)
    sbasis = TensorProductBasis(sdim,polyorder)

    quad1d = ImplicitDomainQuadrature.ReferenceQuadratureRule(quadorder)
    vtpq = tensor_product(quad1d,reference_cell(vdim))
    iquad = tensor_product(quad1d,reference_cell(sdim))

    ftpq = repeat([iquad],nfaces)

    imap = InterpolatingPolynomial(vdim,sbasis)

    vquads = element_quadratures(dgmesh.cellsign,coeffs,levelset,quad1d)
    fquads = face_quadratures(dgmesh.cellsign,coeffs,levelset,facemaps,quad1d)

    icoeffs = interface_coefficients(dgmesh.cellsign,coeffs,levelset,
        sbasis,iquad)
    inormals = interface_normals(dgmesh,icoeffs,imap,coeffs,levelset,
        iquad.points)

    return UniformFunctionSpace(vbasis,sbasis,vtpq,ftpq,vquads,fquads,fnormals,
        icoeffs,iquad,imap,inormals)
end

function interface_normals(dgmesh,icoeffs,imap,lcoeffs,levelset,
    refpoints)

    sdim,nqp = size(refpoints)
    dim = sdim+1

    cellsign = dgmesh.cellsign

    ncells = length(cellsign)
    ft = default_float_type()
    inormals = Vector{Matrix{ft}}(undef,ncells)
    cellmap = CellMap(dim)

    for cellid in 1:ncells
        if cellsign[cellid] == 0
            update!(levelset,lcoeffs[:,cellid])
            update!(imap,icoeffs[cellid])
            update_range!(cellmap,dgmesh.domain[cellid])

            mappedpoints = hcat([imap(refpoints[:,i]) for i in 1:size(refpoints)[2]]...)
            inormals[cellid] = levelset_normal(levelset,mappedpoints,cellmap)
        end
    end
    return inormals
end

function element_quadratures!(equads,cellsign,coeffs,poly,tpq,box,quad1d)

    nf,ncells = size(coeffs)
    @assert length(cellsign) == ncells
    @assert size(equads) == (2,ncells)

    for idx in 1:ncells
        s = cellsign[idx]
        if s == +1
            equads[1,idx] = tpq
        elseif s == -1
            equads[2,idx] = tpq
        elseif s == 0
            update!(poly,coeffs[:,idx])
            equads[1,idx] = quadrature(poly,+1,false,box,quad1d)
            equads[2,idx] = quadrature(poly,-1,false,box,quad1d)
        else
            throw(ArgumentError("Expected `cellsign[idx] ∈ {-1,0,+1}`,
                got `cellsign[$idx] == $s`"))
        end
    end
end

function element_quadratures!(equads,cellsign,coeffs,poly,quad1d)

    nf,ncells = size(coeffs)
    @assert length(cellsign) == ncells
    @assert size(equads) == (2,ncells)

    dim = dimension(poly)
    box = reference_cell(dim)

    tpq = tensor_product(quad1d,box)

    element_quadratures!(equads,cellsign,coeffs,poly,tpq,box,quad1d)
end

function element_quadratures(cellsign,coeffs,poly,quad1d)
    dim = dimension(poly)
    ncells = length(cellsign)
    equads = Array{QuadratureRule{dim}}(undef,2,ncells)
    element_quadratures!(equads,cellsign,coeffs,poly,quad1d)
    return equads
end

function update_face_quadratures!(facequads,levelset,facemaps,
    sign_condition,quad1d)

    nfaces = length(facequads)
    @assert length(facemaps) == nfaces
    for faceid in 1:nfaces
        fmap = facemaps[faceid]
        xiL,xiR = reference_interval(fmap)
        quad = QuadratureRule(quadrature([x->levelset(fmap(x))],
            [sign_condition],xiL,xiR,quad1d))
        facequads[faceid] = quad
    end
end

function assign_all(array,value,ncells)
    for i = 1:ncells
        array[i] = value
    end
end

function face_quadratures!(facequads,cellsign,coeffs,levelset,facemaps,
    quad1d)

    nf,ncells = size(coeffs)
    @assert size(facequads) == (2,ncells)

    nfaces = length(facemaps)

    @assert all([length(facequads[1,i]) == nfaces for i = 1:ncells])
    @assert all([length(facequads[2,i]) == nfaces for i = 1:ncells])

    dim = dimension(levelset)
    reference_face = reference_cell(dim-1)

    tpq = tensor_product(quad1d,reference_face)
    zeroquad = blank_quadrature_rule(dim-1)

    for cellid in 1:ncells
        if cellsign[cellid] == +1
            assign_all(facequads[1,cellid],tpq,nfaces)
            assign_all(facequads[2,cellid],zeroquad,nfaces)
        elseif cellsign[cellid] == -1
            assign_all(facequads[2,cellid],tpq,nfaces)
            assign_all(facequads[1,cellid],zeroquad,nfaces)
        elseif cellsign[cellid] == 0

            update!(levelset,coeffs[:,cellid])

            update_face_quadratures!(facequads[1,cellid],levelset,
                facemaps,+1,quad1d)
            update_face_quadratures!(facequads[2,cellid],levelset,
                facemaps,-1,quad1d)

        end
    end

end

function blank_quadrature_rule(dim)
    ft = default_float_type()
    return QuadratureRule(ImplicitDomainQuadrature.TemporaryQuadrature(ft,dim))
end

function face_quadratures(cellsign,coeffs,levelset,facemaps,quad1d)

    nfuncs,ncells = size(coeffs)
    @assert length(cellsign) == ncells

    nfaces = length(facemaps)

    dim = dimension(levelset)
    facedim = dim-1
    facequads = reshape([Vector{QuadratureRule{facedim}}(undef,nfaces) for i = 1:2ncells],2,ncells)
    face_quadratures!(facequads,cellsign,coeffs,levelset,facemaps,quad1d)
    return facequads

end

function interface_coefficients!(icoeffs,cellsign,coeffs,levelset,basis,quad)

    dim = dimension(levelset)

    sdim = dimension(basis)
    @assert sdim == dim-1
    @assert dimension(quad) == sdim

    ncells = length(cellsign)
    nf = number_of_basis_functions(basis)

    cell = reference_cell(dim)
    mass = lu(mass_matrix(basis,quad,1.0,dim))

    for cellid in 1:ncells
        if cellsign[cellid] == 0
            update!(levelset,coeffs[:,cellid])
            icoeffs[cellid] = fit_zero_levelset(levelset,basis,quad,mass,cell)
        end
    end

end

function interface_coefficients(cellsign,coeffs,levelset,basis,quad)
    nf = number_of_basis_functions(basis)
    dim = dimension(levelset)
    ncells = length(cellsign)
    ft = default_float_type()
    icoeffs = Vector{Vector{ft}}(undef,ncells)
    interface_coefficients!(icoeffs,cellsign,coeffs,levelset,basis,quad)
    return icoeffs
end


text:
abstract type SlidingWindow{TElem,TData,O} <: DataView{TElem,TData} end

Base.parent(A::SlidingWindow) = A.data
Base.length(A::SlidingWindow) = nobs(A)
nobs(A::SlidingWindow) = A.count

allowcontainer(fun, ::SlidingWindow) = true

# compatibility with nested functions
default_obsdim(A::SlidingWindow) = A.obsdim

function _windowsettings(A::SlidingWindow, windowindex::Int)
    offset = 1 + (windowindex-1) * A.stride
    range = offset:offset+A.size-1
    offset, range
end

function _check_windowargs(data, size::Int, stride::Int, obsdim)
    size > 0 || throw(ArgumentError("Specified window size must be strictly greater than 0. Actual: $size"))
    size <= nobs(data,obsdim) || throw(ArgumentError("Specified window size is too large for the given number of observations"))
    stride > 0 || throw(ArgumentError("Specified stride must be strictly greater than 0. Actual: $stride"))
end

# --------------------------------------------------------------------

struct UnlabeledSlidingWindow{TElem,TData,O} <: SlidingWindow{TElem,TData,O}
    data::TData
    size::Int
    stride::Int
    count::Int
    obsdim::O
end

"""
    slidingwindow(data, size, [stride], [obsdim]) -> UnlabeledSlidingWindow

Return a vector-like view of the `data` for which each element is
a fixed size "window" of `size` adjacent observations. By default
these windows are not overlapping. Note that only complete
windows are included in the output, which implies that it is
possible for excess observations to be omitted from the view.

```julia-repl
julia> A = slidingwindow(1:20, 6)
3-element slidingwindow(::UnitRange{Int64}, 6) with element type SubArray{Int64,1,UnitRange{Int64},Tuple{UnitRange{Int64}},true}:
 [1, 2, 3, 4, 5, 6]
 [7, 8, 9, 10, 11, 12]
 [13, 14, 15, 16, 17, 18]
```

Note that the values of `data` itself are not copied. Instead the
function [`datasubset`](@ref) is called when `getindex` is
invoked. To actually get a copy of the data at some window use
the function [`getobs`](@ref).

```julia-repl
julia> A[1]
6-element SubArray{Int64,1,UnitRange{Int64},Tuple{UnitRange{Int64}},true}:
 1
 ⋮
 6

julia> getobs(A, 1)
6-element Array{Int64,1}:
 1
 ⋮
 6
```

The optional parameter `stride` can be used to specify the
distance between the start elements of each adjacent window.
By default the stride is equal to the window size.

```julia-repl
julia> slidingwindow(1:20, 6, stride = 3)
5-element slidingwindow(::UnitRange{Int64}, 6, stride = 3) with element type SubArray{Int64,1,UnitRange{Int64},Tuple{UnitRange{Int64}},true}:
 [1, 2, 3, 4, 5, 6]
 [4, 5, 6, 7, 8, 9]
 [7, 8, 9, 10, 11, 12]
 [10, 11, 12, 13, 14, 15]
 [13, 14, 15, 16, 17, 18]
```

The optional (keyword) parameter `obsdim` allows one to specify
which dimension denotes the observations. see `LearnBase.ObsDim`
for more detail.
"""
function slidingwindow(data::T, size::Int, stride::Int, obsdim::O = default_obsdim(data)) where {T,O}
    _check_windowargs(data, size, stride, obsdim)
    E = typeof(datasubset(data,1:size,obsdim))
    count = floor(Int, (nobs(data,obsdim) - size + stride) / stride)
    UnlabeledSlidingWindow{E,T,O}(data, size, stride, count, obsdim)
end

function slidingwindow(data, size::Int, obsdim::Union{Tuple,ObsDimension})
    # default to using same stride as size
    slidingwindow(data, size, size, obsdim)
end

function slidingwindow(data, size::Int; stride=size, obsdim=default_obsdim(data))
    slidingwindow(data, size, stride, convert(ObsDimension, obsdim))
end

function Base.getindex(A::UnlabeledSlidingWindow, windowindex::Int)
    _, windowrange = _windowsettings(A, windowindex)
    datasubset(A.data, windowrange, A.obsdim)
end

getobs(A::UnlabeledSlidingWindow, indices::AbstractVector) = [getobs(A,i) for i in indices]

# --------------------------------------------------------------------

struct LabeledSlidingWindow{TElem,TData,O,TFun,Exclude} <: SlidingWindow{TElem,TData,O}
    data::TData
    targetfun::TFun
    size::Int
    stride::Int
    count::Int
    offset::Int
    obsdim::O
end

"""
    slidingwindow(f, data, size, [stride], [excludetarget], [obsdim]) -> LabeledSlidingWindow

Return a vector-like view of the `data` for which each element is
a tuple of two elements:

1. A fixed size "window" of `size` adjacent observations. By
   default these windows are not overlapping. This can be
   changed by explicitly specifying a `stride`.

2. A single target (or vector of targets) for the window. The
   content of the target(s) is defined by the label-index
   function `f`.

The label-index function `f` is a unary function that takes the
index of the first observation in the current window and should
return the index (or indices) of the associated target(s) for
that window.

```julia-repl
julia> A = slidingwindow(i->i+6, 1:20, 6)
3-element slidingwindow(::##3#4, ::UnitRange{Int64}, 6) with element type Tuple{...}
 ([1, 2, 3, 4, 5, 6], 7)
 ([7, 8, 9, 10, 11, 12], 13)
 ([13, 14, 15, 16, 17, 18], 19)
```

Note that only complete and in-bound windows are included in the
output, which implies that it is possible for excess observations
to be omitted from the resulting view.

```julia-repl
julia> A = slidingwindow(i->i-1, 1:20, 6)
2-element slidingwindow(::##5#6, ::UnitRange{Int64}, 6) with element type Tuple{...}
 ([7, 8, 9, 10, 11, 12], 6)
 ([13, 14, 15, 16, 17, 18], 12)
```

As hinted above, it is also allowed for `f` to return a vector of
indices. This can be useful for emulating techniques such as
skip-gram.

```julia-repl
julia> data = split("The quick brown fox jumps over the lazy dog")
9-element Array{SubString{String},1}:
 "The"
 "quick"
 ⋮
 "lazy"
 "dog"

julia> A = slidingwindow(i->[i-2:i-1; i+1:i+2], data, 1)
5-element slidingwindow(::##11#12, ::Array{SubString{String},1}, 1) with element type Tuple{...}:
 (["brown"], ["The", "quick", "fox", "jumps"])
 (["fox"], ["quick", "brown", "jumps", "over"])
 (["jumps"], ["brown", "fox", "over", "the"])
 (["over"], ["fox", "jumps", "the", "lazy"])
 (["the"], ["jumps", "over", "lazy", "dog"])
```

Should it so happen that the targets overlap with the features,
then the affected observation(s) will be present in both. To
change this behaviour one can set the optional parameter
`excludetarget = true`. This will remove the target(s) from the
feature window.

```julia-repl
julia> slidingwindow(i->i+2, data, 5, stride = 1, excludetarget = true)
5-element slidingwindow(::##17#18, ::Array{SubString{String},1}, 5, stride = 1) with element type Tuple{...}:
 (["The", "quick", "fox", "jumps"], "brown")
 (["quick", "brown", "jumps", "over"], "fox")
 (["brown", "fox", "over", "the"], "jumps")
 (["fox", "jumps", "the", "lazy"], "over")
 (["jumps", "over", "lazy", "dog"], "the")
```

The optional (keyword) parameter `obsdim` allows one to specify
which dimension denotes the observations. see `LearnBase.ObsDim`
for more detail.
"""
function slidingwindow(f::F, data::T, size::Int, stride::Int, ::Val{Exclude}=Val(false), obsdim::O = default_obsdim(data)) where {F,T,Exclude,O}
    _check_windowargs(data, size, stride, obsdim)
    count = floor(Int, (nobs(data,obsdim) - size + stride) / stride)
    # limit back of sequence until it contains target
    o = 1 + (count-1) * stride
    while maximum(f(o)) > nobs(data,obsdim)
        count -= 1
        o = 1 + (count-1) * stride
    end
    # limit front of sequence until it contains target
    offset = 0
    o = 1 + offset * stride
    while minimum(f(o)) < 1
        offset += 1
        count -= 1
        o = 1 + offset * stride
    end
    # construct view
    E = typeof((datasubset(data,o:o+size-1,obsdim),datasubset(data,f(o),obsdim)))
    LabeledSlidingWindow{E,T,O,F,Exclude}(data, f, size, stride, count, offset, obsdim)
end

function slidingwindow(f, data, size::Int, stride::Int, obsdim::Union{ObsDimension,Tuple}, ::Val{Exclude}=Val(false)) where {Exclude}
    slidingwindow(f, data, size, stride, Val(Exclude), obsdim)
end

function slidingwindow(f::Function, data, size::Int, obsdim::Union{ObsDimension,Tuple}, ::Val{Exclude}=Val(false)) where {Exclude}
    slidingwindow(f, data, size, size, Val(Exclude), obsdim)
end

function slidingwindow(f::Function, data, size::Int, ::Val{Exclude}, obsdim = default_obsdim(data)) where {Exclude}
    slidingwindow(f, data, size, size, Val(Exclude), obsdim)
end

Base.@pure _toVal(::Val{T}) where {T} = Val(T)
Base.@pure _toVal(T) = Val(T)
function slidingwindow(f, data, size::Int; stride=size, excludetarget=Val(false), obsdim=default_obsdim(data))
    slidingwindow(f, data, size, stride, _toVal(excludetarget), convert(ObsDimension, obsdim))
end

function Base.getindex(A::LabeledSlidingWindow{TElem,TData,O,TFun,Exclude}, wi::Int) where {TElem,TData,O,TFun,Exclude}
    windowindex = wi + A.offset
    offset, windowrange = _windowsettings(A, windowindex)
    windowindices = if Exclude
        setdiff(windowrange, A.targetfun(offset))
    else
        windowrange
    end
    X = datasubset(A.data, windowindices, A.obsdim)
    Y = datasubset(A.data, A.targetfun(offset), A.obsdim)
    X, Y
end

function getobs(A::LabeledSlidingWindow, indices::AbstractVector)
    map(i->getobs(A,i), indices)
end

function getobs(A::LabeledSlidingWindow{<:NTuple{<:Any,Tuple}}, i::Int)
    map(a->map(getobs,a), A[i])
end

# --------------------------------------------------------------------

function Base.showarg(io::IO, A::SlidingWindow, toplevel)
    print(io, "slidingwindow(")
    if A isa LabeledSlidingWindow
        Base.showarg(io, A.targetfun, false)
        print(io, ", ")
    end
    Base.showarg(io, parent(A), false)
    print(io, ", ")
    print(io, A.size)
    if A.stride != A.size
        print(io, ", ")
        print(io, "stride = ", A.stride)
    end
    if A.obsdim != default_obsdim(A.data)
        print(io, ", ")
        print(io, "obsdim = ", obsdim_string(A.obsdim))
    end
    print(io, ')')
    toplevel && print(io, " with eltype ", nameof(eltype(A)))
end

# --------------------------------------------------------------------
# batches of sliding windows

const WindowBatchView{TElem,O} = BatchView{TElem,<:SlidingWindow,O}

function BatchView(data::T, size::Int, count::Int, obsdim::O = default_obsdim(data), upto::Bool = false) where {T<:SlidingWindow,O}
    nsize, ncount = _compute_batch_settings(data, size, count, obsdim, upto)
    E = typeof(_getwindowbatch(data, nsize, 1))
    BatchView{E,T,O}(data, nsize, ncount, obsdim)
end

function Base.getindex(A::WindowBatchView, batchindex::Int)
    _getwindowbatch(parent(A), A.size, batchindex)
end

function Base.getindex(A::WindowBatchView, batchindices::AbstractVector)
    [A[bi] for bi in batchindices]
end

# shuffling/resampling like this is broken
allowcontainer(fun, ::BatchView{<:Any,<:SlidingWindow}) = false
allowcontainer(::typeof(shuffleobs), ::BatchView{<:Any,<:SlidingWindow}) = false
allowcontainer(::typeof(splitobs), ::BatchView{<:Any,<:SlidingWindow}) = false

# --------------------------------------------------------------------

getobs(A::BatchView{T,<:UnlabeledSlidingWindow}) where {T} = map(x->getobs.(x), A)
getobs(A::BatchView{T,<:UnlabeledSlidingWindow}, i::AbstractVector) where {T} = map(x->getobs.(x), A)
getobs(A::BatchView{T,<:UnlabeledSlidingWindow}, i::Int) where {T} = getobs.(A[i])

function _getwindowbatch(W::UnlabeledSlidingWindow, batchsize::Int, batchindex::Int)
    batchrange = _batchrange(batchsize, batchindex)
    windowranges = [_windowsettings(W, windowindex)[2] for windowindex in batchrange]
    [datasubset(W.data, [windowranges[wi][oi] for wi in 1:batchsize], W.obsdim) for oi in 1:W.size]
end

# --------------------------------------------------------------------

function _getwindowbatch(A::LabeledSlidingWindow{TElem,TData,O,TFun,Exclude}, batchsize::Int, batchindex::Int) where {TElem,TData,O,TFun,Exclude}
    batchrange = _batchrange(batchsize, batchindex)
    fltr = if Exclude
        s -> (A.targetfun(s[1]), setdiff(s[2], A.targetfun(s[1])))
    else
        s -> (A.targetfun(s[1]), s[2])
    end
    windowsettings = [fltr(_windowsettings(A, windowindex)) for windowindex in batchrange]
    ([datasubset(A.data, [windowsettings[wi][2][oi] for wi in 1:batchsize], A.obsdim)
        for oi in 1:length(windowsettings[1][2])],
     [datasubset(A.data, [windowsettings[wi][1][oi] for wi in 1:batchsize], A.obsdim)
        for oi in 1:length(windowsettings[1][1])])
end


text:
import Plots, SearchBall
const sb = SearchBall
const pl = Plots

pl.gr()

function end_pos(elapsed::Number, v = 1.0, ω = π/16.0, maximumtime::Number=2π/ω)
    init = [1, 0]
    rot = elapsed*ω
    end_pos = (maximumtime - elapsed).*v .* sb.rotate2d(init, rot)
    end_pos
end

path = end_pos.(0:0.01:32)

x = getindex.(path, 1)
y = getindex.(path, 2)

pl.plot(x, y, aspect_ratio=:equal)

function spiral(θ::Number, a::Number=1, b::Number=1)
    a + b*θ
end

θ = 0:0.01:2π
pl.plot(θ, spiral.(θ, 0, 2.0), proj = :polar)

function spiral_param(φ::Number, a::Number=1)
    [a*φ*cos(φ), a*φ*sin(φ)]
end

φ = 0:0.01:2π
vals = spiral_param.(φ)
x_spiral = getindex.(vals, 1)
y_spiral = getindex.(vals, 2)
pl.plot(x_spiral, y_spiral, aspect_ratio=:equal)


text:
module DiffEqFlux

using DiffEqBase, Tracker, DiffResults, DiffEqSensitivity, ForwardDiff,
      Flux, Requires, RecursiveArrayTools, Adapt, LinearAlgebra

gpu_or_cpu(x) = Array
function __init__()
    @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
        gpu_or_cpu(x::CuArrays.CuArray) = CuArrays.CuArray
        gpu_or_cpu(x::TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}) = CuArrays.CuArray
        gpu_or_cpu(x::Transpose{<:Any,<:CuArrays.CuArray}) = CuArrays.CuArray
        gpu_or_cpu(x::Adjoint{<:Any,<:CuArrays.CuArray}) = CuArrays.CuArray
        gpu_or_cpu(x::Adjoint{<:Any,TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}}) = CuArrays.CuArray
        gpu_or_cpu(x::Transpose{<:Any,TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}}) = CuArrays.CuArray
    end
end

include("Flux/layers.jl")
include("Flux/neural_de.jl")
include("Flux/utils.jl")

export diffeq_fd, diffeq_rd, diffeq_adjoint
export neural_ode, neural_ode_rd
export neural_dmsde
end


text:
"""
    test_manifold(
        M::Manifold,
        pts::AbstractVector;
        args,
    )

Test general properties of manifold `M`, given at least three different points
that lie on it (contained in `pts`).

# Arguments
- `basis_has_specialized_diagonalizing_get = false`: if true, assumes that
    [`DiagonalizingOrthonormalBasis`](@ref) given in `basis_types` has
    [`get_coordinates`](@ref) and [`get_vector`](@ref) that work without caching.
- `basis_types_to_from = ()`: basis types that will be tested based on
    [`get_coordinates`](@ref) and [`get_vector`](@ref).
- `basis_types_vecs = ()` : basis types that will be tested based on [`get_vectors`](@ref).
- `default_inverse_retraction_method = ManifoldsBase.LogarithmicInverseRetraction()`:
    default method for inverse retractions ([`log`](@ref)).
- `default_retraction_method = ManifoldsBase.ExponentialRetraction()`: default method for
    retractions ([`exp`](@ref)).
- `exp_log_atol_multiplier = 0`: change absolute tolerance of exp/log tests
    (0 use default, i.e. deactivate atol and use rtol).
- `exp_log_rtol_multiplier = 1`: change the relative tolerance of exp/log tests
    (1 use default). This is deactivated if the `exp_log_atol_multiplier` is nonzero.
- `expected_dimension_type = Integer`: expected type of value returned by
    [`manifold_dimension`](@ref).
- `inverse_retraction_methods = []`: inverse retraction methods that will be tested.
- `is_mutating = true`: whether mutating variants of functions should be tested.
- `is_point_atol_multiplier = 0`: determines atol of `is_manifold_point` checks.
- `is_tangent_atol_multiplier = 0`: determines atol of `is_tangent_vector` checks.
- `mid_point12 = test_exp_log ? shortest_geodesic(M, pts[1], pts[2], 0.5) : nothing`: if not `nothing`, then check
    that `mid_point(M, pts[1], pts[2])` is approximately equal to `mid_point12`. This is
    by default set to `nothing` if `text_exp_log` is set to false.
- `point_distributions = []` : point distributions to test.
- `rand_tvector_atol_multiplier = 0` : chage absolute tolerance in testing random vectors
    (0 use default, i.e. deactivate atol and use rtol) random tangent vectors are tangent
    vectors.
- `retraction_atol_multiplier = 0`: change absolute tolerance of (inverse) retraction tests
    (0 use default, i.e. deactivate atol and use rtol).
- `retraction_rtol_multiplier = 1`: change the relative tolerance of (inverse) retraction
    tests (1 use default). This is deactivated if the `exp_log_atol_multiplier` is nonzero.
- `retraction_methods = []`: retraction methods that will be tested.
- `test_exp_log = true`: if true, checkthat [`exp`](@ref) is the inverse of [`log`](@ref).
- `test_forward_diff = true`: if true, automatic differentiation using
    ForwardDiff is tested.
- `test_injectivity_radius = true`: whether implementation of [`injectivity_radius`](@ref)
    should be tested.
- `test_is_tangent`: if true check that the `default_inverse_retraction_method`
    actually returns valid tangent vectors.
- `test_musical_isomorphisms = false` : test musical isomorphisms.
- `test_mutating_rand = false` : test the mutating random function for points on manifolds.
- `test_project_point = false`: test projections onto the manifold.
- `test_project_tangent = false` : test projections on tangent spaces.
- `test_representation_size = true` : test repersentation size of points/tvectprs.
- `test_reverse_diff = true`: if true, automatic differentiation using
    ReverseDiff is tested.
- `test_tangent_vector_broadcasting = true` : test boradcasting operators on TangentSpace.
- `test_vector_spaces = true` : test Vector bundle of this manifold.
- `test_default_vector_transport = false` : test the default vector transport (usually
    parallel transport).
- `test_vee_hat = false`: test [`vee`](@ref) and [`hat`](@ref) functions.
- `tvector_distributions = []` : tangent vector distributions to test.
- `vector_transport_methods = []`: vector transport methods that should be tested.
"""
function ManifoldTests.test_manifold(
    M::Manifold,
    pts::AbstractVector;
    basis_has_specialized_diagonalizing_get = false,
    basis_types_to_from = (),
    basis_types_vecs = (),
    default_inverse_retraction_method = LogarithmicInverseRetraction(),
    default_retraction_method = ExponentialRetraction(),
    exp_log_atol_multiplier = 0,
    exp_log_rtol_multiplier = 1,
    expected_dimension_type = Integer,
    inverse_retraction_methods = [],
    is_mutating = true,
    is_point_atol_multiplier = 0,
    is_tangent_atol_multiplier = 0,
    point_distributions = [],
    projection_atol_multiplier = 0,
    rand_tvector_atol_multiplier = 0,
    retraction_atol_multiplier = 0,
    retraction_methods = [],
    retraction_rtol_multiplier = 1,
    test_exp_log = true,
    test_forward_diff = true,
    test_is_tangent = true,
    test_injectivity_radius = true,
    test_musical_isomorphisms = false,
    test_mutating_rand = false,
    test_project_point = false,
    test_project_tangent = false,
    test_representation_size = true,
    test_reverse_diff = true,
    test_tangent_vector_broadcasting = true,
    test_default_vector_transport = false,
    test_vector_spaces = true,
    test_vee_hat = false,
    tvector_distributions = [],
    vector_transport_methods = [],
    mid_point12 = test_exp_log ? shortest_geodesic(M, pts[1], pts[2], 0.5) : nothing,
)

    length(pts) ≥ 3 || error("Not enough points (at least three expected)")
    isapprox(M, pts[1], pts[2]) && error("Points 1 and 2 are equal")
    isapprox(M, pts[1], pts[3]) && error("Points 1 and 3 are equal")

    # get a default tangent vector for every of the three tangent spaces
    n = length(pts)
    if default_inverse_retraction_method === nothing
        tv = [zero_tangent_vector(M, pts[i]) for i in 1:n] # no other available
    else
        tv = [
            inverse_retract(
                M,
                pts[i],
                pts[((i + 1) % n) + 1],
                default_inverse_retraction_method,
            ) for i in 1:n
        ]
    end
    Test.Test.@testset "dimension" begin
        Test.@test isa(manifold_dimension(M), expected_dimension_type)
        Test.@test manifold_dimension(M) ≥ 0
        Test.@test manifold_dimension(M) ==
                   vector_space_dimension(Manifolds.VectorBundleFibers(
            Manifolds.TangentSpace,
            M,
        ))
        Test.@test manifold_dimension(M) ==
                   vector_space_dimension(Manifolds.VectorBundleFibers(
            Manifolds.CotangentSpace,
            M,
        ))
    end

    test_representation_size && Test.@testset "representation" begin
        function test_repr(repr)
            Test.@test isa(repr, Tuple)
            for rs in repr
                Test.@test rs > 0
            end
            return nothing
        end

        test_repr(Manifolds.representation_size(M))
        for fiber in (Manifolds.TangentSpace, Manifolds.CotangentSpace)
            test_repr(Manifolds.representation_size(Manifolds.VectorBundleFibers(fiber, M)))
        end
    end

    test_injectivity_radius && Test.@testset "injectivity radius" begin
        Test.@test injectivity_radius(M, pts[1]) > 0
        Test.@test injectivity_radius(M, pts[1]) ≥ injectivity_radius(M)
        for rm in retraction_methods
            Test.@test injectivity_radius(M, rm) > 0
            Test.@test injectivity_radius(M, pts[1], rm) ≥ injectivity_radius(M, rm)
            Test.@test injectivity_radius(M, pts[1], rm) ≤ injectivity_radius(M, pts[1])
        end
    end

    Test.@testset "is_manifold_point" begin
        for pt in pts
            atol = is_point_atol_multiplier * ManifoldTests.find_eps(pt)
            Test.@test is_manifold_point(M, pt; atol = atol)
            Test.@test check_manifold_point(M, pt; atol = atol) === nothing
        end
    end

    test_is_tangent && Test.@testset "is_tangent_vector" begin
        for (p, X) in zip(pts, tv)
            atol = is_tangent_atol_multiplier * ManifoldTests.find_eps(p)
            if !(check_tangent_vector(M, p, X; atol = atol) === nothing)
                print(check_tangent_vector(M, p, X; atol = atol))
            end
            Test.@test is_tangent_vector(M, p, X; atol = atol)
            Test.@test check_tangent_vector(M, p, X; atol = atol) === nothing
        end
    end

    test_exp_log && Test.@testset "log/exp tests" begin
        epsp1p2 = ManifoldTests.find_eps(pts[1], pts[2])
        atolp1p2 = exp_log_atol_multiplier * epsp1p2
        rtolp1p2 =
            exp_log_atol_multiplier == 0.0 ? sqrt(epsp1p2) * exp_log_rtol_multiplier : 0
        X1 = log(M, pts[1], pts[2])
        X2 = log(M, pts[2], pts[1])
        Test.@test isapprox(M, pts[2], exp(M, pts[1], X1); atol = atolp1p2, rtol = rtolp1p2)
        Test.@test isapprox(
            M,
            pts[1],
            exp(M, pts[1], X1, 0);
            atol = atolp1p2,
            rtol = rtolp1p2,
        )
        Test.@test isapprox(
            M,
            pts[2],
            exp(M, pts[1], X1, 1);
            atol = atolp1p2,
            rtol = rtolp1p2,
        )
        if VERSION >= v"1.5" && isa(M, Union{Grassmann,GeneralizedStiefel})
            # TODO: investigate why this is so imprecise on newer Julia versions on CI
            Test.@test isapprox(
                M,
                pts[1],
                exp(M, pts[2], X2);
                # yields 5*10^-8 for the usual 10^-13 we impose on earlier Julia versions
                atol = atolp1p2 * 5 * 10^5,
                rtol = rtolp1p2,
            )
        else
            Test.@test isapprox(M, pts[1], exp(M, pts[2], X2); atol = atolp1p2, rtol = rtolp1p2)
        end
        Test.@test is_manifold_point(
            M,
            exp(M, pts[1], X1);
            atol = atolp1p2,
            rtol = rtolp1p2,
        )
        Test.@test isapprox(
            M,
            pts[1],
            exp(M, pts[1], X1, 0);
            atol = atolp1p2,
            rtol = rtolp1p2,
        )
        for p in pts
            epsx = ManifoldTests.find_eps(p)
            Test.@test isapprox(
                M,
                p,
                zero_tangent_vector(M, p),
                log(M, p, p);
                atol = epsx * exp_log_atol_multiplier,
                rtol = exp_log_atol_multiplier == 0.0 ?
                       sqrt(epsx) * exp_log_rtol_multiplier : 0,
            )
            Test.@test isapprox(
                M,
                p,
                zero_tangent_vector(M, p),
                inverse_retract(M, p, p);
                atol = epsx * exp_log_atol_multiplier,
                rtol = exp_log_atol_multiplier == 0.0 ?
                       sqrt(epsx) * exp_log_rtol_multiplier : 0.0,
            )
        end
        atolp1 = exp_log_atol_multiplier * ManifoldTests.find_eps(pts[1])
        if is_mutating
            zero_tangent_vector!(M, X1, pts[1])
        else
            X1 = zero_tangent_vector(M, pts[1])
        end
        Test.@test isapprox(M, pts[1], X1, zero_tangent_vector(M, pts[1]); atol = atolp1)
        if is_mutating
            log!(M, X1, pts[1], pts[2])
        else
            X1 = log(M, pts[1], pts[2])
        end

        Test.@test isapprox(M, exp(M, pts[1], X1, 1), pts[2]; atol = atolp1)
        Test.@test isapprox(M, exp(M, pts[1], X1, 0), pts[1]; atol = atolp1)

        Test.@test distance(M, pts[1], pts[2]) ≈ norm(M, pts[1], X1)

        X3 = log(M, pts[1], pts[3])

        Test.@test inner(M, pts[1], X1, X3) ≈ conj(inner(M, pts[1], X3, X1))
        Test.@test inner(M, pts[1], X1, X1) ≈ real(inner(M, pts[1], X1, X1))

        Test.@test norm(M, pts[1], X1) isa Real
        Test.@test norm(M, pts[1], X1) ≈ sqrt(inner(M, pts[1], X1, X1))
    end

    Test.@testset "(inverse &) retraction tests" begin
        for (p, X) in zip(pts, tv)
            epsx = ManifoldTests.find_eps(p)
            for retr_method in retraction_methods
                Test.@test is_manifold_point(M, retract(M, p, X, retr_method))
                Test.@test isapprox(
                    M,
                    p,
                    retract(M, p, X, 0, retr_method);
                    atol = epsx * retraction_atol_multiplier,
                    rtol = retraction_atol_multiplier == 0 ?
                           sqrt(epsx) * retraction_rtol_multiplier : 0,
                )
                if is_mutating
                    new_pt = allocate(p)
                    retract!(M, new_pt, p, X, retr_method)
                else
                    new_pt = retract(M, p, X, retr_method)
                end
                Test.@test is_manifold_point(M, new_pt)
            end
        end
        for p in pts
            epsx = ManifoldTests.find_eps(p)
            for inv_retr_method in inverse_retraction_methods
                Test.@test isapprox(
                    M,
                    p,
                    zero_tangent_vector(M, p),
                    inverse_retract(M, p, p, inv_retr_method);
                    atol = epsx * retraction_atol_multiplier,
                    rtol = retraction_atol_multiplier == 0 ?
                           sqrt(epsx) * retraction_rtol_multiplier : 0,
                )
            end
        end
    end

    test_vector_spaces && Test.@testset "vector spaces tests" begin
        for p in pts
            X = zero_tangent_vector(M, p)
            mts = Manifolds.VectorBundleFibers(Manifolds.TangentSpace, M)
            Test.@test isapprox(M, p, X, zero_vector(mts, p))
            if is_mutating
                zero_vector!(mts, X, p)
                Test.@test isapprox(M, p, X, zero_tangent_vector(M, p))
            end
        end
    end

    Test.@testset "basic linear algebra in tangent space" begin
        for (p, X) in zip(pts, tv)
            Test.@test isapprox(
                M,
                p,
                0 * X,
                zero_tangent_vector(M, p);
                atol = ManifoldTests.find_eps(pts[1]),
            )
            Test.@test isapprox(M, p, 2 * X, X + X)
            Test.@test isapprox(M, p, 0 * X, X - X)
            Test.@test isapprox(M, p, (-1) * X, -X)
        end
    end

    test_tangent_vector_broadcasting &&
        Test.@testset "broadcasted linear algebra in tangent space" begin
            for (p, X) in zip(pts, tv)
                Test.@test isapprox(M, p, 3 * X, 2 .* X .+ X)
                Test.@test isapprox(M, p, -X, X .- 2 .* X)
                Test.@test isapprox(M, p, -X, .-X)
                if (isa(X, AbstractArray))
                    Y = allocate(X)
                    Y .= 2 .* X .+ X
                else
                    Y = 2 * X + X
                end
                Test.@test isapprox(M, p, Y, 3 * X)
            end
        end

    test_project_tangent && Test.@testset "project tangent test" begin
        for (p, X) in zip(pts, tv)
            atol = ManifoldTests.find_eps(p) * projection_atol_multiplier
            Test.@test isapprox(M, p, X, project(M, p, X); atol = atol)
            if is_mutating
                X2 = allocate(X)
                project!(M, X2, p, X)
            else
                X2 = project(M, p, X)
            end
            Test.@test isapprox(M, p, X2, X; atol = atol)
        end
    end

    test_project_point && Test.@testset "project point test" begin
        for p in pts
            atol = ManifoldTests.find_eps(p) * projection_atol_multiplier
            Test.@test isapprox(M, p, project(M, p); atol = atol)
            if is_mutating
                p2 = allocate(p)
                project!(M, p2, p)
            else
                p2 = project(M, p)
            end
            Test.@test isapprox(M, p2, p; atol = atol)
        end
    end

    !(
        default_retraction_method === nothing ||
        default_inverse_retraction_method === nothing
    ) && Test.@testset "vector transport" begin
        tvatol = is_tangent_atol_multiplier * ManifoldTests.find_eps(pts[1])
        X1 = inverse_retract(M, pts[1], pts[2], default_inverse_retraction_method)
        X2 = inverse_retract(M, pts[1], pts[3], default_inverse_retraction_method)
        pts32 = retract(M, pts[1], X2, default_retraction_method)
        test_default_vector_transport && Test.@testset "default vector transport" begin
            v1t1 = vector_transport_to(M, pts[1], X1, pts32)
            v1t2 = vector_transport_direction(M, pts[1], X1, X2)
            Test.@test is_tangent_vector(M, pts32, v1t1; atol = tvatol)
            Test.@test is_tangent_vector(M, pts32, v1t2; atol = tvatol)
            Test.@test isapprox(M, pts32, v1t1, v1t2)
            Test.@test isapprox(M, pts[1], vector_transport_to(M, pts[1], X1, pts[1]), X1)

            is_mutating && Test.@testset "mutating variants" begin
                v1t1_m = allocate(v1t1)
                v1t2_m = allocate(v1t2)
                vector_transport_to!(M, v1t1_m, pts[1], X1, pts32)
                vector_transport_direction!(M, v1t2_m, pts[1], X1, X2)
                Test.@test isapprox(M, pts32, v1t1, v1t1_m)
                Test.@test isapprox(M, pts32, v1t2, v1t2_m)
            end
        end

        for vtm in vector_transport_methods
            Test.@testset "vector transport method $(vtm)" begin
                tvatol = is_tangent_atol_multiplier * ManifoldTests.find_eps(pts[1])
                X1 = inverse_retract(M, pts[1], pts[2], default_inverse_retraction_method)
                X2 = inverse_retract(M, pts[1], pts[3], default_inverse_retraction_method)
                pts32 = retract(M, pts[1], X2, default_retraction_method)
                v1t1 = vector_transport_to(M, pts[1], X1, pts32, vtm)
                v1t2 = vector_transport_direction(M, pts[1], X1, X2, vtm)
                Test.@test is_tangent_vector(M, pts32, v1t1; atol = tvatol)
                Test.@test is_tangent_vector(M, pts32, v1t2; atol = tvatol)
                Test.@test isapprox(M, pts32, v1t1, v1t2)
                Test.@test isapprox(
                    M,
                    pts[1],
                    vector_transport_to(M, pts[1], X1, pts[1], vtm),
                    X1,
                )

                is_mutating && Test.@testset "mutating variants" begin
                    v1t1_m = allocate(v1t1)
                    v1t2_m = allocate(v1t2)
                    vector_transport_to!(M, v1t1_m, pts[1], X1, pts32, vtm)
                    vector_transport_direction!(M, v1t2_m, pts[1], X1, X2, vtm)
                    Test.@test isapprox(M, pts32, v1t1, v1t1_m)
                    Test.@test isapprox(M, pts32, v1t2, v1t2_m)
                end
            end
        end
    end

    for btype in basis_types_vecs
        Test.@testset "Basis support for $(btype)" begin
            p = pts[1]
            b = get_basis(M, p, btype)
            Test.@test isa(b, CachedBasis)
            bvectors = get_vectors(M, p, b)
            N = length(bvectors)

            # test orthonormality
            for i in 1:N
                Test.@test norm(M, p, bvectors[i]) ≈ 1
                for j in (i + 1):N
                    Test.@test real(inner(M, p, bvectors[i], bvectors[j])) ≈ 0 atol =
                        sqrt(ManifoldTests.find_eps(p))
                end
            end
            if isa(btype, ProjectedOrthonormalBasis)
                # check projection idempotency
                for i in 1:N
                    Test.@test norm(M, p, bvectors[i]) ≈ 1
                    for j in (i + 1):N
                        Test.@test real(inner(M, p, bvectors[i], bvectors[j])) ≈ 0 atol =
                            sqrt(ManifoldTests.find_eps(p))
                    end
                end
                # check projection idempotency
                for i in 1:N
                    Test.@test isapprox(M, p, project(M, p, bvectors[i]), bvectors[i])
                end
            end
            if !isa(btype, ProjectedOrthonormalBasis) && (
                basis_has_specialized_diagonalizing_get ||
                !isa(btype, DiagonalizingOrthonormalBasis)
            )

                X1 = inverse_retract(M, p, pts[2], default_inverse_retraction_method)
                Xb = get_coordinates(M, p, X1, btype)

                Test.@test get_coordinates(M, p, X1, b) ≈ Xb
                Test.@test isapprox(
                    M,
                    p,
                    get_vector(M, p, Xb, b),
                    get_vector(M, p, Xb, btype),
                )
            end
        end
    end

    for btype in (basis_types_to_from..., basis_types_vecs...)
        p = pts[1]
        N = number_of_coordinates(M, btype)
        if !isa(btype, ProjectedOrthonormalBasis) && (
            basis_has_specialized_diagonalizing_get ||
            !isa(btype, DiagonalizingOrthonormalBasis)
        )

            X1 = inverse_retract(M, p, pts[2], default_inverse_retraction_method)

            Xb = get_coordinates(M, p, X1, btype)
            #Test.@test isa(Xb, AbstractVector{<:Real})
            Test.@test length(Xb) == N
            Xbi = get_vector(M, p, Xb, btype)
            Test.@test isapprox(M, p, X1, Xbi)

            Xs = [[ifelse(i == j, 1, 0) for j in 1:N] for i in 1:N]
            Xs_invs = [get_vector(M, p, Xu, btype) for Xu in Xs]
            # check orthonormality of inverse representation
            for i in 1:N
                Test.@test norm(M, p, Xs_invs[i]) ≈ 1 atol = ManifoldTests.find_eps(p)
                for j in (i + 1):N
                    Test.@test real(inner(M, p, Xs_invs[i], Xs_invs[j])) ≈ 0 atol =
                        sqrt(ManifoldTests.find_eps(p))
                end
            end

            if is_mutating
                Xb_s = allocate(Xb)
                Test.@test get_coordinates!(M, Xb_s, p, X1, btype) === Xb_s
                Test.@test isapprox(Xb_s, Xb; atol = ManifoldTests.find_eps(p))

                Xbi_s = allocate(Xbi)
                Test.@test get_vector!(M, Xbi_s, p, Xb, btype) === Xbi_s
                Test.@test isapprox(M, p, X1, Xbi_s)
            end
        end
    end

    test_vee_hat && Test.@testset "vee and hat" begin
        p = pts[1]
        q = pts[2]
        X = inverse_retract(M, p, q, default_inverse_retraction_method)
        Y = vee(M, p, X)
        Test.@test length(Y) == number_of_coordinates(M, ManifoldsBase.VeeOrthogonalBasis())
        Test.@test isapprox(M, p, X, hat(M, p, Y))
        Y2 = allocate(Y)
        vee_ret = vee!(M, Y2, p, X)
        Test.@test vee_ret === Y2
        Test.@test isapprox(Y, Y2)
        X2 = allocate(X)
        hat_ret = hat!(M, X2, p, Y)
        Test.@test hat_ret === X2
        Test.@test isapprox(M, p, X2, X)
    end

    mid_point12 !== nothing && Test.@testset "midpoint" begin
        epsp1p2 = ManifoldTests.find_eps(pts[1], pts[2])
        atolp1p2 = exp_log_atol_multiplier * epsp1p2
        rtolp1p2 =
            exp_log_atol_multiplier == 0.0 ? sqrt(epsp1p2) * exp_log_rtol_multiplier : 0
        mp = mid_point(M, pts[1], pts[2])
        Test.@test isapprox(M, mp, mid_point12; atol = atolp1p2, rtol = rtolp1p2)
        if is_mutating
            mpm = allocate(mp)
            mid_point!(M, mpm, pts[1], pts[2])
            Test.@test isapprox(M, mpm, mid_point12; atol = atolp1p2, rtol = rtolp1p2)
        end
    end

    test_forward_diff && Test.@testset "ForwardDiff support" begin
        ManifoldTests.test_forwarddiff(M, pts, tv)
    end

    test_reverse_diff && Test.@testset "ReverseDiff support" begin
        ManifoldTests.test_reversediff(M, pts, tv)
    end

    test_musical_isomorphisms && Test.@testset "Musical isomorphisms" begin
        if default_inverse_retraction_method !== nothing
            tv_m = inverse_retract(M, pts[1], pts[2], default_inverse_retraction_method)
        else
            tv_m = zero_tangent_vector(M, pts[1])
        end
        ctv_m = flat(M, pts[1], FVector(TangentSpace, tv_m))
        Test.@test ctv_m.type == CotangentSpace
        tv_m_back = sharp(M, pts[1], ctv_m)
        Test.@test tv_m_back.type == TangentSpace
    end

    Test.@testset "number_eltype" begin
        for (p, X) in zip(pts, tv)
            Test.@test number_eltype(X) == number_eltype(p)
            p = retract(M, p, X, default_retraction_method)
            Test.@test number_eltype(p) == number_eltype(p)
        end
    end

    is_mutating && Test.@testset "copyto!" begin
        for (p, X) in zip(pts, tv)
            p2 = allocate(p)
            copyto!(p2, p)
            Test.@test isapprox(M, p2, p)

            X2 = allocate(X)
            if default_inverse_retraction_method === nothing
                X3 = zero_tangent_vector(M, p)
                copyto!(X2, X3)
                Test.@test isapprox(M, p, X2, zero_tangent_vector(M, p))
            else
                q = retract(M, p, X, default_retraction_method)
                X3 = inverse_retract(M, p, q, default_inverse_retraction_method)
                copyto!(X2, X3)
                Test.@test isapprox(
                    M,
                    p,
                    X2,
                    inverse_retract(M, p, q, default_inverse_retraction_method),
                )
            end
        end
    end

    is_mutating && Test.@testset "point distributions" begin
        for p in pts
            prand = allocate(p)
            for pd in point_distributions
                for _ in 1:10
                    Test.@test is_manifold_point(M, rand(pd))
                    if test_mutating_rand
                        rand!(pd, prand)
                        Test.@test is_manifold_point(M, prand)
                    end
                end
            end
        end
    end

    Test.@testset "tangent vector distributions" begin
        for tvd in tvector_distributions
            supp = Manifolds.support(tvd)
            for _ in 1:10
                randtv = rand(tvd)
                atol = rand_tvector_atol_multiplier * ManifoldTests.find_eps(randtv)
                Test.@test is_tangent_vector(M, supp.point, randtv; atol = atol)
            end
        end
    end
    return nothing
end


text:
module IndexFund

using ..Diderot
using ..Diderot: Arc, Node, Layer

struct Instance
    weights::Vector{Float64}
    similarity::Matrix{Float64}
    number::Int

    function Instance(weights, similarity, number)
        @assert length(weights) == size(similarity, 1)
        @assert size(similarity, 1) == size(similarity, 2)
        @assert number < length(weights)
        return new(weights, similarity, number)
    end
end

# Choose a representative for each stock.
Base.length(instance::Instance) = length(instance.weights)
Diderot.domain_type(instance::Instance) = Int

# Maximize weighted sum of similarities.
Diderot.value_type(instance::Instance) = Float64

# State encodes the stocks already picked (limited cardinality).
Diderot.initial_state(instance::Instance) = BitSet()

function Diderot.transitions(instance::Instance, state, variable)
    results = Dict{Arc{BitSet, Int, Float64}, BitSet}()

    N = length(instance)
    weights = instance.weights[variable]
    similarity = instance.similarity[:, variable]

    # Could we add another stock to the fund?
    if length(state) < instance.number
        for repr in setdiff(1:N, state)
            value = weights * similarity[repr]
            results[Arc(state, repr, value)] = union(state, repr)
        end
    end

    # Pick most similar stock already in the fund.
    if !isempty(state)
        best_repr, best_similarity = 0, -Inf
        for repr in state
            if similarity[repr] > best_similarity
                best_repr = repr
                best_similarity = similarity[repr]
            end
        end
        value = weights * best_similarity
        results[Arc(state, best_repr, value)] = state
    end

    return results
end

# For relaxation:
#  - Merge states by intersection.
merge = intersect
#  - Keep smallest states.
sort_by(pair) = length(pair.first)

relax(max_width) = Diderot.RelaxAllInOne(max_width, intersect, sort_by)

end