text:
module PlotDocs

using Plots, DataFrames, MacroTools, OrderedCollections, Dates
import Plots: _examples

export
    generate_markdown,
    generate_supported_markdown,
    generate_attr_markdown,
    generate_graph_attr_markdown,
    generate_colorschemes_markdown,
    GENDIR

const GENDIR = normpath(@__DIR__, "..", "docs", "src", "generated")
mkpath(GENDIR)

const PLOT_DOCS_URL = "https://github.com/JuliaPlots/PlotDocs.jl/blob/master/src/PlotDocs.jl"

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

recursive_rmlines(x) = x
function recursive_rmlines(x::Expr)
    x = rmlines(x)
    x.args .= recursive_rmlines.(x.args)
    return x
end

function pretty_print_expr(io::IO, expr::Expr)
    if expr.head == :block
        for arg in recursive_rmlines(expr).args
            println(io, arg)
        end
    else
        println(io, recursive_rmlines(expr))
    end
end

function markdown_code_to_string(arr, prefix = "")
    string("`", prefix, join(sort(map(string, arr)), "`, `$prefix"), "`")
end
markdown_symbols_to_string(arr) = isempty(arr) ? "" : markdown_code_to_string(arr, ":")

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

function generate_markdown(pkgname::Symbol; skip = get(Plots._backend_skips, pkgname, Int[]))
    pkg = Plots._backend_instance(pkgname)

    # open the markdown file
    md = open(joinpath(GENDIR, "$(pkgname).md"), "w")

    write(md, """
    ```@meta
    EditURL = "$PLOT_DOCS_URL"
    ```

    ### [Initialize](@id $pkgname-examples)

    ```@example $pkgname
    using Plots
    Plots.reset_defaults() # hide
    $(pkgname)()
    ```
    """)

    for (i,example) in enumerate(_examples)
        i in skip && continue
        # generate animations only for GR
        i in (2, 31) && pkgname != :gr && continue
        # write out the header, description, code block, and image link
        if !isempty(example.header)
            write(md, "### [$(example.header)](@id $pkgname-ref$i)\n")
        end
        write(md, "$(example.desc)\n")
        write(md, """
        ```@example $pkgname
        Plots.reset_defaults() # hide
        """)
        if pkgname ∈ (:unicodeplots, :inspectdr, :gaston)
            write(md, "using Logging; Logging.disable_logging(Logging.Warn) # hide\n")
        end
        for expr in example.exprs
            pretty_print_expr(md, expr)
        end
        if pkgname == :unicodeplots
            write(md, "show(current()) # hide\n")
        end
        if i in (2, 31)
            write(md, "gif(anim, \"anim_$(pkgname)_ex$i.gif\") # hide\n")
        end
        if pkgname ∈ (:plotly, :plotlyjs, :inspectdr, :gaston)
            write(md, "png(\"$(pkgname)_ex$i\") # hide\n")
        end
        write(md, "```\n")
        if pkgname ∈ (:plotly, :plotlyjs, :inspectdr, :gaston)
            write(md, "![]($(pkgname)_ex$i.png)\n")
        end
    end

    write(md, "- Supported arguments: $(markdown_code_to_string(collect(Plots.supported_attrs(pkg))))\n")
    write(md, "- Supported values for linetype: $(markdown_symbols_to_string(Plots.supported_seriestypes(pkg)))\n")
    write(md, "- Supported values for linestyle: $(markdown_symbols_to_string(Plots.supported_styles(pkg)))\n")
    write(md, "- Supported values for marker: $(markdown_symbols_to_string(Plots.supported_markers(pkg)))\n")

    write(md, "(Automatically generated: $(now()))")
    close(md)
end

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

# tables detailing the features that each backend supports

function make_support_df(allvals, func)
    vals = sort(allvals) # rows
    bs = sort(backends())
    bs = filter(be -> !(be in Plots._deprecated_backends), bs) # cols
    df = DataFrames.DataFrame(keys=vals)

    for b in bs
        b_supported_vals = ["" for _ in 1:length(vals)]
        for (i, val) in enumerate(vals)
            if func == Plots.supported_seriestypes
                stype = Plots.seriestype_supported(Plots._backend_instance(b), val)
                b_supported_vals[i] = stype == :native ? "✅" : (stype == :no ? "" : "🔼")
            else
                supported = func(Plots._backend_instance(b))

                b_supported_vals[i] = val in supported ? "✅" : ""
            end
        end
        df[!, b] = b_supported_vals
    end
    df
end

function generate_supported_markdown()
    md = open(joinpath(GENDIR, "supported.md"), "w")

    write(md, """
    ```@meta
    EditURL = "$PLOT_DOCS_URL"
    ```

    ## [Series Types](@id supported)

    Key:

    - ✅ the series type is natively supported by the backend.
    - 🔼 the series type is supported through series recipes.

    """)

    write(md, "```@raw html\n")
    write(md,
        to_html(
            make_support_df(
                Plots.all_seriestypes(),
                Plots.supported_seriestypes,
            )
        )
    )
    write(md, "\n```\n\n")

    supported_args =OrderedDict(
        "Keyword Arguments" => (Plots._all_args, Plots.supported_attrs),
        "Markers" => (Plots._allMarkers, Plots.supported_markers),
        "Line Styles" => (Plots._allStyles,  Plots.supported_styles),
        "Scales" => (Plots._allScales,  Plots.supported_scales)
    )

    for (header, args) in supported_args
        write(md, """

        ## $header

        """)

        write(md, "```@raw html\n")
        write(md, to_html(make_support_df(args...)))
        write(md, "\n```\n\n")
    end

    write(md, "\n(Automatically generated: $(now()))")
    close(md)
end

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

function make_attr_df(ktype::Symbol, defs::KW)
    n = length(defs)
    df = DataFrame(
        Attribute = fill("", n),
        Aliases = fill("", n),
        Default = fill("", n),
        Type = fill("", n),
        Description = fill("", n),
    )
    for (i, (k, def)) in enumerate(defs)
        desc = get(Plots._arg_desc, k, "")
        first_period_idx = findfirst(isequal('.'), desc)

        aliases = sort(collect(keys(filter(p -> p.second == k, Plots._keyAliases))))
        df.Attribute[i] = string(k)
        df.Aliases[i] = join(aliases, ", ")
        if first_period_idx !== nothing
            typedesc = desc[1:first_period_idx-1]
            desc = strip(desc[first_period_idx+1:end])
            df.Default[i] = show_default(def)
            df.Type[i] = string(typedesc)
            df.Description[i] = string(desc)
        end
    end
    sort!(df, [:Attribute])
    return df
end

show_default(x) = string("`", x, "`")
show_default(::Nothing) = "`nothing`"

const ATTRIBUTE_TEXTS = Dict(
    :Series => "These attributes apply to individual series (lines, scatters, heatmaps, etc)",
    :Plot => "These attributes apply to the full Plot. (A Plot contains a tree-like layout of Subplots)",
    :Subplot => "These attributes apply to settings for individual Subplots.",
    :Axis => """
    These attributes apply by default to all Axes in a Subplot (for example the `subplot[:xaxis]`).
    !!! info
        You can also specific the x, y, or z axis for each of these attributes by prefixing the attribute name with x, y, or z
        (for example `xmirror` only sets the mirror attribute for the x axis).
    """,
)

const ATTRIBUTE_DEFAULTS = Dict(
    :Series => Plots._series_defaults,
    :Plot => Plots._plot_defaults,
    :Subplot => Plots._subplot_defaults,
    :Axis => Plots._axis_defaults,
)

function generate_attr_markdown(c)
    # open the markdown file
    cstr = lowercase(string(c))
    attr_text = ATTRIBUTE_TEXTS[c]
    md = open(joinpath(GENDIR, "attributes_$cstr.md"), "w")

    write(md, """
    ```@meta
    EditURL = "$PLOT_DOCS_URL"
    ```
    ### $c

    $attr_text

    """)

    write(md, "```@raw html\n")
    write(md, to_html(make_attr_df(c, ATTRIBUTE_DEFAULTS[c])))
    write(md, "\n```\n\n")

    write(md, "\n(Automatically generated: $(now()))")
    close(md)
end

function generate_attr_markdown()
    for c in (:Series, :Plot, :Subplot, :Axis)
        generate_attr_markdown(c)
    end
end

function generate_graph_attr_markdown()
    md = open(joinpath(GENDIR, "graph_attributes.md"), "w")

    write(md, """
    ```@meta
    EditURL = "$PLOT_DOCS_URL"
    ```
    # [Graph Attributes](@id graph_attributes)

    Where possible, GraphRecipes will adopt attributes from Plots.jl to format visualizations.
    For example, the `linewidth` attribute from Plots.jl has the same effect in GraphRecipes.
    In order to give the user control over the layout of the graph visualization, GraphRecipes
    provides a number of keyword arguments (attributes). Here we describe those attributes
    alongside their default values.

    """)

    df = DataFrame(
        Attribute = [
            "dim",
            "T",
            "curves",
            "curvature_scalar",
            "root",
            "node_weights",
            "names",
            "fontsize",
            "nodeshape",
            "nodesize",
            "nodecolor",
            "x, y, z",
            "method",
            "func",
            "shorten",
            "axis_buffer",
            "layout_kw",
            "edgewidth",
            "edgelabel",
            "edgelabel_offset",
            "self_edge_size",
            "edge_label_box",
        ],
        Aliases = [
            "",
            "",
            "",
            "curvaturescalar, curvature",
            "",
            "nodeweights",
            "",
            "",
            "node_shape",
            "node_size",
            "marker_color",
            "x",
            "",
            "",
            "shorten_edge",
            "axisbuffer",
            "",
            "edge_width, ew",
            "edge_label, el",
            "edgelabeloffset, elo",
            "selfedgesize, ses",
            "edgelabelbox, edgelabel_box, elb",
        ],
        Default = [
            "2",
            "Float64",
            "true",
            "0.05",
            ":top",
            "nothing",
            "[]",
            "7",
            ":hexagon",
            "0.1",
            "1",
            "nothing",
            ":stress",
            "get(_graph_funcs, method, by_axis_local_stress_graph)",
            "0.0",
            "0.2",
            "Dict{Symbol,Any}()",
            "(s, d, w) -> 1",
            "nothing",
            "0.0",
            "0.1",
            "true",
        ],
        Description = [
            "The number of dimensions in the visualization.",
            "The data type for the coordinates of the graph nodes.",
            "Whether or not edges are curved. If `curves == true`, then the edge going from node \$s\$ to node \$d\$ will be defined by a cubic spline passing through three points: (i) node \$s\$, (ii) a point `p` that is distance `curvature_scalar` from the average of node \$s\$ and node \$d\$ and (iii) node \$d\$.",
            "A scalar that defines how much edges curve, see `curves` for more explanation.",
            "For displaying trees, choose from `:top`, `:bottom`, `:left`, `:right`. If you choose `:top`, then the tree will be plotted from the top down.",
            "The weight of the nodes given by a list of numbers. If `node_weights != nothing`, then the size of the nodes will be scaled by the `node_weights` vector.",
            "Names of the nodes given by a list of objects that can be parsed into strings. If the list is smaller than the number of nodes, then GraphRecipes will cycle around the list.",
            "Font size for the node labels and the edge labels.",
            "Shape of the nodes, choose from `:hexagon`, `:circle`, `:ellipse`, `:rect` or `:rectangle`.",
            "The size of nodes in the plot coordinates. Note that if `names` is not empty, then nodes will be scaled to fit the labels inside them.",
            "The color of the nodes. If `nodecolor` is an integer, then it will be taken from the current color pallette. Otherwise, the user can pass any color that would be recognised by the Plots `color` attribute.",
            "The coordinates of the nodes.",
            "The method that GraphRecipes uses to produce an optimal layout, choose from `:spectral`, `:sfdp`, `:circular`, `:shell`, `:stress`, `:spring`, `:tree`, `:buchheim`, `:arcdiagram` or `:chorddiagram`. See [NetworkLayout](https://github.com/JuliaGraphs/NetworkLayout.jl) for further details.",
            "A layout algorithm that can be passed in by the user.",
            "An amount to shorten edges by.",
            "Increase the `xlims` and `ylims`/`zlims` of the plot. Can be useful if part of the graph sits outside of the default view.",
            "A list of keywords to be passed to the layout algorithm, see [NetworkLayout](https://github.com/JuliaGraphs/NetworkLayout.jl) for a list of keyword arguments for each algorithm.",
            "The width of the edge going from \$s\$ to node \$d\$ with weight \$w\$.",
            "A dictionary of `(s, d) => label`, where `s` is an integer for the source node, `d` is an integer for the destiny node and `label` is the desired label for the given edge. Alternatively the user can pass a vector or a matrix describing the edge labels. If you use a vector or matrix, then either `missing`, `false`, `nothing`, `NaN` or `\"\"` values will not be displayed. In the case of multigraphs, triples can be used to define edges.",
            "The distance between edge labels and edges.",
            "The size of self edges.",
            "A box around edge labels that avoids intersections between edge labels and the edges that they are labeling.",
        ]
    )

    write(md, "```@raw html\n")
    write(md, to_html(df))
    write(md, "\n```\n\n")

    write(md, """
    ## Aliases
    Certain keyword arguments have aliases, so GraphRecipes "does what you mean, not
    what you say".

    So for example, `nodeshape=:rect` and `node_shape=:rect` are equivalent. To see the
    available aliases, type `GraphRecipes.graph_aliases`. If you are unhappy with the provided
    aliases, then you can add your own:
    ```julia
    using GraphRecipes, Plots

    push!(GraphRecipes.graph_aliases[:nodecolor],:nc)

    # These two calls produce the same plot, modulo some randomness in the layout.
    plot(graphplot([0 1; 0 0], nodecolor=:red), graphplot([0 1; 0 0], nc=:red))
    ```
    """)

    write(md, "\n(Automatically generated: $(now()))")
    close(md)
end

function generate_colorschemes_markdown()
    md = open(joinpath(GENDIR, "colorschemes.md"), "w")

    write(md, """
    ```@meta
    EditURL = "$PLOT_DOCS_URL"
    ```
    """)

    for line in readlines(normpath(@__DIR__, "..", "docs", "src", "colorschemes.md"))
        write(md, line)
        write(md, "\n")
    end

    write(md, """
    ## misc

    These colorschemes are not defined or provide different colors in ColorSchemes.jl
    They are kept for compatibility with Plots behavior before v1.1.0.

    """)
    write(md, "```@raw html\n")
    write(
        md,
        generate_colorschemes_table(
            [:default; sort(collect(keys(PlotUtils.MISC_COLORSCHEMES)))]
        )
    )
    write(md, "\n```\n\nThe following colorschemes are defined by ColorSchemes.jl.\n\n")
    for cs in ["cmocean", "scientific", "matplotlib", "colorbrewer", "gnuplot", "colorcet", "seaborn", "general"]
        ks = sort([k for (k, v) in PlotUtils.ColorSchemes.colorschemes if occursin(cs, v.category)])
        write(md, "\n\n## $cs\n\n```@raw html\n")
        write(md, generate_colorschemes_table(ks))
        write(md, "\n```\n\n")
    end

    close(md)
end

function colors_svg(cs, w, h)
    n = length(cs)
    ws = min(w / n, h)
    html = """
    <?xml version="1.0" encoding="UTF-8"?>
    <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
     "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
     <svg xmlns="http://www.w3.org/2000/svg" version="1.1"
          width="$(n * ws)mm" height="$(h)mm"
          viewBox="0 0 $n 1" preserveAspectRatio="none"
          shape-rendering="crispEdges" stroke="none">
    """
    for (i, c) in enumerate(cs)
        html *= """
        <rect width="$(ws)mm" height="$(h)mm" x="$(i-1)" y="0" fill="#$(hex(convert(RGB, c)))" />
        """
    end
    html *= "</svg>"
    return html
end

function generate_colorschemes_table(ks)
    extra_dir = get(ENV, "CI", "false") == "true" ? "../" : ""
    html = "<link type=\"text/css\" rel=\"stylesheet\" href=\"$(extra_dir)../assets/tables.css\" />

"
    for header in ["NAME", "palette(NAME)", "cgrad(NAME)"]
        html *= "$header	"
    end
    html *= "
"
    w, h = 60, 5
    for k in ks
        p = palette(k)
        cg = cgrad(k)[range(0, 1, length = 100)]
        cp = length(p) <= 100 ? color_list(p) : cg
        # cp7 = color_list(palette(k, 7))

        html *= "
:$k	"
        html *= colors_svg(cp, w, h)
        html *= "	"
        html *= colors_svg(cg, w, h)
        # html *= "	"
        # html *= colors_svg(cp7, 35, h)
        html *= "
"
    end
    html *= "
"
    return html
end

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

function to_html(df::DataFrames.AbstractDataFrame)
    cnames = DataFrames._names(df)
    extra_dir = get(ENV, "CI", "false") == "true" ? "../" : ""
    html = "<link type=\"text/css\" rel=\"stylesheet\" href=\"$(extra_dir)../assets/tables.css\" />

"
    for column_name in cnames
        html *= "$column_name	"
    end
    html *= "
"
    attrstr = " class=\"attr\""
    for row in 1:size(df,1)
        html *= "
"
        for (i,column_name) in enumerate(cnames)
            data = df[row, i]
            cell = data == nothing ? "" : string(data)
            attrstr = if i == 1
                " class=\"attr\""
            elseif i == length(cnames)
                " class=\"desc\""
            else
                ""
            end
            html *= "$(DataFrames.html_escape(cell))	"
        end
        html *= "
"
    end
    html *= "
"
    return html
end

end # module


text:
export get_independent_variables
export get_harmonic_equations
export is_rearranged
export slow_flow, slow_flow!
export _remove_brackets

show(eom::HarmonicEquation) = show_fields(eom)

"""
    harmonic_ansatz(eom::DifferentialEquation, time::Num; coordinates="Cartesian")

Expand each variable of `diff_eom` using the harmonics assigned to it with `time` as the time variable.
For each harmonic of each variable, instance(s) of `HarmonicVariable` are automatically created and named.

"""
function harmonic_ansatz(diff_eom::DifferentialEquation, time::Num)
    !is_harmonic(diff_eom, time) && error("The differential equation is not harmonic in ", time, " !")
    eqs = collect(values(diff_eom.equations))
    rules, vars = Dict(), []

    # keep count to label new variables
    uv_idx = 1
    a_idx = 1

    for nvar in get_variables(diff_eom) # sum over natural variables
        to_substitute = 0 # combine all the subtitution rules for var
        for ω in diff_eom.harmonics[nvar]
            if !isequal(ω, 0) # nonzero harmonic - create u,v
                rule_u, hvar_u = _create_harmonic_variable(nvar, ω, time, "u", new_symbol="u"*string(uv_idx))
                rule_v, hvar_v = _create_harmonic_variable(nvar, ω, time, "v", new_symbol="v"*string(uv_idx))
                rule = rule_u + rule_v
                uv_idx += 1
                append!(vars, [hvar_u, hvar_v])
            else # zero harmonic - create a
                rule, hvar = _create_harmonic_variable(nvar, ω, time, "a", new_symbol="a"*string(a_idx))
                a_idx += 1
                append!(vars, [hvar])
            end
            to_substitute += rule
        end

        rules[nvar] = to_substitute # total sub rule for nvar
    end
    eqs = substitute_all(eqs, rules)
    HarmonicEquation(eqs, Vector{HarmonicVariable}(vars), diff_eom)
end

function slow_flow!(eom::HarmonicEquation; fast_time::Num, slow_time::Num, degree=2)
    eom.equations = expand_derivatives.(eom.equations) # expand all the derivatives

    # fast_time => slow_time for derivatives up to degree-1
    replace0 = [var => substitute_all(var, fast_time => slow_time) for var in get_variables(eom)]
    replace_degrees = [Differential(fast_time)^deg => Differential(slow_time)^deg for deg in 1:degree-1]
    replace = flatten([replace0, replace_degrees])

    # degree derivatives are removed
    drop = [d(var, fast_time, degree) => 0 for var in get_variables(eom)]

    eom.equations = substitute_all(substitute_all(eom.equations, drop), replace)
    eom.variables = substitute_all(eom.variables, replace)
end

"""
    slow_flow(eom::HarmonicEquation; fast_time::Num, slow_time::Num, degree=2)

Removes all derivatives w.r.t `fast_time` (and their products) in `eom` of power `degree`.
In the remaining derivatives, `fast_time` is replaced by `slow_time`.
"""
function slow_flow(eom::HarmonicEquation; fast_time::Num, slow_time::Num, degree=2)::HarmonicEquation
    new_eq = deepcopy(eom)
    slow_flow!(new_eq, fast_time=fast_time, slow_time=slow_time, degree=degree)
    new_eq
end

#   Drop powers of `var` of degree >= `deg` from the equation set in `eom`.
function drop_powers(eom::HarmonicEquation, var, deg::Integer)
    new_eom = deepcopy(eom)
    new_eom.equations = drop_powers(eom.equations, var, deg)
    new_eom
end

"Rearrange an equation system such that the field equations is equal to the vector specified in new_lhs"
function rearrange!(eom::HarmonicEquation, new_rhs::Vector{Num})
    soln = Symbolics.solve_for(eom.equations, new_rhs, simplify=false,check=true)
    eom.equations = soln .~ new_rhs
    return nothing
end

function rearrange(eom::HarmonicEquation, new_rhs::Vector{Num})
    new_eom = deepcopy(eom)
    rearrange!(new_eom, new_rhs)
    return new_eom
end

"""
$(TYPEDSIGNATURES)
Check if `eom` is rearranged to the standard form, such that the derivatives of the variables are on one side.
"""
function is_rearranged(eom::HarmonicEquation)
    tvar = get_independent_variables(eom)[1]

    # Hopf-containing equations are arranged by construstion (impossible to check)
    isequal(getfield.(eom.equations, :rhs), d(get_variables(eom), tvar)) || in("Hopf", getfield.(eom.variables, :type))
end

"""
$(TYPEDSIGNATURES)
Rearrange `eom` to the standard form, such that the derivatives of the variables are on one side.
"""
function rearrange_standard(eom::HarmonicEquation)
    tvar = get_independent_variables(eom)[1]
    dvars =  d(get_variables(eom), tvar)
    rearrange(eom, dvars)
end

"""
$(TYPEDSIGNATURES)
Get the internal symbols of the independent variables of `eom`.
"""
function get_variables(eom::HarmonicEquation)
    return flatten(get_variables.(eom.variables))
end

"Get the parameters (not time nor variables) of a HarmonicEquation"
function _parameters(eom::HarmonicEquation)
    all_symbols = flatten([cat(get_variables(eq.lhs), get_variables(eq.rhs), dims=1) for eq in eom.equations])
    # subtract the set of independent variables (i.e., time) from all free symbols
    setdiff(all_symbols, get_variables(eom), get_independent_variables(eom))
end

###
# Extending Symbolics.jl's simplify and substitute
###

"Apply `rules` to both `equations` and `variables` field of `eom`"
function substitute_all(eom::HarmonicEquation, rules::Union{Dict, Pair})
    new_eom = deepcopy(eom)
    new_eom.equations = expand_derivatives.(substitute_all(eom.equations, rules))
    new_eom
end

"Simplify the equations in HarmonicEquation."
function simplify!(eom::HarmonicEquation)
    eom.equations = [simplify(eq) for eq in eom.equations]
end

"""
$(TYPEDSIGNATURES)
Return the independent variables (typically time) of `eom`.
"""
function get_independent_variables(eom::HarmonicEquation)::Vector{Num}
    dynamic_vars = flatten(getfield.(eom.variables, Symbol("symbol")))
    return flatten(unique([SymbolicUtils.arguments(var.val) for var in dynamic_vars]))
end

"""
$(TYPEDSIGNATURES)
Extract the Fourier components of `eom` corresponding to the harmonics specified in `eom.variables`.
For each non-zero harmonic of each variable, 2 equations are generated (cos and sin Fourier coefficients).
For each zero (constant) harmonic, 1 equation is generated
`time` does not appear in the resulting equations anymore.

Underlying assumption: all time-dependences are harmonic.
"""
function fourier_transform(eom::HarmonicEquation, time::Num)
    new_eom = deepcopy(eom)
    fourier_transform!(new_eom, time)
    return new_eom
end

function fourier_transform!(eom::HarmonicEquation, time::Num)
    avg_eqs = Vector{Equation}(undef, length(eom.variables))

    # loop over the HarmonicVariables, each generates one equation
    for (i, hvar) in enumerate(eom.variables)
        # find the equation belonging to this variable
        eq_idx = findall(x -> isequal(x, hvar.natural_variable), collect(keys(eom.natural_equation.equations)))
        eq = eom.equations[eq_idx...]
        # "type" is usually "u" or "v" (harmonic) or ["a"] (zero-harmonic)
        if hvar.type == "u"
            avg_eqs[i] = fourier_cos_term(eq, hvar.ω, time)
        elseif hvar.type == "v"
            avg_eqs[i] = fourier_sin_term(eq, hvar.ω, time)
        elseif hvar.type == "a"
            avg_eqs[i] = fourier_cos_term(eq, 0, time) # pick out the constants
        end
    end

    eom.equations = avg_eqs
end

"""
    get_harmonic_equations(diff_eom::DifferentialEquation; fast_time=nothing, slow_time=nothing)

Apply the harmonic ansatz, followed by the slow-flow, Fourier transform and dropping
higher-order derivatives to obtain
a set of ODEs (the harmonic equations) governing the harmonics of `diff_eom`.

The harmonics evolve in `slow_time`, the oscillating terms themselves in `fast_time`.
If no input is used, a variable T is defined for `slow_time` and `fast_time` is taken as the independent variable
of `diff_eom`.

By default, all products of order > 1 of `slow_time`-derivatives are dropped,
which means the equations are linear in the time-derivatives.

# Example
```julia-repl
julia> @variables t, x(t), ω0, ω, F;

# enter the simple harmonic oscillator
julia> diff_eom = DifferentialEquation( d(x,t,2) + ω0^2 * x ~ F *cos(ω*t), x);

# expand x in the harmonic ω
julia> add_harmonic!(diff_eom, x, ω);

# get equations for the harmonics evolving in the slow time T
julia> harmonic_eom = get_harmonic_equations(diff_eom)

A set of 2 harmonic equations
Variables: u1(T), v1(T)
Parameters: ω0, ω, F

Harmonic ansatz:
x(t) = u1*cos(ωt) + v1*sin(ωt)

Harmonic equations:

(ω0^2)*u1(T) + (2//1)*ω*Differential(T)(v1(T)) - (ω^2)*u1(T) ~ F

(ω0^2)*v1(T) - (ω^2)*v1(T) - (2//1)*ω*Differential(T)(u1(T)) ~ 0
```

"""
function get_harmonic_equations(diff_eom::DifferentialEquation; fast_time=nothing, slow_time=nothing)

    slow_time = isnothing(slow_time) ? (@variables T; T) : slow_time
    fast_time = isnothing(fast_time) ? get_independent_variables(diff_eom)[1] : fast_time

    all(isempty.(values(diff_eom.harmonics))) && error("No harmonics specified!")
    eom = harmonic_ansatz(diff_eom, fast_time); # substitute trig functions into the differential equation
    eom = slow_flow(eom, fast_time=fast_time, slow_time=slow_time); # drop 2nd order time derivatives
    fourier_transform!(eom, fast_time); # perform averaging over the frequencies originally specified in dEOM
    ft_eom_simplified = drop_powers(eom, d(get_variables(eom), slow_time), 2); # drop higher powers of the first-order derivatives
    return ft_eom_simplified
end

"Rearrange `eq` to have zero on the right-hand-side."
_set_zero_rhs(eq::Equation) = eq.lhs - eq.rhs ~ 0
_set_zero_rhs(eqs::Vector{Equation}) = [_set_zero_rhs(eq) for eq in eqs]

_remove_brackets(var::Num) = declare_variable(var_name(var))
_remove_brackets(vars::Vector{Num}) = _remove_brackets.(vars)
_remove_brackets(hv::HarmonicVariable) = _remove_brackets(hv.symbol)

"Returns the equation system in `eom`, dropping all argument brackets (i.e., u(T) becomes u)."
function _remove_brackets(eom::HarmonicEquation)
    variable_rules = [var => _remove_brackets(var) for var in get_variables(eom)]
    equations_lhs = Num.(getfield.(eom.equations, :lhs)  - getfield.(eom.equations, :rhs))
    substitute_all(equations_lhs, variable_rules)
end

function drop_powers(eom::HarmonicEquation, terms, deg::Int)
    new_eom = deepcopy(eom)
    new_eom.equations = drop_powers(eom.equations, terms, deg)
    new_eom
end

ext: jl


text:
using Pkg
Pkg.build()
Pkg.test("Coverage"; coverage=true)


text:
@testset "context" begin

ctx = current_context()
dev = current_device()

synchronize(ctx)

let ctx2 = CuContext(dev)
    @test ctx2 == current_context()    # ctor implicitly pushes
    activate(ctx)
    @test ctx == current_context()

    @test device(ctx2) == dev

    CUDA.unsafe_destroy!(ctx2)
end

let global_ctx2 = nothing
    CuContext(dev) do ctx2
        @test ctx2 == current_context()
        @test ctx != ctx2
        global_ctx2 = ctx2
    end
    @test !CUDA.isvalid(global_ctx2)
    @test ctx == current_context()

    @test device(ctx) == dev
    @test current_device() == dev
    device_synchronize()
end

end

# FIXME: these trample over our globally-managed context

# @testset "primary context" begin

# pctx = CuPrimaryContext(device())

# @test !isactive(pctx)
# unsafe_reset!(pctx)
# @test !isactive(pctx)

# @test flags(pctx) == 0
# setflags!(pctx, CUDA.CTX_SCHED_BLOCKING_SYNC)
# @test flags(pctx) == CUDA.CTX_SCHED_BLOCKING_SYNC

# let global_ctx = nothing
#     CuContext(pctx) do ctx
#         @test CUDA.isvalid(ctx)
#         @test isactive(pctx)
#         global_ctx = ctx
#     end
#     @test !isactive(pctx)
#     @test !CUDA.isvalid(global_ctx)
# end

# CuContext(pctx) do ctx
#     @test CUDA.isvalid(ctx)
#     @test isactive(pctx)

#     unsafe_reset!(pctx)

#     @test !isactive(pctx)
#     @test !CUDA.isvalid(ctx)
# end

# let
#     @test !isactive(pctx)

#     ctx1 = CuContext(pctx)
#     @test isactive(pctx)
#     @test CUDA.isvalid(ctx1)

#     unsafe_reset!(pctx)
#     @test !isactive(pctx)
#     @test !CUDA.isvalid(ctx1)
#     CUDA.valid_contexts

#     ctx2 = CuContext(pctx)
#     @test isactive(pctx)
#     @test !CUDA.isvalid(ctx1)
#     @test CUDA.isvalid(ctx2)

#     unsafe_reset!(pctx)
# end

# end

@testset "cache config" begin

config = cache_config()

cache_config!(CUDA.FUNC_CACHE_PREFER_L1)
@test cache_config() == CUDA.FUNC_CACHE_PREFER_L1

cache_config!(config)

end

@testset "shmem config" begin

config = shmem_config()

shmem_config!(CUDA.SHARED_MEM_CONFIG_FOUR_BYTE_BANK_SIZE)
@test shmem_config() == CUDA.SHARED_MEM_CONFIG_FOUR_BYTE_BANK_SIZE

shmem_config!(config)

end

@testset "limits" begin

lim = limit(CUDA.LIMIT_DEV_RUNTIME_SYNC_DEPTH)

lim += 1
limit!(CUDA.LIMIT_DEV_RUNTIME_SYNC_DEPTH, lim)
@test lim == limit(CUDA.LIMIT_DEV_RUNTIME_SYNC_DEPTH)

limit!(CUDA.LIMIT_DEV_RUNTIME_SYNC_DEPTH, lim)

end

############################################################################################

@testset "devices" begin

dev = device()

@test name(dev) isa String
@test uuid(dev) isa Base.UUID
totalmem(dev)
attribute(dev, CUDA.DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK)
@test warpsize(dev) == 32
capability(dev)
@grab_output show(stdout, "text/plain", dev)

@test eval(Meta.parse(repr(dev))) == dev

@test eltype(devices()) == CuDevice
@grab_output show(stdout, "text/plain", CUDA.DEVICE_CPU)
@grab_output show(stdout, "text/plain", CUDA.DEVICE_INVALID)

@test length(devices()) == ndevices()
@grab_output show(stdout, "text/plain", devices())

end

############################################################################################

@testset "errors" begin

let
    ex = CuError(CUDA.SUCCESS)
    @test CUDA.name(ex) == "SUCCESS"
    @test CUDA.description(ex) == "no error"
    @test eval(Meta.parse(repr(ex))) == ex

    io = IOBuffer()
    showerror(io, ex)
    str = String(take!(io))

    @test occursin("0", str)
    @test occursin("no error", str)
end

let
    ex = CuError(CUDA.SUCCESS, "foobar")

    io = IOBuffer()
    showerror(io, ex)
    str = String(take!(io))

    @test occursin("foobar", str)
end

end

############################################################################################

@testset "events" begin

let
    start = CuEvent()
    stop = CuEvent()
    @test start != stop

    record(start)
    record(stop)
    synchronize(stop)

    @test elapsed(start, stop) > 0
end

@test (CUDA.@elapsed begin end) > 0

CuEvent(CUDA.EVENT_BLOCKING_SYNC)
CuEvent(CUDA.EVENT_BLOCKING_SYNC | CUDA.EVENT_DISABLE_TIMING)

@testset "stream wait" begin
    event  = CuEvent()
    stream = CuStream()

    CUDA.record(event, stream)

    CUDA.wait(event)
    synchronize()
end

@testset "event query" begin
    event  = CuEvent()
    @test CUDA.query(event) == true
end

end

############################################################################################

@testset "execution" begin

let
    # test outer CuDim3 constructors
    @test CUDA.CuDim3((Cuint(4),Cuint(3),Cuint(2))) == CUDA.CuDim3(Cuint(4),Cuint(3),Cuint(2))
    @test CUDA.CuDim3((Cuint(3),Cuint(2)))          == CUDA.CuDim3(Cuint(3),Cuint(2),Cuint(1))
    @test CUDA.CuDim3((Cuint(2),))                  == CUDA.CuDim3(Cuint(2),Cuint(1),Cuint(1))
    @test CUDA.CuDim3(Cuint(2))                     == CUDA.CuDim3(Cuint(2),Cuint(1),Cuint(1))

    # outer constructor should type convert
    @test CUDA.CuDim3(2)       == CUDA.CuDim3(Cuint(2),Cuint(1),Cuint(1))
    @test_throws InexactError CUDA.CuDim3(typemax(Int64))

    # CuDim type alias should accept conveniently-typed dimensions
    @test isa(2,        CUDA.CuDim)
    @test isa((2,),     CUDA.CuDim)
    @test isa((2,2),    CUDA.CuDim)
    @test isa((2,2,2),  CUDA.CuDim)
    @test isa(Cuint(2), CUDA.CuDim)
end

@testset "device" begin

let
    md = CuModuleFile(joinpath(@__DIR__, "ptx/dummy.ptx"))
    dummy = CuFunction(md, "dummy")

    # different cudacall syntaxes
    cudacall(dummy, Tuple{})
    cudacall(dummy, Tuple{}; threads=1)
    cudacall(dummy, Tuple{}; threads=1, blocks=1)
    cudacall(dummy, Tuple{}; threads=1, blocks=1, shmem=0)
    cudacall(dummy, Tuple{}; threads=1, blocks=1, shmem=0, stream=stream())
    cudacall(dummy, Tuple{}; threads=1, blocks=1, shmem=0, stream=stream(), cooperative=false)
    cudacall(dummy, ())
    cudacall(dummy, (); threads=1, blocks=1, shmem=0, stream=stream(), cooperative=false)

    # different launch syntaxes
    CUDA.launch(dummy)
    CUDA.launch(dummy; threads=1)
    CUDA.launch(dummy; threads=1, blocks=1)
    CUDA.launch(dummy; threads=1, blocks=1, shmem=0)
    CUDA.launch(dummy; threads=1, blocks=1, shmem=0, stream=stream())
    CUDA.launch(dummy; threads=1, blocks=1, shmem=0, stream=stream(), cooperative=false)
end

let
    md = CuModuleFile(joinpath(@__DIR__, "ptx/vectorops.ptx"))
    vadd = CuFunction(md, "vadd")
    vmul = CuFunction(md, "vmul")
    vsub = CuFunction(md, "vsub")
    vdiv = CuFunction(md, "vdiv")

    a = rand(Float32, 10)
    b = rand(Float32, 10)
    ad = CuArray(a)
    bd = CuArray(b)

    # Addition
    let
        c = zeros(Float32, 10)
        c_d = CuArray(c)
        cudacall(vadd,
                 (CuPtr{Cfloat},CuPtr{Cfloat},CuPtr{Cfloat}), ad, bd, c_d;
                 threads=10)
        c = Array(c_d)
        @test c ≈ a+b
    end

    # Subtraction
    let
        c = zeros(Float32, 10)
        c_d = CuArray(c)
        cudacall(vsub,
                 (CuPtr{Cfloat},CuPtr{Cfloat},CuPtr{Cfloat}), ad, bd, c_d;
                 threads=10)
        c = Array(c_d)
        @test c ≈ a-b
    end

    # Multiplication
    let
        c = zeros(Float32, 10)
        c_d = CuArray(c)
        cudacall(vmul,
                 (CuPtr{Cfloat},CuPtr{Cfloat},CuPtr{Cfloat}), ad, bd, c_d;
                 threads=10)
        c = Array(c_d)
        @test c ≈ a.*b
    end

    # Division
    let
        c = zeros(Float32, 10)
        c_d = CuArray(c)
        cudacall(vdiv,
                 (CuPtr{Cfloat},CuPtr{Cfloat},CuPtr{Cfloat}), ad, bd, c_d;
                 threads=10)
        c = Array(c_d)
        @test c ≈ a./b
    end
end

end

@testset "host" begin
    c = Condition()
    CUDA.launch() do
        notify(c)
    end
    wait(c)
end

@testset "attributes" begin

md = CuModuleFile(joinpath(@__DIR__, "ptx/dummy.ptx"))
dummy = CuFunction(md, "dummy")

val = attributes(dummy)[CUDA.FUNC_ATTRIBUTE_SHARED_SIZE_BYTES]

if CUDA.version() >= v"9.0"
    attributes(dummy)[CUDA.FUNC_ATTRIBUTE_MAX_DYNAMIC_SHARED_SIZE_BYTES] = val
end

end

end

############################################################################################

@testset "graph" begin

let A = CUDA.zeros(Int, 1)
    # ensure compilation
    A .+= 1
    @test Array(A) == [1]

    graph = capture() do
        @test is_capturing()
        A .+= 1
    end
    @test Array(A) == [1]

    exec = instantiate(graph)
    CUDA.launch(exec)
    @test Array(A) == [2]

    graph′ = capture() do
        A .+= 2
    end

    update(exec, graph′)
    CUDA.launch(exec)
    @test Array(A) == [4]
end

let A = CUDA.zeros(Int, 1)
    function iteration(A, val)
        # custom kernel to force compilation on the first iteration
        function kernel(a, val)
            a[] += val
            return
        end
        @cuda kernel(A, val)
        return
    end

    for i in 1:2
        @captured iteration(A, i)
    end
    @test Array(A) == [3]
end

end

############################################################################################

@testset "memory" begin

let
    a,b = Mem.info()
    # NOTE: actually testing this is pretty fragile on CI
    #=@test a == =# CUDA.available_memory()
    #=@test b == =# CUDA.total_memory()
end

# dummy data
T = UInt32
N = 5
data = rand(T, N)
nb = sizeof(data)

# buffers are untyped, so we use a convenience function to get a typed pointer
# we prefer to return a device pointer (for managed buffers) to maximize CUDA coverage
typed_pointer(buf::Union{Mem.Device, Mem.Unified}, T) = convert(CuPtr{T}, buf)
typed_pointer(buf::Mem.Host, T)                       = convert(Ptr{T},   buf)

# allocations and copies
for srcTy in [Mem.Device, Mem.Host, Mem.Unified],
    dstTy in [Mem.Device, Mem.Host, Mem.Unified]

    dummy = Mem.alloc(srcTy, 0)
    Mem.free(dummy)

    src = Mem.alloc(srcTy, nb)
    unsafe_copyto!(typed_pointer(src, T), pointer(data), N)

    dst = Mem.alloc(dstTy, nb)
    unsafe_copyto!(typed_pointer(dst, T), typed_pointer(src, T), N)

    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), typed_pointer(dst, T), N)

    @test data == ref

    if isa(src, Mem.Device) || isa(src, Mem.Unified)
        Mem.set!(typed_pointer(src, T), zero(T), N)
    end

    # test the memory-type attribute
    if isa(src, Mem.Device)
        @test CUDA.memory_type(typed_pointer(src, T)) == CUDA.MEMORYTYPE_DEVICE
    elseif isa(src, Mem.Host)
        @test CUDA.memory_type(convert(Ptr{T}, src)) == CUDA.MEMORYTYPE_HOST
    elseif isa(src, Mem.Unified)
        # unified memory can reside in either place
        # FIXME: does this depend on the current migration, or on the configuration?
        @test CUDA.memory_type(convert(CuPtr{T}, src)) == CUDA.MEMORYTYPE_HOST ||
              CUDA.memory_type(convert(CuPtr{T}, src)) == CUDA.MEMORYTYPE_DEVICE ||
        @test CUDA.memory_type(convert(CuPtr{T}, src)) == CUDA.memory_type(convert(Ptr{T}, src))
    end

    # test device with context in which pointer was allocated.
    @test device(typed_pointer(src, T)) == device()
    if !CUDA.has_stream_ordered(device())
        # NVIDIA bug #3319609
        @test CuContext(typed_pointer(src, T)) == context()
    end

    # test the is-managed attribute
    if isa(src, Mem.Unified)
        @test CUDA.is_managed(convert(Ptr{T}, src))
        @test CUDA.is_managed(convert(CuPtr{T}, src))
    else
        @test !CUDA.is_managed(typed_pointer(src, T))
    end
    # Test conversion to Ptr throwing an error
    if isa(src, Mem.Device)
        @test_throws ArgumentError convert(Ptr, src)
    end

    @grab_output show(stdout, src)
    @grab_output show(stdout, dst)
    Mem.free(src)
    Mem.free(dst)
end

# pointer attributes
let
    src = Mem.alloc(Mem.Device, nb)

    attribute!(typed_pointer(src, T), CUDA.POINTER_ATTRIBUTE_SYNC_MEMOPS, 0)

    Mem.free(src)
end

# asynchronous operations
let
    src = Mem.alloc(Mem.Device, nb)

    unsafe_copyto!(typed_pointer(src, T), pointer(data), N; async=true)

    Mem.set!(typed_pointer(src, T), zero(T), N; stream=stream())

    Mem.free(src)
end

# pinned memory
let
    # create a pinned and mapped buffer
    src = Mem.alloc(Mem.Host, nb, Mem.HOSTALLOC_DEVICEMAP)

    # get the CPU address and copy some data to the buffer
    cpu_ptr = convert(Ptr{T}, src)
    @test CUDA.memory_type(cpu_ptr) == CUDA.MEMORYTYPE_HOST
    unsafe_copyto!(cpu_ptr, pointer(data), N)

    # get the GPU address and copy back the data
    gpu_ptr = convert(CuPtr{T}, src)
    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), gpu_ptr, N)
    @test ref == data

    Mem.free(src)
    # NOTE: don't free dst, it's just a mapped pointer
end

# pinned memory with existing memory
if attribute(device(), CUDA.DEVICE_ATTRIBUTE_HOST_REGISTER_SUPPORTED) != 0
    # register a pinned and mapped buffer
    src = Mem.register(Mem.Host, pointer(data), nb, Mem.HOSTREGISTER_DEVICEMAP)

    # get the GPU address and copy back the data
    gpu_ptr = convert(CuPtr{T}, src)
    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), gpu_ptr, N)
    @test ref == data

    Mem.unregister(src)
end

# unified memory
let
    src = Mem.alloc(Mem.Unified, nb)

    @test_throws BoundsError Mem.prefetch(src, 2*nb; device=CUDA.DEVICE_CPU)
    # FIXME: prefetch doesn't work on some CI devices, unsure why.
    @test_skip Mem.prefetch(src, nb; device=CUDA.DEVICE_CPU)
    Mem.advise(src, Mem.ADVISE_SET_READ_MOSTLY)

    # get the CPU address and copy some data
    cpu_ptr = convert(Ptr{T}, src)
    unsafe_copyto!(cpu_ptr, pointer(data), N)

    # get the GPU address and copy back data
    gpu_ptr = convert(CuPtr{T}, src)
    ref = Array{T}(undef, N)
    unsafe_copyto!(pointer(ref), gpu_ptr, N)
    @test ref == data

    Mem.free(src)
end

# 3d memcpy
let
    # TODO: use cuMemAllocPitch (and put pitch in buffer?) to actually get benefit from this

    data = collect(reshape(1:27, 3, 3, 3))

    dst = Mem.alloc(Mem.Device, sizeof(data))
    Mem.unsafe_copy3d!(typed_pointer(dst, Int), Mem.Device, pointer(data), Mem.Host, length(data))

    check = zeros(Int, size(data))
    Mem.unsafe_copy3d!(pointer(check), Mem.Host, typed_pointer(dst, Int), Mem.Device, length(data))

    @test check == data

    Mem.free(dst)
end
let
    # copying an x-z plane of a 3-D array

    T = Int
    nx, ny, nz = 4, 4, 4
    data = collect(reshape(1:(nx*nz), nx, nz))
    dst = Mem.alloc(Mem.Device, nx*ny*nz*sizeof(data))

    # host to device
    Mem.unsafe_copy3d!(typed_pointer(dst, T), Mem.Device, pointer(data), Mem.Host,
                       nx, 1, nz;
                       dstPos=(1,2,1),
                       srcPitch=nx*sizeof(T), srcHeight=1,
                       dstPitch=nx*sizeof(T), dstHeight=ny)

    # copy back
    check = zeros(T, size(data))
    Mem.unsafe_copy3d!(pointer(check), Mem.Host, typed_pointer(dst, T), Mem.Device,
                       nx, 1, nz;
                       srcPos=(1,2,1),
                       srcPitch=nx*sizeof(T), srcHeight=ny,
                       dstPitch=nx*sizeof(T), dstHeight=1)

    @test check == data

    # copy back into a 3-D array
    check2 = zeros(T, nx, ny, nz)
    Mem.unsafe_copy3d!(pointer(check2), Mem.Host, typed_pointer(dst, T), Mem.Device,
                       nx, 1, nz;
                       srcPos=(1,2,1),
                       dstPos=(1,2,1),
                       srcPitch=nx*sizeof(T), srcHeight=ny,
                       dstPitch=nx*sizeof(T), dstHeight=ny)
    @test check2[:,2,:] == data
end
let
    # copying an y-z plane of a 3-D array

    T = Int
    nx, ny, nz = 4, 4, 4
    data = collect(reshape(1:(ny*nz), ny, nz))
    dst = Mem.alloc(Mem.Device, nx*ny*nz*sizeof(data))

    # host to device
    Mem.unsafe_copy3d!(typed_pointer(dst, T), Mem.Device, pointer(data), Mem.Host,
                       1, ny, nz;
                       dstPos=(2,1,1),
                       srcPitch=1*sizeof(T), srcHeight=ny,
                       dstPitch=nx*sizeof(T), dstHeight=ny)

    # copy back
    check = zeros(T, size(data))
    Mem.unsafe_copy3d!(pointer(check), Mem.Host, typed_pointer(dst, T), Mem.Device,
                       1, ny, nz;
                       srcPos=(2,1,1),
                       srcPitch=nx*sizeof(T), srcHeight=ny,
                       dstPitch=1*sizeof(T), dstHeight=ny)

    @test check == data

    # copy back into a 3-D array
    check2 = zeros(T, nx, ny, nz)
    Mem.unsafe_copy3d!(pointer(check2), Mem.Host, typed_pointer(dst, T), Mem.Device,
                       1, ny, nz;
                       srcPos=(2,1,1),
                       dstPos=(2,1,1),
                       srcPitch=nx*sizeof(T), srcHeight=ny,
                       dstPitch=nx*sizeof(T), dstHeight=ny)
    @test check2[2,:,:] == data
end
let
    # JuliaGPU/CUDA.jl#863: wrong offset calculation
    nx, ny, nz = 1, 2, 1

    A = zeros(Int, nx, nz)
    B = CuArray(reshape([1:(nx*ny*nz)...], (nx, ny, nz)))
    Mem.unsafe_copy3d!(
        pointer(A), Mem.Host, pointer(B), Mem.Device,
        nx, 1, nz;
        srcPos=(1,2,1),
        srcPitch=nx*sizeof(A[1]), srcHeight=ny,
        dstPitch=nx*sizeof(A[1]), dstHeight=1
    )

    @test A == Array(B)[:,2,:]
end

# pinned memory with existing memory
if attribute(device(), CUDA.DEVICE_ATTRIBUTE_HOST_REGISTER_SUPPORTED) != 0
    hA = rand(UInt8, 512)
    @test !CUDA.is_pinned(pointer(hA))
    Mem.pin(hA)
    @test CUDA.is_pinned(pointer(hA))

    # make sure we can double-pin
    Mem.pin(hA)

    # memory copies on pinned memory behave differently, so test that code path
    dA = CUDA.rand(UInt8, 512)
    copyto!(dA, hA)
    copyto!(hA, dA)
end

end

############################################################################################

@testset "module" begin

let
    md = CuModuleFile(joinpath(@__DIR__, "ptx/vadd.ptx"))

    vadd = CuFunction(md, "vadd")
end

# comparisons
let
    md = CuModuleFile(joinpath(@__DIR__, "ptx/vectorops.ptx"))
    vadd = CuFunction(md, "vadd")
    @test vadd == vadd
    @test vadd != CuFunction(md, "vdiv")
end

let
    f = open(joinpath(@__DIR__, "ptx/vadd.ptx"))
    ptx = read(f, String)
    close(f)

    md = CuModule(ptx)
    vadd = CuFunction(md, "vadd")

    md2 = CuModuleFile(joinpath(@__DIR__, "ptx/vadd.ptx"))
    @test md != md2
end

@test_throws_cuerror CUDA.ERROR_INVALID_IMAGE CuModule("foobar")

@testset "globals" begin
    md = CuModuleFile(joinpath(@__DIR__, "ptx/global.ptx"))

    var = CuGlobal{Int32}(md, "foobar")
    @test eltype(var) == Int32
    @test eltype(typeof(var)) == Int32

    @test_throws ArgumentError CuGlobal{Int64}(md, "foobar")

    var[] = Int32(42)
    @test var[] == Int32(42)
end

@testset "linker" begin
    link = CuLink()
    @test link == link
    @test link != CuLink()

    # PTX code
    open(joinpath(@__DIR__, "ptx/empty.ptx")) do f
        add_data!(link, "vadd_parent", read(f, String))
    end
    @test_throws ArgumentError add_data!(link, "vadd_parent", "\0")

    # object code
    # TODO: test with valid object code
    # NOTE: apparently, on Windows cuLinkAddData! _does_ accept object data containing \0
    if !Sys.iswindows()
        @test_throws_cuerror CUDA.ERROR_UNKNOWN add_data!(link, "vadd_parent", UInt8[0])
    end
end

@testset "error log" begin
    @test_throws_message contains("ptxas fatal") CuError CuModule(".version 3.1")

    link = CuLink()
    @test_throws_message contains("ptxas fatal") CuError add_data!(link, "dummy", ".version 3.1")
end

let
    link = CuLink()
    add_file!(link, joinpath(@__DIR__, "ptx/vadd_child.ptx"), CUDA.JIT_INPUT_PTX)
    open(joinpath(@__DIR__, "ptx/vadd_parent.ptx")) do f
        add_data!(link, "vadd_parent", read(f, String))
    end

    obj = complete(link)
    md = CuModule(obj)

    vadd = CuFunction(md, "vadd")

    options = Dict{CUDA.CUjit_option,Any}()
    options[CUDA.JIT_GENERATE_LINE_INFO] = true

    md = CuModule(obj, options)
    vadd = CuFunction(md, "vadd")
end

end

############################################################################################

@testset "occupancy" begin

let md = CuModuleFile(joinpath(@__DIR__, "ptx/dummy.ptx"))
    dummy = CuFunction(md, "dummy")

    active_blocks(dummy, 1)
    active_blocks(dummy, 1; shmem=64)

    occupancy(dummy, 1)
    occupancy(dummy, 1; shmem=64)

    launch_configuration(dummy)
    launch_configuration(dummy; shmem=64)
    launch_configuration(dummy; shmem=64, max_threads=64)

    let cb_calls = 0
        launch_configuration(dummy; shmem=threads->(cb_calls += 1; 0))
        @test cb_calls > 0
    end
end

end

############################################################################################

@testset "pool" begin

dev = device()
if CUDA.version() >= v"11.2" && attribute(dev, CUDA.DEVICE_ATTRIBUTE_MEMORY_POOLS_SUPPORTED) == 1

pool = memory_pool(dev)

pool2 = CuMemoryPool(dev)
@test pool2 != pool
memory_pool!(dev, pool2)
@test pool2 == memory_pool(dev)
@test pool2 != default_memory_pool(dev)

memory_pool!(dev, pool)
@test pool == memory_pool(dev)

@test attribute(UInt64, pool2, CUDA.MEMPOOL_ATTR_RELEASE_THRESHOLD) == 0
attribute!(pool2, CUDA.MEMPOOL_ATTR_RELEASE_THRESHOLD, UInt64(2^30))
@test attribute(UInt64, pool2, CUDA.MEMPOOL_ATTR_RELEASE_THRESHOLD) == 2^30

end

end

############################################################################################

@testset "profile" begin

@test_logs (:warn, r"only informs an external profiler to start") CUDA.Profile.start()
CUDA.Profile.stop()

@test_logs (:warn, r"only informs an external profiler to start") CUDA.@profile begin end

end

############################################################################################

@testset "stream" begin

s = CuStream()
synchronize(s)
@test CUDA.query(s) == true

let s2 = CuStream()
    @test s != s2
    @test !(s == s2)
end

let s3 = CuStream(; flags=CUDA.STREAM_NON_BLOCKING)
    @test s != s3
    @test !(s == s3)
end

prio = priority_range()
let s = CuStream(; priority=first(prio))
    @test priority(s) == first(prio)
end
let s = CuStream(; priority=last(prio))
    @test priority(s) == last(prio)
end

synchronize()
synchronize(; blocking=false)
synchronize(stream())
synchronize(stream(); blocking=false)

@grab_output show(stdout, stream())

end

############################################################################################

@testset "version" begin

@test isa(CUDA.version(), VersionNumber)

@test isa(CUDA.release(), VersionNumber)
@test CUDA.release().patch == 0

end


text:
if parse(Bool, get(ENV, "CI", "false"))
    @info "CI detected: skipping dataset download"
else
    data_dir = datadep"EMNIST"

    n_features = (28, 28)
    n_targets = 1

    for (name, n_obs) in [(:balanced,  (112800, 18800)),
                          (:byclass,  (697932, 116323)),
                          (:bymerge,  (697932, 116323)),
                          (:digits,  (240000, 40000)),
                          (:letters,  (124800, 20800)),
                          (:mnist,  (60000, 10000))]

        d = EMNIST(name)

        @test d.name == name
        @test d.split == :train
        @test extrema(d.features) == (0, 1)
        @test convert2image(d, 1) isa AbstractMatrix{<:Gray}
        @test convert2image(d, 1:2) isa AbstractArray{<:Gray, 3}

        test_supervised_array_dataset(d;
            n_features, n_targets, n_obs=n_obs[1],
            Tx=Float32, Ty=Int,
            conv2img=true)

        d = EMNIST(name, split=:test, Tx=UInt8)

        @test d.name == name
        @test d.split == :test
        @test extrema(d.features) == (0, 255)
        @test convert2image(d, 1) isa AbstractMatrix{<:Gray}

        test_supervised_array_dataset(d;
            n_features, n_targets, n_obs=n_obs[2],
            Tx=UInt8, Ty=Int,
            conv2img=true)
    end
end


text:
#=
# Import a simple triangle or quad mesh from a .msh file.
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("unstruct2dtest");

useLog("unstruct2dlog", level=3)

domain(2, grid=UNSTRUCTURED)
functionSpace(order=2)

# This rectangle covers [0, 0.1]x[0, 0.3]
# Uncomment the desired mesh.
mesh("utriangle.msh")  # Using triangles
#mesh("uquad.msh")     # Using quads

u = variable("u")
testSymbol("v")

boundary(u, 1, DIRICHLET, 0)

coefficient("f", "(-10-(x+1)*200*pi*pi)*sin(10*pi*x)*sin(10*pi*y) + 10*pi*cos(10*pi*x)*sin(10*pi*y)")
coefficient("k", "x+1")
coefficient("C", 10)
weakForm(u, "k*dot(grad(u), grad(v)) + C*u*v+ f*v")

solve(u);

finalize_finch();

# exact solution is sin(10*pi*x)*sin(10*pi*y)
# check error
maxerr = 0;
exact(x,y) = sin(10*pi*x)*sin(10*pi*y);

for i=1:size(Finch.grid_data.allnodes,2)
    x = Finch.grid_data.allnodes[1,i];
    y = Finch.grid_data.allnodes[2,i];
    err = abs(u.values[i] - exact(x,y));
    global maxerr;
    maxerr = max(err,maxerr);
end
println("max error = "*string(maxerr));

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
function read_input(file_path::String)
	file_stream::IOStream = open(file_path, "r")
	file_string::String = readstring(file_stream)
	close(file_stream)

	solid_structure_path::String = read_key(file_string, "solid_structure_path", String)
	diffusion_map_path::String = read_key(file_string, "diffusion_map_path", String)
	D0::Float64 = read_key(file_string, "D0", Float64)
	D1::Float64 = read_key(file_string, "D1", Float64)
	deltat_coarse::Float64 = read_key(file_string, "deltat_coarse", Float64)
	number_of_time_points_coarse::Int64 = read_key(file_string, "number_of_time_points_coarse", Int64)
	number_of_time_points_fine_per_coarse::Int64 = read_key(file_string, "number_of_time_points_fine_per_coarse", Int64)
	number_of_diffusers::Int64 = read_key(file_string, "number_of_diffusers", Int64)
	number_of_cells_x::Int64 = read_key(file_string, "number_of_cells_x", Int64)
	number_of_cells_y::Int64 = read_key(file_string, "number_of_cells_y", Int64)
	number_of_cells_z::Int64 = read_key(file_string, "number_of_cells_z", Int64)
	boundary_condition::String = read_key(file_string, "boundary_condition", String)
	output_file_path::String = read_key(file_string, "output_file_path", String)

	return (
		solid_structure_path,
		diffusion_map_path,
		D0,
		D1,
		deltat_coarse,
		number_of_time_points_coarse,
		number_of_time_points_fine_per_coarse,
		number_of_diffusers,
		number_of_cells_x,
		number_of_cells_y,
		number_of_cells_z,
		boundary_condition,
		output_file_path)
end


text:
# Interface

const Onepb = NamedTuple
const Prodb = AbstractVector{<: NamedTuple}

"""
`AbstractBasisSelector` : object specifying how a finite basis is selected from
the infinite space of symmetric polynomials. This type is pretty superfluous
for now since all basis selector and basis selection algorithms we admit so
far require Downsets.
"""
abstract type AbstractBasisSelector end

"""
`DownsetBasisSelector` : must implemented a non-negative valued `level` function
that is monotonically increasing with respect to the lexiographic ordering of
the basis functions.
"""
abstract type DownsetBasisSelector <: AbstractBasisSelector end

"""
return maximum correlation order of the basis selector. Note that the filter
function could reduce that by a bit.
"""
function maxorder end

"""
* `level(b::Onepb, Bsel::DownsetBasisSelector, basis::OneParticleBasis)`
* `level(bb::Prodb, Bsel::DownsetBasisSelector, basis::OneParticleBasis)`

The first version specifies the level function for
"""
function level end

"""
`maxlevel(Bsel, basis1p)`
"""
function maxlevel end

"""
`level1(b::Onepb, Bsel::DownsetBasisSelector, basis::OneParticleBasis)`

A specialized version of `level` to construct the 1-particle basis (cf.
`init1pspec!`). Fallback is to just call the `level` function, but this gives
some additional flexibility to ensure that the conditions of the `level`
framework are all justified, see docs for more detail.
"""
level1(b::Onepb, Bsel::DownsetBasisSelector, basis::OneParticleBasis) =
      level(b, Bsel, basis)

maxlevel1(Bsel::DownsetBasisSelector, basis::OneParticleBasis) =
      maxlevel(Bsel, basis)

"""
`filter(b_or_bb, Bsel::AbstractBasisSelector, basis::OneParticleBasis)`

After a down-set basis has been constructed, it can still be filtered, which
allows us to construct basis sets that aren't downsets but not too far away from
downsets. The main application is to enfore the constraints on m and l channels
arising from the symmetries.

Fallback implementation always returns `true`.
"""
filter(b_or_bb, Bsel::AbstractBasisSelector, basis::OneParticleBasis) = true

"""
No constraints on the basis - this selects that largest possible basis
subject to additional constraints baked into the one-particle basis.
In practise this should be equivalent to a naive max-norm basis selection,
and likely never used in practise.
"""
struct MaxBasis <: DownsetBasisSelector
   maxorder::Int
end

maxorder(Bsel::MaxBasis) = Bsel.maxorder

level(b::Onepb, Bsel::MaxBasis, basis::OneParticleBasis) =
      degree(b, basis)

level(bb::Prodb, Bsel::MaxBasis, basis::OneParticleBasis) =
      length(bb) == 0 ? 0 : sum(b -> level(b, Bsel, basis), bb)

maxlevel(Bsel::MaxBasis, basis::OneParticleBasis) = Inf

"""
`struct SimpleSparseBasis`:

The most basic form of a sparse basis selection, using the total degree.
Only the maximum correlation order and maximum degree may be specified.
This should primarily be used for testing. Construct it via
```julia
SimpleSparseBasis(maxorder, maxdegree)
```
Note that `maxdegree` is a basic total degree.
"""
struct SimpleSparseBasis <: DownsetBasisSelector
   maxorder::Int
   maxlevel::Float64
end

level(b::Onepb, Bsel::SimpleSparseBasis, basis::OneParticleBasis) =
      degree(b, basis)

level(bb::Prodb, Bsel::SimpleSparseBasis, basis::OneParticleBasis) =
      length(bb) == 0 ? 0 : sum( degree(b, basis) for b in bb )

maxlevel(bb, Bsel::SimpleSparseBasis, args...) =
      Bsel.maxlevel

maxlevel(Bsel::SimpleSparseBasis, args...) =
      Bsel.maxlevel

maxorder(Bsel::SimpleSparseBasis, args...) =
      Bsel.maxorder

"""
`AbstractSparseBasis`: Super-type for sparse basis selection as sub-levelsets of
the levelset function `level` and corresponding (possibly order-dependent)
levels provided in the dictionary `maxlevels::Dict{Any, Float64}`. In the default
implementation the levelset function and the degree function are identical.

Basis functions are selected in two steps. First, "admissible" basis
specifications are generated as a sub-levelset of the leveset function using the
implementation of the function `gensparse`.  After that, basis functions that do
not satisfy the conditons implemented in the function `filter` are removed from
the basis set.
"""
abstract type AbstractSparseBasis <: DownsetBasisSelector end

maxorder(Bsel::AbstractSparseBasis) = Bsel.maxorder

level(b::Onepb, Bsel::AbstractSparseBasis, basis::OneParticleBasis) =
      degree(b, basis, Bsel.weight)

level(bb::Prodb, Bsel::AbstractSparseBasis, basis::OneParticleBasis)  =  (
      length(bb) == 0 ? 0.0
                      : norm( level.(bb, Ref(Bsel), Ref(basis)), Bsel.p ) )

maxlevel(bb::Prodb, Bsel::AbstractSparseBasis, basis::OneParticleBasis) =
      maxlevel(length(bb), Bsel, basis)

maxlevel(ord::Integer, Bsel::AbstractSparseBasis, basis::OneParticleBasis) =
      ( haskey(Bsel.maxlevels, ord) ? Bsel.maxlevels[ord]
                                    : Bsel.maxlevels["default"] )

maxlevel(Bsel::AbstractSparseBasis, basis::OneParticleBasis) =
      maximum( maxlevel(ord, Bsel, basis) for ord = 1:maxorder(Bsel) )

"""
`SparseBasis`: basic implementation of an `AbstractSparseBasis`.
- `maxorder` : maximum correlation order
- `p` : degree (level) is computed via a weighted p-norm of the 1p basis
      function degrees
- `maxlevels` : for each order `ord`, `maxlevels[ord]` is the maximum level
      for basis functions of that order. if `maxlevels[ord]` is not
      specified, then `maxlevels["default"]` is the fallback.
- `weight` : specify weights for the different channels of the 1p basis
      functions, e.g., for an `Rn * Ylm` basis, one might specify
      `weight[:n] = 1` and `weight[:l] = 2` to have more radial and
      fewer angular basis functions.
"""
struct SparseBasis <: AbstractSparseBasis
   maxorder::Int
   weight::Dict{Symbol, Float64}
   maxlevels::Dict{Any, Float64}
   p::Float64
end

function SparseBasis(;  maxorder::Integer = nothing,
                        p = 1,
                        weight = Dict(:l => 1.0, :n => 1.0),
                        default_maxdeg = nothing,
                        maxlevels = nothing )
   if (default_maxdeg != nothing) && (maxlevels == nothing )
      return SparseBasis(maxorder, weight,
                         Dict{Any, Float64}("default" => default_maxdeg),
                         p)
   elseif (default_maxdeg == nothing) && (maxlevels != nothing)
      SparseBasis(maxorder, weight, maxlevels, p)
   else
      error("""Either both or neither optional arguments `maxlevels` and
               `default_maxdeg` were provided. To avoid ambiguity ensure that
               exactly one of these arguments is provided.""")
   end
end

"""
`CategorySparseBasis`: extension of `SparseBasis` that implements a
constraint on the correlation orders for within-category correlations, i.e., for
each category `cat` contained in OneParticleBasis associated with the symbol
`isym`, it is required that the within-category correlation order `corr[cat]`
satisfies
```julia
minorder_dict[cat] <= corr[cat] <= maxorder_dict[cat].
```
"""
struct CategorySparseBasis <: AbstractSparseBasis
   maxorder::Int
   weight::Dict{Symbol, Float64}
   maxlevels::Dict{Any, Float64}
   p::Float64
   isym::Symbol
   minorder_dict::Dict{Any, Int}
   maxorder_dict::Dict{Any, Int}
   weight_cat::Dict{Any, Float64}
end

function CategorySparseBasis(isym::Symbol, categories;
                             maxorder::Integer = nothing,
                             p = 1,
                             weight = Dict{Symbol, Float64}(),
                             default_maxdeg = nothing,
                             maxlevels::Dict{Any, Float64} = nothing,
                             minorder_dict = Dict{Any, Float64}(),
                             maxorder_dict = Dict{Any, Float64}(),
                             weight_cat = Dict(c => 1.0 for c in categories),
                             )
   if (default_maxdeg != nothing) && (maxlevels == nothing )
      CategorySparseBasis(maxorder, weight,
                          Dict{Any, Float64}("default" => default_maxdeg),
                          p, isym, minorder_dict, maxorder_dict, weight_cat)
   elseif (default_maxdeg == nothing) && (maxlevels != nothing)
      CategorySparseBasis(maxorder, weight, maxlevels, p, isym, minorder_dict,
                          maxorder_dict, weight_cat)
   else
      @error """Either both or neither optional arguments `maxlevels` and
                `default_maxdeg` were provided. To avoid ambiguity ensure that
                exactly one of these arguments is provided."""
   end
end

maxorder(Bsel::CategorySparseBasis, category) = Bsel.maxorder_dict[category]

minorder(Bsel::CategorySparseBasis, category) = Bsel.minorder_dict[category]

function filter(bb, Bsel::CategorySparseBasis, basis::OneParticleBasis)
   # auxiliary function to count the number of 1pbasis functions in bb
   # for which b.isym == s.
   num_b_is_(s) = sum([(getproperty(b, Bsel.isym) == s) for b in bb])

   # Within category min correlation order constaint:
   cond_ord_cats_min = all( num_b_is_(s) >= minorder(Bsel, s)
                            for s in keys(Bsel.minorder_dict) )
   # Within category max correlation order constaint:
   cond_ord_cats_max = all( num_b_is_(s) <= maxorder(Bsel, s)
                            for s in keys(Bsel.minorder_dict) )

   return cond_ord_cats_min && cond_ord_cats_max
end

# maxorder and maxlevel are inherited from the abstract interface

level(b::Onepb, Bsel::CategorySparseBasis, basis::OneParticleBasis) =
      cat_weighted_degree(b, Bsel, basis)

level(bb::Prodb, Bsel::CategorySparseBasis, basis::OneParticleBasis)  =
      cat_weighted_degree(bb, Bsel, basis)

# Category-weighted degree function
cat_weighted_degree(b::Onepb, Bsel::CategorySparseBasis, basis::OneParticleBasis) =
      degree(b, basis, Bsel.weight) * Bsel.weight_cat[getproperty(b, Bsel.isym)]

cat_weighted_degree(bb::Prodb, Bsel::CategorySparseBasis, basis::OneParticleBasis) = (
      length(bb) == 0 ? 0.0
                      : norm(cat_weighted_degree.(bb, Ref(Bsel), Ref(basis)), Bsel.p)
      )

# ---------------------------
# Some useful filters

struct NoConstant
end

(::NoConstant)(bb) = (length(bb) > 0)

"""
`EvenL`: selects all basis functions where the sum `L = sum_i l_i` of the degrees `l_i` of the spherical harmonics is even.
"""
struct EvenL
      isym::Symbol
      categories
end

function (f::ACE.EvenL)(bb)
      if isempty(bb)
            return true
      else
            suml(s) = sum( [getl(O3(), b) for b in bb if getproperty(b, f.isym) == s])
            return all(iseven(suml(s)) for s in f.categories)
      end
end

#=
"""
`DownsetIntersection`: Basis selector whose set of admissible specifications is the intersection
of the sets of admissible specifications of the sparse basis selectors contained in the lists `DBsels` and `ABsels`.
"""
struct DownsetIntersection <: DownsetBasisSelector
   DBsels::Vector{DownsetBasisSelector}
   ABsels::Vector{AbstractBasisSelector}
   maxorder::Int
end

maxlevel(Bsel::DownsetIntersection, basis::OneParticleBasis) = 1.0
maxorder(Bsel::DownsetIntersection) = Bsel.maxorder

function Base.intersect(Bsel1::DownsetIntersection,Bsel2::DownsetBasisSelector)
      return DownsetIntersection(vcat(Bsel1.DBsels,[Bsel2]), BSel1.ABsels, minimum([maxorder(Bsel1),maxorder(Bsel2)]))
end

function Base.intersect(Bsel1::DownsetIntersection,Bsel2::AbstractBasisSelector)
      return DownsetIntersection(Bsel1.DBsels, vcat(Bsel1.ABsels,[Bsel2]), minimum([maxorder(Bsel1),maxorder(Bsel2)]))
end

function Base.intersect(Bsel1::DownsetBasisSelector, Bsel2::AbstractBasisSelector)
      Bsel = DownsetIntersection(Vector{DownsetBasisSelector}([]), Vector{AbstractBasisSelector}([]), Inf)
      return intersect(intersect(Bsel, Bsel1), Bsel2)
end

Base.intersect(Bsel1::AbstractBasisSelector, Bsel2::DownsetBasisSelector)  = intersect(Bsel2,Bsel1)
Base.intersect(Bsel1::AbstractBasisSelector, Bsel2::DownsetIntersection) = intersect(Bsel2,Bsel1)

function Base.intersect(Bsel1::DownsetIntersection, Bsel2::DownsetIntersection)
      Bsel = deepcopy(Bsel1)
      for b in Bsel2
            Bsel = intersect(Bsel, b)
      end
      return Bsel
end

function level(bb, Bsel::DownsetIntersection, basis::OneParticleBasis)
   return maximum([ level(bb, bsel, basis)/maxlevel(bsel,basis) for bsel in Bsel.DBsels])
   #return maximum([ level(bb, bsel, basis)/maxlevel(length(bb),bsel,basis) for bsel in Bsel.DBsels])
end

function filter(bb, Bsel::DownsetIntersection, basis::OneParticleBasis)
   return all([filter(bb, bsel, basis) for bsel in Bsel.DBsels]) && all([filter(bb, bsel, basis) for bsel in Bsel.ABsels])
end
=#


text:
"""
    BnBSolver

Stores solver specific functions used to solve BnB problem. Has the following fields:
* `Lower_Prob::Any`:        Stores lower problem function (default = [])
* `Upper_Prob::Any`:        Stores upper problem function (default = [])
* `Preprocess::Any`:        Stores preprocessing function (default = [])
* `Postprocess::Any`:       Stores postprocessing function (default = [])
* `Term_Check::Any`:        Stores termination check function (default = 'Term_Check')
* `Branch_Sto::Any`:        Stores branching function (default = 'BM_depth_best!')
* `Node_Select::Any`:       Stores node selection function (default = 'NS_best')
* `Bisect_Func::Any`:       Stores branching function (default = 'Bisect_Rel')
* `Verbosity::String`:      Stores output selection (default = "Normal")
* `max_iter::Number`:       max number of iterations (default = "Inf")
* `iter_lim::Bool`:         determines if iteration limit is checked (default = false)
* `max_nodes::Int64`:       max number of nodes to store in memory (default = 1E6)
* `BnB_atol::Float64`:      absolute tolerance for BnB (default = 1E-4)
* `BnB_rtol::Float64`:      relative tolerance for BnB (default = 1E-4)
* `itr_intv::Int64`:        number of iterations to skip between printing iteration summary (default = 20)
* `hdr_intv::Int64`:        number of iterations to skip between printing header (default = 1)
* `converged::Any`:         convergence criterion (default = Conv_Check)
* `BnB_digits::Int64`:      digits displayed before decimal (default = 3)
* `hist_return::Bool`:      returns LBD, UBD array and time vector (default = false)
* `opt::Any`:               optional storage array (default = [])
* `exhaust::Bool`:          exhaustive search? (default = false)
* `target_upper::Float64`:  required upper bound (default = -Inf)
"""
mutable struct BnBSolver
  Lower_Prob::Any    # Stores lower problem function
  Upper_Prob::Any    # Stores upper problem function
  Preprocess::Any    # Preprocessing function
  Postprocess::Any   # Post processing function
  Repeat_Chk::Any
  Term_Check::Any    # Stores termination check function
  Branch_Sto::Any    # Stores branching function
  Node_Select::Any   # Stores node selection function
  Bisect_Func::Any   # Stores branching function
  Verbosity::String  # Stores output selection
  max_iter::Number   # max number of iterations
  iter_lim::Bool     # determines if iteration limit is checked
  max_nodes::Int64   # max number of nodes to store in memory
  BnB_atol::Float64   # absolute tolerance for BnB
  BnB_rtol::Float64  # relative tolerance for BnB
  itr_intv::Int64    # number of iterations to skip between printing iteration summary
  hdr_intv::Int64    # number of iterations to skip between printing header
  converged::Any     # convergence criterion
  BnB_digits::Int64  # digits displayed before decimal
  hist_return::Bool  # returns LBD, UBD array and time vector
  opt::Any           # optional storage array
  exhaust::Bool      # exhaustive search: find all solns or find first
  target_upper::Float64
end

"""
    BnBSolver()

Initializes solver with default parameters: best-first search, relative-width
bisection, no iteration limit, 1E6 node limit, 1E-4 absolute and relative
tolerances, no target upper bound for termination. No pre, or post
processing nodes and no repetition.
"""
BnBSolver() =  BnBSolver([],
                         [],
                         default_pre,
                         default_post,
                         Repeat_Node_Default,
                         Term_Check,
                         BM_depth_best!,
                         NS_best,
                         Bisect_Rel,
                         "Normal",
                         Inf,
                         false,
                         1E6,
                         1E-4,
                         1E-4,
                         1,
                         20,
                         Conv_Check,
                         3,
                         false,
                         [],
                         false,
                         -Inf)


text:
## Interoperability with Distributions.jl random variables
## ========================================================

## Conversion between Sigma random variaables and Distribution.jl random variables
## ===============================================================================
# Generates conversion functions of form:
# convert(::Type{Distributions.Poisson}, X::PoissonRandVar) =
#   Distributions.Poisson(X.λ.val)

# get_params(X::ElementaryRandVar, param_names) = [X.(p).val for p in param_names]

# # Generate functions
# for erv_type in subtypes(ElementaryRandVar)
#   # param 1 is dimension, ignore
#   param_names = erv_type.names[2:end]
#   dist_name = interop_name(erv_type)

#   p = Expr(:., :Distributions, QuoteNode(dist_name))
#   eval(:(convert(::Type{Distributions.$dist_name}, X::$(erv_type.name.name)) = $p(get_params(X, $param_names)...)))
# end

concretize(args::Tuple{Vararg{ConstantRandVar}}) = [arg.val for arg in args]
convert{T <: Distributions.Distribution}(D::Type{T}, X::ElementaryRandVar) =
  T(concretize(args(X))...)

interop_name{T <: ElementaryRandVar}(t::Type{T}) = Symbol(string(t.name.name)[1:end-7])
distribution_type{T <: ElementaryRandVar}(t::Type{T}) = getfield(Distributions, interop_name(T))

## Elementary Random Variables with Constant Parameters should be compatible
## with all of Distributios.jl fnuctions
## =========================================================================

## Distributions.jl interface
import Distributions: params, succprob, failprob, scale, location, shape,
  rate, ncategories, ntrials, dof

import Distributions: mean, var, std, median, modes, mode, skewness,
  kurtosis, isplatykurtic, isleptokurtic, ismesokurtic, entropy, entropy,
  mgf, cf

# Probability Evaluation
import Distributions: insupport, pdf, logpdf, loglikelihood, logcdf, ccdf, logccdf,
  cquantile, invlogcdf, invlogccdf

params(X::ElementaryRandVar) = params(convert(Distributions.Distribution, X))
succprob(X::ElementaryRandVar) = succprob(convert(Distributions.Distribution, X))
failprob(X::ElementaryRandVar) = failprob(convert(Distributions.Distribution, X))
scale(X::ElementaryRandVar) = scale(convert(Distributions.Distribution, X))
location(X::ElementaryRandVar) = location(convert(Distributions.Distribution, X))
shape(X::ElementaryRandVar) = shape(convert(Distributions.Distribution, X))
rate(X::ElementaryRandVar) = rate(convert(Distributions.Distribution, X))
ncategories(X::ElementaryRandVar) = ncategories(convert(Distributions.Distribution, X))
ntrials(X::ElementaryRandVar) = ntrials(convert(Distributions.Distribution, X))


text:
# ---
# title: 144. Binary Tree Preorder Traversal
# id: problem144
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Stack, Tree
# link: <https://leetcode.com/problems/binary-tree-preorder-traversal/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, return _the preorder traversal of its nodes
# ' values_.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/15/inorder_1.jpg)
#
#
#
#     Input: root = [1,null,2,3]
#     Output: [1,2,3]
#
#
# **Example 2:**
#
#
#
#     Input: root = []
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: root = [1]
#     Output: [1]
#
#
# **Example 4:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/15/inorder_5.jpg)
#
#
#
#     Input: root = [1,2]
#     Output: [1,2]
#
#
# **Example 5:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/15/inorder_4.jpg)
#
#
#
#     Input: root = [1,null,2]
#     Output: [1,2]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[0, 100]`.
#   * `-100 <= Node.val <= 100`
#
#
#
# **Follow up:**
#
# Recursive solution is trivial, could you do it iteratively?
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Documenter, GrowthMaps, Weave, IJulia

basedir = dirname(@__FILE__)

jmdpath = joinpath(basedir, "src/example.jmd")

# Generate examples markdown and images
mdpath = joinpath(basedir, "src/example.md")
weave(jmdpath, out_path=mdpath, doctype="github")

# Generate examples notebook
notebookdir = joinpath(basedir, "build/notebook")
notebookpath = joinpath(notebookdir, "example.ipynb")
mkpath(notebookdir)
convert_doc(jmdpath, notebookpath)

# Generate HTML docs
makedocs(
    modules = [GrowthMaps],
    sitename = "GrowthMaps.jl",
    pages = [
        "Home" => "index.md",
        "Examples" => "example.md",
    ],
    clean = false,
)

deploydocs(
    repo = "github.com/cesaraustralia/GrowthMaps.jl.git",
)


text:
include("./../src/Lar.jl")
using SparseArrays, ViewerGL, LinearAlgebraicRepresentation
LAR = LinearAlgebraicRepresentation
LarA = Lar.Arrangement
GL = ViewerGL

#=== LAURENTINA DATASET ===#
#  pre-crop E, V = 24043×22808
#  cropped  E, V = 21454×20389
#  arranged E, V =
filename = "examples/osm/laurentina.osm"
filename = "examples/osm/laurentina2.osm"

#=== MARCONI DATASET ===#
#  pre-crop V =     ; E =     ; F =   0
#  cropped  V = 5153; E = 5489; F =   0
#  arranged V = 4600; E = 5062; F = 838
#filename = "examples/osm/marconi.osm"

#=== ALBANO LAZIALE DATASET ===#
#  pre-crop E, V =
#  cropped  E, V = 9205×8879
#  arranged E, V = 9305×8892
#filename = "examples/osm/albano.osm"

model = Lar.osm2lar(filename, true, false)

GL.VIEW([GL.GLGrid(model.G, LAR.cop2lar(model.T[1]), GL.COLORS[3])])

model, edge_map = LarA.planar_arrangement(model, spzeros(Int8, 0), true, false)

GL.VIEW([GL.GLGrid(model.G, LAR.cop2lar(model.T[1]), GL.COLORS[3])])

FVs = Lar.triangulate2D(model)
GL.VIEW(GL.GLExplode(model.G,FVs, 1., 1., 1.,99,1));
GL.VIEW(GL.GLExplode(model.G,FVs,1.2,1.2,1.2,99,1));

GL.VIEW( GL.GLExplode(
    model.G,
    [f for f in FVs if length(f) > 20],
    1.5, 1.5, 1.5, 99, 1)
);

#==

#=== NOW IT LOOPS ===#
LarA = Lar.Arrangement

sigma           = SparseArrays.spzeros(Int8, 0)
return_edge_map = false
multiproc       = false

# model, edge_map = LarA.planar_arrangement_1(model, sigma, true, multiproc)

edgenum = size(model, 1, 1)
edge_map = Array{Array{Int, 1}, 1}(undef,edgenum)
rV = Lar.Points(zeros(0, 2))
rEV = SparseArrays.spzeros(Int8, 0, 0)
finalcells_num = 0

# spaceindex computation
bigPI = LarA.spaceIndex(model)
...

bicon_comps = LarA.biconnected_components(model.T[1])
V           = convert(Lar.Points, model.G')
copEV       = model.T[1]

# arrangement of isolated components
n = size(bicon_comps, 1)
shells = Array{Lar.Chain, 1}(undef, n)
boundaries = Array{Lar.ChainOp, 1}(undef, n)
EVs = Array{Lar.ChainOp, 1}(undef, n)
# for each component
for p in 2 : n
    ev = copEV[sort(bicon_comps[p]), :]
    # computation of 2-cells
    fe = LarA.minimal_2cycles(V, ev)
    # exterior cycle
    shell_num = LarA.get_external_cycle(V, ev, fe)
    # decompose each fe (co-boundary local to component)
    EVs[p] = ev
    tokeep = setdiff(1:fe.m, shell_num)
    boundaries[p] = fe[tokeep, :]
    shells[p] = fe[shell_num, :]
end

p = 1
ev = copEV[sort(bicon_comps[p]), :]
# computation of 2-cells
fe = LarA.minimal_2cycles(V, ev)                         # <---| Here it loops
# exterior cycle
shell_num = LarA.get_external_cycle(V, ev, fe)
# decompose each fe (co-boundary local to component)
EVs[p] = ev
tokeep = setdiff(1:fe.m, shell_num)
boundaries[p] = fe[tokeep, :]
shells[p] = fe[shell_num, :]

==#


text:
################################################################################
# Architecture
################################################################################

# We want to be able to control the selection of the color for each tile to
# allow different patterns to be generated easily.
#
# Control this assignment with an ArchitectureColor that will control dispatch
# to methods throughout architecture construction.

"""
    abstract type ArchitectureColor end

Abstract supertype for color assignment rules.

API:

    `pick(::ArchitectureColor, addr::Address) :: SquareColor`

Allow selection of the color for the tile at `addr`.
"""
abstract type ArchitectureColor end

# Color squares in a chessboard pattern.
struct ChessboardColor <: ArchitectureColor end
# Determine color based on the sum of the components of the address.
pick(::ChessboardColor, addr) = isodd(sum(addr.I)) ? White : Black

struct HashColor <: ArchitectureColor end
pick(::HashColor, addr) = isodd(hash(addr)) ? White : Black

abstract type ArchitectureStyle end

struct Rectangle2D <: ArchitectureStyle end
portnames(::Rectangle2D) = ("north", "east", "south", "west")
nd(::Rectangle2D) = 2
function rules(::Rectangle2D)
    return tuple(
        ConnectionRule([
            Offset((-1, 0), "north_out", "south_in"),
            Offset((1, 0), "south_out", "north_in"),
            Offset((0, 1), "east_out", "west_in"),
            Offset((0, -1), "west_out", "east_in"),
        ]),
    )
end

struct Rectangle3D <: ArchitectureStyle end
portnames(::Rectangle3D) = ("north", "east", "south", "west", "up", "down")
nd(::Rectangle3D) = 3
function rules(::Rectangle3D)
    return tuple(
        ConnectionRule([
            Offset((-1, 0, 0), "north_out", "south_in"),
            Offset((1, 0, 0), "south_out", "north_in"),
            Offset((0, 1, 0), "east_out", "west_in"),
            Offset((0, -1, 0), "west_out", "east_in"),
            Offset((0, 0, 1), "up_out", "down_in"),
            Offset((0, 0, -1), "down_out", "up_in"),
        ]),
    )
end

#=
          ____          ____
         /    \        /    \
        /      \      /      \
   ____/  (0,1) \____/ (0,3)  \
  /    \        /    \        /
 /      \      /      \      /
/ (0,0)  \____/ (0,2)  \____/
\        /    \        /    \
 \      /      \      /      \
  \____/ (1,1)  \____/ (1,3)  \
  /    \        /    \        /
 /      \      /      \      /
/ (1,0)  \____/ (1,2)  \____/
\        /    \        /
 \      /      \      /
  \____/        \____/
=#

struct Hexagonal2D <: ArchitectureStyle end
portnames(::Hexagonal2D) = ("30", "90", "150", "210", "270", "330")
nd(::Hexagonal2D) = 2
function rules(::Hexagonal2D)
    return (
        # Rule to apply if the column address is even.
        ConnectionRule(
            [
                Offset((0, 1), "30_out", "210_in"),
                Offset((-1, 0), "90_out", "270_in"),
                Offset((0, -1), "150_out", "330_in"),
                Offset((1, -1), "210_out", "30_in"),
                Offset((1, 0), "270_out", "90_in"),
                Offset((1, 1), "330_out", "150_in"),
            ];
            address_filter = x -> iseven(x.I[2]),
        ),
        # Rule to apply if the column address is odd.
        ConnectionRule(
            [
                Offset((-1, 1), "30_out", "210_in"),
                Offset((-1, 0), "90_out", "270_in"),
                Offset((-1, -1), "150_out", "330_in"),
                Offset((0, -1), "210_out", "30_in"),
                Offset((1, 0), "270_out", "90_in"),
                Offset((0, 1), "330_out", "150_in"),
            ];
            address_filter = x -> isodd(x.I[2]),
        ),
    )
end

# Build mappable component
function mappable(color::SquareColor)
    component = Component("square"; metadata = Dict("color" => color))
    add_port(component, "in", Input)
    add_port(component, "out", Output)

    return component
end

# Since Crossbars are Mapper primitives, we are now ready to build our tile.
function build_tile(color::SquareColor, style::ArchitectureStyle)
    tile = Component("tile")

    # Ports
    directions = portnames(style)
    for direction in directions
        add_port(tile, "$(direction)_in", Input)
        add_port(tile, "$(direction)_out", Output)
    end

    # Sub-components
    add_child(tile, mappable(color), "square")

    # Add 1 to the numer of IO directions to allow a connection to be made to
    # the mappable core.
    nports = length(directions) + 1
    add_child(tile, build_mux(nports, nports), "mux")

    # Links

    # Mux to Square
    add_link(tile, "square.out", "mux.in[0]")
    add_link(tile, "mux.out[0]", "square.in")

    # Connect IO based on the provided directions.
    for (idx, direction) in enumerate(directions)
        add_link(tile, "mux.out[$idx]", "$(direction)_out")
        add_link(tile, "$(direction)_in", "mux.in[$idx]")
    end

    return tile
end

function architecture(
    tiles_per_side, style::ArchitectureStyle, color_rule::ArchitectureColor
)

    # Get the number of dimensions for this architecture.
    N = nd(style)

    arch = TopLevel{N}("architecture")

    # Instantiate squares
    for address in CartesianIndices(Tuple(tiles_per_side for _ in 1:N))
        color = pick(color_rule, address)
        add_child(arch, build_tile(color, style), address)
    end

    # Connect cores.
    for rule in rules(style)
        connection_rule(arch, rule)
    end

    return arch
end


text:
module Thunks

include("core.jl")
include("macros.jl")

export Thunk, Unevaluated, Reversible, Think, Checkpointable, reify, thunk, @lazy,
    @noeval, @reversible, undo
end


text:
import HoTRG_lemon.SpinModule:
	getTrotterparameter,
	getTrotteriteration,
	getTrotterlayers,
	getBeta,
	setTrotterparameter!,
	setTrotterIteration!

abstract QuantumLattice{T} <: Lattice{T}

function buildQuantumLattice(latticeinfo::LatticeInfo, spinmodel::SpinModel)
	if isClassical(spinmodel)
		error("wrong lattice")
	end
	if isTwoDimension(latticeinfo) & isSquareLattice(latticeinfo)
		return Quantum2dSquareLattice(spinmodel)
	elseif isTwoDimension(latticeinfo) & isFractalLattice(latticeinfo)
		return Quantum2dFractalLattice(spinmodel)
	else
		error("wrong latticeinfo")
	end
end

#---
# functions inherit from QuantumSpinModel
get_functions = [
				:getTrotterparameter,
				:getTrotteriteration,
				:getTrotterlayers,
				:getBeta
				]
for f in get_functions
	eval(quote
	  function $(f){T}(lattice::QuantumLattice{T})
	  	$(f)(lattice.spinmodel)
		end
	end)
end

set_functions = [
				:setTrotterparameter!,
				:setTrotterIteration!
				]

for f in set_functions
	eval(quote

		  function $(f){T}(lattice::QuantumLattice{T}, parameter)
		  $(f)(lattice.spinmodel, parameter)
			initialize!(lattice)
		end

	end)
end


text:
# http://www.analytictech.com/networks/dataentry.htm
# http://www.analytictech.com/ucinet/help/hs5000.htm

function read_dl(io::IO)
    g = Graph(16)
    readline(io)
    readline(io)
    readline(io)
    readline(io)
    readline(io)
    readline(io)
    while !eof(io)
        l = readline(io)
        a, b, c = split(l, ' ')
        addedge!(g, parse(Int, a), parse(Int, b))
    end
    return g
end


text:
#=
Implementation of a metric.

Using the Colors module.

You might not be able to include this,
if you also define another metric having the same signature.

I initially used this code
because the Colors module exists
and I thought it would work since it does compute a "difference" in some sense.

But the Colors module has different concerns:
the vagaries of human perception.
We don't need that, instead we need fast.
Experimentally: too slow because of issues with allocations in a call to colordiff(),
because of issues with Julia precompiling failing to infer properly.
=#

using Colors

#=
More generally, call it a "value."
=#
function pointDifference(targetColor, corpusColor )
    result = 0

    result = colordiff(
        targetColor,
        corpusColor,
        metric )
    # WAS metric=DE_94() i.e. computing the metric deep in the loops
    return result
end

#=
The maximum value of the metric.
Approximate constants from the literature.
According to Wolfram (but not the Julia site)
the "typical" max for CIE2000 colordiff
1.16839  # CIE2000
1.48916  # CIE84 DE_94

!!! Depends on metric passed to colorDiff
=#

const maxPointDifference = 1.16839

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
module CUDA

using Reexport

@reexport using CUDAdrv

@eval $(Symbol("@elapsed")) = $(getfield(CUDAdrv, Symbol("@elapsed")))
@eval $(Symbol("@profile")) = $(getfield(CUDAdrv, Symbol("@profile")))

@reexport using CUDAnative

## math intrinsics
for intr in [:cos, :cos_fast, :cospi, :sin, :sin_fast, :sinpi, :tan, :tan_fast, :acos,
             :asin, :atan, :atan2, :angle, :cosh, :sinh, :tanh, :acosh, :asinh, :atanh,
             :log, :log_fast, :log, :log_fast, :log10, :log10_fast, :log1p, :log2,
             :log2_fast, :logb, :ilogb, :exp, :exp_fast, :exp2, :exp10, :exp10_fast, :expm1,
             :ldexp, :exp, :exp_fast, :erf, :erfinv, :erfc, :erfcinv, :erfcx, :brev, :clz,
             :ffs, :byte_perm, :popc, :isfinite, :isinf, :isnan, :nearbyint, :nextafter,
             :signbit, :copysign, :abs, :sqrt, :rsqrt, :cbrt, :rcbrt, :pow, :pow_fast, :pow,
             :abs2, :ceil, :floor, :min, :max, :saturate, :mod, :rem, :div_fast, :lgamma,
             :tgamma, :j0, :j1, :jn, :y0, :y1, :yn, :normcdf, :normcdfinv, :hypot, :fma,
             :sad, :dim, :mul24, :mul64hi, :mulhi, :hadd, :rhadd, :scalbn]
    @eval const $intr = $(getfield(CUDAnative, intr))
end

## atomics
const atomic_xchg!  = CUDAnative.atomic_xchg!
const atomic_add!   = CUDAnative.atomic_add!
const atomic_and!   = CUDAnative.atomic_and!
const atomic_or!    = CUDAnative.atomic_or!
const atomic_xor!   = CUDAnative.atomic_xor!
const atomic_min!   = CUDAnative.atomic_min!
const atomic_max!   = CUDAnative.atomic_max!
const atomic_inc!   = CUDAnative.atomic_inc!
const atomic_dec!   = CUDAnative.atomic_dec!

@reexport using CuArrays

## array constructors
const zeros = CuArrays.zeros
const ones  = CuArrays.ones
const fill  = CuArrays.fill

## random numbers
const fill          = CuArrays.fill
const seed!         = CuArrays.seed!
const rand          = CuArrays.rand
const randn         = CuArrays.randn
const rand_logn     = CuArrays.rand_logn
const rand_poisson  = CuArrays.rand_poisson

@eval $(Symbol("@sync"))        = $(getfield(CuArrays, Symbol("@sync")))
@eval $(Symbol("@time"))        = $(getfield(CuArrays, Symbol("@time")))
@eval $(Symbol("@allocated"))   = $(getfield(CuArrays, Symbol("@allocated")))

end


text:
export randompolicy, leftdownpolicy

"""
    Chooses a direction randomly
"""
function randompolicy(_)
    Float32[0.25, 0.25, 0.25, 0.25]
end

"""
    Chooses left or down, and the other direction rarely.
    Intended for use with the Greedy player
"""
function leftdownpolicy(_)
    Float32[0.99-2eps(Float32), eps(Float32), eps(Float32), 0.01]
end


text:
# using JuMP
# using LazySets
#TODO: You could use Polyhedra to get the feasible set (poly = polyhedron(m) m=Model())
#TODO: Naming of the constraints => in order to have better debugging

# Used for DenseAxisArray ========================================
# x ∈ ℝn×N+1 and 𝒳 ⊂ ℝⁿ
function add_constraint!(m::Model, x::AbstractArray{VariableRef,2}, 𝒳::AbstractPolyhedron)
    a, b, s = get_constraints(𝒳)
    N = size(x,2)-1 # Since the first element of the JumpVariable is 0
    for j=0:N
        @constraint(m, [i=1:s], a[i]'*x[:, j] <= b[i])
    end
end

# x ∈ ℝⁿ and 𝒳 ⊂ ℝⁿ
function add_constraint!(m::Model, x::AbstractVector{VariableRef}, 𝒳::AbstractPolyhedron)
    a, b, s = get_constraints(𝒳)
    @constraint(m, [i=1:s], a[i]'*x <= b[i])
end

# x ∈ ℝⁿ and 𝒳 = 𝓍 ∈ ℝⁿ
function add_constraint!(m::Model, x, 𝓍::AbstractVector{<:Real})
    n = length(𝓍)
    @constraint(m, [i=1:n], x[i] == 𝓍[i])
end

# Used for Array{VariableRef, 2} ========================================
# x ∈ ℝⁿ and 𝒳 = 𝓍 ∈ ℝⁿ
function add_constraint!(m::Model, x::Vector{VariableRef}, 𝓍::AbstractVector{<:Real})
    n = length(𝓍)
    @constraint(m, x .== 𝓍)
end

function add_constraint!(m::Model, x::Array{VariableRef,2}, 𝒳::AbstractPolyhedron)
    a, b, s = get_constraints(𝒳)
    N = size(x,2) # Since, here, the first element of the JumpVariable is 1
    for j=1:N
        @constraint(m, [i=1:s], a[i]'*x[:, j] <= b[i])
    end
end

# # x ∈ ℝⁿ and 𝒳 = 𝓍 ± ϵ ∈ ℝⁿ (relaxed terminal constraint)
# function add_constraint!(m::Model, x, 𝓍::AbstractVector, ϵ=1e-5)
#     n = length(𝓍)
#     @constraint(m, [i=1:n], x[i] <= 𝓍[i] + ϵ)
#     @constraint(m, [i=1:n], x[i] <= 𝓍[i] - ϵ)
# end

# # TODO: make it a macro! such that the constr name can be added to the constraint
# function constrain_variable_by_set!(model::Model, x::AbstractVector{VariableRef},
#                                     𝒳::Interval{Float64,LazySets.IntervalArithmetic.Interval{Float64}})
#     N = length(x)
#     a = [1.0, -1.0]
#     b = [𝒳.dat.hi, -𝒳.dat.lo]
#     for j=0:N-1
#         @constraint(model, [i=1:2], a[i]'*x[j] <= b[i])
#     end
# end
#
# # 1D variable x and 1D Interval set 𝒳
# function constrain_variable_by_set!(model::Model, x::VariableRef,
#                                     𝒳::Interval{Float64,LazySets.IntervalArithmetic.Interval{Float64}})
#     a = [1.0, -1.0]
#     b = [𝒳.dat.hi, -𝒳.dat.lo]
#     @constraint(model, [i=1:2], a[i]'*x <= b[i])
# end

""" Checks JuMP.MOI.TerminationStatusCode """
function nonlin_opt_okay(model::Model)
    term_status = termination_status(model)
    if Int(term_status) == 1
        # everything is fine
    elseif Int(term_status) == 4
        @warn("No global Optimum found: $term_status")
    elseif (Int(term_status) == 7 || Int(term_status) == 10)
        @warn("Solved to relaxed tolerances: $term_status")
    else
        error("Optimization not a local or global Optimum: $term_status")
    end
    nothing
end

function get_ineq_constr(model::Model)
    return all_constraints(model, GenericAffExpr{Float64, VariableRef}, MOI.LessThan{Float64})
end

function get_eq_constr(model::Model)
    return all_constraints(model, GenericAffExpr{Float64, VariableRef}, MOI.EqualTo{Float64})
end

# list_of_constraint_types(model)
# TODO: print name of constraints!
# print(model)
# Min f(l(z + -0.200000000000002, v), l(z + 0.200000000000002, v))
# Subject to
#  0.5 z - v == 0.0
#  z <= 4.799999999999998
#  -z <= 4.799999999999998
#  v <= 2.0
#  -v <= 2.0


text:
function string2cmd(str::String)

ind = findall(" ", str);
n = length(ind);
cmdstr = Array{String,1}(undef,n+1);
cmdstr[1] = str[1:(ind[1][1]-1)];
for i = 2:n
    cmdstr[i] = str[(ind[i-1][1]+1):(ind[i][1]-1)];
end
i = n+1;
cmdstr[i] = str[(ind[i-1][1]+1):end];

return Cmd(cmdstr)

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
module KUparser

using Compat, DistributedArrays

include("util.jl");	  export @date, evalparse #, testnet
include("types.jl");      export Sentence, Corpus, Vocab, WordId, DepRel, PosTag, Position, Cost, Move, SFtype, WVtype, Pvec, Dvec, wdim, wcnt, wtype
include("parser.jl");     export ArcEager13, ArcEagerR1, ArcHybrid13, ArcHybridR1, Parser, reset!
include("features.jl");   export Fvec, flen, xsize, ysize
include("flist.jl");      export Flist
include("rparser.jl");    export rparse
include("oparser.jl");    export oparse
include("gparser.jl");    export gparse
#include("bparser.jl");    export bparse, bparse_pmap

end # module

# include("resetworkers.jl");


text:
function bubblechart(; x::AbstractVector = Int[],
					   y::AbstractVector = Int[],
					   group::AbstractVector = [],
					   pointSize::AbstractVector = Int[],
					   pointShape::AbstractString = "circle")

    v = scatterplot(x = x, y = y, group = group)
    #Get name before destroying
    table = v.data[1].name

    empty!(v.data)
    add_data!(v, name = table, x = x, y = y, group = group, y2 = pointSize)

    v.marks[1].properties.enter.size = VegaValueRef(field = "y2", mult = 30)
    v.marks[1].properties.enter.shape = VegaValueRef(value = pointShape)

    return v

end

max_stars_repo_name: jpthor/vega