text:
## FUNCTIONS TO INSPECT METADATA OF REGISTERED MODELS AND TO
## FACILITATE MODEL SEARCH

# sort and add to the model trait names:
const property_names = sort(MODEL_TRAITS_IN_REGISTRY)
const alpha = [:name, :package_name, :is_supervised]
const omega = [:input_scitype, :target_scitype, :output_scitype]
const both = vcat(alpha, omega)
filter!(!in(both), property_names)
prepend!(property_names, alpha)
append!(property_names, omega)
const PROPERTY_NAMES = Tuple(property_names)

ModelProxy = NamedTuple{PROPERTY_NAMES}

function Base.isless(p1::ModelProxy, p2::ModelProxy)
    if isless(p1.name, p2.name)
        return true
    elseif p1.name == p2.name
        return isless(p1.package_name, p2.package_name)
    else
        return false
    end
end

import MLJModelInterface.==
function ==(m1::ModelProxy, m2::ModelProxy)
    m1.name == m2.name && m1.package_name == m2.package_name
    # tests = map(keys(m1)) do k
    #     v1 = getproperty(m1, k)
    #     v2 = getproperty(m2, k)
    #     if k isa AbstractVector
    #         Set(v1) == Set(v2)
    #     else
    #         v1 == v2
    #     end
    # end
    # return all(tests)
end

Base.show(stream::IO, p::ModelProxy) =
    print(stream, "(name = $(p.name), package_name = $(p.package_name), "*
          "... )")

function Base.show(stream::IO, ::MIME"text/plain", p::ModelProxy)
    printstyled(stream, p.docstring, bold=false, color=:magenta)
    println(stream)
    PrettyPrinting.pprint(stream, p)
end

# returns named tuple version of the dictionary i=info_dict(SomeModelType):
function info_as_named_tuple(i)
    property_values = Tuple(i[property] for property in PROPERTY_NAMES)
    return NamedTuple{PROPERTY_NAMES}(property_values)
end

## INFO

StatisticalTraits.info(handle::Handle) =
    info_as_named_tuple(INFO_GIVEN_HANDLE[handle])

"""
    info(name::String; pkg=nothing, interactive=false)

Returns the metadata for the registered model type with specified
`name`. The key-word argument `pkg` is required in the case of
duplicate names.

"""
function StatisticalTraits.info(name::String; pkg=nothing, interactive=false)
    name in NAMES ||
        throw(ArgumentError("There is no model named \"$name\" in "*
                            "the registry. \n Run `models()` to view all "*
                            "registered models."))
    # get the handle:
    if pkg == nothing
        handle  = Handle(name) # returns (name=..., pkg=missing) if ambiguous
        if ismissing(handle.pkg)
            pkgs = PKGS_GIVEN_NAME[name]
            if interactive
                choice = request(
                    "Multiple packages provide $name. Choose a package: ",
                    pkgs...)
                handle = Handle(name, pkgs[choice])
            else
                message = "Ambiguous model type name. Use pkg=... .\n"*
                    "The model $name is provided by these packages:\n $pkgs.\n"
                throw(ArgumentError(message))
            end
        end
    else
        handle = Handle(name, pkg)
        haskey(INFO_GIVEN_HANDLE, handle) ||
            throw(ArgumentError("The package \"$pkg\" does not appear to "*
                                "provide the model \"$name\". \n"*
                                "Use models() to list all models. "))
    end
    return info(handle)

end

"""
   info(model::Model)

Return the traits associated with the specified `model`. Equivalent to
`info(name; pkg=pkg)` where `name::String` is the name of the model type, and
`pkg::String` the name of the package containing it.

"""
function StatisticalTraits.info(M::Type{<:MMI.Model})
    values =
        tuple([eval(:($trait($M))) for trait in PROPERTY_NAMES]...)
    info_as_named_tuple(info_dict(M))
    return NamedTuple{PROPERTY_NAMES}(values)
end
StatisticalTraits.info(model::MMI.Model) = info(typeof(model))

## MATCHING

# Note. `ModelProxy` is the type of a model's metadata entry (a named
# tuple). So, `info("PCA")` has this type, for example.

# Basic idea

if false

    matching(model::MLJModels.ModelProxy, X) =
        !(model.is_supervised) && scitype(X) <: model.input_scitype

    matching(model::MLJModels.ModelProxy, X, y) =
        model.is_supervised &&
        scitype(X) <: model.input_scitype &&
        scitype(y) <: model.target_scitype

    matching(model::MLJModels.ModelProxy, X, y, w::AbstractVector{<:Real}) =
        model.is_supervised &&
        model.supports_weights &&
        scitype(X) <: model.input_scitype &&
        scitype(y) <: model.target_scitype

end

# Implementation

struct Checker{is_supervised,
               supports_weights,
               supports_class_weights,
               input_scitype,
               target_scitype} end

function Base.getproperty(::Checker{is_supervised,
                                    supports_weights,
                                    supports_class_weights,
                                    input_scitype,
                                    target_scitype},
                          field::Symbol) where {is_supervised,
                                                supports_weights,
                                                supports_class_weights,
                                                input_scitype,
                                                target_scitype}
    if field === :is_supervised
        return is_supervised
    elseif field === :supports_weights
        return supports_weights
    elseif field === :supports_class_weights
        return supports_class_weights
    elseif field === :input_scitype
        return input_scitype
    elseif field === :target_scitype
        return target_scitype
    else
        throw(ArgumentError("Unsupported property. "))
    end
end

Base.propertynames(::Checker) =
    (:is_supervised,
     :supports_weights,
     :supports_class_weights,
     :input_scitype,
     :target_scitype)

function _as_named_tuple(s::Checker)
    names = propertynames(s)
    NamedTuple{names}(Tuple(getproperty(s, p) for p in names))
end

# function Base.show(io::IO, ::MIME"text/plain", S::Checker)
#     show(io, MIME("text/plain"), _as_named_tuple(S))
# end

not_missing_and_true(x) = !ismissing(x) && x

matching(X)       = Checker{false,missing,missing,scitype(X),missing}()
matching(X, y)    = Checker{true,missing,missing,scitype(X),scitype(y)}()
matching(X, y, w) = Checker{true,true,false,scitype(X),scitype(y)}()
matching(X, y, w::AbstractDict) =
    Checker{true,false,true,scitype(X),scitype(y)}()

(f::Checker{false,
            missing,
            missing,
            XS,
            missing})(model::MLJModels.ModelProxy) where XS =
    !(model.is_supervised) &&
    XS <: model.input_scitype

(f::Checker{true,
            missing,
            missing,
            XS,
            yS})(model::MLJModels.ModelProxy) where {XS,yS} =
    model.is_supervised &&
    XS <: model.input_scitype &&
    yS <: model.target_scitype

(f::Checker{true,
            true,
            false,
            XS,
            yS})(model::MLJModels.ModelProxy) where {XS,yS} =
    model.is_supervised &&
    not_missing_and_true(model.supports_weights) &&
    XS <: model.input_scitype &&
    yS <: model.target_scitype

(f::Checker{true,
            false,
            true,
            XS,
            yS})(model::MLJModels.ModelProxy) where {XS,yS} =
    model.is_supervised &&
    not_missing_and_true(model.supports_class_weights) &&
    XS <: model.input_scitype &&
    yS <: model.target_scitype

(f::Checker)(name::String; pkg=nothing) = f(info(name, pkg=pkg))
(f::Checker)(realmodel::Model) = f(info(realmodel))

matching(model::MLJModels.ModelProxy, args...) = matching(args...)(model)
matching(name::String, args...; pkg=nothing) =
    matching(info(name, pkg=pkg), args...)
matching(realmodel::Model, args...) = matching(info(realmodel), args...)

## MODEL QUERY

"""
    models()

List all models in the MLJ registry. Here and below *model* means the
registry metadata entry for a genuine model type (a proxy for types
whose defining code may not be loaded).

    models(filters..)

List all models `m` for which `filter(m)` is true, for each `filter`
in `filters`.

    models(matching(X, y))

List all supervised models compatible with training data `X`, `y`.

    models(matching(X))

List all unsupervised models compatible with training data `X`.

Excluded in the listings are the built-in model-wraps, like `EnsembleModel`,
`TunedModel`, and `IteratedModel`.

### Example

If

    task(model) = model.is_supervised && model.is_probabilistic

then `models(task)` lists all supervised models making probabilistic
predictions.

See also: [`localmodels`](@ref).

"""
function models(conditions...)
    unsorted = filter(info.(keys(INFO_GIVEN_HANDLE))) do model
        all(c(model) for c in conditions)
    end
    return sort!(unsorted)
end

"""
    models(needle::Union{AbstractString,Regex})

List all models whole `name` or `docstring` matches a given `needle`.
"""
function models(needle::Union{AbstractString,Regex})
    f = model ->
        occursin(needle, model.name) || occursin(needle, model.docstring)
    return models(f)
end

# get the model types in top-level of given module's namespace:
function localmodeltypes(modl; toplevel=false)
    ft = finaltypes(Model)
    return filter!(ft) do M
        name = MLJModelInterface.name(M)
        test1 = !toplevel || isdefined(modl, Symbol(name))
        !MLJModelInterface.is_wrapper(M) && test1
    end
end

"""
    localmodels(; modl=Main)
    localmodels(filters...; modl=Main)
    localmodels(needle::Union{AbstractString,Regex}; modl=Main)

List all models currently available to the user from the module `modl`
without importing a package, and which additional pass through the
specified filters. Here a *filter* is a `Bool`-valued function on
models.

Use `load_path` to get the path to some model returned, as in these
examples:

    ms = localmodels()
    model = ms[1]
    load_path(model)

See also [`models`](@ref), [`load_path`](@ref).

"""
function localmodels(args...; modl=Main, toplevel=false)
    modeltypes = localmodeltypes(modl, toplevel=toplevel)
    handles = map(modeltypes) do M
        Handle(MMI.name(M), MMI.package_name(M))
    end
    return filter(models(args...)) do model
        Handle(model.name, model.package_name) in handles
    end
end


text:
#
# Calculates requests per second from client to server without response, but
# using only naked sockets.
# Results: rps: 29303.170665295427, time: 34.12599992752075
#
server = listen(2000)
sock   = accept(server)

function run()
  for i=1:1000000 write(sock, UInt8(rand(1:100))) end
end

ext: jl

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
function updateλs(λⱼ, λⱼ₋₁, N2Fᵢ, N2Fⱼ, N2Fⱼ₋₁)
    # Fit a parabola for a line search of the minimum of the norm of F(xᵢ + λ δxᵢ) = p(λ)
    # modified from C.T.Kelley, 2003
    #
    # input:
    #       λⱼ = current λ
    #       λⱼ₋₁ = previous λ
    #       N2Fᵢ = |F(xᵢ)|²
    #       N2Fⱼ = |F(xᵢ + λⱼ d)|²
    #       N2Fⱼ₋₁ = |F(xᵢ + λⱼ₋₁ d)|²
    #
    # output:
    #       λⱼ₊₁, λⱼ = new values of λ's assuming a parabola
    #
    # internal parameters:
    #       σ₀ = 0.1 and σ₁ = 0.5, safeguard bounds for the linesearch

    # Set internal parameters.
    σ₀ = 0.1
    σ₁ = 0.5
    # Compute coefficients of interpolation polynomial.
    #   p(λ) = |F(xᵢ)|² + (c₁ λ + c₂ λ²) / d
    # with d = (λⱼ - λⱼ₋₁) λⱼ λⱼ₋₁ < 0.
    # If c₂ > 0 we have negative curvature and default to
    #   λⱼ₊₁ = σ₁ λ.
    # Otherwise return minimizes p(λ) but ensure that σ₀ < λⱼ₊₁ / λⱼ < σ₁
    c₂ = λⱼ₋₁ * (N2Fⱼ - N2Fᵢ) - λⱼ * (N2Fⱼ₋₁ - N2Fᵢ)
    if c₂ >= 0
        λⱼ₊₁ = σ₁ * λⱼ
        return λⱼ₊₁, λⱼ
    else
        c₁ = λⱼ^2 * (N2Fⱼ₋₁ - N2Fᵢ) - λⱼ₋₁^2 * (N2Fⱼ - N2Fᵢ)
        λⱼ₊₁ = - 0.5 * c₁ / c₂
        λⱼ₊₁ = min(max(λⱼ₊₁, σ₀ * λⱼ), σ₁ * λⱼ)
        return λⱼ₊₁, λⱼ
    end

end

function searchLineArmijo!(δxᵢ, xᵢ, Fᵢ, F, maxItArmijo, nrm, preprint = "")
    preprint ≠ "" ? preprint = preprint * "    │" : nothing
    j = 0 # Armijo iteration counter
    α = 1e-4 # Armijo parameter from C. T. Kelley [2003] (K03 hereafter)
    λⱼ₋₁, λⱼ = 1.0, 1.0 # Relative step size
    # Norm of F(xᵢ) never changes
    NFᵢ = nrm(Fᵢ)
    N2Fᵢ = NFᵢ^2
    # If norm at Newton step is already good enough, or if Newton step is already to small, accept the Newton step
    if typeof(xᵢ) ≠ typeof(δxᵢ)
        xᵢ = convert(typeof(δxᵢ), xᵢ) # required for type
    end
    xⱼ = xᵢ .+ δxᵢ # Newton step
    Fⱼ = F(xⱼ)
    NFⱼ = nrm(Fⱼ)
    if (NFⱼ < (1 - α * λⱼ) * NFᵢ) || (nrm(δxᵢ) / nrm(xᵢ) < 1e-10) || eltype(Fᵢ) ≠ Float64
        # Update xᵢ, and Fᵢ
        xᵢ .= xⱼ
        Fᵢ .= Fⱼ
        return δxᵢ, xᵢ, Fᵢ, j, (j > maxItArmijo)
    else
        # Otherwise, start with at the middle
        λⱼ = 0.5
        xⱼ .= xᵢ .+ λⱼ .* δxᵢ
        Fⱼ .= F(xⱼ)
        NFⱼ, NFⱼ₋₁ = nrm(Fⱼ),  NFⱼ
        # Initialize the squares of norms
        N2Fⱼ, N2Fⱼ₋₁ = NFⱼ^2, NFⱼ₋₁^2
        while (NFⱼ >= (1 - α * λⱼ) * NFᵢ) && (j <= maxItArmijo)
            # Update λ's using assuming a parabola (see K03)
            λⱼ, λⱼ₋₁ = updateλs(λⱼ, λⱼ₋₁, N2Fᵢ, N2Fⱼ, N2Fⱼ₋₁)
            # Update xⱼ, F(xⱼ), and norms
            xⱼ .= xᵢ .+ λⱼ .* δxᵢ
            Fⱼ .= F(xⱼ)
            NFⱼ = nrm(Fⱼ)
            if preprint ≠ ""
                j == 0 ? println("(Armijo line search)") : nothing
                print(preprint)
                RNδxⱼ = nrm(λⱼ .* δxᵢ) / nrm(xᵢ)
                @printf "%3d    %8.1e   %8.1e\n" j NFⱼ RNδxⱼ
            end
            N2Fⱼ, N2Fⱼ₋₁ = NFⱼ^2, N2Fⱼ
            j += 1
        end
        # Update δxᵢ, xᵢ, and Fᵢ
        δxᵢ .= λⱼ .* δxᵢ
        xᵢ .= xⱼ
        Fᵢ .= Fⱼ
        (preprint == "" || j ≤ maxItArmijo) ? nothing : println(preprint * "x ─> Armijo failure!")
        return δxᵢ, xᵢ, Fᵢ, j, (j > maxItArmijo)
    end
end

function updateJacobian!(JF, rShamᵢ, rSham₀, ArmijoFail, ∇ₓF, xᵢ; update_factors=true)
    if update_factors && rShamᵢ ≤ rSham₀ && ~ArmijoFail
        JF.age += 1
        JF.facage +=1
    else
        JF.fac = factorize(JF.fac, ∇ₓF(xᵢ), update_factors=update_factors)
        JF.age = 0
        update_factors ? JF.facage = 0 : JF.facage += 1
    end
    return JF
end

# Added inplace factorize for real-valued case
# Not usre this is the best way to do all this...
# Do not re-export this odd inplace factorize
import LinearAlgebra: factorize
function factorize(Mf::SuiteSparse.UMFPACK.UmfpackLU{T,<:Int}, M::SparseMatrixCSC{T,<:Int}; update_factors=true) where T<:Union{Float64,Complex{Float64}}
    update_factors ? factorize(M) : Mf
end

function NewtonChordShamanskii(F, ∇ₓF, nrm, xinit, τstop; preprint="", maxItNewton=50)
    if preprint ≠ ""
        println(preprint * "(No initial Jacobian factors fed to Newton solver)")
    end
    # Initial Jacobian
    J = ∇ₓF(xinit)
    Jfinit = factorize(J) # construct JF
    return NewtonChordShamanskii(F, ∇ₓF, nrm, xinit, τstop, Jfinit; preprint=preprint, maxItNewton=maxItNewton)
end

function initialize_AgedJacobianFactor(firstJ, Jfinit, update_factors)
    if update_factors # varying A
        JF = AgedJacobianFactors(Jfinit, 0, 0) # construct JF
    else              # constant A
        JF = AgedJacobianFactors(udpate_J_nonrealparts(firstJ, Jfinit), 0, 0)
    end
end

function udpate_J_nonrealparts(firstJ::SparseMatrixCSC{<:Dual, <:Int}, Jfinit)
    return DualFactors(Jfinit, dualpart.(firstJ))
end
function udpate_J_nonrealparts(firstJ::SparseMatrixCSC{<:Hyper, <:Int}, Jfinit)
    return HyperDualFactors(Jfinit, ε₁part.(firstJ), ε₂part.(firstJ), ε₁ε₂part.(firstJ))
end

function NewtonChordShamanskii(F, ∇ₓF, nrm, xinit, τstop, Jfinit; preprint="", maxItNewton=50, update_factors=true, firstJ=false)
    if preprint ≠ ""
        println(preprint * "Solving F(x) = 0 (using Shamanskii Method)")
        preprint_end = preprint * "└─> "
        preprint = preprint * "│   "
    end

    maxItArmijo = 20 # Max iterations in Armijo line search
     # Max quasi Newton iterations
    rSham₀ = 0.5     # Shamanskii minimum reduction to keep old J

    # Initial values for while loop
    Fᵢ = F(xinit)      ;  Fᵢ₋₁ = copy(Fᵢ)
    solType = eltype(Fᵢ) # The type of the solution will be that of J (or f)
    xᵢ₋₁ = convert(Vector{solType}, xinit)
    xᵢ = copy(xᵢ₋₁)
    Nxᵢ = nrm(xᵢ)   ;  Nxᵢ₋₁ = Nxᵢ
    NFᵢ = nrm(Fᵢ)   ;  NFᵢ₋₁ = NFᵢ
    δxᵢ = copy(xᵢ₋₁)

    # Initialize the first Jacobian
    if firstJ
        JF = AgedJacobianFactors(Jfinit, 10, 10)
    else
        JF = initialize_AgedJacobianFactor(∇ₓF(xinit), Jfinit, update_factors)
    end

    τᵢ = Nxᵢ / NFᵢ
    rShamᵢ = 0

    i = 0 # counter of quasi-Newton δxⱼs
    if preprint ≠ ""
        print(preprint)
        @printf "%3s   %8s   %8s   %7s   %7s" "iteration" "|F(x)|" "|δx|/|x|" "Jac age" "fac age"
        if solType == Complex{Float64}
            @printf "%12s" "|δix|/|ix|"
        elseif solType == Dual{Float64}
            @printf "%12s" "|δεx|/|εx|"
        elseif solType == Hyper{Float64}
            @printf "%12s" "|δhx|/|hx|"
        end
        println("")
        print(preprint)
        @printf "%5d       %8.1e   %8s   %7s   %7s\n" i NFᵢ "" "" ""
    end

    # Condition on relative tolerance for imaginary part of xᵢ
    nonrealtolx = true
    ArmijoFail = false

    # main iteration loop
    while (τᵢ < τstop || nonrealtolx) && (i < maxItNewton)
        i += 1
        if preprint ≠ ""
            print(preprint)
            @printf "%5d       " i
        end
        # fill containers with current iteration values
        xᵢ₋₁ .= xᵢ
        Fᵢ₋₁ .= Fᵢ
        Nxᵢ₋₁ = Nxᵢ
        NFᵢ₋₁ = NFᵢ

        # Update Jacobian (or rather its factors, JF)
        # Shamanski: only update J if |F(xᵢ)| / |F(xᵢ₋₁)| > rSham₀
        JF = updateJacobian!(JF, rShamᵢ, rSham₀, ArmijoFail, ∇ₓF, xᵢ, update_factors=update_factors)

        # Newton Step
        δxᵢ .= JF.fac \ -Fᵢ

        # Update δxᵢ, xᵢ, and Fᵢ with an Armijo line search
        δxᵢ, xᵢ, Fᵢ, nArmijo, ArmijoFail = searchLineArmijo!(δxᵢ, xᵢ, Fᵢ, F, maxItArmijo, nrm, preprint)

        Nxᵢ = nrm(xᵢ)
        NFᵢ = nrm(Fᵢ)
        τᵢ = Nxᵢ / NFᵢ
        rShamᵢ = NFᵢ / NFᵢ₋₁

        if ArmijoFail
            if JF.age > 0 # If Armijo fails and J is old, come back and update J
                xᵢ .= xᵢ₋₁
                Fᵢ .= Fᵢ₋₁
            else # else if J is fresh it's a complete failure
                preprint ≠ "" ? println("Complete Failure") : nothing
                return xᵢ
            end
        end

        Nδxᵢ = nrm(δxᵢ)
        RNδxᵢ = Nδxᵢ / Nxᵢ₋₁

        # Check that the non-real part is converging
        # Here I chose to only look at the relative step size, |δx|/|x|
        if solType == Float64
            nonrealtolx = false
        elseif solType == Complex{Float64}
            RNimδxᵢ = nrm(imag.(δxᵢ)) / nrm(imag.(xᵢ₋₁))
            nonrealtolx = (RNimδxᵢ > 1e-7)
        elseif solType == Dual{Float64}
            RNεδxᵢ = nrm(dualpart.(δxᵢ)) / nrm(dualpart.(xᵢ₋₁))
            nonrealtolx = (RNεδxᵢ > 1e-7)
        elseif solType == Hyper{Float64}
            RNhδxᵢ = nrm(ε₁part.(δxᵢ)) / nrm(ε₁part.(xᵢ₋₁)) +
                     nrm(ε₂part.(δxᵢ)) / nrm(ε₂part.(xᵢ₋₁)) +
                     nrm(ε₁ε₂part.(δxᵢ)) / nrm(ε₁ε₂part.(xᵢ₋₁))
            nonrealtolx = (RNhδxᵢ > 1e-7)
        end

        if preprint ≠ ""
            if ArmijoFail
                if JF.age > 0
                    print(preprint * "Armijo Failure, but old J")
                else
                    print(preprint * "Complete Armijo failure (new J): ")
                    @printf("|x|/|F(x)| = %.2g years\n", τᵢ/(365*24*60*60))
                end
            end
            nArmijo == 0 ? nothing : print(preprint * "    └─────> ") # alignment thing
            @printf "%8.1e   %8.1e   " NFᵢ RNδxᵢ
            print_marker(JF.age)
            print_marker(JF.facage)
            if solType == Complex{Float64}
                @printf "%8.1e" RNimδxᵢ
            elseif solType == Dual{Float64}
                @printf "%8.1e" RNεδxᵢ
            elseif solType == Hyper{Float64}
                @printf "%8.1e" RNhδxᵢ
            end
            println("")
        end
    end

    if preprint ≠ ""
        (τᵢ < τstop) && print(preprint_end * "Newton has reached max iterations, ")
        (τᵢ > τstop) && print(preprint_end * "Newton has converged, ")
        @printf("|x|/|F(x)| = %.2g years\n", τᵢ/(365*24*60*60))
    end

    return xᵢ
end

function print_marker(i)
    if i==0
        @printf " (!)%2d    " i
    else
        @printf "%6d    " i
    end
end


text:
using Documenter, MagneticReadHead

makedocs(;
    modules=[MagneticReadHead],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/oxinabox/MagneticReadHead.jl/blob/{commit}{path}#L{line}",
    sitename="MagneticReadHead.jl",
    authors="Lyndon White",
    assets=[],
)

deploydocs(;
    repo="github.com/oxinabox/MagneticReadHead.jl",
)


text:
using Distances

"""
    init_dist_matrix(data)
    init_dist_matrix(data, training_data)

initialize a `D_out` object for `dist_matrix!()`.
"""
function init_dist_matrix(data::AbstractArray{tp, N}) where {tp, N}
  T = size(data, 1)
  dat = zeros(tp, T, size(data, N))
  tdat = zeros(tp, size(data, N), T)
  D = zeros(tp, T, T)
  D_out = (D, dat, tdat)
  return(D_out)
end

function init_dist_matrix(data::AbstractArray{tp, N}, training_data::AbstractArray{tp, N}) where {tp, N}
  T = size(data, 1)
  Ttrain = size(training_data, 1)
  dat = zeros(tp, T, size(data, N))
  traindat = zeros(tp, Ttrain, size(training_data, N))
  tdat = zeros(tp, size(data, N), T)
  ttraindat = zeros(tp, size(training_data, N), Ttrain)
  D = zeros(tp, Ttrain, T)
  D_out = (D, dat, tdat, traindat, ttraindat)
  return(D_out)
end

"""
    dist_matrix!(D_out, data, ...)

compute the distance matrix of `data`, similar to `dist_matrix()`. `D_out` object has to be preallocated, i.e. with `init_dist_matrix`.

# Examples
```
julia> dc = randn(10,4, 4,3)
julia> D_out = init_dist_matrix(dc)
julia> dist_matrix!(D_out, dc, lat = 2, lon = 2)
julia> D_out[1]
```
"""
function dist_matrix!(D_out::Tuple{Array{tp,2},Array{tp,2},Array{tp,2}}, data::AbstractArray{tp, N}; dist::String = "Euclidean", space::Int = 0, lat::Int = 0, lon::Int = 0, Q = 0, dims = 2) where {tp, N}
  #@assert N == 2 || N == 3 || N  = 4
  (D, dat, tdat) = D_out
  if N == 2 copyto!(dat, data) end
  if N == 3 copyto!(dat, view(data, :, space, :)) end
  if N == 4 copyto!(dat, view(data, :, lat, lon, :))  end
  transpose!(tdat, dat)
  if(dist == "Euclidean")         pairwise!(D, Euclidean(), tdat, dims = dims)
  elseif(dist == "SqEuclidean")   pairwise!(D, SqEuclidean(), tdat, dims = dims)
  elseif(dist == "Chebyshev")     pairwise!(D, Chebyshev(), tdat, dims = dims)
  elseif(dist == "Cityblock")     pairwise!(D, Cityblock(), tdat, dims = dims)
  elseif(dist == "JSDivergence")  pairwise!(D, JSDivergence(), tdat, dims = dims)
  elseif(dist == "Mahalanobis")   pairwise!(D, Mahalanobis(Q), tdat, dims = dims)
  elseif(dist == "SqMahalanobis") pairwise!(D, SqMahalanobis(Q), tdat, dims = dims)
  else print("$dist is not a defined distance metric, has to be one of 'Euclidean', 'SqEuclidean', 'Chebyshev', 'Cityblock' or 'JSDivergence'")
  end
  return(D_out[1])
end

function dist_matrix!(D_out::Tuple{Array{tp,2},Array{tp,2},Array{tp,2},Array{tp,2},Array{tp,2}},
                             data::AbstractArray{tp, N}, training_data::AbstractArray{tp, N}; dist::String = "Euclidean", space::Int = 0, lat::Int = 0, lon::Int = 0, Q = 0, dims = 2) where {tp, N}
  #@assert N == 2 || N == 3 || N  = 4
  (D, dat, tdat, traindat, ttraindat) = D_out
  if N == 2 copyto!(dat, data) end
  if N == 3 copyto!(dat, view(data, :, space, :)) end
  if N == 4 copyto!(dat, view(data, :, lat, lon, :))  end
  if N == 2 copyto!(traindat, training_data) end
  if N == 3 copyto!(traindat, view(training_data, :, space, :)) end
  if N == 4 copyto!(traindat, view(training_data, :, lat, lon, :))  end
  transpose!(tdat, dat)
  transpose!(ttraindat, traindat)
  if(dist == "Euclidean")         pairwise!(D, Euclidean(), ttraindat, tdat, dims = dims)
  elseif(dist == "SqEuclidean")   pairwise!(D, SqEuclidean(), ttraindat, tdat, dims = dims)
  elseif(dist == "Chebyshev")     pairwise!(D, Chebyshev(), ttraindat, tdat, dims = dims)
  elseif(dist == "Cityblock")     pairwise!(D, Cityblock(), ttraindat, tdat, dims = dims)
  elseif(dist == "JSDivergence")  pairwise!(D, JSDivergence(), ttraindat, tdat, dims = dims)
  elseif(dist == "Mahalanobis")   pairwise!(D, Mahalanobis(Q), ttraindat, tdat, dims = dims)
  elseif(dist == "SqMahalanobis") pairwise!(D, SqMahalanobis(Q), ttraindat, tdat, dims = dims)
  else print("$dist is not a defined distance metric, has to be one of 'Euclidean', 'SqEuclidean', 'Chebyshev', 'Cityblock' or 'JSDivergence'")
  end
  return(D_out[1])
end

"""
    dist_matrix(data::AbstractArray{tp, N}; dist::String = "Euclidean", space::Int = 0, lat::Int = 0, lon::Int = 0, Q = 0) where {tp, N}
    dist_matrix(data::AbstractArray{tp, N}, training_data; dist::String = "Euclidean", space::Int = 0, lat::Int = 0, lon::Int = 0, Q = 0) where {tp, N}

compute the distance matrix of `data` or the distance matrix between data and training data i.e. the pairwise distances along the first dimension of data, using the last dimension as variables.
`dist` is a distance metric, currently `Euclidean`(default), `SqEuclidean`, `Chebyshev`, `Cityblock`, `JSDivergence`, `Mahalanobis` and `SqMahalanobis` are supported.
The latter two need a covariance matrix `Q` as input argument.

# Examples
```
julia> dc = randn(10, 4,3)
julia> D = dist_matrix(dc, space = 2)
```
"""
function dist_matrix(data::AbstractArray{tp, N}; dist::String = "Euclidean", space::Int = 0, lat::Int = 0, lon::Int = 0, Q = 0, dims::Int = 2) where {tp, N}
  D_out = init_dist_matrix(data)
  dist_matrix!(D_out, data, dist = dist, space = space, lat = lat, lon = lon ,Q = Q, dims = dims)
  return(D_out[1])
end

function dist_matrix(data::AbstractArray{tp, N}, training_data::AbstractArray{tp, N}; dist::String = "Euclidean", space::Int = 0, lat::Int = 0, lon::Int = 0, Q = 0, dims::Int = 2) where {tp, N}
  D_out = init_dist_matrix(data, training_data)
  dist_matrix!(D_out, data, training_data, dist = dist, space = space, lat = lat, lon = lon ,Q = Q, dims = dims)
  return(D_out[1])
end

"""
    knn_dists(D, k::Int, temp_excl::Int = 0)

returns the k-nearest neighbors of a distance matrix `D`. Excludes `temp_excl` (default: `temp_excl = 0`) distances
from the main diagonal of `D` to be also nearest neighbors.

# Examples
```
julia> dc = randn(20, 4,3)
julia> D = dist_matrix(dc, space = 2)
julia> knn_dists_out = knn_dists(D, 3, 1)
julia> knn_dists_out[5] # distances
julia> knn_dists_out[4] # indices
```
"""
function knn_dists(D::AbstractArray, k::Int, temp_excl::Int = 0)
    T = size(D,1)
    if ((k + temp_excl) > T-1) print("k has to be smaller size(D,1)") end
    knn_dists_out = init_knn_dists(T, k)
    knn_dists!(knn_dists_out, D, temp_excl)
    return(knn_dists_out)
end

"""
    init_knn_dists(T::Int, k::Int)
    init_knn_dists(datacube::AbstractArray, k::Int)

initialize a preallocated `knn_dists_out` object. `k`is the number of nerarest neighbors, `T` the number of time steps (i.e. size of the first dimension) or a multidimensional `datacube`.
"""
function init_knn_dists(T::Int, k::Int)
    ix = zeros(Int64, T)
    v = zeros(Float64, T)
    indices = zeros(Int64, T, k)
    nndists = zeros(Float64, T, k)
    knn_dists_out = (k, ix, v, indices, nndists)
    return(knn_dists_out)
end

function init_knn_dists(datacube::AbstractArray, k::Int)
    T = size(datacube, 1)
    ix = zeros(Int64, T)
    v = zeros(Float64, T)
    indices = zeros(Int64, T, k)
    nndists = zeros(Float64, T, k)
    knn_dists_out = (k, ix, v, indices, nndists)
    return(knn_dists_out)
end

"""
    knn_dists!(knn_dists_out, D, temp_excl::Int = 0)

returns the k-nearest neighbors of a distance matrix `D`. Similar to `knn_dists()`, but uses preallocated input object `knn_dists_out`, initialized with `init_knn_dists()`.
Please note that the number of nearest neighbors `k` is not necessary, as it is already determined by the `knn_dists_out` object.

# Examples
```
julia> dc = randn(20, 4,3)
julia> D = dist_matrix(dc, space = 2)
julia> knn_dists_out = init_knn_dists(dc, 3)
julia> knn_dists!(knn_dists_out, D)
julia> knn_dists_out[5] # distances
julia> knn_dists_out[4] # indices
```
"""
function knn_dists!(knn_dists_out::Tuple{Int64,Array{Int64,1},Array{Float64,1},Array{Int64,2},Array{Float64,2}}, D::AbstractArray, temp_excl::Int = 0)
    (k, ix, v, indices, nndists) = knn_dists_out
    T = size(D,1)
    if ((k + temp_excl) > T-1) print("k has to be smaller size(D,1)") end
    maxD = maximum(D)
    for i = 1:T
        copyto!(v, view(D,:,i))
        for excl = -temp_excl:temp_excl
          if(i+excl > 0 && i+excl <= T)
            v[i+excl]= maxD
          end
        end
        sortperm!(ix, v)
        for j = 1:k
            indices[i,j] = ix[j]
            nndists[i,j] = v[ix[j]]
        end
    end
    return(knn_dists_out)
end

# compute kernel matrix from distance matrix
"""
    kernel_matrix(D::AbstractArray, σ::Float64 = 1.0[, kernel::String = "gauss", dimension::Int64 = 1])

compute a kernel matrix out of distance matrix `D`, given `σ`. Optionally normalized by the `dimension`, if `kernel = "normalized_gauss"`.
compute `D` with `dist_matrix()`.

# Examples
```
julia> dc = randn(20, 4,3)
julia> D = dist_matrix(dc, space = 2)
julia> K = kernel_matrix(D, 2.0)
```
"""
function kernel_matrix(D::AbstractArray, σ::Float64 = 1.0, kernel::String = "gauss", dimension::Int64 = 1)
  K = similar(D)
  kernel_matrix!(K, D, σ, kernel, dimension)
  return K
end

# compute kernel matrix from distance matrix
"""
    kernel_matrix!(K, D::AbstractArray, σ::Float64 = 1.0[, kernel::String = "gauss", dimension::Int64 = 1])

compute a kernel matrix out of distance matrix `D`. Similar to `kernel_matrix()`, but with preallocated Array K (`K = similar(D)`) for output.

# Examples
```
julia> dc = randn(20, 4,3)
julia> D = dist_matrix(dc, space = 2)
julia> kernel_matrix!(D, D, 2.0) # overwrites distance matrix
```
"""
function kernel_matrix!(K::AbstractArray{T,N}, D::AbstractArray{T,N}, σ::Real = 1.0, kernel::String = "gauss", dimension::Int64 = 10) where {T,N}
    #if(size(D, 1) != size(D, 2)) print("D is not a distance matrix with equal dimensions")
    σ = convert(T, σ)
    if(kernel == "normalized_gauss") # k integral gets one
    for i in eachindex(K)
      @inbounds K[i] = exp.(-0.5 * D[i]./(σ*σ))./(σ * (2 *pi).^(dimension/2))#exp.(-0.5 * D[i]./(σ*σ))./((2 *pi*σ*σ).^(dimension/2))
    end
    elseif (kernel == "gauss")
    for i in eachindex(K)
        @inbounds K[i] = exp.(-0.5 * D[i]./(σ*σ))
    end
    end
  return(K)
end

###################################
#end


text:
#Subject
"""
    applylimitrule!(data::PKSubject, rule::LimitRule)

Apply rule to PK subject .

* STEP 1 (NaN step): replace all `NaN` and `missing` values with nan keyword value (if `nan` not NaN);
* STEP 2 (LLOQ step): replace values below `lloq` with `btmax` value if this value befor Tmax or with atmax if this value after Tmax (if `lloq` not NaN);
* STEP 3 (remove NaN): `rm` == true, then remove all `NaN` and `missing` values.
"""
function applylimitrule!(data::PKSubject, rule::LimitRule)
    applylimitrule!(data.time, data.obs, rule)
    data
end
"""
    applylimitrule!(f::Function, data::DataSet{T}, rule::LimitRule) where T <: PKSubject
"""
function applylimitrule!(f::Function, data::DataSet{T}, rule::LimitRule) where T <: PKSubject
    for i in data
        if f(i) applylimitrule!(i, rule) end
    end
    data
end
#DS ind Int
"""
    applylimitrule!(data::DataSet{T}, rule::LimitRule, ind::Int) where T <: PKSubject
"""
function applylimitrule!(data::DataSet{T}, rule::LimitRule, ind::Int) where T <: PKSubject
    applylimitrule!(data[ind], rule)
    data
end
#DS iter Int
"""
    applylimitrule!(data::DataSet{T}, rule::LimitRule, inds::Union{Vector{Int}, UnitRange{Int}, Tuple{Vararg{Int}}}) where T <: PKSubject
"""
function applylimitrule!(data::DataSet{T}, rule::LimitRule, inds::Union{Vector{Int}, UnitRange{Int}, Tuple{Vararg{Int}}}) where T <: PKSubject
    for i in inds
        applylimitrule!(data[i], rule)
    end
    data
end
#DS all
"""
    applylimitrule!(data::DataSet{T}, rule::LimitRule) where T <: PKSubject
"""
function applylimitrule!(data::DataSet{T}, rule::LimitRule) where T <: PKSubject
    for i = 1:length(data)
        applylimitrule!(data[i], rule)
    end
    data
end
#DS Dict
"""
    applylimitrule!(data::DataSet{T}, rule::LimitRule, sort::Dict) where T <: PKSubject
"""
function applylimitrule!(data::DataSet{T}, rule::LimitRule, sort::Dict) where T <: PKSubject
    for i = 1:length(data)
        if sort ⊆ data[i].id applylimitrule!(data[i], rule) end
    end
    data
end
"""
    applylimitrule!(time, obs, rule::LimitRule)
"""
function applylimitrule!(time, obs, rule::LimitRule)
    cmax, tmax, tmaxn = ctmax(time, obs)
    #NaN Rule
    obsn = length(obs)
    if !isnan(rule.nan)
        for i = 1:obsn
            if isnanormissing(obs[i])
                obs[i] = rule.nan
            end
        end
    end
    #LLOQ rule
    if !isnan(rule.lloq)
        for i = 1:obsn
            if !isnanormissing(obs[i]) && obs[i] <= rule.lloq
                if i <= tmaxn
                    obs[i] = rule.btmax
                else
                    obs[i] = rule.atmax
                end
            end
        end
    end
    #NaN Remove rule
    if rule.rm
        inds = findall(isnanormissing, obs)
        deleteat!(time, inds)
        deleteat!(obs, inds)
    end
    time, obs
end


text:
import Pkg
import UUIDs

struct AlwaysAssertionError <: Exception
    msg::String
end

abstract type CapStrategy end

struct UpperBound <: CapStrategy
end

struct MonotonicUpperBound <: CapStrategy
end

abstract type LatestVersionOption end

struct ExcludeLatestVersion <: LatestVersionOption
end

struct CapLatestVersion <: LatestVersionOption
end

struct Package
    name::String
end


text:
function getUpdatedFactoredVal_3D(HO::Bool,kappaSquaredLoc::Float64,T::Array{EikTtype,3},Done::Array{Bool,3},loc1::Int64,
									loc2::Int64, loc3::Int64, h::Array{Float64,1}, hinv::Array{Float64,1}, n::Array{Int64,1},
									src::Array{Int64,1}, T0loc::Float64, G0loc::Array{Float64,1}, OpLoc::Array{Int8,1},A::Array{Float64,1},B::Array{Float64,1})
# op: -1 for backward. 1 for forward.
# A and B are just containers for answers so we won't allocate tuples all the time....
if HO
	(t1,t1a) = calcDirectionHO_3D(1, T, Done, loc1, loc2, loc3, h, hinv, n, src, T0loc, G0loc, OpLoc,A,B);
	(t2,t2a) = calcDirectionHO_3D(2, T, Done, loc1, loc2, loc3, h, hinv, n, src, T0loc, G0loc, OpLoc,A,B);
	(t3,t3a) = calcDirectionHO_3D(3, T, Done, loc1, loc2, loc3, h, hinv, n, src, T0loc, G0loc, OpLoc,A,B);
else
	t1 = calcDirectionFO_3D(1, T, Done, loc1, loc2, loc3, h, hinv, n, src, T0loc, G0loc, OpLoc,A,B);
	t2 = calcDirectionFO_3D(2, T, Done, loc1, loc2, loc3, h, hinv, n, src, T0loc, G0loc, OpLoc,A,B);
	t3 = calcDirectionFO_3D(3, T, Done, loc1, loc2, loc3, h, hinv, n, src, T0loc, G0loc, OpLoc,A,B);
end

#### Piecewise Quadratic Solve
@inbounds A1 = A[1]; A2 = A[2]; A3 = A[3]; B1 = B[1]; B2 = B[2]; B3 = B[3];
T_loc = solvePiecewiseQuadratic(A1,B1,A2,B2,A3,B3,kappaSquaredLoc,OpLoc);

if FactoredEikonalFastMarching.imposeMonotonicity
	if (t1 > 0.0 && T_loc*T0loc < t1+1e-16)
		if OpLoc[1]==-1 || OpLoc[1]==1
			A1 = T0loc/h[1];
			B1 = t1/h[1];
		else
			A1 = (3.0*T0loc)/(2.0*h[1]);
			B1 = (4.0*t1 - t1a)/(2.0*h[1]);
		end
		# A1 = 0.0;
		# B1 = 0.0;
		# OpLoc[1] = 0;
	end
	if (t2 > 0.0 && T_loc*T0loc < t2+1e-16)
		if OpLoc[2]==-1 || OpLoc[2]==1
			A2 = T0loc/h[2];
			B2 = t2/h[2];
		else
			A2 = (3.0*T0loc)/(2.0*h[2]);
			B2 = (4.0*t2 - t2a)/(2.0*h[2]);
		end
		# A2 = 0.0;
		# B2 = 0.0;
		# OpLoc[2] = 0;
	end
	if (t3 > 0.0 && T_loc*T0loc < t3+1e-16)
		if OpLoc[3]==-1 || OpLoc[3]==1
			A3 = T0loc/h[3];
			B3 = t3/h[3];
		else
			A3 = (3.0*T0loc)/(2.0*h[3]);
			B3 = (4.0*t3 - t3a)/(2.0*h[3]);
		end
		#A3 = 0.0;
		#B3 = 0.0;
		#OpLoc[3] = 0;
	end
	T_loc = solvePiecewiseQuadratic(A1,B1,A2,B2,A3,B3,kappaSquaredLoc,OpLoc);
end
return T_loc;
end

function solvePiecewiseQuadratic(A1::Float64,B1::Float64,A2::Float64,B2::Float64,A3::Float64,B3::Float64,kappaSquaredLoc::Float64,OpLoc::Array{Int8,1})
a1 = A1*A1; a2 = A2*A2; a3 = A3*A3; b1 = -A1*B1; b2 = -A2*B2; b3 = -A3*B3; c1 = B1*B1; c2 = B2*B2; c3 = B3*B3;
T_loc = solveQuadratic(a1+a2+a3,b1+b2+b3,c1+c2+c3,kappaSquaredLoc);

g1 = (A1*T_loc - B1);
g2 = (A2*T_loc - B2);
g3 = (A3*T_loc - B3);

if g1<0.0 || g2<0.0 || g3<0.0
	val1 = (OpLoc[1]==0) ? Inf : B1/A1;
	val2 = (OpLoc[2]==0) ? Inf : B2/A2;
	val3 = (OpLoc[3]==0) ? Inf : B3/A3;
	# Below, solvePiecewiseQuadratic uses only the first two cells in OpLoc.
	if val1 < val2 # val1 is in for two.
		if val2 < val3 # solve for val1 and val2.
			T_loc = solvePiecewiseQuadratic(A1,B1,A2,B2,kappaSquaredLoc,a1+a2,b1+b2,c1+c2,OpLoc);
			OpLoc[3] = 0; #G_loc3 = 0.0;
		else # solve for val1 and val3.
			OpLoc[2] = OpLoc[3];
			T_loc = solvePiecewiseQuadratic(A1,B1,A3,B3,kappaSquaredLoc,a1+a3,b1+b3,c1+c3,OpLoc);
			OpLoc[3] = OpLoc[2]; # this is a hack to handle dimension 1 and 3 by a 2D code.
			OpLoc[2] = 0; #G_loc2 = 0.0;
		end
	else # val2 is in for two
		if val1 < val3 # solve for val2 and val1
 			T_loc = solvePiecewiseQuadratic(A1,B1,A2,B2,kappaSquaredLoc,a1+a2,b1+b2,c1+c2,OpLoc);
			OpLoc[3] = 0; #G_loc3 = 0.0;
		else # solve for val2 and val3.
		    OpLoc[1] = OpLoc[3];
			T_loc = solvePiecewiseQuadratic(A3,B3,A2,B2,kappaSquaredLoc,a3+a2,b3+b2,c3+c2,OpLoc);
			OpLoc[3] = OpLoc[1];
			OpLoc[1] = 0; #G_loc1 = 0.0;
		end
	end
end
return T_loc;
end

function calcDirectionHO_3D(dim::Int64,T::Array{EikTtype,3},Done::Array{Bool,3},
	loc1::Int64,loc2::Int64,loc3::Int64,h::Array{Float64,1},hinv::Array{Float64,1},n::Array{Int64,1},src::Array{Int64,1},
	T0loc::Float64,G0loc::Array{Float64,1},OpLoc::Array{Int8,1},A::Array{Float64,1},B::Array{Float64,1})
# OpLoc is used as an output to avoid tuple allocation.
if dim == 1
	loc1p1 = loc1 + 1;
	loc1m1 = loc1 - 1;
	loc1p2 = loc1 + 2;
	loc1m2 = loc1 - 2;
	loc2p1 = loc2m1 = loc2p2 = loc2m2 = loc2;
	loc3p1 = loc3m1 = loc3p2 = loc3m2 = loc3;
elseif dim == 2
	loc2p1 = loc2 + 1;
	loc2m1 = loc2 - 1;
	loc2p2 = loc2 + 2;
	loc2m2 = loc2 - 2;
	loc1p1 = loc1m1 = loc1p2 = loc1m2 = loc1;
	loc3p1 = loc3m1 = loc3p2 = loc3m2 = loc3;
else
	loc3p1 = loc3 + 1;
	loc3m1 = loc3 - 1;
	loc3p2 = loc3 + 2;
	loc3m2 = loc3 - 2;
	loc1p1 = loc1m1 = loc1p2 = loc1m2 = loc1;
	loc2p1 = loc2m1 = loc2p2 = loc2m2 = loc2;
end

tBWD = 0.0;
tFWD = 0.0;

ChooseBWD = getPaddedArr3D(Done,loc1m1,loc2m1,loc3m1);
if ChooseBWD
	tBWD = T[loc1m1,loc2m1,loc3m1]*T03D(n,h,src,loc1m1,loc2m1,loc3m1);
end
ChooseFWD = getPaddedArr3D(Done,loc1p1,loc2p1,loc3p1);
if ChooseFWD
	tFWD = T[loc1p1,loc2p1,loc3p1]*T03D(n,h,src,loc1p1,loc2p1,loc3p1);
end
if ChooseBWD && ChooseFWD
	ChooseFWD = tFWD < tBWD;
	ChooseBWD = !ChooseFWD;
end

A[dim] = 0.0;
B[dim] = 0.0;
OpLoc[dim] = 0;
t = T0loc*hinv[dim];
if ChooseBWD
	@inbounds if getPaddedArr3D(Done,loc1m2,loc2m2,loc3m2) && tBWD > T[loc1m2,loc2m2,loc3m2]*T03D(n,h,src,loc1m2,loc2m2,loc3m2)
		@inbounds A[dim] = 1.5*t + G0loc[dim]; #going for second order!!!
        @inbounds B[dim] = (2*T[loc1m1,loc2m1,loc3m1] - 0.5*T[loc1m2,loc2m2,loc3m2])*t;
        @inbounds OpLoc[dim] = -2;
		return (tBWD,T[loc1m2,loc2m2,loc3m2]*T03D(n,h,src,loc1m2,loc2m2,loc3m2))
	else # first order.
		@inbounds A[dim] = t + G0loc[dim];
		@inbounds B[dim] = T[loc1m1,loc2m1,loc3m1]*t;
		@inbounds OpLoc[dim] = -1; # Choosing N means backward derivative.
		return (tBWD,0.0)
	end
elseif ChooseFWD
	@inbounds if getPaddedArr3D(Done,loc1p2,loc2p2,loc3p2) && tFWD > T[loc1p2,loc2p2,loc3p2]*T03D(n,h,src,loc1p2,loc2p2,loc3p2)
		@inbounds A[dim] = 1.5*t - G0loc[dim];
        @inbounds B[dim] = (2*T[loc1p1,loc2p1,loc3p1] - 0.5*T[loc1p2,loc2p2,loc3p2])*t;
        @inbounds OpLoc[dim] = 2;
		return (tFWD,T[loc1p2,loc2p2,loc3p2]*T03D(n,h,src,loc1p2,loc2p2,loc3p2));
	else
		@inbounds A[dim] = t - G0loc[dim];
		@inbounds B[dim] = T[loc1p1,loc2p1,loc3p1]*t;
		@inbounds OpLoc[dim] = 1;
		return (tFWD,0.0);
	end
end
return (0.0,0.0);
end

function calcDirectionFO_3D(dim::Int64,T::Array{EikTtype,3},Done::Array{Bool,3},loc1::Int64,loc2::Int64,loc3::Int64,
h::Array{Float64,1},hinv::Array{Float64,1},n::Array{Int64,1},src::Array{Int64,1},T0loc::Float64,G0loc::Array{Float64,1},OpLoc::Array{Int8,1},A::Array{Float64,1},B::Array{Float64,1})
loc1p1 = loc1m1 = loc1;
loc2p1 = loc2m1 = loc2;
loc3p1 = loc3m1 = loc3;
if dim == 1
	loc1p1 += 1;
	loc1m1 -= 1;
elseif dim == 2 # dim is assumed to be 2
	loc2p1 += 1;
	loc2m1 -= 1;
else
	loc3p1 += 1;
	loc3m1 -= 1;
end

ChooseBWD = getPaddedArr3D(Done,loc1m1,loc2m1,loc3m1);
ChooseFWD = getPaddedArr3D(Done,loc1p1,loc2p1,loc3p1);
if ChooseBWD && ChooseFWD
	@inbounds ChooseFWD = T03D(n,h,src,loc1p1,loc2p1,loc3p1)*T[loc1p1,loc2p1,loc3p1] < T[loc1m1,loc2m1,loc3m1]*T03D(n,h,src,loc1m1,loc2m1,loc3m1);
	ChooseBWD = !ChooseFWD;
end

tBWD = 0.0;
tFWD = 0.0;

ChooseBWD = getPaddedArr3D(Done,loc1m1,loc2m1,loc3m1);
ChooseFWD = getPaddedArr3D(Done,loc1p1,loc2p1,loc3p1);

if ChooseBWD && ChooseFWD
	tBWD = T[loc1m1,loc2m1,loc3m1]*T03D(n,h,src,loc1m1,loc2m1,loc3m1);
	tFWD = T03D(n,h,src,loc1p1,loc2p1,loc3p1)*T[loc1p1,loc2p1,loc3p1];
	@inbounds ChooseFWD = tFWD < tBWD;
	ChooseBWD = !ChooseFWD;
elseif FactoredEikonalFastMarching.imposeMonotonicity
	if ChooseBWD
		tBWD = T[loc1m1,loc2m1,loc3m1]*T03D(n,h,src,loc1m1,loc2m1,loc3m1);
	end
	if ChooseFWD
		tFWD =  T03D(n,h,src,loc1p1,loc2p1,loc3p1)*T[loc1p1,loc2p1,loc3p1];
	end
end

A[dim] = 0.0;
B[dim] = 0.0;
OpLoc[dim] = 0;
t = T0loc*hinv[dim];
if ChooseBWD
	@inbounds A[dim] = t + G0loc[dim];
	@inbounds B[dim] = T[loc1m1,loc2m1,loc3m1]*t;
	@inbounds OpLoc[dim] = -1; # Choosing N means backward derivative.
	return tBWD;
elseif ChooseFWD
	@inbounds A[dim] = t - G0loc[dim];
	@inbounds B[dim] = T[loc1p1,loc2p1,loc3p1]*t;
	@inbounds OpLoc[dim] = 1;
	return tFWD;
end
return 0.0;
end


text:
@testset "CIC" begin

    nx = 64 # nombre de pts suivant x
    ny = 16   # nombre de pts suivant y

    alpha = 0.1
    kx = 0.5
    dimx = 2 * pi / kx
    dimy = 1

    mesh = TwoDGrid(dimx, nx, dimy, ny)

    dx = mesh.dx
    dy = mesh.dy

    rho = zeros(nx, ny)

    for i = 1:nx
        aux2 = alpha * cos(kx * mesh.x[i])
        for j = 1:ny
            rho[i, j] = aux2
        end
    end

    nbpart = 1000 * nx * ny

    p = ParticleGroup{2,2}(nbpart, charge = 1.0, mass = 1.0, n_weights = 1)

    sampler = LandauDamping(alpha, kx)

    ParticleInCell.sample!(p, mesh, sampler)

    kernel = CloudInCell()

    ρ = compute_rho(p, kernel, mesh)

    @test ρ ≈ rho atol=1e-2

end


text:
# Let block - makes local scope

a = let
  i = 1
  i += 3
end

println(a)
# println(i) - i not defined here

# Begin block - no local scope
c = begin
  i = 3
end

println(c)
println(i)

# Constants

const ONE = 1
println(ONE)

# Warnings thrown for each of these
# ONE = 1
# ONE = 2

# Julia does not allow this
# ONE = 1.0

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
using ExactD, Test, LinearAlgebra

@testset "parity operator" begin
    @test ExactD.parity(10, 1, 2) == 1.
    @test ExactD.parity(10, 2, 0) == -1.
    @test ExactD.parity(10, 1, 3) == 1.
    @test ExactD.parity(14, 1, 3) == -1.
end

@testset "hard-core boson many body operators" begin
@testset "dense many-body operator" begin
    L = 4
    N = 2
    J::Array{Float64, 2} = reshape(collect(1:L^2), (L, L))
    V::Array{Float64, 2} = reshape(collect(1:L^2), (L, L))
    Op = build_many_body_op(L, N, J, V)

    # Diagonal terms.
    @test Op[1, 1] ≈ 7. + 7.
    @test Op[4, 4] ≈ 17. + 17.
    @test Op[6, 6] ≈ 27. + 27.

    # Off-diagonal terms.
    @test Op[1, 2] ≈ 10.
    @test Op[2, 5] ≈ 0.
    @test Op[4, 2] ≈ 12.
    @test Op[6, 1] ≈ 0.
    @test Op[3, 5] ≈ 15.

    C = 1.3
    Op2 = build_many_body_op(L, N, J, V, C)

    # Diagonal terms.
    @test Op2[1, 1] ≈ 8.3 + 7.
    @test Op2[4, 4] ≈ 18.3 + 17.
    @test Op2[6, 6] ≈ 28.3 + 27.
end

@testset "sparse many-body operator" begin
    L = 6
    N = 3
    basis = make_LN_basis(L, N)
    J::Array{Float64, 2} = reshape(sin.(1:L^2), (L, L))
    V::Array{Float64, 2} = reshape(tan.(1:L^2), (L, L))
    C = 1.3
    sp_Op = build_sparse_many_body_op(L, basis, J, V, C)
    de_Op = build_many_body_op(L, basis, J, V, C)
    @test norm(Array(sp_Op) - de_Op) ≈ 0.0 atol=1e-12
end
end  # testset "hard-core boson many body operators"

@testset "dense many-body operator with L spins s=1" begin
    L = 4
    Sz = 2
    J = reshape(collect(1.:L^2), (L, L))
    W = reshape(collect(1.:L^2), (L, L))
    Jz = collect(1.:L)/2
    Jx = zeros(L)
    C = 1.33
    Op = build_spin1_many_body_op(L, Sz, J, W, Jz, Jx, C)

    # Diagonal terms: S^z + (S^z_i)^2 + S^z_i*S^z_j.
    @test Op[1, 1] ≈ 1.5 + 7 + 5 + C
    @test Op[4, 4] ≈ 2.5 + 17 + 13 + C
    @test Op[7, 7] ≈ 4. + 34 + 12 + C
    @test Op[10, 10] ≈ 1. + 34 - 18 + C

    # Off-diagonal terms.
    @test Op[1, 2] ≈ 20.
    @test Op[1, 4] ≈ 28.
    @test Op[1, 6] ≈ 0.
    @test Op[7, 6] ≈ 4.
    @test Op[6, 7] ≈ 10.
    @test Op[1, 9] ≈ 30.
    @test Op[9, 1] ≈ 24.
    @test Op[1, 10] ≈ 24.

    L = 3
    Sz = 0
    J = reshape(collect(1.:L^2), (L, L))
    W = reshape(collect(1.:L^2), (L, L))
    Jz = collect(1.:L)/4
    Jx = zeros(L)
    Op = build_spin1_many_body_op(L, Sz, J, W, Jz, Jx)
    # Off-diagonal terms.
    @test Op[4, 2] ≈ 12.
    @test Op[4, 5] ≈ 4.

    # Test only Sx terms.
    L = 3
    basis = make_spin1_L_basis(L)
    J = zeros(L, L)
    W = zeros(L, L)
    Jz = zeros(L)
    Jx = collect(1.:L)
    Op = build_spin1_many_body_op(L, basis, J, W, Jz, Jx)

    @test all(Op .≈ transpose(Op))
    @test Op[1, 3] ≈ 1/sqrt(2)
    @test Op[1, 2] ≈ 1/sqrt(2)
    @test Op[12, 15] ≈ 2/sqrt(2)
    @test Op[12, 18] ≈ 2/sqrt(2)
    @test Op[27, 9] ≈ 3/sqrt(2)
end


text:
@defcomp geography begin
    regions = Index()

    area = Variable(index=[time,regions])

    landloss = Parameter(index=[time,regions])
    area0 = Parameter(index=[regions])

    function run_timestep(p, v, d, t);

        if is_first(t)
            for r in d.regions
                v.area[t, r] = p.area0[r]
            end
        else
            for r in d.regions
                v.area[t, r] = v.area[t - 1, r] - p.landloss[t - 1, r]
            end
        end
    end
end


text:
# This script downloads and bundles CA certificates from Amazon Web Services (AWS).
# Currently this makes the following three .pem files available:
# - Amazon root certificate: to use SSL/TLS with various AWS services, including serverless RDS
# - RDS root certificate: to use SSL/TLS with non-serverless RDS (https://docs.aws.amazon.com/AmazonRDS/latest/UserGuide/UsingWithRDS.SSL.html)
# - A combined .pem file of the two certificates above, to make it easier to use SSL/TLS with either serverless or non-serverless RDS

# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "AmazonCACerts"

amazon_root_certificate_filename = "AmazonRootCA1.pem"
rds_root_certificate_version = "2019"
rds_root_certificate_filename = "rds-ca-$rds_root_certificate_version-root.pem"
combined_root_certificates_filename = "combined-amazon-root-ca.pem"

tarball_version = VersionNumber("$rds_root_certificate_version.1")

# Collection of sources required to build AmazonCACerts
sources = [
    FileSource("https://www.amazontrust.com/repository/$amazon_root_certificate_filename",
    "2c43952ee9e000ff2acc4e2ed0897c0a72ad5fa72c3d934e81741cbd54f05bd1",
    filename="$amazon_root_certificate_filename"),

    FileSource("https://s3.amazonaws.com/rds-downloads/$rds_root_certificate_filename",
    "d464378fbb8b981d2b28a1deafffd0113554e6adfb34535134f411bf3c689e73",
    filename="$rds_root_certificate_filename")
]

# Bash recipe for building across all platforms
script = """
cd \$WORKSPACE/srcdir/
mkdir -p \$prefix/share
cp $amazon_root_certificate_filename \$prefix/share/$amazon_root_certificate_filename
cp $rds_root_certificate_filename \$prefix/share/$rds_root_certificate_filename
cat $rds_root_certificate_filename $amazon_root_certificate_filename > \$prefix/share/$combined_root_certificates_filename
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [AnyPlatform()]

# The products that we will ensure are always built
products = [
    FileProduct("share/$amazon_root_certificate_filename", :amazon_ca_root_certificate),
    FileProduct("share/$rds_root_certificate_filename", :rds_ca_root_certificate),
    FileProduct("share/$combined_root_certificates_filename", :combined_ca_root_certificates),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, tarball_version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module Model

# parameter types
export Image, TiledImage, ImageTile,
       SkyPatch, PsfComponent,
       GalaxyComponent, GalaxyPrototype,
       PriorParams, UnconstrainedParams,
       CanonicalParams, BrightnessParams, StarPosParams,
       GalaxyPosParams, GalaxyShapeParams,
       VariationalParams, FreeVariationalParams,
       PsfParams, RawPSF, CatalogEntry,
       init_source

# functions
export align

# constants
export band_letters, D, Ia, B, psf_K,
       galaxy_prototypes, prior,
       shape_standard_alignment,
       brightness_standard_alignment,
       gal_shape_alignment,
       ids_names, ids_free_names,
       ids, ids_free, star_ids, gal_ids, gal_shape_ids, psf_ids, bids

import Base.convert
import Base.+
import Distributions
import FITSIO, WCS
import WCS.WCSTransform
import ForwardDiff
import ..Log

import Base.length

const cfgdir = joinpath(Pkg.dir("Celeste"), "cfg")

const band_letters = ['u', 'g', 'r', 'i', 'z']

# The number of bands (colors).
const B = length(band_letters)

include("model/light_source_model.jl")
include("model/psf_model.jl")
include("model/image_model.jl")
include("model/param_set.jl")
include("model/imaged_sources.jl")
include("model/wcs_utils.jl")

# A vector of variational parameters.  The outer index is
# of celestial objects, and the inner index is over individual
# parameters for that object (referenced using ParamIndex).

typealias VariationalParams{NumType <: Number} Vector{Vector{NumType}}
typealias FreeVariationalParams{NumType <: Number} Vector{Vector{NumType}}

end  # module


text:
using Test, AdvancedHMC
using AdvancedHMC: DualValue, PhasePoint
using LinearAlgebra: dot, diagm

include("common.jl")

@testset "PhasePoint" begin
    for T in [Float32, Float64]
        init_z1() = PhasePoint(
            [T(NaN)],
            [T(NaN)],
            DualValue(zero(T), [zero(T)]),
            DualValue(zero(T), [zero(T)])
        )
        init_z2() = PhasePoint(
            [T(Inf)],
            [T(Inf)],
            DualValue(zero(T), [zero(T)]),
            DualValue(zero(T), [zero(T)])
        )

        @test_logs (:warn, "The current proposal will be rejected due to numerical error(s).") init_z1()
        @test_logs (:warn, "The current proposal will be rejected due to numerical error(s).") init_z2()

        z1 = init_z1()
        z2 = init_z2()

        @test z1.ℓπ.value == z2.ℓπ.value
        @test z1.ℓκ.value == z2.ℓκ.value
    end
end

@testset "Metric" begin
    n_tests = 10

    for T in [Float32, Float64]
        for _ in 1:n_tests
            θ_init = randn(T, D)
            r_init = randn(T, D)

            h = Hamiltonian(UnitEuclideanMetric(T, D), ℓπ, ∂ℓπ∂θ)
            @test -AdvancedHMC.neg_energy(h, r_init, θ_init) == sum(abs2, r_init) / 2
            @test AdvancedHMC.∂H∂r(h, r_init) == r_init

            M⁻¹ = ones(T, D) + abs.(randn(T, D))
            h = Hamiltonian(DiagEuclideanMetric(M⁻¹), ℓπ, ∂ℓπ∂θ)
            @test -AdvancedHMC.neg_energy(h, r_init, θ_init) ≈ r_init' * diagm(0 => M⁻¹) * r_init / 2
            @test AdvancedHMC.∂H∂r(h, r_init) == M⁻¹ .* r_init

            m = randn(T, D, D)
            M⁻¹ = m' * m
            h = Hamiltonian(DenseEuclideanMetric(M⁻¹), ℓπ, ∂ℓπ∂θ)
            @test -AdvancedHMC.neg_energy(h, r_init, θ_init) ≈ r_init' * M⁻¹ * r_init / 2
            @test AdvancedHMC.∂H∂r(h, r_init) == M⁻¹ * r_init
        end
    end
end


text:
# This code is derived from the Modia project and is licensed as follows:
# https://github.com/ModiaSim/Modia.jl/blob/b61daad643ef7edd0c1ccce6bf462c6acfb4ad1a/LICENSE

function try_assign_eq!(ict::IncrementalCycleTracker, vj::Integer, eq::Integer)
    G = ict.graph
    add_edge_checked!(ict, Iterators.filter(!=(vj), 𝑠neighbors(G.graph, eq)), vj) do G
        G.matching[vj] = eq
        G.ne += length(𝑠neighbors(G.graph, eq)) - 1
    end
end

function tearEquations!(ict::IncrementalCycleTracker, Gsolvable, es::Vector{Int},
                        vs::Vector{Int})
    G = ict.graph
    vActive = BitSet(vs)

    for eq in es  # iterate only over equations that are not in eSolvedFixed
        for vj in Gsolvable[eq]
            if G.matching[vj] === unassigned && (vj in vActive)
                r = try_assign_eq!(ict, vj, eq)
                r && break
            end
        end
    end

    return ict
end

function tear_graph_block_modia!(var_eq_matching, graph, solvable_graph, eqs, vars)
    ict = IncrementalCycleTracker(DiCMOBiGraph{true}(graph); dir = :in)
    tearEquations!(ict, solvable_graph.fadjlist, eqs, vars)
    for var in vars
        var_eq_matching[var] = ict.graph.matching[var]
    end
    return nothing
end

function tear_graph_modia(structure::SystemStructure; varfilter = v -> true,
                          eqfilter = eq -> true)
    @unpack graph, solvable_graph = structure
    var_eq_matching = complete(maximal_matching(graph, eqfilter, varfilter))
    var_sccs::Vector{Union{Vector{Int}, Int}} = find_var_sccs(graph, var_eq_matching)

    for vars in var_sccs
        filtered_vars = filter(varfilter, vars)
        ieqs = Int[var_eq_matching[v]
                   for v in filtered_vars if var_eq_matching[v] !== unassigned]
        for var in vars
            var_eq_matching[var] = unassigned
        end
        tear_graph_block_modia!(var_eq_matching, graph, solvable_graph, ieqs, filtered_vars)
    end

    return var_eq_matching
end


text:
function get_sensors_ring(; num_pts=180, φ_inicial=0, φ_final=2π, kR=1, θ=π/2)
    ϕ_range = range(φ_inicial, φ_final; length=num_pts)
    θ_range = ones(length(ϕ_range)) .* θ

    x = kR .* cos.(ϕ_range) .* sin.(θ_range)
    y = kR .* sin.(ϕ_range) .* sin.(θ_range)
    z = kR .* cos.(θ_range)
    return Array(transpose([x y z]))
end


text:
using helloIBox
using Test

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


text:
function main(r::Robot)
    for i in (North, East, South, West)
        path = go_to_border_and_return_path!(r, i; go_around_barriers = true)
        putmarker!(r)
        go_by_path!(r, path)
    end
end

max_stars_repo_name: Savko-fokin/mirea-progs

max_issues_repo_name: Savko-fokin/mirea-progs

max_forks_repo_name: Savko-fokin/mirea-progs


text:
###  only single scale in SVG[^JS]

struct LabelGeometry <: Gadfly.GeometryElement
    # One of :dynamic, :left, :right, :above, :below, :centered
    position::Symbol

    # If true, hide labels that can't be made to not-overlap during dynamic
    # label layout.
    hide_overlaps::Bool

    tag::Symbol
end
LabelGeometry(; position=:dynamic, hide_overlaps=true, tag=empty_tag) =
        LabelGeometry(position, hide_overlaps, tag)

element_aesthetics(::LabelGeometry) = [:x, :y, :label]

default_statistic(::LabelGeometry) = Gadfly.Stat.identity()

"""
    Geom.label[(; position=:dynamic, hide_overlaps=true)]

Place the text strings in the `label` aesthetic at the `x` and `y` coordinates
on the plot frame.  Offset the text according to `position`, which can be
`:left`, `:right`, `:above`, `:below`, `:centered`, or `:dynamic`.  The latter
tries a variety of positions for each label to minimize the number that overlap.
"""
const label = LabelGeometry

# A deferred context function for labeling points in a plot. Optimizing label
# placement depends on knowing the absolute size of the containing context.
function deferred_label_context(geom::LabelGeometry,
                                aes::Gadfly.Aesthetics,
                                theme::Gadfly.Theme,
                                drawctx::ParentDrawContext)

    # Label layout is non-trivial problem. Quite a few papers and at least one
    # Phd thesis has been written on the topic. The approach here is pretty
    # simple. A label may be placed anywhere surrounding a point or hidden.
    # Simulated annealing is used to try to minimize a penalty, which is equal
    # to the number of overlapping or out fo bounds labels, plus terms
    # penalizing hidden labels.
    #
    # TODO:
    # Penalize to prefer certain label positions over others.

    parent_transform = drawctx.t
    units = drawctx.units
    parent_box = drawctx.box

    num_labels = length(aes.label)

    if aes.size == nothing
        padding = fill(theme.point_size, num_labels) .+ theme.label_padding
    else
        padding = aes.size .+ theme.label_padding
    end

    point_positions = Array{AbsoluteVec2}(undef, 0)
    for (x, y) in zip(aes.x, aes.y)
        x = Compose.resolve_position(parent_box, units, parent_transform, Compose.x_measure(x))
        y = Compose.resolve_position(parent_box, units, parent_transform, Compose.y_measure(y))
        push!(point_positions, (x, y))
    end

    label_extents = text_extents(theme.point_label_font, theme.point_label_font_size, aes.label...)

    # the first n values in label_point_{boxes,extents} correspond to the labels
    label_point_boxes = Absolute2DBox[]
    label_point_extents = AbsoluteVec2[]
    for (point_position, pad, (text_width, text_height)) in
            zip(point_positions, padding, label_extents)
        push!(label_point_boxes, Absolute2DBox(
                (point_position[1]-text_width/8, point_position[2]+pad),
                (text_width, text_height)))
        push!(label_point_extents, AbsoluteVec2((text_width, text_height)))
    end

    # the second n values in label_point_{boxes,extents} correspond to the points
    for ((x,y),pad) in zip(point_positions, padding)
        push!(label_point_boxes, Absolute2DBox((x-pad, y-pad), (2pad,2pad)))
        push!(label_point_extents, (2pad,2pad))
    end

    "Return a box containing every point that label `i` could possibly overlap."
    max_extents(i) = Absolute2DBox(
            (point_positions[i][1] - padding[i] - label_point_extents[i][1],
             point_positions[i][2] - padding[i] - label_point_extents[i][2]),
            (2*(padding[i] + label_point_extents[i][1]),
             2*(padding[i] + label_point_extents[i][2])))

    "Returns `true` if boxes `a` and `b` overlap."
    function overlaps(a::Absolute2DBox, b::Absolute2DBox)
        a.x0[1] + a.a[1] >= b.x0[1] && a.x0[1] <= b.x0[1] + b.a[1] &&
        a.x0[2] + a.a[2] >= b.x0[2] && a.x0[2] <= b.x0[2] + b.a[2]
    end

    # Checking for label overlaps is O(n^2). To mitigate these costs, we build a
    # sparse overlap matrix. This also costs O(n^2), but we only have to do it
    # once, rather than every iteration of annealing.
    possible_overlaps = [Array{Int}(undef, 0) for _ in 1:length(label_point_boxes)]

    # TODO: this whole thing would be much more effecient if we forbid from
    # the start labels that overlap points. We should be able to precompute
    # that, since they're static.

    for j in 1:num_labels
        for i in (j+1):num_labels  # automatically skips i==j
            if overlaps(max_extents(i), max_extents(j))
                push!(possible_overlaps[i], j)
                push!(possible_overlaps[j], i)
            end
        end

        for i in (num_labels+1):2*num_labels
            # skip the point box corresponding to label
            i==j+num_labels && continue

            if overlaps(label_point_boxes[i], max_extents(j))
                push!(possible_overlaps[i], j)
                push!(possible_overlaps[j], i)
            end
        end
    end

    "Returns `true` if `a` is fully contained in plot window."
    plot_contains(a::Absolute2DBox) =
            0mm < a.x0[1] && a.x0[1] + a.a[1] < parent_box.a[1] &&
            0mm < a.x0[2] && a.x0[2] + a.a[2] < parent_box.a[2]

    # This variable holds the value of the objective function we wish to
    # minimize. A label overlap is a penalty of 1. Other penaties (out of bounds
    # labels, hidden labels) or calibrated to that.
    total_penalty = 0.0

    for i in 1:num_labels
        if !plot_contains(label_point_boxes[i])
            total_penalty += theme.label_out_of_bounds_penalty
        end
    end

    for j in 1:num_labels
        for i in possible_overlaps[j]
            if i > j && overlaps(label_point_boxes[i], label_point_boxes[j])
                total_penalty += 1
            end
        end
    end

    label_visibility = fill(true, length(label_point_boxes))

    num_iterations = num_labels * theme.label_placement_iterations
    for k in 1:num_iterations
        total_penalty==0 && break

        j = rand(1:num_labels)

        new_total_penalty = total_penalty
        candidate_box = BoundingBox(0mm, 0mm, 0mm, 0mm)
        propose_hide = false

        # Propose flipping the visibility of the label.
        if geom.hide_overlaps && label_visibility[j] && rand() < theme.label_visibility_flip_pr
            propose_hide = true
            new_total_penalty += theme.label_hidden_penalty

        # Propose a change to label placement.
        else
            if !label_visibility[j]
                new_total_penalty -= theme.label_hidden_penalty
            end

            r = rand()
            point_x, point_y = point_positions[j]
            xspan = label_point_extents[j][1]
            yspan = label_point_extents[j][2]

            # TODO: it doesn't know what xspan, yspan are. That's the major
            # source of slowness.
            if rand() < 0.5
                xpos = Gadfly.lerp(rand(), (point_x - 7xspan/8), (point_x - 6xspan/8))
            else
                xpos = Gadfly.lerp(rand(), (point_x - 2xspan/8), (point_x - 1xspan/8))
            end

            ypos = Gadfly.lerp(rand(), (point_y - 3yspan/4), (point_y - 1yspan/4))

            # choose a side
            if r < 0.25 # above
                candidate_box = Absolute2DBox(
                        (xpos, point_y - label_point_extents[j][2] - padding[j]),
                        (label_point_extents[j][1], label_point_extents[j][2]))
            elseif 0.25 <= r < 0.5 # right
                candidate_box = Absolute2DBox(
                        (point_x + padding[j], ypos),
                        (label_point_extents[j][1], label_point_extents[j][2]))
            elseif 0.5 <= r < 0.75 # below
                candidate_box = Absolute2DBox(
                        (xpos, point_y + padding[j]),
                        (label_point_extents[j][1], label_point_extents[j][2]))
            else # left
                candidate_box = Absolute2DBox(
                        (point_x - label_point_extents[j][1] - padding[j], ypos),
                        (label_point_extents[j][1], label_point_extents[j][2]))
            end
        end

        if !plot_contains(label_point_boxes[j]) && label_visibility[j]
            new_total_penalty -= theme.label_out_of_bounds_penalty
        end

        if !propose_hide && !plot_contains(candidate_box)
            new_total_penalty += theme.label_out_of_bounds_penalty
        end

        for i in possible_overlaps[j]
            if label_visibility[i] && label_visibility[j] &&
                        overlaps(label_point_boxes[i], label_point_boxes[j])
                new_total_penalty -= 1
            end

            if label_visibility[i] && overlaps(label_point_boxes[i], candidate_box)
                new_total_penalty += 1
            end
        end

        improvement = total_penalty - new_total_penalty

        T = 0.1 * (1.0 - (k / (1 + num_iterations)))
        if improvement >= 0 || rand() < exp(improvement / T)
            if propose_hide
                label_visibility[j] = false
            else
                label_visibility[j] = true
                label_point_boxes[j] = candidate_box
            end
            total_penalty = new_total_penalty
        end
    end

    return compose!(context(),
        (context(), text([point_positions[i][1] for i in 1:num_labels],
             [point_positions[i][2] for i in 1:num_labels],
             aes.label,
             [hcenter], [vcenter], [Rotation()],
             [(label_point_boxes[i].x0[1] - point_positions[i][1] + label_point_extents[i][1]/2,
               label_point_boxes[i].x0[2] - point_positions[i][2] + label_point_extents[i][2]/2)
              for i in 1:num_labels],
             tag=geom.tag), svgclass("marker")),
        visible(label_visibility),
        font(theme.point_label_font),
        fontsize(theme.point_label_font_size),
        fill(theme.point_label_color),
        stroke(nothing),
        svgclass("geometry"))
end

const label_layouts = Dict(
    :left     => (hright,  vcenter, -1,  0),
    :right    => (hleft,   vcenter,  1,  0),
    :above    => (hcenter, vbottom,  0, -1),
    :below    => (hcenter, vtop,     0,  1),
    :centered => (hcenter, vcenter,  0,  0)
)

function render(geom::LabelGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.Label", aes, :label, :x, :y)

    if geom.position == :dynamic
        return ctxpromise(drawctx -> deferred_label_context(geom, aes, theme, drawctx))
    else
        if !in(geom.position, [:left, :right, :above, :below, :centered])
            error("""
                The position argument of Geom.label must be one of :dynamic,
                :left, :right, :above, :below, :centered
                """)
        end

        hpos, vpos, xoff, yoff = label_layouts[geom.position]
        if aes.size == nothing
            offsets = [(xoff*(theme.point_size + theme.label_padding),
                        yoff*(theme.point_size + theme.label_padding))]
        else
            offsets = [(xoff*(point_size+theme.label_padding),
                        yoff*(point_size+theme.label_padding)) for point_size in aes.size]
        end

        return compose!(context(),
            (context(), text([Compose.x_measure(x) for x in aes.x],
                 [Compose.y_measure(y) for y in aes.y],
                 aes.label,
                 [hpos], [vpos], [Rotation()], offsets, tag=geom.tag),
                 svgclass("marker")),
            font(theme.point_label_font),
            fontsize(theme.point_label_font_size),
            fill(theme.point_label_color),
            stroke(nothing),
            svgclass("geometry"))
    end
end


text:
# function _swap_rows(A::AbstractMatrix{NB}, i::T, nlines::T) where {T<:Int, NB<:Number}
    # for n ∈ (i+1):nlines        # iterate in lines below to check if could be swap
        # if A[n,i] ≠ 0.0         # condition to swap row
            # L = copy(A[i,:])    # copy line to swap
            # A[i,:] = A[n,:]     # swap occur
            # A[n,:] = L
            # break
        # end
    # end
    # return A
# end

"""
    gauss_jordan(A::AbstractMatrix{T}) where T<:Number

Gaussian elimination, also known as row reduction, is an algorithm for solving systems of linear equations.
It consists of a sequence of operations performed on the corresponding matrix of coefficients.
This method can also be used to compute the rank of a matrix, the determinant of a square matrix, and the inverse of an invertible matrix.
https://en.wikipedia.org/wiki/Gaussian_elimination

# Examples/Tests
```julia
julia> M1 = [1 2 3; 4 5 6];
julia> M2 = [1 2 3; 4 8 12];

julia> @test gauss_jordan(M1) == [1 0 -1; 0 1 2]        # Test Passed
julia> @test_throws AssertionError gauss_jordan(M2)     # Test Passed - Thrown: AssertionError
```

# Contributed By:- [AugustoCL](https://github.com/AugustoCL)
"""
function gauss_jordan(A::AbstractMatrix{T}) where {T<:Number}

    # convert to float to avoid InexactError: Int64()
    (T <: Integer) && (A = convert.(Float64, A))

    # check if matrix is singular
    m, n = size(A)
    if m == n
        @assert det(A) ≠ 0.0 "Must insert a non-singular matrix"
    else
        @assert det(A[:,1:end-1]) ≠ 0.0 "Must insert a non-singular matrix or a system matrix [A b]"
    end

    # execute the gauss-jordan elimination
    for i ∈ axes(A, 1)
        (A[i,i] == 0.0) && map([A,i,m]) do (x,y,z)
            for n ∈ (y+1):z                             # check if need swap rows -> this came from function `_swap_rows` (commented above)
                if x[n,y] ≠ 0.0                         # Since it is used only once, an anonymous function with do-block should suffice.
                    L = copy(x[y,:])
                    x[y,:] = x[n,:]
                    x[n,:] = L
                    break
                end
            end
            return x
        end

        @. A[i,:] = A[i,:] ./ A[i,i]                    # divide pivot line by pivot element

        for j ∈ axes(A, 1)                              # iterate each line for each pivot column, except pivot line
            if j ≠ i                                    # jump pivot line
                @. A[j,:] = A[j,:] - A[i,:]*A[j,i]      # apply gauss jordan in each line
            end
        end
    end

    return A
end

max_stars_repo_name: eliascarv/Julia

max_issues_repo_name: eliascarv/Julia

max_forks_repo_name: eliascarv/Julia


text:
"""
Model to approximate cross products of node voltages

```
wdcr[(i,j)] <= wdc[i]*wdc[j]
```
"""
function constraint_voltage_dc(pm::_PM.AbstractWRMModel, n::Int)
    wdc = _PM.var(pm, n, :wdc)
    wdcr = _PM.var(pm, n, :wdcr)

    for (i,j) in _PM.ids(pm, n, :buspairsdc)
        JuMP.@constraint(pm.model, [ wdc[i]/sqrt(2), wdc[j]/sqrt(2), wdcr[(i,j)]/sqrt(2), wdcr[(i,j)]/sqrt(2)] in JuMP.RotatedSecondOrderCone() )
    end
end

"""
Limits dc branch current

```
p[f_idx] <= wdc[f_bus] * Imax
```
"""
function constraint_dc_branch_current(pm::_PM.AbstractWRMModel, n::Int, f_bus, f_idx, ccm_max, p)
    p_dc_fr = _PM.var(pm, n, :p_dcgrid, f_idx)
    wdc_fr = _PM.var(pm, n, :wdc, f_bus)

    JuMP.@constraint(pm.model, p_dc_fr <= wdc_fr * ccm_max * p^2)
end

############## TNEP constraints ###################
"""
Model to approximate cross products of node voltages

```
wdcr[(i,j)] <= wdc[i]*wdc[j]
```
"""
function constraint_voltage_dc_ne(pm::_PM.AbstractWRMModel, n::Int)
    wdc = _PM.var(pm, n, :wdc_ne)
    wdc_frto = _PM.var(pm, n, :wdcr_ne)
    wdc_du_frto = _PM.var(pm, n, :wdcr_du)
    wdc_du_to = _PM.var(pm, n, :wdc_du_to)
    wdc_du_fr = _PM.var(pm, n, :wdc_du_fr)
    z  = _PM.var(pm, n, :branchdc_ne)

        for (l,i,j) in pm.ref[:it][:pm][:nw][n][:arcs_dcgrid_from_ne]
        wdc_to = []
        wdc_fr = []
        wdc_to, wdc_fr = contraint_ohms_dc_branch_busvoltage_structure_W(pm, n, i, j, wdc_du_to, wdc_du_fr)
        JuMP.@constraint(pm.model, [ wdc_du_to[l]/sqrt(2), wdc_du_fr[l]/sqrt(2), wdc_du_frto[l]/sqrt(2), wdc_du_frto[l]/sqrt(2)] in JuMP.RotatedSecondOrderCone() )
        end
end


text:
## GPUArrays interfaces

GPUArrays.device(x::MtlArray) = x.dev

## execution

struct mtlArrayBackend <: AbstractGPUBackend end

struct mtlKernelContext <: AbstractKernelContext end

@inline function GPUArrays.launch_heuristic(::mtlArrayBackend, f::F, args::Vararg{Any,N};
                                             elements::Int, elements_per_thread::Int) where {F,N}
    kernel = @metal launch=false f(mtlKernelContext(), args...)

    # TODO: Be more efficient about creating pipeline states as they're expensive
    # The pipeline state automatically computes occupancy stats, so we just need to parse its fields
    pipeline = MtlComputePipelineState(kernel.fun.lib.device, kernel.fun)
    threads_needed = cld(elements, elements_per_thread)
    # Limit the threadgroup size
    _threads = min(threads_needed, pipeline.maxTotalThreadsPerThreadgroup)
    _blocks  = cld(threads_needed, _threads)
    return (threads=_threads, blocks=_blocks)
end

function GPUArrays.gpu_call(::mtlArrayBackend, f, args, threads::Int, blocks::Int;
                            name::Union{String,Nothing})
    @metal threads=threads grid=blocks name=name f(mtlKernelContext(), args...)
end

## on-device

# indexing
GPUArrays.blockidx(ctx::mtlKernelContext)  = Metal.threadgroup_position_in_grid_1d()
GPUArrays.blockdim(ctx::mtlKernelContext)  = Metal.threads_per_threadgroup_1d()
GPUArrays.threadidx(ctx::mtlKernelContext) = Metal.thread_position_in_threadgroup_1d()
GPUArrays.griddim(ctx::mtlKernelContext)   = Metal.threadgroups_per_grid_1d()

# math

@inline GPUArrays.cos(ctx::mtlKernelContext, x)  = cos(x)
@inline GPUArrays.sin(ctx::mtlKernelContext, x)  = sin(x)
@inline GPUArrays.sqrt(ctx::mtlKernelContext, x) = sqrt(x)
@inline GPUArrays.log(ctx::mtlKernelContext, x)  = log(x)

# memory

# @inline function GPUArrays.LocalMemory(::oneKernelContext, ::Type{T}, ::Val{dims}, ::Val{id}
#                                       ) where {T, dims, id}
#     ptr = oneAPI.emit_localmemory(Val(id), T, Val(prod(dims)))
#     oneDeviceArray(dims, LLVMPtr{T, onePI.AS.Local}(ptr))
# end

# synchronization

@inline GPUArrays.synchronize_threads(::mtlKernelContext) = Metal.threadgroup_barrier()

#
# Host abstractions
#

GPUArrays.backend(::Type{<:MtlArray}) = mtlArrayBackend()

const GLOBAL_RNGs = Dict{MtlDevice,GPUArrays.RNG}()
function GPUArrays.default_rng(::Type{<:MtlArray})
    dev = MtlDevice(1)
    get!(GLOBAL_RNGs, dev) do
        N = 128 # Size of default oneAPI working group with barrier, so should be good for Metal
        state = MtlArray{NTuple{4, UInt32}}(undef, N)
        rng = GPUArrays.RNG(state)
        Random.seed!(rng)
        rng
    end
end


text:
using Bonsai, JSON3, StructTypes
using Test
using StructTypes: @Struct
using Bonsai: open_api_parameters, ParameterObject,
	ResponseObject,  RequestBodyObject,
	handler_writes, HttpParameter, handler_reads
using CodeInfoTools: code_inferred
using Bonsai: PathItemObject, MediaTypeObject, ParameterObject, OperationObject, OpenAPI
# https://raw.githubusercontent.com/OAI/OpenAPI-Specification/main/examples/v3.0/petstore.json
# https://blog.stoplight.io/openapi-json-schema#:~:text=You%20can%20use%20JSON%20Schema,generate%20an%20entire%20mock%20server.

@Struct struct Id
	id::String
end

@Struct struct Pet1
	id::Int64
	name::String
	tag::String
end

@Struct struct Limit1
	limit::Int
	offset::Int
end

@Struct struct Offset
	start::Int
	n::Int
end

@Struct struct AuthHeaders
	x_pass::String
	x_user::String
end

@testset "Query" begin
	@test length(open_api_parameters(Query{Limit1})) == 2
	@test length(open_api_parameters(Query{Offset})) == 2
end

@testset "Body" begin
	b1 = Body(Pet)
	@test RequestBodyObject(typeof(b1)) isa RequestBodyObject
end

@testset "OpenAPI" begin
	app = App()

	app.post("/pets/") do stream
		body = Bonsai.read(stream, Body(Id))
	end

	app.get("/pets/") do stream
		pet = Pet1(body.id,"bob", "cat")
		Bonsai.write(stream , pet)
	end

	app.get("/pets/{id}") do stream
		pets = [ Pet1(body.id, "bob", "cat") for i in 1:10]
		Bonsai.write(pets)
	end

	app.get("/pets/{id}") do stream, next
		pets = [ Pet1(body.id, "bob", "cat") for i in 1:10]
		Bonsai.write(pets)
	end

	get_pets = match(app.paths,  Params(), "GET", ["pets"], 1)
	create_pets = match(app.paths,  Params(), "POST", ["pets"], 1)

	get_pets = Bonsai.match_middleware(app.paths,  Params(), "GET", ["pets"], 1)

	@test Bonsai.RequestBodyObject(
		Bonsai.handler_reads(create_pets.fn)[1]
	) isa Bonsai.RequestBodyObject

	@test OpenAPI(app) isa OpenAPI
end

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
using PATHSolver
using Test

using MathOptInterface
const MOI = MathOptInterface

@testset "MOI.Name" begin
    model = PATHSolver.Optimizer()
    @test MOI.get(model, MOI.SolverName()) == PATHSolver.c_api_Path_Version()
end

@testset "MOI.AbstractOptimizer" begin
    @test PATHSolver.Optimizer() isa MOI.AbstractOptimizer
end

@testset "RawParameter" begin
    model = PATHSolver.Optimizer()
    @test MOI.get(model, MOI.RawParameter("output")) === nothing
    MOI.set(model, MOI.RawParameter("output"), "no")
    @test MOI.get(model, MOI.RawParameter("output")) == "no"
end

@testset "Invalid models" begin
    @testset "Infeasible" begin
        model = PATHSolver.Optimizer()
        x = MOI.add_variable(model)
        MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Interval(0.0, -1.0))
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.INVALID_MODEL
        @test MOI.get(model, MOI.RawStatusString()) == "Problem has a bound error"
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.UNKNOWN_RESULT_STATUS
    end
    @testset "Binary" begin
        model = PATHSolver.Optimizer()
        x = MOI.add_variable(model)
        @test_throws(
            MOI.UnsupportedConstraint,
            MOI.add_constraint(model, MOI.SingleVariable(x), MOI.ZeroOne())
        )
    end
    @testset "wrong dimension" begin
        model = PATHSolver.Optimizer()
        x = MOI.add_variable(model)
        MOI.add_constraint(
            model,
            MOI.VectorAffineFunction(
                MOI.VectorAffineTerm{Float64}[],
                [0.0]
            ),
            MOI.Complements(1)
        )
        @test_throws(
            ErrorException(
                "Dimension of constant vector 1 does not match the required dimension of " *
                "the complementarity set 2."
            ),
            MOI.optimize!(model)
        )
    end
    @testset "non-zero variable offset" begin
        model = PATHSolver.Optimizer()
        x = MOI.add_variable(model)
        MOI.add_constraint(
            model,
            MOI.VectorAffineFunction(
                [
                    MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x)),
                    MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x))
                ],
                [0.0, 1.0]
            ),
            MOI.Complements(1)
        )
        @test_throws(
            ErrorException(
                "VectorAffineFunction malformed: a constant associated with a " *
                "complemented variable is not zero: [1.0]."
            ),
            MOI.optimize!(model)
        )
    end
    @testset "output dimension too large" begin
        model = PATHSolver.Optimizer()
        x = MOI.add_variable(model)
        MOI.add_constraint(
            model,
            MOI.VectorAffineFunction(
                [
                    MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x)),
                    MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(1.0, x))
                ],
                [0.0, 0.0]
            ),
            MOI.Complements(1)
        )
        @test_throws(
            ErrorException("VectorAffineFunction malformed: output_index 3 is too large."),
            MOI.optimize!(model)
        )
    end
    @testset "missing complement variable" begin
        model = PATHSolver.Optimizer()
        x = MOI.add_variable(model)
        MOI.add_constraint(
            model,
            MOI.VectorAffineFunction(
                [
                    MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))
                ],
                [0.0, 0.0]
            ),
            MOI.Complements(1)
        )
        @test_throws(
            ErrorException("VectorAffineFunction malformed: expected variable in row 2."),
            MOI.optimize!(model)
        )
    end
    @testset "output dimension too large" begin
        model = PATHSolver.Optimizer()
        x = MOI.add_variable(model)
        MOI.add_constraint(
            model,
            MOI.VectorAffineFunction(
                [
                    MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x)),
                    MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(2.0, x))
                ],
                [0.0, 0.0]
            ),
            MOI.Complements(1)
        )
        @test_throws(
            ErrorException(
                "VectorAffineFunction malformed: variable $(x) has a coefficient that is " *
                "not 1 in row 2 of the VectorAffineFunction."
            ),
            MOI.optimize!(model)
        )
    end
    @testset "variable in multiple constraints" begin
        model = PATHSolver.Optimizer()
        x = MOI.add_variable(model)
        MOI.add_constraint(
            model,
            MOI.VectorAffineFunction(
                [
                    MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x)),
                    MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x)),
                    MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(1.0, x)),
                    MOI.VectorAffineTerm(4, MOI.ScalarAffineTerm(1.0, x))
                ],
                [0.0, 0.0, 0.0, 0.0]
            ),
            MOI.Complements(2)
        )
        @test_throws(
            ErrorException(
                "The variable $(x) appears in more than one complementarity constraint."
            ),
            MOI.optimize!(model)
        )
    end
end

@testset "Example 1" begin
    model = PATHSolver.Optimizer()
    MOI.set(model, MOI.RawParameter("time_limit"), 60)
    @test MOI.supports(model, MOI.Silent()) == true
    @test MOI.get(model, MOI.Silent()) == false
    MOI.set(model, MOI.Silent(), true)
    @test MOI.get(model, MOI.Silent()) == true
    x = MOI.add_variables(model, 4)
    @test MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex)
    @test MOI.get(model, MOI.VariablePrimalStart(), x[1]) === nothing
    MOI.add_constraint.(model, MOI.SingleVariable.(x), MOI.Interval(0.0, 10.0))
    MOI.set.(model, MOI.VariablePrimalStart(), x, 0.0)
    M = Float64[
        0  0 -1 -1;
        0  0  1 -2;
        1 -1  2 -2;
        1  2 -2  4
    ]
    q = [2; 2; -2; -6]
    for i in 1:4
        terms = [
            MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x[i]))
        ]
        for j in 1:4
            iszero(M[i, j]) && continue
            push!(
                terms,
                MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(M[i, j], x[j]))
            )
        end
        MOI.add_constraint(
            model,
            MOI.VectorAffineFunction(terms, [q[i], 0.0]),
            MOI.Complements(1)
        )
    end
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.NO_SOLUTION
    @test MOI.get(model, MOI.RawStatusString()) == "MOI.optimize! was not called yet"
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.LOCALLY_SOLVED
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.RawStatusString()) == "The problem was solved"
    x_val = MOI.get.(model, MOI.VariablePrimal(), x)
    @test isapprox(x_val, [2.8, 0.0, 0.8, 1.2])
end

# An implementation of the GAMS model transmcp.gms: transporation model
# as an equilibrium problem.
#
# https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_transmcp.html
#
# Transportation model as equilibrium problem (TRANSMCP,SEQ=126)
#
#    Dantzig's original transportation model (TRNSPORT) is
#    reformulated as a linear complementarity problem.  We first
#    solve the model with fixed demand and supply quantities, and
#    then we incorporate price-responsiveness on both sides of the
#    market.
#
# Dantzig, G B, Chapter 3.3. In Linear Programming and Extensions.
# Princeton University Press, Princeton, New Jersey, 1963.
@testset "transmcp" begin
    plants = ["seattle", "san-diego"]
    P = length(plants)
    capacity = Dict("seattle" => 350, "san-diego" => 600)
    markets = ["new-york", "chicago", "topeka"]
    M = length(markets)
    demand = Dict("new-york" => 325, "chicago"  => 300, "topeka" => 275)
    distance = Dict(
        ("seattle" => "new-york")   => 2.5,
        ("seattle" => "chicago")    => 1.7,
        ("seattle" => "topeka")     => 1.8,
        ("san-diego" => "new-york") => 2.5,
        ("san-diego" => "chicago")  => 1.8,
        ("san-diego" => "topeka")   => 1.4
    )

    model = PATHSolver.Optimizer()
    MOI.set(model, MOI.Silent(), true)

    # w[i in plants] >= 0
    w = MOI.add_variables(model, P)
    MOI.add_constraint.(model, MOI.SingleVariable.(w), MOI.GreaterThan(0.0))

    # p[j in markets] >= 0
    p = MOI.add_variables(model, M)
    MOI.add_constraint.(model, MOI.SingleVariable.(p), MOI.GreaterThan(0.0))

    # x[i in plants, j in markets] >= 0
    x = reshape(MOI.add_variables(model, P * M), P, M)
    MOI.add_constraint.(model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0))

    # w[i] + 90 * distance[i => j] / 1000 - p[j] ⟂ x[i, j]
    for (i, plant) in enumerate(plants), (j, market) in enumerate(markets)
        MOI.add_constraint(
            model,
            MOI.VectorAffineFunction(
                [
                    MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, w[i])),
                    MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(-1.0, p[j])),
                    MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x[i, j]))
                ],
                [90 * distance[plant => market] / 1000, 0.0]
            ),
            MOI.Complements(1)
        )
    end

    # capacity[i] - sum(x[i, :]) ⟂ w[i]
    terms = MOI.VectorAffineTerm{Float64}[]
    for (i, plant) in enumerate(plants)
        for j = 1:M
            push!(terms, MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(-1.0, x[i, j])))
        end
        push!(terms, MOI.VectorAffineTerm(P + i, MOI.ScalarAffineTerm(1.0, w[i])))
    end
    q = vcat([capacity[p] for p in plants], zeros(P))
    MOI.add_constraint(
        model, MOI.VectorAffineFunction(terms, q), MOI.Complements(P)
    )

    # sum(x[:, j]) - demand[j] ⟂ p[j]
    terms = MOI.VectorAffineTerm{Float64}[]
    for (j, market) in enumerate(markets)
        for i = 1:P
            push!(terms, MOI.VectorAffineTerm(j, MOI.ScalarAffineTerm(1.0, x[i, j])))
        end
        push!(terms, MOI.VectorAffineTerm(M + j, MOI.ScalarAffineTerm(1.0, p[j])))
    end
    q = vcat([-demand[m] for m in markets], zeros(M))
    MOI.add_constraint(
        model, MOI.VectorAffineFunction(terms, q), MOI.Complements(M)
    )

    MOI.optimize!(model)

    @test isapprox(
        MOI.get.(model, MOI.VariablePrimal(), p),
        [0.225, 0.153, 0.126],
        atol = 1e-3
    )
end


text:
using LinearAlgebra
function MSE(y, t; mean_f = true) # this is Mean Square Error
    if size(y) != size(t)
        y = y'
    end
    mse = sum((y - t) .^ 2)
    if mean_f
        mse /= length(y) #mse is Mean Square error
    end
    return mse
end
function CEE(p, t; mean_f = true, t_f = false) #this is Cross entropy Error
    function safe_log(x, miniv = 0.00000000001)
        return map(log, clamp.(x, miniv, Inf))
    end
    if t_f
        loss = -sum(t .* safe_log(p) + (1 .- t) .* safe_log(1 .- p))
    else
        loss = -sum(t .* safe_log(p))
    end
    if mean_f
        loss = loss / length(t)
    end
    return loss
end


text:
function facs_binning_snr(filepath; debug=false, quiet=false)
    if !debug
        parameters = Dict{Symbol, Vector}(
            :representation => [10, 100, 1000],
            :bottleneck_representation => [10, 100, 1000],
            :seq_depth => [10, 100, 1000],
            :σ => [1.0, 1.0, 0.5],
            :bin_info => [OrderedDict(:bin1 => (0.0, p), :bin2 => (1.0-p, 1.0)) for p in linspace(0.5, 0.025, 30)]
        )
        num_runs = 100
    else
        parameters = Dict{Symbol, Vector}(
            :representation => [100],
            :bottleneck_representation => [100],
            :seq_depth => [100],
            :σ => [1.0],
            :bin_info => [OrderedDict(:bin1 => (0.0, 0.25), :bin2 => (0.75, 1.0))]
        )
        num_runs = 1
    end

    runs = grouped_param_space(FacsScreen(), parameters, [:bin_info], num_runs)

    before = time()
    results = pmap(args -> run_exp(args[1], Library(CRISPRi()), compute_snr; run_idx=args[2]), runs)
    (!quiet) && println("$(time() - before) seconds")
    results = DataFrame(permutedims(hcat(results...), [2, 1]))
    results[:crisprtype] = "CRISPRi"
    before = time()
    results2 = pmap(args -> run_exp(args[1], Library(CRISPRn()), compute_snr; run_idx=args[2]), runs)
    (!quiet) && println("$(time() - before) seconds")
    results2 = DataFrame(permutedims(hcat(results2...), [2, 1]))
    results2[:crisprtype] = "CRISPRn"
    results = vcat(results, results2)

    hierarchy = reshape([:snr, :signal, :noise], (3, 1))
    new_names = [[:technique, :score]...; fieldnames(FacsScreen)...; :run_idx; :crisprtype;]

    results = construct_hierarchical_label(hierarchy, results, new_names)
    results[:bin_info] = Float64[el[:bin1][2] for el in results[:bin_info]]
    CSV.write(filepath, results)
end


text:
module day17

include("inputs.jl")

ex = inputs.ex_d17()
input = inputs.input_d17()

function parse_input(input, dim)
    input = split(input, "\n") .|> String
    a = [char=='#' for line in input for char in line]
    @assert rem(length(a), sqrt(length(a)))==0.0 # square-able dim
    side = Int(sqrt(length(a)))
    input2d = reshape(a, (side, side)) |> permutedims
    if dim == 2
        return input2d
    elseif dim == 3
        return reshape(input2d, size(input2d)..., 1)
    elseif dim == 4
        return reshape(input2d, size(input2d)..., 1, 1)
    end
end

function count_actives(input, dim)
    cycles = 6
    parsed_input = parse_input(input, dim)

    grid_size = size(parsed_input) .+ (2*cycles)
    grid = zeros(Bool, grid_size)

    seed_loc = UnitRange.(cycles+1, cycles .+ size(parsed_input))
    grid[seed_loc...] = parsed_input

    adj = CartesianIndices((repeat([-1:1], dim)...,))

    next_grid = copy(grid)
    for cycle in 1:cycles
        for coord ∈ CartesianIndices(grid)
            n_adj = sum(get(grid, coord + a, 0) for a in adj) - grid[coord]
            if grid[coord] && n_adj ∈ [2, 3]
                # println("keep ", coord, "true")
                next_grid[coord] = true
            elseif grid[coord] && n_adj ∉ [2, 3]
                 next_grid[coord] = false
            elseif !grid[coord] && n_adj == 3
                # println("set ", coord, "from false to true")
                next_grid[coord] = true
            elseif !grid[coord] && n_adj != 3
                 next_grid[coord] = false
            end
        end
        grid .= next_grid
    end
    return sum(grid)
end

part1(input) = count_actives(input, 3)

@assert part1(ex) == 112
@info "Day 17, Part 1 answer: $(part1(input))"

part2(input) = count_actives(input, 4)
@assert part2(ex) == 848
@info "Day 17, Part 2 answer: $(part2(input))"

end


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

@doc raw"""JSON represents any valid JSON value. These types are supported: bool, int64, float64, string, []interface{}, map[string]interface{} and nil.

    IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON(;
    )

"""
mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON <: SwaggerModel

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON(;)
        o = new()
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON = Dict{Symbol,Symbol}()
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON = Dict{Symbol,String}()
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSON)
    true
end

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


text:
# test_with_actors2.jl

# Basically the same as with test_local_channels.jl

using Pkg
Pkg.activate(".")
Pkg.precompile()

using Rocket
using DataFrames

# start docker containers
cmd = `docker start test_sshd`
run(cmd)

cmd = `docker start test_sshd2`
run(cmd)

cmd = `docker ps`
run(cmd)

sleep(5)

# enable distrbuted computing
using Distributed

#addprocs(4; exeflags=`--project=$(Base.active_project())`)
addprocs([("rob@172.17.0.2", 1), ("rob@172.17.0.3", 1)]; exeflags=`--project=$(Base.active_project())`, tunnel=true, dir="/home/rob")

#addprocs([("pi@192.168.2.3", 1)]; exename=`/home/pi/julia/julia-1.3.1/bin/julia`, dir="/home/pi")

# remove processes > 3
while length(procs()) ≥ 4
    rmprocs(procs()[length(procs())])
end

gl_pid = procs()[2] # general ledger
ar_pid = procs()[3] # accounts receivable (orders/bankstatements)

# activate the packages
@everywhere begin
    using AppliSales
    using AppliGeneralLedger
    using AppliAR
    using Query
end;

struct StmActor <: Actor{String}
    ar_pid::Int64
    StmActor(ar_pid::Int) = new(ar_pid)
end
Rocket.on_next!(actor::StmActor, data::String) = begin
    if data == "READ_STMS"
        stms = AppliAR.read_bank_statements("./bank.csv")
        unpaid_inv = @fetchfrom actor.ar_pid retrieve_unpaid_invoices()
        entries = @fetchfrom actor.ar_pid AppliAR.process(unpaid_inv, stms)
        subscribe!(from(entries), gl_actor)
    end
end
Rocket.on_complete!(actor::StmActor) = @info("StmActor completed!")
Rocket.on_error!(actor::StmActor, err) = @info(error(err))

struct SalesActor <: Actor{String} end
Rocket.on_next!(actor::SalesActor, data::String) = begin
    if data == "START"
        #ar_actor = ARActor()
        orders = @fetch AppliSales.process()
        subscribe!(from(orders), ar_actor)
    end
end
Rocket.on_complete!(actor::SalesActor) = @info("SalesActor completed!")
Rocket.on_error!(actor::SalesActor, err) = @info(error(err))

struct ARActor <: Actor{AppliSales.Order}
    values::Vector{AppliGeneralLedger.JournalEntry}
    ar_pid::Int64
    gl_pid::Int64
    ARActor(ar_pid, gl_pid) = new(Vector{AppliGeneralLedger.JournalEntry}(), ar_pid, gl_pid)
end
Rocket.on_next!(actor::ARActor, data::AppliSales.Order) = begin
        d = @fetchfrom actor.ar_pid AppliAR.process([data])
        push!(actor.values, d[1])
end
Rocket.on_complete!(actor::ARActor) = begin
    #println(actor.values)
    @info(typeof(actor.values))
    subscribe!(from(actor.values), gl_actor)
    @info("ARActor Completed!")
end
Rocket.on_error!(actor::ARActor, err) = @info(error(err))

struct GLActor <: Actor{Any}
    gl_pid::Int64
    GLActor(gl_pid) = new(gl_pid)
end
Rocket.on_next!(actor::GLActor, data::Any) = begin
    if data isa AppliGeneralLedger.JournalEntry
        result = @fetchfrom actor.gl_pid AppliGeneralLedger.process([data])
    end
end
Rocket.on_complete!(actor::GLActor) = @info("GLActor completed!")
Rocket.on_error!(actor::GLActor, err) = @info(error(err))

sales_actor = SalesActor()
ar_actor = ARActor(ar_pid, gl_pid)
gl_actor = GLActor(gl_pid)
stm_actor = StmActor(ar_pid)

subscribe!(from(["START"]), sales_actor)

subscribe!(from(["READ_STMS"]), stm_actor)

# print aging report
r1 = @fetchfrom ar_pid report()
result = DataFrame(r1)
println("\nUnpaid invoices\n===============")
show(result)

# print general ledger
r2 = @fetchfrom gl_pid AppliGeneralLedger.read_from_file("./test_ledger.txt")
df = DataFrame(r2)
println("\nGeneral Ledger mutations\n========================")
show(df)

df2 = r2 |> @filter(_.accountid == 1300) |> DataFrame
println("\nAccounts receivable\n===================")
show(df2)

println("\nBalance Accounts Receivable: $(sum(df2.debit) - sum(df2.credit))")

# open shell in container
cmd = `ssh rob@172.17.0.2`
@info("after run(cmd) is activated: goto console, press Enter, and rm test3_* files. Leave the container with Ctrl-D")
run(cmd)

# open shell in container
cmd = `ssh rob@172.17.0.3`
@info("after run(cmd) is activated: goto console, press Enter, and rm test3_invoicing.sqlite*. Leave the container with Ctrl-D")
run(cmd)
@info("Ctrl-L to clean the consule. Close julia with Ctrl-D.")

# end