text:
"A result from a PostgreSQL database query"
mutable struct Result
    "A pointer to a libpq PGresult object (C_NULL if cleared)"
    result::Ptr{libpq_c.PGresult}

    "True if the PG object has been cleaned up"
    closed::Atomic{Bool}

    "PostgreSQL Oids for each column in the result"
    column_oids::Vector{Oid}

    "Julia types for each column in the result"
    column_types::Vector{Type}

    "Whether to expect NULL for each column (whether output data can have `missing`)"
    not_null::Vector{Bool}

    "Conversions from PostgreSQL data to Julia types for each column in the result"
    column_funcs::Vector{Base.Callable}

    "Name of each column in the result"
    column_names::Vector{String}

    # TODO: attach encoding per https://wiki.postgresql.org/wiki/Driver_development#Result_object_and_client_encoding
    function Result(
        result::Ptr{libpq_c.PGresult},
        jl_conn::Connection;
        column_types::Union{AbstractDict, AbstractVector}=ColumnTypeMap(),
        type_map::AbstractDict=PQTypeMap(),
        conversions::AbstractDict=PQConversions(),
        not_null=false,
    )
        jl_result = new(result, Atomic{Bool}(result == C_NULL))

        type_lookup = LayerDict(
            PQTypeMap(type_map),
            jl_conn.type_map,
            LIBPQ_TYPE_MAP,
            _DEFAULT_TYPE_MAP,
        )

        func_lookup = LayerDict(
            PQConversions(conversions),
            jl_conn.func_map,
            LIBPQ_CONVERSIONS,
            _DEFAULT_CONVERSIONS,
            _FALLBACK_CONVERSION,
        )

        jl_result.column_oids = col_oids = map(1:num_columns(jl_result)) do col_num
            libpq_c.PQftype(jl_result.result, col_num - 1)
        end

        jl_result.column_names = map(1:num_columns(jl_result)) do col_num
            unsafe_string(libpq_c.PQfname(jl_result.result, col_num - 1))
        end

        column_type_map = ColumnTypeMap()
        for (k, v) in pairs(column_types)
            column_type_map[column_number(jl_result, k)] = v
        end

        jl_result.column_types = col_types = collect(Type, imap(enumerate(col_oids)) do itr
            col_num, col_oid = itr
            get(column_type_map, col_num) do
                get(type_lookup, col_oid, String)
            end
        end)

        jl_result.column_funcs = collect(Base.Callable, imap(col_oids, col_types) do oid, typ
            func_lookup[(oid, typ)]
        end)

        # figure out which columns the user says may contain nulls
        if not_null isa Bool
            jl_result.not_null = fill(not_null, size(col_types))
        elseif not_null isa AbstractArray
            if eltype(not_null) === Bool
                if length(not_null) != length(col_types)
                    throw(ArgumentError(
                        "The length of keyword argument not_null, when an array, must be equal to the number of columns"
                    ))
                end

                jl_result.not_null = not_null
            else
                # assume array of column names
                jl_result.not_null = fill(false, size(col_types))

                for col_name in not_null
                    col_num = column_number(jl_result, col_name)
                    jl_result.not_null[col_num] = true
                end
            end
        else
            throw(ArgumentError(
                "Unsupported type $(typeof(not_null)) for keyword argument not_null"
            ))
        end

        finalizer(close, jl_result)

        return jl_result
    end
end

"""
    show(io::IO, jl_result::Result)

Show a PostgreSQL result and whether it has been cleared.
"""
function Base.show(io::IO, jl_result::Result)
    print(io, "PostgreSQL result")

    if !isopen(jl_result)
        print(io, " (cleared)")
    end
end

"""
    status(jl_result::Result) -> libpq_c.ExecStatusType

Return the status of a result's corresponding database query according to libpq.

See also: [`error_message`](@ref)
"""
status(jl_result::Result) = libpq_c.PQresultStatus(jl_result.result)

"""
    error_message(jl_result::Result; verbose=false) -> String

Return the error message associated with the result, or an empty string if there was no
error.
If `verbose`, have libpq generate a more verbose version of the error message if possible.
Includes a trailing newline.
"""
function error_message(jl_result::Result; verbose=false)
    return verbose ? _verbose_error_message(jl_result) : _error_message(jl_result)
end

function _error_message(jl_result::Result)
    return unsafe_string(libpq_c.PQresultErrorMessage(jl_result.result))
end

function _verbose_error_message(jl_result::Result)
    msg_ptr = libpq_c.PQresultVerboseErrorMessage(
        jl_result.result,
        libpq_c.PQERRORS_VERBOSE,
        libpq_c.PQSHOW_CONTEXT_ALWAYS,
    )

    if msg_ptr == C_NULL
        error(LOGGER, Errors.JLResultError(
            "libpq could not allocate memory for the result error message"
        ))
    end

    msg = unsafe_string(msg_ptr)
    libpq_c.PQfreemem(msg_ptr)
    return msg
end

"""
    error_field(jl_result::Result, field_code::Char) -> Union{String, Nothing}

Get an individual field from the error report in a [`Result`](@ref).
Returns `nothing` if that field is not provided for this error, or if there is no error or
warning in this `Result`.

See [](https://www.postgresql.org/docs/10/libpq-exec.html#LIBPQ-PQRESULTERRORFIELD)
for all available fields.

## Example

```
julia> LibPQ.error_field(result, LibPQ.libpq_c.PG_DIAG_SEVERITY)
"ERROR"
```
"""
function error_field(jl_result::Result, field_code::Union{Char, Integer})
    ret = libpq_c.PQresultErrorField(jl_result.result, field_code)
    return ret == C_NULL ? nothing : unsafe_string(ret)
end

"""
    close(jl_result::Result)

Clean up the memory used by the `PGresult` object.
The `Result` will no longer be usable.
"""
function Base.close(jl_result::Result)
    if !atomic_cas!(jl_result.closed, false, true)
        ptr, jl_result.result = jl_result.result, C_NULL
        libpq_c.PQclear(ptr)
    end
    return nothing
end

"""
    isopen(jl_result::Result)

Determine whether the given `Result` has been `close`d, i.e. whether the memory
associated with the underlying `PGresult` object has been cleared.
"""
Base.isopen(jl_result::Result) = !jl_result.closed[]

"""
    handle_result(jl_result::Result; throw_error::Bool=true) -> Result

Check status and handle errors for newly-created result objects.

If `throw_error` is `true`, throw an error and clear the result if the query results in a
fatal error or unreadable response.
Otherwise a warning is shown.

Also print an info message about the result.
"""
function handle_result(jl_result::Result; throw_error::Bool=true)
    result_status = status(jl_result)

    if result_status in (libpq_c.PGRES_BAD_RESPONSE, libpq_c.PGRES_FATAL_ERROR)
        err = Errors.PQResultError(jl_result)

        if throw_error
            close(jl_result)
            error(LOGGER, err)
        else
            warn(LOGGER, err)
        end
    else
        if result_status == libpq_c.PGRES_NONFATAL_ERROR
            warn(LOGGER, Errors.PQResultError(jl_result))
        end

        debug(LOGGER, unsafe_string(libpq_c.PQcmdStatus(jl_result.result)))
    end

    return jl_result
end

"""
    execute(
        {jl_conn::Connection, query::AbstractString | stmt::Statement},
        [parameters::Union{AbstractVector, Tuple},]
        throw_error::Bool=true,
        column_types::AbstractDict=ColumnTypeMap(),
        type_map::AbstractDict=LibPQ.PQTypeMap(),
        conversions::AbstractDict=LibPQ.PQConversions(),
    ) -> Result

Run a query on the PostgreSQL database and return a `Result`.
If `throw_error` is `true`, throw an error and clear the result if the query results in a
fatal error or unreadable response.

The query may be passed as `Connection` and `AbstractString` (SQL) arguments, or as a
`Statement`.

`execute` optionally takes a `parameters` vector which passes query parameters as strings to
PostgreSQL.

`column_types` accepts type overrides for columns in the result which take priority over
those in `type_map`.
For information on the `column_types`, `type_map`, and `conversions` arguments, see
[Type Conversions](@ref typeconv).
"""
function execute end

function execute(
    jl_conn::Connection,
    query::AbstractString;
    throw_error::Bool=true,
    kwargs...
)
    result = lock(jl_conn) do
        _execute(jl_conn.conn, query)
    end

    return handle_result(Result(result, jl_conn; kwargs...); throw_error=throw_error)
end

function execute(
    jl_conn::Connection,
    query::AbstractString,
    parameters::Union{AbstractVector, Tuple};
    throw_error::Bool=true,
    kwargs...
)
    string_params = string_parameters(parameters)
    pointer_params = parameter_pointers(string_params)

    result = lock(jl_conn) do
        _execute(jl_conn.conn, query, pointer_params)
    end

    return handle_result(Result(result, jl_conn; kwargs...); throw_error=throw_error)
end

function _execute(conn_ptr::Ptr{libpq_c.PGconn}, query::AbstractString)
    return libpq_c.PQexec(conn_ptr, query)
end

function _execute(
    conn_ptr::Ptr{libpq_c.PGconn},
    query::AbstractString,
    parameters::Vector{Ptr{UInt8}},
)
    num_params = length(parameters)

    return libpq_c.PQexecParams(
        conn_ptr,
        query,
        num_params,
        C_NULL,  # set paramTypes to C_NULL to have the server infer a type
        parameters,
        C_NULL,  # paramLengths is ignored for text format parameters
        zeros(Cint, num_params),  # all parameters in text format
        zero(Cint),  # return result in text format
    )
end

"""
    string_parameters(parameters::AbstractVector) -> Vector{Union{String, Missing}}

Convert parameters to strings which can be passed to libpq, propagating `missing`.
"""
function string_parameters end

string_parameters(parameters::AbstractVector{<:Parameter}) = parameters

# Tuples of parameters
string_parameters(parameters::Tuple) = string_parameters(collect(parameters))

# vector which can't contain missing
string_parameters(parameters::AbstractVector) = map(string_parameter, parameters)

# vector which might contain missings
function string_parameters(parameters::AbstractVector{>:Missing})
    collect(
        Union{String, Missing},
        imap(parameters) do parameter
            ismissing(parameter) ? missing : string_parameter(parameter)
        end
    )
end

string_parameter(parameter) = string(parameter)

function string_parameter(parameter::AbstractVector)
    io = IOBuffer()
    print(io, "{")
    join(io, (_array_element(el) for el in parameter), ",")
    print(io, "}")
    String(take!(io))
end

_array_element(el::AbstractString) = "\"$el\""
_array_element(el::Missing) = "NULL"
_array_element(el) = string_parameter(el)

"""
    parameter_pointers(parameters::AbstractVector{<:Parameter}) -> Vector{Ptr{UInt8}}

Given a vector of parameters, returns a vector of pointers to either the string bytes in the
original or `C_NULL` if the element is `missing`.
"""
function parameter_pointers(parameters::AbstractVector{<:Parameter})
    pointers = Vector{Ptr{UInt8}}(undef, length(parameters))

    map!(pointers, parameters) do parameter
        ismissing(parameter) ? C_NULL : pointer(parameter)
    end

    return pointers
end

"""
    num_params(jl_result::Result) -> Int

Return the number of parameters in a prepared statement.
If this result did not come from the description of a prepared statement, return 0.
"""
function num_params(jl_result::Result)::Int
    # todo: check cleared?
    libpq_c.PQnparams(jl_result.result)
end

"""
    num_rows(jl_result::Result) -> Int

Return the number of rows in the query result.
This will be 0 if the query would never return data.
"""
function num_rows(jl_result::Result)::Int
    # todo: check cleared?
    libpq_c.PQntuples(jl_result.result)
end

"""
    num_affected_rows(jl_result::Result) -> Int

Return the number of rows affected by the command returning the result.
This is useful for counting the rows affected by operations such as INSERT,
UPDATE and DELETE that do not return rows but affect them.
This will be 0 if the query does not affect any row.
"""
function num_affected_rows(jl_result::Result)::Int
    # todo: check cleared?
    str = unsafe_string(libpq_c.PQcmdTuples(jl_result.result))
    if isempty(str)
        throw(ArgumentError("Result generated by an incompatible command"))
    else
        return parse(Int, str)
    end
end

"""
    num_columns(jl_result::Result) -> Int

Return the number of columns in the query result.
This will be 0 if the query would never return data.
"""
function num_columns(jl_result::Result)::Int
    # todo: check cleared?
    libpq_c.PQnfields(jl_result.result)
end

"""
    column_name(jl_result::Result, column_number::Integer) -> String

Return the name of the column at index `column_number` (1-based).
"""
function column_name(jl_result::Result, column_number::Integer)
    return jl_result.column_names[column_number]
end

"""
    column_names(jl_result::Result) -> Vector{String}

Return the names of all the columns in the query result.
"""
column_names(jl_result::Result) = copy(jl_result.column_names)

"""
    column_number(jl_result::Result, column_name::Union{AbstractString, Symbol}) -> Int

Return the index (1-based) of the column named `column_name`.
"""
function column_number(jl_result::Result, column_name::Union{AbstractString, Symbol})::Int
    return something(findfirst(isequal(String(column_name)), jl_result.column_names), 0)
end

"""
    column_number(jl_result::Result, column_idx::Integer) -> Int

Return the index of the column if it is valid, or error.
"""
function column_number(jl_result::Result, column_idx::Integer)::Int
    @boundscheck if !checkindex(Bool, Base.OneTo(num_columns(jl_result)), column_idx)
        throw(BoundsError(column_names(jl_result), column_idx))
    end

    return column_idx
end

"""
    column_oids(jl_result::Result) -> Vector{LibPQ.Oid}

Return the PostgreSQL oids for each column in the result.
"""
column_oids(jl_result::Result) = jl_result.column_oids

"""
    column_types(jl_result::Result) -> Vector{Type}

Return the corresponding Julia types for each column in the result.
"""
column_types(jl_result::Result) = jl_result.column_types

"""
    getindex(jl_result::Result, row::Integer, col::Integer) -> Union{_, Missing}

Return the parsed value of the result at the row and column specified (1-indexed).
The returned value will be `missing` if `NULL`, or will be of the type specified in
[`column_types`](@ref).
"""
function Base.getindex(jl_result::Result, row::Integer, col::Integer)
    if isnull(jl_result, row, col)
        return missing
    else
        oid = column_oids(jl_result)[col]
        T = column_types(jl_result)[col]
        return jl_result.column_funcs[col](PQValue{oid}(jl_result, row, col))::T
    end
end


text:
include("Misc/acb_root_ctx.jl")
include("Misc/fmpr.jl")
include("Misc/FactoredElem.jl")
include("Misc/Integer.jl")
include("Misc/coprime.jl")
include("Misc/Magma.jl")
include("Misc/Matrix.jl")
include("Misc/NumberField.jl")
include("Misc/Poly.jl")
include("Misc/Residue.jl")
include("Misc/ResidueRingPoly.jl")
include("Misc/pAdic.jl")
include("Misc/FiniteField.jl")
include("Misc/Plesken.jl")
include("Misc/Arb.jl")
include("Misc/psibound.jl")
include("Misc/AbelianGrp.jl")
include("Misc/UnitsModM.jl")
include("Misc/ZZModUInt.jl")

ext: jl


text:
module MaxHelpingHandSidewaysJumpThru

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/SidewaysJumpThru" SidewaysJumpThru(x::Integer, y::Integer, height::Integer=Maple.defaultBlockHeight,
    left::Bool=true, texture::String="wood", animationDelay::Number=0.0, letSeekersThrough::Bool=false, surfaceIndex::Int=-1, pushPlayer::Bool=false)
@mapdef Entity "MaxHelpingHand/AttachedSidewaysJumpThru" AttachedSidewaysJumpThru(x::Integer, y::Integer, height::Integer=Maple.defaultBlockHeight,
    left::Bool=true, texture::String="wood", animationDelay::Number=0.0, letSeekersThrough::Bool=false, surfaceIndex::Int=-1, pushPlayer::Bool=false)

textures = ["wood", "dream", "temple", "templeB", "cliffside", "reflection", "core", "moon"]
const placements = Ahorn.PlacementDict(
    "Sideways Jump Through (Left) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        SidewaysJumpThru,
        "rectangle",
        Dict{String, Any}(
            "texture" => "wood",
            "left" => true
        )
    ),
    "Sideways Jump Through (Right) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        SidewaysJumpThru,
        "rectangle",
        Dict{String, Any}(
            "texture" => "wood",
            "left" => false
        )
    ),
    "Sideways Jump Through (Attached, Left)\n(max480's Helping Hand)" => Ahorn.EntityPlacement(
        AttachedSidewaysJumpThru,
        "rectangle",
        Dict{String, Any}(
            "texture" => "wood",
            "left" => true
        )
    ),
    "Sideways Jump Through (Attached, Right)\n(max480's Helping Hand)" => Ahorn.EntityPlacement(
        AttachedSidewaysJumpThru,
        "rectangle",
        Dict{String, Any}(
            "texture" => "wood",
            "left" => false
        )
    )
)

# Sideways Jumpthrus and Attached ones look the same and have the same options, so we're going to handle them together!
const jumpthruUnion = Union{SidewaysJumpThru, AttachedSidewaysJumpThru}

quads = Tuple{Integer, Integer, Integer, Integer}[
    (0, 0, 8, 7) (8, 0, 8, 7) (16, 0, 8, 7);
    (0, 8, 8, 5) (8, 8, 8, 5) (16, 8, 8, 5)
]

Ahorn.editingOptions(entity::jumpthruUnion) = Dict{String, Any}(
    "texture" => textures,
    "surfaceIndex" => Maple.tileset_sound_ids
)

Ahorn.minimumSize(entity::jumpthruUnion) = 0, 8
Ahorn.resizable(entity::jumpthruUnion) = false, true

function Ahorn.selection(entity::jumpthruUnion)
    x, y = Ahorn.position(entity)
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, 8, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::jumpthruUnion, room::Maple.Room)
    texture = get(entity.data, "texture", "wood")
    texture = texture == "default" ? "wood" : texture

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    height = Int(get(entity.data, "height", 8))
    left = get(entity.data, "left", true)

    startX = div(x, 8) + 1
    startY = div(y, 8) + 1
    stopY = startY + div(height, 8) - 1
    animated = Number(get(entity.data, "animationDelay", 0)) > 0

    Ahorn.Cairo.save(ctx)

    Ahorn.rotate(ctx, pi / 2)

    if left
        Ahorn.scale(ctx, 1, -1)
    end

    len = stopY - startY
    for i in 0:len
        if animated
            Ahorn.drawImage(ctx, "objects/jumpthru/$(texture)00", 8 * i, left ? 0 : -8)
        else
            connected = false
            qx = 2
            if i == 0
                connected = get(room.fgTiles.data, (startY - 1, startX), false) != '0'
                qx = 1

            elseif i == len
                connected = get(room.fgTiles.data, (stopY + 1, startX), false) != '0'
                qx = 3
            end

            quad = quads[2 - connected, qx]
            Ahorn.drawImage(ctx, "objects/jumpthru/$(texture)", 8 * i, left ? 0 : -8, quad...)
        end
    end

    Ahorn.Cairo.restore(ctx)
end

end

max_stars_repo_name: swoolcock/MaxHelpingHand

max_issues_repo_name: swoolcock/MaxHelpingHand

max_forks_repo_name: swoolcock/MaxHelpingHand


text:
## FUNTIONS TO LOAD MODEL IMPLEMENTATION CODE

# helper:

"""
    load(name::String; pkg=nothing, modl=Main, verbosity=1, name=nothing)

Load the model implementation code for the model type with specified
`name` (first argument) into the module `modl`, specifying `pkg` if
necesssary to resolve duplicate names. Return an instance of the model type
with default hyper-parameters.

To bind the name of the model type to a name different from the first
argument, specify `name=...`.

    load(proxy; kwargs...1)

In the case that `proxy` is a return value of `traits` (ie, has the
form `(name = ..., package_name = ..., etc)`) this is equivalent to
`load(name=proxy.name, pkg=proxy.package_name; kwargs...)`.

See also [`@load`](@ref)

### Examples

    load("RidgeRegressor",
         pkg="MLJLinearModels",
         modl=MyPackage,
         name="NewRidge",
         verbosity=0)

    load(localmodels()[1])

See also [`@load`](@ref)

"""
function load(proxy::ModelProxy; modl=Main, verbosity=1, name=nothing)

    # decide what to print
    toprint = verbosity > 0

    # get name, package and load path:
    model_name = proxy.name
    pkg = proxy.package_name
    handle = (name=model_name, pkg=pkg)

    # see if the model has already been loaded:
    type_already_loaded = handle in map(localmodels(modl=modl)) do m
        (name=m.name, pkg=m.package_name)
    end

    if type_already_loaded && name == nothing
        toprint && @info "Model code for $model_name already loaded"
        # return an instance of the type:
        for M in localmodeltypes(modl)
             i = info(M)
            if i.name == model_name && i.package_name == pkg
                return M()
            end
        end
    end

    # from now on work with symbols not strings:
    path = INFO_GIVEN_HANDLE[handle][:load_path]
    path_components = Symbol.(split(path, '.') )
    model_name = Symbol(model_name)
    pkg = Symbol(pkg)

    verbosity > 0 && @info "Loading into module \"$modl\": "

    # create `name`, the actual name (symbol) to be bound to the
    # new type in the global namespace of `modl`. Note two packages
    # might provide models with the same `model_name`.
    if name === nothing
        name = MLJBase.available_name(modl, Symbol(model_name))
        name == model_name || verbosity < 0 ||
            @warn "New model type being bound to "*
            "`$name` to avoid conflict with an existing name. "
    else
        name = Symbol(name)
    end

    # the package providing the implementation of the MLJ model
    # interface could be different from `pkg`, which is the package
    # name exposed to the user (and where the core alogrithms live,
    # but that this not actually relevant here).
    api_pkg = path_components[1]

    # if needed, put MLJModels in the calling module's namespace:
    if api_pkg == :MLJModels
        toprint && print("import MLJModels ")
        # the latter pass exists as a fallback, and shouldn't happen
        # for end users, but we need that for testing, etc
        load_ex =
            isdefined(modl, :MLJ) ? :(import MLJ.MLJModels) :
            :(import MLJModels)
        modl.eval(load_ex)
        toprint && println('\u2714')
    end

    # load `api_pkg`, unless this is MLJModels, in which case load
    # `pkg` to trigger lazy-loading of implementation code:
    pkg_ex = api_pkg == :MLJModels ? pkg : api_pkg
    toprint && print("import $pkg_ex ")
    modl.eval(:(import $pkg_ex))

    toprint && println('\u2714')

    # load the model:
    # the latter pass exists as a fallback, and shouldn't happen for
    # end users, but we need that for testing, etc
    if api_pkg == "MLJModels" && isdefined(modl, :MLJ)
        pushfirst!(path_components, :MLJ)
    end
    root_components = path_components[1:(end - 1)]
    root_str = join(string.(root_components), '.')
    import_ex = Expr(:import, Expr(:(.), root_components...))
    path_str = join(string.(path_components), '.')
    path_ex = path_str |> Meta.parse
    program =
        quote
            $import_ex
            const $name = $path_ex
        end

    toprint && print("import ", root_str, " ")
    modl.eval(program)
    toprint && println('\u2714')

    return modl.eval(:($name()))

end

load(name::String; pkg=nothing, kwargs...) =
    load(info(name; pkg=pkg); kwargs...)

"""
    @load name pkg=nothing verbosity=0 name=nothing

Load the model implementation code for the model named in the first
argument into the calling module, specfying `pkg` in the case of
duplicate names. Return a model instance with default
hyper-parameters.

To bind the new model type to a name different from the first
argument, specify `name=...`.

### Examples

    @load DecisionTreeeRegressor
    @load PCA verbosity=1
    @load SVC pkg=LIBSVM name=MyOtherSVC

See also [`load`](@ref)

"""
macro load(name_ex, kw_exs...)
    name_ = string(name_ex)

    # parse kwargs:
    warning = "Invalid @load syntax.\n "*
    "Sample usage: @load PCA pkg=\"MultivariateStats\" verbosity=1"
    for ex in kw_exs
        ex.head == :(=) || throw(ArgumentError(warning))
        variable_ex = ex.args[1]
        value_ex = ex.args[2]
        if variable_ex == :pkg
            pkg = string(value_ex)
        elseif variable_ex == :verbosity
            verbosity = value_ex
        elseif variable_ex == :name
	    model_name = value_ex
	else
            throw(ArgumentError(warning))
        end
    end
    (@isdefined pkg) || (pkg = nothing)
    (@isdefined verbosity) || (verbosity = 0)
    (@isdefined model_name) || (model_name = nothing)
    # get rid brackets in name_, as in
    # "(MLJModels.Clustering).KMedoids":
    name = filter(name_) do c !(c in ['(',')']) end

    load(name, modl=__module__, pkg=pkg, verbosity=verbosity, name=model_name)

end


text:
# Many derivative free methods of different orders
#
# TODO: rework Order5 #https://pdfs.semanticscholar.org/ce50/3210d96f653a14b28da96600d5990d2abe97.pdf
# https://content.sciendo.com/view/journals/tmj/10/4/article-p103.xml 7 and 8
# order8: https://www.hindawi.com/journals/ijmms/2012/493456/ref/

##################################################
## Guard against non-robust algorithms
##
## By default, we do this by deciding if we
## should take a secant step instead of the algorithm For example, for
## Steffensen which is quadratically convergent and the Secant method
## which is only superlinear,
## the error, e_{n+1} = x_{n+1} - alpha, may be smaller after a secant
## step than a Steffensen step. (It is only once x_n is close enough
## to alpha that the method is quadratically convergent.
## The Steffensen error is
## Δn+1 = f[x,x+fx, alpha]/f[x, x+fx] * (1 - f[x, alpha]) (x-alpha)^2
##      ≈ f''/(2f') * ( 1 + f') Δn^2
## The Secant error is
## Δn+1 = f[x,x_{-1},alpha] / f[x,x_{-1}] * (x-alpha) * (x_{-1} - alpha)
##      ≈  f''/(2f')  Δn ⋅ Δn-1
## The ratio is ≈ (1 + f')(Δn / Δn-1)
## It seems reasonable, that a Steffensen step is preferred when
## the ratio satisfies -1 < (1+f') ⋅ Δn /Δn-1 < 1
## We could use f' ~ fp = (fx1-fx0)/(x1-x0); but our proxy for
## Δn/Δn-1 is problematic, as we don't know alpha, and using xn-x_{n-1}
## can be an issue when only x1 and not x0 is specified. This needs
## working around.
##
## Instead, as Steffensen is related to Newton as much as
## (f(x+fx) - fx)/fx  ≈ f'(x), we take a Steffensen step if |fx|
## is small enough. For this we use |fx| <= x/1000; which
## seems to work reasonably well over several different test cases.
@inline function do_guarded_step(M::AbstractSecant, o::UnivariateZeroState{T,S}) where {T, S}
    x, fx = o.xn1, o.fxn1
    1000 * abs(fx) >  max(oneunit(S), abs(x) * oneunit(S) /oneunit(T)) * one(T)
end

# check if we should guard against step for method M; call N if yes, P if not
function update_state_guarded(M::AbstractSecant,N::AbstractUnivariateZeroMethod, P::AbstractUnivariateZeroMethod, fs, o, options)
    if do_guarded_step(M, o)
        #@debug "do secant step"
        return update_state(N, fs, o, options)
    else
        #@debug "do $N step"
        update_state(P, fs, o, options)
    end
end

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

function init_state(method::AbstractSecant, fs, x::Number)
    x1 = float(x)
    x0 = _default_secant_step(x1)
    init_state(method, fs, (x0, x1))
end

function init_state(method::AbstractSecant, fs, x::Union{Tuple, Vector})
    x0, x1 = promote(float(x[1]), float(x[2]))
    fx0, fx1 = fs(x0), fs(x1)
    state = UnivariateZeroState(x1, x0, zero(x1)/zero(x1)*oneunit(x1), eltype(x1)[],
                                fx1, fx0, fx1, eltype(fx1)[],
                                0, 2,
                                false, false, false, false, "")

    state
end

function init_state!(state::UnivariateZeroState{T, S}, method::AbstractSecant, fs, x::Number) where {T, S}
    x1::T = float(x)
    x0::T = _default_secant_step(x1)
    init_state!(state, method, fs, (x0, x1))
end

function init_state!(state::UnivariateZeroState{T, S}, ::AbstractSecant, f, x::Union{Tuple, Vector}) where {T, S}
    x0,x1 = promote(float.(x)...)
    fx0, fx1 = promote(f(x0), f(x1))
    init_state!(state, x1, x0, T[], fx1, fx0, S[])
    state.fnevals = 2
    nothing
end

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

## Order0 and Secant are related
"""
    Order0()

The `Order0` method is engineered to be a more robust, though possibly
slower, alternative to the other derivative-free root-finding
methods. The implementation roughly follows the algorithm described in
*Personal Calculator Has Key to Solve Any Equation f(x) = 0*, the
SOLVE button from the
[HP-34C](http://www.hpl.hp.com/hpjournal/pdfs/IssuePDFs/1979-12.pdf).
The basic idea is to use a secant step. If along the way a bracket is
found, switch to bisection, using `AlefeldPotraShi`.  If the secant
step fails to decrease the function value, a quadratic step is used up
to 4 times.

This is not really 0-order: the secant method has order
1.6...[Wikipedia](https://en.wikipedia.org/wiki/Secant_method#Comparison_with_other_root-finding_methods_
and the the bracketing method has order
1.6180...[Wikipedia](http://www.ams.org/journals/mcom/1993-61-204/S0025-5718-1993-1192965-2/S0025-5718-1993-1192965-2.pdf)
so for reasonable starting points, this algorithm should be
superlinear, and relatively robust to non-reasonable starting points.

"""
struct Order0 <: AbstractSecant end

function find_zero(fs, x0, method::Order0;
                   tracks::AbstractTracks=NullTracks(),
                   verbose=false,
                   kwargs...)
    M = Order1()
    N = AlefeldPotraShi()

    find_zero(callable_function(fs), x0, M, N; tracks=tracks,verbose=verbose, kwargs...)
end

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

## Secant
## https://en.wikipedia.org/wiki/Secant_method
"""
    Order1()
    Roots.Secant()

The `Order1()` method is an alias for `Secant`. It specifies the
[secant method](https://en.wikipedia.org/wiki/Secant_method).
This method keeps two values in its state, `x_n` and `x_n1`. The
updated point is the intersection point of x axis with the secant line
formed from the two points. The secant method uses 1 function
evaluation per step and has order `(1+sqrt(5))/2`.

The error, `eᵢ = xᵢ - α`, satisfies
`eᵢ₊₂ = f[xᵢ₊₁,xᵢ,α] / f[xᵢ₊₁,xᵢ] * (xᵢ₊₁-α) * (xᵢ - α)`.

"""
struct Secant <: AbstractSecant end
const Order1 = Secant

function update_state(method::Order1, fs, o::UnivariateZeroState{T,S}, options) where {T, S}

    xn0, xn1 = o.xn0, o.xn1
    fxn0, fxn1 = o.fxn0, o.fxn1

    delta = fxn1 * (xn1 - xn0) / (fxn1 - fxn0)

    if isinf(delta) || isnan(delta)
        o.stopped = true
        o.message = "Increment `Δx` has issues. "
         return
     end

    o.xn0 = xn1
    o.xn1 -= delta
    o.fxn0 = fxn1
    o.fxn1 = (tmp::S = fs(o.xn1))
    incfn(o)

    nothing

end

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

"""
    Roots.Order1B()
    Roots.King()

A superlinear (order 1.6...) modification of the secant method for multiple roots.
Presented in A SECANT METHOD FOR MULTIPLE ROOTS, by RICHARD F. KING, BIT 17 (1977), 321-328

The basic idea is similar to Schroder's method: apply the secant method
to  `f/f'`. However, this uses `f' ~ fp = (fx - f(x-fx))/fx` (a Steffensen step). In
this implementation, `Order1B`, when `fx` is too big, a single secant step of `f`
is used.

The *asymptotic* error, `eᵢ = xᵢ - α`, is given by
`eᵢ₊₂ = 1/2⋅G''/G'⋅ eᵢ⋅eᵢ₊₁ + (1/6⋅G'''/G' - (1/2⋅G''/G'))^2⋅eᵢ⋅eᵢ₊₁⋅(eᵢ+eᵢ₊₁)`.

"""
struct King <: AbstractSecant end
struct Order1B <: AbstractSecant end

function update_state(method::Order1B, fs, o::UnivariateZeroState{T,S}, options)  where {T, S}
    update_state_guarded(method, Secant(), King(), fs, o, options)
end

function update_state(method::King, fs,
                      o::UnivariateZeroState{T,S}, options)  where {T, S}

    x0, x1 = o.xn0, o.xn1
    fx0, fx1 = o.fxn0, o.fxn1

    # G(x,f0,f_1) = -fx^2/(f_1 - f0)
    f0 = fx1
    f_1::S = fs(x1 - f0*oneunit(T)/oneunit(S))
    incfn(o, 1)

    G1 = -f0^2 / (f_1 - f0)

    if isempty(o.fm)
        f0 = fx0
        tmp::S = fs(x0 - f0*oneunit(T)/oneunit(S))
        f_1 = tmp
        incfn(o,1)
        G0 = -f0^2 / (f_1 - f0)
    else
        G0 = first(o.fm)
    end

    m = (x1-x0)/(G1-G0) # approximate value of `m`, the multiplicity
    if abs(m) <= 1e-2 * oneunit(m)
        #@info "small m estimate, stopping"
        o.stopped  = true
        o.message = "Estimate for multiplicity has issues. "
        return
    end

    delta = G1 * (x1 - x0) / (G1 - G0)
    empty!(o.fm); push!(o.fm, G1)

    if isissue(delta)
        o.stopped  = true
        o.message = "Increment `Δx` has issues. "
        return
    end

    o.xn0, o.fxn0 = o.xn1, o.fxn1
    o.xn1 -= delta
    o.fxn1 = fs(o.xn1)
    incfn(o)

    nothing
end

### Steffensen
## https://en.wikipedia.org/wiki/Steffensen's_method#Simple_description

"""
    Order2()
    Roots.Steffensen()

The quadratically converging
[Steffensen](https://en.wikipedia.org/wiki/Steffensen's_method#Simple_description)
method is used for the derivative-free `Order2()` algorithm. Unlike
the quadratically converging Newton's method, no derivative is
necessary, though like Newton's method, two function calls per step
are. Steffensen's algorithm is more sensitive than Newton's method to
poor initial guesses when `f(x)` is large, due to how `f'(x)` is
approximated. This algorithm, `Order2`, replaces a Steffensen step with a secant
step when `f(x)` is large.

The error, `eᵢ - α`, satisfies
`eᵢ₊₁ = f[xᵢ, xᵢ+fᵢ, α] / f[xᵢ,xᵢ+fᵢ] ⋅ (1 - f[xᵢ,α] ⋅ eᵢ²`
"""
struct Order2 <: AbstractSecant end
struct Steffensen <: AbstractSecant end

function update_state(method::Order2, fs, o::UnivariateZeroState{T,S}, options)  where {T, S}
     update_state_guarded(method, Secant(), Steffensen(), fs, o, options)
end

function update_state(method::Steffensen, fs, o::UnivariateZeroState{T,S}, options)  where {T, S}

    x0, x1 = o.xn0, o.xn1
    fx0, fx1 = o.fxn0, o.fxn1

    sgn = sign((fx1 - fx0) / (x1 - x0))
    x2 = x1 - sgn * fx1 / oneunit(S) * oneunit(T)

    f0 = fx1
    f1::S = fs(x2)
    incfn(o, 1)

    delta = -sgn * f0 * f0 / (f1 - f0) * oneunit(T) / oneunit(S)

    if isissue(delta)
        o.stopped  = true
        o.message = "Increment `Δx` has issues. "
        return
    end

    o.xn0, o.fxn0 = o.xn1, o.fxn1
    o.xn1 -= delta
    o.fxn1 = fs(o.xn1)
    incfn(o)

    nothing
end

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

### Order2B() Esser method
"""
    Roots.Order2B()
    Roots.Esser()

Esser's method. This is a quadratically convergent method that, like
Schroder's method, does not depend on the multiplicity of the
zero. Schroder's method has update step `x - r2/(r2-r1) * r1`, where `ri =
f^(i-1)/f^(i)`. Esser approximates `f' ~ f[x-h, x+h], f'' ~
f[x-h,x,x+h]`, where `h = fx`, as with Steffensen's method, Requiring 3
function calls per step. The implementation `Order2B` uses a secant
step when `|fx|` is considered too large.

Esser, H. Computing (1975) 14: 367. https://doi.org/10.1007/BF02253547
Eine stets quadratisch konvergente Modifikation des Steffensen-Verfahrens

Example
```
f(x) = cos(x) - x
g(x) = f(x)^2
x0 = pi/4
find_zero(f, x0, Order2(), verbose=true)        #  3 steps / 7 function calls
find_zero(f, x0, Roots.Order2B(), verbose=true) #  4 / 9
find_zero(g, x0, Order2(), verbose=true)        #  22 / 45
find_zero(g, x0, Roots.Order2B(), verbose=true) #  4 / 10
```
"""
struct Esser <: AbstractSecant end
struct Order2B <: AbstractSecant end

function update_state(method::Order2B, fs, o::UnivariateZeroState{T,S}, options)  where {T, S}
     update_state_guarded(method, Secant(), Esser(), fs, o, options)
end

function update_state(method::Esser, fs,
                      o::UnivariateZeroState{T,S}, options)  where {T, S}

    x1, fx1 = o.xn1, o.fxn1

    f0 = fx1

    f1  = fs(x1 + f0 * oneunit(T) / oneunit(S))
    f_1 = fs(x1 - f0 * oneunit(T) / oneunit(S))
    incfn(o, 2)

    # h = f0
    # r1 = f/f' ~ f/f[x+h,x-h]
    # r2 = f'/f'' ~ f[x+h, x-h]/f[x-h,x,x+h]
    r1 = f0 * 2*f0 / (f1 - f_1) * oneunit(T) / oneunit(S)
    r2 = (f1 - f_1)/(f1 - 2*f0 + f_1) * f0/2 * oneunit(T) / oneunit(S)

    k = r2/(r2-r1)  # ~ m

    if abs(k) <= 1e-2 * oneunit(k)
        #@info "estimate for m is *too* small"
        o.stopped  = true
        o.message = "Estimate for multiplicity had issues. "
        return
    end

    delta = k * r1

    if isissue(delta)
        o.stopped  = true
        o.message = "Increment `Δx` has issues. "
        return
    end

    o.xn0, o.fxn0 = o.xn1, o.fxn1
    o.xn1 -= delta

    o.fxn1 = fs(o.xn1)
    incfn(o)

    nothing
end

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

"""
    Order5()

Implements an order 5 algorithm from *A New Fifth Order Derivative
Free Newton-Type Method for Solving Nonlinear Equations* by Manoj
Kumar, Akhilesh Kumar Singh, and Akanksha, Appl. Math. Inf. Sci. 9,
No. 3, 1507-1513 (2015), DOI: [10.12785/amis/090346](https://doi.org/10.12785/amis/090346). Four function
calls per step are needed.

The error, `eᵢ = xᵢ - α`, satisfies
`eᵢ₊₁ = K₁ ⋅ K₅ ⋅ M ⋅ eᵢ⁵ + O(eᵢ⁶)`

"""
struct Order5 <: AbstractSecant end
struct KumarSinghAkanksha <: AbstractSecant end

function update_state(M::Union{Order5, KumarSinghAkanksha}, fs, o::UnivariateZeroState{T,S},
                      options) where {T, S}

    xn = o.xn1
    fxn = o.fxn1

    wn::T = steff_step(M, o.xn1, o.fxn1)

    fwn::S = fs(wn)
    incfn(o)

    fp, issue = _fbracket(o.xn1, wn, o.fxn1, fwn)
    if issue
        o.xn0, o.xn1 = o.xn1, wn
        o.fxn0, o.fxn1 = o.fxn1, fwn
        o.message = "Issue with divided difference f[xn, wn]. "
        o.stopped  = true
        return
    end

    yn::T = o.xn1 - o.fxn1 / fp
    fyn::S = fs(yn)
    incfn(o)

    zn::T = xn - (fxn + fyn) / fp
    fzn::S = fs(zn)
    incfn(o)

    fp, issue = _fbracket_ratio(yn, o.xn1, wn, fyn, o.fxn1, fwn)
    if issue
        o.xn0, o.xn1 = o.xn1, yn
        o.fxn0, o.fxn1 = o.fxn1, fyn
        o.message = "Issue with f[xn,yn]*f[yn,wn] / f[xn, wn]"
        o.stopped = true
        return
    end

    o.xn0 = o.xn1
    o.fxn0 = o.fxn1
    o.xn1 = zn  - fzn  / fp
    o.fxn1 = fs(o.xn1)
    incfn(o)

    nothing
end

## If we have a derivative, we have this. (Deprecate?)
function update_state(method::Order5, fs::Union{FirstDerivative,SecondDerivative},
                      o::UnivariateZeroState{T,S}, options)  where {T, S}

    xn, fxn = o.xn1, o.fxn1
    a::T, b::S = fΔx(fs, xn)
    fpxn = a/b
    incfn(o)

    if isissue(fpxn)
        o.stopped  = true
        return
    end

    yn::T = xn - fxn / fpxn
    fyn::S, Δyn::T = fΔx(fs, yn)
    fpyn = fyn / Δyn
    incfn(o, 2)

    if isissue(fpyn)
        o.xn0, o.xn1 = xn, yn
        o.fxn0, o.fxn1 = fxn, fyn
        o.stopped  = true
        return
    end

    zn::T = xn  - (fxn + fyn) / fpxn
    fzn::S = fs(zn)
    incfn(o, 1)

    xn1 = zn - fzn / fpyn
    fxn1 = fs(xn1)
    incfn(o, 1)

    o.xn0, o.xn1 = xn, xn1
    o.fxn0, o.fxn1 = fxn, fxn1

    nothing
end

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

## cf also: https://doi.org/10.1515/tmj-2017-0049
"""
    Order8()

Implements an eighth-order algorithm from *New Eighth-Order
Derivative-Free Methods for Solving Nonlinear Equations* by Rajinder
Thukral, International Journal of Mathematics and Mathematical
Sciences Volume 2012 (2012), Article ID 493456, 12 pages DOI:
[10.1155/2012/493456](https://doi.org/10.1155/2012/493456). Four function calls per step are required.

The error, `eᵢ = xᵢ - α`, is expressed as `eᵢ₊₁ = K ⋅ eᵢ⁸` in
(2.25) of the paper for an explicit `K`.

"""
struct Order8 <: AbstractSecant end
struct Thukral8 <: AbstractSecant end

function update_state(M::Union{Thukral8, Order8}, fs, o::UnivariateZeroState{T,S}, options) where {T, S}

    xn = o.xn1
    fxn = o.fxn1

    wn::T = steff_step(M, xn, fxn)
    fwn::S = fs(wn)
    incfn(o)

    if isissue(fwn)
        o.xn0,o.xn1 = xn, wn
        o.fxn0,o.fxn1 = fxn, fwn
        o.stopped = true
        o.message = "issue with Steffensen step fwn"
        return
    end

    fp, issue = _fbracket(xn, wn, fxn, fwn)

    if issue
        o.stopped = true
        o.message = "issue with divided difference f[xn, wn]. "
        return
    end

    yn::T = xn - fxn / fp
    fyn::S = fs(yn)
    incfn(o)

    fp, issue = _fbracket(yn, xn, fyn, fxn)
    if issue #fp
        o.xn0,o.xn1 = xn, yn
        o.fxn0,o.fxn1 = fxn, fyn
        o.stopped = true
        o.message = "issue with divided difference f[xn, yn]. "
        return
    end

    phi = (1 + fyn / fwn)           # pick one of options
    zn =  yn - phi * fyn / fp
    fzn::S = fs(zn)
    incfn(o)

    fp, issue =  _fbracket_diff(xn, yn, zn, fxn, fyn, fzn)
    if issue
        o.xn0,o.xn1 = xn, zn
        o.fxn0,o.fxn1 = fxn, fzn
        o.message = "issue with divided difference  f[y,z] - f[x,y] + f[x,z]. "
        o.stopped = true
        return
    end

    w = 1 / (1 - fzn/fwn)

    xi = (1 - 2fyn*fyn*fyn / (fwn * fwn * fxn))

    xn1 = zn - w * xi * fzn / fp
    fxn1::S = fs(xn1)
    incfn(o)

    o.xn0,o.xn1 = xn, xn1
    o.fxn0,o.fxn1 = fxn, fxn1

    nothing
end

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

"""
    Order16()

Implements the order 16 algorithm from
*New Sixteenth-Order Derivative-Free Methods for Solving Nonlinear Equations*
by R. Thukral,
American Journal of Computational and Applied Mathematics
p-ISSN: 2165-8935;    e-ISSN: 2165-8943; 2012;  2(3): 112-118
doi: 10.5923/j.ajcam.20120203.08.

Five function calls per step are required. Though rapidly converging,
this method generally isn't faster (fewer function calls/steps) over
other methods when using `Float64` values, but may be useful for
solving over `BigFloat`.

The error, `eᵢ = xᵢ - α`, is expressed as `eᵢ₊₁ = K⋅eᵢ¹⁶` for an explicit `K`
in equation (50) of the paper.

"""
struct Order16 <: AbstractSecant end
struct Thukral16 <: AbstractSecant end

function update_state(M::Union{Thukral16, Order16}, fs, o::UnivariateZeroState{T,S}, options) where {T, S}
    xn = o.xn1
    fxn = o.fxn1

    wn::T = steff_step(M, xn, fxn)
    fwn::S = fs(wn)
    incfn(o)

    fp, issue = _fbracket(xn, wn, fxn, fwn)

    if issue
        o.xn0, o.xn1 = xn, wn
        o.fxn0, o.fxn1 = fxn, fwn
        o.message = "issue with f[xn,wn]"
        o.stopped = true
        return
    end

    yn::T = xn - fxn / fp
    fyn::S = fs(yn)
    incfn(o)

    fp, issue = _fbracket_ratio(yn, xn, wn, fyn, fxn, fwn)
    if issue
        o.xn0, o.xn1 = xn, yn
        o.fxn0, o.fxn1 = fxn, fyn
        o.message = "issue with f[xn,yn]*f[yn,wn]/f[xn,wn]"
        o.stopped = true
        return
    end

    zn = yn - fyn / fp
    fzn::S = fs(zn)
    incfn(o)

    fp, issue = _fbracket_diff(xn, yn, zn, fxn, fyn, fzn)
    u2, u3, u4 = fzn/fwn, fyn/fxn, fyn/fwn

    eta = 1 / (1 + 2*u3*u4^2) / (1 - u2)
    if issue
        o.xn0, o.xn1 = xn, zn
        o.fxn0, o.fxn1 = fxn, fzn
        o.stopped = true
        o.message = "Approximate derivative failed"
        return
    end

    an = zn - eta * fzn / fp
    fan::S = fs(an)
    incfn(o)

    fp, issue = _fbracket_ratio(an, yn, zn, fan, fyn, fzn)
    if issue
        o.xn0, o.xn1 = xn, an
        o.fxn0, o.fxn1 = fxn, fan
        o.stopped = true
        o.message = "Approximate derivative failed"
        return
    end

    u1, u5, u6 = fzn/fxn, fan/fxn, fan/fwn

    fp1, issue = _fbracket(xn,yn, fxn, fyn)

    sigma =  1 + u1*u2 - u1*u3*u4^2 + u5 + u6 + u1^2*u4 +
    u2^2*u3 + 3*u1*u4^2*(u3^2 - u4^2)/(fp1/oneunit(fp1))

    xn1 = an - sigma * fan / fp
    fxn1::S = fs(xn1)
    incfn(o)

    o.xn0, o.xn1 = xn, xn1
    o.fxn0, o.fxn1 = fxn, fxn1

    nothing
end

##################################################
## some means of guarding against large fx when taking a secant step
## TODO: rework this
function steff_step(M::Union{Order5, Order8, Order16}, x::S, fx::T) where {S, T}

    xbar, fxbar = real(x/oneunit(x)), fx/oneunit(fx)
    thresh =  max(1, abs(xbar)) * sqrt(eps(one(xbar))) #^(1/2) # max(1, sqrt(abs(x/fx))) * 1e-6

    out = abs(fxbar) <= thresh ? fxbar  : sign(fx) * thresh
    x + out * oneunit(x)

end


text:
using SphereQuadratures
using Test

@testset "SphereQuadratures.jl" begin
    # Write your tests here.
end


text:
"""
    complement(g)

Return the [graph complement](https://en.wikipedia.org/wiki/Complement_graph)
of a graph

### Implementation Notes
Preserves the `eltype` of the input graph.

# Examples
```jldoctest
julia> g = SimpleDiGraph([0 1 0 0 0; 0 0 1 0 0; 1 0 0 1 0; 0 0 0 0 1; 0 0 0 1 0]);

julia> foreach(println, edges(complement(g)))
Edge 1 => 3
Edge 1 => 4
Edge 1 => 5
Edge 2 => 1
Edge 2 => 4
Edge 2 => 5
Edge 3 => 2
Edge 3 => 5
Edge 4 => 1
Edge 4 => 2
Edge 4 => 3
Edge 5 => 1
Edge 5 => 2
Edge 5 => 3
```
"""
function complement(g::Graph)
    gnv = nv(g)
    h = SimpleGraph(gnv)
    for i = 1:gnv
        for j = (i + 1):gnv
            if !has_edge(g, i, j)
                add_edge!(h, i, j)
            end
        end
    end
    return h
end

function complement(g::DiGraph)
    gnv = nv(g)
    h = SimpleDiGraph(gnv)
    for i in vertices(g), j in vertices(g)
        if i != j && !has_edge(g, i, j)
            add_edge!(h, i, j)
        end
    end
    return h
end

"""
    reverse(g)

Return a directed graph where all edges are reversed from the
original directed graph.

### Implementation Notes
Preserves the eltype of the input graph.

# Examples
```jldoctest
julia> g = SimpleDiGraph([0 1 0 0 0; 0 0 1 0 0; 1 0 0 1 0; 0 0 0 0 1; 0 0 0 1 0]);

julia> foreach(println, edges(reverse(g)))
Edge 1 => 3
Edge 2 => 1
Edge 3 => 2
Edge 4 => 3
Edge 4 => 5
Edge 5 => 4
```
"""
function reverse end
@traitfn function reverse(g::G::IsDirected) where G<:AbstractSimpleGraph
    gnv = nv(g)
    gne = ne(g)
    h = SimpleDiGraph(gnv)
    h.fadjlist = deepcopy_adjlist(g.badjlist)
    h.badjlist = deepcopy_adjlist(g.fadjlist)
    h.ne = gne
    return h
end

"""
    reverse!(g)

In-place reverse of a directed graph (modifies the original graph).
See [`reverse`](@ref) for a non-modifying version.
"""
function reverse! end
@traitfn function reverse!(g::G::IsDirected) where G<:AbstractSimpleGraph
    g.fadjlist, g.badjlist = g.badjlist, g.fadjlist
    return g
end

"""
    blockdiag(g, h)

Return a graph with ``|V(g)| + |V(h)|`` vertices and ``|E(g)| + |E(h)|``
edges where the vertices and edges from graph `h` are appended to graph `g`.

### Implementation Notes
Preserves the eltype of the input graph. Will error if the
number of vertices in the generated graph exceeds the `eltype`.

# Examples
```jldoctest
julia> g1 = SimpleDiGraph([0 1 0 0 0; 0 0 1 0 0; 1 0 0 1 0; 0 0 0 0 1; 0 0 0 1 0]);

julia> g2 = SimpleDiGraph([0 1 0; 0 0 1; 1 0 0]);

julia> blockdiag(g1, g2)
{8, 9} directed simple Int64 graph

julia> foreach(println, edges(blockdiag(g1, g2)))
Edge 1 => 2
Edge 2 => 3
Edge 3 => 1
Edge 3 => 4
Edge 4 => 5
Edge 5 => 4
Edge 6 => 7
Edge 7 => 8
Edge 8 => 6
```
"""
function blockdiag(g::T, h::T) where T <: AbstractGraph
    gnv = nv(g)
    r = T(gnv + nv(h))
    for e in edges(g)
        add_edge!(r, e)
    end
    for e in edges(h)
        add_edge!(r, gnv + src(e), gnv + dst(e))
    end
    return r
end

"""
    intersect(g, h)

Return a graph with edges that are only in both graph `g` and graph `h`.

### Implementation Notes
This function may produce a graph with 0-degree vertices.
Preserves the eltype of the input graph.

# Examples
```jldoctest
julia> g1 = SimpleDiGraph([0 1 0 0 0; 0 0 1 0 0; 1 0 0 1 0; 0 0 0 0 1; 0 0 0 1 0]);

julia> g2 = SimpleDiGraph([0 1 0; 0 0 1; 1 0 0]);

julia> foreach(println, edges(intersect(g1, g2)))
Edge 1 => 2
Edge 2 => 3
Edge 3 => 1
```
"""
function intersect(g::T, h::T) where T <: AbstractGraph
    gnv = nv(g)
    hnv = nv(h)

    r = T(min(gnv, hnv))
    for e in intersect(edges(g), edges(h))
        add_edge!(r, e)
    end
    return r
end

"""
    difference(g, h)

Return a graph with edges in graph `g` that are not in graph `h`.

### Implementation Notes
Note that this function may produce a graph with 0-degree vertices.
Preserves the `eltype` of the input graph.

# Examples
```jldoctest
julia> g1 = SimpleDiGraph([0 1 0 0 0; 0 0 1 0 0; 1 0 0 1 0; 0 0 0 0 1; 0 0 0 1 0]);

julia> g2 = SimpleDiGraph([0 1 0; 0 0 1; 1 0 0]);

julia> foreach(println, edges(difference(g1, g2)))
Edge 3 => 4
Edge 4 => 5
Edge 5 => 4
```
"""
function difference(g::T, h::T) where T <: AbstractGraph
    gnv = nv(g)
    hnv = nv(h)

    r = T(gnv)
    for e in edges(g)
        !has_edge(h, e) && add_edge!(r, e)
    end
    return r
end

"""
    symmetric_difference(g, h)

Return a graph with edges from graph `g` that do not exist in graph `h`,
and vice versa.

### Implementation Notes
Note that this function may produce a graph with 0-degree vertices.
Preserves the eltype of the input graph. Will error if the
number of vertices in the generated graph exceeds the eltype.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleGraph(3); h = SimpleGraph(3);

julia> add_edge!(g, 1, 2);

julia> add_edge!(h, 1, 3);

julia> add_edge!(h, 2, 3);

julia> f = symmetric_difference(g, h);

julia> collect(edges(f))
3-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 1 => 3
 Edge 2 => 3
```
"""
function symmetric_difference(g::T, h::T) where T <: AbstractGraph
    gnv = nv(g)
    hnv = nv(h)

    r = T(max(gnv, hnv))
    for e in edges(g)
        !has_edge(h, e) && add_edge!(r, e)
    end
    for e in edges(h)
        !has_edge(g, e) && add_edge!(r, e)
    end
    return r
end

"""
    union(g, h)

Return a graph that combines graphs `g` and `h` by taking the set union
of all vertices and edges.

### Implementation Notes
Preserves the eltype of the input graph. Will error if the
number of vertices in the generated graph exceeds the eltype.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = SimpleGraph(3); h = SimpleGraph(5);

julia> add_edge!(g, 1, 2);

julia> add_edge!(g, 1, 3);

julia> add_edge!(h, 3, 4);

julia> add_edge!(h, 3, 5);

julia> add_edge!(h, 4, 5);

julia> f = union(g, h);

julia> collect(edges(f))
5-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 1 => 3
 Edge 3 => 4
 Edge 3 => 5
 Edge 4 => 5
```
"""
function union(g::T, h::T) where T <: AbstractSimpleGraph
    gnv = nv(g)
    hnv = nv(h)

    r = T(max(gnv, hnv))
    r.ne = ne(g)
    for i in vertices(g)
        r.fadjlist[i] = deepcopy(g.fadjlist[i])
        if is_directed(g)
            r.badjlist[i] = deepcopy(g.badjlist[i])
        end
    end
    for e in edges(h)
        add_edge!(r, e)
    end
    return r
end

"""
    join(g, h)

Return a graph that combines graphs `g` and `h` using `blockdiag` and then
adds all the edges between the vertices in `g` and those in `h`.

### Implementation Notes
Preserves the eltype of the input graph. Will error if the number of vertices
in the generated graph exceeds the eltype.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = join(StarGraph(3), PathGraph(2))
{5, 9} undirected simple Int64 graph

julia> collect(edges(g))
9-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 1 => 3
 Edge 1 => 4
 Edge 1 => 5
 Edge 2 => 4
 Edge 2 => 5
 Edge 3 => 4
 Edge 3 => 5
 Edge 4 => 5
```
"""
function join(g::T, h::T) where T <: AbstractGraph
    r = blockdiag(g, h)
    for i in vertices(g)
        for j = (nv(g) + 1):(nv(g) + nv(h))
            add_edge!(r, i, j)
        end
    end
    return r
end

"""
    crosspath(len::Integer, g::Graph)

Return a graph that duplicates `g` `len` times and connects each vertex
with its copies in a path.

### Implementation Notes
Preserves the eltype of the input graph. Will error if the number of vertices
in the generated graph exceeds the eltype.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = crosspath(3, PathGraph(3))
{9, 12} undirected simple Int64 graph

julia> collect(edges(g))
12-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 1 => 4
 Edge 2 => 3
 Edge 2 => 5
 Edge 3 => 6
 Edge 4 => 5
 Edge 4 => 7
 Edge 5 => 6
 Edge 5 => 8
 Edge 6 => 9
 Edge 7 => 8
 Edge 8 => 9
```
"""
function crosspath end
# see https://github.com/mauro3/SimpleTraits.jl/issues/47#issuecomment-327880153 for syntax
@traitfn function crosspath(len::Integer, g::AG::(!IsDirected)) where {T, AG <: AbstractGraph{T}}
    p = PathGraph(len)
    h = Graph{T}(p)
    return cartesian_product(h, g)
end

# The following operators allow one to use a LightGraphs.Graph as a matrix in eigensolvers for spectral ranking and partitioning.
# """Provides multiplication of a graph `g` by a vector `v` such that spectral
# graph functions in [GraphMatrices.jl](https://github.com/jpfairbanks/GraphMatrices.jl) can utilize LightGraphs natively.
# """
function *(g::Graph, v::Vector{T}) where T <: Real
    length(v) == nv(g) || throw(ArgumentError("Vector size must equal number of vertices"))
    y = zeros(T, nv(g))
    for e in edges(g)
        i = src(e)
        j = dst(e)
        y[i] += v[j]
        y[j] += v[i]
    end
    return y
end

function *(g::DiGraph, v::Vector{T}) where T <: Real
    length(v) == nv(g) || throw(ArgumentError("Vector size must equal number of vertices"))
    y = zeros(T, nv(g))
    for e in edges(g)
        i = src(e)
        j = dst(e)
        y[i] += v[j]
    end
    return y
end

"""
    sum(g, i)

Return a vector of indegree (`i`=1) or outdegree (`i`=2) values for graph `g`.

# Examples
```jldoctest
julia> g = SimpleDiGraph([0 1 0 0 0; 0 0 1 0 0; 1 0 0 1 0; 0 0 0 0 1; 0 0 0 1 0]);

julia> sum(g, 2)
5-element Array{Int64,1}:
 1
 1
 2
 1
 1

julia> sum(g, 1)
5-element Array{Int64,1}:
 1
 1
 1
 2
 1
```
"""
function sum(g::AbstractGraph, dim::Int)
    dim == 1 && return indegree(g, vertices(g))
    dim == 2 && return outdegree(g, vertices(g))
    throw(ArgumentError("dimension must be <= 2"))
end

size(g::AbstractGraph) = (nv(g), nv(g))
"""
    size(g, i)

Return the number of vertices in `g` if `i`=1 or `i`=2, or `1` otherwise.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = CycleGraph(4);

julia> size(g, 1)
4

julia> size(g, 2)
4

julia> size(g, 3)
1
```
"""
size(g::Graph, dim::Int) = (dim == 1 || dim == 2) ? nv(g) : 1

"""
    sum(g)

Return the number of edges in `g`.

# Examples
```jldoctest
julia> g = SimpleGraph([0 1 0; 1 0 1; 0 1 0]);

julia> sum(g)
2
```
"""
sum(g::AbstractGraph) = ne(g)

"""
    sparse(g)

Return the default adjacency matrix of `g`.
"""
sparse(g::AbstractGraph) = adjacency_matrix(g)

length(g::AbstractGraph) = nv(g) * nv(g)
ndims(g::AbstractGraph) = 2
issymmetric(g::AbstractGraph) = !is_directed(g)

"""
    cartesian_product(g, h)

Return the (cartesian product)[https://en.wikipedia.org/wiki/Cartesian_product_of_graphs]
of `g` and `h`.

### Implementation Notes
Preserves the eltype of the input graph. Will error if the number of vertices
in the generated graph exceeds the eltype.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = cartesian_product(StarGraph(3), PathGraph(3))
{9, 12} undirected simple Int64 graph

julia> collect(edges(g))
12-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 1 => 4
 Edge 1 => 7
 Edge 2 => 3
 Edge 2 => 5
 Edge 2 => 8
 Edge 3 => 6
 Edge 3 => 9
 Edge 4 => 5
 Edge 5 => 6
 Edge 7 => 8
 Edge 8 => 9
```
"""
function cartesian_product(g::G, h::G) where G <: AbstractGraph
    z = G(nv(g) * nv(h))
    id(i, j) = (i - 1) * nv(h) + j
    for e in edges(g)
        i1, i2 = Tuple(e)
        for j = 1:nv(h)
            add_edge!(z, id(i1, j), id(i2, j))
        end
    end

    for e in edges(h)
        j1, j2 = Tuple(e)
        for i in vertices(g)
            add_edge!(z, id(i, j1), id(i, j2))
        end
    end
    return z
end

"""
    tensor_product(g, h)

Return the (tensor product)[https://en.wikipedia.org/wiki/Tensor_product_of_graphs]
of `g` and `h`.

### Implementation Notes
Preserves the eltype of the input graph. Will error if the number of vertices
in the generated graph exceeds the eltype.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = tensor_product(StarGraph(3), PathGraph(3))
{9, 8} undirected simple Int64 graph

julia> collect(edges(g))
8-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 5
 Edge 1 => 8
 Edge 2 => 4
 Edge 2 => 6
 Edge 2 => 7
 Edge 2 => 9
 Edge 3 => 5
 Edge 3 => 8
```
"""
function tensor_product(g::G, h::G) where G <: AbstractGraph
    z = G(nv(g) * nv(h))
    id(i, j) = (i - 1) * nv(h) + j
    undirected = !is_directed(g)
    for e1 in edges(g)
        i1, i2 = Tuple(e1)
        for e2 in edges(h)
            j1, j2 = Tuple(e2)
            add_edge!(z, id(i1, j1), id(i2, j2))
            if undirected
                add_edge!(z, id(i1, j2), id(i2, j1))
            end
        end
    end
    return z
end

## subgraphs ###

"""
    induced_subgraph(g, vlist)
    induced_subgraph(g, elist)

Return the subgraph of `g` induced by the vertices in  `vlist` or edges in `elist`
along with a vector mapping the new vertices to the old ones
(the  vertex `i` in the subgraph corresponds to the vertex `vmap[i]` in `g`.)

The returned graph has `length(vlist)` vertices, with the new vertex `i`
corresponding to the vertex of the original graph in the `i`-th position
of `vlist`.

### Usage Examples
```doctestjl
julia> g = CompleteGraph(10)

julia> sg, vmap = induced_subgraph(g, 5:8)

julia> @assert g[5:8] == sg

julia> @assert nv(sg) == 4

julia> @assert ne(sg) == 6

julia> @assert vm[4] == 8

julia> sg, vmap = induced_subgraph(g, [2,8,3,4])

julia> @assert sg == g[[2,8,3,4]]

julia> elist = [Edge(1,2), Edge(3,4), Edge(4,8)]

julia> sg, vmap = induced_subgraph(g, elist)

julia> @assert sg == g[elist]
```
"""
function induced_subgraph(g::T, vlist::AbstractVector{U}) where T <: AbstractGraph where U <: Integer
    allunique(vlist) || throw(ArgumentError("Vertices in subgraph list must be unique"))
    h = T(length(vlist))
    newvid = Dict{U,U}()
    vmap = Vector{U}(undef, length(vlist))
    for (i, v) in enumerate(vlist)
        newvid[v] = U(i)
        vmap[i] = v
    end

    vset = Set(vlist)
    for s in vlist
        for d in outneighbors(g, s)
            # println("s = $s, d = $d")
            if d in vset && has_edge(g, s, d)
                newe = Edge(newvid[s], newvid[d])
                add_edge!(h, newe)
            end
        end
    end
    return h, vmap
end

function induced_subgraph(g::AG, elist::AbstractVector{U}) where AG <: AbstractGraph{T} where T where U <: AbstractEdge
    h = zero(g)
    newvid = Dict{T,T}()
    vmap = Vector{T}()

    for e in elist
        u, v = Tuple(e)
        for i in (u, v)
            if !haskey(newvid, i)
                add_vertex!(h)
                newvid[i] = nv(h)
                push!(vmap, i)
            end
        end
        add_edge!(h, newvid[u], newvid[v])
    end
    return h, vmap
end

"""
    g[iter]

Return the subgraph induced by `iter`.
Equivalent to [`induced_subgraph`](@ref)`(g, iter)[1]`.
"""
getindex(g::AbstractGraph, iter) = induced_subgraph(g, iter)[1]

"""
    egonet(g, v, d, distmx=weights(g))

Return the subgraph of `g` induced by the neighbors of `v` up to distance
`d`, using weights (optionally) provided by `distmx`.
This is equivalent to [`induced_subgraph`](@ref)`(g, neighborhood(g, v, d, dir=dir))[1].`

### Optional Arguments
- `dir=:out`: if `g` is directed, this argument specifies the edge direction
with respect to `v` (i.e. `:in` or `:out`).
"""
egonet(g::AbstractGraph{T}, v::Integer, d::Integer, distmx::AbstractMatrix{U}=weights(g); dir=:out) where T <: Integer where U <: Real =
    g[neighborhood(g, v, d, distmx, dir=dir)]

"""
    compute_shifts(n::Int, x::AbstractArray)

Determine how many elements of `x` are less than `i` for all `i` in `1:n`.
"""
function compute_shifts(n::Integer, x::AbstractArray)
    tmp = zeros(eltype(x), n)
    tmp[x[2:end]] .= 1
    return cumsum!(tmp, tmp)
end

"""
    merge_vertices(g::AbstractGraph, vs)

Create a new graph where all vertices in `vs` have been aliased to the same vertex `minimum(vs)`.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = PathGraph(5);

julia> collect(edges(g))
4-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 2 => 3
 Edge 3 => 4
 Edge 4 => 5

julia> h = merge_vertices(g, [2, 3]);

julia> collect(edges(h))
3-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 2 => 3
 Edge 3 => 4
```
"""
function merge_vertices(g::AbstractGraph, vs)
    labels = collect(1:nv(g))
    # Use lowest value as new vertex id.
    sort!(vs)
    nvnew = nv(g) - length(unique(vs)) + 1
    nvnew <= nv(g) || return g
    (v0, vm) = extrema(vs)
    v0 > 0 || throw(ArgumentError("invalid vertex ID: $v0 in list of vertices to be merged"))
    vm <= nv(g) || throw(ArgumentError("vertex $vm not found in graph")) # TODO 0.7: change to DomainError?
    labels[vs] .= v0
    shifts = compute_shifts(nv(g), vs[2:end])
    for v in vertices(g)
        if labels[v] != v0
            labels[v] -= shifts[v]
        end
    end

    #if v in vs then labels[v] == v0 else labels[v] == v
    newg = SimpleGraph(nvnew)
    for e in edges(g)
        u, w = src(e), dst(e)
        if labels[u] != labels[w] #not a new self loop
            add_edge!(newg, labels[u], labels[w])
        end
    end
    return newg
end

"""
    merge_vertices!(g, vs)

Combine vertices specified in `vs` into single vertex whose
index will be the lowest value in `vs`. All edges connected to vertices in `vs`
connect to the new merged vertex.

Return a vector with new vertex values are indexed by the original vertex indices.

### Implementation Notes
Supports [`SimpleGraph`](@ref) only.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = PathGraph(5);

julia> collect(edges(g))
4-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 2 => 3
 Edge 3 => 4
 Edge 4 => 5

julia> merge_vertices!(g, [2, 3])
5-element Array{Int64,1}:
 1
 2
 2
 3
 4

julia> collect(edges(g))
3-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 2 => 3
 Edge 3 => 4
```
"""
function merge_vertices!(g::Graph{T}, vs::Vector{U} where U <: Integer) where T
    vs = sort!(unique(vs))
    merged_vertex = popfirst!(vs)

    x = zeros(Int, nv(g))
    x[vs] .= 1
    new_vertex_ids = collect(1:nv(g)) .- cumsum(x)
    new_vertex_ids[vs] .= merged_vertex

    for i in vertices(g)
        # Adjust connections to merged vertices
        if (i != merged_vertex) && !insorted(i, vs)
            nbrs_to_rewire = Set{T}()
            for j in outneighbors(g, i)
                if insorted(j, vs)
                    push!(nbrs_to_rewire, merged_vertex)
                else
                    push!(nbrs_to_rewire, new_vertex_ids[j])
                end
            end
            g.fadjlist[new_vertex_ids[i]] = sort(collect(nbrs_to_rewire))

        # Collect connections to new merged vertex
        else
            nbrs_to_merge = Set{T}()
            for element in filter(x -> !(insorted(x, vs)) && (x != merged_vertex), g.fadjlist[i])
                push!(nbrs_to_merge, new_vertex_ids[element])
            end

            for j in vs, e in outneighbors(g, j)
                if new_vertex_ids[e] != merged_vertex
                    push!(nbrs_to_merge, new_vertex_ids[e])
                end
            end
            g.fadjlist[i] = sort(collect(nbrs_to_merge))
        end
    end

    # Drop excess vertices
    g.fadjlist = g.fadjlist[1:(end - length(vs))]

    # Correct edge counts
    g.ne = sum(degree(g, i) for i in vertices(g)) / 2

    return new_vertex_ids
end


text:
using Optim

# I made a pull request for this change to be made in Optim. We'll see where that goes. If
#   it doesn't get merged, I'll have to make Optim a dependency and put this in src.
function Optim.initial_state(method::IPNewton, options, d::TwiceDifferentiable, constraints::TwiceDifferentiableConstraints, initial_x::ModelingStruct{T}) where T
    # Check feasibility of the initial state
    mc = Optim.nconstraints(constraints)
    constr_c = Array{T}(undef, mc)
    # TODO: When we change to `value!` from NLSolversBase instead of c!
    # we can also update `initial_convergence` for ConstrainedOptimizer in interior.jl
    constraints.c!(constr_c, initial_x)
    if !Optim.isinterior(constraints, initial_x, constr_c)
        @warn("Initial guess is not an interior point")
        Base.show_backtrace(stderr, backtrace())
        println(stderr)
    end
    # Allocate fields for the objective function
    n = length(initial_x)
    g = similar(initial_x)
    s = similar(initial_x)
    f_x_previous = NaN
    f_x, g_x = Optim.value_gradient!(d, initial_x)
    g .= g_x # needs to be a separate copy of g_x
    H = Matrix{T}(undef, n, n)
    Hd = Vector{Int8}(undef, n)
    Optim.hessian!(d, initial_x)
    copyto!(H, Optim.hessian(d))

    # More constraints
    constr_J = Array{T}(undef, mc, n)
    gtilde = copy(g)
    constraints.jacobian!(constr_J, initial_x)
    μ = T(1)
    bstate = Optim.BarrierStateVars(constraints.bounds, initial_x, constr_c)
    bgrad = copy(bstate)
    bstep = copy(bstate)
    # b_ls = BarrierLineSearch(similar(constr_c), similar(bstate))
    b_ls = Optim.BarrierLineSearchGrad(copy(constr_c), copy(constr_J), copy(bstate), copy(bstate))

    state = Optim.IPNewtonState(
        copy(initial_x), # Maintain current state in state.x
        f_x, # Store current f in state.f_x
        copy(initial_x), # Maintain previous state in state.x_previous
        g, # Store current gradient in state.g (TODO: includes Lagrangian calculation?)
        T(NaN), # Store previous f in state.f_x_previous
        H,
        0,    # will be replaced
        Hd,
        similar(initial_x), # Maintain current x-search direction in state.s
        μ,
        μ,
        T(NaN),
        T(NaN),
        bstate,
        bgrad,
        bstep,
        constr_c,
        constr_J,
        T(NaN),
        Optim.@initial_linesearch()..., # Maintain a cache for line search results in state.lsr
        b_ls,
        gtilde,
        0)

    Hinfo = (state.H, Optim.hessianI(initial_x, constraints, 1 ./ bstate.slack_c, 1))
    Optim.initialize_μ_λ!(state, constraints.bounds, Hinfo, method.μ0)
    Optim.update_fg!(d, constraints, state, method)
    Optim.update_h!(d, constraints, state, method)

    state
end

max_stars_repo_name: jonniedie/NamedViewVectors


text:
function parametercolors()
    α_color = :red
    ωₚ_color = :navy
    γ_color = :orange
    r_color = :purple
    ϕₘ_color = :blue
    β_color = :green
    ν_color = :magenta
    σₗ_color = :maroon
    σᵣ_color = :teal
    return [α_color,
            ωₚ_color,
            γ_color,
            r_color,
            ϕₘ_color,
            β_color,
            ν_color,
            σₗ_color,
            σᵣ_color]
end
function makeparameternames()
    return [
        "\\alpha",
        "\\omega_p",
        "\\gamma",
        "r",
        "\\phi_m",
        "\\beta",
        "\\nu",
        "\\sigma_l",
        "\\sigma_r"
    ]
end

mod2pi_nan(x) = isnan(x) ? x : mod2pi(x)


text:
# The composite operators are built using basic operators (scalar, array, and
# derivative) using arithmetic or other operator compositions. The composite
# operator types are lazy and maintain the structure used to build them.

# Recursive routines that use `getops`
function update_coefficients!(L::AbstractDiffEqCompositeOperator,u,p,t)
  for op in getops(L)
    update_coefficients!(op,u,p,t)
  end
  L
end
isconstant(L::AbstractDiffEqCompositeOperator) = all(isconstant, getops(L))

# Scaled operator (α * A)
struct DiffEqScaledOperator{T,F,OpType<:AbstractDiffEqLinearOperator{T}} <: AbstractDiffEqCompositeOperator{T}
  coeff::DiffEqScalar{T,F}
  op::OpType
end
*(α::DiffEqScalar{T,F}, L::AbstractDiffEqLinearOperator{T}) where {T,F} = DiffEqScaledOperator(α, L)
*(α::Number, L::AbstractDiffEqLinearOperator{T}) where T = DiffEqScaledOperator(DiffEqScalar(convert(T,α)), L)
-(L::AbstractDiffEqLinearOperator{T}) where {T} = DiffEqScalar(-one(T)) * L
getops(L::DiffEqScaledOperator) = (L.coeff, L.op)
Matrix(L::DiffEqScaledOperator) = L.coeff * Matrix(L.op)
convert(::Type{AbstractMatrix}, L::DiffEqScaledOperator) = L.coeff * convert(AbstractMatrix, L.op)

size(L::DiffEqScaledOperator, args...) = size(L.op, args...)
opnorm(L::DiffEqScaledOperator, p::Real=2) = abs(L.coeff) * opnorm(L.op, p)
getindex(L::DiffEqScaledOperator, i::Int) = L.coeff * L.op[i]
getindex(L::DiffEqScaledOperator, I::Vararg{Int, N}) where {N} =
  L.coeff * L.op[I...]
*(L::DiffEqScaledOperator, x::AbstractArray) = L.coeff * (L.op * x)
*(x::AbstractArray, L::DiffEqScaledOperator) = (L.op * x) * L.coeff
/(L::DiffEqScaledOperator, x::AbstractArray) = L.coeff * (L.op / x)
/(x::AbstractArray, L::DiffEqScaledOperator) = 1/L.coeff * (x / L.op)
\(L::DiffEqScaledOperator, x::AbstractArray) = 1/L.coeff * (L.op \ x)
\(x::AbstractArray, L::DiffEqScaledOperator) = L.coeff * (x \ L)
for N in (2,3)
  @eval begin
    mul!(Y::AbstractArray{T,$N}, L::DiffEqScaledOperator{T}, B::AbstractArray{T,$N}) where {T} =
        lmul!(Y, L.coeff, mul!(Y, L.op, B))
  end
end
ldiv!(Y::AbstractArray, L::DiffEqScaledOperator, B::AbstractArray) =
  lmul!(1/L.coeff, ldiv!(Y, L.op, B))
factorize(L::DiffEqScaledOperator) = L.coeff * factorize(L.op)
for fact in (:lu, :lu!, :qr, :qr!, :cholesky, :cholesky!, :ldlt, :ldlt!,
  :bunchkaufman, :bunchkaufman!, :lq, :lq!, :svd, :svd!)
  @eval LinearAlgebra.$fact(L::DiffEqScaledOperator, args...) =
    L.coeff * fact(L.op, args...)
end

# Linear Combination
struct DiffEqOperatorCombination{T,O<:Tuple{Vararg{AbstractDiffEqLinearOperator{T}}},
    C<:AbstractVector{T}} <: AbstractDiffEqCompositeOperator{T}
    ops::O
    cache::C
    function DiffEqOperatorCombination(ops; cache=nothing)
        T = eltype(ops[1])
        for i in 2:length(ops)
            @assert size(ops[i]) == size(ops[1]) "Operators must be of the same size to be combined! Mismatch between $(ops[i]) and $(ops[i-1]), which are operators $i and $(i-1) respectively"
        end
        if cache == nothing
            cache = zeros(T, size(ops[1], 1))
        end
        new{T,typeof(ops),typeof(cache)}(ops, cache)
    end
end
+(ops::AbstractDiffEqLinearOperator...) = DiffEqOperatorCombination(ops)
+(L1::DiffEqOperatorCombination, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorCombination((L1.ops..., L2))
+(L1::AbstractDiffEqLinearOperator, L2::DiffEqOperatorCombination) = DiffEqOperatorCombination((L1, L2.ops...))
+(L1::DiffEqOperatorCombination, L2::DiffEqOperatorCombination) = DiffEqOperatorCombination((L1.ops..., L2.ops...))
-(L1::AbstractDiffEqLinearOperator, L2::AbstractDiffEqLinearOperator) = L1 + (-L2)
getops(L::DiffEqOperatorCombination) = L.ops
Matrix(L::DiffEqOperatorCombination) = sum(Matrix, L.ops)
convert(::Type{AbstractMatrix}, L::DiffEqOperatorCombination) =
  sum(op -> convert(AbstractMatrix, op), L.ops)

size(L::DiffEqOperatorCombination, args...) = size(L.ops[1], args...)
getindex(L::DiffEqOperatorCombination, i::Int) = sum(op -> op[i], L.ops)
getindex(L::DiffEqOperatorCombination, I::Vararg{Int, N}) where {N} = sum(op -> op[I...], L.ops)
*(L::DiffEqOperatorCombination, x::AbstractArray) = sum(op -> op * x, L.ops)
*(x::AbstractArray, L::DiffEqOperatorCombination) = sum(op -> x * op, L.ops)
/(L::DiffEqOperatorCombination, x::AbstractArray) = sum(op -> op / x, L.ops)
\(x::AbstractArray, L::DiffEqOperatorCombination) = sum(op -> x \ op, L.ops)
function mul!(y::AbstractVector, L::DiffEqOperatorCombination, b::AbstractVector)
  mul!(y, L.ops[1], b)
  for op in L.ops[2:end]
    mul!(L.cache, op, b)
    y .+= L.cache
  end
  return y
end

# Composition (A ∘ B)
struct DiffEqOperatorComposition{T,O<:Tuple{Vararg{AbstractDiffEqLinearOperator{T}}},
  C<:Tuple{Vararg{AbstractVector{T}}}} <: AbstractDiffEqCompositeOperator{T}
  ops::O # stored in the order of application
  caches::C
  function DiffEqOperatorComposition(ops; caches=nothing)
    T = eltype(ops[1])
    for i in 2:length(ops)
      @assert size(ops[i-1], 1) == size(ops[i], 2) "Operations do not have compatible sizes! Mismatch between $(ops[i]) and $(ops[i-1]), which are operators $i and $(i-1) respectively."
    end

    if caches == nothing
      # Construct a list of caches to be used by mul! and ldiv!
      caches = []
      for op in ops[1:end-1]
        tmp = Vector{T}(undef, size(op, 1))
        fill!(tmp,0)
        push!(caches,tmp)
      end
      caches = tuple(caches...)
    end
    new{T,typeof(ops),typeof(caches)}(ops, caches)
  end
end
*(ops::AbstractDiffEqLinearOperator...) = DiffEqOperatorComposition(reverse(ops))
∘(L1::AbstractDiffEqLinearOperator, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorComposition((L2, L1))
*(L1::DiffEqOperatorComposition, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorComposition((L2, L1.ops...))
∘(L1::DiffEqOperatorComposition, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorComposition((L2, L1.ops...))
*(L1::AbstractDiffEqLinearOperator, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1))
∘(L1::AbstractDiffEqLinearOperator, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1))
*(L1::DiffEqOperatorComposition, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1.ops...))
∘(L1::DiffEqOperatorComposition, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1.ops...))
getops(L::DiffEqOperatorComposition) = L.ops
Matrix(L::DiffEqOperatorComposition) = prod(Matrix, reverse(L.ops))
convert(::Type{AbstractMatrix}, L::DiffEqOperatorComposition) =
  prod(op -> convert(AbstractMatrix, op), reverse(L.ops))

size(L::DiffEqOperatorComposition) = (size(L.ops[end], 1), size(L.ops[1], 2))
size(L::DiffEqOperatorComposition, m::Integer) = size(L)[m]
opnorm(L::DiffEqOperatorComposition) = prod(opnorm, L.ops)
*(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op*acc, L.ops; init=x)
*(x::AbstractArray, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc*op, reverse(L.ops); init=x)
/(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op/acc, L.ops; init=x)
/(x::AbstractArray, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc/op, L.ops; init=x)
\(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op\acc, reverse(L.ops); init=x)
\(x::AbstractArray, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc\op, reverse(L.ops); init=x)
function mul!(y::AbstractVector, L::DiffEqOperatorComposition, b::AbstractVector)
  mul!(L.caches[1], L.ops[1], b)
  for i in 2:length(L.ops) - 1
    mul!(L.caches[i], L.ops[i], L.caches[i-1])
  end
  mul!(y, L.ops[end], L.caches[end])
end
function ldiv!(y::AbstractVector, L::DiffEqOperatorComposition, b::AbstractVector)
  ldiv!(L.caches[end], L.ops[end], b)
  for i in length(L.ops) - 1:-1:2
    ldiv!(L.caches[i-1], L.ops[i], L.caches[i])
  end
  ldiv!(y, L.ops[1], L.caches[1])
end
factorize(L::DiffEqOperatorComposition) = prod(factorize, reverse(L.ops))
for fact in (:lu, :lu!, :qr, :qr!, :cholesky, :cholesky!, :ldlt, :ldlt!,
  :bunchkaufman, :bunchkaufman!, :lq, :lq!, :svd, :svd!)
  @eval LinearAlgebra.$fact(L::DiffEqOperatorComposition, args...) =
    prod(op -> $fact(op, args...), reverse(L.ops))
end


text:
#ДАНО: Робот - в произвольной клетке, где-то между внутренними перегородкаи прямоугольной формы и
#внешней рамкой поля (прямоугольники могут быть и вырожденными, т.е. представлять собой отрезки)
#РЕЗУЛЬТАТ: Робот - в исходном положении и в 4-х приграничных клетках,
#2-е из которых имеют ту же широту, а 2-е - ту же долготу, что и Робот, стоят маркеры.
function mark_centers(r)
    num_steps = through_rectangles_into_angle(r,(Sud,West))
    num_steps_to_ost = sum(num_steps[1:2:end])
    num_steps_to_nord = sum(num_steps[2:2:end])

    movements!(r,Nord,num_steps_to_nord)
    putmarker!(r)
    num_steps_to_sud = movements!(r,Nord)

    movements!(r,Ost,num_steps_to_ost)
    putmarker!(r)
    num_steps_to_west = movements!(r,Ost)

    movements!(r,Sud,num_steps_to_sud)
    putmarker!(r)
    movements!(r,Sud)

    movements!(r,Sud,num_steps_to_west)
    putmarker!(r)
    movements!(r,Sud)

    movements!(r,(Ost,Nord),num_steps)

end
function through_rectangles_into_angle(r,angle::NTuple{2,HorizonSide})
    num_steps=[]
    while isborder(r,angle[1])==false || isborder(r,angle[2]) # Робот - не в юго-западном углу
        push!(num_steps, movements!(r, angle[2]))
        push!(num_steps, movements!(r, angle[1]))
    end
    return num_steps
end
function movements!(r,sides,num_steps::Vector{Int})
    for (i,n) in enumerate(num_steps)
        movements!(r, sides[mod(i-1, length(sides))+1], n)
    end
end