text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

isexpr(@nospecialize(ex), heads) = isa(ex, Expr) && in(ex.head, heads)
isexpr(@nospecialize(ex), heads, n::Int) = isa(ex, Expr) && in(ex.head, heads) && length(ex.args) == n
const is_expr = isexpr

## symbols ##

"""
    gensym([tag])

Generates a symbol which will not conflict with other variable names (in the same module).
"""
gensym() = ccall(:jl_gensym, Ref{Symbol}, ())

gensym(s::String) = ccall(:jl_tagged_gensym, Ref{Symbol}, (Ptr{UInt8}, Csize_t), s, sizeof(s))

gensym(ss::String...) = map(gensym, ss)
gensym(s::Symbol) = ccall(:jl_tagged_gensym, Ref{Symbol}, (Ptr{UInt8}, Csize_t), s, -1 % Csize_t)

"""
    @gensym

Generates a gensym symbol for a variable. For example, `@gensym x y` is transformed into
`x = gensym("x"); y = gensym("y")`.
"""
macro gensym(names...)
    blk = Expr(:block)
    for name in names
        push!(blk.args, :($(esc(name)) = gensym($(string(name)))))
    end
    push!(blk.args, :nothing)
    return blk
end

## expressions ##

isexpr(@nospecialize(ex), head::Symbol) = isa(ex, Expr) && ex.head === head
isexpr(@nospecialize(ex), head::Symbol, n::Int) = isa(ex, Expr) && ex.head === head && length(ex.args) == n

copy(e::Expr) = exprarray(e.head, copy_exprargs(e.args))

# copy parts of an AST that the compiler mutates
function copy_exprs(@nospecialize(x))
    if isa(x, Expr)
        return copy(x)
    elseif isa(x, PhiNode)
        values = x.values
        nvalues = length(values)
        new_values = Vector{Any}(undef, nvalues)
        @inbounds for i = 1:nvalues
            isassigned(values, i) || continue
            new_values[i] = copy_exprs(values[i])
        end
        return PhiNode(copy(x.edges), new_values)
    elseif isa(x, PhiCNode)
        values = x.values
        nvalues = length(values)
        new_values = Vector{Any}(undef, nvalues)
        @inbounds for i = 1:nvalues
            isassigned(values, i) || continue
            new_values[i] = copy_exprs(values[i])
        end
        return PhiCNode(new_values)
    end
    return x
end
copy_exprargs(x::Array{Any,1}) = Any[copy_exprs(@inbounds x[i]) for i in 1:length(x)]

@eval exprarray(head::Symbol, arg::Array{Any,1}) = $(Expr(:new, :Expr, :head, :arg))

# create copies of the CodeInfo definition, and any mutable fields
function copy(c::CodeInfo)
    cnew = ccall(:jl_copy_code_info, Ref{CodeInfo}, (Any,), c)
    cnew.code = copy_exprargs(cnew.code)
    cnew.slotnames = copy(cnew.slotnames)
    cnew.slotflags = copy(cnew.slotflags)
    cnew.codelocs  = copy(cnew.codelocs)
    cnew.linetable = copy(cnew.linetable::Union{Vector{Any},Vector{Core.LineInfoNode}})
    cnew.ssaflags  = copy(cnew.ssaflags)
    cnew.edges     = cnew.edges === nothing ? nothing : copy(cnew.edges::Vector)
    ssavaluetypes  = cnew.ssavaluetypes
    ssavaluetypes isa Vector{Any} && (cnew.ssavaluetypes = copy(ssavaluetypes))
    return cnew
end

==(x::Expr, y::Expr) = x.head === y.head && isequal(x.args, y.args)
==(x::QuoteNode, y::QuoteNode) = isequal(x.value, y.value)
==(stmt1::Core.PhiNode, stmt2::Core.PhiNode) = stmt1.edges == stmt2.edges && stmt1.values == stmt2.values

"""
    macroexpand(m::Module, x; recursive=true)

Take the expression `x` and return an equivalent expression with all macros removed (expanded)
for executing in module `m`.
The `recursive` keyword controls whether deeper levels of nested macros are also expanded.
This is demonstrated in the example below:
```julia-repl
julia> module M
           macro m1()
               42
           end
           macro m2()
               :(@m1())
           end
       end
M

julia> macroexpand(M, :(@m2()), recursive=true)
42

julia> macroexpand(M, :(@m2()), recursive=false)
:(#= REPL[16]:6 =# M.@m1)
```
"""
function macroexpand(m::Module, @nospecialize(x); recursive=true)
    if recursive
        ccall(:jl_macroexpand, Any, (Any, Any), x, m)
    else
        ccall(:jl_macroexpand1, Any, (Any, Any), x, m)
    end
end

"""
    @macroexpand

Return equivalent expression with all macros removed (expanded).

There are differences between `@macroexpand` and [`macroexpand`](@ref).

* While [`macroexpand`](@ref) takes a keyword argument `recursive`, `@macroexpand`
  is always recursive. For a non recursive macro version, see [`@macroexpand1`](@ref).

* While [`macroexpand`](@ref) has an explicit `module` argument, `@macroexpand` always
  expands with respect to the module in which it is called.

This is best seen in the following example:
```julia-repl
julia> module M
           macro m()
               1
           end
           function f()
               (@macroexpand(@m),
                macroexpand(M, :(@m)),
                macroexpand(Main, :(@m))
               )
           end
       end
M

julia> macro m()
           2
       end
@m (macro with 1 method)

julia> M.f()
(1, 1, 2)
```
With `@macroexpand` the expression expands where `@macroexpand` appears in the code (module `M` in the example).
With `macroexpand` the expression expands in the module given as the first argument.
"""
macro macroexpand(code)
    return :(macroexpand($__module__, $(QuoteNode(code)), recursive=true))
end

"""
    @macroexpand1

Non recursive version of [`@macroexpand`](@ref).
"""
macro macroexpand1(code)
    return :(macroexpand($__module__, $(QuoteNode(code)), recursive=false))
end

## misc syntax ##

"""
    Core.eval(m::Module, expr)

Evaluate an expression in the given module and return the result.
"""
Core.eval

"""
    @inline

Give a hint to the compiler that this function is worth inlining.

Small functions typically do not need the `@inline` annotation,
as the compiler does it automatically. By using `@inline` on bigger functions,
an extra nudge can be given to the compiler to inline it.

`@inline` can be applied immediately before the definition or in its function body.

```julia
# annotate long-form definition
@inline function longdef(x)
    ...
end

# annotate short-form definition
@inline shortdef(x) = ...

# annotate anonymous function that a `do` block creates
f() do
    @inline
    ...
end
```

!!! compat "Julia 1.8"
    The usage within a function body requires at least Julia 1.8.

---
    @inline block

Give a hint to the compiler that calls within `block` are worth inlining.

```julia
# The compiler will try to inline `f`
@inline f(...)

# The compiler will try to inline `f`, `g` and `+`
@inline f(...) + g(...)
```

!!! note
    A callsite annotation always has the precedence over the annotation applied to the
    definition of the called function:
    ```julia
    @noinline function explicit_noinline(args...)
        # body
    end

    let
        @inline explicit_noinline(args...) # will be inlined
    end
    ```

!!! note
    When there are nested callsite annotations, the innermost annotation has the precedence:
    ```julia
    @noinline let a0, b0 = ...
        a = @inline f(a0)  # the compiler will try to inline this call
        b = f(b0)          # the compiler will NOT try to inline this call
        return a, b
    end
    ```

!!! warning
    Although a callsite annotation will try to force inlining in regardless of the cost model,
    there are still chances it can't succeed in it. Especially, recursive calls can not be
    inlined even if they are annotated as `@inline`d.

!!! compat "Julia 1.8"
    The callsite annotation requires at least Julia 1.8.
"""
macro inline(x)
    return annotate_meta_def_or_block(x, :inline)
end

"""
    @noinline

Give a hint to the compiler that it should not inline a function.

Small functions are typically inlined automatically.
By using `@noinline` on small functions, auto-inlining can be
prevented.

`@noinline` can be applied immediately before the definition or in its function body.

```julia
# annotate long-form definition
@noinline function longdef(x)
    ...
end

# annotate short-form definition
@noinline shortdef(x) = ...

# annotate anonymous function that a `do` block creates
f() do
    @noinline
    ...
end
```

!!! compat "Julia 1.8"
    The usage within a function body requires at least Julia 1.8.

---
    @noinline block

Give a hint to the compiler that it should not inline the calls within `block`.

```julia
# The compiler will try to not inline `f`
@noinline f(...)

# The compiler will try to not inline `f`, `g` and `+`
@noinline f(...) + g(...)
```

!!! note
    A callsite annotation always has the precedence over the annotation applied to the
    definition of the called function:
    ```julia
    @inline function explicit_inline(args...)
        # body
    end

    let
        @noinline explicit_inline(args...) # will not be inlined
    end
    ```

!!! note
    When there are nested callsite annotations, the innermost annotation has the precedence:
    ```julia
    @inline let a0, b0 = ...
        a = @noinline f(a0)  # the compiler will NOT try to inline this call
        b = f(b0)            # the compiler will try to inline this call
        return a, b
    end
    ```

!!! compat "Julia 1.8"
    The callsite annotation requires at least Julia 1.8.

---
!!! note
    If the function is trivial (for example returning a constant) it might get inlined anyway.
"""
macro noinline(x)
    return annotate_meta_def_or_block(x, :noinline)
end

"""
    @pure ex

`@pure` gives the compiler a hint for the definition of a pure function,
helping for type inference.

!!! warning
    This macro is intended for internal compiler use and may be subject to changes.

!!! warning
    In Julia 1.8 and higher, it is favorable to use [`@assume_effects`](@ref) instead of `@pure`.
    This is because `@assume_effects` allows a finer grained control over Julia's purity
    modeling and the effect system enables a wider range of optimizations.
"""
macro pure(ex)
    esc(isa(ex, Expr) ? pushmeta!(ex, :pure) : ex)
end

"""
    @constprop setting ex

`@constprop` controls the mode of interprocedural constant propagation for the
annotated function. Two `setting`s are supported:

- `@constprop :aggressive ex`: apply constant propagation aggressively.
  For a method where the return type depends on the value of the arguments,
  this can yield improved inference results at the cost of additional compile time.
- `@constprop :none ex`: disable constant propagation. This can reduce compile
  times for functions that Julia might otherwise deem worthy of constant-propagation.
  Common cases are for functions with `Bool`- or `Symbol`-valued arguments or keyword arguments.
"""
macro constprop(setting, ex)
    if isa(setting, QuoteNode)
        setting = setting.value
    end
    setting === :aggressive && return esc(isa(ex, Expr) ? pushmeta!(ex, :aggressive_constprop) : ex)
    setting === :none && return esc(isa(ex, Expr) ? pushmeta!(ex, :no_constprop) : ex)
    throw(ArgumentError("@constprop $setting not supported"))
end

"""
    @assume_effects setting... ex

`@assume_effects` overrides the compiler's effect modeling for the given method.
`ex` must be a method definition or `@ccall` expression.

!!! compat "Julia 1.8"
    Using `Base.@assume_effects` requires Julia version 1.8.

```jldoctest
julia> Base.@assume_effects :terminates_locally function pow(x)
           # this :terminates_locally allows `pow` to be constant-folded
           res = 1
           1 < x < 20 || error("bad pow")
           while x > 1
               res *= x
               x -= 1
           end
           return res
       end
pow (generic function with 1 method)

julia> code_typed() do
           pow(12)
       end
1-element Vector{Any}:
 CodeInfo(
1 ─     return 479001600
) => Int64

julia> Base.@assume_effects :total !:nothrow @ccall jl_type_intersection(Vector{Int}::Any, Vector{<:Integer}::Any)::Any
Vector{Int64} (alias for Array{Int64, 1})
```

!!! warning
    Improper use of this macro causes undefined behavior (including crashes,
    incorrect answers, or other hard to track bugs). Use with care and only if
    absolutely required.

In general, each `setting` value makes an assertion about the behavior of the
function, without requiring the compiler to prove that this behavior is indeed
true. These assertions are made for all world ages. It is thus advisable to limit
the use of generic functions that may later be extended to invalidate the
assumption (which would cause undefined behavior).

The following `setting`s are supported.
- `:consistent`
- `:effect_free`
- `:nothrow`
- `:terminates_globally`
- `:terminates_locally`
- `:notaskstate`
- `:foldable`
- `:total`

# Extended help

---
# `:consistent`

The `:consistent` setting asserts that for egal (`===`) inputs:
- The manner of termination (return value, exception, non-termination) will always be the same.
- If the method returns, the results will always be egal.

!!! note
    This in particular implies that the return value of the method must be
    immutable. Multiple allocations of mutable objects (even with identical
    contents) are not egal.

!!! note
    The `:consistent`-cy assertion is made world-age wise. More formally, write
    ``fᵢ`` for the evaluation of ``f`` in world-age ``i``, then we require:
    ```math
    ∀ i, x, y: x ≡ y → fᵢ(x) ≡ fᵢ(y)
    ```
    However, for two world ages ``i``, ``j`` s.t. ``i ≠ j``, we may have ``fᵢ(x) ≢ fⱼ(y)``.

    A further implication is that `:consistent` functions may not make their
    return value dependent on the state of the heap or any other global state
    that is not constant for a given world age.

!!! note
    The `:consistent`-cy includes all legal rewrites performed by the optimizer.
    For example, floating-point fastmath operations are not considered `:consistent`,
    because the optimizer may rewrite them causing the output to not be `:consistent`,
    even for the same world age (e.g. because one ran in the interpreter, while
    the other was optimized).

!!! note
    If `:consistent` functions terminate by throwing an exception, that exception
    itself is not required to meet the egality requirement specified above.

---
# `:effect_free`

The `:effect_free` setting asserts that the method is free of externally semantically
visible side effects. The following is an incomplete list of externally semantically
visible side effects:
- Changing the value of a global variable.
- Mutating the heap (e.g. an array or mutable value), except as noted below
- Changing the method table (e.g. through calls to eval)
- File/Network/etc. I/O
- Task switching

However, the following are explicitly not semantically visible, even if they
may be observable:
- Memory allocations (both mutable and immutable)
- Elapsed time
- Garbage collection
- Heap mutations of objects whose lifetime does not exceed the method (i.e.
  were allocated in the method and do not escape).
- The returned value (which is externally visible, but not a side effect)

The rule of thumb here is that an externally visible side effect is anything
that would affect the execution of the remainder of the program if the function
were not executed.

!!! note
    The `:effect_free` assertion is made both for the method itself and any code
    that is executed by the method. Keep in mind that the assertion must be
    valid for all world ages and limit use of this assertion accordingly.

---
# `:nothrow`

The `:nothrow` settings asserts that this method does not terminate abnormally
(i.e. will either always return a value or never return).

!!! note
    It is permissible for `:nothrow` annotated methods to make use of exception
    handling internally as long as the exception is not rethrown out of the
    method itself.

!!! note
    `MethodErrors` and similar exceptions count as abnormal termination.

---
# `:terminates_globally`

The `:terminates_globally` settings asserts that this method will eventually terminate
(either normally or abnormally), i.e. does not loop indefinitely.

!!! note
    This `:terminates_globally` assertion covers any other methods called by the annotated method.

!!! note
    The compiler will consider this a strong indication that the method will
    terminate relatively *quickly* and may (if otherwise legal), call this
    method at compile time. I.e. it is a bad idea to annotate this setting
    on a method that *technically*, but not *practically*, terminates.

---
# `:terminates_locally`

The `:terminates_locally` setting is like `:terminates_globally`, except that it only
applies to syntactic control flow *within* the annotated method. It is thus
a much weaker (and thus safer) assertion that allows for the possibility of
non-termination if the method calls some other method that does not terminate.

!!! note
    `:terminates_globally` implies `:terminates_locally`.

---
# `:notaskstate`

The `:notaskstate` setting asserts that the method does not use or modify the
local task state (task local storage, RNG state, etc.) and may thus be safely
moved between tasks without observable results.

!!! note
    The implementation of exception handling makes use of state stored in the
    task object. However, this state is currently not considered to be within
    the scope of `:notaskstate` and is tracked separately using the `:nothrow`
    effect.

!!! note
    The `:notaskstate` assertion concerns the state of the *currently running task*.
    If a reference to a `Task` object is obtained by some other means that
    does not consider which task is *currently* running, the `:notaskstate`
    effect need not be tainted. This is true, even if said task object happens
    to be `===` to the currently running task.

!!! note
    Access to task state usually also results in the tainting of other effects,
    such as `:effect_free` (if task state is modified) or `:consistent` (if
    task state is used in the computation of the result). In particular,
    code that is not `:notaskstate`, but is `:effect_free` and `:consistent`
    may still be dead-code-eliminated and thus promoted to `:total`.

---
# `:foldable`

This setting is a convenient shortcut for the set of effects that the compiler
requires to be guaranteed to constant fold a call at compile time. It is
currently equivalent to the following `setting`s:

- `:consistent`
- `:effect_free`
- `:terminates_globally`

!!! note
    This list in particular does not include `:nothrow`. The compiler will still
    attempt constant propagation and note any thrown error at compile time. Note
    however, that by the `:consistent`-cy requirements, any such annotated call
    must consistently throw given the same argument values.

---
# `:total`

This `setting` is the maximum possible set of effects. It currently implies
the following other `setting`s:
- `:consistent`
- `:effect_free`
- `:nothrow`
- `:terminates_globally`
- `:notaskstate`

!!! warning
    `:total` is a very strong assertion and will likely gain additional semantics
    in future versions of julia (e.g. if additional effects are added and included
    in the definition of `:total`). As a result, it should be used with care.
    Whenever possible, prefer to use the minimum possible set of specific effect
    assertions required for a particular application. In cases where a large
    number of effect overrides apply to a set of functions, a custom macro is
    recommended over the use of `:total`.

---

## Negated effects

Effect names may be prefixed by `!` to indicate that the effect should be removed
from an earlier meta effect. For example, `:total !:nothrow` indicates that while
the call is generally total, it may however throw.

## Comparison to `@pure`

`@assume_effects :foldable` is similar to [`@pure`](@ref) with the primary
distinction that the `:consistent`-cy requirement applies world-age wise rather
than globally as described above. However, in particular, a method annotated
`@pure` should always be at least `:foldable`.
Another advantage is that effects introduced by `@assume_effects` are propagated to
callers interprocedurally while a purity defined by `@pure` is not.
"""
macro assume_effects(args...)
    (consistent, effect_free, nothrow, terminates_globally, terminates_locally, notaskstate) =
        (false, false, false, false, false, false, false)
    for org_setting in args[1:end-1]
        (setting, val) = compute_assumed_setting(org_setting)
        if setting === :consistent
            consistent = val
        elseif setting === :effect_free
            effect_free = val
        elseif setting === :nothrow
            nothrow = val
        elseif setting === :terminates_globally
            terminates_globally = val
        elseif setting === :terminates_locally
            terminates_locally = val
        elseif setting === :notaskstate
            notaskstate = val
        elseif setting === :foldable
            consistent = effect_free = terminates_globally = val
        elseif setting === :total
            consistent = effect_free = nothrow = terminates_globally = notaskstate = val
        else
            throw(ArgumentError("@assume_effects $org_setting not supported"))
        end
    end
    ex = args[end]
    isa(ex, Expr) || throw(ArgumentError("Bad expression `$ex` in `@assume_effects [settings] ex`"))
    if ex.head === :macrocall && ex.args[1] == Symbol("@ccall")
        ex.args[1] = GlobalRef(Base, Symbol("@ccall_effects"))
        insert!(ex.args, 3, Core.Compiler.encode_effects_override(Core.Compiler.EffectsOverride(
            consistent, effect_free, nothrow, terminates_globally, terminates_locally, notaskstate
        )))
        return esc(ex)
    end
    return esc(pushmeta!(ex, :purity, consistent, effect_free, nothrow, terminates_globally, terminates_locally, notaskstate))
end

function compute_assumed_setting(@nospecialize(setting), val::Bool=true)
    if isexpr(setting, :call) && setting.args[1] === :(!)
        return compute_assumed_setting(setting.args[2], !val)
    elseif isa(setting, QuoteNode)
        return compute_assumed_setting(setting.value, val)
    else
        return (setting, val)
    end
end

"""
    @propagate_inbounds

Tells the compiler to inline a function while retaining the caller's inbounds context.
"""
macro propagate_inbounds(ex)
    if isa(ex, Expr)
        pushmeta!(ex, :inline)
        pushmeta!(ex, :propagate_inbounds)
    end
    esc(ex)
end

"""
    @polly

Tells the compiler to apply the polyhedral optimizer Polly to a function.
"""
macro polly(ex)
    esc(isa(ex, Expr) ? pushmeta!(ex, :polly) : ex)
end

## some macro utilities ##

unwrap_macrocalls(@nospecialize(x)) = x
function unwrap_macrocalls(ex::Expr)
    inner = ex
    while inner.head === :macrocall
        inner = inner.args[end]::Expr
    end
    return inner
end

function pushmeta!(ex::Expr, sym::Symbol, args::Any...)
    if isempty(args)
        tag = sym
    else
        tag = Expr(sym, args...)::Expr
    end

    inner = unwrap_macrocalls(ex)

    idx, exargs = findmeta(inner)
    if idx != 0
        push!(exargs[idx].args, tag)
    else
        body = inner.args[2]::Expr
        pushfirst!(body.args, Expr(:meta, tag))
    end
    ex
end

popmeta!(body, sym) = _getmeta(body, sym, true)
peekmeta(body, sym) = _getmeta(body, sym, false)

function _getmeta(body::Expr, sym::Symbol, delete::Bool)
    body.head === :block || return false, []
    _getmeta(body.args, sym, delete)
end
_getmeta(arg, sym, delete::Bool) = (false, [])
function _getmeta(body::Array{Any,1}, sym::Symbol, delete::Bool)
    idx, blockargs = findmeta_block(body, args -> findmetaarg(args,sym)!=0)
    if idx == 0
        return false, []
    end
    metaargs = blockargs[idx].args
    i = findmetaarg(blockargs[idx].args, sym)
    if i == 0
        return false, []
    end
    ret = isa(metaargs[i], Expr) ? (metaargs[i]::Expr).args : []
    if delete
        deleteat!(metaargs, i)
        isempty(metaargs) && deleteat!(blockargs, idx)
    end
    true, ret
end

# Find index of `sym` in a meta expression argument list, or 0.
function findmetaarg(metaargs, sym)
    for i = 1:length(metaargs)
        arg = metaargs[i]
        if (isa(arg, Symbol) && (arg::Symbol)    == sym) ||
           (isa(arg, Expr)   && (arg::Expr).head == sym)
            return i
        end
    end
    return 0
end

function annotate_meta_def_or_block(@nospecialize(ex), meta::Symbol)
    inner = unwrap_macrocalls(ex)
    if is_function_def(inner)
        # annotation on a definition
        return esc(pushmeta!(ex, meta))
    else
        # annotation on a block
        return Expr(:block,
                    Expr(meta, true),
                    Expr(:local, Expr(:(=), :val, esc(ex))),
                    Expr(meta, false),
                    :val)
    end
end

function is_short_function_def(@nospecialize(ex))
    isexpr(ex, :(=)) || return false
    while length(ex.args) >= 1 && isa(ex.args[1], Expr)
        (ex.args[1].head === :call) && return true
        (ex.args[1].head === :where || ex.args[1].head === :(::)) || return false
        ex = ex.args[1]
    end
    return false
end
is_function_def(@nospecialize(ex)) =
    return isexpr(ex, :function) || is_short_function_def(ex) || isexpr(ex, :->)

function findmeta(ex::Expr)
    if is_function_def(ex)
        body = ex.args[2]::Expr
        body.head === :block || error(body, " is not a block expression")
        return findmeta_block(ex.args)
    end
    error(ex, " is not a function expression")
end

findmeta(ex::Array{Any,1}) = findmeta_block(ex)

function findmeta_block(exargs, argsmatch=args->true)
    for i = 1:length(exargs)
        a = exargs[i]
        if isa(a, Expr)
            if a.head === :meta && argsmatch(a.args)
                return i, exargs
            elseif a.head === :block
                idx, exa = findmeta_block(a.args, argsmatch)
                if idx != 0
                    return idx, exa
                end
            end
        end
    end
    return 0, []
end

remove_linenums!(ex) = ex
function remove_linenums!(ex::Expr)
    if ex.head === :block || ex.head === :quote
        # remove line number expressions from metadata (not argument literal or inert) position
        filter!(ex.args) do x
            isa(x, Expr) && x.head === :line && return false
            isa(x, LineNumberNode) && return false
            return true
        end
    end
    for subex in ex.args
        subex isa Expr && remove_linenums!(subex)
    end
    return ex
end
function remove_linenums!(src::CodeInfo)
    src.codelocs .= 0
    length(src.linetable) > 1 && resize!(src.linetable, 1)
    return src
end

macro generated()
    return Expr(:generated)
end

"""
    @generated f

`@generated` is used to annotate a function which will be generated.
In the body of the generated function, only types of arguments can be read
(not the values). The function returns a quoted expression evaluated when the
function is called. The `@generated` macro should not be used on functions mutating
the global scope or depending on mutable elements.

See [Metaprogramming](@ref) for further details.

## Example:
```jldoctest
julia> @generated function bar(x)
           if x <: Integer
               return :(x ^ 2)
           else
               return :(x)
           end
       end
bar (generic function with 1 method)

julia> bar(4)
16

julia> bar("baz")
"baz"
```
"""
macro generated(f)
    if isa(f, Expr) && (f.head === :function || is_short_function_def(f))
        body = f.args[2]
        lno = body.args[1]
        tmp = gensym("tmp")
        return Expr(:escape,
                    Expr(f.head, f.args[1],
                         Expr(:block,
                              lno,
                              Expr(:if, Expr(:generated),
                                   # https://github.com/JuliaLang/julia/issues/25678
                                   Expr(:block,
                                        :(local $tmp = $body),
                                        :(if $tmp isa $(GlobalRef(Core, :CodeInfo)); return $tmp; else $tmp; end)),
                                   Expr(:block,
                                        Expr(:meta, :generated_only),
                                        Expr(:return, nothing))))))
    else
        error("invalid syntax; @generated must be used with a function definition")
    end
end

"""
    @atomic var
    @atomic order ex

Mark `var` or `ex` as being performed atomically, if `ex` is a supported expression.

    @atomic a.b.x = new
    @atomic a.b.x += addend
    @atomic :release a.b.x = new
    @atomic :acquire_release a.b.x += addend

Perform the store operation expressed on the right atomically and return the
new value.

With `=`, this operation translates to a `setproperty!(a.b, :x, new)` call.
With any operator also, this operation translates to a `modifyproperty!(a.b,
:x, +, addend)[2]` call.

    @atomic a.b.x max arg2
    @atomic a.b.x + arg2
    @atomic max(a.b.x, arg2)
    @atomic :acquire_release max(a.b.x, arg2)
    @atomic :acquire_release a.b.x + arg2
    @atomic :acquire_release a.b.x max arg2

Perform the binary operation expressed on the right atomically. Store the
result into the field in the first argument and return the values `(old, new)`.

This operation translates to a `modifyproperty!(a.b, :x, func, arg2)` call.

See [Per-field atomics](@ref man-atomics) section in the manual for more details.

```jldoctest
julia> mutable struct Atomic{T}; @atomic x::T; end

julia> a = Atomic(1)
Atomic{Int64}(1)

julia> @atomic a.x # fetch field x of a, with sequential consistency
1

julia> @atomic :sequentially_consistent a.x = 2 # set field x of a, with sequential consistency
2

julia> @atomic a.x += 1 # increment field x of a, with sequential consistency
3

julia> @atomic a.x + 1 # increment field x of a, with sequential consistency
3 => 4

julia> @atomic a.x # fetch field x of a, with sequential consistency
4

julia> @atomic max(a.x, 10) # change field x of a to the max value, with sequential consistency
4 => 10

julia> @atomic a.x max 5 # again change field x of a to the max value, with sequential consistency
10 => 10
```

!!! compat "Julia 1.7"
    This functionality requires at least Julia 1.7.
"""
macro atomic(ex)
    if !isa(ex, Symbol) && !is_expr(ex, :(::))
        return make_atomic(QuoteNode(:sequentially_consistent), ex)
    end
    return esc(Expr(:atomic, ex))
end
macro atomic(order, ex)
    order isa QuoteNode || (order = esc(order))
    return make_atomic(order, ex)
end
macro atomic(a1, op, a2)
    return make_atomic(QuoteNode(:sequentially_consistent), a1, op, a2)
end
macro atomic(order, a1, op, a2)
    order isa QuoteNode || (order = esc(order))
    return make_atomic(order, a1, op, a2)
end
function make_atomic(order, ex)
    @nospecialize
    if ex isa Expr
        if isexpr(ex, :., 2)
            l, r = esc(ex.args[1]), esc(ex.args[2])
            return :(getproperty($l, $r, $order))
        elseif isexpr(ex, :call, 3)
            return make_atomic(order, ex.args[2], ex.args[1], ex.args[3])
        elseif ex.head === :(=)
            l, r = ex.args[1], esc(ex.args[2])
            if is_expr(l, :., 2)
                ll, lr = esc(l.args[1]), esc(l.args[2])
                return :(setproperty!($ll, $lr, $r, $order))
            end
        end
        if length(ex.args) == 2
            if ex.head === :(+=)
                op = :+
            elseif ex.head === :(-=)
                op = :-
            elseif @isdefined string
                shead = string(ex.head)
                if endswith(shead, '=')
                    op = Symbol(shead[1:prevind(shead, end)])
                end
            end
            if @isdefined(op)
                return Expr(:ref, make_atomic(order, ex.args[1], op, ex.args[2]), 2)
            end
        end
    end
    error("could not parse @atomic expression $ex")
end
function make_atomic(order, a1, op, a2)
    @nospecialize
    is_expr(a1, :., 2) || error("@atomic modify expression missing field access")
    a1l, a1r, op, a2 = esc(a1.args[1]), esc(a1.args[2]), esc(op), esc(a2)
    return :(modifyproperty!($a1l, $a1r, $op, $a2, $order))
end

"""
    @atomicswap a.b.x = new
    @atomicswap :sequentially_consistent a.b.x = new

Stores `new` into `a.b.x` and returns the old value of `a.b.x`.

This operation translates to a `swapproperty!(a.b, :x, new)` call.

See [Per-field atomics](@ref man-atomics) section in the manual for more details.

```jldoctest
julia> mutable struct Atomic{T}; @atomic x::T; end

julia> a = Atomic(1)
Atomic{Int64}(1)

julia> @atomicswap a.x = 2+2 # replace field x of a with 4, with sequential consistency
1

julia> @atomic a.x # fetch field x of a, with sequential consistency
4
```

!!! compat "Julia 1.7"
    This functionality requires at least Julia 1.7.
"""
macro atomicswap(order, ex)
    order isa QuoteNode || (order = esc(order))
    return make_atomicswap(order, ex)
end
macro atomicswap(ex)
    return make_atomicswap(QuoteNode(:sequentially_consistent), ex)
end
function make_atomicswap(order, ex)
    @nospecialize
    is_expr(ex, :(=), 2) || error("@atomicswap expression missing assignment")
    l, val = ex.args[1], esc(ex.args[2])
    is_expr(l, :., 2) || error("@atomicswap expression missing field access")
    ll, lr = esc(l.args[1]), esc(l.args[2])
    return :(swapproperty!($ll, $lr, $val, $order))
end

"""
    @atomicreplace a.b.x expected => desired
    @atomicreplace :sequentially_consistent a.b.x expected => desired
    @atomicreplace :sequentially_consistent :monotonic a.b.x expected => desired

Perform the conditional replacement expressed by the pair atomically, returning
the values `(old, success::Bool)`. Where `success` indicates whether the
replacement was completed.

This operation translates to a `replaceproperty!(a.b, :x, expected, desired)` call.

See [Per-field atomics](@ref man-atomics) section in the manual for more details.

```jldoctest
julia> mutable struct Atomic{T}; @atomic x::T; end

julia> a = Atomic(1)
Atomic{Int64}(1)

julia> @atomicreplace a.x 1 => 2 # replace field x of a with 2 if it was 1, with sequential consistency
(old = 1, success = true)

julia> @atomic a.x # fetch field x of a, with sequential consistency
2

julia> @atomicreplace a.x 1 => 2 # replace field x of a with 2 if it was 1, with sequential consistency
(old = 2, success = false)

julia> xchg = 2 => 0; # replace field x of a with 0 if it was 1, with sequential consistency

julia> @atomicreplace a.x xchg
(old = 2, success = true)

julia> @atomic a.x # fetch field x of a, with sequential consistency
0
```

!!! compat "Julia 1.7"
    This functionality requires at least Julia 1.7.
"""
macro atomicreplace(success_order, fail_order, ex, old_new)
    fail_order isa QuoteNode || (fail_order = esc(fail_order))
    success_order isa QuoteNode || (success_order = esc(success_order))
    return make_atomicreplace(success_order, fail_order, ex, old_new)
end
macro atomicreplace(order, ex, old_new)
    order isa QuoteNode || (order = esc(order))
    return make_atomicreplace(order, order, ex, old_new)
end
macro atomicreplace(ex, old_new)
    return make_atomicreplace(QuoteNode(:sequentially_consistent), QuoteNode(:sequentially_consistent), ex, old_new)
end
function make_atomicreplace(success_order, fail_order, ex, old_new)
    @nospecialize
    is_expr(ex, :., 2) || error("@atomicreplace expression missing field access")
    ll, lr = esc(ex.args[1]), esc(ex.args[2])
    if is_expr(old_new, :call, 3) && old_new.args[1] === :(=>)
        exp, rep = esc(old_new.args[2]), esc(old_new.args[3])
        return :(replaceproperty!($ll, $lr, $exp, $rep, $success_order, $fail_order))
    else
        old_new = esc(old_new)
        return :(replaceproperty!($ll, $lr, $old_new::Pair..., $success_order, $fail_order))
    end
end

max_stars_repo_name: cffbots/julia

max_issues_repo_name: Eurus-Holmes/julia

max_forks_repo_name: Eurus-Holmes/julia


text:
#
println("Fa) Symbolic evaluation by means of special values & recursion relations of the Wigner n-j symbols.")
#
using SymEngine

if      false
    # Use recursion relations for Wigner 3j symbols
    ja  = Basic(:ja);    jb = Basic(:jb);    jc = Basic(:jc);    ma = Basic(:ma);    mb = Basic(:mb);    mc = Basic(:mc)
    w3j = W3j(ja, jb, jc, ma, mb, mc)
    wa  = RacahAlgebra.recursionW3j(w3j, RacahAlgebra.RecursionW3jMagnetic())
    ## wa  = RacahAlgebra.recursionW3j(w3j, RacahAlgebra.RecursionW3jOneStep())
    ## wa  = RacahAlgebra.recursionW3j(w3j, RacahAlgebra.RecursionW3jHalfStep())
    ## wa  = RacahAlgebra.recursionW3j(w3j, RacahAlgebra.RecursionW3jLouck())
elseif  false
    # Look for special value for Wigner 3j symbol
    j   = Basic(:j);     m = Basic(:m)
    w3j = W3j(j+3//2, j, 3//2, m, -m-3//2, 3//2)
    rex = RacahAlgebra.evaluate(w3j);           @show rex
    rex = RacahAlgebra.equivalentForm(w3j);     @show rex
    RacahAlgebra.evaluate(rex, special=true)
elseif  false
    # Look for special value for Wigner 6j symbol
    a   = Basic(:a);     b  = Basic(:b);    c  = Basic(:c)
    w6j = W6j( a, b, c, 2, c-2, b-2)
    rex = RacahAlgebra.evaluate(w6j);           @show rex
    rex = RacahAlgebra.equivalentForm(w6j);     @show rex
    RacahAlgebra.evaluate(rex, special=true)
elseif  true
    # Look for special value for Wigner 9j symbol
    a   = Basic(:a);    b = Basic(:b);    c = Basic(:c);    d = Basic(:d);    ee = Basic(:ee);    f = Basic(:f)
    w9j = W9j(a, b, ee, c, d, ee, f, f, Basic(0))
    rex = RacahAlgebra.evaluate(w9j);           @show rex
    rex = RacahAlgebra.equivalentForm(w9j);     @show rex
    RacahAlgebra.evaluate(rex, special=true)
elseif  false
    # Evaluate the recoupling coefficient for 3 angular momenta explicitly from Racah expression; from Racah-02.ps
    j1   = Basic(:j1);    j2 = Basic(:j2);    j3 = Basic(:j3);    J12 = Basic(:J12);    J23 = Basic(:J23);    J = Basic(:J)
    m1   = Basic(:m1);    m2 = Basic(:m2);    m3 = Basic(:m3);    M12 = Basic(:M12);    M23 = Basic(:M23);    M = Basic(:M)
    w3ja = W3j(J12, j3, J, M12, m3, -M);        w3jb = W3j(j1, j2, J12, m1, m2, -M12)
    w3jc = W3j(j2, j3, J23, m2, m3, -M23);      w3jd = W3j(j1, J23, J, m1, M23, -M)
    rex  = RacahExpression( [m1, m2, m3, M12, M23], -J12 + 2*j3 - 2*M - 2*j1 - M12 - M23 + J23,
                            (2*J+1) * sqrt( (2*J12+1)*(2*J23+1) ), Kronecker[], Triangle[], [w3ja, w3jb, w3jc, w3jd], W6j[], W9j[] )
    rex  = RacahAlgebra.evaluate(rex);   @show rex
    rex  = RacahAlgebra.evaluate(rex)
elseif  false
    # Evaluate the Racah expression for the Coulomb interaction strength; from Racah-02.ps
    ja   = Basic(:ja);    jb = Basic(:jb);    jc = Basic(:jc);    jd = Basic(:jd);    Jab = Basic(:Jab);    Jcd = Basic(:Jcd);    L = Basic(:L)
    ma   = Basic(:ma);    mb = Basic(:mb);    mc = Basic(:mc);    md = Basic(:md);    Mab = Basic(:Mab);    Mcd = Basic(:Mcd);    M = Basic(:M)
    rexa = RacahAlgebra.ClebschGordan(ja, ma, jb, mb, Jab, Mab);    rexb = RacahAlgebra.ClebschGordan(jc, mc, jd, md, Jcd, Mcd)
    w3ja = W3j(ja, L, jc, -ma, M, mc);                              w3jb = W3j(jb, L, jd, -mb, -M, md)
    rex  = rexa * rexb * RacahExpression( [ma, mb, mc, md, M], L - M + ja - ma + jb - mb, Basic(1), Kronecker[], Triangle[], [w3ja, w3jb], W6j[], W9j[] )
    rex  = RacahAlgebra.evaluate(rex);   @show rex
    rex  = RacahAlgebra.evaluate(rex)
elseif  true
    # Evaluate the Racah expression for the Feynman-Goldstone diagram; from Racah-02.ps
    ja   = Basic(:ja);    jb = Basic(:jb);    jr = Basic(:jr);    js = Basic(:js);    L1 = Basic(:L1);    L2 = Basic(:L2)
    ma   = Basic(:ma);    mb = Basic(:mb);    mr = Basic(:mr);    ms = Basic(:ms);    M1 = Basic(:M1);    M2 = Basic(:M2)
    w3ja = W3j(jr, L1, ja, -mr, M1, ma);      w3jb = W3j(js, L1, jb, -ms, -M1, mb)
    w3jc = W3j(ja, L2, jr, -ma, M2, mr);      w3jd = W3j(jb, L2, js, -mb, -M2, ms)
    rex  = RacahExpression( [ma, mb, mr, ms, M1, M2], L1 - M1 + L2 - M2 + ja - ma + jb - mb + jr - mr + js - ms,
                            Basic(1), Kronecker[], Triangle[], [w3ja, w3jb, w3jc, w3jd], W6j[], W9j[] )
    rex  = RacahAlgebra.evaluate(rex);   @show rex
    rex  = RacahAlgebra.evaluate(rex)
elseif  false
    # Test for special values of Wigner 3j symbols; selectW3j(1..16)
    w3j = RacahAlgebra.selectW3j(16);                   println("w3j-original        = $w3j")
    rex = RacahAlgebra.equivalentForm(w3j);             println("rex-equivalent      = $rex")
    rex = RacahAlgebra.evaluate(rex, special=true);     println("rex-special value   = $rex")
elseif  false
    # Test for special values of the Wigner 6j symbols; selectW6j(1..22)
    w6j = RacahAlgebra.selectW6j(22);                   println("w6j-original        = $w6j")
    rex = RacahAlgebra.equivalentForm(w6j);             println("rex-equivalent      = $rex")
    rex = RacahAlgebra.evaluate(rex, special=true);     println("rex-special value   = $rex")
    ## ws  = [:ja => 2//1, :jb => 2, :jc => 1]
    ## wb  = RacahAlgebra.subs(w6j, ws);                   println("\nw6j-substitution  = $wb")
    ## wc  = RacahAlgebra.subs(wa, ws);                    println("\nrhs-substitution  = $wc")
    ## wd  = RacahAlgebra.evaluateNumerical(wb);           println("\nlhs-numerical     = $wd")
    ## we  = RacahAlgebra.evaluateNumerical(wc);           println("\nrhs-numerical     = $we")
elseif  true
    # Test for special values of the Wigner 9-j symbols; selectW6j(1..2)
    w9j = RacahAlgebra.selectW9j(2);                    println("w9j-original        = $w9j")
    rex = RacahAlgebra.equivalentForm(w9j);             println("rex-equivalent      = $rex")
    rex = RacahAlgebra.evaluate(rex, special=true);     println("rex-special value   = $rex")
end

ext: jl


text:
"""
    MLFlow

Base type which defines location and version for MLFlow API service.

# Fields
- `baseuri::String`: base MLFlow tracking URI, e.g. `http://localhost:5000`
- `apiversion`: used API version, e.g. `2.0`

# Constructors

- `MLFlow(baseuri; apiversion=2.0)`
- `MLFlow()` - defaults to `MLFlow("http://localhost:5000")`

# Examples

```@example
mlf = MLFlow()
```

"""
struct MLFlow
    baseuri::String
    apiversion
end
MLFlow(baseuri; apiversion=2.0) = MLFlow(baseuri, apiversion)
MLFlow() = MLFlow("http://localhost:5000", 2.0)
Base.show(io::IO, t::MLFlow) = show(io, ShowCase(t, new_lines=true))

"""
    MLFlowExperiment

Represents an MLFlow experiment.

# Fields
- `name::String`: experiment name.
- `lifecycle_stage::String`: life cycle stage, one of ["active", "deleted"]
- `experiment_id::Integer`: experiment identifier.
- `tags::Any`: list of tags.
- `artifact_location::String`: where are experiment artifacts stored.

# Constructors

- `MLFlowExperiment(name, lifecycle_stage, experiment_id, tags, artifact_location)`
- `MLFlowExperiment(exp::Dict{String,Any})`

"""
struct MLFlowExperiment
    name::String
    lifecycle_stage::String
    experiment_id::Integer
    tags::Any
    artifact_location::String
end
function MLFlowExperiment(exp::Dict{String,Any})
    name = get(exp, "name", missing)
    lifecycle_stage = get(exp, "lifecycle_stage", missing)
    experiment_id = parse(Int, get(exp, "experiment_id", missing))
    tags = get(exp, "tags", missing)
    artifact_location = get(exp, "artifact_location", missing)
    MLFlowExperiment(name, lifecycle_stage, experiment_id, tags, artifact_location)
end
Base.show(io::IO, t::MLFlowExperiment) = show(io, ShowCase(t, new_lines=true))

"""
    MLFlowRunStatus

Represents the status of an MLFlow Run.

# Fields
- `status::String`: one of RUNNING/SCHEDULED/FINISHED/FAILED/KILLED

# Constructors

- `MLFlowRunStatus(status::String)`
"""
struct MLFlowRunStatus
    status::String
    function MLFlowRunStatus(status::String)
        acceptable_statuses = ["RUNNING", "SCHEDULED", "FINISHED", "FAILED", "KILLED"]
        status ∈ acceptable_statuses || error("Invalid status $status - choose one of $acceptable_statuses")
        new(status)
    end
end
Base.show(io::IO, t::MLFlowRunStatus) = show(io, ShowCase(t, new_lines=true))

"""
    MLFlowRunInfo

Represents run metadata.

# Fields
- `run_id::String`: run identifier.
- `experiment_id::Integer`: experiment identifier.
- `status::MLFlowRunStatus`: run status.
- `start_time::Union{Int64,Missing}`: when was the run started, UNIX time in milliseconds.
- `end_time::Union{Int64,Missing}`: when did the run end, UNIX time in milliseconds.
- `artifact_uri::String`: where are artifacts from this run stored.
- `lifecycle_stage::String`: one of `active` or `deleted`.

# Constructors

- `MLFlowRunInfo(run_id, experiment_id, status, start_time, end_time, artifact_uri, lifecycle_stage)`
- `MLFlowRunInfo(info::Dict{String,Any})`
"""
struct MLFlowRunInfo
    run_id::String
    experiment_id::Integer
    status::MLFlowRunStatus
    start_time::Union{Int64,Missing}
    end_time::Union{Int64,Missing}
    artifact_uri::String
    lifecycle_stage::String
end
function MLFlowRunInfo(info::Dict{String,Any})
    run_id = get(info, "run_id", missing)
    experiment_id = get(info, "experiment_id", missing)
    status = get(info, "status", missing)
    start_time = get(info, "start_time", missing)
    end_time = get(info, "end_time", missing)
    artifact_uri = get(info, "artifact_uri", "")
    lifecycle_stage = get(info, "lifecycle_stage", "")

    experiment_id = ismissing(experiment_id) ? experiment_id : parse(Int64, experiment_id)
    status = ismissing(status) ? status : MLFlowRunStatus(status)

    # support for mlflow 1.21.0
    if !ismissing(start_time) && !(typeof(start_time) <: Int)
        start_time = parse(Int64, start_time)
    end
    if !ismissing(end_time) && !(typeof(end_time) <: Int)
        end_time = parse(Int64, end_time)
    end
    MLFlowRunInfo(run_id, experiment_id, status, start_time, end_time, artifact_uri, lifecycle_stage)
end
Base.show(io::IO, t::MLFlowRunInfo) = show(io, ShowCase(t, new_lines=true))
get_run_id(runinfo::MLFlowRunInfo) = runinfo.run_id

"""
    MLFlowRunDataMetric

Represents a metric.

# Fields
- `key::String`: metric identifier.
- `value::Float64`: metric value.
- `step::Int64`: step.
- `timestamp::Int64`: timestamp in UNIX time in milliseconds.

# Constructors

- `MLFlowRunDataMetric(d::Dict{String,Any})`

"""
struct MLFlowRunDataMetric
    key::String
    value::Float64
    step::Int64
    timestamp::Int64
end
function MLFlowRunDataMetric(d::Dict{String,Any})
    key = d["key"]
    value = d["value"]
    if typeof(d["step"]) <: Int
        step = d["step"]
    else
        step = parse(Int64, d["step"])
    end
    if typeof(d["timestamp"]) <: Int
        timestamp = d["timestamp"]
    else
        timestamp = parse(Int64, d["timestamp"])
    end
    MLFlowRunDataMetric(key, value, step, timestamp)
end
Base.show(io::IO, t::MLFlowRunDataMetric) = show(io, ShowCase(t, new_lines=true))

"""
    MLFlowRunData

Represents run data.

# Fields
- `metrics::Dict{String,MLFlowRunDataMetric}`: run metrics.
- `params::Dict{String,String}`: run parameters.
- `tags`: list of run tags.

# Constructors

- `MLFlowRunData(data::Dict{String,Any})`

"""
struct MLFlowRunData
    metrics::Dict{String,MLFlowRunDataMetric}
    params::Union{Dict{String,String},Missing}
    tags
end
function MLFlowRunData(data::Dict{String,Any})
    metrics = Dict{String,MLFlowRunDataMetric}()
    if haskey(data, "metrics")
        for metric in data["metrics"]
            v = MLFlowRunDataMetric(metric)
            metrics[v.key] = v
        end
    end
    if haskey(data, "params")
        params = Dict{String,String}()
        for p in data["params"]
            params[p["key"]] = p["value"]
        end
    else
        params = Dict{String,String}()
    end
    tags = haskey(data, "tags") ? data["tags"] : missing
    MLFlowRunData(metrics, params, tags)
end
Base.show(io::IO, t::MLFlowRunData) = show(io, ShowCase(t, new_lines=true))
get_params(rundata::MLFlowRunData) = rundata.params

"""
    MLFlowRun

Represents an MLFlow run.

# Fields
- `info::MLFlowRunInfo`: Run metadata.
- `data::MLFlowRunData`: Run data.

# Constructors

- `MLFlowRun(rundata::MLFlowRunData)`
- `MLFlowRun(runinfo::MLFlowRunInfo)`
- `MLFlowRun(info::Dict{String,Any})`
- `MLFlowRun(info::Dict{String,Any}, data::Dict{String,Any})`

"""
struct MLFlowRun
    info::Union{MLFlowRunInfo,Missing}
    data::Union{MLFlowRunData,Missing}
end
MLFlowRun(rundata::MLFlowRunData) =
    MLFlowRun(missing, rundata)
MLFlowRun(runinfo::MLFlowRunInfo) =
    MLFlowRun(runinfo, missing)
MLFlowRun(info::Dict{String,Any}) =
    MLFlowRun(MLFlowRunInfo(info), missing)
MLFlowRun(info::Dict{String,Any}, data::Dict{String,Any}) =
    MLFlowRun(MLFlowRunInfo(info), MLFlowRunData(data))
Base.show(io::IO, t::MLFlowRun) = show(io, ShowCase(t, new_lines=true))
get_info(run::MLFlowRun) = run.info
get_data(run::MLFlowRun) = run.data
get_run_id(run::MLFlowRun) = get_run_id(run.info)
get_params(run::MLFlowRun) = get_params(run.data)

"""
    MLFlowArtifactFileInfo

Metadata of a single artifact file -- result of [`listartifacts`](@ref).

# Fields
- `filepath::String`: File path, including the root artifact directory of a run.
- `filesize::Int64`: Size in bytes.
"""
struct MLFlowArtifactFileInfo
    filepath::String
    filesize::Int64
end
Base.show(io::IO, t::MLFlowArtifactFileInfo) = show(io, ShowCase(t, new_lines=true))
get_path(mlfafi::MLFlowArtifactFileInfo) = mlfafi.filepath
get_size(mlfafi::MLFlowArtifactFileInfo) = mlfafi.filesize

"""
    MLFlowArtifactDirInfo

Metadata of a single artifact directory -- result of [`listartifacts`](@ref).

# Fields
- `dirpath::String`: Directory path, including the root artifact directory of a run.
"""
struct MLFlowArtifactDirInfo
    dirpath::String
end
Base.show(io::IO, t::MLFlowArtifactDirInfo) = show(io, ShowCase(t, new_lines=true))
get_path(mlfadi::MLFlowArtifactDirInfo) = mlfadi.dirpath
get_size(mlfadi::MLFlowArtifactDirInfo) = 0

max_stars_repo_name: JuliaAI/MLFlowClient

max_issues_repo_name: JuliaAI/MLFlowClient

max_forks_repo_name: JuliaAI/MLFlowClient


text:
using CoordinateSplittingPTrees
import AbstractTrees
using Test, Random, LinearAlgebra, SparseArrays
using GLPK
include("functions.jl")

struct DummyMeta
    val::Float64
end

@testset "World" begin
    world = @inferred(World([0], [1], [0], 10))
    @test eltype(world.lower) == Float64
    world = @inferred(World([0.0], [1], [0], 10))
    @test eltype(world.lower) == Float64
    world = @inferred(World(zeros(5), ones(5), rand(5), 10))
    @test eltype(world.lower) == Float64
    @test_throws ArgumentError World(zeros(2), ones(2), [7.0,0.5], 0)
    world = @inferred(World(zeros(2), [Inf,1], [7.0,0.5], 0))
    @test CoordinateSplittingPTrees.baseposition(world) == [7.0, 0.5]
    world = @inferred(World([1,2], 0))
    @test world.lower == [-Inf,-Inf]
    @test world.upper == [Inf,Inf]
    @test CoordinateSplittingPTrees.baseposition(world) == [1.0,2.0]
    world = @inferred(World([1,2], x->3.0))
    @test world.meta == 3.0

    # Behavior reserved for use by QuadSplit (2-tuple splits)
    world = @inferred(World(zeros(3), fill(3, 3), fill((1,2), 3), DummyMeta(0)))
    @test eltype(world.lower) == Float64
    @test CoordinateSplittingPTrees.baseposition(world) == ones(3)
    world = World(zeros(3), fill(3, 3), fill((2,1), 3), DummyMeta(0))
    @test CoordinateSplittingPTrees.baseposition(world) == fill(2, 3)
    @test_throws ArgumentError World(zeros(3), fill(3, 3), fill((1,1), 3), DummyMeta(0))
    world = World(zeros(3), fill(3, 3), fill((2,1), 3), DummyMeta(0))
    @test CoordinateSplittingPTrees.baseposition(world) == fill(2, 3)
    world = World(fill((2,1), 3), DummyMeta(0))
    @test world.lower == [-Inf,-Inf,-Inf]
    @test world.upper == [Inf,Inf,Inf]
    @test CoordinateSplittingPTrees.baseposition(world) == fill(2, 3)

    world = World(fill(-Inf, 3), fill(Inf, 3), fill((1,2.5), 3), 0)
    root = Box{2}(world)
    @test position(root) == ones(3)
    box = CoordinateSplittingPTrees.addpoint_distinct!(root, [5,1.0,3], x->0)
    @test position(box) == [5,2.5,3]
end

@testset "Geometry and iteration, CS1" begin
    # For comparing boxbounds
    myapproxeq(t1::Tuple{Real,Real}, t2::Tuple{Real,Real}) = t1[1] ≈ t2[1] && t1[2] ≈ t2[2]
    myapproxeq(v1::Vector, v2::Vector) = all(myapproxeq.(v1, v2))

    # 1-d
    world = World([0], [1], [0], 10)
    @test eltype(world.lower) == Float64
    root = @inferred(Box{1}(world))
    @test_throws ErrorException Box(root, 1, -0.5, 20)  # out-of-bounds
    @test_throws ErrorException Box(root, 1,  1.5, 20)  # out-of-bounds
    @test_throws ErrorException Box(root, 1,    0, 20)  # coincides with parent
    box = @inferred(Box(root, 1, 1, 20))[1]
    @test @inferred(meta(box)) == 20
    @test @inferred(position(box)) == [1]
    @test @inferred(position(box, 1)) == 1
    @test_throws BoundsError position(box, 2)
    @test @inferred(boxbounds(box)) == [(0.5,1)]
    @test @inferred(boxbounds(box, 1)) == (0.5,1)
    @test_throws BoundsError boxbounds(box, 2)
    @test CoordinateSplittingPTrees.epswidth(boxbounds(box, 1)) == eps()
    @test CoordinateSplittingPTrees.epswidth(boxbounds(getleaf(root), 1)) == eps(0.5)
    # evaluation point is at the boundary of parents
    @test_throws AssertionError Box(root, 1, 0.5, 30)  # can only split leaves
    @test_throws ErrorException Box(getleaf(root), 1, 0.5, 30)  # can't eval at upper edge
    boxes = collect(root)
    @test boxes == [root, root.split.self, box]
    boxes = collect(box)
    @test boxes == [box]
    boxes = collect(leaves(root))
    @test boxes == [root.split.self, box]
    box1 = Box(box, 1, 0.5, 30)[1]
    @test boxbounds(box1, 1) == (0.5, 0.75)
    @test boxbounds(getleaf(box), 1) == (0.75, 1.0)
    @test boxbounds(getleaf(root), 1) == (0.0, 0.5)
    @test boxbounds(box1) == [(0.5, 0.75)]
    @test boxbounds(getleaf(box)) == [(0.75, 1.0)]
    @test boxbounds(getleaf(root)) == [(0.0, 0.5)]
    leaf = find_leaf_at(root, [0.3])
    @test isleaf(leaf)
    @test leaf == getleaf(root)
    leaf = find_leaf_at(root, [0.5])
    @test isleaf(leaf)
    @test leaf === box1
    leaf = find_leaf_at(root, [0.8])
    @test isleaf(leaf)
    @test leaf === getleaf(box)
    @test_throws ErrorException find_leaf_at(root, [1.2])
    @test length(root) == 5
    @test length(leaves(root)) == 3
    @test position(leaf) ∈ root
    @test position(leaf) ∈ leaf

    world = World([0], [Inf], [0], 10) # with infinite size
    root = Box{1}(world)
    box = Box(root, 1, 1, 20)[1]
    @test boxbounds(box) == [(0.5,Inf)]
    @test boxbounds(box, 1) == (0.5,Inf)
    @test CoordinateSplittingPTrees.epswidth(boxbounds(box, 1)) == eps(0.5)
    @test CoordinateSplittingPTrees.epswidth(boxbounds(getleaf(root), 1)) == eps(0.5)

    world = World([0], [1], [1], 10)
    root = Box{1}(world)
    box1 = Box(root, 1, 0, 20)[1]
    leaf = find_leaf_at(root, 0.5)
    @test leaf == getleaf(root)

    geom = Dict()
    root = generate_randboxes(Box{1}, 1, 10, (args...)->record_geometry!(geom, args...))
    for leaf in leaves(root)
        pos, bbs = geom[leaf]
        @test position(leaf) == pos
        @test boxbounds(leaf) == bbs
    end

    # 2-d
    world = World([0,-Inf], [Inf,Inf], [1,1], nothing)
    root = @inferred(Box{1}(world))
    box1 = @inferred(Box(root, 1, 2, nothing))[1]
    box2 = Box(getleaf(root), 2, 2, nothing)[1]

    b = getleaf(root)
    @test @inferred(position(b, 1)) == position(b, 2) == 1
    @test @inferred(position(b)) == [1,1]
    @test @inferred(boxbounds(b, 1)) == (0,1.5)
    @test boxbounds(b, 2) == (-Inf,1.5)
    @test @inferred(boxbounds(b)) == [(0,1.5), (-Inf,1.5)]

    b = box2
    @test position(b, 1) == 1
    @test position(b, 2) == 2
    @test position(b) == [1,2]
    @test boxbounds(b, 1) == (0,1.5)
    @test boxbounds(b, 2) == (1.5,Inf)
    @test boxbounds(b) == [(0,1.5), (1.5,Inf)]

    b = box1
    @test position(b, 1) == 2
    @test position(b, 2) == 1
    @test position(b) == [2,1]
    @test boxbounds(b, 1) == (1.5,Inf)
    @test boxbounds(b, 2) == (-Inf,Inf)  # never split along 2
    @test boxbounds(b) == [(1.5,Inf), (-Inf,Inf)]

    b = getleaf(root)
    @test b.parent.split.dims == (2,)
    b2 = Box(b, 1, 1/3, nothing)[1]
    @test position(b2, 1) ≈ 1/3
    @test position(b2, 2) == 1
    @test position(b2) ≈ [1/3, 1]
    @test myapproxeq(boxbounds(b2, 1), (0, 2/3))
    @test boxbounds(b2, 2) == (-Inf,1.5)
    @test myapproxeq(boxbounds(b2), [(0, 2/3), (-Inf,1.5)])
    b1 = b.split.self
    @test position(b1, 1) == 1
    @test position(b1, 2) == 1
    @test position(b1) == [1, 1]
    @test myapproxeq(boxbounds(b1, 1), (2/3, 1.5))
    @test boxbounds(b1, 2) == (-Inf,1.5)
    @test myapproxeq(boxbounds(b1), [(2/3, 1.5), (-Inf,1.5)])
    b2_2 = Box(b2, 2, -1, nothing)[1]
    @test position(b2_2, 1) ≈ 1/3
    @test position(b2_2, 2) ≈ -1
    @test position(b2_2) ≈ [1/3, -1]
    @test myapproxeq(boxbounds(b2_2, 1), (0, 2/3))
    @test boxbounds(b2_2, 2) == (-Inf,0)
    @test myapproxeq(boxbounds(b2_2), [(0, 2/3), (-Inf,0)])
    b2_1 = b2.split.self
    @test position(b2_1, 1) == position(b2, 1)
    @test position(b2_1, 2) == position(b2, 2)
    @test position(b2_1) == position(b2)
    @test myapproxeq(boxbounds(b2_1, 1), (0, 2/3))
    @test boxbounds(b2_1, 2) == (0,1.5)
    @test myapproxeq(boxbounds(b2_1), [(0, 2/3), (0,1.5)])
    b1_2 = Box(b1, 2, -1, nothing)[1]
    @test position(b1_2, 1) == 1
    @test position(b1_2, 2) == -1
    @test position(b1_2) == [1, -1]
    @test myapproxeq(boxbounds(b1_2, 1), (2/3, 1.5))
    @test boxbounds(b1_2, 2) == (-Inf, 0)
    @test myapproxeq(boxbounds(b1_2), [(2/3, 1.5), (-Inf,0)])
    b1_1 = b1.split.self
    @test position(b1_1, 1) == position(b1, 1)
    @test position(b1_1, 2) == position(b1, 2)
    @test position(b1_1) == position(b1)
    @test myapproxeq(boxbounds(b1_1, 1), (2/3, 1.5))
    @test boxbounds(b1_1, 2) == (0, 1.5)
    @test myapproxeq(boxbounds(b1_1), [(2/3, 1.5), (0, 1.5)])

    @test all(CoordinateSplittingPTrees.splittable, boxbounds(b1_1))

    geom = Dict()
    root = generate_randboxes(Box{1}, 2, 10, (args...)->record_geometry!(geom, args...))
    for leaf in leaves(root)
        pos, bbs = geom[leaf]
        @test position(leaf) == pos
        @test boxbounds(leaf) == bbs
    end

    # 5-d
    geom = Dict()
    root = generate_randboxes(Box{1}, 5, 20, (args...)->record_geometry!(geom, args...))
    for leaf in leaves(root)
        pos, bbs = geom[leaf]
        @test position(leaf) == pos
        @test boxbounds(leaf) == bbs
    end
end

@testset "Geometry and iteration, CS2" begin
    myapproxeq(t1::Tuple{Real,Real}, t2::Tuple{Real,Real}) = t1[1] ≈ t2[1] && t1[2] ≈ t2[2]
    myapproxeq(v1::Vector, v2::Vector) = all(myapproxeq.(v1, v2))

    # 1-d
    world = World([-Inf], [Inf], [1], 0)
    root = Box{2}(world)
    addpoint!(root, [1.2], x->0)
    # test for displacement along fictive dimension
    chldrn = collect(AbstractTrees.children(root))
    @test !CoordinateSplittingPTrees.isfake(chldrn[1])
    @test !CoordinateSplittingPTrees.isfake(chldrn[2])
    @test  CoordinateSplittingPTrees.isfake(chldrn[3])
    @test  CoordinateSplittingPTrees.isfake(chldrn[4])
    @test length(leaves(root)) == 2
    box1 = addpoint!(root, [0.75], x->0)
    @test length(leaves(root)) == 3
    box2 = addpoint!(root, [0.5], x->0)
    @test length(leaves(root)) == 4
    @test boxbounds(box2, 1) == (-Inf, 5/8)
    @test boxbounds(box1, 1) == (-Inf, 7/8)
    @test boxbounds(getleaf(box1), 1) == (5/8, 7/8)

    world = World([0], [1], [0], 10)
    root = Box{1}(world)
    box1 = Box(root, 1, 1, 20)[1]
    leaf = find_leaf_at(root, 0.5)
    @test leaf == box1
    world = World([0], [1], [1], 10)
    root = Box{1}(world)
    box1 = Box(root, 1, 0, 20)[1]
    leaf = find_leaf_at(root, 0.5)
    @test leaf == getleaf(root)

    # 2-d
    world = World([0,-Inf], [Inf,Inf], [1,1], nothing)
    root = @inferred(Box{2}(world))
    boxes1 = @inferred(Box(root, (1,2), (2,2), (nothing, nothing, nothing)))

    b = getleaf(root)
    @test @inferred(position(b, 1)) == position(b, 2) == 1
    @test @inferred(position(b)) == [1,1]
    @test @inferred(boxbounds(b, 1)) == (0,1.5)
    @test boxbounds(b, 2) == (-Inf,1.5)
    @test @inferred(boxbounds(b)) == [(0,1.5), (-Inf,1.5)]

    b = boxes1[2]
    @test position(b, 1) == 1
    @test position(b, 2) == 2
    @test position(b) == [1,2]
    @test boxbounds(b, 1) == (0,1.5)
    @test boxbounds(b, 2) == (1.5,Inf)
    @test boxbounds(b) == [(0,1.5), (1.5,Inf)]

    b = boxes1[1]
    @test position(b, 1) == 2
    @test position(b, 2) == 1
    @test position(b) == [2,1]
    @test boxbounds(b, 1) == (1.5,Inf)
    @test boxbounds(b, 2) == (-Inf,1.5)
    @test boxbounds(b) == [(1.5,Inf), (-Inf,1.5)]

    # 3-d (fake dimensions and iteration)
    world = World(fill(-Inf,3), fill(Inf,3), fill(0,3), 1)
    root = Box{2}(world)
    @test length(root) == 1
    @test length(leaves(root)) == 1
    boxes1 = Box(root, (1,2), (1,1), (2, 3, 4)) # split along dims 1 & 2
    @test length(root) == 5
    @test length(leaves(root)) == 4
    boxes2 = Box(boxes1[end], (3,4), (1,1), (5, 6, 7)) # split along dims 1 & 2
    @test length(leaves(root)) == 5
    v = [value(leaf) for leaf in leaves(root)]
    @test v == [1:5;]

    # 4-d
    world = World([0,-Inf,-5,-Inf], [Inf,Inf,50,20], fill(1, 4), 1)
    root = Box{2}(world)
    @test length(root) == 1
    @test length(leaves(root)) == 1
    boxes1 = Box(root, (1,2), (2,2), (2, 3, 4)) # split along dims 1 & 2
    @test length(root) == 5
    @test length(leaves(root)) == 4

    b = getleaf(root)
    @test position(b, 1) == position(b, 2) == position(b, 3) == position(b, 4) == 1
    @test position(b) == [1,1,1,1]
    @test boxbounds(b, 1) == (0,1.5)
    @test boxbounds(b, 2) == (-Inf,1.5)
    @test boxbounds(b, 3) == (-5,50)
    @test boxbounds(b, 4) == (-Inf,20)
    @test boxbounds(b) == [(0,1.5), (-Inf,1.5), (-5,50), (-Inf,20)]
    @test meta(b) == 1

    b = boxes1[1]
    @test position(b, 1) == 2
    @test position(b, 2) == position(b, 3) == position(b, 4) == 1
    @test position(b) == [2,1,1,1]
    @test boxbounds(b, 1) == (1.5,Inf)
    @test boxbounds(b, 2) == (-Inf,1.5)
    @test boxbounds(b, 3) == (-5,50)
    @test boxbounds(b, 4) == (-Inf,20)
    @test boxbounds(b) == [(1.5,Inf), (-Inf,1.5), (-5,50), (-Inf,20)]
    @test meta(b) == 2

    b = boxes1[2]
    @test position(b, 1) == position(b, 3) == position(b, 4) == 1
    @test position(b, 2) == 2
    @test position(b) == [1,2,1,1]
    @test boxbounds(b, 1) == (0,1.5)
    @test boxbounds(b, 2) == (1.5,Inf)
    @test boxbounds(b, 3) == (-5,50)
    @test boxbounds(b, 4) == (-Inf,20)
    @test boxbounds(b) == [(0,1.5), (1.5,Inf), (-5,50), (-Inf,20)]
    @test meta(b) == 3

    b = boxes1[3]
    @test position(b, 1) == 2
    @test position(b, 2) == 2
    @test position(b) == [2,2,1,1]
    @test boxbounds(b, 1) == (1.5,Inf)
    @test boxbounds(b, 2) == (1.5,Inf)
    @test boxbounds(b, 3) == (-5,50)
    @test boxbounds(b, 4) == (-Inf,20)
    @test boxbounds(b) == [(1.5,Inf), (1.5,Inf), (-5,50), (-Inf,20)]
    @test meta(b) == 4

    boxes2 = Box(boxes1[2], (3,4), (2,3), (5, 6, 7)) # split box@[1,2,1,1] along dims 3 & 4
    @test length(root) == 9
    @test length(leaves(root)) == 7
    b = getleaf(boxes1[2])
    @test position(b) == [1,2,1,1]
    @test boxbounds(b, 3) == (-5,1.5)
    @test boxbounds(b, 4) == (-Inf,2)
    @test boxbounds(b) == [(0,1.5), (1.5,Inf), (-5,1.5), (-Inf,2)]
    @test meta(b) == 3

    b = boxes2[1]
    @test position(b, 2) == position(b, 3) == 2
    @test position(b, 1) == position(b, 4) == 1
    @test position(b) == [1,2,2,1]
    @test boxbounds(b, 1) == (0,1.5)
    @test boxbounds(b, 2) == (1.5,Inf)
    @test boxbounds(b, 3) == (1.5,50)
    @test boxbounds(b, 4) == (-Inf,2)
    @test boxbounds(b) == [(0,1.5), (1.5,Inf), (1.5,50), (-Inf,2)]
    @test meta(b) == 5

    b = boxes2[2]
    @test position(b, 1) == position(b, 3) == 1
    @test position(b, 2) == 2
    @test position(b, 4) == 3
    @test position(b) == [1,2,1,3]
    @test boxbounds(b, 1) == (0,1.5)
    @test boxbounds(b, 2) == (1.5,Inf)
    @test boxbounds(b, 3) == (-5,1.5)
    @test boxbounds(b, 4) == (2,20)
    @test boxbounds(b) == [(0,1.5), (1.5,Inf), (-5,1.5), (2,20)]
    @test meta(b) == 6

    b = boxes2[3]
    @test position(b, 2) == position(b, 3) == 2
    @test position(b, 1) == 1
    @test position(b, 4) == 3
    @test position(b) == [1,2,2,3]
    @test boxbounds(b, 1) == (0,1.5)
    @test boxbounds(b, 2) == (1.5,Inf)
    @test boxbounds(b, 3) == (1.5,50)
    @test boxbounds(b, 4) == (2,20)
    @test boxbounds(b) == [(0,1.5), (1.5,Inf), (1.5,50), (2,20)]
    @test meta(b) == 7

    boxes3 = Box(boxes2[3], (3,1), (1.8,1.3), (8, 9, 10)) # split box@[1,2,2,3] along dims 3 & 1
    @test length(root) == 13
    @test length(leaves(root)) == 10
    b = getleaf(boxes2[3])
    @test position(b) == [1,2,2,3]
    @test boxbounds(b, 1) == (0,1.15)
    @test boxbounds(b, 3) == (1.9,50)
    @test boxbounds(b) == [(0,1.15), (1.5,Inf), (1.9,50), (2,20)]
    @test meta(b) == 7

    b = boxes3[1]
    @test position(b, 1) == 1
    @test position(b, 2) == 2
    @test position(b, 3) == 1.8
    @test position(b, 4) == 3
    @test position(b) == [1,2,1.8,3]
    @test boxbounds(b, 1) == (0,1.15)
    @test boxbounds(b, 2) == (1.5,Inf)
    @test boxbounds(b, 3) == (1.5,1.9)
    @test boxbounds(b, 4) == (2,20)
    @test boxbounds(b) == [(0,1.15), (1.5,Inf), (1.5,1.9), (2,20)]
    @test meta(b) == 8

    b = boxes3[2]
    @test position(b, 1) == 1.3
    @test position(b, 2) == position(b, 3) == 2
    @test position(b, 4) == 3
    @test position(b) == [1.3,2,2,3]
    @test boxbounds(b, 1) == (1.15,1.5)
    @test boxbounds(b, 2) == (1.5,Inf)
    @test boxbounds(b, 3) == (1.9,50)
    @test boxbounds(b, 4) == (2,20)
    @test boxbounds(b) == [(1.15,1.5), (1.5,Inf), (1.9,50), (2,20)]
    @test meta(b) == 9

    b = boxes3[3]
    @test position(b, 1) == 1.3
    @test position(b, 2) == 2
    @test position(b, 3) == 1.8
    @test position(b, 4) == 3
    @test position(b) == [1.3,2,1.8,3]
    @test boxbounds(b, 1) == (1.15,1.5)
    @test boxbounds(b, 2) == (1.5,Inf)
    @test boxbounds(b, 3) == (1.5,1.9)
    @test boxbounds(b, 4) == (2,20)
    @test boxbounds(b) == [(1.15,1.5), (1.5,Inf), (1.5,1.9), (2,20)]
    @test meta(b) == 10

    # 5-d
    geom = Dict()
    root = generate_randboxes(Box{2}, 5, 100, (args...)->record_geometry!(geom, args...))
    for leaf in leaves(root)
        pos, bbs = geom[leaf]
        @test position(leaf) == pos
        @test boxbounds(leaf) == bbs
    end
end

@testset "Split iteration" begin
    f(x) = rand()

    # Odd dimensionality for CSp with p even (fictive dimensions)
    n = 3
    world = World(fill(-Inf, n), fill(Inf, n), fill(0, n), rand())
    root = Box{2}(world)
    x = ones(n)
    box = addpoint!(root, x, [(1,2), (3,)], f)
    @test position(box) == x
    x = [0.6, 2, 2]
    box = addpoint!(root, x, [(2,3), (1,)], f)
    @test position(box) == [0.6,2,2]
    @test boxbounds(box) == [(0.5, 0.8), (1.5,Inf), (1.5,Inf)]

    s = [split.dims for split in splits(box)]
    @test s == [(1,4), (2,3), (3,4), (1,2)]
    box = box.parent
    s = [split.dims for split in splits(box)]
    @test s == [(1,4), (2,3), (3,4), (1,2)]
    box = box.parent
    s = [split.dims for split in splits(box)]
    @test s == [(2,3), (1,4), (3,4), (1,2)]
    box = box.parent
    s = [split.dims for split in splits(box)]
    @test s == [(3,4), (2,3), (1,4), (1,2)]
    box = box.parent
    s = [split.dims for split in splits(box)]
    @test s == [(1,2), (3,4), (2,3), (1,4)]

    n = 6
    world = World(fill(-Inf, n), fill(Inf, n), fill(0, n), rand())
    root = Box{2}(world)
    x = ones(n); box = addpoint!(root, x, [(1,2), (3,4), (5,6)], f)
    position(box) == x
    x = [0.8; 0.8; 0.8; 0.2; 0.2; 0.2]; addpoint!(root, x, [(1,3), (2,5), (4,6)], f)
    x = fill(0.2, n); addpoint!(root, x, [(1,5), (3,6), (2,4)], f)
    box = find_leaf_at(root, [0.2,0,0.2,0,0.2,0])
    s = [split.dims for split in splits(box)]
    @test s == [(2,4), (3,6), (1,5), (3,4), (1,3), (2,5), (4,6), (5,6), (1,2)]
    box = box.parent
    s = [split.dims for split in splits(box)]
    @test s == [(3,6), (2,4), (1,5), (3,4), (1,3), (2,5), (4,6), (5,6), (1,2)]
    box = find_leaf_at(root, [0.8,0.8,0.8,0.2,0.2,0.2])
    s = [split.dims for split in splits(box)]
    @test s == [(4,6), (2,5), (1,3), (5,6), (3,4), (1,5), (3,6), (2,4), (1,2)]
    box = find_leaf_at(root, ones(6))
    s = [split.dims for split in splits(box)]
    @test s == [(5,6), (1,3), (2,5), (4,6), (3,4), (1,5), (3,6), (2,4), (1,2)]
    box = find_leaf_at(root, [1,1,0,0,0,0])
    s = [split.dims for split in splits(box)]
    @test s == [(1,3), (2,5), (4,6), (5,6), (3,4), (1,5), (3,6), (2,4), (1,2)]
    s = [split.dims for split in splits(root)]
    @test s == [(1,2), (1,5), (3,6), (2,4), (3,4), (1,3), (2,5), (4,6), (5,6)]

    n = 8
    world = World(fill(-Inf, n), fill(Inf, n), fill(0, n), rand())
    for i = 1:20
        root = Box{2}(world)
        x = randn(n); addpoint!(root, x, [(1,2), (3,4), (5,6), (7,8)], f)
        x = randn(n); addpoint!(root, x, [(3,7), (4,8), (1,5), (2,6)], f)
        x = randn(n); addpoint!(root, x, [(6,8), (1,3), (2,4), (5,7)], f)
        x = randn(n); addpoint!(root, x, [(4,5), (1,6), (2,7), (3,8)], f)
        box = minimum(root)
        s = [split.dims for split in splits(box)]
        @test length(s) == 16 && length(unique(s)) == 16
    end
end

@testset "Collinearity" begin
    for i = 1:20
        root = generate_randboxes(Box{1}, 5, 20)
        X = collect_positions(root)
        for leaf in leaves(root)
            @test CoordinateSplittingPTrees.ncollinear(leaf) == collinear_dims(X, position(leaf))
        end
    end

    for i = 1:20
        root = generate_randboxes(Box{2}, 5, 20)
        X = collect_positions(root)
        for leaf in leaves(root)
            @test CoordinateSplittingPTrees.ncollinear(leaf) == collinear_dims(X, position(leaf))
        end
    end
end

@testset "Chains" begin
    # Cases with and without fictive dimensions
    f(x) = rand()
    for (n, dimlistss) in ((3, ([(1,2), (3,)],  [(3,1), (2,)])),
                           (4, ([(1,2), (3,4)], [(3,1), (2,4)])),
                           (4, ([(1,2), (3,4)], [(1,2), (3,4)])))
        x0 = randn(n)
        world = World(fill(-Inf, n), fill(Inf, n), x0, f(x0))
        root = box = Box{2}(world)
        chaintops = typeof(root)[]
        repeats = dimlistss[2] == dimlistss[1]
        for dimlists in dimlistss
            x = position(box) .+ 1
            bx = box
            box = addpoint!(box, x, dimlists, f)
            for i = 1:ceil(Int, n/2)
                push!(chaintops, bx)
                if repeats
                    # If the dimlists repeat, then any split box is
                    # top of some chain
                    bx = bx.split.others.children[end]
                end
            end
        end
        if repeats
            chaintops[end] = chaintops[end].parent # need full chain
        end
        i = length(chaintops)
        while !isroot(box)
            top, success = chaintop(box)
            @test success
            @test top == chaintops[i]
            ntests = 0
            for (j, split) in enumerate(chain(top))
                @test j < 3
                if j == 1
                    @test split == top.split
                    ntests += 1
                elseif j == 2
                    @test split == top.split.others.children[end].split
                    ntests += 1
                end
            end
            @test ntests == 2
            for child in box.parent.split.others.children
                top, success = CoordinateSplittingPTrees.chaintop(child)
                @test success
                @test top == chaintops[i]
            end
            i -= 1
            box = box.parent
        end
    end
end

@testset "p-th order direct coefficients" begin
    bprod(x, B::AbstractVector) = B'*x
    bprod(x, B::AbstractMatrix) = x'*B*x/2
    function bprod(x, B::AbstractArray{T,3}) where T
        s = zero(T)
        for I in CartesianIndices(size(B))
            p = oneunit(T)
            for d in I.I
                p *= x[d]
            end
            s += p*B[I]
        end
        return s/6
    end

    for p in (1, 2, 3)
        for n in (6, 7)
            x0 = zeros(n)
            world = World(fill(-Inf, n), fill(Inf, n), x0, 0.0)
            sz = ntuple(d->n, p)
            B = CoordinateSplittingPTrees.SymmetricArray(randn(sz))
            f(x) = bprod(x, B)
            root = Box{p}(world)
            filled = CoordinateSplittingPTrees.SymmetricArray(falses(sz))
            if p > 1
                for I in CartesianIndices(sz)
                    if anydups(I.I)
                        filled[I] = true
                    end
                end
            end
            # Generate splits
            while !all(filled)
                local pairings
                while true
                    pairings = dimpair(randperm(n), p)
                    isnovel = true
                    for pr in pairings
                        length(pr) == p || continue
                        isnovel = !filled[pr...]
                    end
                    isnovel && break
                end
                for pr in pairings
                    length(pr) == p || break
                    filled[pr...] = true
                end
                addpoint!(root, randn(n), pairings, f)
            end
            for box in leaves(root)
                Cp = CoordinateSplittingPTrees.coefficients_p(box)
                testequal_offdiag(Cp, B)
            end
        end
    end
end

hesscount = 0

@testset "Sparse Hessians" begin
    function fsparse(x)
        dx = diff(x)       # only (i,i+1) Q terms are nonzero
        return sum(abs2, dx)/2
    end
    x0 = [0.1,0.2,0.3,0.4]
    world = World(x0, fsparse)
    root = Box{2}(world)
    box = addpoint!(root, 2*x0, [(1,3), (2,4)], fsparse) # neither split has nonzero in Hessian
    box = addpoint!(root, 3*x0, [(2,3), (1,4)], fsparse)
    box = addpoint!(root, 4*x0, [(1,4), (2,3)], fsparse)
    box = addpoint!(root, 5*x0, [(1,2), (3,4)], fsparse)
    callback = (box,val)->(global hesscount; hesscount+=1)
    Q = zeros(4, 4)
    @test_throws MethodError CoordinateSplittingPTrees.coefficients_p!(Q, box, callback)  # type must enforce symmetry
    Q = CoordinateSplittingPTrees.SymmetricArray(Q)
    CoordinateSplittingPTrees.coefficients_p!(Q, box, callback)
    for i = 1:3
        @test Q[i,i+1] ≈ -1
        @test Q[i+1,i] ≈ -1
        @test isnan(Q[i,i])
    end
    @test abs((Q[1,3])) < 1e-8
    @test abs((Q[2,4])) < 1e-8
    @test abs((Q[1,4])) < 1e-8
    global hesscount
    @test hesscount == 24
    Q = SymTridiagonal(zeros(4), zeros(3))
    hesscount = 0
    CoordinateSplittingPTrees.coefficients_p!(Q, box, callback)
    for i = 1:3
        @test Q[i,i+1] ≈ -1
        @test Q[i+1,i] ≈ -1
        @test isnan(Q[i,i])
    end
    @test hesscount == 12
    Q = sparse([1,1,2,2,3,3,4,1], [1,2,2,3,3,4,4,4], zeros(8))
    @test_throws MethodError CoordinateSplittingPTrees.coefficients_p!(Q, box, callback)
    Q = CoordinateSplittingPTrees.SymmetricArray(Q)
    @test CoordinateSplittingPTrees.nnz_offdiag_sym(Q) == 4
    hesscount = 0
    CoordinateSplittingPTrees.coefficients_p!(Q, box, callback)
    for i = 1:3
        @test Q[i,i+1] ≈ -1
        @test Q[i+1,i] ≈ -1
        @test isnan(Q[i,i])
    end
    @test hesscount == 16
end

@testset "Choosing split dimensions" begin
    f(x) = rand()
    n = 2
    x0 = randn(n)
    world = World(fill(-Inf, n), fill(Inf, n), x0, f(x0))
    root = Box{2}(world)
    pairs = CoordinateSplittingPTrees.choose_dimensions(root)
    @test pairs == [(1,2)]
    box = addpoint!(root, randn(n), pairs, f)
    pairs = CoordinateSplittingPTrees.choose_dimensions(box)
    @test pairs == [(1,2)]

    n = 3
    x0 = randn(n)
    world = World(fill(-Inf, n), fill(Inf, n), x0, f(x0))
    pairseq = []
    root = Box{2}(world)
    pairs = CoordinateSplittingPTrees.choose_dimensions(root)
    @test pairs == [(1,2), (3,)]
    push!(pairseq, pairs)
    box = addpoint!(root, randn(n), pairs, f)
    pairs = CoordinateSplittingPTrees.choose_dimensions(box)
    push!(pairseq, pairs)
    @test pairs[1] ∈ ((1,3), (2,3))
    box = addpoint!(root, randpoint_inside(box), pairs, f)
    pairs2 = CoordinateSplittingPTrees.choose_dimensions(box)
    push!(pairseq, pairs2)
    @test pairs2[1] ∈ ((1,3), (2,3))
    @test pairs2[1] != pairs[1]
    box = addpoint!(root, randpoint_inside(box), pairs2, f)
    for i = 1:3
        pairs = CoordinateSplittingPTrees.choose_dimensions(box)
        @test pairs == pairseq[i]
        box = addpoint!(root, randpoint_inside(box), pairs, f)
    end

    root = Box{2}(world)
    pairs = CoordinateSplittingPTrees.choose_dimensions(root)
    box = addpoint!(root, randn(n), pairs, f)
    pairs = CoordinateSplittingPTrees.choose_dimensions(box)
    boxd = addpoint!(root, randpoint_inside(box), pairs, f)
    boxp = box.parent.split.others.children[2]
    pairs2 = CoordinateSplittingPTrees.choose_dimensions(boxp)
    @test pairs2[1] ∈ ((1,3), (2,3))
    @test pairs2[1] != pairs[1]
    boxpd = addpoint!(root, randpoint_inside(boxp), pairs2, f)
    @test CoordinateSplittingPTrees.choose_dimensions(boxpd) ==
          CoordinateSplittingPTrees.choose_dimensions(boxd)  == [(1,2), (3,)]
end

@testset "Qcoef" begin
    using CoordinateSplittingPTrees: Qcoef_value, Qcoef_lineq, chi
    for n in (7, 8)
        dims = randperm(n)
        seen = falses(n)
        prec = falses(n, n)
        for i = 1:2:n
            d1 = d2 = dims[i]
            prec[:,d1] = seen
            if i < n
                d2 = dims[i+1]
                prec[:,d2] = seen
            end
            seen[d1] = true
            if i < n
                seen[d2] = true
            end
        end
        b = randn(n)
        y = randn(n)
        for i = 1:n  # ensure y is distinct from b in all coordinates
            while y[i] == b[i]
                y[i] = randn()
            end
        end
        x = copy(b)
        for i = 1:2:n
            d1 = d2 = dims[i]
            Δx1 = Δx2 = y[d1] - b[d1]
            if i < n
                d2 = dims[i+1]
                Δx2 = y[d2] - b[d2]
            end
            for k = 1:n
                (k == d1 || k == d2) && continue
                @test Qcoef_lineq(d1, k, Δx1, x[k], b[k], y[k], prec[d1,k]) == 0
                @test Qcoef_lineq(d2, k, Δx2, x[k], b[k], y[k], prec[d2,k]) == 0
            end
            x[d1] = y[d1]
            if i < n
                x[d2] = y[d2]
            end
            for k = 1:n
                (k == d1 || k == d2) && continue
                @test Qcoef_value(k, d1, x, b, y, prec) == 0
                @test Qcoef_value(d1, k, x, b, y, prec) == 0
                @test Qcoef_value(k, d2, x, b, y, prec) == 0
                @test Qcoef_value(d2, k, x, b, y, prec) == 0
                @test abs((x[d1] - b[d1]) * (x[k] - chi(k, d1, b, y, prec))/2 +
                          (x[k] - b[k]) * (x[d1] - chi(d1, k, b, y, prec))/2) < 1e-12
                @test abs((x[d2] - b[d2]) * (x[k] - chi(k, d2, b, y, prec))/2 +
                          (x[k] - b[k]) * (x[d2] - chi(d2, k, b, y, prec))/2) < 1e-12
            end
            @test Qcoef_value(d1, d1, x, b, y, prec) == 0
            @test Qcoef_value(d2, d2, x, b, y, prec) == 0
            @test abs((x[d1] - b[d1]) * (x[d1] - chi(d1, d1, b, y, prec))/2 +
                      (x[d1] - b[d1]) * (x[d1] - chi(d1, d1, b, y, prec))/2) < 1e-12
            @test abs((x[d2] - b[d2]) * (x[d2] - chi(d2, d2, b, y, prec))/2 +
                      (x[d2] - b[d2]) * (x[d2] - chi(d2, d2, b, y, prec))/2) < 1e-12
            @test Qcoef_value(d1, d2, x, b, y, prec) ≈ (
                (x[d1] - b[d1]) * (x[d2] - chi(d2, d1, b, y, prec))/4 +
                (x[d2] - b[d2]) * (x[d1] - chi(d1, d2, b, y, prec))/4)
        end
        x = randn(n)
        for j = 1:n, i = 1:n
            @test Qcoef_value(i, j, x, b, y, prec) ≈ (
                (x[i] - b[i]) * (x[j] - chi(j, i, b, y, prec))/4 +
                (x[j] - b[j]) * (x[i] - chi(i, j, b, y, prec))/4)
        end
    end
end

@testset "CS2 full models" begin
    for n = (2, 3, 7, 8)
        B = CoordinateSplittingPTrees.SymmetricArray(randn(n, n))
        f(x) = (x'*B*x)/2
        x0 = randn(n)
        world = World(fill(-Inf, n), fill(Inf, n), x0, f(x0))
        for iter = 1:10  # ensure a variety of splitting patterns
            root = box = Box{2}(world)
            while true
                box = addpoint!(root, randn(n), f)
                Cp = CoordinateSplittingPTrees.coefficients_p(box)
                !hasnan_offdiag(Cp) && break
            end
            # One extra to make sure we can determine the diagonal components
            box = addpoint!(root, randn(n), f)
            c, g, Q, b, y, prec, firstmatch = CoordinateSplittingPTrees.fit_quadratic_native(box)
            for leaf in leaves(root)
                x = position(leaf)
                @test modelvalue(c, g, Q, b, y, prec, x) ≈ f(x) rtol=1e-7
                @test modelvalue_chi(c, g, Q, b, y, prec, x) ≈ f(x) rtol=1e-7
            end
            c, g, Q, b = CoordinateSplittingPTrees.fit_quadratic(box)
            @test Q ≈ B
            @test g ≈ Q*b
            x = position(box)
            @test c + g'*(x - b) + (x-b)'*Q*(x-b)/2 ≈ f(x) rtol=1e-8
            @test abs(c - g'*b + (b'*Q*b)/2) < 1e-8
            # gdiag approach (which can yield different answers for non-quadratic functions,
            # but should be identical for quadratic)
            Q = CoordinateSplittingPTrees.coefficients_p(box)
            c, g, Q, b = CoordinateSplittingPTrees.fit_quadratic_gdiag!(Q, box)
            @test Q ≈ B
            @test g ≈ Q*b
        end
    end
end

@testset "Lowerbound models" begin
    f(x) = sum(abs2, x)/2 + x[1]*x[2]/5
    world = World([1.0,1.0,1.0,1.0], f)
    n = ndims(world)
    root = Box{2}(world)
    for i = 1:n+2
        addpoint!(root, randn(n), f)
    end
    box = minimum(root)

    updater, model, gp, dp, Q = CoordinateSplittingPTrees.lowerbound_model(box, GLPK.Optimizer())
    c, g, Q, b = CoordinateSplittingPTrees.fit_quadratic_lowerbound!(Q, updater, model, gp, dp, box)
    Qtarget = Matrix(Diagonal(ones(n)))
    Qtarget[1,2] = Qtarget[2,1] = 0.2
    @test Q ≈ Qtarget rtol=0.01
    @test g ≈ Q*b rtol=0.01

    # Check that it works with an excess of points
    updater, model, gp, dp, Q = CoordinateSplittingPTrees.lowerbound_model(box, GLPK.Optimizer(), 3n)
    c, g, Q, b = CoordinateSplittingPTrees.fit_quadratic_lowerbound!(Q, updater, model, gp, dp, box)
    Qtarget = Matrix(Diagonal(ones(n)))
    Qtarget[1,2] = Qtarget[2,1] = 0.2
    @test Q ≈ Qtarget rtol=0.01
    @test g ≈ Q*b rtol=0.01

    # Check that nothing breaks if we ask for more constraints than we have boxes
    updater, model, gp, dp, Q = CoordinateSplittingPTrees.lowerbound_model(box, GLPK.Optimizer(), length(leaves(root))+1)
    c, g, Q, b = CoordinateSplittingPTrees.fit_quadratic_lowerbound!(Q, updater, model, gp, dp, box)
    Qtarget = Matrix(Diagonal(ones(n)))
    Qtarget[1,2] = Qtarget[2,1] = 0.2
    @test Q ≈ Qtarget rtol=0.01
    @test g ≈ Q*b rtol=0.01
end

@testset "Gauss elim" begin
    ige = CoordinateSplittingPTrees.IGE{Float64}(3)
    insert!(ige, [0,1,0], 1)
    insert!(ige, [1,0,0], 2)
    insert!(ige, [0,0,1], 3)
    x = CoordinateSplittingPTrees.solve(ige)
    @test x ≈ [2,1,3]
    empty!(ige)
    U, _ = qr(randn(3,3))  # random unitary matrix
    A = U*Diagonal(range(1, stop=3, length=3))*U'  # random posdef with eigvals 1,2,3
    rhs = randn(3)
    xtrue = A\rhs
    for i = 1:3
        insert!(ige, A[i,:], rhs[i])
        if i < 3
            @test !iscomplete(ige)
            @test_throws LinearAlgebra.LAPACKException CoordinateSplittingPTrees.solve(ige)
        end
    end
    @test iscomplete(ige)
    x = CoordinateSplittingPTrees.solve(ige)
    @test x ≈ xtrue
    insert!(ige, A[3,:], 0)
    x = CoordinateSplittingPTrees.solve(ige)
    @test x ≈ xtrue
    A = U*Diagonal([1.0,1.0,1e-12])*U'  # random posdef one small eigval
    empty!(ige)
    for i = 1:3
        insert!(ige, A[i,:], rhs[i])
    end
    x = CoordinateSplittingPTrees.solve(ige)
    @test x ≈ A \ rhs rtol=1e-4
end

@testset "posdef" begin
    Q = Float64[4 1; 1 4]
    Qp = CoordinateSplittingPTrees.possemidef(Q)
    @test Qp == Q
    Q = Float64[1 4; 4 1]
    Qp = CoordinateSplittingPTrees.possemidef(Q)
    @test Qp == [4 4; 4 4]
    Q = Float64[4.5 3; 3 0.5]
    Qp = CoordinateSplittingPTrees.possemidef(Q)
    @test Qp == [9 3; 3 1]
    # This case is deliberately weaker in what it tests, since the
    # desired behavior for 0 diagonal entries is less clear.
    Q = Float64[0 2; 2 1]
    Qp = CoordinateSplittingPTrees.possemidef(Q)
    @test all(isfinite, Qp) && Qp[1,2] == Q[1,2]
    cholesky(Qp) # test only that it's positive-definite
end

@testset "Display" begin
    io = IOBuffer()
    io2 = IOBuffer()
    # CS1, 1-d
    world = World([0], [1], [0], 10)
    root = Box{1}(world)
    box = Box(root, 1, 1, 20)[1]
    box1 = Box(box, 1, 0.5, 30)[1]
    splitprint(io, root)
    @test String(take!(io)) == "(1,)[l, (1,)[l, l]]"
    splitprint_colored(io, root, box1)
    printstyled(io2, "(1,)", color=:cyan)
    print(io2, "[l, ")
    printstyled(io2, "(1,)", color=:cyan)
    print(io2, "[l, ")
    printstyled(io2, "l", color=:light_red)
    print(io2, "]]")
    @test String(take!(io)) == String(take!(io2))
    print_tree(io, root)
    @test String(take!(io)) == """
(1,)
├─ Box10@[0.0]
└─ (1,)
   ├─ Box20@[1.0]
   └─ Box30@[0.5]
"""

    # CS1, 2-d
    world = World([0,-Inf], [Inf,Inf], [1,1], 'A')
    root = Box{1}(world)
    box1 = Box(root, 1, 2, 'B')[1]
    box2 = Box(getleaf(root), 2, 2, 'C')[1]
    print_tree(io, root)
    @test String(take!(io)) == """
(1,)
├─ (2,)
│  ├─ Box'A'@[1.0, 1.0]
│  └─ Box'C'@[1.0, 2.0]
└─ Box'B'@[2.0, 1.0]
"""

    # CS2, 3-d
    world = World(fill(-Inf,3), fill(Inf,3), fill(1, 3), 100)
    root = Box{2}(world)
    boxes1 = Box(root, (1,2), (10,20), (200,300,400))
    boxes2 = Box(boxes1[2], (3,1), (-10,3.5), (500,600,700))
    print_tree(io, root)
    @test String(take!(io)) == """
(1, 2)
├─ Box100@[1.0, 1.0, 1.0]
├─ Box200@[10.0, 1.0, 1.0]
├─ (3, 1)
│  ├─ Box300@[1.0, 20.0, 1.0]
│  ├─ Box500@[1.0, 20.0, -10.0]
│  ├─ Box600@[3.5, 20.0, 1.0]
│  └─ Box700@[3.5, 20.0, -10.0]
└─ Box400@[10.0, 20.0, 1.0]
"""
end

@testset "Box metadata" begin
    world = World([1, 1], 10.0)
    root = Box{2}(world, (-1, 1))
    local incrementor
    let counter = 1
        incrementor(box) = (counter += 1; (-counter,counter))
    end
    addpoint!(root, [2.0, 1.8], x->sum(x), incrementor)
    @test root.metabox == (-1, 1)
    i = 1
    for leaf in leaves(root)
        i += 1
        @test leaf.metabox == (-i, i)
    end
    io = IOBuffer()
    print_metabox(io, root)
    @test String(take!(io)) == "Box10.0@[1.0, 1.0] with metadata (-1, 1)"
end


text:
mutable struct UpdatableCholesky{T, UT <: AbstractMatrix{T}} <: Factorization{T}
    U_full::UT # pre-allocated U matrix
    n::Int # size of factorization
    info::Int
end
const UCholesky = UpdatableCholesky

"""
```
updatable_cholesky(A::AbstractMatrix, m::Int = 2size(A, 1); check::Bool = true)
```
Computes an updatable Cholesky factorization starting with the `n by n` matrix `A`,
and pre-allocates an `m x m` matrix for future additions to the matrix (defaults to twice the size `2n`).
"""
function updatable_cholesky(A::AbstractMatrix, m::Int = 2size(A, 1); check::Bool = true)
    n = checksquare(A)
    m ≥ n || throw(DimensionMismatch("$m < size(A, 1) == $n"))
    U_full = zeros(eltype(A), m, m)
    @. U_full[1:n, 1:n] = A
    updatable_cholesky!(U_full, n, check = check)
end
const ucholesky = updatable_cholesky

"""
```
updatable_cholesky!(U_full::AbstractMatrix, n::Int; check::Bool = true)
```
Computes an updatable Cholesky factorization. Assumes that the currently available
data `A` is in the upper left `n x n` submatrix.
"""
function updatable_cholesky!(U_full::AbstractMatrix, n::Int; check::Bool = true)
    m = checksquare(U_full)
    A = @view U_full[1:n, 1:n]
    C = cholesky!(A, check = check)
    UpdatableCholesky(U_full, n, C.info)
end
const ucholesky! = updatable_cholesky!

# constructing cholesky from UpdatableCholesky
function LinearAlgebra.Cholesky(F::UpdatableCholesky)
    Cholesky(F.U, :U, F.info)
end

# converting Cholesky factorization to UpdatableCholesky
function updatable_cholesky(C::Union{Cholesky, CholeskyPivoted}, m::Int = size(C, 1); check::Bool = true)
	U_full = zeros(eltype(C), m, m) # in principle, could do this in place
	n = size(C, 1)
	U = C.U
	if C isa CholeskyPivoted
		U = @view C.U[:, invperm(C.p)]
	end
	@. U_full[1:n, 1:n] = U
	UpdatableCholesky(U_full, n, C.info)
end

# fall-back for arbitrary factorization object
function updatable_cholesky(F::Factorization, m::Int = size(F, 1); check::Bool = true)
	updatable_cholesky(Hermitian(Matrix(F)), m, check = check)
end

function Base.getproperty(F::UpdatableCholesky, s::Symbol)
    if s == :U
        U = @view F.U_full[1:F.n, 1:F.n]
        UpperTriangular(U)
    elseif s == :L
        F.U'
    else
        getfield(F, s)
    end
end
Base.eltype(F::UpdatableCholesky{T}) where T = T
Base.size(F::UpdatableCholesky) = (F.n, F.n)
Base.size(F::UpdatableCholesky, i::Int) = i > 2 ? 1 : size(F)[i]
Base.AbstractMatrix(F::UpdatableCholesky) = Matrix(F)
Base.adjoint(C::UpdatableCholesky) = C # since it is Hermitian
LinearAlgebra.issuccess(C::UpdatableCholesky) = C.info == 0
function Base.Matrix(F::UpdatableCholesky)
    U = F.U
    U'U
end

function Base.copy(F::UpdatableCholesky)
    UpdatableCholesky(copy(F.U_full), F.n, F.info)
end

# linear solves with updatable cholesky
LinearAlgebra.:\(C::UpdatableCholesky, x::AbstractVector) = ldiv!(C, copy(x))
LinearAlgebra.:\(C::UpdatableCholesky, X::AbstractMatrix) = ldiv!(C, copy(X))
LinearAlgebra.:/(X::AbstractMatrix, C::UpdatableCholesky) = rdiv!(copy(X), C)
LinearAlgebra.rdiv!(X::AbstractMatrix, F::UpdatableCholesky) = ldiv!(F, X')'
function LinearAlgebra.ldiv!(F::UpdatableCholesky, x::AbstractVecOrMat)
    U = F.U
    ldiv!(U', x)
    ldiv!(U, x)
    return x
end

"""
```
add_column!(C::UpdatableCholesky, A::AbstractMatrix)
```
appending columns in A - and their corresponding rows A' - to `UpdatableCholesky`
 factorization `C`. Complexity is `O(m³ + n²m)` where `n = size(C, 1)` and `m = size(A, 2)`.
"""
function add_column!(C::UpdatableCholesky, A::AbstractMatrix)
    ensure_enough_memory!(C, A)
    n, m = size(C, 1), size(A, 2)
    new_n = n + m
    if size(A, 1) ≠ new_n
        throw(DimensionMismatch("New columns `A` of $(size(A)) do not have the required dimensionality $new_n to add to UpdatableCholesky of size $n."))
    end
    A12 = @view A[1:n, :]
    A22 = @view A[n+1:end, :]

    S12 = @view C.U_full[1:n, n+1:new_n]
    ldiv!(S12, C.U', A12) # this is squared in the number of old indices and linear in m: O(n²m)

    # S22 = cholesky!(A22 - S12' * S12).U
    S22 = @view C.U_full[n+1:new_n, n+1:new_n]
    @. S22 = A22
    mul!(S22, S12', S12, -1, 1)
    C22 = cholesky!(S22) # this is cubic in the number of new indices only `O(m³)`
    C.n = new_n
    return C
end
function add_column!(C::UpdatableCholesky, a::AbstractVector)
    add_column!(C, reshape(a, :, 1))
end

"""
```
remove_column!(C::UpdatableCholesky, i::Int)
```
updating `UpdatableCholesky` factorization `C` corresponding to removal of the
`i`th row and column of the original matrix. Complexity is `O(n²)`.
"""
function remove_column!(C::UpdatableCholesky, i::Int)
    n = C.n
    if i == n # removing the last row / column is trivial
        C.n -= 1
        return C
    end

    # Create vector for rank update
    t = zeros(eltype(C), n)
    # t = @view C.U[:, i] # IDEA: can we use existing memory? Problem: Base's lowrankupdate! does not accept views.
    @. t[1:i] = 0
    @. t[i+1:n] = conj.(C.U[i, i+1:n])
    # Do rank update
    D = Cholesky(C)
    lowrankupdate!(D, t)

    @. C.U[:, i] = 0 # shift memory to the left, IDEA: could instead update view, but this will be more complex with subsequent additions
    for j in i+1:n
        @. C.U[1:i-1, j-1] = C.U[1:i-1, j]
        @. C.U[i:j-1, j-1] = C.U[i+1:j, j]
    end
    C.n -= 1
    return C
end

const addcol! = add_column!
const remcol! = remove_column!

"""
```
ensure_enough_memory!(C::UpdatableCholesky, A::AbstractMatrix, memory_expansion_factor::Int = 2)
```
makes sure we have enough memory pre-allocated in U_full to accomodate
additions of columns in `A`. `memory_expansion_factor` determines how much memory
is requested to accomodate future column additions, if we don't have enough space
to add `A`.
"""
function ensure_enough_memory!(C::UpdatableCholesky, A::AbstractMatrix, memory_expansion_factor::Int = 2)
    if size(C.U_full, 1) < C.n + size(A, 2) # create a bigger U_full, and store it in C.U_full
        n = C.n
        m = memory_expansion_factor * n # new size
        U_full = zeros(eltype(C), m, m)
        @. U_full[1:n, 1:n] = C.U
        C.U_full = U_full
    end
    return C
end


text:
origin(c::Cylinder{N, T}) where {N, T} = c.origin
extremity(c::Cylinder{N, T}) where {N, T} = c.extremity
radius(c::Cylinder{N, T}) where {N, T} = c.r
height(c::Cylinder{N, T}) where {N, T} = norm(c.extremity - c.origin)
direction(c::Cylinder{N, T}) where {N, T} = (c.extremity .- c.origin) ./ height(c)

function rotation(c::Cylinder{2, T}) where T
    d2 = direction(c); u = @SVector [d2[1], d2[2], T(0)]
    v = @MVector [u[2], -u[1], T(0)]
    normalize!(v)
    return hcat(v, u, @SVector T[0, 0, 1])
end
function rotation(c::Cylinder{3, T}) where T
    d3 = direction(c); u = @SVector [d3[1], d3[2], d3[3]]
    if abs(u[1]) > 0 || abs(u[2]) > 0
        v = @MVector [u[2], -u[1], T(0)]
    else
        v = @MVector [T(0), -u[3], u[2]]
    end
    normalize!(v)
    w = @SVector [u[2] * v[3] - u[3] * v[2], -u[1] * v[3] + u[3] * v[1], u[1] * v[2] - u[2] * v[1]]
    return hcat(v, w, u)
end


text:
struct Triangle{T}
    vertices ::NTuple{3, Int}
    r        ::T
end

Triangle(vs, r::T) where T = Triangle{T}(vs, r)

function Base.show(io::IO, Δ::Triangle{T}) where T
    if T <: Number
        print(io, "Δ($(Δ.vertices), $(round(Δ.r, 5)))")
    else
        print(io, "Δ($(Δ.vertices), $(Δ.r))")
    end
end

# Return "empty" triangle. Used as a placeholder value when the iterator is done.
nulltriangle(::Type{T}) where T = Triangle{T}((0, 0, 0), zero(T))

# ============================================================================ #
# Equilateral Iterator
# ============================================================================ #

# Return a vector that maps ordering to indices in distance matrix, for example,
# ordertoedge(M)[3] returns the third shortest edge (=third smallest entry) in M
function ordertoedge(M)
    n = size(M, 1)
    res = Vector{Tuple{Int, Int}}(n * (n - 1) ÷ 2)
    i = 1
    for j in 2:n, k in 1:j-1
        res[i] = (j, k)
        i += 1
    end
    sort!(res, lt = (e1, e2) -> isless(M[e1...], M[e2...]))
    res
end

struct EquilateralIterator{T, M<:AbstractMatrix{T}}
    M     ::M
    edges ::Vector{Tuple{Int, Int}}
    tol   ::T
    upper ::T
end

"""
    equilaterals(M; tol = sampledensity(M), upper = Inf)

Iterator that returns approximately (up to `tol`) equilateral triangles in
distance matrix `M`, sorted by longest side length. Skips triangles with
longest side length above `upper`.
"""
function equilaterals(M::AbstractMatrix{T};
                      tol   = sampledensity(M),
                      upper = typemax(T)) where T
    t = convert(T, tol)
    u = convert(T, upper)
    edges = ordertoedge(M)
    EquilateralIterator{T, typeof(M)}(M, edges, t, u)
end

function Base.start(ei::EquilateralIterator{T}) where T
    M     = ei.M
    edges = ei.edges
    tol   = ei.tol
    upper = ei.upper

    @inbounds for eindex in eachindex(edges)
        i, j = edges[eindex]
        a = M[i, j]
        a > upper && return (true, nulltriangle(T), 0, 0)
        a ≤ 0 && continue

        for k in 1:size(M, 1)
            (k == i || k == j || i == j) && continue
            b = M[k, i]
            b ≤ 0 && continue
            c = M[k, j]
            c ≤ 0 && continue
            if 0 ≤ a - b ≤ tol && 0 ≤ a - c ≤ tol && abs(b - c) ≤ tol
                return (false, Triangle((i, j, k), a), eindex, k + 1)
            end
        end
    end
    (true, nulltriangle(T), 0, 0)
end

function Base.next(ei::EquilateralIterator{T}, st) where T
    M     = ei.M
    edges = ei.edges
    tol   = ei.tol
    upper = ei.upper

    _, Δ, eindex, kstart = st
    n = size(M, 1)

    if kstart > n
        eindex += 1
        kstart  = 1
    end

    if eindex > length(edges) || M[edges[eindex]...] > upper
        return Δ, (true, nulltriangle(T), 0, 0)
    end

    @inbounds while eindex ≤ length(edges)
        i, j = edges[eindex]
        a = M[i, j]
        a > upper && return Δ, (true, nulltriangle(T), 0, 0)
        a ≤ 0 && continue

        for k in kstart:n
            (k == i || k == j || i == j) && continue
            b = M[k, i]
            b ≤ 0 && continue
            c = M[k, j]
            c ≤ 0 && continue
            if 0 ≤ a - b ≤ tol && 0 ≤ a - c ≤ tol && abs(b - c) ≤ tol
                return Δ, (false, Triangle((i, j, k), a), eindex, k + 1)
            end
        end

        eindex += 1
        kstart  = 1
    end

    Δ, (true, nulltriangle(T), 0, 0)
end

Base.done(ei::EquilateralIterator, st) = st[1]

Base.iteratorsize(::EquilateralIterator)   = Base.SizeUnknown()
Base.iteratoreltype(::EquilateralIterator) = Base.HasEltype()

Base.eltype(::EquilateralIterator{T}) where T       = Triangle{T}
Base.eltype(::Type{EquilateralIterator{T}}) where T = Triangle{T}