text:
"""
Central type for routing. Not meant for extending.

API
---
* [`allroutes`](@ref)
* [`getroute`](@ref)
* [`alllinks`](@ref)
* [`getlink`](@ref)
* [`start_vertices`](@ref)
* [`stop_vertices`](@ref)
* [`getchannel`](@ref)
* [`getmap`](@ref)
* [`getgraph`](@ref)
* [`iscongested`](@ref)
* [`clear_route`](@ref)
* [`setroute`](@ref)
"""
struct RoutingStruct{L<:RoutingLink,C<:RoutingChannel}
    "Base [`RoutingGraph`](@ref) for the underlying routing architecture."
    architecture_graph::RoutingGraph

    """
    Annotating [`RoutingLink`](@ref) for each routing element in
    `architecture_graph`.
    """
    graph_vertex_annotations::Vector{L}

    """
    Graphs of routing resources used by each channel.
    """
    routings::Vector{SparseDiGraph{Int}}

    """
    `Vector{RoutingChannel}` containing the channel information for the
    taskgraph.
    """
    channels::Vector{C}

    """
    Convenience structure mapping local channel indices back to edge indices
    in the parent taskgraph.
    """
    channel_index_to_taskgraph_index::Dict{Int,Int}
end

function RoutingStruct(map::Map)
    # Unpack some fields from the map
    toplevel = map.toplevel
    taskgraph = map.taskgraph
    @debug "Building Resource Graph"

    architecture_graph = routing_graph(toplevel)
    # Annotate the links in the routing graph with the custom structure defined
    # by the architecture type.
    graph_vertex_annotations = annotate(toplevel, rules(map), architecture_graph)
    # Get start and stop nodes for each taskgraph.
    channels, channel_dict = build_routing_taskgraph(map, architecture_graph)
    # Initialize the paths variable.
    routings = [SparseDiGraph{Int}() for i in 1:length(channels)]

    return RoutingStruct(
        architecture_graph, graph_vertex_annotations, routings, channels, channel_dict
    )
end

#---------#
# Methods #
#---------#
"""
$(SIGNATURES)

Return all routings in `routing_struct`.
"""
allroutes(routing_struct::RoutingStruct) = routing_struct.routings

"""
$(SIGNATURES)

Return route for channel `index`.
"""
getroute(routing_struct::RoutingStruct, index) = routing_struct.routings[index]

"""
$(SIGNATURES)

Return `Vector{RoutingLink}` for all links in `routing_struct.
"""
alllinks(routing_struct::RoutingStruct) = routing_struct.graph_vertex_annotations

"""
$(SIGNATURES)

Return `<:RoutingLink` for link in `routing_struct` with indes `i`.
"""
function getlink(routing_struct::RoutingStruct, i::Integer)
    return routing_struct.graph_vertex_annotations[i]
end

"""
$(SIGNATURES)

Return `Vector{PortVertices}` of start vertices for channel index `i`.

    start_vertices(channel::RoutingChannel) :: Vector{PortVertices}

Return [`Vector{PortVertices}`](@ref PortVertices) of start vertices for `channel`.
"""
function start_vertices(routing_struct::RoutingStruct, i::ChannelIndex)
    return start_vertices(routing_struct.channels[i])
end

"""
$(SIGNATURES)

Return `Vector{PortVertices}` of stop vertices for channel index `i`.

    stop_vertices(channel::RoutingChannel) :: Vector{PortVertices}

Return [`Vector{PortVertices}`](@ref PortVertices) of stop vertices for `channel`.
"""
function stop_vertices(routing_struct::RoutingStruct, i::ChannelIndex)
    return stop_vertices(routing_struct.channels[i])
end

"""
$(SIGNATURES)

Return `<:RoutingChannel` with indesx `i`.
"""
getchannel(routing_struct::RoutingStruct, i::ChannelIndex) = routing_struct.channels[i]

getmap(routing_struct::RoutingStruct) = getmap(routing_struct.architecture_graph)

"""
$(SIGNATURES)

Return the [`RoutingGraph`](@ref) member of `routing_struct`.
"""
getgraph(routing_struct::RoutingStruct) = routing_struct.architecture_graph

"""
    iscongested(routing_struct, [path])

Return `true` if routing congestion exists in `routing_struct`. If path is given
either as a [`ChannelIndex`](@ref) or `SparseDiGraph`, return `true` if just
the specified path is congested.

Method List
-----------
$(METHODLIST)
"""
function iscongested(routing_struct::RoutingStruct)
    return iscongested(routing_struct.graph_vertex_annotations)
end

function iscongested(routing_struct::RoutingStruct, path)
    return iscongested(routing_struct, getroute(routing_struct, path))
end

function iscongested(routing_struct::RoutingStruct, graph::SparseDiGraph{Int})
    for i in vertices(graph)
        iscongested(getlink(routing_struct, i)) && return true
    end
    return false
end

"""
    clear_route(rs::RoutingStruct, channel::ChannelIndex)

Rip up the current routing for the given link.
"""
function clear_route(rs::RoutingStruct, channel::ChannelIndex)
    #=
    1. Get the path for the link.
    2. Step through each architecture link on that path, remove the link index
        from that link info.
    3. Set the path to an empty set.
    =#
    path = getroute(rs, channel)
    for i in vertices(path)
        remchannel(getlink(rs, i), channel)
    end
    # Clear the path variable
    rs.routings[channel] = SparseDiGraph{Int}()
    return nothing
end

"""
$(SIGNATURES)

Assign `route` to `channel`.
"""
function setroute(
    routing_struct::RoutingStruct, route::SparseDiGraph, channel::ChannelIndex
)
    # This should always be the case - this assertion is to catch bugs.
    @assert nv(getroute(routing_struct, channel)) == 0
    for i in vertices(route)
        addchannel(getlink(routing_struct, i), channel)
    end
    return routing_struct.routings[channel] = route
end

function record(m::Map, r::RoutingStruct)
    # Get the dictionary mapping channel indices back to taskgraph edge indices.
    channel_dict = r.channel_index_to_taskgraph_index
    routes = translate_routes(getgraph(r), allroutes(r))
    for (i, route) in enumerate(routes)
        m.mapping[channel_dict[i]] = route
    end
end

# function get_routing_path(toplevel::TopLevel, g, maprev)
#     routing_path = SparseDiGraph{Any}()
#     # Add vertices
#     for v in vertices(g)
#         path = maprev[v]
#         add_vertex!(routing_path, path)
#     end
#     # Add edges
#     for src in vertices(g), dst in outneighbors(g, src)
#         add_edge!(routing_path, maprev[src], maprev[dst])
#     end
#     return routing_path
# end

################################################################################
# Routing Link Annotation
################################################################################
function annotate(toplevel::TopLevel, ruleset::RuleSet, rg::RoutingGraph)
    @debug "Annotating Graph Links"
    maprev = rev_dict(rg.map)

    # Initialize this to any. We'll clean it up later.
    routing_links_any = Vector{Any}(undef, nv(rg.graph))

    for (i, path) in maprev
        routing_links_any[i] = annotate(ruleset, toplevel[path])
    end
    # Clean up types
    routing_links = typeunion(routing_links_any)

    @debug "Type of Routing Link Annotations: $(typeof(routing_links))"
    return routing_links
end

################################################################################
# Routing Taskgraph Constructor
################################################################################
function build_routing_taskgraph(m::Map, r::RoutingGraph)
    # Debug printing
    @debug "Building Default Routing Taskgraph"
    # Unpack map
    taskgraph = m.taskgraph
    toplevel = m.toplevel
    ruleset = rules(m)

    # Get the list of channels that need routing.
    edges = getedges(taskgraph)
    edge_indices_to_route = [i for (i, e) in enumerate(edges) if needsrouting(ruleset, e)]

    # Create routing channels for all edges that need routing.
    channels = map(edge_indices_to_route) do index
        # Unpack the edge
        edge = edges[index]
        # Get source and destination nodes paths
        #
        # Wrap map in "Ref" to get scalar broadcasting behavior.
        sources = MapperCore.getpath.(Ref(m), getsources(edge))
        sinks = MapperCore.getpath.(Ref(m), getsinks(edge))
        # Convert these to indices in the routing graph
        start = collect_nodes(toplevel, ruleset, r.map, edge, sources, MapperCore.Source)
        stop = collect_nodes(toplevel, ruleset, r.map, edge, sinks, MapperCore.Sink)
        # Build the routing channel type
        return routing_channel(ruleset, start, stop, edge)
    end

    # Create a dictionary mapping indices in the "channels" vector to indices
    # in the original vector of edges.
    channel_dict = Dict(i => edge_indices_to_route[i] for i in 1:length(channels))

    # Return the collection of channels
    return channels, channel_dict
end

function collect_nodes(
    toplevel::TopLevel, ruleset::RuleSet, pathmap, edge::TaskgraphEdge, paths, dir
)
    nodes = Vector{PortVertices}()
    # Iterate through the source paths - get the port names.
    for path in paths
        # Get the component from the architecture
        component = toplevel[path]
        ports = get_routing_ports(ruleset, edge, component, dir)

        # The ports are just a collection of strings.
        #
        # 1. Create port path types from these ports using the path for the
        #    component that they belong to.
        #
        # 2. Use these full paths to index into the portmap dictionary to
        #    get the numbers in the routing graph.
        port_paths = [catpath(path, Path{Port}(port)) for port in ports]
        port_indices = PortVertices([pathmap[pp] for pp in port_paths])

        # Add this to the collection of nodes
        push!(nodes, port_indices)
    end
    return nodes
end

# Get vector of ports of component "c" that can serve as source/sink for the
# taskgraph edge.
function get_routing_ports(
    ruleset::RuleSet, edge::TaskgraphEdge, component::Component, dir::MapperCore.Direction
)
    if dir == MapperCore.Source
        return get_routing_ports(ruleset, edge, component, dir, is_source_port)
    else
        return get_routing_ports(ruleset, edge, component, dir, is_sink_port)
    end
end

function get_routing_ports(ruleset, edge, component, dir, f::Function)
    return [
        path for (path, port) in component.ports if
        checkclass(invert(port), dir) && f(ruleset, port, edge)
    ]
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

include("bat_sample.jl")
include("mcmc/mcmc.jl")
include("sampled_density.jl")
include("importance/importance_sampler.jl")
include("partitioned_sampling/partitioned_sampling.jl")

ext: jl


text:
module FileTools

using DataFrames
using Printf
using Dates
import Base.parse
using DelimitedFiles

include("tsffile.jl")
include("asciifile.jl")
include("campbellfile.jl")
include("dygraphsfile.jl")
include("eopfile.jl")
include("atmacsfile.jl")
include("ggpfile.jl")
include("gpcpdfile.jl")
include("hydrusfile.jl")
include("baytapfile.jl")
include("igetsfile.jl")
include("otherfile.jl")
include("dwdfile.jl")
include("tidefile.jl")
include("gravityeffectfile.jl")
include("et0calc.jl")

export readtsf, writetsf, loadascii, writeascii, readcampbell
export readdygraphs, writedygraphs
export loadeop
export stackframes, loadatmacs
export readggp, writeggp, ggpdata2blocks
export readgpcpd, readgpcpd_head, readgpcpd_time, readgpcpd_lonlat, coorShift_lon
export writeatmosph, writeprofile1d, readhydrus1d_obsnode, readhydrus1d_nodinf, readatmosph
export writebaytap, baytap2tsoft
export igetsimport, igetsexport
export dwdclimateraw, readgssd, convertgssd
export vav2tsoft, eterna2tsoft
export read_layerResponse
export ET0calc_dsc, ET0calc_dta, ET0calc_read

end #module


text:
import HTTP
import Markdown: htmlesc
import UUIDs: UUID
import JSON

# Serve everything from `/frontend`, and create HTTP endpoints to open notebooks.

to_symbol_dict(d) = Dict(Symbol(k) => v for (k, v) in d)

function build_file_tree(session, root_directory, parent_id="")
    tree = Dict()

    tree["id"] = parent_id * (parent_id == "" ? "" : "-") * basename(root_directory)
    tree["name"] = basename(root_directory)
    tree["type"] = "directory"
    tree["children"] = []

    files = []
    dirs = []
    for node = readdir(root_directory)
        nodepath = joinpath(root_directory, node)
        if isfile(nodepath)
            push!(files, node)
        else
            push!(dirs, node)
        end
    end

    for dir = dirs
        push!(tree["children"], build_file_tree(session, joinpath(root_directory, dir), tree["id"]))
    end
    for file = files
        push!(tree["children"], Dict(
            "id" => tree["id"] * "-" * makeid(8),
            "name" => file,
            "type" => "file",
            "running" => joinpath(root_directory, file) ∈ (x -> x.path).(values(session.notebooks))
        ))
    end

    tree
end

to_notebook_root_relative(session::ServerSession, path::String) = replace(replace(path, session.options.server.notebook_root => ""), "\\" => "/")

"Attempts to find the MIME pair corresponding to the extension of a filename. Defaults to `text/plain`."
function mime_fromfilename(filename)
    # This bad boy is from: https://developer.mozilla.org/en-US/docs/Web/HTTP/Basics_of_HTTP/MIME_types/Common_types
    mimepairs = Dict(".aac" => "audio/aac", ".bin" => "application/octet-stream", ".bmp" => "image/bmp", ".css" => "text/css", ".csv" => "text/csv", ".eot" => "application/vnd.ms-fontobject", ".gz" => "application/gzip", ".gif" => "image/gif", ".htm" => "text/html", ".html" => "text/html", ".ico" => "image/vnd.microsoft.icon", ".jpeg" => "image/jpeg", ".jpg" => "image/jpeg", ".js" => "text/javascript", ".json" => "application/json", ".jsonld" => "application/ld+json", ".mjs" => "text/javascript", ".mp3" => "audio/mpeg", ".mp4" => "video/mp4", ".mpeg" => "video/mpeg", ".oga" => "audio/ogg", ".ogv" => "video/ogg", ".ogx" => "application/ogg", ".opus" => "audio/opus", ".otf" => "font/otf", ".png" => "image/png", ".pdf" => "application/pdf", ".rtf" => "application/rtf", ".sh" => "application/x-sh", ".svg" => "image/svg+xml", ".tar" => "application/x-tar", ".tif" => "image/tiff", ".tiff" => "image/tiff", ".ttf" => "font/ttf", ".txt" => "text/plain", ".wav" => "audio/wav", ".weba" => "audio/webm", ".webm" => "video/webm", ".webp" => "image/webp", ".woff" => "font/woff", ".woff2" => "font/woff2", ".xhtml" => "application/xhtml+xml", ".xml" => "application/xml", ".xul" => "application/vnd.mozilla.xul+xml", ".zip" => "application/zip")
    file_extension = getkey(mimepairs, '.' * split(filename, '.')[end], ".txt")
    MIME(mimepairs[file_extension])
end

trial_extensions = ["html", "js", "css"]
function asset_response(path)
    default_asset = project_relative_path("frontend/lab", "index.html")
    if !isfile(path)
        for ext = trial_extensions
            ext_path = path * "." * ext
            if isfile(ext_path)
                return asset_response(ext_path)
            end
        end
        return asset_response(default_asset)
    end
    try
        @assert isfile(path)
        response = HTTP.Response(200, read(path, String))
        push!(response.headers, "Content-Type" => string(mime_fromfilename(path)))
        push!(response.headers, "Access-Control-Allow-Origin" => "*")
        response
    catch e
        asset_response(default_asset)
    end
end

function json_response(obj)
    response = HTTP.Response(200, JSON.json(obj));
    push!(response.headers, "Content-Type" => "text/json")
    response
end

function error_response(status_code::Integer, title, advice, body="")
    template = read(project_relative_path("frontend", "error.jl.html"), String)

    body_title = body == "" ? "" : "Error message:"
    filled_in = replace(replace(replace(replace(template, "\$TITLE" => title), "\$ADVICE" => advice), "\$BODYTITLE" => body_title), "\$BODY" => htmlesc(body))

    response = HTTP.Response(status_code, filled_in)
    push!(response.headers, "Content-Type" => string(mime_fromfilename(".html")))
    response
end

function notebook_redirect_response(notebook; home_url="./")
    response = HTTP.Response(302, "")
    push!(response.headers, "Location" => home_url * "edit?id=" * string(notebook.notebook_id))
    return response
end

"""
Return whether the `request` was authenticated in one of two ways:
1. the session's `secret` was included in the URL as a search parameter, or
2. the session's `secret` was included in a cookie.
"""
function is_authenticated(session::ServerSession, request::HTTP.Request)
    (
        secret_in_url = try
            uri = HTTP.URI(request.target)
            query = HTTP.queryparams(uri)
            get(query, "secret", "") == session.secret
        catch e
            @warn "Failed to authenticate request using URL" exception = (e, catch_backtrace())
            false
        end
    ) || (
        secret_in_cookie = try
            cookies = HTTP.cookies(request)
            any(cookies) do cookie
                cookie.name == "secret" && cookie.value == session.secret
            end
        catch e
            @warn "Failed to authenticate request using cookies" exception = (e, catch_backtrace())
            false
        end
    )
    # that ) || ( kind of looks like Krabs from spongebob
end

function http_router_for(session::ServerSession)
    router = HTTP.Router()
    security = session.options.security

    function add_set_secret_cookie!(response::HTTP.Response)
        push!(response.headers, "Set-Cookie" => "secret=$(session.secret); SameSite=Strict; HttpOnly")
        response
    end

    """
        with_authentication(f::Function)

    Returns a function `HTTP.Request → HTTP.Response` which does three things:
    1. Check whether the request is authenticated (by calling `is_authenticated`), if not, return a 403 error.
    2. Call your `f(request)` to create the response message.
    3. Add a `Set-Cookie` header to the response with the session's `secret`.
    """
    function with_authentication(f::Function; required::Bool)::Function
        return function (request::HTTP.Request)
            if !required || is_authenticated(session, request)
                response = f(request)
                add_set_secret_cookie!(response)
                response
            else
                error_response(403, "Not yet authenticated", "Open the link that was printed in the terminal where you launched Pluto. It includes a secret, which is needed to access this server.

If you are running the server yourself and want to change this configuration, have a look at the keyword arguments to Pluto.run.

Please report this error if you did not expect it!")
            end
        end
    end

    function create_serve_onefile(path)
        return request::HTTP.Request -> asset_response(normpath(path))
    end

    # / does not need security.require_secret_for_open_links
    # because this is how we handle the case:
    #    require_secret_for_open_links == true
    #    require_secret_for_access == false
    # Access to all 'risky' endpoints is still restricted to requests that have the secret cookie, but visiting `/` is allowed, and it will set the cookie. From then on the security situation is identical to
    #    secret_for_access == true
    HTTP.@register(router, "GET", "/", with_authentication(
        create_serve_onefile(project_relative_path("frontend/lab", "index.html"));
        required=security.require_secret_for_access
        ))
    HTTP.@register(router, "GET", "/edit", with_authentication(
        create_serve_onefile(project_relative_path("frontend", "editor.html"));
        required=security.require_secret_for_access ||
        security.require_secret_for_open_links,
    ))
    # the /edit page also uses with_authentication, but this is not how access to notebooks is secured: this is done by requiring the WS connection to be authenticated.
    # we still use it for /edit to do the cookie stuff, and show a more helpful error, instead of the WS never connecting.

    HTTP.@register(router, "GET", "/ping", r -> HTTP.Response(200, "OK!"))
    HTTP.@register(router, "GET", "/possible_binder_token_please", r -> HTTP.Response(200, session.binder_token === nothing ? "" : session.binder_token))

    function try_launch_notebook_response(action::Function, path_or_url::AbstractString; title="", advice="", home_url="./")
        try
            nb = action(session, path_or_url)
            notebook_redirect_response(nb; home_url=home_url)
        catch e
            if e isa SessionActions.NotebookIsRunningException
                notebook_redirect_response(e.notebook; home_url=home_url)
            else
                error_response(500, title, advice, sprint(showerror, e, stacktrace(catch_backtrace())))
            end
        end
    end

    serve_newfile = with_authentication(;
        required=security.require_secret_for_access ||
        security.require_secret_for_open_links
    ) do request::HTTP.Request
        notebook = SessionActions.new(session; to_symbol_dict(HTTP.queryparams(HTTP.URI(request.target)))...)
        json_response(Dict(
            "id" => string(notebook.notebook_id),
            "path" => to_notebook_root_relative(session, notebook.path)
        ))
    end
    HTTP.@register(router, "GET", "/new", serve_newfile)

    serve_openfile = with_authentication(;
        required=security.require_secret_for_access ||
        security.require_secret_for_open_links
    ) do request::HTTP.Request
        try
            uri = HTTP.URI(request.target)
            query = HTTP.queryparams(uri)
            if haskey(query, "path")
                path = joinpath(session.options.server.notebook_root, query["path"])
                if isfile(path)
                    return try_launch_notebook_response(SessionActions.open, path, title="Failed to load notebook", advice="The file $(htmlesc(path)) could not be loaded. Please report this error!")
                else
                    return error_response(404, "Can't find a file here", "Please check whether $(htmlesc(path)) exists.")
                end
            elseif haskey(query, "url")
                url = query["url"]
                return try_launch_notebook_response(SessionActions.open_url, url, title="Failed to load notebook", advice="The notebook from $(htmlesc(url)) could not be loaded. Please report this error!")
            else
                error("Empty request")
            end
        catch e
            return error_response(400, "Bad query", "Please report this error!", sprint(showerror, e, stacktrace(catch_backtrace())))
        end
    end

    HTTP.@register(router, "GET", "/open", serve_openfile)

    serve_sample = with_authentication(;
        required=security.require_secret_for_access ||
        security.require_secret_for_open_links
    ) do request::HTTP.Request
        uri = HTTP.URI(request.target)
        sample_path = split(HTTP.unescapeuri(uri.path), "sample/")[2]
        sample_path_without_dotjl = "sample " * sample_path[1:end - 3]

        path = numbered_until_new(joinpath(new_notebooks_directory(), sample_path_without_dotjl))
        readwrite(project_relative_path("sample", sample_path), path)

        try_launch_notebook_response(SessionActions.open, path, home_url="../", title="Failed to load sample", advice="Please report this error!")
    end
    HTTP.@register(router, "GET", "/sample/*", serve_sample)

    serve_notebookfile = with_authentication(;
        required=security.require_secret_for_access ||
        security.require_secret_for_open_links
    ) do request::HTTP.Request
        try
            uri = HTTP.URI(request.target)
            query = HTTP.queryparams(uri)
            id = UUID(query["id"])
            notebook = session.notebooks[id]

            response = HTTP.Response(200, sprint(save_notebook, notebook))
            push!(response.headers, "Content-Type" => "text/plain; charset=utf-8")
            push!(response.headers, "Content-Disposition" => "inline; filename=\"$(basename(notebook.path))\"")
            response
        catch e
            return error_response(400, "Bad query", "Please report this error!", sprint(showerror, e, stacktrace(catch_backtrace())))
        end
    end
    HTTP.@register(router, "GET", "/notebookfile", serve_notebookfile)

    serve_file_structure = with_authentication(;
        required=security.require_secret_for_access ||
        security.require_secret_for_open_links
    ) do request::HTTP.Request
        json_response(build_file_tree(session, session.options.server.notebook_root))
    end
    HTTP.@register(router, "GET", "/tree", serve_file_structure)

    serve_file_move = with_authentication(;
        required=security.require_secret_for_access ||
        security.require_secret_for_open_links
    ) do request::HTTP.Request
        params = HTTP.queryparams(HTTP.URI(request.target))
        success_response() = json_response(Dict("success" => true))
        error_response(msg) = json_response(Dict("success" => false, "error" => msg))
        try
            src_path = normpath(joinpath(session.options.server.notebook_root, params["src"]))
            dst_path = normpath(joinpath(session.options.server.notebook_root, params["dst"]))

            if joinpath(src_path, "") == joinpath(dst_path, "")
                return success_response()
            end

            notebook_paths = (x -> x.path).(values(session.notebooks))
            nb_src_occurs = occursin.(repeat(src_path, length(notebook_paths)), notebook_paths)
            if isdir(src_path) && length(nb_src_occurs) > 0 && sum(nb_src_occurs) > 0
                return error_response("Folders cannot be moved while a Pluto notebook is running inside it")
            end

            # TODO: Possibly use SessionActions.move(session, notebook, newpath)
            mv(src_path, dst_path)

            # Change notebook path to match destination path
            notebook_running = findall(x -> x == src_path, notebook_paths)
            if length(notebook_running) > 0
                nb = session.notebooks[collect(keys(session.notebooks))[notebook_running[1]]]
                nb.path = dst_path
            end

            success_response()
        catch e
            error_response(string(e))
        end

    end
    HTTP.@register(router, "GET", "/fileMove", serve_file_move)

    function serve_asset(request::HTTP.Request)
        uri = HTTP.URI(request.target)

        filepath = project_relative_path("frontend", relpath(HTTP.unescapeuri(uri.path), "/"))
        asset_response(filepath)
    end
    HTTP.@register(router, "GET", "/*", serve_asset)
    HTTP.@register(router, "GET", "/favicon.ico", create_serve_onefile(project_relative_path("frontend", "img", "favicon.ico")))

    return router
end


text:
##############################################################################################################
#                                      Load Packages
##############################################################################################################
cd(@__DIR__)
using Pkg
Pkg.activate("../../")
using Revise, TwoDots, TwoDotsModels
import TwoDotsModels: decide!
include("model_1.jl")
##############################################################################################################
#                                         Model Setup
##############################################################################################################
model = Model(;real_time = true)
game = Game()
gui = GUI(;game)
run!(model, game, gui)


text:
# basis.jl
#
# Types and methods for defining and working with basis vectors of vector spaces

# BasisVector and Basis
#-----------------------
# immutable type for specifying basis vectors: every vector space should have a canonical basis, where the basis
# vectors can be specified using an identifier of some custom type T (e.g. user defined)
immutable BasisVector{S<:VectorSpace,T}
    space::S
    identifier::T
end

space(b::BasisVector) = b.space
Base.in(b::BasisVector,V::VectorSpace) = space(b) == V

Base.show(io::IO, b::BasisVector) = print(io, "BasisVector($(b.space),$(b.identifier))")

# a basis is just an iterator over the canonical basis vectors of the space
immutable Basis{S<:VectorSpace}
    space::S
end

space(B::Basis) = B.space
basis(V::VectorSpace) = Basis(V) # lowercase convenience method

Base.show(io::IO, B::Basis) = print(io, "Basis($(B.space))")


text:
function initL{Tboundary <:BoundaryStrategy }(
  pf :: PartialD_OP{PartialD{1}, StencilsnN{3,3}, Tboundary},
  data :: Array{Float64,1}
  )
  a = zeros(data);
  b = zeros(data);
  c = zeros(data);
  a[2:end-1] = 2/3.0
  c[2:end-1] = 1/6.0
  b[2:end-1] = 1/6.0
  a[1] = 1.0; a[end] = 1.0;
  c[1] = 3.0; b[end] = 3.0;
  return(a,b,c)
end

function initL{Tboundary <:BoundaryStrategy }(
  pf :: PartialD_OP{PartialD{2}, StencilsnN{3,3}, Tboundary},
  data :: Array{Float64,1}
  )
  a = zeros(data);
  b = zeros(data);
  c = zeros(data);
  a[2:end-1] = 5/6.0
  c[2:end-1] = 1/12.0
  b[2:end-1] = 1/12.0
  a[1] = 1.0; a[end] = 1.0;
  c[1] = 11.0; b[end] = 11.0;
  return(a,b,c)
end

function initR{Tboundary <:BoundaryStrategy }(
  pf :: PartialD_OP{PartialD{1}, StencilsnN{3,3}, Tboundary},
  data :: Array{Float64,1}
  )
  res = zeros(data)
  res[2:end-1] = (data[3:end] - data[1:end-2])/2
  pboundary = [-17/6, 3/2 , 3/2, -1/6]
  res[1] = sum( pboundary .* data[1:4])
  res[end] =  - sum( pboundary .* data[end:-1:end-3])
  return res
end

function initR{Tboundary <:BoundaryStrategy }(
  pf :: PartialD_OP{PartialD{2}, StencilsnN{3,3}, Tboundary},
  data :: Array{Float64,1}
  )
  res = zeros(data)
  res[2:end-1] = data[3:end] + data[1:end-2] - 2*data[2:end-1]
  pboundary = [13, -27 , 15, -1]
  res[1] = sum( pboundary .* data[1:4])
  res[end] =   sum( pboundary .* data[end:-1:end-3])
  return res
end


text:
module LightFSM

export FiniteStateMachine,
  fire!,
  finished,
  can

type FiniteStateMachine{StateT<:Enum,EventT<:Enum}
  map::Vector{Vector{Pair{StateT,StateT}}}
  current::StateT
  terminal::StateT
  FiniteStateMachine(map,current,terminal) =
    new(map,current,terminal)
end

FiniteStateMachine{StateT<:Enum,EventT<:Enum}(::Type{StateT},::Type{EventT}) =
  FiniteStateMachine{StateT,EventT}(Vector{Vector{Pair{StateT,StateT}}}(),instances(StateT)[1],instances(StateT)[end])

function FiniteStateMachine{StateT<:Enum,EventT<:Enum}(map::Dict{EventT,Vector{Pair{StateT,StateT}}},
  current::StateT,terminal::StateT)
  fsm = FiniteStateMachine(StateT,EventT)
  m = Vector{Vector{Pair{StateT,StateT}}}(length(instances(EventT)))
  for kv in map
    m[Int(kv[1])] = copy(kv[2])
  end
  fsm.map = m
  fsm.current = current
  fsm.terminal = terminal
  return fsm
end

function fire!{StateT<:Enum,EventT<:Enum}(fsm::FiniteStateMachine{StateT,EventT},
  event::EventT)
  found = false
  if isdefined(fsm.map,Int(event))
    @inbounds for i in 1:length(fsm.map[Int(event)])
      if fsm.current == fsm.map[Int(event)][i][1]
        fsm.current = fsm.map[Int(event)][i][2]
        found = true
        break
      end
    end
  end
  return found
end

finished{StateT<:Enum,EventT<:Enum}(fsm::FiniteStateMachine{StateT,EventT}) =
  fsm.current == fsm.terminal

function can{StateT<:Enum,EventT<:Enum}(fsm::FiniteStateMachine{StateT,EventT},
  event::EventT)
  found = false
  if isdefined(fsm.map,Int(event))
    @inbounds for i in 1:length(fsm.map[Int(event)])
      if fsm.current == fsm.map[Int(event)][i][1]
        found = true
        break
      end
    end
  end
  return found
end

end # module


text:
# lbm that calculates eigenvectors/eigenvalues
# cuts off eigenvalues after a certain drop in orders of magnitude (input?)
using Plots
using LinearAlgebra
using SparseArrays
using Arpack

include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/experimental/LBM.jl")

function lattice_boltzmann_basis(xres, yres, f, N, tau, deltax, deltat, threshmag, count=0, rhofull=zeros(xres*yres, N))
    # returns the eigenvalues/eigenvectors that define the basis of the density function over the given time period
    c = deltax/deltat

    ftransform = zeros(xres*yres, 9)
    for j=1:yres
        for i=1:xres
            ftransform[(j-1)*xres+i,:] = f[i,j,:]
        end
    end

    stensor = PeriodicStreamTensor(xres, yres)
    fnew = stensor * ftransform

    rho = zeros((xres*yres))
    for i=1:xres*yres
        rho[i] = sum(fnew[i,:])
    end

    u = zeros(xres*yres, 2)
    for i=1:xres*yres
        ux = (fnew[i,2] - fnew[i,4] + 1/sqrt(2) * (fnew[i,6] + fnew[i,9] - fnew[i,7] - fnew[i,8])) / rho[i]
        uy = (fnew[i,3] - fnew[i,5] + 1/sqrt(2) * (fnew[i,6] + fnew[i,8] - fnew[i,9] - fnew[i,8])) / rho[i]
        u[i,1] = ux
        u[i,2] = uy
    end

    unitvecs = zeros(9,2)
    unitvecs[:,1] = [0, c, 0, -c, 0, 1/sqrt(2)*c, -1/sqrt(2)*c, -1/sqrt(2)*c, 1/sqrt(2)*c]
    unitvecs[:,2] = [0, 0, c, 0, -c, 1/sqrt(2)*c, 1/sqrt(2)*c, -1/sqrt(2)*c, -1/sqrt(2)*c]

    w = zeros(9)
    w[1] = 4/9
    w[2] = 1/9
    w[3] = 1/9
    w[4] = 1/9
    w[5] = 1/9
    w[6] = 1/36
    w[7] = 1/36
    w[8] = 1/36
    w[9] = 1/36

    s = zeros((xres*yres, 9))
    usquared = zeros((xres*yres))
    for i=1:xres*yres
        usquared[i] = dot(u[i,:],u[i,:])
        s[i,:] = w[:] .* (3*(unitvecs*u[i,:])/c^2 .+ 9/2*((unitvecs*u[i,:])).^2/c^3 .- 3/2*usquared[i]/c^2)
    end

    feq = zeros((xres*yres, 9))
    for i=1:xres*yres
        feq[i,:] = w[:].*rho[i] .+ rho[i].*s[i,:]
    end

    fnext = zeros((xres*yres, 9))
    for i=1:xres*yres
        fnext[i,:] = fnew[i,:] .- (deltat/(tau+.5*deltat)).*(fnew[i,:] - feq[i,:])
    end

    fn = zeros(xres, yres, 9)
    for j=1:yres
        for i=1:xres
            fn[i,j,:] = fnext[(j-1)*xres+i,:]
        end
    end

    if count == 0
        rhofull = rho
    else
        rhofull = hcat(rhofull, rho)
    end

    count = count + 1

    if count < N
        lattice_boltzmann_basis(xres, yres, fn, N, tau, deltax, deltat, threshmag, count, rhofull)
    else
        p = size(rhofull)[2]
        covarmatrix = 1/(p-1) * rhofull * transpose(rhofull)
        nevv = 100
        evals, evecs = eigs(covarmatrix, nev=nevv, which=:LM)

        for i=1:nevv-1
            if real(evals[i])/real(evals[1]) < threshmag
                evalsnew = zeros(i)
                evecsnew = zeros(xres*yres,i)
                evalsnew[:] .= evals[1:i,1]
                evecsnew[:,:] .= evecs[:,1:i]
                coeff = transpose(evecsnew) * rhofull
                return evalsnew, evecsnew, coeff
            end
        end

    end
end

xres = 100
yres = 100
N = 100
tau = 10000
deltax = 1
deltat = 1
threshmag = 10^-8

fsine = zeros(xres, yres, 9)
for i=1:xres
    for j=1:yres
        for k=1:9
            fsine[i,j,k] = sin(i*pi/xres)*sin(j*pi/yres)
        end
    end
end

evals, evecs, coefficients = lattice_boltzmann_basis(xres, yres, fsine, N, tau, deltax, deltat, threshmag)
println(size(evals))


text:
using Test, MPI, Impero, LinearAlgebra
import Impero: compute

using ClimateMachine
using ClimateMachine.MPIStateArrays

ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()
const mpicomm = MPI.COMM_WORLD

Q = MPIStateArray{Float64}(mpicomm, ArrayType, 4, 6, 8)

compute(Q::MPIStateArray) = Q.realdata
@wrapper q=Q
compute(2*q)
compute(2*q) - 2 * Q

@testset "Impero Algebra" begin
    @test norm(compute(2*q) - 2 * Q) == 0.0
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
#=
O módulo de datas oferece dois caminhos para lidar com datas: Date e Datetime. Um
representa o tempo com precisão diária, enquanto o outro representa o tempo com
precisão na casa dos milisegundos, respectivamente, mas ambos são subtipos do tipo
abstrato TimeType.

O motivo para a existência de dois tipos é simples: objetos Date são mais simples,
pois dispensam a necessidade de lidar com horas, minutos, segundos e todas as
considerações usuais, como fusohorário, horário de verão e afins. Ou seja: se
o necessário for apenas lidar com dias, anos ou meses, use Date. Se a granularidade
necessária formaior, use Datetime.

Ambos Date e Datetime são, essencialmente, wrappers imutáveis de objetos Int64.
O Datetime ignora timezones por natureza (naive, no vocabulário Python).
Funcionalidades envolvendo timezones podem ser adquiridas através do módulo
TimeZones.jl

Documentação original: https://docs.julialang.org/en/v1/stdlib/Dates/#Adjuster-Functions
=#

import Dates

# =============================================================
# Constructors
# =============================================================

Dates.DateTime(2021)
Dates.DateTime(2021, 1)
Dates.DateTime(2021, 1, 1)
Dates.DateTime(2021, 1, 1, 12)
Dates.DateTime(2021, 1, 1, 12, 30)
Dates.DateTime(2021, 1, 1, 12, 30, 59)
Dates.DateTime(2021, 1, 1, 12, 30, 59, 1)

Dates.Date(2021)
Dates.Date(2021, 1)
Dates.Date(2021, 1, 1)
Dates.Date(Dates.Year(2013), Dates.Month(7), Dates.Day(1))
Dates.Date(Dates.Month(7), Dates.Year(2013))

# Objetos também podem ser instanciados com strings
Dates.Date("2021-01-01")

# A biblioteca possui alguns constructors utilitários
Dates.now()
Dates.today()

# =============================================================
# Formatação
# =============================================================

df = Dates.DateFormat("y-m-d")
Dates.Date("2021-01-01", df)

# Também é possível formatar os objetos por broadcasting
anos = ["2020", "2021"]
Dates.Date.(anos, Dates.DateFormat("yyyy"))

# =============================================================
# Durações e comparações
# =============================================================

dt1 = Dates.Date(1900, 1, 1)
dt2 = Dates.Date(2000, 1, 2)

# subtrair Dates retornará a diferença em dias.
# subtrair Datetimes retornará a diferença em milisegundos.
dt1 - dt2
dt2 - dt1

dt1 < dt2
dt1 > dt2

# Nota: operações envolvendo +, * e / não são suportadas.

# =============================================================
# Accessor functions
# =============================================================

Dates.year(dt1)
Dates.month(dt1)
Dates.day(dt1)

# Funções de acesso composto são oferecidas por conveniência
Dates.yearmonth(dt1)
Dates.monthday(dt1)
Dates.yearmonthday(dt1)

# =============================================================
# Query functions
# Funções de query oferecem informações de calendário sobre objetos TimeType
# =============================================================

Dates.dayofweek(dt1)
Dates.dayname(dt1)
Dates.monthname(dt1)
Dates.daysinmonth(dt1)
Dates.dayofyear(dt1)
Dates.isleapyear(dt1)

# Dayname e monthname podem receber um parâmetro de locale opcional com nomes
# de dias e meses de outras linguagens/ locales, assim como as abreviações.

# namely dayabbr and monthabbr. First the mapping is loaded into the LOCALES variable:

meses = [
    "janeiro", "fevereiro", "março", "abriç", "maio", "junho",
    "julho", "agosto", "setembro", "outubro", "novembro", "dezembro"
];

abreviacao_meses = [
    "jan","feb","mar","abr","mai","jun",
    "jul","ago","set","out","nov","dec"
];

dias = ["segunda", "terça", "quarta", "quinta", "sexta", "sábado", "domingo"];

abreviacoa_dias = ["seg", "ter", "qua", "qui", "sex", "sáb", "dom"];

Dates.LOCALES["brasil"] = Dates.DateLocale(
    meses, abreviacao_meses, dias, abreviacoa_dias
);

Dates.dayabbr(dt1; locale="brasil")
Dates.dayname(dt1; locale="brasil")
Dates.monthabbr(dt1; locale="brasil")
Dates.monthname(dt1; locale="brasil")

# =============================================================
# Aritimética com datas
# =============================================================

# A princípio, a ordem das operações não afeta o resultado:
Dates.Date(2000, 1, 29) + Dates.Day(1) + Dates.Month(1)
Dates.Date(2000, 1, 29) + Dates.Month(1) + Dates.Day(1)

# Contudo, a ordem pode se tornar muito importante, se alguma ordem for forçada:
(Dates.Date(2000, 1, 29) + Dates.Day(1)) + Dates.Month(1)
(Dates.Date(2000, 1, 29) + Dates.Month(1)) + Dates.Day(1)

# O que está acontecendo no exemplo acima?
# Quando solicitamos que 1 mês seja adicionado à data, o valor do dia é verificado.
# Se o dia for maior que o máximo do mês, o dia é aproximado (para baixo)
# para o máximo do mês. Com esse exemplo, podemos perceber que o design da biblioteca
# é sempre somar os "tempos" ordenadamente: primeiro anos, depois meses, depois dias
# e assim por diante. Se nenhuma ordem for forçada, o comportamento
# será bastante canônico.

# As operações podem ser feitas com ranges!
dr1 = Dates.Date(2021, 1, 29):Dates.Day(1):Dates.Date(2021, 2, 3);
collect(dr1)

dr2 = Dates.Date(2021, 1, 29):Dates.Month(1):Dates.Date(2021, 07, 29);
collect(dr2)

# =============================================================
# Adjuster Functions
# Para auxiliar nas operações, a biblioteca oferece funções de ajuste
# que facilitam cálculos mais complexos.
# =============================================================

dt1 = Dates.Date(2000, 6, 3)

# ajusta a data pra a segunda-feira da data da semana inputada
Dates.firstdayofweek(dt1)

# ajusta para o último dia do mês da data inputada
Dates.lastdayofmonth(dt1)

# Retorna true se o dia da semana de X for terça
e_terca = x -> Dates.dayofweek(x) == Dates.Tuesday;
# # 2021-12-25 foi um sábado
Dates.tonext(e_terca, Dates.Date(2021, 12, 25))

# O mesmo acima, porém escrito de uma forma mais conveniente:
Dates.tonext(Dates.Date(2021, 12, 25), Dates.Tuesday)

# A formulação acima pode ser útil para expressões mais complexas:
Dates.tonext(Dates.Date(2021, 12, 25)) do x
    # Return true on the 4th Thursday of November (Thanksgiving)
    Dates.dayofweek(x) == Dates.Thursday &&
    Dates.dayofweekofmonth(x) == 4 &&
    Dates.month(x) == Dates.November
end

# A função filter do pacote Base pode ser usada para filtrar datas válidas
# dadas certas condições:
dr = Dates.Date(2021):Dates.Day(1):Dates.Date(2022);
filter(dr) do x
    Dates.dayofweek(x) == Dates.Tue &&
    Dates.April <= Dates.month(x) <= Dates.Nov &&
    Dates.dayofweekofmonth(x) == 2
end

max_stars_repo_name: ericmiguel/julia-notes

max_issues_repo_name: ericmiguel/julia-notes

max_forks_repo_name: ericmiguel/julia-notes


text:
file = joinpath(dirname(@__FILE__), "reference/smc_cloud_fix=true.jld2")
cloud = load(file, "cloud")
split_cloud(file, 2)
rejoined_cloud = join_cloud(file, 2)

@testset "Test split and join clouds" begin
    @test cloud.particles           == rejoined_cloud.particles
    @test SMC.get_vals(cloud)       == SMC.get_vals(rejoined_cloud)
    @test SMC.get_loglh(cloud)      == SMC.get_loglh(rejoined_cloud)
    @test SMC.get_old_loglh(cloud)   == SMC.get_old_loglh(rejoined_cloud)
    @test SMC.get_logpost(cloud)    == SMC.get_logpost(rejoined_cloud)
    @test cloud.ESS                 == rejoined_cloud.ESS
    @test cloud.c                   == rejoined_cloud.c
    @test cloud.stage_index         == rejoined_cloud.stage_index
    @test cloud.total_sampling_time == rejoined_cloud.total_sampling_time
    @test cloud.accept              == rejoined_cloud.accept
    @test cloud.n_Φ                 == rejoined_cloud.n_Φ
    @test cloud.resamples           == rejoined_cloud.resamples
    @test cloud.tempering_schedule  == rejoined_cloud.tempering_schedule
end


text:
import ..@import_interface
@import_interface

#one step of gradient optimization.  This is done as a closure so that memory
#operations can be conserved.
@generated function gradientoptimize{F}(mlmodel::MLAlgorithm{F}, input::Matrix, answers::Matrix, cost::Function)
  parameter_count = parameters(mlmodel)
  output_count = outputs(mlmodel)
  quote
    #hyperparameters:
    alpha_parameter = 0.01
    delta_x = 0.0001
    #return this as the stepping function.
    #generate a storage array.
    storage_array = Vector{F}($parameter_count)
    delta_array = Vector{F}($parameter_count)
    slope_array = Vector{F}($parameter_count)

    const input_width = size(input, 2)
    result_array = Matrix{F}($output_count, input_width)

    magfactor = alpha_parameter / delta_x

    for round = 1:300
      #flatten the mlmodel into the storage array
      flatten!(storage_array, mlmodel)
      #calculate the base cost.
      batch_evaluate!(result_array, mlmodel, input, Val{input_width})

      base_cost = cost(answers, result_array)

      @inbounds delta_array[:] = storage_array[:]

      @inbounds for idx = 1:$parameter_count
        old_value = delta_array[idx]
        #nudge the delta array over a bit.
        delta_array[idx] += delta_x
        #deserialize the delta array into the model, then execute the model.
        unflatten!(mlmodel, delta_array)

        batch_evaluate!(result_array, mlmodel, input, Val{input_width})
        #calculate the cost of this execution round.
        delta_cost = cost(answers, result_array)
        #take the delta to generate the slopes.
        slope_array[idx] = base_cost - delta_cost
        #restore the old value
        delta_array[idx] = old_value
      end
      #update the storage array.
      storage_array += magfactor * slope_array
      #unflatten the new storage array into the xornet.
      unflatten!(mlmodel, storage_array)
    end
  end
end

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
using DFControl, Test, Suppressor

testdir = @__DIR__
@time begin
    test_server = Server("localhost_test", "ponet", "localhost", 8123, Servers.Bash, "", "julia", homedir(), 0)
    if Servers.isalive(test_server)
        try
            Servers.kill_server(test_server)
        catch
            nothing
        end
    end
    Servers.save(test_server)
    DFControl.Servers.start(test_server)
    if Servers.isalive(test_server)
        try
            Servers.kill_server(test_server)
        catch
            nothing
        end
    end
    @async DFControl.Resource.run(8123)
    @time @testset "constants" begin
        include("constant_tests.jl")
    end
    @time @testset "documenation" begin
        include("documentation_tests.jl")
    end
    @time @testset "Setting defaults" begin
        include("defaults_tests.jl")
    end
    @time @testset "Job from CIF file" begin
        include("jobfromcif_tests.jl")
    end
    @time @testset "Job control tests" begin
        include("job_control_tests.jl")
    end
    @time @testset "Remove defaults" begin
        include("rmdefaults_tests.jl")
    end
    rm(DFControl.config_path("servers/localhost_test.json"))
    rm(DFControl.config_path("environments/test_default.json"))
end


text:
# http://numbers.computation.free.fr/Constants/Algorithms/splitting.html

binarysplitting(a, b) = begin
    # directly compute P(a, a + 1), Q(a, a + 1) and T(a, a + 1)
    (b - a) == 1 && return compute(a, b)

    # recursively compute P(a, b), Q(a, b) and T(a, b)

    # m is the midpoint of a and b
    m = (a + b) >>> 1

    # recursively calculate P(a, m), Q(a, m) and T(a, m)
    am = binarysplitting(a, m)
    # recursively calculate P(m, b), Q(m, b) and T(m, b)
    mb = binarysplitting(m, b)

    # now combine
    combine(am, mb)
end

binarysplitting!(s, a, b) = begin
    # directly compute P(a, a + 1), Q(a, a + 1) and T(a, a + 1)
    (b - a) == 1 && (compute!(s, a, b); return)

    # recursively compute P(a, b), Q(a, b) and T(a, b)

    # m is the midpoint of a and b
    m = (a + b) >>> 1

    # recursively calculate P(a, m), Q(a, m) and T(a, m)
    binarysplitting!(s, a, m)

    # recursively calculate P(m, b), Q(m, b) and T(m, b)
    push!(s)
    binarysplitting!(s, m, b)
    pop!(s)

    # now combine
    combine!(s)
    return
end


text:
export AbstractField, Field, Fields, TextField, NumericField, GeoField, TagField, redis_args

"""
    Fields

BareModel for Fields available in the RediSearch Framework.

# Members
- `NUMERIC`
- `TEXT`
- `NOSTEM`
- `PHONETIC`
- `WEIGHT`
- `GEO`
- `TAG`
- `SEPARATOR`
- `SORTABLE`
- `NOINDEX`
- `AS`
"""
baremodule Fields
    const NUMERIC = "NUMERIC"
    const TEXT = "TEXT"
    const NOSTEM = "NOSTEM"
    const PHONETIC = "PHONETIC"
    const WEIGHT = "WEIGHT"
    const GEO = "GEO"
    const TAG = "TAG"
    const SEPARATOR = "SEPARATOR"
    const SORTABLE = "SORTABLE"
    const NOINDEX = "NOINDEX"
    const AS = "AS"
end

abstract type AbstractField end

"""
    Field::AbstractField(
        name::String
        args::Vector{Union{String,Number}}
        sortable::Bool
        no_index::Bool
        as_name::Union{String,Nothing}
    )

Field object to include in RediSearch schema.

# Arguments
- `name::String`: Name of the field.

# KeyWords
- `args::Vector{Union{String,Number}}`: List of Arguments fot the field.
- `sortable::Bool`: Sortable status.
- `no_index::Bool`: Index Status.
- `as_name::Union{String,Nothing}`: Name to reference field as.

# Returns
- `Field::AbstractField`
"""
mutable struct Field <: AbstractField
    name::String
    args::Vector{Union{String,Number}}
    sortable::Bool
    no_index::Bool
    as_name::Union{String,Nothing}
end

function Field(
    name::String;
    args=Vector{Union{String,Number}}(),
    sortable=false,
    no_index=false,
    as_name=nothing,
)
    if !isnothing(as_name)
        args = vcat(args, [Fields.AS, as_name])
    end

    if sortable
        append!(args, [Fields.SORTABLE])
    end

    if no_index
        sortable && throw(ErrorException("Cannot Allow for Sortable and No Index on Field"))
        append!(args, [Fields.NOINDEX])
    end

    return Field(name, args, sortable, no_index, as_name)
end

"""
    TextField(name; weight=1.0, no_stem=false, phonetic=nothing, kwargs...)

Creates a TextField.

# Arguments
- `name`: Name of the field.

# KeyWords
- `weight`: Weight of the field.
- `no_stem`: Stem bool.
- `phonetic`: Phonetic reference.
- `kwargs...`: Additional fields defined in the `Field` object.

# Returns
- `TextField::Field`
"""
function TextField(name; weight=1.0, no_stem=false, phonetic="", kwargs...)
    field = Field(name; args=[Fields.TEXT, Fields.WEIGHT, weight], kwargs...)

    if !isempty(phonetic) && phonetic in ["dm:en", "dm:fr", "dm:pt", "dm:es"]
        append!(field.args, [Fields.PHONETIC, phonetic])
    end

    if no_stem
        append!(field.args, [Fields.NOSTEM])
    end

    return field
end

"""
    NumericField(name; kwargs...)

Creates a NumericField object.

# Arguments
- `name`: Name of the field.

# KeyWords
- `kwargs...`: Additional fields defined in the `Field` object.

# Returns
- `NumericField::Field`
"""
function NumericField(name; kwargs...)
    return Field(name; args=[Fields.NUMERIC], kwargs...)
end

"""
    GeoField(name; kwargs...)

Creates a GeoField object.

# Arguments
- `name`: Name of the field.

# KeyWords
- `kwargs...`: Additional fields defined in the `Field` object.

# Returns
- `GeoField::Field`
"""
function GeoField(name; kwargs...)
    return Field(name; args=[Fields.GEO], kwargs...)
end

"""
    TagField(name; separator=",", kwargs...)

Creates a TagField object.

# Arguments
- `name`: Name of the field.

# KeyWords
- `separator`: Separator string.
- `kwargs...`: Additional fields defined in the `Field` object.

# Returns
- `TagField::Field`
"""
function TagField(name; separator=",", kwargs...)
    return Field(name; args=[Fields.TAG, Fields.SEPARATOR, separator], kwargs...)
end

"""
    redis_args(field::AbstractField)

Create an array of strings to represent the RedisArgs of a Field.

# Arguments
- `field::AbstractField`: Field to fetch arguements for.

# Returns
- `AbstractArray{String}`: List of strings that represent Redis arguments.
"""
function redis_args(field::AbstractField)
    return vcat(field.name, field.args)
end


text:
"""
    mean(A::AbstractArray, Estimator::EmpiricalMean)

The usual empirical mean estimator. Nothing fancy.
"""
struct EmpiricalMean <: MeanEstimator end
function mean(A::AbstractArray, Estimator::EmpiricalMean)
    return mean(A)
end

"""
MedianOfMean(x::AbstractArray, k::Integer)

Compute the Median of Mean with `k` groups (it does not permute the samples)
"""
function MoM(x::AbstractArray, k::Integer)
    n = length(x)
    @assert n % k == 0 "Check that n=$n is a mutliple of k=$k"
    m = n ÷ k
    return median(mean(@view(x[(1+(l-1)*m):(l*m)])) for l = 1:k)
end
struct MedianOfMean{T<:Integer} <: RobustMean
    k::T # number of blocks
end
"""
mean(A::AbstractArray, Estimator::MoM)

The Median of Mean estimator.
"""
function mean(A::AbstractArray, Estimator::MedianOfMean)
    return MoM(A, Estimator.k)
end

struct TrimmedMean{T<:Integer} <: RobustMean
    α::T #
    β::T
end
"""
    TrimmedMean(x::AbstractArray, k::Integer)

Compute the Trimmed Mean thresolding data smaller `α` or larger than `β`
"""
TrimMean(x::AbstractArray, α, β) = mean(φ(x[i], α, β) for i in eachindex(x))
function φ(x, α, β)
    if x < α
        return α
    elseif x > β
        return β
    else
        x
    end
end
"""
mean(A::AbstractArray, Estimator::TrimmedMean)

The Trimmed Mean estimator. (p)
"""
function mean(A::AbstractArray, Estimator::TrimmedMean)
    return TrimMean(A, Estimator.α, Estimator.β)
end

"""
    Z_estimator(x::AbstractArray, α, ψ::Function; ini = median(x))

Implement Z estimators given a
"""
function Z_estimator(x::AbstractArray, α, ψ::Function; ini=median(x))
    #! For some reason it seems that `find_zero` only accepts `Float64`
    f(y) = Rα(y, x, α, ψ)
    return find_zero(f, ini)
end

Rα(y, x::AbstractArray, α, ψ::Function) = sum(ψ(α * (x[i] - y)) for i in eachindex(x))

struct Z_Estimator{F,T} <: MeanEstimator
    α::T # scaling parameter
    ψ::F # Influence function
end
"""
mean(A::AbstractArray, Estimator::Z_Estimator)

A Z estimator given an influence function `x->ψ(x)` and a scaling parameter `α`.
"""
function mean(A::AbstractArray, Estimator::Z_Estimator; kwargs...)
    return Z_estimator(A, Estimator.α, Estimator.ψ; kwargs...)
end

"""
    ψ_Catoni(x)

Catoni's influence function
"""
ψ_Catoni(x) = x < 0 ? -log(1 - x + x^2 / 2) : log(1 + x + x^2 / 2)

"""
    ψ_Huber(x)

ψ_Huber's influence function
"""
ψ_Huber(x) = abs(x) ≤ 1 ? x : sign(x)


text:
mutable struct Route
    copyrights::String
    startaddress::String
    endaddress::String
    travelmode::Vector{String}
    traveltime::Vector{Int} # seconds
    traveldist::Vector{Int} # kilometers
    lonlat::Vector{Tuple{Float64,Float64}}
    bbox::Tuple{Float64,Float64,Float64,Float64} # minlon,minlat,maxlon,maxlat
    instructions::Vector{String}
    misc::Dict{String,Any}
    polylines::Vector{Vector{Tuple{Float64,Float64}}}

    Route() = new()
end

function directions(;kwargs...)
    resultjson = requestjson("directions"; kwargs...)
    if resultjson["status"] == "OK"
        return getroutes(resultjson["routes"])
    else
        return Route[]
    end
end

"""
Returns a vector of transit routes under the GoogleMaps Directions API.

## Parameters

* `origin`: the origin location.
* `dest`: the destination location.
* [optional] `datetime`: defaults to the current time

See https://developers.google.com/maps/documentation/directions/intro
for details on the origin and destination specification.

## Keyword Arguments

* `key::String`: Your Google Maps API Key.
* `mode`: includes `driving` (default), `walking`, `bicycling`, `transit`.
* `alternatives::Bool` (default: `true`): specifies whether the Directions
    service may provide more than one route alternative in the response.
"""
function directions(
        origin::String,
        dest::String,
        datetime::DateTime;
        key::String = "",
        mode::String = "transit",
        alternatives::Bool = true,
        kwargs...
    )
    key != "" && push!(kwargs, (:key, key))
    push!(kwargs, (:departure_time, string(round(Int,Dates.datetime2unix(datetime)))))
    push!(kwargs, (:mode, mode))
    push!(kwargs, (:alternatives, string(alternatives)))
    push!(kwargs, (:origin, origin))
    push!(kwargs, (:destination, dest))
    directions(; kwargs...)
end

directions(origin::String, dest::String; kwargs...) =
    directions(origin, dest, now(); kwargs...)

function directions(
        lonlat_origin::Tuple{Float64,Float64},
        lonlat_dest::Tuple{Float64,Float64},
        args...;
        kwargs...
    )
    directions(
        string(lonlat_origin[2], ",", lonlat_origin[1]),
        string(lonlat_dest[2], ",", lonlat_dest[1]),
        args...; kwargs...
    )
end

function Base.show(io::IO, r::Route)
    pprint(s::String, width::Int=50) =
        length(s) > width ? string(s[1:(width-3)], "...") : s
    println(io, "Route (", r.copyrights, ")")
    println(io, "=====")
    println(io, " *         from: ", pprint(r.startaddress))
    println(io, " *           to: ", pprint(r.endaddress))
    println(io, " * instructions: ", isempty(r.instructions) ? "no" : "yes")
    println(io, " *   travelmode: ", r.travelmode)
end

"""
Returns a list of routes for the result from Google Maps Directions API.

## Parameters

* `polyline::Bool` (default: `true`): whether or not to include the polyline
    corresponding to the route
* `instructions` (default: `true`): whether or not to return the
    instructions for the route.

"""
function getroutes(
        routes::Vector;
        instructions::Bool = true,
        stepcallback::Function = (route, s) -> nothing
    )
    result = Route[]
    sizehint!(result, length(routes))
    for r in routes
        route = Route()
        route.copyrights = r["copyrights"]; route.misc = Dict()

        minlon = r["bounds"]["southwest"]["lng"]
        minlat = r["bounds"]["southwest"]["lat"]
        maxlon = r["bounds"]["northeast"]["lng"]
        maxlat = r["bounds"]["northeast"]["lat"]
        route.bbox = (minlon,minlat,maxlon,maxlat)

        @assert length(r["legs"]) == 1 # no waypoints
        leg = r["legs"][1]

        route.startaddress = leg["start_address"]
        route.endaddress = leg["end_address"]

        steps = leg["steps"]
        @assert length(steps) >= 1
        route.travelmode = [s["travel_mode"] for s in steps]
        route.traveltime = [s["duration"]["value"] for s in steps]
        route.traveldist = [s["distance"]["value"] for s in steps]

        route.instructions = instructions ? [s["html_instructions"] for s in steps] : []
        route.polylines = [decodepolyline(s["polyline"]["points"]) for s in steps]

        startlon = steps[1]["start_location"]["lng"]
        startlat = steps[1]["start_location"]["lat"]
        route.lonlat = [(startlon,startlat)]
        sizehint!(route.lonlat, length(steps))
        for s in steps
            stepcallback(route, s)
            lonlat = (s["end_location"]["lng"], s["end_location"]["lat"])
            push!(route.lonlat, lonlat)
        end
        push!(result, route)
    end
    result
end


text:
include("tests.jl")
include("test-matrix.jl")
include("test-ode.jl")
include("test-logical.jl")


text:
# This file is part of TaylorSeries.jl, MIT licensed
#

using TaylorSeries

using Test
using LinearAlgebra, SparseArrays
eeuler = Base.MathConstants.e

@testset "Tests for Taylor1 expansions" begin
    ta(a) = Taylor1([a,one(a)],15)
    t = Taylor1(Int,15)
    tim = im*t
    zt = zero(t)
    ot = 1.0*one(t)
    tol1 = eps(1.0)

    @test Taylor1 <: AbstractSeries
    @test Taylor1{Float64} <: AbstractSeries{Float64}

    @test Taylor1([1,2,3,4,5], 2) == Taylor1([1,2,3])
    @test get_order(Taylor1([1,2,3,4,5], 2)) == 2

    @test size(t) == (16,)
    @test firstindex(t) == 0
    @test lastindex(t) == 15
    @test eachindex(t) == 0:15
    @test iterate(t) == (0.0, 1)
    @test iterate(t, 1) == (1.0, 2)
    @test iterate(t, 16) == nothing
    @test axes(t) == ()
    @test axes([t]) == (Base.OneTo(1),)

    v = [1,2]
    @test typeof(TaylorSeries.resize_coeffs1!(v,3)) == Nothing
    @test v == [1,2,0,0]
    TaylorSeries.resize_coeffs1!(v,0)
    @test v == [1]
    TaylorSeries.resize_coeffs1!(v,3)
    setindex!(Taylor1(v),3,2)
    @test v == [1,0,3,0]
    pol_int = Taylor1(v)
    @test pol_int[:] == [1,0,3,0]
    @test pol_int[:] == pol_int.coeffs[:]
    @test pol_int[1:2:3] == pol_int.coeffs[2:2:4]
    setindex!(pol_int,0,0:2)
    @test v == zero(v)
    setindex!(pol_int,1,:)
    @test v == ones(Int, 4)
    setindex!(pol_int, v, :)
    @test v == ones(Int, 4)
    setindex!(pol_int, zeros(Int, 4), 0:3)
    @test v == zeros(Int, 4)
    pol_int[:] .= 0
    @test v == zero(v)
    pol_int[0:2:end] = 2
    @test all(v[1:2:end] .== 2)
    pol_int[0:2:3] = [0, 1]
    @test all(v[1:2:3] .== [0, 1])
    rv = [rand(0:3) for i in 1:4]
    @test Taylor1(rv)[:] == rv
    y = sin(Taylor1(16))
    @test y[:] == y.coeffs
    y[:] .= cos(Taylor1(16))[:]
    @test y == cos(Taylor1(16))
    @test y[:] == cos(Taylor1(16))[:]
    y = sin(Taylor1(16))
    rv = rand(5)
    y[0:4] .= rv
    @test y[0:4] == rv
    @test y[5:end] == y.coeffs[6:end]
    rv = rand( length(y.coeffs) )
    y[:] .= rv
    @test y[:] == rv
    y[:] .= cos(Taylor1(16)).coeffs
    @test y == cos(Taylor1(16))
    @test y[:] == cos(Taylor1(16))[:]
    y[:] .= 0.0
    @test y[:] == zero(y[:])
    y = sin(Taylor1(16))
    rv = rand.(length(0:4))
    y[0:4] .= rv
    @test y[0:4] == rv
    @test y[6:end] == sin(Taylor1(16))[6:end]
    rv = rand.(length(y))
    y[:] .= rv
    @test y[:] == rv
    y[0:4:end] .= 1.0
    @test all(y.coeffs[1:4:end] .== 1.0)
    y[0:2:8] .= rv[1:5]
    @test y.coeffs[1:2:9] == rv[1:5]
    @test_throws AssertionError y[0:2:3] = rv

    @test Taylor1([0,1,0,0]) == Taylor1(3)
    @test getcoeff(Taylor1(Complex{Float64},3),1) == complex(1.0,0.0)
    @test Taylor1(Complex{Float64},3)[1] == complex(1.0,0.0)
    @test getindex(Taylor1(3),1) == 1.0
    @inferred convert(Taylor1{Complex{Float64}},ot) == Taylor1{Complex{Float64}}
    @test eltype(convert(Taylor1{Complex{Float64}},ot)) == Complex{Float64}
    @test eltype(convert(Taylor1{Complex{Float64}},1)) == Complex{Float64}
    @test eltype(convert(Taylor1, 1im)) == Complex{Int}
    @test convert(Taylor1, 1im) == Taylor1(1im, 0)
    @test convert(Taylor1{Int},[0,2]) == 2*t
    @test convert(Taylor1{Complex{Int}},[0,2]) == (2+0im)*t
    @test convert(Taylor1{BigFloat},[0.0, 1.0]) == ta(big(0.0))
    @test promote(t,Taylor1(1.0,0)) == (ta(0.0),ot)
    @test promote(0,Taylor1(1.0,0)) == (zt,ot)
    @test eltype(promote(ta(0),zeros(Int,2))[2]) == Int
    @test eltype(promote(ta(0.0),zeros(Int,2))[2]) == Float64
    @test eltype(promote(0,Taylor1(ot))[1]) == Float64
    @test eltype(promote(1.0+im, zt)[1]) == Complex{Float64}

    @test length(Taylor1(10)) == 11
    @test length.( TaylorSeries.fixorder(zt, Taylor1([1])) ) == (16, 16)
    @test length.( TaylorSeries.fixorder(zt, Taylor1([1], 1)) ) == (2, 2)
    @test eltype(TaylorSeries.fixorder(zt,Taylor1([1]))[1]) == Int
    @test TaylorSeries.findfirst(t) == 1
    @test TaylorSeries.findfirst(t^2) == 2
    @test TaylorSeries.findfirst(ot) == 0
    @test TaylorSeries.findfirst(zt) == -1
    @test TaylorSeries.findlast(t) == 1
    @test TaylorSeries.findlast(t^2) == 2
    @test TaylorSeries.findlast(ot) == 0
    @test TaylorSeries.findlast(zt) == -1
    @test iszero(zero(t))
    @test !iszero(one(t))
    @test isinf(Taylor1([typemax(1.0)]))
    @test isnan(Taylor1([typemax(1.0), NaN]))

    @test constant_term(2.0) == 2.0
    @test constant_term(t) == 0
    @test constant_term(tim) == complex(0, 0)
    @test constant_term([zt, t]) == [0, 0]
    @test linear_polynomial(2) == 2
    @test linear_polynomial(t) == t
    @test linear_polynomial(tim^2) == zero(tim)
    @test linear_polynomial([zero(tim), tim]) == [zero(tim), tim]

    @test ot == 1
    @test 0.0 == zt
    @test getcoeff(tim,1) == complex(0,1)
    @test zt+1.0 == ot
    @test 1.0-ot == zt
    @test t+t == 2t
    @test t-t == zt
    @test +t == -(-t)

    tsquare = Taylor1([0,0,1],15)
    @test t * true == t
    @test false * t == zero(t)
    @test t^0 == t^0.0 == one(t)
    @test t*t == tsquare
    @test t*1 == t
    @test 0*t == zt
    @test (-t)^2 == tsquare
    @test t^3 == tsquare*t
    @test zero(t)/t == zero(t)
    @test one(t)/one(t) == 1.0
    @test tsquare/t == t
    @test t/(t*3) == (1/3)*ot
    @test t/3im == -tim/3
    @test 1/(1-t) == Taylor1(ones(t.order+1))
    @test Taylor1([0,1,1])/t == t+1
    @test (t+im)^2 == tsquare+2im*t-1
    @test (t+im)^3 == Taylor1([-1im,-3,3im,1],15)
    @test (t+im)^4 == Taylor1([1,-4im,-6,4im,1],15)
    @test imag(tsquare+2im*t-1) == 2t
    @test (Rational(1,2)*tsquare)[2] == 1//2
    @test t^2/tsquare == ot
    @test ((1+t)^(1/3))[2]+1/9 ≤ tol1
    @test (1.0-tsquare)^3 == (1.0-t)^3*(1.0+t)^3
    @test (1-tsquare)^2 == (1+t)^2.0 * (1-t)^2.0
    @test (sqrt(1+t))[2] == -1/8
    @test ((1-tsquare)^(1//2))^2 == 1-tsquare
    @test ((1-t)^(1//4))[14] == -4188908511//549755813888
    @test abs(((1+t)^3.2)[13] + 5.4021062656e-5) < tol1
    @test Taylor1(BigFloat,5)/6 == 1im*Taylor1(5)/complex(0,BigInt(6))
    @test Taylor1(BigFloat,5)/(6*Taylor1(3)) == 1/BigInt(6)
    @test Taylor1(BigFloat,5)/(6im*Taylor1(3)) == -1im/BigInt(6)

    # These tests involve some sort of factorization
    @test t/(t+t^2) == 1/(1+t)
    @test sqrt(t^2+t^3) == t*sqrt(1+t)
    @test (t^3+t^4)^(1/3) ≈ t*(1+t)^(1/3)
    @test norm((t^3+t^4)^(1/3) - t*(1+t)^(1/3), Inf) < eps()
    @test ((t^3+t^4)^(1/3))[15] ≈ -8617640/1162261467

    trational = ta(0//1)
    @inferred ta(0//1) == Taylor1{Rational{Int}}
    @test eltype(trational) == Rational{Int}
    @test trational + 1//3 == Taylor1([1//3,1],15)
    @test complex(3,1)*trational^2 == Taylor1([0//1,0//1,complex(3,1)//1],15)
    @test trational^2/3 == Taylor1([0//1,0//1,1//3],15)
    @test trational^3/complex(7,1) == Taylor1([0,0,0,complex(7//50,-1//50)],15)
    @test sqrt(zero(t)) == zero(t)

    @test isapprox( rem(4.1 + t,4)[0], 0.1 )
    @test isapprox( mod(4.1 + t,4)[0], 0.1 )
    @test isapprox( rem(1+Taylor1(Int,4),4.0)[0], 1.0 )
    @test isapprox( mod(1+Taylor1(Int,4),4.0)[0], 1.0 )
    @test isapprox( mod2pi(2pi+0.1+t)[0], 0.1 )

    @test abs(ta(1)) == ta(1)
    @test abs(ta(-1.0)) == -ta(-1.0)

    @test taylor_expand(x->2x,order=10) == 2*Taylor1(10)
    @test taylor_expand(x->x^2+1) == Taylor1(15)*Taylor1(15) + 1
    @test evaluate(taylor_expand(cos,0.)) == cos(0.)
    @test evaluate(taylor_expand(tan,pi/4)) == tan(pi/4)
    @test eltype(taylor_expand(x->x^2+1,1)) == eltype(1)
    tsq = t^2
    update!(tsq,2.0)
    @test tsq == (t+2.0)^2
    update!(tsq,-2.0)
    @test tsq == t^2

    @test log(exp(tsquare)) == tsquare
    @test exp(log(1-tsquare)) == 1-tsquare
    @test log((1-t)^2) == 2*log(1-t)
    @test real(exp(tim)) == cos(t)
    @test imag(exp(tim)) == sin(t)
    @test exp(conj(tim)) == cos(t)-im*sin(t) == exp(tim')
    @test (exp(t))^(2im) == cos(2t)+im*sin(2t)
    @test (exp(t))^Taylor1([-5.2im]) == cos(5.2t)-im*sin(5.2t)
    @test getcoeff(convert(Taylor1{Rational{Int}},cos(t)),8) == 1//factorial(8)
    @test abs((tan(t))[7]- 17/315) < tol1
    @test abs((tan(t))[13]- 21844/6081075) < tol1
    @test evaluate(exp(Taylor1([0,1],17)),1.0) == 1.0*eeuler
    @test evaluate(exp(Taylor1([0,1],1))) == 1.0
    @test evaluate(exp(t),t^2) == exp(t^2)
    #Test function-like behavior for Taylor1s
    t17 = Taylor1([0,1],17)
    myexpfun = exp(t17)
    @test myexpfun(1.0) == 1.0*eeuler
    @test myexpfun() == 1.0
    @test myexpfun(t17^2) == exp(t17^2)
    @test exp(t17^2)(t17) == exp(t17^2)
    p = cos(t17)
    q = sin(t17)
    @test cos(-im*t)(1)+im*sin(-im*t)(1) == exp(-im*t)(im)
    @test p(-im*t17)(1)+im*q(-im*t17)(1) == exp(-im*t17)(im)
    cossin1 = x->p(q(x))
    @test evaluate(p, evaluate(q, pi/4)) == cossin1(pi/4)
    cossin2 = p(q)
    @test evaluate(evaluate(p,q), pi/4) == cossin2(pi/4)
    @test evaluate(p, q) == cossin2
    @test p(q)() == evaluate(evaluate(p, q))
    @test evaluate(p, q) == p(q)
    @test evaluate(q, p) == q(p)
    cs = x->cos(sin(x))
    csdiff = (cs(t17)-cossin2(t17)).(-2:0.1:2)
    @test norm(csdiff, 1) < 5e-15
    a = [p,q]
    @test a(0.1) == evaluate.([p,q],0.1)
    @test a.(0.1) == a(0.1)
    @test a.() == evaluate.([p, q])
    @test a.() == [p(), q()]
    @test a.() == a()
    vr = rand(2)
    @test p.(vr) == evaluate.([p], vr)
    Mr = rand(3,3,3)
    @test p.(Mr) == evaluate.([p], Mr)
    mytaylor1 = Taylor1(rand(20))
    vr = rand(5)
    @test p(vr) == p.(vr)
    @test p(vr) == evaluate.([p],vr)
    @test p(Mr) == p.(Mr)
    @test p(Mr) == evaluate.([p], Mr)
    taylor_a = Taylor1(Int,10)
    taylor_x = exp(Taylor1(Float64,13))
    @test taylor_x(taylor_a) == evaluate(taylor_x, taylor_a)
    A_T1 = [t 2t 3t; 4t 5t 6t ]
    @test evaluate(A_T1,1.0) == [1.0  2.0  3.0; 4.0  5.0  6.0]
    @test evaluate(A_T1,1.0) == A_T1(1.0)
    @test evaluate(A_T1) == A_T1()
    @test A_T1(tsquare) == [tsquare 2tsquare 3tsquare; 4tsquare 5tsquare 6tsquare]
    @test view(A_T1, :, :)(1.0) == A_T1(1.0)
    @test view(A_T1, :, 1)(1.0) == A_T1[:,1](1.0)

    @test sin(asin(tsquare)) == tsquare
    @test tan(atan(tsquare)) == tsquare
    @test atan(tan(tsquare)) == tsquare
    @test asin(t) + acos(t) == pi/2
    @test derivative(acos(t)) == - 1/sqrt(1-t^2)

    @test - sinh(t) + cosh(t) == exp(-t)
    @test  sinh(t) + cosh(t) == exp(t)
    @test evaluate(- sinh(t)^2 + cosh(t)^2 , rand()) == 1
    @test evaluate(- sinh(t)^2 + cosh(t)^2 , 0) == 1
    @test tanh(t + 0im) == -1im * tan(t*1im)
    @test  evaluate(tanh(t/2),1.5) == evaluate(sinh(t) / (cosh(t) + 1),1.5)
    @test cosh(t) == real(cos(im*t))
    @test sinh(t) == imag(sin(im*t))

    ut = 1.0*t
    tt = zero(ut)
    TaylorSeries.one!(tt, ut, 0)
    @test tt[0] == 1.0
    TaylorSeries.one!(tt, ut, 1)
    @test tt[1] == 0.0
    TaylorSeries.abs!(tt, 1.0+ut, 0)
    @test tt[0] == 1.0
    TaylorSeries.add!(tt, ut, ut, 1)
    @test tt[1] == 2.0
    TaylorSeries.add!(tt, -3.0, 0)
    @test tt[0] == -3.0
    TaylorSeries.add!(tt, -3.0, 1)
    @test tt[1] == 0.0
    TaylorSeries.subst!(tt, ut, ut, 1)
    @test tt[1] == 0.0
    TaylorSeries.subst!(tt, -3.0, 0)
    @test tt[0] == 3.0
    TaylorSeries.subst!(tt, -2.5, 1)
    @test tt[1] == 0.0
    iind, cind = TaylorSeries.divfactorization(ut, ut)
    @test iind == 1
    @test cind == 1.0
    TaylorSeries.div!(tt, ut, ut, 0, iind)
    @test tt[0] == cind
    TaylorSeries.div!(tt, 1+ut, 1+ut, 0)
    @test tt[0] == 1.0
    TaylorSeries.div!(tt, 1, 1+ut, 0)
    @test tt[0] == 1.0
    TaylorSeries.pow!(tt, 1.0+t, 1.5, 0, 0)
    @test tt[0] == 1.0
    TaylorSeries.pow!(tt, 1.0+t, 1.5, 0)
    @test tt[0] == 1.0
    TaylorSeries.pow!(tt, 1.0+t, 0.5, 1)
    @test tt[1] == 0.5
    TaylorSeries.pow!(tt, 1.0+t, 0, 0)
    @test tt[0] == 1.0
    TaylorSeries.pow!(tt, 1.0+t, 1, 1)
    @test tt[1] == 1.0
    tt = zero(ut)
    TaylorSeries.pow!(tt, 1.0+t, 2, 0)
    @test tt[0] == 1.0
    TaylorSeries.pow!(tt, 1.0+t, 2, 1)
    @test tt[1] == 2.0
    TaylorSeries.pow!(tt, 1.0+t, 2, 2)
    @test tt[2] == 1.0
    TaylorSeries.sqrt!(tt, 1.0+t, 0, 0)
    @test tt[0] == 1.0
    TaylorSeries.sqrt!(tt, 1.0+t, 0)
    @test tt[0] == 1.0
    TaylorSeries.exp!(tt, 1.0*t, 0)
    @test tt[0] == exp(t[0])
    TaylorSeries.log!(tt, 1.0+t, 0)
    @test tt[0] == 0.0
    ct = zero(ut)
    TaylorSeries.sincos!(tt, ct, 1.0*t, 0)
    @test tt[0] == sin(t[0])
    @test ct[0] == cos(t[0])
    TaylorSeries.tan!(tt, 1.0*t, ct, 0)
    @test tt[0] == tan(t[0])
    @test ct[0] == tan(t[0])^2
    TaylorSeries.asin!(tt, 1.0*t, ct, 0)
    @test tt[0] == asin(t[0])
    @test ct[0] == sqrt(1.0-t[0]^2)
    TaylorSeries.acos!(tt, 1.0*t, ct, 0)
    @test tt[0] == acos(t[0])
    @test ct[0] == sqrt(1.0-t[0]^2)
    TaylorSeries.atan!(tt, ut, ct, 0)
    @test tt[0] == atan(t[0])
    @test ct[0] == 1.0+t[0]^2
    TaylorSeries.sinhcosh!(tt, ct, ut, 0)
    @test tt[0] == sinh(t[0])
    @test ct[0] == cosh(t[0])
    TaylorSeries.tanh!(tt, ut, ct, 0)
    @test tt[0] == tanh(t[0])
    @test ct[0] == tanh(t[0])^2

    v = [sin(t), exp(-t)]
    vv = Vector{Float64}(undef, 2)
    @test evaluate!(v, zero(Int), vv) == nothing
    @test vv == [0.0,1.0]
    @test evaluate(v) == vv
    @test evaluate(v, complex(0.0,0.2)) ==
        [complex(0.0,sinh(0.2)),complex(cos(0.2),sin(-0.2))]

    @test differentiate(exp(ta(1.0)), 0) == exp(ta(1.0))
    expected_result_approx = Taylor1(convert(Vector{Float64},exp(ta(1.0))[0:10]))
    @test derivative(exp(ta(1.0)), 5) ≈ expected_result_approx atol=eps() rtol=0.0
    expected_result_approx = Taylor1(convert(Vector{Float64},exp(ta(1.0pi))[0:12]),15)
    @test derivative(exp(ta(1.0pi)), 3) ≈ expected_result_approx atol=eps(16.0) rtol=0.0
    expected_result_approx = Taylor1(convert(Vector{Float64},exp(ta(1.0pi))[0:5]),15)
    @test derivative(exp(ta(1.0pi)), 10) ≈ expected_result_approx atol=eps(64.0) rtol=0.0

    @test derivative(exp(ta(1.0)), 5)() == exp(1.0)
    @test derivative(exp(ta(1.0pi)), 3)() == exp(1.0pi)
    @test isapprox(derivative(exp(ta(1.0pi)), 10)() , exp(1.0pi) )

    @test derivative(5, exp(ta(1.0))) == exp(1.0)
    @test derivative(3, exp(ta(1.0pi))) == exp(1.0pi)
    @test isapprox(derivative(10, exp(ta(1.0pi))) , exp(1.0pi) )

    @test integrate(derivative(exp(t)),1) == exp(t)
    @test integrate(cos(t)) == sin(t)

    @test promote(ta(0.0), t) == (ta(0.0),ta(0.0))

    @test norm((inverse(exp(t)-1) - log(1+t)).coeffs) < 2tol1
    cfs = [(-n)^(n-1)/factorial(n) for n = 1:15]
    @test norm(inverse(t*exp(t))[1:end]./cfs .- 1) < 4tol1

    @test_throws ArgumentError Taylor1([1,2,3], -2)
    @test_throws ArgumentError abs(ta(big(0)))
    @test_throws ArgumentError 1/t
    @test_throws ArgumentError zt/zt
    @test_throws ArgumentError t^1.5
    @test_throws ArgumentError t^(-2)
    @test_throws ArgumentError sqrt(t)
    @test_throws ArgumentError log(t)
    @test_throws ArgumentError cos(t)/sin(t)
    @test_throws AssertionError derivative(30, exp(ta(1.0pi)))
    @test_throws ArgumentError inverse(exp(t))
    @test_throws ArgumentError abs(t)

    use_show_default(true)
    aa = sqrt(2)+Taylor1(2)
    @test string(aa) == "Taylor1{Float64}([1.4142135623730951, 1.0, 0.0], 2)"
    @test string([aa, aa]) ==
        "Taylor1{Float64}[Taylor1{Float64}([1.4142135623730951, 1.0, 0.0], 2), " *
        "Taylor1{Float64}([1.4142135623730951, 1.0, 0.0], 2)]"
    use_show_default(false)
    @test string(aa) == " 1.4142135623730951 + 1.0 t + 𝒪(t³)"
    set_taylor1_varname(" x ")
    @test string(aa) == " 1.4142135623730951 + 1.0 x + 𝒪(x³)"
    set_taylor1_varname("t")
    displayBigO(false)
    @test string(ta(-3)) == " - 3 + 1 t "
    @test string(ta(0)^3-3) == " - 3 + 1 t³ "
    @test TaylorSeries.pretty_print(ta(3im)) == " ( 3 im )  + ( 1 ) t "
    @test string(Taylor1([1,2,3,4,5], 2)) == string(Taylor1([1,2,3]))
    displayBigO(true)
    @test string(ta(-3)) == " - 3 + 1 t + 𝒪(t¹⁶)"
    @test string(ta(0)^3-3) == " - 3 + 1 t³ + 𝒪(t¹⁶)"
    @test TaylorSeries.pretty_print(ta(3im)) == " ( 3 im )  + ( 1 ) t + 𝒪(t¹⁶)"
    @test string(Taylor1([1,2,3,4,5], 2)) == string(Taylor1([1,2,3]))

    a = collect(1:12)
    t_a = Taylor1(a,15)
    t_C = complex(3.0,4.0) * t_a
    rnd = rand(10)
    @test typeof( norm(Taylor1(rnd)) ) == Float64
    @test norm(Taylor1(rnd)) > 0
    @test norm(t_a) == norm(a)
    @test norm(Taylor1(a,15),3) == sum((a.^3))^(1/3)
    @test norm(t_a,Inf) == 12
    @test norm(t_C) == norm(complex(3.0,4.0)*a)

    @test TaylorSeries.rtoldefault(Taylor1{Int}) == 0
    @test TaylorSeries.rtoldefault(Taylor1{Float64}) == sqrt(eps(Float64))
    @test TaylorSeries.rtoldefault(Taylor1{BigFloat}) == sqrt(eps(BigFloat))
    @test TaylorSeries.real(Taylor1{Float64}) == Taylor1{Float64}
    @test TaylorSeries.real(Taylor1{Complex{Float64}}) == Taylor1{Float64}
    @test isfinite(t_C)
    @test isfinite(t_a)
    @test !isfinite( Taylor1([0, Inf]) )
    @test !isfinite( Taylor1([NaN, 0]) )
    b = convert(Vector{Float64}, a)
    b[3] += eps(10.0)
    b[5] -= eps(10.0)
    t_b = Taylor1(b,15)
    t_C2 = t_C+eps(100.0)
    t_C3 = t_C+eps(100.0)*im
    @test isapprox(t_C, t_C)
    @test t_a ≈ t_a
    @test t_a ≈ t_b
    @test t_C ≈ t_C2
    @test t_C ≈ t_C3
    @test t_C3 ≈ t_C2
    t = Taylor1(25)
    p = sin(t)
    q = sin(t+eps())
    @test t ≈ t
    @test t ≈ t+sqrt(eps())
    @test isapprox(p, q, atol=eps())

    t = Taylor1(35)
    @test Taylor1([180.0, rad2deg(1.0)], 35) == rad2deg(pi+t)
    @test sin(pi/2+deg2rad(1.0)t) == sin(deg2rad(90+t))
    a = Taylor1(rand(10))
    b = Taylor1(rand(10))
    c = deepcopy(a)
    TaylorSeries.deg2rad!(b, a, 0)
    @test a == c
    @test a[0]*(pi/180) == b[0]
    # TaylorSeries.deg2rad!.(b, a, [0,1,2])
    # @test a == c
    # for i in 0:2
    #     @test a[i]*(pi/180) == b[i]
    # end
    a = Taylor1(rand(10))
    b = Taylor1(rand(10))
    c = deepcopy(a)
    TaylorSeries.rad2deg!(b, a, 0)
    @test a == c
    @test a[0]*(180/pi) == b[0]
    # TaylorSeries.rad2deg!.(b, a, [0,1,2])
    # @test a == c
    # for i in 0:2
    #     @test a[i]*(180/pi) == b[i]
    # end

    # Test additional Taylor1 constructors
    @test Taylor1{Float64}(true) == Taylor1([1.0])
    @test Taylor1{Float64}(false) == Taylor1([0.0])
    @test Taylor1{Int}(true) == Taylor1([1])
    @test Taylor1{Int}(false) == Taylor1([0])
end

@testset "Test `inv` for `Matrix{Taylor1{Float64}}``" begin
    t = Taylor1(5)
    a = Diagonal(rand(0:10,3)) + rand(3, 3)
    ainv = inv(a)
    b = Taylor1.(a, 5)
    binv = inv(b)
    tol = 1.0e-11

    for its = 1:10
        a .= Diagonal(rand(2:12,3)) + rand(3, 3)
        ainv .= inv(a)
        b .= Taylor1.(a, 5)
        binv .= inv(b)
        @test norm(binv - ainv, Inf) ≤ tol
        @test norm(b*binv - I, Inf) ≤ tol
        @test norm(binv*b - I, Inf) ≤ tol
        @test norm(triu(b)*inv(UpperTriangular(b)) - I, Inf) ≤ tol
        @test norm(inv(LowerTriangular(b))*tril(b) - I, Inf) ≤ tol

        b .= b .+ t
        binv .= inv(b)
        @test norm(b*binv - I, Inf) ≤ tol
        @test norm(binv*b - I, Inf) ≤ tol
        @test norm(triu(b)*inv(triu(b)) - I, Inf) ≤ tol
        @test norm(inv(tril(b))*tril(b) - I, Inf) ≤ tol
    end
end

@testset "Matrix multiplication for Taylor1" begin
    order = 30
    n1 = 100
    k1 = 90

    order = max(n1,k1)
    B1 = randn(n1,order)
    Y1 = randn(k1,order)

    A1  = randn(k1,n1)

    for A in (A1,sparse(A1))
        # B and Y contain elements of different orders
        B  = Taylor1{Float64}[Taylor1(collect(B1[i,1:i]),i) for i=1:n1]
        Y  = Taylor1{Float64}[Taylor1(collect(Y1[k,1:k]),k) for k=1:k1]
        Bcopy = deepcopy(B)
        mul!(Y,A,B)

        # do we get the same result when using the `A*B` form?
        @test A*B≈Y
        # Y should be extended after the multilpication
        @test reduce(&, [y1.order for y1 in Y] .== Y[1].order)
        # B should be unchanged
        @test B==Bcopy

        # is the result compatible with the matrix multiplication?  We
        # only check the zeroth order of the Taylor series.
        y1=sum(Y)[0]
        Y=A*B1[:,1]
        y2=sum(Y)

        # There is a small numerical error when comparing the generic
        # multiplication and the specialized version
        @test abs(y1-y2) < n1*(eps(y1)+eps(y2))

        @test_throws DimensionMismatch mul!(Y,A[:,1:end-1],B)
        @test_throws DimensionMismatch mul!(Y,A[1:end-1,:],B)
        @test_throws DimensionMismatch mul!(Y,A,B[1:end-1])
        @test_throws DimensionMismatch mul!(Y[1:end-1],A,B)
    end
end


text:
#####
##### Loss Plot
#####

function plot_losses(getlosses, range, title, xlabel=:none)
  fields = fieldnames(Report.Loss)
  labels = [string(f) for _ in 1:1, f in fields]
  data = [(getlosses(i)..., f) for i in range, f in fields]
  xs = map(d -> d[1], data)
  ys = map(d -> getfield(d[2], d[3]), data)
  xlims = (minimum(xs), maximum(xs))
  return Plots.plot(xs, ys,
    label=labels, title=title, xlims=xlims, ylims=(0, Inf), xlabel=xlabel)
end

#####
##### Iteration summary plots
#####

function learning_iter_plot(rep::Report.Learning, params::Params)
  n = length(rep.checkpoints)
  nbatches = rep.checkpoints[end].batch_id
  losses_plot = plot_losses(0:n, "Losses") do i
    if i == 0
      (0, rep.initial_status.loss)
    else
      (rep.checkpoints[i].batch_id, rep.checkpoints[i].status_after.loss)
    end
  end
  checkpoints_plot = Plots.hline(
    [0, params.arena.update_threshold],
    title="Checkpoints")
  Plots.plot!(checkpoints_plot,
    [c.batch_id for c in rep.checkpoints],
    [c.reward for c in rep.checkpoints],
    ylims=(-1.0, 1.0),
    t=:scatter,
    legend=:none)
  Plots.xlims!(losses_plot, (0, nbatches))
  Plots.xlims!(checkpoints_plot, (0, nbatches))
  return Plots.plot(losses_plot, checkpoints_plot, layout=(2, 1))
end

function performances_plot(rep::Report.Iteration)
  # Global
  global_labels = []
  global_content = []
  push!(global_labels, "Self Play")
  push!(global_content, rep.perfs_self_play.time)
  if !isnothing(rep.memory)
    push!(global_labels, "Memory Analysis")
    push!(global_content, rep.perfs_memory_analysis.time)
  end
  push!(global_labels, "Learning")
  push!(global_content, rep.perfs_learning.time)
  glob = Plots.pie(global_labels, global_content,
    title="Global",
    legend=:right)
  # Self-play details
  self_play =
    let gcratio =
      rep.perfs_self_play.gc_time / rep.perfs_self_play.time
    let itratio = rep.self_play.inference_time_ratio
      Plots.pie(
        ["MCTS (inference)", "MCTS (other)", "GC"],
        [(1 - gcratio) * itratio, (1 - gcratio) * (1 - itratio), gcratio],
        title="Self Play") end end
  # Learning details
  learning = Plots.pie(
    ["Samples conversion", "Loss computation", "Optimization", "Arena (MCTS)"],
    [ rep.learning.time_convert,
      rep.learning.time_loss,
      rep.learning.time_train,
      rep.learning.time_eval],
    title="Learning")
  return Plots.plot(glob, self_play, learning)
end

function plot_iteration(
    report::Report.Iteration,
    params::Params,
    dir::String,
    itc::Int)
  # Summary plot
  splot = learning_iter_plot(report.learning, params)
  # Performances plot
  pplot = performances_plot(report)
  # Losses plot
  losses = Util.momentum_smoothing(report.learning.losses, 0.1)
  lplot = Plots.plot(collect(eachindex(losses)), losses,
    title="Loss on Minibatches",
    ylims=(0, Inf),
    legend=nothing,
    xlabel="Batch number")
  # Saving everything
  plots = [splot, pplot, lplot]
  names = ["iter_summary", "iter_perfs", "iter_loss"]
  for (plot, name) in zip(plots, names)
    pdir = joinpath(dir, name)
    isdir(pdir) || mkdir(pdir)
    Plots.savefig(plot, joinpath(pdir, "$itc"))
  end
end

#####
##### Training summary plots
#####

function plot_benchmark(
    params::Params,
    benchs::Vector{Benchmark.Report},
    dir::String)
  isempty(benchs) && return
  n = length(benchs) - 1
  nduels = length(benchs[1])
  nduels >= 1 || return
  @assert all(length(b) == nduels for b in benchs)
  isdir(dir) || mkpath(dir)
  labels = ["$(d.player) / $(d.baseline)" for _ in 1:1, d in benchs[1]]
  # Average reward
  avgz_data = [[b[i].avgz for b in benchs] for i in 1:nduels]
  avgz = Plots.plot(0:n,
    avgz_data,
    title="Average Reward",
    ylims=(-1.0, 1.0),
    legend=:bottomright,
    label=labels,
    xlabel="Iteration number")
  Plots.savefig(avgz, joinpath(dir, "benchmark_reward"))
  if params.ternary_rewards
    function compute_percentage(b, f)
      stats = Benchmark.TernaryOutcomeStatistics(b)
      return ceil(Int, 100 * (f(stats) / length(b.rewards)))
    end
    function compute_data(f)
      return [[compute_percentage(b[i], f) for b in benchs] for i in 1:nduels]
    end
    # Percentage of lost games
    ploss = Plots.plot(0:n,
      compute_data(s -> s.num_lost),
      title="Percentage of Lost Games",
      ylims=(0, 100),
      legend=:topright,
      label=labels,
      xlabel="Iteration number")
    Plots.savefig(ploss, joinpath(dir, "benchmark_lost_games"))
    # Percentage of won games
    pwin = Plots.plot(0:n,
      compute_data(s -> s.num_won),
      title="Percentage of Won Games",
      ylims=(0, 100),
      legend=:bottomright,
      label=labels,
      xlabel="Iteration number")
    Plots.savefig(pwin, joinpath(dir, "benchmark_won_games"))
  end
end

function plot_training(
    params::Params,
    iterations::Vector{Report.Iteration},
    dir::String)
  n = length(iterations)
  iszero(n) && return
  isdir(dir) || mkpath(dir)
  plots, files = [], []
  # Exploration depth
  expdepth = Plots.plot(1:n,
    [it.self_play.average_exploration_depth for it in iterations],
    title="Average Exploration Depth",
    ylims=(0, Inf),
    legend=:none,
    xlabel="Iteration number")
  # Number of samples
  nsamples = Plots.plot(0:n,
    [0;[it.self_play.memory_size for it in iterations]],
    title="Experience Buffer Size",
    label="Number of samples",
    xlabel="Iteration number")
  Plots.plot!(nsamples, 0:n,
    [0;[it.self_play.memory_num_distinct_boards for it in iterations]],
    label="Number of distinct boards")
  # Performances during evaluation
  arena = Plots.plot(1:n, [
    maximum(c.reward for c in it.learning.checkpoints)
    for it in iterations],
    title="Arena Results",
    ylims=(-1, 1),
    t=:bar,
    legend=:none,
    xlabel="Iteration number")
  Plots.hline!(arena, [0, params.arena.update_threshold])
  # Loss on the full memory after an iteration
  lfmt = "Loss on Full Memory"
  losses_fullmem = plot_losses(1:n, lfmt, "Iteration number") do i
    (i, iterations[i].learning.initial_status.loss)
  end
  # Plots related to the memory analysis
  if all(it -> !isnothing(it.memory), iterations)
    # Loss on last batch
    losses_last = plot_losses(
      1:n, "Loss on Last Batch", "Iteration number") do i
      (i, iterations[i].memory.latest_batch.status.loss)
    end
    # Loss per game stage
    nstages = minimum(length(it.memory.per_game_stage) for it in iterations)
    colors = range(colorant"blue", stop=colorant"red", length=nstages)
    losses_ps = Plots.plot(
      title="Loss per Game Stage", ylims=(0, Inf), xlabel="Iteration number")
    for s in 1:nstages
      tmin = minimum([
        it.memory.per_game_stage[s].min_remaining_length
        for it in iterations])
      tmax = maximum([
        it.memory.per_game_stage[s].max_remaining_length
        for it in iterations])
      Plots.plot!(losses_ps, 1:n, [
          it.memory.per_game_stage[s].samples_stats.status.loss.L
          for it in iterations],
        label="$tmin to $tmax turns left",
        color=colors[s])
    end
    append!(plots, [losses_last, losses_ps])
    append!(files, ["loss_last_batch", "loss_per_stage"])
  end
  # Policies entropy
  entropies = Plots.plot(1:n,
    [it.learning.initial_status.Hp for it in iterations],
    ylims=(0, Inf),
    title="Policy Entropy",
    label="MCTS",
    xlabel="Iteration number")
  Plots.plot!(entropies, 1:n,
    [it.learning.initial_status.Hpnet for it in iterations],
    label="Network")
  # Assembling everything together
  append!(plots,
    [arena, entropies, nsamples, expdepth, losses_fullmem])
  append!(files,
    ["arena", "entropies", "nsamples", "exploration_depth", "loss"])
  for (file, plot) in zip(files, plots)
    Plots.savefig(plot, joinpath(dir, file))
  end
end
# To test:
# AlphaZero.plot_training("sessions/connect-four")

max_stars_repo_name: btheodorou/Bridge-NN-Research

max_issues_repo_name: btheodorou/Bridge-NN-Research

max_forks_repo_name: btheodorou/Bridge-NN-Research


text:
using DataFrames, FixedEffectModels, Random
N = 10000000
K = 100
id1 = rand(1:div(N, K), N)
id2 = rand(1:K, N)
w = cos.(id1)

x1 = 5 * cos.(id1) + 5 * sin.(id2) + randn(N)
x2 =  cos.(id1) +  sin.(id2) + randn(N)
y= 3 .* x1 .+ 5 .* x2 .+ cos.(id1) .+ cos.(id2).^2 .+ randn(N)

df = DataFrame(id1 = categorical(id1), id2 = categorical(id2), x1 = x1, x2 = x2, w = w, y = y)
df.x3 =  cos.(id1) + sin.(id2) + randn(N)
df.x4 =  cos.(id1) + sin.(id2) + randn(N)
df.x5 =  cos.(id1) + sin.(id2) + randn(N)
df.x6 =  cos.(id1) + sin.(id2) + randn(N)
df.x7 =  cos.(id1) + sin.(id2) + randn(N)

@time reg(df, @model(y ~ x1 + x2))
# 0.582029 seconds (852 allocations: 535.311 MiB, 18.28% gc time)
@time reg(df, @model(y ~ x1 + x2, vcov = cluster(id2)))
# 0.809652 seconds (1.55 k allocations: 649.821 MiB, 14.40% gc time)
@time reg(df, @model(y ~ x1 + x2 + fe(id1)))
# 1.655732 seconds (1.21 k allocations: 734.353 MiB, 16.88% gc time)
@time reg(df, @model(y ~ x1 + x2 + fe(id1), vcov = cluster(id1)))
#  2.113248 seconds (499.66 k allocations: 811.751 MiB, 15.08% gc time)
@time reg(df, @model(y ~ x1 + x2 + fe(id1) + fe(id2)))
# 3.553678 seconds (1.36 k allocations: 1005.101 MiB, 10.55% gc time))
@time reg(df, @model(y ~ x1 + x2 + fe(id1) + fe(id2), vcov = cluster(id1 + id2)))

# more regressors
@time reg(df, @model(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + fe(id1))
#  5.654325 seconds (3.42 k allocations: 3.993 GiB, 18.19% gc time)
@time reg(df, @model(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + fe(id1) + fe(id2))
# 10.520039 seconds (3.50 k allocations: 4.334 GiB, 6.19% gc time)
@time reg(df, @model(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + fe(id1)+ fe(id2)), method = :lsmr_threads)
#  9.591568 seconds (3.94 k allocations: 6.372 GiB, 5.24% gc time)