text:
const LOCATIONS = Dict( k => i-1 for (i, k) in enumerate((
    "none", "upper right", "upper left", "lower left", "lower right",
    "right", "center left", "center right", "lower center", "upper center", "center",
    "outer upper right", "outer center right", "outer lower right"
)))

# Legend
function legend!(p::PlotObject, args...; location=1, kinds=Tuple{}(), kwargs...)
    location = lookup(location, LOCATIONS)
    # Reset main viewport if there was a legend
    if haskey(p.attributes, :location) && p.attributes[:location] ∈ LEGEND_LOCATIONS[:right_out]
        p.viewport.inner[2] += p.legend.size[1]
    end
    chosen = choosegeoms(p, kinds)
    for i = 1:min(length(args), length(chosen))
        j = chosen[i]
        p.geoms[j] = Geometry(p.geoms[j], label=args[i])
    end
    maxrows = Int(get(kwargs, :maxrows, length(p.geoms)))
    p.legend = Legend(p.geoms, p.viewport.inner, maxrows)
    # Redefine viewport if legend is set outside
    if p.legend.size ≠ NULLPAIR && location ∈ LEGEND_LOCATIONS[:right_out]
        p.viewport.inner[2] -= p.legend.size[1]
    end
    p.attributes[:location] = location
end

legend!(f::Figure, args...; kwargs...) = legend!(currentplot(f), args...; kwargs...)

"""
    legend(labels...; kwargs...)

Set the legend of the plot, using a series of `labels` (strings).

In addition to the legend strings, the keyword argument
`location` can be used to define the location of the legend with
respect to the plot axes and the keyword argument `maxrows`
to distribute the legend labels in a grid with a maximum number of rows.

Locations are defined as a number or a string, as indicated
in the following table --- based on the convention of
[Matplotlib legends](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.legend.html):

|#  | String                |
|--:|:----------------------|
|  0| `"none"`              |
|  1| `"upper right"`       |
|  2| `"upper left"`        |
|  3| `"lower left"`        |
|  4| `"lower right"`       |
|  5| `"right"`             |
|  6| `"center left"`       |
|  7| `"center right"`      |
|  8| `"lower center"`      |
|  9| `"upper center"`      |
| 10| `"center"`            |
| 11| `"outer upper right"` |
| 12| `"outer center right"`|
| 13| `"outer lower right"` |

The labels are assigned to the geometries contained in the plot,
in the same order as they were created. The assignment can be
restricted to specific kinds of geometries through the keyword argument
`kinds`, which can take a `Symbol` or a collection of `Symbol`s
that identify the kinds. Use the helper function [`geometrykinds`](@ref)
to see the list of kinds available in the current plot.

Only geometries with non-empty labels and an available guide for legends
will be presented in the legend.

# Examples

```julia
# Set the legends to "a" and "b"
legend("a", "b")
```
"""
function legend(args::AbstractString...; kwargs...)
    f = gcf()
    legend!(currentplot(f), args...; kwargs...)
    return f
end

"""
    geometrykinds([p])

Return a list with symbols that represent the kind
of the geometries included in the given plot or figure `p`.

If no argument is given, it takes the current plot of
the current figure.

# Examples

```julia
julia> # Plot a set of points at values `(x, y)`
julia> # and a regression line passing through `(x, ŷ)`
julia> scatter(x, y)
julia> plot(x, ŷ)
julia> geometrykinds()
2-element Array{Symbol,1}:
 :scatter
 :line
```
"""
geometrykinds(p::PlotObject) = [g.kind for g in p.geoms]
geometrykinds(f::Figure=gcf()) = geometrykinds(currentplot(f))

function choosegeoms(p::PlotObject, kinds=Tuple{}())
    isempty(kinds) && return collect(1:length(p.geoms))
    gk = geometrykinds(p)
    findall(k -> k ∈ kinds, gk)
end

choosegeoms(p::PlotObject, kinds::Symbol) = choosegeoms(p, (kinds,))

# Hold
hold!(p::PlotObject, state::Bool) = (p.attributes[:hold] = state)

hold!(f::Figure, state) = hold!(currentplot(f), state)

"""
    hold(flag::Bool)

Set the hold flag for combining multiple plots.

`hold(true)` prevents clearing previous plots, so that next plots
will be drawn on top of the previous one until `hold(false)` is called.

Use the keyword argument `hold=<true/false>` in plotting functions, to
set the hold flag during the creation of plots.
"""
hold(state) = hold!(currentplot(gcf()), state)

# Title
function title!(p::PlotObject, s)
    if isempty(s)
        delete!(p.attributes, :title)
    else
        p.attributes[:title] = s
    end
    return nothing
end

title!(f::Figure, s) = title!(currentplot(f), s)

"""
    title(s)

Set the plot title as the string `s`.

Use the keyword argument `title=s` in plotting functions, to
set the title during the creation of plots.

# Examples

```julia
# Set the plot title to "Example Plot"
title("Example Plot")
# Clear the plot title
title("")
```
"""
function title(s::AbstractString)
    f = gcf()
    title!(currentplot(f), s)
    return f
end

const AXISLABEL_DOC = """
    xlabel(s)
    ylabel(s)
    zlabel(s)

Set the X, Y or Z axis labels as the string `s`.

Use the keyword argument `xlab=s`, etc. in plotting functions, to
set the axis labels during the creation of plots.

# Examples

```julia
# Set the x-axis label to "x"
xlabel("x")
# Clear the y-axis label
ylabel("")
```
"""

const TICKS_DOC = """
    xticks(minor[, major = 1])
    yticks(minor[, major = 1])
    zticks(minor[, major = 1])

Set the `minor`intervals of the ticks for the X, Y or Z axis,
and (optionally) the number of minor ticks between `major` ticks.

Use the keyword argument `xticks=(minor, major)`, etc. in plotting functions, to
set the tick intervals during the creation of plots (both the minor and major
values are required in this case).

# Examples

```julia
# Minor ticks every 0.2 units in the X axis
xticks(0.2)
# Major ticks every 1 unit (5 minor ticks) in the Y axis
yticks(0.2, 5)
```
"""

# Attributes for axes

const AXISLIM_DOC = """
    xlim(inf, sup [, adjust::Bool = false])
    xlim((inf, sup), ...)
    ylim(inf, sup, ...)
    ylim((inf, sup), ...)
    zlim(inf, sup, ...)
    zlim((inf, sup), ...)

Set the limits for the plot axes.

The axis limits can either be passed as individual arguments or as a
tuple of `(inf, sup)` values. Setting either limit to `nothing` will
cause it to be automatically determined based on the data, which is the
default behavior.

Additionally to the limits, the flag `adjust` can be used to
tell whether or not the limits have to be adjusted.

Use the keyword argument `xlim=(inf, sup)`, etc. in plotting functions, to
set the axis limits during the creation of plots.

# Examples

```julia
# Set the x-axis limits to -1 and 1
xlim((-1, 1))
# Reset the x-axis limits to be determined automatically
xlim()
# Set the y-axis upper limit and set the lower limit to 0
ylim((0, nothing))
# Reset the y-axis lower limit and set the upper limit to 1
ylim((nothing, 1))
```
"""

const AXISLOG_DOC = """
    xlog(flag)
    ylog(flag)
    zlog(flag)

Set the X-, Y- or Z-axis to be drawn in logarithmic scale (`flag == true`),
or in linear scale (`flag == false`).

Use the keyword argument `xlog=<true/false>`, etc. in plotting functions, to
set the logarithmic axes during the creation of plots.

!!! note

    When the axis is set to logarithmic scale, its lower limit is adjusted
    to represent only positive values, even if the data of the plot contain
    zero or negative values. The aspect of logarithmic axes with limits
    explicitly set to contain negative values (with [`xlim`](@ref), etc.)
    is undefined.

# Examples

```julia
# Set the x-axis limits to log scale
xlog(true)
# Ensure that the y-axis is in linear scale
ylog(false)
```
"""

const AXISFLIP_DOC = """
    xflip(flag)
    yflip(flag)
    zflip(flag)

Reverse the direction of the X-, Y- or Z-axis (`flag == true`),
or set them back to their normal direction (`flag == false` ).

Use the keyword argument `xflip=<true/false>`, etc. in plotting functions, to
set reversed axes during the creation of plots.

# Examples

```julia
# Reverse the x-axis
xflip(true)
# Ensure that the y-axis is not reversed
yflip(false)
```
"""

@eval function _config_axislimits!(ax, p, (minval, maxval), adjust)
    data_limits = minmax(p.geoms, p.axes.options[:scale])[ax]
    limits = set_limits((minval, maxval), data_limits)
    adjust && (limits = GR.adjustlimits(limits...))
    p.axes.ranges[ax] = limits
    tickdata = p.axes.tickdata
    if haskey(tickdata, ax)
        axisticks = tickdata[ax]
        if get(p.attributes, Symbol(ax,:flip), false)
            limits = reverse(limits)
        end
        tickdata[ax] = (axisticks[1], limits, axisticks[3])
    end
    return nothing
end

for ax = ("x", "y", "z")
    # xlabel, etc.
    fname! = Symbol(ax, :label!)
    fname = Symbol(ax, :label)
    @eval function $fname!(p::PlotObject, s)
        if isempty(s)
            delete!(p.attributes, Symbol($ax, :label))
        else
            p.attributes[Symbol($ax, :label)] = s
        end
        return nothing
    end
    @eval $fname!(f::Figure, s) = $fname!(currentplot(f), s)
    @eval @doc AXISLABEL_DOC function $fname(s::AbstractString)
        f = gcf()
        $fname!(currentplot(f), s)
        return f
    end

    # xticks, etc.
    fname! = Symbol(ax, :ticks!)
    fname = Symbol(ax, :ticks)
    @eval function $fname!(p::PlotObject, minor, major=1)
        tickdata = p.axes.tickdata
        if haskey(tickdata, Symbol($ax))
            tickdata[Symbol($ax)] = (float(minor), tickdata[Symbol($ax)][2], Int(major))
        end
        p.attributes[Symbol($ax, :ticks)] = (minor, major)
        return nothing
    end
    @eval $fname!(f::Figure, args...) = $fname!(currentplot(f), args...)
    @eval @doc TICKS_DOC function $fname(args...)
        f = gcf()
        $fname!(currentplot(f), args...)
        return f
    end

    # xlim, etc.
    fname! = Symbol(ax, :lim!)
    fname = Symbol(ax, :lim)
    @eval function $fname!(p::PlotObject, limits, adjust=false)
        _config_axislimits!(Symbol($ax), p, limits, adjust)
        p.attributes[Symbol($ax, :lim)] = limits
    end
    @eval function $fname!(p::PlotObject, minval::Union{Nothing, Number}, maxval::Union{Nothing, Number}, adjust=false)
        $fname!(p, (minval, maxval), adjust)
    end
    @eval $fname!(p::PlotObject) = $fname!(p, (nothing, nothing))
    @eval $fname!(f::Figure, args...) = $fname!(currentplot(f), args...)
    @eval @doc AXISLIM_DOC function $fname(args...)
        f = gcf()
        $fname!(currentplot(f), args...)
        return f
    end

    # xlog, xflip, etc.
    for (attr, docstr) ∈ (("log", :AXISLOG_DOC), ("flip", :AXISFLIP_DOC))
        fname! = Symbol(ax, attr, :!)
        fname = Symbol(ax, attr)
        @eval function $fname!(p::PlotObject, flag)
            if p.axes.kind ∈ (:axes2d, :axes3d)
                p.attributes[Symbol($ax, $attr)] = flag
                newscale = set_scale(; p.attributes...)
                if p.axes.options[:scale] != newscale
                    p.axes.options[:scale] = newscale
                    axlimits = get(p.attributes, Symbol($ax, :lim), (nothing, nothing))
                    adjust = !get(p.attributes, Symbol($ax, :log), false)
                    _config_axislimits!(Symbol($ax), p, axlimits, adjust)
                end
            end
            return nothing
        end
        @eval $fname!(f::Figure, args...) = $fname!(currentplot(f), args...)
        @eval @doc $docstr function $fname(args...)
            f = gcf()
            $fname!(currentplot(f), args...)
            return f
        end
    end
end

const TICKLABELS_DOC = """
    xticklabels(f)
    yticklabels(f)

Customize the string of the X and Y axes tick labels.

The labels of the tick axis can be defined by a function
with one argument (the numeric value of the tick position) that
returns a string, or by an array of strings that are located
sequentially at X = 1, 2, etc.

Use the keyword argument `xticklabels=s`, etc. in plotting functions, to
set the axis tick labels during the creation of plots.

# Examples

```julia
# Label the range (0-1) of the Y-axis as percent values
yticklabels(p -> Base.Printf.@sprintf("%0.0f%%", 100p))
# Label the X-axis with a sequence of strings
xticklabels(["first", "second", "third"])
```
"""

for ax = ("x", "y")
    fname! = Symbol(ax, :ticklabels!)
    fname = Symbol(ax, :ticklabels)
    @eval function $fname!(p::PlotObject, s)
        set_ticklabels!(p.axes.ticklabels; $fname = s)
        p.attributes[Symbol($ax, :ticklabels)] = s
    end
    @eval $fname!(f::Figure, s) = $fname!(currentplot(f), s)
    @eval @doc TICKLABELS_DOC function $fname(s)
        f = gcf()
        $fname!(currentplot(f), s)
        return f
    end
end

# Grid
function grid!(p::PlotObject, flag)
    p.axes.options[:grid] = Int(flag)
    p.attributes[:grid] = flag
end

grid!(f::Figure, flag) = grid!(currentplot(f), flag)

"""
    grid(flag::Bool)

Draw or disable the grid of the current plot axes.

Use the keyword argument `grid=<true/false>`, etc. in plotting functions, to
set the grid during the creation of plots.
"""
function grid(flag)
    f = gcf()
    grid!(currentplot(f), flag)
    return f
end

# Colorbar
colorbar!(p::PlotObject, flag::Bool) = (p.attributes[:colorbar] = flag)

function colorbar!(p::PlotObject, levels::Integer)
    p.colorbar = Colorbar(p.axes, levels)
    colorbar!(p, true)
end

colorbar!(f::Figure, flag) = colorbar!(currentplot(f), flag)

"""
    colorbar(flag::Bool)
    colorbar(levels::Integer)

Set the color bar of the current plot.

The input argument can be a `Bool` (`true` or `false`) to show or hide
the colorbar -- if it is available, or an `Integer` to set the number
of levels shown in the color bar (256 levels by default).

Color bars are only presented when there is actual color data in the plot,
regardless of the usage of this function.

Use the keyword argument `colorbar=<true/false>`, etc. in plotting functions, to
enable or disable the color bar during the creation of plots.
"""
function colorbar(flag)
    f = gcf()
    colorbar!(currentplot(f), flag)
    return f
end

# Aspect ratio
function aspectratio!(p::PlotObject, r)
    margins = plotmargins(p.legend, p.colorbar; p.attributes...)
    set_ratio!(p.viewport.inner, r, margins)
    p.attributes[:ratio] = r
end

aspectratio!(f::Figure, r) = aspectratio!(currentplot(f), r)

"""
    aspectratio(r)

Set the aspect of the current plot to a given width : height ratio.

Use the keyword argument `aspectratio=r`, etc. in plotting functions, to
set the aspect ratio during the creation of plots.

# Examples

```julia
$(_example("aspectratio"))
```
"""
function aspectratio(r)
    f = gcf()
    aspectratio!(currentplot(f), r)
    return f
end

# Radians in polar axes

function radians!(p::PlotObject, flag)
    if p.axes.kind ≠ :polar
        return nothing
    end
    p.axes.options[:radians] = Int(flag)
    p.attributes[:radians] = flag
end

radians!(f::Figure, flag) = radians!(currentplot(f), flag)

"""
    radians(flag::Bool)

Set the scale of angles in polar plots.

Use `radians(true)` to represent angles in radians (default setting),
and `radians(false)` to represent them in degrees.

This operation only modifies the guides of the polar plot grid lines.
The existing geometries are left without changes

Use the keyword argument `radians=<true/false>`, etc. in plotting functions, to
set the scale of angles during the creation of polar plots.

# Example

```julia
# Example data
θ = LinRange(0, 2π, 40)
r = sin.(θ)
# Draw the polar plot (by default in radians)
polar(θ, r)
# Change the angula scale
radians(false)
```
"""
function radians(flag)
    f = gcf()
    radians!(currentplot(f), flag)
    return f
end

# Pan and zoom

function panzoom!(p::PlotObject, x, y, r = 0.0)
    GR.savestate()
    GR.setviewport(p.viewport.inner...)
    GR.setwindow(p.axes.ranges[:x]..., p.axes.ranges[:y]...)
    xmin, xmax, ymin, ymax = GR.panzoom(x, y, r)
    GR.restorestate()
    xlim!(p, (xmin, xmax))
    ylim!(p, (ymin, ymax))
    return nothing
end

panzoom!(f::Figure, args...) = panzoom!(currentplot(f), args...)

"""
    panzoom(x, y[, s = 0])

Pan/zoom the axes of the current plot.

The focus of the zoom is set at a point with an offset of `(x, y)` units
in normalized device coordinates (NDC) from the center of the current axes.
The corners of the axes are linearly displaced towards that point,
such that the size of the new axes is `s` times their original size.

If `s` is set to 0 (the default value), the center of the axes
is displaced at the focus, without resizing-

# Example

```julia
# Move the center 1 unit right and 0.2 up (NDC)
panzoom(1, 0.2)
# Reduce the focus of the axes to half their size
# focusing on the previous point
panzoom(1, 0.2, 0.5)
```
"""
function panzoom(args...)
    f = gcf()
    panzoom!(currentplot(f), args...)
    return f
end

# zoom for axes2d and axes3d with gr3
zoom2d!(p, r) = panzoom!(p, 0.0, 0.0, r)

function zoomgr3!(p::PlotObject, r)
    p.axes.camera[1:3] ./= r
    return nothing
end

function zoom!(p::PlotObject, r)
    if p.axes.kind == :axes2d
        zoom2d!(p, r)
    elseif get(p.axes.options, :render3d, 0) == 2
        zoomgr3!(p, r)
    end
end

zoom!(f::Figure, r) = zoom!(currentplot(f), r)

"""
    zoom(r)

Zoom the plot by the ratio indicated by `r`.

In two-dimensional plots, the "zoomed" axes are centered around the same point,
but proportionally resized to `r` times the original size.

In three-dimensional scenes defined with "camera" settings
(e.g. in [`isosurface`](@ref) plots), the camera distance is divided by `r`.

# Examples

```julia
# Reduce the axes to half their size
zoom(0.5)
```
"""
function zoom(r)
    f = gcf()
    zoom!(currentplot(f), r)
    return f
end

# 3-D perspectives

function viewpoint!(p::PlotObject, rotation, tilt)
    p.axes.perspective .= [rotation, tilt]
    if get(p.axes.options, :render3d, 0) == 2
        distance = norm(view(p.axes.camera, 1:3))
        p.axes.camera .= set_camera(distance, rotation, tilt)
    end
    return nothing
end

viewpoint!(f::Figure, rotation, tilt) = viewpoint!(currentplot(f), rotation, tilt)

"""
    viewpoint(rotation, tilt)

Set the viewpoint of three-dimensional plots.

`rotation` and `tilt` must be integer values that indicate the
"azimuth" and "elevation" angles of the line of sight (in degrees).

If both angles are zero, the plot is viewed in the direction of the Y axis
(i.e. the X-Z plane is seen). Positive `rotation` values mean a
counterclockwise rotation of the line of sight (or a clockwise rotation of the scene)
around the vertical (Z) axis. Positive `tilt` values mean an ascension
of the view point.

# Examples

```julia
# Reset the view to the X-Y plane
# (rotation=0, tilt=90)
viewpoint(0, 90)
```
"""
function viewpoint(rotation, tilt)
    f = gcf()
    viewpoint!(currentplot(f), rotation, tilt)
    return f
end

function rotate!(p::PlotObject, angle)
    p.axes.perspective[1] += angle
    if get(p.axes.options, :render3d, 0) == 2
        _rotate!(view(p.axes.camera, 1:3), angle)
        _rotate!(view(p.axes.camera, 7:9), angle)
    end
    return nothing
end

rotate!(f::Figure, angle) = rotate!(currentfigure(f), angle)

function tilt!(p::PlotObject, angle)
    p.axes.perspective[2] += angle
    if get(p.axes.options, :render3d, 0) == 2
        rotation = p.axes.perspective[1]
        camera_position = view(p.axes.camera, 1:3)
        _rotate!(camera_position, -rotation)
        _tilt!(camera_position, angle)
        _rotate!(camera_position, rotation)
        up_vector = view(p.axes.camera, 7:9)
        _rotate!(up_vector, -rotation)
        _tilt!(up_vector, angle)
        _rotate!(up_vector, rotation)
    end
    return nothing
end

tilt!(f::Figure, angle) = tilt!(currentfigure(f), angle)

"""
    rotate(angle::Int)

Rotate the viewpoint of the current plot
by `angle` degrees around the vertical axis of the scene,
with respect to its current position.

# Examples

```julia
# Rotate 10 degrees to the right
rotate(10)
```
"""
function rotate(angle)
    f = gcf()
    rotate!(currentplot(f), angle)
    return f
end

"""
    tilt(angle::Int)

Tilt (elevate) the viewpoint of the current plot
by `angle` degrees over the horizontal plane,
with respect to its current position.

# Examples

```julia
# Tilt 10 degrees up
tilt(10)
```
"""
function tilt(angle)
    f = gcf()
    tilt!(currentplot(f), angle)
    return f
end

# Only for 3-D scenes with gr3

movefocus!(p::PlotObject, target) = _focus!(p.axes.camera, target)
movefocus!(f::Figure, target) = movefocus!(currentplot(f), target)

"""
    movefocus(target)

Rotate the camera view axis, moving the focus to the `target` point.

This only affects 3-D scenes created with camera settings, e.g.
[`isosurface`](@ref) plots. Moving the focus point rotates the camera
without changing its position; in order to rotate the camera around
the center of the scene, use the functions
[`rotate`](@ref), [`tilt`](@ref) or [`viewpoint`](@ref).

# Examples

```julia
# Move the focus to the point (1.0, 0.5, 0.0)
movefocus([1.0, 0.5, 0.0])
```
"""
function movefocus(target)
    f = gcf()
    movefocus!(currentplot(f), target)
    return f
end

function turncamera!(p::PlotObject, angle)
    params = p.axes.camera
    # Rotate up vector towards right vector around axis
    axis = normalize([params[4]-params[1], params[5]-params[2], params[6]-params[3]])
    up_vector = params[7:9]
    right_vector = axis × up_vector
    up_vector .= cosd(angle).*up_vector .+ sind(angle).*right_vector
    return nothing
end

turncamera!(f::Figure, angle) = turncamera!(currentplot(f), angle)

"""
    turncamera(angle)

Turn the orientation of the camera by `angle` degrees
around its view axis (only for 3-D scenes created with camera settings).

# Examples

```julia
# Turn the perspective 10 degrees
turncamera(10)
```
"""
function turncamera(angle)
    f = gcf()
    turncamera!(currentplot(f), angle)
    return f
end

"""
    colormap!(p, cmap)

Apply a colormap `cmap` to the given plot `p`, which can be a `PlotObject`,
or a `Figure` (in such case the colormap is applied to all the plots contained in it).

The value of `cmap` can be the number or the name of any of the
[GR built-in colormaps](https://gr-framework.org/colormaps.html)
(see [`colormap`](@ref) for more details).

Use the keyword argument `colormap` in plotting functions, to set a particular
colormap during the creation of plots (in this case it can only be identified
by its number).

# Examples

```julia
# Create a surface plot with the "grayscale" colormap (2)
surface(x, y, z, colormap=2)
# Change it to the "viridis" colormap
colormap!(gcf(), "viridis")
```
"""
function colormap!(p::PlotObject, cmap)
    p.attributes[:colormap] = Int(cmap)
    return nothing
end

function colormap!(p::PlotObject, cmap::AbstractString)
    cmap = lowercase(replace(cmap, (' ', '_') => ""))
    colormap!(p, COLORMAPS[cmap])
end

function colormap!(f::Figure, cmap)
    for p in f.plots
        colormap!(p, cmap)
    end
    return f
end

"""
    colorscheme!(p, scheme)

Apply a color `scheme` to the given plot `p`, which can be a `PlotObject`,
or a `Figure` (in such case the scheme is applied to all the plots contained in it).

The value of `scheme` can be the number or the name of any available
color scheme (see [`colorscheme`](@ref) for more details).

Use the keyword argument `scheme` in plotting functions, to set a particular
color scheme during the creation of plots (in this case only the number of
an already exisiting scheme is allowed).

# Examples

```julia
# Create a plot with a dark scheme (2)
plot(x, y, scheme=2)
# Change it to the standard light scheme
colorscheme!(currentplot(), "light")
```
"""
function colorscheme!(p::PlotObject, scheme)
    p.attributes[:scheme] = Int(scheme)
    return nothing
end

function colorscheme!(p::PlotObject, scheme::AbstractString)
    scheme = replace(scheme, " " => "")
    scheme = lowercase(scheme)
    scheme_dict = Dict("none" => 0, "light" => 1, "dark" => 2,
        "solarizedlight" => 3, "solarizeddark" => 4)
    colorscheme!(p, scheme_dict[scheme])
end

function colorscheme!(f::Figure, scheme)
    for p in f.plots
        colorscheme!(p, scheme)
    end
    return f
end

# Custom background

"""
    background!(p, bgcolor[, alpha])

Add a custom background color to the given plot object or to all the plots
inside the given figure. See [`background`](@ref) for more details.
"""
function background!(p::PlotObject, bgcolor)
    p.attributes[:backgroundcolor] = Int(bgcolor)
    return nothing
end

function background!(p::PlotObject, bgcolor, alpha)
    p.attributes[:backgroundcolor] = Int(bgcolor)
    p.attributes[:backgroundalpha] = alpha
    return nothing
end

background!(p::PlotObject, ::Nothing) = background!(p, -1)

function background!(f::Figure, args...)
    for p in f.plots
        background!(p, args...)
    end
    return f
end

"""
    background(color[, alpha])

Add a custom background color to the current figure.

The argument can be an hexadecimal color code or `nothing` for a transparent
background. A partially transparent color can be defined adding the alpha
value between 0 and 1 as second argument.

Use the keyword arguments `backgroundcolor` and `backgroundalpha`
in plotting functions, to set a particular background color configuration
during the creation of plots.

This overrides the default background defined by the [`colorscheme`](@ref) for
the area outside the axes and legends of all the plots contained in the figure.
Use [`background!`](@ref) to modify the background of individual subplots.

# Examples
```julia
# Create a plot with light blue background
plot(x, y, backgroundcolor=0x88ccff)
# Remove the background
background(nothing)
```
"""
background(args...) = background!(gcf(), args...)


text:
module Hiop
using Libdl
using LinearAlgebra

export solveProblem
export HiopProblem

const libhiop = "$(ENV["JULIA_HIOP_LIBRARY_PATH"])/libhiop.$(Libdl.dlext)"
const HIOP_INFINITY = 1e+20

mutable struct cHiopProblem
    refcppHiop::Ptr{Cvoid}
    hiopinterface::Ptr{Cvoid}
    jprob::Ptr{Cvoid}
    solution::Ptr{Cdouble}
    obj_value::Cdouble
    get_starting_point::Ptr{Cvoid}
    get_prob_sizes::Ptr{Cvoid}
    get_vars_info::Ptr{Cvoid}
    get_cons_info::Ptr{Cvoid}
    eval_f::Ptr{Cvoid}
    eval_grad_f::Ptr{Cvoid}
    eval_cons::Ptr{Cvoid}
    get_sparse_dense_blocks_info::Ptr{Cvoid}
    eval_Jac_cons::Ptr{Cvoid}
    eval_Hess_Lagr::Ptr{Cvoid}
end

# Forward declarations coming after the constructor
function get_starting_point_wrapper end
function get_sparse_dense_blocks_info_wrapper end
function get_prob_sizes_wrapper end
function get_vars_info_wrapper end
function get_cons_info_wrapper end
function eval_f_wrapper end
function eval_grad_f_wrapper end
function eval_g_wrapper end
function eval_jac_g_wrapper end
function eval_h_wrapper end
function freeProblem end

mutable struct HiopProblem
    cprob::cHiopProblem  # Reference to the C data structure
    n::Int64  # Num vars
    m::Int64  # Num cons
    ns::Int64 # Sparse whatever
    nd::Int64 # Dense whatever
    nx_sparse::Int32
    nx_dense::Int32
    nnz_sparse_Jaceq::Int32
    nnz_sparse_Jacineq::Int32
    jacidxlist::Array{Int64}
    nnz_sparse_Hess_Lagr_SS::Int32
    nnz_sparse_Hess_Lagr_SD::Int32
    hesidxlist::Array{Int64}
    x::Vector{Float64}  # Starting and final solution
    x_L::Vector{Float64}  # Starting and final solution
    x_U::Vector{Float64}  # Starting and final solution
    g::Vector{Float64}  # Final constraint values
    g_L::Vector{Float64}  # Final constraint values
    g_U::Vector{Float64}  # Final constraint values
    mult_g::Vector{Float64} # lagrange multipliers on constraints
    mult_x_L::Vector{Float64} # lagrange multipliers on lower bounds
    mult_x_U::Vector{Float64} # lagrange multipliers on upper bounds
    obj_val::Float64  # Final objective
    status::Int  # Final status

    # Callbacks
    eval_f::Function
    eval_g::Function
    eval_grad_f::Function
    eval_jac_g::Function
    eval_h  # Can be nothing
    x0::Vector{Float64}
    user_data::Any

    function HiopProblem(ns::Int, nd::Int, nx_sparse::Int32, nx_dense::Int32,
                         nnz_sparse_Jaceq::Int32, nnz_sparse_Jacineq::Int32, nnz_sparse_Hess_Lagr_SS::Int32, nnz_sparse_Hess_Lagr_SD::Int32,
                         n::Int, x_L::Vector{Float64}, x_U::Vector{Float64},
                         m::Int, g_L::Vector{Float64}, g_U::Vector{Float64},
                         eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h = nothing, user_data = nothing)
        # Wrap callbacks
        prob = new(cHiopProblem(C_NULL, C_NULL, C_NULL, C_NULL, 0.0, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL),
                   n, m, ns, nd,
                   nx_sparse, nx_dense,
                   nnz_sparse_Jaceq, nnz_sparse_Jacineq, Vector{Int64}(),
        nnz_sparse_Hess_Lagr_SS, nnz_sparse_Hess_Lagr_SD, Vector{Int64}(),
        zeros(Float64, n), x_L, x_U,
        zeros(Float64, m), g_L, g_U,
        zeros(Float64,m),
        zeros(Float64,n), zeros(Float64,n), 0.0, 0,
        eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h, zeros(Float64, n), user_data)

        prob.cprob.get_starting_point = @cfunction(get_starting_point_wrapper, Cint,
                                                   (Clonglong, Ptr{Cdouble}, Ptr{Cvoid}))
        prob.cprob.get_sparse_dense_blocks_info = @cfunction(get_sparse_dense_blocks_info_wrapper, Cint,
                                                             (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cvoid}))
        prob.cprob.get_prob_sizes = @cfunction(get_prob_sizes_wrapper, Cint,
                                               (Ptr{Clonglong}, Ptr{Clonglong}, Ptr{Cvoid}))
        prob.cprob.get_vars_info = @cfunction(get_vars_info_wrapper, Cint,
                                              (Ptr{Clonglong}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
        prob.cprob.get_cons_info = @cfunction(get_cons_info_wrapper, Cint,
                                              (Ptr{Clonglong}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
        prob.cprob.eval_f = @cfunction(eval_f_wrapper, Cint,
                                       (Clonglong, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cvoid}))
        prob.cprob.eval_grad_f = @cfunction(eval_grad_f_wrapper, Cint,
                                            (Clonglong, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cvoid}))
        prob.cprob.eval_cons = @cfunction(eval_g_wrapper, Cint,
                                          (Clonglong, Clonglong, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cvoid}))
        prob.cprob.eval_Jac_cons = @cfunction(eval_jac_g_wrapper, Cint,
                                              (Clonglong, Clonglong, Ptr{Cdouble}, Cint, Clonglong, Clonglong,
                                               Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
        prob.cprob.eval_Hess_Lagr = @cfunction(eval_h_wrapper, Cint,
                                               (Clonglong, Clonglong,
                                                Ptr{Cdouble}, Cint, Cdouble,
                                                Ptr{Cdouble}, Cint,
                                                Clonglong, Clonglong,
                                                Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
                                                Ptr{Cdouble},
                                                Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cvoid}))
        prob.cprob.jprob = pointer_from_objref(prob)
        ret = ccall((:hiop_createProblem, libhiop),
                    Cint, (Ptr{cHiopProblem},),
                    pointer_from_objref(prob.cprob))
        if ret != 0
            error("HiOp: Failed to construct problem.")
        end
        return prob
    end
end

# From Ipopt/src/Interfaces/IpReturnCodes_inc.h
ApplicationReturnStatus = Dict(
    0=>:Solve_Succeeded,
    1=>:Solved_To_Acceptable_Level,
    2=>:Infeasible_Problem_Detected,
    3=>:Search_Direction_Becomes_Too_Small,
    4=>:Diverging_Iterates,
    5=>:User_Requested_Stop,
    6=>:Feasible_Point_Found,
    -1=>:Maximum_Iterations_Exceeded,
    -2=>:Restoration_Failed,
    -3=>:Error_In_Step_Computation,
    -4=>:Maximum_CpuTime_Exceeded,
    -10=>:Not_Enough_Degrees_Of_Freedom,
    -11=>:Invalid_Problem_Definition,
    -12=>:Invalid_Option,
    -13=>:Invalid_Number_Detected,
    -100=>:Unrecoverable_Exception,
    -101=>:NonIpopt_Exception_Thrown,
    -102=>:Insufficient_Memory,
    -199=>:Internal_Error,
)

###########################################################################
# Callback wrappers
###########################################################################

function get_starting_point_wrapper(n::Clonglong, x0_::Ptr{Cdouble}, prob_::Ptr{Cvoid})
    prob = unsafe_pointer_to_objref(prob_)::HiopProblem
    x0 = unsafe_wrap(Array{Float64}, x0_, prob.n)
    x0 .= prob.x0
    return Int32(1)
end
function get_prob_sizes_wrapper(n_::Ptr{Clonglong}, m_::Ptr{Clonglong}, prob_::Ptr{Cvoid})
    prob = unsafe_pointer_to_objref(prob_)::HiopProblem
    n = unsafe_wrap(Array{Int64}, n_, 1)
    m = unsafe_wrap(Array{Int64}, m_, 1)
    n[1] = prob.n
    m[1] = prob.m
    return Int32(1)
end

function get_vars_info_wrapper(n_::Ptr{Clonglong}, xlow_::Ptr{Cdouble}, xupp_::Ptr{Cdouble}, prob_::Ptr{Cvoid})
    prob = unsafe_pointer_to_objref(prob_)::HiopProblem
    xlow = unsafe_wrap(Array{Float64}, xlow_, prob.n)
    xupp = unsafe_wrap(Array{Float64}, xupp_, prob.n)
    xlow .= prob.x_L
    xupp .= prob.x_U
    return Int32(1)
end

function get_cons_info_wrapper(m_::Ptr{Clonglong}, clow_::Ptr{Cdouble}, cupp_::Ptr{Cdouble}, prob_::Ptr{Cvoid})
    prob = unsafe_pointer_to_objref(prob_)::HiopProblem
    clow = unsafe_wrap(Array{Float64}, clow_, prob.m)
    cupp = unsafe_wrap(Array{Float64}, cupp_, prob.m)
    clow .= prob.g_L
    cupp .= prob.g_U
    return Int32(1)
end

function get_sparse_dense_blocks_info_wrapper(nx_sparse_::Ptr{Cint}, nx_dense_::Ptr{Cint},
                                              nnz_sparse_Jaceq_::Ptr{Cint},
                                              nnz_sparse_Jacineq_::Ptr{Cint},
                                              nnz_sparse_Hess_Lagr_SS_::Ptr{Cint},
                                              nnz_sparse_Hess_Lagr_SD_::Ptr{Cint},
                                              prob_::Ptr{Cvoid})
    prob = unsafe_pointer_to_objref(prob_)::HiopProblem
    nx_sparse = unsafe_wrap(Array{Int32}, nx_sparse_, 1)
    nx_dense = unsafe_wrap(Array{Int32}, nx_dense_, 1)
    nnz_sparse_Jaceq = unsafe_wrap(Array{Int32}, nnz_sparse_Jaceq_, 1)
    nnz_sparse_Jacineq = unsafe_wrap(Array{Int32}, nnz_sparse_Jacineq_, 1)
    nnz_sparse_Hess_Lagr_SS = unsafe_wrap(Array{Int32}, nnz_sparse_Hess_Lagr_SS_, 1)
    nnz_sparse_Hess_Lagr_SD = unsafe_wrap(Array{Int32}, nnz_sparse_Hess_Lagr_SD_, 1)

    nx_sparse[1] = prob.nx_sparse
    nx_dense[1] = prob.nx_dense
    nnz_sparse_Jaceq[1] = prob.nnz_sparse_Jaceq
    nnz_sparse_Jacineq[1] = prob.nnz_sparse_Jacineq
    nnz_sparse_Hess_Lagr_SS[1] = prob.nnz_sparse_Hess_Lagr_SS
    nnz_sparse_Hess_Lagr_SD[1] = prob.nnz_sparse_Hess_Lagr_SD
    return Int32(1)
end

function eval_f_wrapper(n::Clonglong, x_ptr::Ptr{Float64},
                        new_x::Cint, obj_ptr::Ptr{Float64}, prob_::Ptr{Cvoid})
    # Extract Julia the problem from the pointer
    prob = unsafe_pointer_to_objref(prob_)::HiopProblem
    # Calculate the new objective
    prob.obj_val = convert(Float64, prob.eval_f(unsafe_wrap(Array,x_ptr, Int(n)), prob))::Float64
    # Fill out the pointer
    unsafe_store!(obj_ptr, prob.obj_val)
    # Done
    return Int32(1)
end

# Constraints (eval_g)
function eval_g_wrapper(n::Clonglong, m::Clonglong,
                        x_ptr::Ptr{Float64}, new_x::Cint, g_ptr::Ptr{Float64},
                        user_data::Ptr{Cvoid})
    # Extract Julia the problem from the pointer
    prob = unsafe_pointer_to_objref(user_data)::HiopProblem
    x = unsafe_wrap(Array{Float64}, x_ptr, prob.n)
    # Calculate the new constraint values
    new_g = unsafe_wrap(Array, g_ptr, Int(m))
    prob.eval_g(x, new_g, prob)
    # Done
    return Int32(1)
end

# Objective gradient (eval_grad_f)
function eval_grad_f_wrapper(n::Clonglong, x_ptr::Ptr{Float64},
                             new_x::Cint, grad_f_ptr::Ptr{Float64},
                             user_data::Ptr{Cvoid})
    # Extract Julia the problem from the pointer
    prob = unsafe_pointer_to_objref(user_data)::HiopProblem
    # Calculate the gradient
    new_grad_f = unsafe_wrap(Array, grad_f_ptr, Int(n))
    prob.eval_grad_f(unsafe_wrap(Array,x_ptr, Int(n)), new_grad_f, prob)
    # Done
    return Int32(1)
end

# Jacobian (eval_jac_g)
function eval_jac_g_wrapper(n::Clonglong, m::Clonglong,
                            x_ptr::Ptr{Cdouble}, new_x::Cint,
                            nsparse::Clonglong, ndense::Clonglong,
                            nnzJacS::Cint, iJacS_ptr::Ptr{Cint},
                            jJacS_ptr::Ptr{Cint}, MJacS_ptr::Ptr{Cdouble},
                            JacD_ptr::Ptr{Cdouble}, user_data::Ptr{Cvoid})
    # Extract Julia the problem from the pointer
    prob = unsafe_pointer_to_objref(user_data)::HiopProblem
    # Determine mode
    mode = Vector{Symbol}()
    if iJacS_ptr != C_NULL && jJacS_ptr != C_NULL
        push!(mode, :Structure)
    end
    if MJacS_ptr != C_NULL
        push!(mode, :Sparse)
    end
    if JacD_ptr != C_NULL
        push!(mode, :Dense)
    end
    x = unsafe_wrap(Array, x_ptr, Int(n))
    iJacS = unsafe_wrap(Array, iJacS_ptr, Int(nnzJacS))
    jJacS = unsafe_wrap(Array, jJacS_ptr, Int(nnzJacS))
    MJacS = unsafe_wrap(Array, MJacS_ptr, Int(nnzJacS))
    # @show n, m, nsparse, ndense, prob.ns, prob.nd
    JacD = unsafe_wrap(Array, JacD_ptr, prob.m * prob.nd)
    prob.eval_jac_g(mode, x, iJacS, jJacS, MJacS, JacD, prob)
    # Done
    return Int32(1)
end

# Hessian
function eval_h_wrapper(n::Clonglong, m::Clonglong,
                        x_ptr::Ptr{Cdouble}, new_x::Cint, obj_factor::Cdouble,
                        lambda_ptr::Ptr{Cdouble}, new_lambda::Cint,
                        nsparse::Clonglong, ndense::Clonglong,
                        nnzHSS::Cint, iHSS_ptr::Ptr{Cint},
                        jHSS_ptr::Ptr{Cint}, MHSS_ptr::Ptr{Cdouble},
                        HDD_ptr::Ptr{Cdouble},
                        nnzHSD::Cint, iHSD_ptr::Ptr{Cint},
                        jHSD_ptr::Ptr{Cint}, MHSD_ptr::Ptr{Cdouble},
                        user_data::Ptr{Cvoid})
    # Extract Julia the problem from the pointer
    prob = unsafe_pointer_to_objref(user_data)::HiopProblem
    # Determine mode
    mode = Vector{Symbol}()
    if iHSS_ptr != C_NULL && jHSS_ptr != C_NULL
        push!(mode, :Structure)
    end
    if MHSS_ptr != C_NULL
        push!(mode, :Sparse)
    end
    if HDD_ptr != C_NULL
        push!(mode, :Dense)
    end
    x = unsafe_wrap(Array, x_ptr, Int(n))
    lambda = unsafe_wrap(Array, lambda_ptr, Int(m))
    iHSS = unsafe_wrap(Array, iHSS_ptr, Int(nnzHSS))
    jHSS = unsafe_wrap(Array, jHSS_ptr, Int(nnzHSS))
    MHSS = unsafe_wrap(Array, MHSS_ptr, Int(nnzHSS))
    HDD = unsafe_wrap(Array, HDD_ptr, Int(prob.nd * prob.nd))
    iHSD = unsafe_wrap(Array, iHSD_ptr, Int(nnzHSD))
    jHSD = unsafe_wrap(Array, jHSD_ptr, Int(nnzHSD))
    MHSD = unsafe_wrap(Array, MHSD_ptr, Int(nnzHSD))
    if prob.eval_h != nothing
        prob.eval_h(mode, x, obj_factor, lambda, iHSS, jHSS, MHSS, HDD, iHSD, jHSD, MHSD, prob)
        # else
        #   MHSS .= 0.0
    end
    # Done
    return Int32(1)
end

###########################################################################
# C function wrappers
###########################################################################

function freeProblem(prob::HiopProblem)
    if prob.cprob.refcppHiop != C_NULL
        ccall((:hiop_destroyProblem, libhiop),
              Cint, (Ptr{cHiopProblem},),
              pointer_from_objref(prob.cprob))
        prob.cprob.refcppHiop = C_NULL
    end
end

function solveProblem(prob::HiopProblem)
    final_objval = [0.0]
    prob.cprob.solution = pointer(prob.x)
    ret = ccall((:hiop_solveProblem, libhiop),
                Cint, (Ptr{cHiopProblem}, ),
                pointer_from_objref(prob.cprob))
    prob.obj_val = prob.cprob.obj_value
    return Int(ret)
end

include("MOI_wrapper.jl")

end # module

ext: jl


text:
struct SuggestItem{T <: AbstractString}
    term::T
    distance::Int
    count::Int
end

SuggestItem(term::String32, distance::Int, count::Int) = SuggestItem{String}(String(term.s), distance, count)
term(si::SuggestItem) = si.term

Base.:isless(si1::SuggestItem, si2::SuggestItem) = (si1.distance < si2.distance) ||
    ((si1.distance == si2.distance) && (si1.count > si2.count)) ||
    ((si1.distance == si2.distance) && (si1.count == si2.count) && (si1.term < si2.term))

function set_options!(dict::SymSpell; kwargs...)
    for (opt, v) in kwargs
        if opt in fieldnames(LookupOptions)
            setfield!(dict.lookup, opt, v)
        else
            @warn "Lookup option '$opt' is not supported"
        end
    end
end

"""
    lookup(dict, phrase, max_edit_distance, include_unknown, ignore_token, transfer_casing)

Find suggested spellings for a given phrase word.
Parameters
----------
phrase : str
    The word being spell checked.
verbosity : :class:`Verbosity`
    The value controlling the quantity/closeness of the
    returned suggestions.
max_edit_distance : int, optional
    The maximum edit distance between phrase and suggested
    words. Set to :attr:`_max_dictionary_edit_distance` by
    default
include_unknown : bool, optional
    A flag to determine whether to include phrase word in
    suggestions, if no words within edit distance found.
ignore_token : regex pattern, optional
    A regex pattern describing what words/phrases to ignore and
    leave unchanged
transfer_casing : bool, optional
    A flag to determine whether the casing --- i.e., uppercase
    vs lowercase --- should be carried over from `phrase`.
Returns
-------
suggestions : list
    suggestions is a list of :class:`SuggestItem` objects
    representing suggested correct spellings for the phrase
    word, sorted by edit distance, and secondarily by count
    frequency.
Raises
------
ValueError
    If `max_edit_distance` is greater than
    :attr:`_max_dictionary_edit_distance`
"""
function lookup(dict, phrase::S; include_unknown = dict.lookup.include_unknown,
        ignore_token = dict.lookup.ignore_token, transfer_casing = dict.lookup.transfer_casing,
        verbosity::Union{String, Verbosity} = dict.lookup.verbosity,
        max_edit_distance = dict.lookup.max_edit_distance) where {S <: AbstractString}

        if verbosity isa String
            verbosity = verbosity == "top" ? VerbosityTOP : verbosity == "closest" ? VerbosityCLOSEST : VerbosityALL
        end
        lookup(dict, String32(phrase), include_unknown, ignore_token,
            transfer_casing, verbosity, max_edit_distance)
end

Base.:getindex(dict, phrase) = term.(lookup(dict, phrase))

# True Damerau-Levenshtein Edit Distance: adjust distance, if both distances>0
# We allow simultaneous edits (deletes) of max_edit_distance on on both the dictionary and
# the phrase term. For replaces and adjacent transposes the resulting edit distance stays
# <= max_edit_distance. For inserts and deletes the resulting edit distance might exceed
# max_edit_distance. To prevent suggestions of a higher edit distance, we need to calculate the
# resulting edit distance, if there are simultaneous edits on both sides.
# Example: (bank==bnak and bank==bink, but bank!=kanb and bank!=xban and bank!=baxn for
# max_edit_distance=1). Two deletes on each side of a pair makes them all equal, but the first two
# pairs have edit distance=1, the others edit distance=2.
function lookup(dict::SymSpell{T, K}, phrase::String32, include_unknown::Bool, ignore_token,
                transfer_casing::Bool, verbosity::Verbosity,
                max_edit_distance::Int) where {T, K}
    if max_edit_distance > dict.max_dictionary_edit_distance
        throw(ArgumentError("Distance $max_edit_distance larger than dictionary threshold $(dict.max_dictionary_edit_distance)"))
    end

    suggestions = SuggestItem{String}[]

    phrase_original = phrase
    if transfer_casing
        phrase = lowercase(phrase)
    end
    phrase_len = length(phrase)

    early_exit = let phrase = phrase_original, suggestions = suggestions, include_unknown = include_unknown, max_edit_distance = max_edit_distance
        function early_exit()
            if include_unknown && isempty(suggestions)
                push!(suggestions, SuggestItem(phrase, max_edit_distance + 1, 0))
            end
            suggestions
        end
    end

    # early exit - word is too big to possibly match any words
    if phrase_len - max_edit_distance > dict.max_length
        return early_exit()
    end

    # quick look for exact match
    if phrase_original in keys(dict.words_idx)
        push!(suggestions, SuggestItem(phrase, 0, dict.words[dict.words_idx[phrase_original]].count))
        if verbosity != VerbosityALL
            return suggestions
        end
        idx = dict.words_idx[phrase]
    else
        idx = K(0)
    end

    if ignore_token != nothing && occursin(ignore_token, phrase)
        suggestion_cnt = 1
        push!(suggestions, SuggestItem(phrase, 0, suggestion_cnt))
        # early exit - return exact match, unless caller wants all matches
        if verbosity != VerbosityALL
            return early_exit()
        end
    end

    # early termination, if we only want to check if word in
    # dictionary or get its frequency e.g. for word segmentation
    if max_edit_distance == 0
        return early_exit()
    end

    considered_deletes = Set{String}()
    considered_suggestions = Set{K}()

    # we considered the phrase already in the
    # 'phrase in keys(dict.words)' above
    if idx > 0
        push!(considered_suggestions, idx)
    end

    max_edit_distance_2 = max_edit_distance
    candidates = String[]

    # add original prefix
    phrase_prefix_len = min(phrase_len, dict.prefix_length)
    push!(candidates, String(phrase.s[1:phrase_prefix_len]))

    is_first = true

    local v2
    local v0

    while !isempty(candidates)
        candidate = popfirst!(candidates)
        candidate_len = length(candidate)
        len_diff = phrase_prefix_len - candidate_len
        # early termination: if candidate distance is already
        # higher than suggestion distance, than there are no better
        # suggestions to be expected
        if len_diff > max_edit_distance_2
            # skip to next candidate if Verbosity.ALL, look no
            # further if Verbosity.TOP or CLOSEST (candidates are
            # ordered by delete distance, so none are closer than
            # current)
            verbosity == VerbosityALL && continue
            break
        end

        if candidate in keys(dict.deletes)
            for suggestion_id in dict.deletes[candidate]
                @inbounds dict_word = dict.words[suggestion_id]
                suggestion = dict_word.word
                suggestion_len = length(dict_word.word)
                suggestion_cnt = dict_word.count

                !simple_filter(phrase_len, phrase_prefix_len,
                    suggestion_len, candidate_len,
                    max_edit_distance_2, dict.prefix_length) && continue

                # number of edits in prefix == maxeditdistance AND no identical suffix, then
                # editdistance > max_edit_distance and no need for Levenshtein calculation
                # (phraseLen >= prefixLength) && (suggestionLen >= prefixLength)
                if dict.prefix_length - max_edit_distance == candidate_len
                    min_distance = min(phrase_len, suggestion_len) - dict.prefix_length

                    min_distance > 1 &&
                        equal_suffixes(phrase, suggestion, min_distance) && continue
                    if min_distance > 0
                        p1 = phrase_len + 1 - min_distance
                        s1 = suggestion_len + 1 - min_distance
                        @inbounds if phrase[p1] != suggestion[s1]
                            @inbounds (phrase[p1 - 1] != suggestion[s1] || phrase[p1] != suggestion[s1 - 1]) && continue
                        end
                    end
                end

                # DL and considered_suggestions compete with each othere
                # In ALL mode it's usually better to test for already
                # verified suggestions
                if verbosity == VerbosityALL
                    suggestion_id in considered_suggestions && continue
                    push!(considered_suggestions, suggestion_id)
                end
                if is_first
                    v2 = Vector{Int}(undef, phrase_len + max_edit_distance)
                    v0 = similar(v2)
                    is_first = false
                end
                distance = evaluate!(phrase, suggestion, max_edit_distance_2, v0, v2)
                distance > max_edit_distance_2 && continue

                if verbosity != VerbosityALL
                    suggestion_id in considered_suggestions && continue
                    push!(considered_suggestions, suggestion_id)
                end

                # do not process higher distances than those already found, if verbosity<ALL (note:
                # max_edit_distance_2 will always equal max_edit_distance when Verbosity.ALL)
                if distance <= max_edit_distance_2
                    si = SuggestItem(suggestion, distance, suggestion_cnt)
                    if !isempty(suggestions)
                        if verbosity == VerbosityCLOSEST
                            # we will calculate DamLev distance
                            # only to the smallest found distance so far
                            if distance < max_edit_distance_2
                                empty!(suggestions)
                            end
                        elseif verbosity == VerbosityTOP
                            suggestions[1] = si < suggestions[1] ? si : suggestions[1]
                            max_edit_distance_2 = distance
                            continue
                        end
                    end

                    if verbosity != VerbosityALL
                        max_edit_distance_2 = distance
                    end
                    push!(suggestions, si)
                end
            end
        end

        # add edits: derive edits (deletes) from candidate (phrase)
        # and add them to candidates list. this is a recursive
        # process until the maximum edit distance has been reached
        if len_diff < max_edit_distance && candidate_len <= dict.prefix_length
            # do not create edits with edit distance smaller than
            # suggestions already found
            if verbosity != VerbosityALL && len_diff >= max_edit_distance_2
                continue
            end
            add_edits!(considered_deletes, candidates, candidate, candidate_len)
        end
    end

    if length(suggestions) > 1
        sort!(suggestions)
    end

    if transfer_casing
        for (i, s) in enumerate(suggestions)
            suggestions[i] = SuggestItem(transfer_casing_for_similar_text(phrase_original, s.term), s.distance, s.count)
        end
    end

    early_exit()

    return suggestions
end

@inline function simple_filter(phrase_len::Int, phrase_prefix_len::Int, suggestion_len::Int,
     candidate_len::Int, max_edit_distance_2::Int, prefix_length::Int)

    phrase_len - suggestion_len > max_edit_distance_2 && return false
    suggestion_len - phrase_len > max_edit_distance_2 && return false
    suggestion_len < candidate_len && return false

    suggestion_prefix_len = suggestion_len > prefix_length ? prefix_length : suggestion_len

    ((suggestion_prefix_len > phrase_prefix_len) &
        (suggestion_prefix_len - candidate_len > max_edit_distance_2)) &&
        return false

    if candidate_len == 0
        # suggestions which have no common chars with phrase (phrase_len<=max_edit_distance &&
        # suggestion_len<=max_edit_distance)

        distance = max(phrase_len, suggestion_len)
        distance > max_edit_distance_2 && return false
    end

    return true
end

@inline function equal_suffixes(phrase::String32, suggestion::String32, min_distance)
    id1 = length(phrase) + 2 - min_distance
    id2 = length(suggestion) + 2 - min_distance
    for i in 0:(length(phrase) - id1)
        @inbounds phrase[id1 + i] != suggestion[id2 + i] && return true
    end

    return false
end

@inline function add_edits!(considered_deletes, candidates, candidate, candidate_len)
    sz = ncodeunits(candidate)
    idx1 = 0
    idx2 = nextind(candidate, idx1)
    idx3 = idx2
    for i in eachindex(0:candidate_len-1)
        idx3 = nextind(candidate, idx3)
        # black magic, fast string with deleted character
        delete = Base._string_n(sz - idx3 + idx2)
        unsafe_copyto!(pointer(delete), pointer(candidate), idx2 - 1)
        unsafe_copyto!(pointer(delete, idx2), pointer(candidate, idx3), sz - idx3 + 1)
        idx1 = idx2
        idx2 = idx3
        if !(delete in considered_deletes)
            push!(considered_deletes, delete)
            push!(candidates, delete)
        end
    end
end


text:
using FANN
using Base.Test

net = ANN([1, 1, 1], [:sigmoid_symmetric, :linear], b=0.1, errorfunc=:linear)

nw = FANN.nweights(net)
@test nw == 4

FANN.weights!(net, Float64[1, 2, 3, 4])
w = FANN.weights(net)
@test w == Float64[1, 2, 3, 4]


text:
######################################################################
# lq.jl
#
# Solves LQ descent direction problems
######################################################################
export LQ, apply_LQ_gains

# TODO:
#  I'm pretty sure I can speed this up by removing following arrays:
#    P, G, r
#  I think I can just store the previous one?
function LQ(A::MF, B::MF, a::MF, b::MF, Q::MF, R::MF, N::Int)
	# Also needed for LQR
	P = Array(Matrix{Float64}, N+1)
	G = Array(Matrix{Float64}, N)
	K = Array(Matrix{Float64}, N)
	P[N+1] = Q

	# Solely for LQ
	r = Array(Vector{Float64}, N+1)
	#r[N+1] = zeros(size(B,1))
	r[N+1] = 0.5*a[:,N+1]
	println("RARARAR")
	C = Array(Vector{Float64}, N)

	# Sweep from the back
	for n = (N-1):-1:0
		G[n+1] = R + (B' * P[n+1+1] * B)
		K[n+1] = inv(G[n+1]) * B' * P[n+1+1] * A
		P[n+1] = Q + (A' * P[n+1+1] * A) - (K[n+1]' * G[n+1] * K[n+1])
		r[n+1] = (A'-K[n+1]'*B')r[n+1+1] + .5*a[:,n+1] - .5*K[n+1]'*b[:,n+1]
		C[n+1] = inv(G[n+1]) * (B'*r[n+1+1] + .5*b[:,n+1])
	end

	return K, C
end

# This appears to be the one currrently in use
function LQ(A::VMF, B::VMF, a::MF, b::MF, Q::MF, R::MF, N::Int)
	# Also needed for LQR
	P = Array{Matrix{Float64}}(N+1)
	G = Array{Matrix{Float64}}(N)
	K = Array{Matrix{Float64}}(N)
	P[N+1] = Q

	# Solely for LQ
	r = Array{Vector{Float64}}(N+1)
	r[N+1] = 0.5*a[:,N+1]
	C = Array{Vector{Float64}}(N)

	# Sweep from the back
	# really n = (N-1):-1:0, but then all indices need an extra +1
	#  this is annoying so I just do n = N:-1:1
	#for n = (N-1):-1:0
	for n = N:-1:1
		# Temporary... only for testing example for paper
		#h = .6
		#A[n] = eye(4)
		#A[n][3,1] = -h; A[n][4,2] = -h; A[n][1,3] = h; A[n][2,4] = h
		#B[n] = zeros(4,2)
		#B[n][3,1] = h; B[n][4,2] = h
		# End Temporary

		G[n] = R + (B[n]' * P[n+1] * B[n])
		Ginv = inv(G[n])
		K[n] = Ginv * B[n]' * P[n+1] * A[n]
		Kp = K[n]'
		P[n] = Q + (A[n]' * P[n+1] * A[n]) - (Kp * G[n] * K[n])
		r[n] = (A[n]'-Kp*B[n]')r[n+1] + .5*a[:,n] - .5*Kp*b[:,n]
		C[n] = Ginv * (B[n]'*r[n+1] + .5*b[:,n])
	end

	return K, C
end

# Experimental...
# used in linear pto
function LQ(K::VMF, Ginv::VMF, A::VMF, B::VMF, a::MF, b::MF, Q::MF, R::MF, N::Int)

	# Solely for LQ
    r = Array{Vector{Float64}}(N+1)
	r[N+1] = 0.5*a[:,N+1]
    C = Array{Vector{Float64}}(N)

	# Sweep from the back
	# really n = (N-1):-1:0, but then all indices need an extra +1
	#  this is annoying so I just do n = N:-1:1
	for n = N:-1:1
		Kp = K[n]'
		r[n] = (A[n]'-Kp*B[n]')r[n+1] + .5*a[:,n] - .5*Kp*b[:,n]
		C[n] = Ginv[n] * (B[n]'*r[n+1] + .5*b[:,n])
	end

	return C
end

# used in linear pto
function LQ2(A::VMF, B::VMF, Q::MF, R::MF, N::Int)
	# Also needed for LQR
    P = Array{Matrix{Float64}}(N+1)
    G = Array{Matrix{Float64}}(N)
    Ginv = Array{Matrix{Float64}}(N)
    K = Array{Matrix{Float64}}(N)
	P[N+1] = Q

	# Sweep from the back
	# really n = (N-1):-1:0, but then all indices need an extra +1
	#  this is annoying so I just do n = N:-1:1
	#for n = (N-1):-1:0
	for n = N:-1:1
		G[n] = R + (B[n]' * P[n+1] * B[n])
		Ginv[n] = inv(G[n])
		K[n] = Ginv[n] * B[n]' * P[n+1] * A[n]
		Kp = K[n]'
		P[n] = Q + (A[n]' * P[n+1] * A[n]) - (Kp * G[n] * K[n])
	end

	return K, Ginv
end

function apply_LQ_gains(A::MF, B::MF, K::Vector{MF}, C::VVF)
	N = length(K)

	z = [zeros(size(B,1))]
    v = Array{Vector{Float64}}(0)
	for n = 1:N
		push!(v, -K[n]*z[n] - C[n])
		push!(z, A*z[n] + B*v[n])
	end
	return z, v
end

function apply_LQ_gains(A::VMF, B::VMF, K::VMF, C::VVF)
	N = length(K)

	z = [zeros(size(B[1],1))]
	v = Vector{Float64}[]
	for n = 1:N
		push!(v, -K[n]*z[n] - C[n])
		push!(z, A[n]*z[n] + B[n]*v[n])
	end
	return z, v
end

# Nope, not here
#function apply_LQ_gains(d::Dynamics, K::Vector{MF}, C::VVF)
#	N = length(K)
#
#	z = [zeros(size(B,1))]
#	v = Array(Vector{Float64}, 0)
#	for n = 1:N
#		push!(v, -K[n]*z[n] - C[n])
#		push!(z, A*z[n] + B*v[n])
#	end
#	return z, v
#end


text:
"""
    lis_length(str_1::String, str_2::String)

The longest increasing subsequence (LIS) algorithm is dedicated to finding the longest
increasing subsequence for a given array.  This subsequence does not have to be contiguous
or unique.

    1. Figure out the number of possible different subsequences of a string with length `m` and `n`.
    2. Recontracted `for-loops` for defining the outer layer of the storing matrix `C`.
    3. Checking if `str_1[i] == str_2[j]` and updating the count in  matrix `C`.
    4.Return the latest element of `C[m,n]`.

For more information see: [https://en.wikipedia.org/wiki/Longest_increasing_subsequence](https://en.wikipedia.org/wiki/Longest_increasing_subsequence)

# Arguments
- `array::Array{Int64,1}`: array of integers unsortetd or sorted

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> arr = [10, 22, 9, 33, 21, 50, 41, 60]
julia> ClassicAlgorithmsCollections.lis_length(arr)
4
```
"""
function lis_length(array::Array{Int64,1})
    n = length(array)
    list = fill(1, n)
    for i in 2:n
        for j in 1:i
            if array[i] > array[j] && list[i] < list[j] + 1
                list[i] = list[j] + 1
            end
        end
    end

    result = 0
    for i in 1:n
        if list[i] > result
            result = list[i]
        end
    end
    return result
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
function Mapper2.ismappable(::Test3d, c::Component)
    return search_metadata(c, "mappable", true)
end

function Mapper2.canmap(T::Test3d, t::TaskgraphNode, c::Component)
    return Mapper2.ismappable(T, c)
end


text:
@testset "support" begin

@testset "command-line options" begin

code = """
    using LLVM
    LLVM.clopts("-version")
"""

out = Pipe()
run(pipeline(julia_cmd(`-e $code`), stdout=out, stderr=out))
close(out.in)

@test occursin("LLVM (http://llvm.org/)", read(out, String))

end

end


text:
const READ_SCOPE = "smartRead"
const WRITE_SCOPE = "smartWrite"


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAuthenticationV1BoundObjectReference <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    uid::Any # spec type: Union{ Nothing, String } # spec name: uid

    function IoK8sApiAuthenticationV1BoundObjectReference(;apiVersion=nothing, kind=nothing, name=nothing, uid=nothing)
        o = new()
        validate_property(IoK8sApiAuthenticationV1BoundObjectReference, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiAuthenticationV1BoundObjectReference, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiAuthenticationV1BoundObjectReference, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiAuthenticationV1BoundObjectReference, Symbol("uid"), uid)
        setfield!(o, Symbol("uid"), uid)
        o
    end
end # type IoK8sApiAuthenticationV1BoundObjectReference

const _property_map_IoK8sApiAuthenticationV1BoundObjectReference = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("name")=>Symbol("name"), Symbol("uid")=>Symbol("uid"))
const _property_types_IoK8sApiAuthenticationV1BoundObjectReference = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("name")=>"String", Symbol("uid")=>"String")
Base.propertynames(::Type{ IoK8sApiAuthenticationV1BoundObjectReference }) = collect(keys(_property_map_IoK8sApiAuthenticationV1BoundObjectReference))
Swagger.property_type(::Type{ IoK8sApiAuthenticationV1BoundObjectReference }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAuthenticationV1BoundObjectReference[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthenticationV1BoundObjectReference }, property_name::Symbol) =  _property_map_IoK8sApiAuthenticationV1BoundObjectReference[property_name]

function check_required(o::IoK8sApiAuthenticationV1BoundObjectReference)
    true
end

function validate_property(::Type{ IoK8sApiAuthenticationV1BoundObjectReference }, name::Symbol, val)
end


text:
# FASTA
# =====
#
# Reader and writer of the FASTA file format.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"Metadata for FASTA sequence records containing just a `description` field"
type FASTAMetadata
    description::StringField
end

function FASTAMetadata()
    return FASTAMetadata(StringField())
end

function Base.:(==)(a::FASTAMetadata, b::FASTAMetadata)
    return a.description == b.description
end

function Base.copy(metadata::FASTAMetadata)
    return FASTAMetadata(copy(metadata.description))
end

"""
    FASTASeqRecord(name, seq[, description=""])

Create a sequence record for the FASTA file format.
"""
typealias FASTASeqRecord{S} SeqRecord{S,FASTAMetadata}

function (::Type{FASTASeqRecord})(name::AbstractString,
                                  seq::Sequence,
                                  description::AbstractString="")
    return SeqRecord(name, seq, FASTAMetadata(description))
end

include("fai.jl")
include("reader.jl")
include("writer.jl")


text:
using PushGP: BoolOr, Dup

@testset "Examples from PushGP docs" begin

@testset "Example 1 from Push3 docs" begin
    # Simple example from the page https://faculty.hampshire.edu/lspector/push3-description.html
    # ( 2 3 INTEGER.* 4.1 5.2 FLOAT.+ TRUE FALSE BOOLEAN.OR )
    interp = PushInterpreter()
    p = [2, 3, Mul{Int}(), 4.1, 5.2, Add{Float64}(), true, false, BoolOr()]
    execute!(interp, p)
    @test length(stack(interp, Bool)) == 1
    @test last(stack(interp, Bool)) == true
    @test length(stack(interp, Float64)) == 1
    @test last(stack(interp, Float64)) == 9.3
    @test length(stack(interp, Int)) == 1
    @test last(stack(interp, Int)) == 6
end # @testset "Example 1 from Push3 docs" begin

@testset "Example 2 from Push3 docs" begin
    # Example 2 from the page https://faculty.hampshire.edu/lspector/push3-description.html
    # ( 5 1.23 INTEGER.+ ( 4 ) INTEGER.- 5.67 FLOAT.* )
    interp = PushInterpreter()
    p = [5, 1.23, Add{Int}(), [4], Sub{Int}(), 5.67, Mul{Float64}()]
    execute!(interp, p)
    @test length(stack(interp, Int)) == 1
    @test last(stack(interp, Int)) == 1
    @test length(stack(interp, Float64)) == 1
    @test last(stack(interp, Float64)) == 6.9741
end # @testset "Example 2 from Push3 docs" begin

@testset "Example 3 from Push3 docs" begin
    # Example 3 from the page https://faculty.hampshire.edu/lspector/push3-description.html
    # ( 5 INTEGER.DUP INTEGER.+ )
    interp = PushInterpreter()
    p = [5, Dup{Int}(), Add{Int}()]
    execute!(interp, p)
    @test length(stack(interp, Int)) == 1
    @test last(stack(interp, Int)) == 10
end

end # @testset "Examples from PushGP docs" begin


text:
module ReadWriteLocks

export ReadWriteLock, read_lock, write_lock, lock!, unlock!

if isdefined(Base, :lock!)
    import Base: lock!
end

if isdefined(Base, :unlock!)
    import Base: unlock!
end

if !isdefined(Base, :Mutex)
    typealias Mutex ReentrantLock

    lock!(x::Mutex) = lock(x)
    unlock!(x::Mutex) = unlock(x)
end

abstract _Lock

immutable ReadLock{T<:_Lock}
    rwlock::T
end

immutable WriteLock{T<:_Lock}
    rwlock::T
end

type ReadWriteLock <: _Lock
    readers::Int
    writer::Bool
    lock::Mutex  # reentrant mutex
    condition::Condition
    read_lock::ReadLock
    write_lock::WriteLock

    function ReadWriteLock()
        rwlock = new(false, 0, Mutex(), Condition())
        rwlock.read_lock = ReadLock(rwlock)
        rwlock.write_lock = WriteLock(rwlock)

        return rwlock
    end
end

read_lock(rwlock::ReadWriteLock) = rwlock.read_lock
write_lock(rwlock::ReadWriteLock) = rwlock.write_lock

function lock!(read_lock::ReadLock)
    rwlock = read_lock.rwlock
    lock!(rwlock.lock)

    try
        while rwlock.writer
            wait(rwlock.condition)
        end

        rwlock.readers += 1
    finally
        unlock!(rwlock.lock)
    end

    return nothing
end

function unlock!(read_lock::ReadLock)
    rwlock = read_lock.rwlock
    lock!(rwlock.lock)

    try
        rwlock.readers -= 1
        if rwlock.readers == 0
            notify(rwlock.condition; all=true)
        end
    finally
        unlock!(rwlock.lock)
    end

    return nothing
end

function lock!(write_lock::WriteLock)
    rwlock = write_lock.rwlock
    lock!(rwlock.lock)

    try
        while rwlock.readers > 0 || rwlock.writer
            wait(rwlock.condition)
        end

        rwlock.writer = true
    finally
        unlock!(rwlock.lock)
    end

    return nothing
end

function unlock!(write_lock::WriteLock)
    rwlock = write_lock.rwlock
    lock!(rwlock.lock)

    try
        rwlock.writer = false
        notify(rwlock.condition; all=true)
    finally
        unlock!(rwlock.lock)
    end

    return nothing
end

end # module


text:
"""
Selenium 動作確認用サーバ
"""

using Genie, Genie.Router, Genie.Renderer.Json, Genie.Requests
using HTTP

route("/") do
    read("./public/index.html") |> String
end

# POST /api/server => {server: "apache"|"nginx", cache: [100..500:100]}
route("/api/server", method = POST) do
    server = postpayload(:server, "nginx")
    cache = postpayload(:cache, "100")
    (:server => server, :cache => cache) |> json
end

# Genie設定: http://localhost:8080
Genie.config.run_as_server = true
Genie.config.server_port = 8080
Genie.config.server_host = "0.0.0.0"
# Nginxのように静的ファイルを配信: ./public/*
# Genie.config.server_handle_static_files = true
# Genie.config.server_document_root = "./public/"

# Genie起動
Genie.startup()

max_stars_repo_name: amenoyoya/ssl-rpa

max_issues_repo_name: amenoyoya/ssl-rpa

max_forks_repo_name: amenoyoya/ssl-rpa


text:
module StreamStats

using OffsetArrays

function compute_shanon(data)
    """
    Calculate Shannon entropy value for a given byte array.
        Keyword arguments:
        data -- data bytes
    """
    entropy::Float16 = 0
    d = Array([(count(x->x==i,data)) for i in data])
    for i in d
        it = float(i/length(data))
        if i > 0
            entropy += - it * log(it ,2)
        end
    end
    return entropy
end

function compute_monte_carlo(data)
    """
    Calculate Monte Carlo Pi approximation for a given byte array.
    Keyword arguments:
    data -- data bytes
    """
    set_length::Int32 = (length(data) / 2) - 1
    r_square::Int16 = 128^2
    circle_surface::Float16 = 0

    for i in 1:set_length
        (data[i*2] - 128)
        if ((data[i*2] - 128) ^ 2 + (data[i*2+1] - 128) ^ 2) <= r_square
            circle_surface += 1
        end
    end
    return 4 * circle_surface / set_length
end

function get_pi_deviation(pi_value)
    """
    Returns an absolute percentage of difference between the provided Pi
    value and canonic.
    Keyword arguments:
    pi_value -- Pi value to be tested for difference
    """
    return abs(100 - (pi_value * 100 / pi))
end

function chi_squared(data)
    """Calculate Chi-Squared value for a given byte array.
    Keyword arguments:
    data -- data bytes
    """
    expected = length(data) / 256
    observed = OffsetVector([0 for i in 0:255], 0:255)

    for b in data
        observed[b] += 1
    end

    chi_squared = 0
    for o in observed
        chi_squared += (o - expected) ^ 2 / expected
    end

    return chi_squared
end

end


text:
# using Profile

function rand_poly(deg, vars)
    result = 0
    indices = vcat(collect(1:length(vars)), collect(1:length(vars)))
    monomials = []
    for d in 0:deg
        for subs in StructuralIdentifiability.IterTools.subsets(indices, d)
            push!(monomials, subs)
        end
    end

    for subs in monomials
        monom = rand(-50:50)
        for v_ind in subs
            monom *= vars[v_ind]
        end
        result += monom
    end

    return result
end

@testset "Power series solution for an ODE system" begin

    R, (x, x_dot) = Nemo.PolynomialRing(Nemo.QQ, ["x", "x_dot"])
    exp_t = ps_ode_solution([x_dot - x], Dict{fmpq_mpoly,fmpq}(x => 1), Dict{fmpq_mpoly,Array{fmpq,1}}(), 20)[x]
    @test valuation(ps_diff(exp_t) - exp_t) == 19

    R, (x, y, x_dot, y_dot, u) = Nemo.PolynomialRing(Nemo.QQ, ["x", "y", "x_dot", "y_dot", "u"])
    prec = 100
    eqs = [
        x_dot - x + 3 * x * y - u,
        y_dot + 2 * y - 4 * x * y
    ]
    u_coeff = [rand(1:5) for i in 1:prec]
    sol = ps_ode_solution(eqs, Dict(x => 0, y => -2), Dict(u => u_coeff), prec)
    @test map(e -> valuation(evaluate(e, [sol[v] for v in gens(R)])), eqs) == [prec - 1, prec - 1]

    F = Nemo.GF(2^31 - 1)
    prec = 100

    # Testing the function ps_ode_solution by itself
    for i in 1:30
        # Setting up the ring
        NUMX = 5
        NUMU = 3
        varnames = vcat(
            ["x_$(i)_dot" for i in 1:NUMX],
            ["x_$i" for i in 1:NUMX],
            ["u_$i" for i in 1:NUMU],
        )
        R, vars = Nemo.PolynomialRing(F, varnames)

        # Generating the initial conditions and inputs
        ic = Dict(vars[i + NUMX] => F(rand(-5:5)) for i in 1:NUMX)
        inputs = Dict(u => [F(rand(-3:3)) for i in 1:prec] for u in vars[(2 * NUMX + 1):end])
        # a dictionary for evaluation at zero (to avoid singularities)
        eval_at_zero = merge(ic, Dict(u => val[1] for (u, val) in inputs))

        # Generating denominators not vanishing at t = 0
        denominators = [rand_poly(1, vars[(NUMX + 1):end]) for i in 1:NUMX]
        while any([eval_at_dict(p, eval_at_zero) == 0 for p in denominators])
            denominators = [rand_poly(1, vars[(NUMX + 1):end]) for i in 1:NUMX]
        end

        eqs = [denominators[i] * vars[i] - rand_poly(2, vars[(NUMX + 1):end]) for i in 1:NUMX]
        @time sol = ps_ode_solution(eqs, ic, inputs, prec)
        evals = map(e -> valuation(evaluate(e, [sol[v] for v in gens(R)])), eqs)
        for e in evals
            @test e >= prec - 1
        end
    end

    # Testing ps_ode_solution in conjuntion with the ODE class
    for i in 1:30
        # Setting up the ring
        NUMX = 3
        NUMP = 3
        NUMU = 2
        varnames = vcat(
            ["x_$i" for i in 1:NUMX],
            ["p_$i" for i in 1:NUMP],
            ["u_$i" for i in 1:NUMU],
        )
        R, vars = Nemo.PolynomialRing(F, varnames)
        PType = gfp_mpoly
        TDict = Dict{PType,Union{PType,Generic.Frac{PType}}}

        # Generating the intial conditions etc
        ic = Dict(vars[i] => F(rand(-5:5)) for i in 1:NUMX)
        param_vals = Dict(vars[i + NUMX] => F(rand(-5:5)) for i in 1:NUMP)
        inputs = Dict(u => [F(rand(-3:3)) for i in 1:prec] for u in vars[(NUMX + NUMP + 1):end])
        eval_at_zero = merge(ic, param_vals, Dict(u => val[1] for (u, val) in inputs))

        # Generating denominators not vanishing at t = 0
        denominators = [rand_poly(1, vars) for i in 1:NUMX]
        while any([eval_at_dict(p, eval_at_zero) == 0 for p in denominators])
            denominators = [rand_poly(1, vars) for i in 1:NUMX]
        end

        eqs = TDict(vars[i] => rand_poly(2, vars) // denominators[i] for i in 1:NUMX)
        ode = ODE{PType}(eqs, TDict(), vars[(NUMX + NUMP + 1):end])
        @time sol = power_series_solution(ode, param_vals, ic, inputs, prec)
        ps_ring = parent(collect(values(sol))[1])
        for (p, val) in param_vals
            sol[p] = ps_ring(val)
        end
        eval_point = [sol[v] for v in vars]
        for (xd, f) in eqs
            num, den = unpack_fraction(f)
            lhs = divexact(evaluate(num, eval_point), evaluate(den, eval_point))
            rhs = ps_diff(sol[xd])
            @test valuation(lhs - rhs) >= prec - 1
        end
    end

end


text:
function sp_s(xg,Gf,phit,phib,Vsbstar)
    inv_phit1 = 1/phit
    oneSixth = 1/6.0
    invSqrt2 = 1/sqrt(2.0)
    xi           =  1.0 + Gf * invSqrt2
    inv_xi       =  1.0 /  xi
    Ux           =  Vsbstar * inv_phit1
    xn         =  phib * inv_phit1 + Ux
    if (xn < se)
        delta     =  exp(-xn)
    else
        delta     = ke / P3(xn - se)
    end
    margin       =  1.0e-5 * xi

    inv_xi = 1/xi
    Gf2 = Gf^2
    inv_Gf2 = 1/Gf2
    if (abs(xg) <= margin)
        SP_S_temp1 =  inv_xi * inv_xi * oneSixth * invSqrt2
        sp         =  xg * inv_xi * (1.0 + xg * (1.0 - (delta)) * Gf * SP_S_temp1)
    else
        if (xg < -margin)
            SP_S_yg     = -xg
            SP_S_ysub   = 1.25 * (SP_S_yg * inv_xi)
            SP_S_eta    = 0.5 * (SP_S_ysub + 10 - sqrt((SP_S_ysub - 6.0) * (SP_S_ysub - 6.0) + 64.0))
            SP_S_temp   = SP_S_yg - SP_S_eta
            SP_S_a      = SP_S_temp * SP_S_temp + Gf2*(SP_S_eta + 1.0)
            SP_S_c      = 2.0 * SP_S_temp - Gf2
            SP_S_tau    = -SP_S_eta + log(SP_S_a * inv_Gf2)
            SP_S_y0 = sigma(SP_S_a, SP_S_c, SP_S_tau, SP_S_eta)
            SP_S_delta0 = expl_high(SP_S_y0)
            SP_S_delta1 = 1.0 / SP_S_delta0
            SP_S_temp   = 1.0 / (2.0 + SP_S_y0 * SP_S_y0)
            SP_S_xi0    = SP_S_y0 * SP_S_y0 * SP_S_temp
            SP_S_xi1    = 4.0 * (SP_S_y0 * SP_S_temp * SP_S_temp)
            SP_S_xi2    = (8.0 * SP_S_temp - 12.0 * SP_S_xi0) * SP_S_temp * SP_S_temp
            SP_S_temp   = SP_S_yg - SP_S_y0
            SP_S_temp1  = (delta) * SP_S_delta1
            SP_S_pC     = 2.0 * SP_S_temp + Gf2 * (SP_S_delta0 - 1.0 - SP_S_temp1 + (delta) * (1.0 - SP_S_xi1))
            SP_S_qC     = SP_S_temp * SP_S_temp - Gf2 * (SP_S_delta0 - SP_S_y0 - 1.0 + SP_S_temp1 + (delta) * (SP_S_y0 - 1.0 - SP_S_xi0))
            SP_S_temp   = 2.0 - Gf2 * (SP_S_delta0 + SP_S_temp1 - (delta) * SP_S_xi2)
            SP_S_temp   = SP_S_pC * SP_S_pC - 2.0 * (SP_S_qC * SP_S_temp)
            sp          = -SP_S_y0 - 2.0 * (SP_S_qC / (SP_S_pC + sqrt(SP_S_temp)))
        else
            SP_xg1    = 1.0 / (1.25 + Gf * 7.324648775608221e-001)
            SP_S_A_fac= (xi * 1.25 * SP_xg1 - 1.0) * SP_xg1
            SP_S_xbar = xg * inv_xi * (1.0 + SP_S_A_fac * xg)
            SP_S_temp = expl_low(-SP_S_xbar)
            SP_S_w    = 1.0 - SP_S_temp
            SP_S_x1   = xg + Gf2 * 0.5 - Gf * sqrt(xg + Gf2 * 0.25 - SP_S_w)
            SP_S_bx   = (xn) + 3.0
            SP_S_eta  = MINA(SP_S_x1, SP_S_bx, 5.0) - 0.5 * (SP_S_bx - sqrt(SP_S_bx * SP_S_bx + 5.0))
            SP_S_temp = xg - SP_S_eta
            SP_S_temp1= exp(-SP_S_eta)
            SP_S_temp2= 1.0 / (2.0 + SP_S_eta * SP_S_eta)
            SP_S_xi0  = SP_S_eta * SP_S_eta * SP_S_temp2
            SP_S_xi1  = 4.0 * (SP_S_eta * SP_S_temp2 * SP_S_temp2)
            SP_S_xi2  = (8.0 * SP_S_temp2 - 12.0 * SP_S_xi0) * SP_S_temp2 * SP_S_temp2
            SP_S_a    = max(1.0e-40, SP_S_temp * SP_S_temp - Gf2 * (SP_S_temp1 + SP_S_eta - 1.0 - (delta) * (SP_S_eta + 1.0 + SP_S_xi0)))
            SP_S_b    = 1.0 - 0.5 * (Gf2 * (SP_S_temp1 - (delta) * SP_S_xi2))
            SP_S_c    = 2.0 * SP_S_temp + Gf2 * (1.0 - SP_S_temp1 - (delta) * (1.0 + SP_S_xi1))
            SP_S_tau  = (xn) - SP_S_eta + log(SP_S_a / Gf2)
            SP_S_x0 = sigma2(SP_S_a, SP_S_b, SP_S_c, SP_S_tau, SP_S_eta)
            if (SP_S_x0 < se05)
                SP_S_delta0 = exp(SP_S_x0)
                SP_S_delta1 = 1.0 / SP_S_delta0
                SP_S_delta0 = (delta) * SP_S_delta0
            else
                if (SP_S_x0 > (xn) - se05)
                    SP_S_delta0 = exp(SP_S_x0 - (xn))
                    SP_S_delta1 = (delta) / SP_S_delta0
                else
                    SP_S_delta0 = ke05 / P3((xn) - SP_S_x0 - se05)
                    SP_S_delta1 = ke05 / P3(SP_S_x0 - se05)
                end
            end
            SP_S_temp   = 1.0 / (2.0 + SP_S_x0 * SP_S_x0)
            SP_S_xi0    = SP_S_x0 * SP_S_x0 * SP_S_temp
            SP_S_xi1    = 4.0 * (SP_S_x0 * SP_S_temp * SP_S_temp)
            SP_S_xi2    = (8.0 * SP_S_temp - 12.0 * SP_S_xi0) * SP_S_temp * SP_S_temp
            SP_S_temp   = xg - SP_S_x0
            SP_S_pC     = 2.0 * SP_S_temp + Gf2 * (1.0 - SP_S_delta1 + SP_S_delta0 - (delta) * (1.0 + SP_S_xi1))
            SP_S_qC     = SP_S_temp * SP_S_temp - Gf2 * (SP_S_delta1 + SP_S_x0 - 1.0 + SP_S_delta0 - (delta) * (SP_S_x0 + 1.0 + SP_S_xi0))
            SP_S_temp   = 2.0 - Gf2 * (SP_S_delta1 + SP_S_delta0 - (delta) * SP_S_xi2)
            SP_S_temp   = SP_S_pC * SP_S_pC - 2.0 * (SP_S_qC * SP_S_temp)
            sp          = SP_S_x0 + 2.0 * (SP_S_qC / (SP_S_pC + sqrt(SP_S_temp)))
        end
    end
end

P3(u) =  (1.0 + (u) * (1.0 + 0.5 * ((u) * (1.0 + (u)/3.0))))

function expl_high(x)
    if ((x) < se05)
            res       = exp(x)
    else
            res       =  ke05inv * P3((x) - se05)
    end
end

function expl_low(x)
    if ((x) > -se05)
            res       =  exp(x)
    else
            res       = ke05 / P3(-se05 - (x))
    end
    return res
end

function sigma2(a,b,c,tau,eta)
    nu           =  (a) + (c)
    mutau        =  (nu) * (nu) + (tau) * (0.5 * ((c) * (c)) - (a) * (b))
    y            =  (eta) + (a) * nu * (tau) / (mutau + (nu / mutau) * (tau) * (tau) * (c) * ((c) * (c) * 1/3.0 - (a) * (b)))
    return y
end

function sigma(a,c,tau,eta)
    nu           =  (a) + (c)
    mutau        =  nu * nu + (tau) * (0.5 * ((c) * (c)) - (a))
    y            =  (eta) + (a) * nu * (tau) / (mutau + (nu / mutau) * (tau) * (tau) * (c) * ((c) * (c) * 1/3.0 - (a)))
    return y
end

function Vp_EKV(Vg,Vfb,Na,tox,T,phif,sigma =0.0 )
    GAMMA = γ(Na,tox)
    PHI = 2*phif
    VTO = Vfb+PHI+GAMMA*sqrt(PHI)
    VGprime = Vg - VTO + PHI + GAMMA * sqrt(PHI)
    VP = VGprime - PHI - GAMMA* (sqrt(VGprime+(GAMMA/2.0)*(GAMMA/2.0))-(GAMMA/2.0));
    return VP
end

function n_EKV(Vg,Vfb,Na,tox,T,phif,sigma =0.0 )
    GAMMA = γ(Na,tox)
    PHI = 2*phif
    VP = Vp_EKV(Vg,Vfb,Na,tox,T,phif)
    n = 1.0 + GAMMA / (2.0*sqrt(PHI + VP + 4.0*$vt));
    return n
end

function q_EKV(Vg,V,Vfb,Na,tox,T,phif,sigma =0.0 )
    VP = Vp_EKV(Vg,Vfb,Na,tox,T,phif)
    x=(VP-VS)/kb/T
    iff = (ln(1.0+exp( x /2.0)))*(ln(1.0+exp( x /2.0)));
    return iff
end

function id_EKV(Vg,Vd,Vs,Vfb,Na,tox,T,phif,sigma =0.0 )
    pt = kb*T
    Ispec = 2*n_EKV(Vg,Vfb,Na,tox,T,phif)*pt*pt
    return Ispec*(q_EKV(Vg,Vs,Vfb,Na,tox,T,phif,sigma)-q_EKV(Vg,Vd,Vfb,Na,tox,T,phif,sigma))

end


text:
"""
        greedy_star2_coloring

    Find a coloring of a given input graph such that
    no two vertices connected by an edge have the same
    color using greedy approach. The number of colors
    used may be equal or greater than the chromatic
    number `χ(G)` of the graph.

    A star coloring is a special type of distance - 1  coloring,
    For a coloring to be called a star coloring, it must satisfy
    two conditions:

    1. every pair of adjacent vertices receives distinct  colors
    (a distance-1 coloring)

    2. For any vertex v, any color that leads to a two-colored path
    involving v and three other vertices  is  impermissible  for  v.
    In other words, every path on four vertices uses at least three
    colors.

    Reference: Gebremedhin AH, Manne F, Pothen A. **What color is your Jacobian? Graph coloring for computing derivatives.** SIAM review. 2005;47(4):629-705.

    TODO: add text explaining the difference between star1 and
    star2
"""
function color_graph(g::LightGraphs.AbstractGraph, :: GreedyStar2Color)
    v = nv(g)
    colorvec = zeros(Int, v)

    forbidden_colors = zeros(Int, v+1)

    for vertex_i = vertices(g)

        for w in inneighbors(g, vertex_i)
            if colorvec[w] != 0
                forbidden_colors[colorvec[w]] = vertex_i
            end

            for x in inneighbors(g, w)
                if colorvec[x] != 0
                    if colorvec[w] == 0
                        forbidden_colors[colorvec[x]] = vertex_i
                    else
                        if colorvec[x] < colorvec[w]
                            forbidden_colors[colorvec[x]] = vertex_i
                        end
                    end
                end
            end
        end

        colorvec[vertex_i] = find_min_color(forbidden_colors, vertex_i)
    end

    colorvec
end


text:
# Tests for the random subspace ensemble
function t_randomsubspace()

	## Test the random subspace ensemble on toy data

	X = [1 2 3; -1 -2 3; 0 0 3] # data
	y = [1,1,0]		    # labels

	f_exec(model,X) = X.+model # the execution function adds to the input the value contained in the model

	# Create an ensemble that for each variable (i.e. line) of X, takes the maximum, which is '3' in all cases. 'y' is not used as
	# the maximum function only takes the first argument of the input `x` (`x` is represented in the training function by Tuple(X,y)
	ensemble = j4pr.RandomSubspace.randomsubspace_train(X, y, 5, 1, x->maximum(x[1],2), f_exec, j4pr.ClassifierCombiner.NoCombiner(), true)

	# Test that the `trained` ensemble members are matrices of a single element of value 3
	Test.@test all((ensemble.members[i] ==Matrix([3]) for i in eachindex(ensemble.members)))

	# Run the ensemble execution (we expect that for the lines specified by ensemble.idx the corresponding values form ensemble.members are added)
	Xt = [1 1 1;2 2 2;3 3 3]
	ensemble_results = j4pr.RandomSubspace.randomsubspace_exec(ensemble,Xt)

	Test.@test all((ensemble_results[i] == f_exec(ensemble.members[i], Xt[ensemble.idx[i],:]) for i in eachindex(ensemble.members)))

	## Test the random subspace ensemble on the fish dataset and a classifier (the test is just a functionality assertion, does not verify values)
	A = j4pr.DataGenerator.fish(20) 			# load data
	tr,ts = j4pr.splitobs(j4pr.shuffleobs(A), 0.5)		# shuffle and split data
	w=j4pr.parzen(0.1)

	C = 2; # number of classes (for the 'fish dataset')
	L = 5; # ensemble size
	M = 2; # variables per ensemble memeber

	# Loop through combiners (must be non-trainable), train subspace on `tr` and run it on `ts`
	for comb in [j4pr.ClassifierCombiner.GeneralizedMeanCombiner(L,C,1.0), j4pr.ClassifierCombiner.ProductCombiner(L,C)]
		Test.@test try
			wes=j4pr.randomsubspace(w, L, M, comb, parallel_execution=false) # ensemble of 5 members, 2 variables, serial execution
			wep=j4pr.randomsubspace(w, L, M, comb, parallel_execution=true)  # ensemble of 5 members, 2 variables, parallel execution

			# Train
			west = wes(tr)
			wept = wep(tr)

			# Execute
			west(ts);
			west(ts);

			# If we're here, everything worked
			true
		catch
			false
		end

	end

	# Test variable sampling without replacing
	Test.@test try
		wes2=j4pr.randomsubspace(w, 2, 1, j4pr.ClassifierCombiner.NoCombiner(),false, parallel_execution=false)
		west2 = wes2(tr)
		west2(+ts);
		true
	catch
		false
	end

end


text:
ms = map(measures()) do m
    m.name
end
@test "cross_entropy" in ms
@test "rms"  in ms

ms = map(measures(matching(categorical(1:3)))) do m
         m.name
end
@test "cross_entropy" in ms
@test !("rms"  in ms)

ms = map(measures(matching(rand(3)))) do m
         m.name
end
@test !("cross_entropy" in ms)
@test "rms"  in ms

true


text:
export ncross, grow, center

Polygon(points::Vararg{T, N}) where {T <: Point, N} = Polygon([points...])

needed_attributes(::Type{Polygon}) = (Visible, Linewidth, Stroke, Linestyle, Fill)

Base.convert(::Type{Polygon}, ps::Vector{Point}) = Polygon(ps)

function ncross(pos::Point, n::Int, r::Real, thickness::Real, angle::Angle=deg(0))
    if n < 3
        error("ncross needs at least three arms")
    end
    angles = deg.(range(0, 360, length=n+1)[1:end-1]) .+ angle
    routcorners = rotate.(Point(r, -0.5*thickness * r), angles)
    loutcorners = rotate.(Point(r,  0.5*thickness * r), angles)
    inner_point1 = intersection(
        Line(loutcorners[1], Point(0, loutcorners[1].y)),
        Line(routcorners[2], routcorners[2] - Point(angles[2]))
    )
    inner_points = rotate.(inner_point1, angles)
    collect(Iterators.flatten((r, l, i) for (r, l, i) in zip(routcorners, loutcorners, inner_points))) .+ pos
end

# function nstar(pos::Point)

function center(pol::Polygon)
    Point(sum([p.xy for p in pol.points]) / length(pol.points))
end

function grow(p::Polygon, factor::Real; from::Union{Point, Nothing}=nothing)
    c = isnothing(from) ? center(p) : from
    dists = from_to.(c, p.points)
    Polygon(dists .* factor .+ c)
end


text:
## activate package environment
cd(@__DIR__)
using Pkg
Pkg.activate("Project.toml")

## parameters (calculation for perpendicular incident)
λ = 520 # nm
N = [1.0; 3.812 + 0.1903im; 4.1891+0.3633im] # air - window layer - absorber layer
D = 100 # thickness of window layer
# (refractive indices are from https://refractiveindex.info/, data by Aspnes 1986)

## function for fresnel coefficients (for perpendicular incident)
function fresnel_coeff(n,m,pol)
    if pol == "s"
        r = (n-m)/(n+m)
        t = (2*n)/(n+m)
    elseif pol == "p"
        r = (m-n)/(m+n)
        t = (2*n)/(m+n)
    end
    return r,t
end

## function for interface propagation matrix
function interface_prop(n,m,pol)
    r,t = fresnel_coeff(n,m,pol)
    I = 1/t * [1 r;r 1]
    return I
end

## function for propagation through film Matrix
function film_prop(d,n)
    beta = 2*pi*d*n/λ
    L = [exp(-1im*beta) 0; 0 exp(1im*beta)]
    return L
end

## calcuation of scattering matrix
# Emanuele Centurioni: Generalized matrix method for calculation of internal light energy flux
# in mixed coherent and incoherent multilayers
Ss = interface_prop(N[1],N[2],"s") * film_prop(D,N[2]) * interface_prop(N[2],N[3],"s")
Sp = interface_prop(N[1],N[2],"p") * film_prop(D,N[2]) * interface_prop(N[2],N[3],"p")

## backside transmittance
Ts = abs2(1/Ss[1,1])*real(N[3])/real(N[1])
Tp = abs2(1/Ss[1,1])*real(conj(N[3]))/real(N[1])
T  = round((Ts+Tp)/2,sigdigits=4)

## ouput result
println("")
println("transmission into absorber T = $T")

max_stars_repo_name: MaxGrossmann/TRPL

max_issues_repo_name: MaxGrossmann/TRPL

max_forks_repo_name: MaxGrossmann/TRPL


text:
# Segment and survey types

"""
`SegmentType`: Type of a segment in a Bottom Hole Assembly (BHA), currently defined:\n
  `Collar`, `Pipe`, `Stabilizer`, `Bit`
"""
abstract type SegmentType end

"""
`Collar(material, length, id, od, fc)`: Definition of a Collar segment.

where:\n
  `material`::Symbol, e.g. :Steel or :Composite\n
  `length`  ::Float64, length of segment, in [m] or [ft]\n
  `id`      ::Float74, internal diameter, in [m] or [inch]\n
  `od`      ::Float74, external diameter, in [m] or [inch]\n
  `fc`      ::Float74, friction coefficient, in [-]
"""
mutable struct Collar <: SegmentType
  material::Symbol
  length::Float64
  id::Float64
  od::Float64
  fc::Float64
end

"""
`Pipe(material, length, id, od, fc)`: Definition of a Pipe segment.

where:\n
  `material`::Symbol, e.g. :Steel or :Composite\n
  `length`  ::Float64, length of segment, in [m] or [ft]\n
  `id`      ::Float74, internal diameter, in [m] or [inch]\n
  `od`      ::Float74, external diameter, in [m] or [inch]\n
  `fc`      ::Float74, friction coefficient, in [-]
"""
mutable struct Pipe <: SegmentType
  material::Symbol
  length::Float64
  id::Float64
  od::Float64
  fc::Float64
end

"""
`Stabilizer(material, id, od, fc)`: Definition of a Stabilizer segment.

where:\n
  `material`::Symbol, e.g. :Steel or :Composite\n
  `length`  ::Float64, length of segment, in [m] or [ft]\n
  `id`      ::Float74, internal diameter, in [m] or [inch]\n
  `od`      ::Float74, external diameter, in [m] or [inch]\n
  `fc`      ::Float74, friction coefficient, in [-]
"""
mutable struct Stabilizer <: SegmentType
  material::Symbol
  length::Float64
  id::Float64
  od::Float64
  fc::Float64
end

"""
`Rig(id, od)`: Definition of a Rig segment. This element must
be the topmost element of the BHA and will fix the position of
the BHA in the center.

where:\n
  `material`::Symbol, e.g. :Steel or :Composite\n
  `length`  ::Float64, length of segment, in [m] or [ft]\n
  `id`      ::Float74, internal diameter, in [m] or [inch]\n
  `od`      ::Float74, external diameter, in [m] or [inch]\n
  `fc`      ::Float74, friction coefficient, in [-]
"""
mutable struct Rig <: SegmentType
  material::Symbol
  length::Float64
  id::Float64
  od::Float64
  fc::Float64
end

"""
`Bit(material, id, od, fc)`: Definition of a Bit segment.

where:\n
  `material`::Symbol, e.g. :Steel or :Composite\n
  `length`  ::Float64, length of segment, in [m] or [ft]\n
  `id`      ::Float74, internal diameter, in [m] or [inch]\n
  `od`      ::Float74, external diameter, in [m] or [inch]\n
  `fc`      ::Float74, friction coefficient, in [-]
"""
mutable struct Bit <: SegmentType
  material::Symbol
  length::Float64
  id::Float64
  od::Float64
  fc::Float64
end
"""
`Bentsub(material, id, od, fc, angle)`: Definition of the Bentsub.

where:\n
  `angle`   ::Bentsub angle [degrees]
"""
mutable struct Bentsub <: SegmentType
  material::Symbol
  angle::Float64
end

"""
# ExtForce

ExtForce: Segment with no length that handles an external force applied to the assemby.

The forces fx and fy are assumed to be in the local coordinate system.

### Tpe
```julia
ExtForce(depth, inclination, heading, Fx, Fy)
```
### Arguments
```julia
* `fx`          : Force in x direction
* `fy`          : Force in y direction
```
"""
mutable struct ExtForce <: SegmentType
  fx::Float64
  fy::Float64
end

# Segment and survey types

"""
`SurveyType`: Type of a survey record, currently defined:\n
  `RigSurvey`, `HoleSurvey`, `BitSurvey`
"""
abstract type SurveyType end
mutable struct RigSurvey <: SurveyType
  depth::Float64
  inclination::Float64          # [Degrees]
  heading::Float64              # [Degrees]
  diameter::Float64
end

mutable struct HoleSurvey <: SurveyType
  depth::Float64
  inclination::Float64          # [Degrees]
  heading::Float64              # [Degrees]
  diameter::Float64
end

mutable struct BitSurvey <: SurveyType
  depth::Float64
  inclination::Float64          # [Degrees]
  heading::Float64              # [Degrees]
  diameter::Float64
end

mutable struct TpRecord
  dalpha::Float64
  wob::Float64
  inclination::Float64
  fbit::Float64
  firsttouch::Float64
  f_firsttouch::Float64
  solution::DataFrame
end

function tprecord_show(io::IO, tprecord::TpRecord, compact::Bool)
  println(io, Float64[getfield(tprecord, f) for f in fieldnames(tprecord)[1:6]])
end

show(io::IO, tprecord::TpRecord) = tprecord_show(io, tprecord, false)
showcompact(io::IO, tprecord::TpRecord) = tprecord_show(io, tprecord, true)

# Trajectory type

"""
`Trajectory(inclination, heading, diameter)`: Definition of the well trajectory.

where:\n
  `inclination` ::Float64, inclination from vertical, in [radians]\n
  `heading`     ::Float64, direction from North, in [radians]\n
  `diameter`    ::Float74, internal diameter, in [inch]\n
"""
mutable struct Trajectory
  heading::Float64          # Used, but irrelevant
  diameter::Float64
end
Trajectory() = Trajectory(0.0, 0.0)

# BHA type

"""
`BHA(material, length, id, od, fc)`: Definition of a BHA type.

where:\n
  `ratio`       ::Float64, used to generate the FE mesh, in []\n
  `medium`      ::Symbol, currently only :mud is supported\n
  `maxiters`    ::Int, maximum number of FE up/down iterations\n
  `dlerrormax`  ::Float74, maximum dl difference between 2 iterations\n
  `drerrormax`  ::Float74, maximum dr difference between 2 iterations\n
  `segments`    ::Vector of Segments (see SegmentType)\n
  `surveys`     ::Vector of Survey records (see SurveyType)\n
"""
mutable struct BHAJ{T}
  ratio::Float64
  medium::Symbol
  maxiters::Int
  dlerrormax::Float64
  drerrormax::Float64
  segments::Array{SegmentType, 1}
  surveys::Array{SurveyType, 1}
  trajectory::Trajectory
  tprecords::Array{TpRecord, 1}
  name::T
  pdir::String
  debuglevel::Int
end

# External constructors
BHAJ(name::T, pdir::String, debuglevel::Int=0) where T = BHAJ(
  0.6,
  :mud,
  Int(1e5),
  1e-7,
  1e-7,
  Array{SegmentType, 1}(),
  Array{SurveyType, 1}(),
  Trajectory(),
  Array{TpRecord, 1}(),
  name,
  pdir,
  debuglevel
)

"""
# Finite Element object

Element: Finite element object

### Tpe
```julia
Fe(length, id, od, radius, fc, inclination, heading, mass, ea, eix, eiy, gl)
```
### Arguments
```julia
* `length`      : Length of element
* `id`          : Inside diameter
* `od`          : Outside diameter
* `radius`      : Radius
* `fc`          : Friction coefficient
* `inclination` : Inclination [Radians]
* `heading`     : Heading [Radians]
* `mass`        : Element mass
                  Young moduli
* `ea`          : E*A in z direction
* `er`          : E*I in r direction
* `gj`          : Shear modulus
* `nu`          : Poisson constant
```
"""
mutable struct Element
  length::Float64
  id::Float64                 # Inside diameter
  od::Float64                 # Outside diameter
  radius::Float64
  fc::Float64                 # Friction coefficient
  inclination::Float64
  heading::Float64
  mass::Float64
  ea::Float64
  er::Float64
  gj::Float64
  nu::Float64
end
Element() = Element(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

"""
# Node object

Node: Node object

### Constructor
```julia
Node(z, x, y, tvd, xvert, dispew, dispns, inclination, deltainclination, dip,
  heading, deltaheading, holeradius, stringradius, fc, clearance,
  weight, mass, massinertia, axialinertia, radialinertia)

* dip  : Angle [Radians] from horizon
* inclination : Angle [Radians] from vertical
```
### Fields
```julia
* `x`                 : x coordinate
* `y`                 : y coordinate
* `z`                 : z coordinate
* `tvd`               : Measured depth
* `xvert`             : True vertical depth
* `dispew`            : Angle from east to west [Radians]
* `dispns`            : Angle from north to south [Radians]
* `inclination`       : Inclination [Radians]
* `deltainclination`  : Add'l inclinationination angle (bentsub) [Radians]
* `dip`               : Dip at node [Radians]
* `heading`           : Heading [Radians]
* `deltaheading`         : Add'l heading angle (bentsub) [Radians]
* `holeradius`        : Hole radius
* `stringradius`      : String radius
* `fc`                : Friction coefficient
* `clearance`         : Clearance element in hole
* `weight`            : Weight (force)
* `mass`              : Mass
* `massinertia`       : Second moment of inertia
* `axialinertia`      : Moment of inertia
* `radialinertia`     : Radial inertia
* `externalxforce`    : External force in x direction applied at node
* `externalyforce`    : External force in y direction applied at node
* `fixedxcoord`       : Fixed location in x direction at node
* `fixedycoord`       : Fixed location in y direction at node
* `fixedxtheta`       : Fixed angle from z towards x direction at node
* `fixedytheta`       : Fixed angle from z towards y direction at node
```
"""
mutable struct Node
  x::Float64
  y::Float64
  z::Float64
  tvd::Float64
  xvert::Float64
  dispew::Float64
  dispns::Float64
  inclination::Float64
  deltainclination::Float64
  dip::Float64
  heading::Float64
  deltaheading::Float64
  holeradius::Float64
  stringradius::Float64
  fc::Float64
  clearance::Float64
  weight::Float64
  mass::Float64
  massinertia::Float64
  axialinertia::Float64
  radialinertia::Float64
  externalxforce::Float64
  externalyforce::Float64
  fixedxcoord::Float64
  fixedycoord::Float64
  fixedxtheta::Float64
  fixedytheta::Float64
end
Node() = Node(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)


text:
module TestUniqueID

import Historic
using Test

using ..Utils

function test_single_thread()
    ids = [Historic.uniqueid() for _ in 1:100]
    @test allunique(ids)
end

function test_multi_thread()
    outputs = [UInt64[] for _ in 1:Threads.nthreads()]
    Utils.foreachthread(outputs) do ids
        for _ in 1:1000
            push!(ids, Historic.uniqueid())
        end
    end
    ids = reduce(vcat, outputs)
    @test allunique(ids)
end

end  # module