text:
# IJulia redirects stdout and stderr into "stream" messages sent to the
# Jupyter front-end.

# create a wrapper type around redirected stdio streams,
# both for overloading things like `flush` and so that we
# can set properties like `color`.
struct IJuliaStdio{IO_t <: IO} <: Base.AbstractPipe
    io::IOContext{IO_t}
end
IJuliaStdio(io::IO, stream::AbstractString="unknown") =
    IJuliaStdio{typeof(io)}(IOContext(io, :color=>Base.have_color,
                            :jupyter_stream=>stream,
                            :displaysize=>displaysize()))
Base.pipe_reader(io::IJuliaStdio) = io.io.io
Base.pipe_writer(io::IJuliaStdio) = io.io.io
Base.lock(io::IJuliaStdio) = lock(io.io.io)
Base.unlock(io::IJuliaStdio) = unlock(io.io.io)
Base.in(key_value::Pair, io::IJuliaStdio) = in(key_value, io.io)
Base.haskey(io::IJuliaStdio, key) = haskey(io.io, key)
Base.getindex(io::IJuliaStdio, key) = getindex(io.io, key)
Base.get(io::IJuliaStdio, key, default) = get(io.io, key, default)
Base.displaysize(io::IJuliaStdio) = displaysize(io.io)
Base.unwrapcontext(io::IJuliaStdio) = Base.unwrapcontext(io.io)
Base.setup_stdio(io::IJuliaStdio, readable::Bool) = Base.setup_stdio(io.io.io, readable)

for s in ("stdout", "stderr", "stdin")
    f = Symbol("redirect_", s)
    Sq = QuoteNode(Symbol(uppercase(s)))
    sq = QuoteNode(Symbol(s))
    @eval function Base.$f(io::IJuliaStdio)
        io[:jupyter_stream] != $s && throw(ArgumentError(string("expecting ", $s, " stream")))
        if isdefined(Base, :stdout)
            Core.eval(Base, Expr(:(=), $sq, io))
        else
            # On Julia 0.6-, the variables are called Base.STDIO, not Base.stdio
            Core.eval(Base, Expr(:(=), $Sq, io))
            # We also need to change Compat.stdxxx for packages that use it
            # ref: https://github.com/JuliaLang/IJulia.jl/issues/690
            Core.eval(Compat, Expr(:(=), $sq, io))
        end
        return io
    end
end

# logging in verbose mode goes to original stdio streams.  Use macros
# so that we do not even evaluate the arguments in no-verbose modes

using Printf
function get_log_preface()
    t = now()
    taskname = get(task_local_storage(), :IJulia_task, "")
    @sprintf("%02d:%02d:%02d(%s): ", Dates.hour(t),Dates.minute(t),Dates.second(t),taskname)
end

macro vprintln(x...)
    quote
        if verbose::Bool
            println(orig_stdout[], get_log_preface(), $(map(esc, x)...))
        end
    end
end

macro verror_show(e, bt)
    quote
        if verbose::Bool
            showerror(orig_stderr[], $(esc(e)), $(esc(bt)))
        end
    end
end

#name=>iobuffer for each stream ("stdout","stderr") so they can be sent in flush
const bufs = Dict{String,IOBuffer}()
const stream_interval = 0.1
# maximum number of bytes in libuv/os buffer before emptying
const max_bytes = 10*1024
# max output per code cell is 512 kb by default
const max_output_per_request = Ref(1 << 19)

"""
Continually read from (size limited) Libuv/OS buffer into an `IObuffer` to avoid problems when
the Libuv/OS buffer gets full (https://github.com/JuliaLang/julia/issues/8789). Send data immediately
when buffer contains more than `max_bytes` bytes. Otherwise, if data is available it will be sent every
`stream_interval` seconds (see the Timers set up in watch_stdio). Truncate the output to `max_output_per_request`
bytes per execution request since excessive output can bring browsers to a grinding halt.
"""
function watch_stream(rd::IO, name::AbstractString)
    task_local_storage(:IJulia_task, "read $name task")
    try
        buf = IOBuffer()
        bufs[name] = buf
        while !eof(rd) # blocks until something is available
            nb = bytesavailable(rd)
            if nb > 0
                stdio_bytes[] += nb
                # if this stream has surpassed the maximum output limit then ignore future bytes
                if stdio_bytes[] >= max_output_per_request[]
                    read(rd, nb) # read from libuv/os buffer and discard
                    if stdio_bytes[] - nb < max_output_per_request[]
                        send_ipython(publish[], msg_pub(execute_msg, "stream",
                                     Dict("name" => "stderr", "text" => "Excessive output truncated after $(stdio_bytes[]) bytes.")))
                    end
                else
                    write(buf, read(rd, nb))
                end
            end
            if buf.size > 0
                if buf.size >= max_bytes
                    #send immediately
                    send_stream(name)
                end
            end
        end
    catch e
        # the IPython manager may send us a SIGINT if the user
        # chooses to interrupt the kernel; don't crash on this
        if isa(e, InterruptException)
            watch_stream(rd, name)
        else
            rethrow()
        end
    end
end

function send_stdio(name)
    if verbose::Bool && !haskey(task_local_storage(), :IJulia_task)
        task_local_storage(:IJulia_task, "send $name task")
    end
    send_stream(name)
end

send_stdout(t::Timer) = send_stdio("stdout")
send_stderr(t::Timer) = send_stdio("stderr")

"""
Jupyter associates cells with message headers. Once a cell's execution state has
been set as to idle, it will silently drop stream messages (i.e. output to
stdout and stderr) - see https://github.com/jupyter/notebook/issues/518.
When using Interact, and a widget's state changes, a new
message header is sent to the IJulia kernel, and while Reactive
is updating Signal graph state, it's execution state is busy, meaning Jupyter
will not drop stream messages if Interact can set the header message under which
the stream messages will be sent. Hence the need for this function.
"""
function set_cur_msg(msg)
    global execute_msg = msg
end

function send_stream(name::AbstractString)
    buf = bufs[name]
    if buf.size > 0
        d = take!(buf)
        n = num_utf8_trailing(d)
        dextra = d[end-(n-1):end]
        resize!(d, length(d) - n)
        s = String(copy(d))
        if isvalid(String, s)
            write(buf, dextra) # assume that the rest of the string will be written later
            length(d) == 0 && return
        else
            # fallback: base64-encode non-UTF8 binary data
            sbuf = IOBuffer()
            print(sbuf, "base64 binary data: ")
            b64 = Base64EncodePipe(sbuf)
            write(b64, d)
            write(b64, dextra)
            close(b64)
            print(sbuf, '\n')
            s = String(take!(sbuf))
        end
        send_ipython(publish[],
             msg_pub(execute_msg, "stream",
                     Dict("name" => name, "text" => s)))
    end
end

"""
If `d` ends with an incomplete UTF8-encoded character, return the number of trailing incomplete bytes.
Otherwise, return `0`.
"""
function num_utf8_trailing(d::Vector{UInt8})
    i = length(d)
    # find last non-continuation byte in d:
    while i >= 1 && ((d[i] & 0xc0) == 0x80)
        i -= 1
    end
    i < 1 && return 0
    c = d[i]
    # compute number of expected UTF-8 bytes starting at i:
    n = c <= 0x7f ? 1 : c < 0xe0 ? 2 : c < 0xf0 ? 3 : 4
    nend = length(d) + 1 - i # num bytes from i to end
    return nend == n ? 0 : nend
end

"""
    readprompt(prompt::AbstractString; password::Bool=false)

Display the `prompt` string, request user input,
and return the string entered by the user.  If `password`
is `true`, the user's input is not displayed during typing.
"""
function readprompt(prompt::AbstractString; password::Bool=false)
    if !execute_msg.content["allow_stdin"]
        error("IJulia: this front-end does not implement stdin")
    end
    send_ipython(raw_input[],
                 msg_reply(execute_msg, "input_request",
                           Dict("prompt"=>prompt, "password"=>password)))
    while true
        msg = recv_ipython(raw_input[])
        if msg.header["msg_type"] == "input_reply"
            return msg.content["value"]
        else
            error("IJulia error: unknown stdin reply")
        end
    end
end

# override prompts using julia#28038 in 0.7
if isdefined(Base, :prompt) && isdefined(Base, :getpass)
    function check_prompt_streams(input::IJuliaStdio, output::IJuliaStdio)
        if get(input,:jupyter_stream,"unknown") != "stdin" ||
            get(output,:jupyter_stream,"unknown") != "stdout"
            throw(ArgumentError("prompt on IJulia stdio streams only works for stdin/stdout"))
         end
     end
    function Base.prompt(input::IJuliaStdio, output::IJuliaStdio, message::AbstractString; default::AbstractString="")
        check_prompt_streams(input, output)
        val = chomp(readprompt(message * ": "))
        return isempty(val) ? default : val
    end
    function Base.getpass(input::IJuliaStdio, output::IJuliaStdio, message::AbstractString)
        check_prompt_streams(input, output)
        # fixme: should we do more to zero memory associated with the password?
        #        doing this properly might require working with the raw ZMQ message buffer here
        return Base.SecretBuffer!(Vector{UInt8}(codeunits(readprompt(message * ": ", password=true))))
    end
end

# IJulia issue #42: there doesn't seem to be a good way to make a task
# that blocks until there is a read request from STDIN ... this makes
# it very hard to properly redirect all reads from STDIN to pyin messages.
# In the meantime, however, we can just hack it so that readline works:
import Base.readline
function readline(io::IJuliaStdio)
    if get(io,:jupyter_stream,"unknown") == "stdin"
        return readprompt("stdin> ")
    else
        readline(io.io)
    end
end

function watch_stdio()
    task_local_storage(:IJulia_task, "init task")
    if capture_stdout
        read_task = @async watch_stream(read_stdout[], "stdout")
        #send stdout stream msgs every stream_interval secs (if there is output to send)
        Timer(send_stdout, stream_interval, interval=stream_interval)
    end
    if capture_stderr
        readerr_task = @async watch_stream(read_stderr[], "stderr")
        #send STDERR stream msgs every stream_interval secs (if there is output to send)
        Timer(send_stderr, stream_interval, interval=stream_interval)
    end
end

function flush_all()
    flush_cstdio() # flush writes to stdout/stderr by external C code
    flush(stdout)
    flush(stderr)
end

function oslibuv_flush()
    #refs: https://github.com/JuliaLang/IJulia.jl/issues/347#issuecomment-144505862
    #      https://github.com/JuliaLang/IJulia.jl/issues/347#issuecomment-144605024
    @static if Sys.iswindows()
        ccall(:SwitchToThread, stdcall, Cvoid, ())
    end
    yield()
    yield()
end

import Base.flush
function flush(io::IJuliaStdio)
    flush(io.io)
    oslibuv_flush()
    send_stream(get(io,:jupyter_stream,"unknown"))
end


text:
import LinearAlgebra
import LinearAlgebra: UpperTriangular, istril
import Base.GMP.MPZ

# Algorithm 3 from George Havas, Bohdan S. Majewski, and Keith R. Matthews,
# "Extended GCD and Hermite Normal Form Algorithms via Lattice Basis Reduction"

function gcd_reduce1!(k, i, m, λ, D, a, B, tmpa, tmpD, tmpB, tmpλ)
    q = @inbounds begin
        if !iszero(a[i])
            div(a[k], a[i], RoundNearest)
        elseif 2*λ[i,k] > D[i+1]
            div(λ[i,k], D[i+1], RoundNearest)
        else
            zero(BigInt)
        end
    end

    @inbounds if !iszero(q)
        MPZ.mul!(tmpa, q, a[i])
        MPZ.sub!(a[k], tmpa)
        @simd for j in 1:m
            MPZ.mul!(tmpB[j], q, B[j,i])
            MPZ.sub!(B[j,k], tmpB[j])
        end
        MPZ.mul!(tmpD, q, D[i+1])
        MPZ.sub!(λ[i,k], tmpD)
        @simd for j in 1:i-1
            MPZ.mul!(tmpλ[j], q, λ[j,i])
            MPZ.sub!(λ[j,k], tmpλ[j])
        end
    end
    nothing
end

function gcd_swap!(k, m, λ, D, a, B, tmpa, tmpD, tmpB, tmpλ, tmpt)
    @inbounds begin
        a[k], a[k-1] = a[k-1], a[k]
        B[:,k], B[:,k-1] = B[:,k-1], B[:,k]
        @simd for j in 1:k-2
            λ[j,k], λ[j,k-1] = λ[j,k-1], λ[j,k]
        end
        @simd for i in k+1:m
            MPZ.mul!(tmpB[i], λ[k-1,i], D[k+1])
            MPZ.mul!(tmpλ[i], λ[k,i], λ[k-1,k])
            MPZ.sub!(tmpt[i], tmpB[i], tmpλ[i])
            MPZ.mul!(tmpB[i], λ[k-1,i], λ[k-1,k])
            MPZ.mul!(tmpλ[i], λ[k,i], D[k-1])
            MPZ.add!(tmpλ[i], tmpB[i])
            MPZ.cdiv_q!(λ[k-1,i], tmpλ[i], D[k])
            MPZ.cdiv_q!(λ[k,i], tmpt[i], D[k])
        end
        MPZ.mul!(tmpD, D[k-1], D[k+1])
        MPZ.pow_ui!(tmpa, λ[k-1,k], UInt(2))
        MPZ.add!(tmpD, tmpa)
        MPZ.cdiv_q!(D[k], tmpD, D[k])
    end
    nothing
end

function extended_gcd(s)
    m = length(s)
    B = BigInt[i==j for i in 1:m, j in 1:m] # LinearAlgebra.I defaults to fill
    λ = BigInt[0 for _ in 1:m, _ in 1:m] # zeros defaults to fill
    # λ = zeros(Int, m, m)
    D = BigInt[1 for _ in 1:m+1]
    a = BigInt.(abs.(s))

    tmpa = BigInt()
    tmpD = BigInt()
    tmpB = [BigInt() for _ in 1:m]
    tmpλ = [BigInt() for _ in 1:m]
    tmpt = [BigInt() for _ in 1:m]

    k = 2
    @inbounds while k <= m
        gcd_reduce1!(k, k-1, m, λ, D, a, B, tmpa, tmpD, tmpB, tmpλ)
        if !iszero(a[k-1]) || ((iszero(a[k-1]) & iszero(a[k])) &&
                               4*(D[k-1]*D[k+1] + λ[k-1,k]^2) < 3*D[k]^2)
            gcd_swap!(k, m, λ, D, a, B, tmpa, tmpD, tmpB, tmpλ, tmpt)
            k -= (k > 2)
        else
            for i in k-2:-1:1
                gcd_reduce1!(k, i, m, λ, D, a, B, tmpa, tmpD, tmpB, tmpλ)
            end
            k += 1
        end
    end
    @inbounds if a[end] < 0
        MPZ.neg!(a[end])
        foreach(MPZ.neg!, @view B[:,end])
    end
    i = 1
    @inbounds for x in s
        if x < 0
            MPZ.neg!(B[i,end])
        end
        i+=1
    end
    return @inbounds (a[end], B[:,end])
end

# Normal form for a list of vectors in ℚ³

function nf3D(list::AbstractVector{<:AbstractVector{<:Rational{T}}}) where T
    U = widen(T)
    n = length(list)
    ratbasis = zero(SizedMatrix{3,3,Rational{U}})
    # We start by computing a basis for the i first vectors of v where i (≤ n) is minimal
    i = 1
    @inbounds while i <= n
        if !iszero(list[i])
            ratbasis[:,3] .= list[i]
            i += 1
            break
        end
        i += 1
    end
    @inbounds while i <= n
        v = list[i]
        valid = false
        # v is valid if it is linearly independant of ratbasis[:,3]
        if iszero(v[1])
            valid = (!iszero(ratbasis[1,3]) & (!iszero(v[2]) | !iszero(v[3]))) ||
                    (v[2] * ratbasis[3,3] != v[3] * ratbasis[2,3])
        else
            x = ratbasis[1,3] // v[1]
            valid = ratbasis[2,3] != x * v[2] || ratbasis[3,3] != x * v[3]
        end
        if valid
            ratbasis[:,2] .= v
            i += 1
            break
        end
        i += 1
    end
    d = zero(Rational{U})
    @inbounds while i <= n
        v = list[i]
        ratbasis[:,1] .= v
        d = det(ratbasis)
        iszero(d) || break
        i += 1
    end
    @assert !iszero(d) # equivalent to @assert !iszero(det(ratbasis))
    # At this point ratbasis is a basis in which the i-1 first coordinates are
    # [0, 0, 1] followed by a (possibly empty) series of [0, 0, x], then
    # [0, 1, 0] followed by a (possibly empty) series of [0, x, y], then [1, 0, 0]
    @inbounds if d < 0
        ratbasis[:,1] .= .-ratbasis[:,1] # Ensure a direct transformation
    end

    lcms = SizedVector{3,T}([1,1,1])
    # per-column lcm of the vectors expressed in ratbasis
    expressed_in_ratbasis = Vector{SVector{3,Rational{T}}}(undef, n)
    invratbasis = inv(SMatrix{3,3,Rational{U}}(ratbasis))
    @inbounds for j in 1:n
        expressed_in_ratbasis[j] = invratbasis * list[j]
        lcms .= lcm.(lcms, denominator.(expressed_in_ratbasis[j]))
    end

    expressed_in_basis = [MVector{3,Int}(undef) for _ in 1:n]
    @inbounds for j in 1:n
        x = expressed_in_ratbasis[j]
        expressed_in_basis[j] .= div.(numerator.(x).*lcms, denominator.(x))
    end
    #= lcms contains the lcm of all denominators of the coordinates of translations
       expressed in ratbasis.
       We now work in an artificial basis in which the coordinates are all integers.
       The integer coordinates are stored in expressed_in_basis.
       expressed_in_basis is now going to be reduced to its normal form.
    =#

    # intcoords = [MVector{3,Int}(undef) for _ in 1:n]
    intbasis = zero(MMatrix{3, 3, Int}) # The artificial integer basis
    @inbounds for j in 1:3
        d, coefs = Tuple{Int, Vector{Int}}(extended_gcd(x[j] for x in expressed_in_basis))
        for i in 1:n
            intbasis[:,j] .+= coefs[i] .* expressed_in_basis[i]
        end
        @assert intbasis[j,j] > 0
        for i in 1:n
            coord = expressed_in_basis[i][j] .÷ d
            # intcoords[i][j] = coord
            expressed_in_basis[i] .-= coord .* intbasis[:,j]
        end
    end
    # normalization = MMatrix{3,3,Int}(LinearAlgebra.I)
    @inbounds for (i,j) in ((2,1), (3,1), (3,2))
        # normalization[i,j] = (intbasis[i,j] ÷ intbasis[i,i]) - signbit(intbasis[i,j])
        # intbasis[:,j] .-= normalization[i,j] * intbasis[:,i]
        x = signbit(intbasis[i,j])
        intbasis[:,j] .-= (((intbasis[i,j] + x) ÷ intbasis[i,i]) - x) * intbasis[:,i]
    end

    # We finally go back from the integer normal form to a rational basis that expresses
    # all the original vectors with integer coordinates
    return SMatrix{3,3,Rational{U}}(ratbasis * [@inbounds(intbasis[i,j] // lcms[i]) for i in 1:3, j in 1:3])

    # newcoords = [SVector{3,Int}(normalization * x) for x in intcoords]
    # @assert (newbasis,) .* newcoords == list

    # return newbasis, newcoords
    nothing
end

ext: jl


text:
using ..ChemistryFeaturization.Utils.ElementFeatureUtils
using DataFrames
using ..ChemistryFeaturization: elements
using ..ChemistryFeaturization.AbstractType: AbstractCodec, AbstractAtoms
using ..ChemistryFeaturization.Codec: OneHotOneCold

include("abstractfeatures.jl")

# TODO - consider edge cases in constructor. add this stuff into modulify.

# TODO: figure out what scheme would look like that is flexible to direct-value encoding (may just need a different feature type since it'll have to handle normalization, etc. too)
"""
    ElementFeatureDescriptor

Describe features associated with individual atoms that depend only upon their elemental identity

## Fields
- `name::String`: Name of the feature
- `encoder_decoder::AbstractCodec`: Codec defined which handles the feature's encoding and decoding logic
- `categorical::Bool`: flag for whether the feature is categorical or continuous-valued
- `lookup_table::DataFrame`: table containing values of feature for every encodable element
"""
struct ElementFeatureDescriptor <: AbstractAtomFeatureDescriptor
    name::String
    encoder_decoder::AbstractCodec
    categorical::Bool
    lookup_table::DataFrame
end

function ElementFeatureDescriptor(feature_name::String, encoder_decoder::AbstractCodec)
    lookup_table = atom_data_df

    colnames = names(lookup_table)
    @assert feature_name in colnames && "Symbol" in colnames "Your lookup table must have a column called :Symbol and one with the same name as your feature to be usable!"

    lookup_table = lookup_table[:, ["Symbol", feature_name]]
    dropmissing!(lookup_table)

    ElementFeatureDescriptor(
        feature_name,
        encoder_decoder,
        default_categorical(feature_name, lookup_table),
        lookup_table,
    )
end

"""
    ElementFeatureDescriptor(feature_name, lookup_table, categorical, contextual, length, encodable_elements)

Construct a feature object that encodes features associated with individual atoms that depend only upon their elemental identity.
If a Codec isn't explicity specified, [OneHotOneCold](@ref) with [default_efd_encode](@ref) and [default_efd_decode](@ref)
as the encoding and decoding functions respectively is the default choice.

## Arguments
- `name::String`: the name of the feature
- `lookup_table::DataFrame`: table containing values of feature for every encodable element
- `nbins::Integer`: Number of bins to use for one-cold decoding of continuous-valued features
- `logspaced::Bool`: whether onehot-style bins should be logarithmically spaced or not
- `categorical::Bool`: flag for whether the feature is categorical or continuous-valued
"""
function ElementFeatureDescriptor(
    feature_name::String,
    lookup_table::DataFrame = atom_data_df;
    nbins::Integer = default_nbins,
    logspaced::Bool = default_log(feature_name, lookup_table),
    categorical::Bool = default_categorical(feature_name, lookup_table),
)
    colnames = names(lookup_table)
    @assert feature_name in colnames && "Symbol" in colnames "Your lookup table must have a column called :Symbol and one with the same name as your feature to be usable!"

    lookup_table = lookup_table[:, ["Symbol", feature_name]]
    dropmissing!(lookup_table)

    ElementFeatureDescriptor(
        feature_name,
        OneHotOneCold(default_efd_encode, default_efd_decode, nbins, logspaced),
        categorical,
        lookup_table,
    )
end

# pretty printing, short version
Base.show(io::IO, efd::ElementFeatureDescriptor) = print(io, "ElementFeature $(efd.name)")

# pretty printing, long version
function Base.show(io::IO, ::MIME"text/plain", efd::ElementFeatureDescriptor)
    st = "ElementFeature $(efd.name):\n   categorical: $(efd.categorical)\n   encoded length: $(output_shape(efd))"
    print(io, st)
end

encodable_elements(efd::ElementFeatureDescriptor) = efd.lookup_table[:, :Symbol]

function encodable_elements(feature_name::String, lookup_table::DataFrame = atom_data_df)
    info = lookup_table[:, [Symbol(feature_name), :Symbol]]
    return info[
        findall(x -> !ismissing(x), getproperty(info, Symbol(feature_name))),
        :Symbol,
    ]
end

function (efd::ElementFeatureDescriptor)(a::AbstractAtoms)
    @assert all([el in encodable_elements(efd) for el in elements(a)]) "Feature $(efd.name) cannot encode some element(s) in this structure!"
    efd.encoder_decoder(efd, a)
end

(ed::OneHotOneCold)(efd::ElementFeatureDescriptor, a::AbstractAtoms) =
    ed.encode_f(efd, a, ed.nbins, ed.logspaced)

(ed::OneHotOneCold)(efd::ElementFeatureDescriptor, encoded_feature) =
    ed.decode_f(efd, encoded_feature, ed.nbins, ed.logspaced)

"""
    output_shape(efd::ElementFeatureDescriptor)

Get the output-shape for an ElementFeatureDescriptor object using the logic assoicated with its
Codec.
"""
function output_shape(efd::ElementFeatureDescriptor, ed::OneHotOneCold)
    return efd.categorical ? length(unique(efd.lookup_table[:, Symbol(efd.name)])) :
           ed.nbins
end

"""
    default_efd_encode(efd::ElementFeatureDescriptor, a::AbstractAtoms, nbins::Integer, logspaced::Bool)

Default one-hot encoding function for an ElementFeatureDescriptor object.
"""
function default_efd_encode(
    efd::ElementFeatureDescriptor,
    a::AbstractAtoms,
    nbins::Integer,
    logspaced::Bool,
)
    reduce(
        hcat,
        map(
            e -> onehot_lookup_encoder(
                e,
                efd.name,
                efd.lookup_table;
                nbins = nbins,
                logspaced = logspaced,
                categorical = efd.categorical,
            ),
            elements(a),
        ),
    )
end

"""
    default_efd_decode(efd::ElementFeatureDescriptor, encoded_feature, nbins::Integer, logspaced::Bool)

Default one-cold decoding logic for an ElementFeatureDescriptor object.
"""
default_efd_decode(efd::ElementFeatureDescriptor, encoded_feature, nbins, logspaced) =
    onecold_decoder(
        encoded_feature,
        efd.name,
        efd.lookup_table;
        nbins = nbins,
        logspaced = logspaced,
        categorical = efd.categorical,
    )


text:
export Finite

struct Finite{T <: Real} <: ConvexSet
    x::Matrix{T}
end

using Distances

function project!(s::Finite{<:Real}, v::Vector{T}, y::Vector{T}) where {T <: Real}
    dist = Distances.colwise(Distances.SqEuclidean(), s.x, y)
    j = argmin(dist)
    v = s.x[:, j]
end

function project(s::Finite{<:Real}, y::Vector{<:Real})
    project!(s, similar(y), y)
end


text:
"""
    half-cheetah
    q = (x, z, t, q1, ..., q8)
        x - lateral position
        z - vertical position
        t - body orientation
        q4 - leg 1 hip angle   (absolute)
		q5 - leg 1 knee angle  (absolute)
        q6 - leg 1 ankle angle (absolute)
		q7 - leg 3 hip angle   (absolute)
        q8 - leg 3 knee angle  (absolute)
        q9 - leg 3 ankle angle (absolute)
"""
mutable struct Halfcheetah4{T} <: Model{T}
    # dimensions
	nq::Int # generalized coordinates
    nu::Int # controls
    nw::Int # parameters
    nc::Int # contact points

    # parameters

    # torso
    l_torso1::T
    l_torso2::T
    m_torso::T
    J_torso::T

    # leg 1
        # thigh
    l_thigh1::T
    lc_thigh1::T
    m_thigh1::T
    J_thigh1::T
        # calf
    l_calf1::T
    lc_calf1::T
    m_calf1::T
    J_calf1::T
		# foot
	l_foot1::T
    lc_foot1::T
    m_foot1::T
    J_foot1::T

	# leg 3
        # thigh
    l_thigh3::T
    lc_thigh3::T
    m_thigh3::T
    J_thigh3::T
        # calf
    l_calf3::T
    lc_calf3::T
    m_calf3::T
    J_calf3::T
		# foot
	l_foot3::T
	lc_foot3::T
	m_foot3::T
	J_foot3::T

    # joint friction
	friction_joint::Vector{T}
	joint_spring_stiffness::Vector{T}
	joint_spring_offset::Vector{T}

    # environment
    friction_body_world::Vector{T}
    friction_foot_world::Vector{T}
    gravity::T
end

function kinematics_hip(model::Halfcheetah4, q; hip=:none)
	x = q[1]
	z = q[2]
    t_torso = q[3]

    if hip == :hip1
        l = model.l_torso1
    elseif hip == :hip2
        l = -model.l_torso2
    end

    return [
            x + l * cos(t_torso);
            z + l * sin(t_torso)
           ]
end

function kinematics_jacobian_hip(model::Halfcheetah4, q; hip=:none)
	x = q[1]
	z = q[2]
    t_torso = q[3]

    if hip == :hip1
        l = model.l_torso1
    elseif hip == :hip2
        l = -model.l_torso2
    end

    jac = zeros(eltype(q), 2, model.nq)

    jac[1, 1] = 1.0
    jac[1, 3] = -l * sin(t_torso)
    jac[2, 2] = 1.0
    jac[2, 3] = l * cos(t_torso)

    return jac
end

function kinematics_thigh(model::Halfcheetah4, q; leg=:none, mode = :ee)
	x = q[1]
	z = q[2]
    t_torso = q[3]

    if leg == :leg1
        t_hip = q[4]
        l_torso = model.l_torso1
        le_thigh = model.l_thigh1
        lc_thigh = model.lc_thigh1
        hip = :hip1
    elseif leg == :leg3
        t_hip = q[7]
        l_torso = -model.l_torso2
        le_thigh = model.l_thigh3
        lc_thigh = model.lc_thigh3
        hip = :hip2
    else
        @error "incorrect leg specified"
    end

    if mode == :ee
        l_thigh = le_thigh
    elseif mode == :com
        l_thigh = lc_thigh
    else
        @error "incorrect mode specified"
    end
    k_hip = kinematics_hip(model, q, hip=hip)

    return k_hip + [l_thigh * sin(t_hip);
                   -l_thigh * cos(t_hip)]
end

function kinematics_jacobian_thigh(model::Halfcheetah4, q; leg=:none, mode = :ee)
	x = q[1]
	z = q[2]
    t_torso = q[3]

    if leg == :leg1
        idx = 4
        l_torso = model.l_torso1
        le_thigh = model.l_thigh1
        lc_thigh = model.lc_thigh1
        hip = :hip1
    elseif leg == :leg3
        idx = 7
        l_torso = -model.l_torso2
        le_thigh = model.l_thigh3
        lc_thigh = model.lc_thigh3
        hip = :hip2
    else
        @error "incorrect leg specified"
    end

    t_hip = q[idx]

    if mode == :ee
        l_thigh = le_thigh
    elseif mode == :com
        l_thigh = lc_thigh
    else
        @error "incorrect mode specified"
    end

    jac = kinematics_jacobian_hip(model, q, hip=hip)

    jac[1, idx] += l_thigh * cos(t_hip)
    jac[2, idx] += l_thigh * sin(t_hip)

    return jac
end

function kinematics_calf(model::Halfcheetah4, q; leg=:none, mode = :none)
	x = q[1]
	z = q[2]
    t_torso = q[3]

    if leg == :leg1
        idx = 4
        l_torso = model.l_torso1
        l_thigh = model.l_thigh1
        le_calf = model.l_calf1
        lc_calf = model.lc_calf1
    elseif leg == :leg3
        idx = 7
        l_torso = -model.l_torso2
        l_thigh = model.l_thigh3
        le_calf = model.l_calf3
        lc_calf = model.lc_calf3
    else
        @error "incorrect leg specified"
    end

    t_hip = q[idx]
    t_calf = q[idx + 1]

    if mode == :ee
        l_calf = le_calf
    elseif mode == :com
        l_calf = lc_calf
    else
        @error "incorrect mode specified"
    end

    k_thigh = kinematics_thigh(model, q, leg=leg, mode=:ee)

    return k_thigh + [l_calf * sin(t_calf);
                     -l_calf * cos(t_calf)]
end

function kinematics_jacobian_calf(model::Halfcheetah4, q; leg=:none, mode = :none)
	x = q[1]
	z = q[2]
    t_torso = q[3]

    if leg == :leg1
        idx = 4
        l_torso = model.l_torso1
        l_thigh = model.l_thigh1
        le_calf = model.l_calf1
        lc_calf = model.lc_calf1
    elseif leg == :leg3
        idx = 7
        l_torso = -model.l_torso2
        l_thigh = model.l_thigh3
        le_calf = model.l_calf3
        lc_calf = model.lc_calf3
    else
        @error "incorrect leg specified"
    end

    t_hip = q[idx]
    t_calf = q[idx + 1]

    if mode == :ee
        l_calf = le_calf
    elseif mode == :com
        l_calf = lc_calf
    else
        @error "incorrect mode specified"
    end

    jac = kinematics_jacobian_thigh(model, q, leg=leg, mode=:ee)
    jac[1, idx + 1] += l_calf * cos(t_calf)
    jac[2, idx + 1] += l_calf * sin(t_calf)

    return jac
end

function kinematics_foot(model::Halfcheetah4, q; leg=:none, mode = :none)
	x = q[1]
	z = q[2]
    t_torso = q[3]

    if leg == :leg1
        idx = 4
        l_torso = model.l_torso1
        l_thigh = model.l_thigh1
		l_calf = model.l_calf1
        le_foot = model.l_foot1
        lc_foot = model.lc_foot1
    elseif leg == :leg3
        idx = 7
        l_torso = -model.l_torso2
        l_thigh = model.l_thigh3
        l_calf = model.l_calf3
		le_foot = model.l_foot3
		lc_foot = model.lc_foot3
    else
        @error "incorrect leg specified"
    end

    t_hip = q[idx]
	t_calf = q[idx + 1]
    t_foot = q[idx + 2]

    if mode == :ee
        l_foot = le_foot
    elseif mode == :com
        l_foot = lc_foot
    else
        @error "incorrect mode specified"
    end

    k_calf = kinematics_calf(model, q, leg=leg, mode=:ee)

    return k_calf + [l_foot * sin(t_foot);
                     -l_foot * cos(t_foot)]
end

function kinematics_jacobian_foot(model::Halfcheetah4, q; leg=:none, mode = :none)
	x = q[1]
	z = q[2]
    t_torso = q[3]

    if leg == :leg1
        idx = 4
		l_torso = model.l_torso1
        l_thigh = model.l_thigh1
		l_calf = model.l_calf1
        le_foot = model.l_foot1
        lc_foot = model.lc_foot1
    elseif leg == :leg3
        idx = 7
		l_torso = -model.l_torso2
        l_thigh = model.l_thigh3
        l_calf = model.l_calf3
		le_foot = model.l_foot3
		lc_foot = model.lc_foot3
    else
        @error "incorrect leg specified"
    end

    t_hip = q[idx]
	t_calf = q[idx + 1]
    t_foot = q[idx + 2]

    if mode == :ee
        l_foot = le_foot
    elseif mode == :com
        l_foot = lc_foot
    else
        @error "incorrect mode specified"
    end

    jac = kinematics_jacobian_calf(model, q, leg=leg, mode=:ee)
    jac[1, idx + 2] += l_foot * cos(t_foot)
    jac[2, idx + 2] += l_foot * sin(t_foot)

    return jac
end

function kinematics(model::Halfcheetah4, q)
	p_toe_1 = kinematics_foot(model, q, leg=:leg1, mode=:ee) # toe 1
	p_toe_3 = kinematics_foot(model, q, leg=:leg3, mode=:ee) # toe 3

	p_heel_1 = kinematics_calf(model, q, leg=:leg1, mode=:ee) # heel 1
	p_heel_3 = kinematics_calf(model, q, leg=:leg3, mode=:ee) # heel 3

    # p_knee_1 = kinematics_thigh(model, q, leg=:leg1, mode=:ee) # knee 1
	# p_knee_3 = kinematics_thigh(model, q, leg=:leg3, mode=:ee) # knee 3

    # p_hip_1 = kinematics_hip(model, q, hip=:hip1) # hip 1
	# p_hip_2 = kinematics_hip(model, q, hip=:hip2) # hip 2

	[
	 p_toe_1; p_toe_3;
     p_heel_1; p_heel_3;
    #  p_knee_1; p_knee_3;
    #  p_hip_1; p_hip_2
    ]
end

# Lagrangian
function lagrangian(model::Halfcheetah4, q, q̇)
	L = 0.0

	# torso
	p_torso = q[1:2]
	v_torso = q̇[1:2]

	L += 0.5 * model.m_torso * transpose(v_torso) * v_torso
	L += 0.5 * model.J_torso * q̇[3]^2.0
	L -= model.m_torso * model.gravity * p_torso[2]

	# thigh 1
	p_thigh_1 = kinematics_thigh(model, q, leg=:leg1, mode=:com)
	J_thigh_1 = kinematics_jacobian_thigh(model, q, leg=:leg1, mode=:com)
	v_thigh_1 = J_thigh_1 * q̇

	L += 0.5 * model.m_thigh1 * transpose(v_thigh_1) * v_thigh_1
	L += 0.5 * model.J_thigh1 * q̇[4]^2.0
	L -= model.m_thigh1 * model.gravity * p_thigh_1[2]

	# calf 1
	p_calf_1 = kinematics_calf(model, q, leg=:leg1, mode=:com)
	J_calf_1 = kinematics_jacobian_calf(model, q, leg=:leg1, mode=:com)
	v_calf_1 = J_calf_1 * q̇

	L += 0.5 * model.m_calf1 * transpose(v_calf_1) * v_calf_1
	L += 0.5 * model.J_calf1 * q̇[5]^2.0
	L -= model.m_calf1 * model.gravity * p_calf_1[2]

	# foot 1
	p_foot_1 = kinematics_foot(model, q, leg=:leg1, mode=:com)
	J_foot_1 = kinematics_jacobian_foot(model, q, leg=:leg1, mode=:com)
	v_foot_1 = J_foot_1 * q̇

	L += 0.5 * model.m_foot1 * transpose(v_foot_1) * v_foot_1
	L += 0.5 * model.J_foot1 * q̇[6]^2.0
	L -= model.m_foot1 * model.gravity * p_foot_1[2]

	# thigh 3
	p_thigh_3 = kinematics_thigh(model, q, leg=:leg3, mode=:com)
	J_thigh_3 = kinematics_jacobian_thigh(model, q, leg=:leg3, mode=:com)
	v_thigh_3 = J_thigh_3 * q̇

	L += 0.5 * model.m_thigh3 * transpose(v_thigh_3) * v_thigh_3
	L += 0.5 * model.J_thigh3 * q̇[7]^2.0
	L -= model.m_thigh3 * model.gravity * p_thigh_3[2]

	# calf 3
	p_calf_3 = kinematics_calf(model, q, leg=:leg3, mode=:com)
	J_calf_3 = kinematics_jacobian_calf(model, q, leg=:leg3, mode=:com)
	v_calf_3 = J_calf_3 * q̇

	L += 0.5 * model.m_calf3 * transpose(v_calf_3) * v_calf_3
	L += 0.5 * model.J_calf3 * q̇[8]^2.0
	L -= model.m_calf3 * model.gravity * p_calf_3[2]

	# foot 3
	p_foot_3 = kinematics_foot(model, q, leg=:leg3, mode=:com)
	J_foot_3 = kinematics_jacobian_foot(model, q, leg=:leg3, mode=:com)
	v_foot_3 = J_foot_3 * q̇

	L += 0.5 * model.m_foot3 * transpose(v_foot_3) * v_foot_3
	L += 0.5 * model.J_foot3 * q̇[8]^2.0
	L -= model.m_foot3 * model.gravity * p_foot_3[2]

	t_torso = q[3]
	t_hip1 = q[4]
	t_knee1 = q[5]
	t_foot1 = q[6]
	t_hip2 = q[7]
	t_knee3 = q[8]
	t_foot3 = q[9]
	L -= 0.5 * joint_spring_stiffness[1] * (-joint_spring_offset[1] + t_torso - t_hip1)^2
	L -= 0.5 * joint_spring_stiffness[2] * (-joint_spring_offset[2] + t_hip1  - t_knee1)^2
	L -= 0.5 * joint_spring_stiffness[3] * (-joint_spring_offset[3] + t_knee1 - t_foot1)^2
	L -= 0.5 * joint_spring_stiffness[4] * (-joint_spring_offset[4] + t_torso - t_hip2)^2
	L -= 0.5 * joint_spring_stiffness[5] * (-joint_spring_offset[5] + t_hip2  - t_knee3)^2
	L -= 0.5 * joint_spring_stiffness[6] * (-joint_spring_offset[6] + t_knee3 - t_foot3)^2
	return L
end

function signed_distance(model::Halfcheetah4, q)
	k = kinematics(model, q)
    return [k[2], k[4], k[6], k[8]]#, k[10], k[12], k[14], k[16]]
end

function input_jacobian(::Halfcheetah4, q)
	# x    z    θ    h1   k1   a1   h2   k2   a2
	[1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0; # infeasible controls
	 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0; # infeasible controls
	 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0; # infeasible controls
     0.0  0.0 -1.0  1.0  0.0  0.0  0.0  0.0  0.0;
	 0.0  0.0  0.0 -1.0  1.0  0.0  0.0  0.0  0.0;
     0.0  0.0  0.0  0.0 -1.0  1.0  0.0  0.0  0.0;
	 0.0  0.0 -1.0  0.0  0.0  0.0  1.0  0.0  0.0;
     0.0  0.0  0.0  0.0  0.0  0.0 -1.0  1.0  0.0;
     0.0  0.0  0.0  0.0  0.0  0.0  0.0 -1.0  1.0]
end

function contact_jacobian(model::Halfcheetah4, q)
	J1 = kinematics_jacobian_foot(model, q, leg=:leg1, mode=:ee)
	J2 = kinematics_jacobian_foot(model, q, leg=:leg3, mode=:ee)

    J3 = kinematics_jacobian_calf(model, q, leg=:leg1, mode=:ee)
	J4 = kinematics_jacobian_calf(model, q, leg=:leg3, mode=:ee)

    # J5 = kinematics_jacobian_thigh(model, q, leg=:leg1, mode=:ee)
	# J6 = kinematics_jacobian_thigh(model, q, leg=:leg3, mode=:ee)

    # J7 = kinematics_jacobian_hip(model, q, hip=:hip1)
    # J8 = kinematics_jacobian_hip(model, q, hip=:hip2)

    [
     J1;
     J2;
     J3;
	 J4;
	#  J5;
	#  J6;
    #  J7;
    #  J8;
    ]
end

# nominal configuration
function nominal_configuration(model::Halfcheetah4)
	[0.0, 0.6, 0.0, -0.52, +0.59159, 0.69159, +0.65800, -1.042, 0.458]
end
# nominal state
function nominal_state(model::Halfcheetah4)
	[[0.0, 0.6, 0.0, -0.52, +0.59159, 0.69159, +0.65800, -1.042, 0.458]; zeros(model.nq)]
end

# friction coefficients
function friction_coefficients(model::Halfcheetah4)
    [model.friction_foot_world; model.friction_body_world]
end

# dimensions
nq = 3 + 2 * 3            # generalized coordinates
nu = 3 + 2 * 3            # controls
nw = 0                    # parameters
nc = 4                    # contact points

# parameters
gravity = 9.81            # gravity
friction_body_world = zeros(0)#0.4; 0.4; 0.4; 0.4] # coefficient of friction
friction_foot_world = [0.4; 0.4; 0.4; 0.4] # coefficient of friction
friction_joint = [4.5, 3.0, 1.5, 6.0, 4.5, 3.0] ./ 100 # fthigh, fcalf, ffoot, bthigh, bcalf, bfoot
joint_spring_stiffness = [180.0, 120.0, 60.0, 240.0, 180.0, 120.0] ./ 5 # fthigh, fcalf, ffoot, bthigh, bcalf, bfoot
joint_spring_offset = [0.52, -1.11159, -0.10000, -0.65800, 1.7, -1.5] # fthigh, fcalf, ffoot, bthigh, bcalf, bfoot

# similar to Unitree A1
scale = 10.0
m_torso = 6.36031 / scale
m_thigh1 = 1.42558 / scale
m_calf1 = 1.17885 / scale
m_foot1 = 0.84986 / scale
m_thigh3 = 1.5352 / scale
m_calf3 = 1.5809 / scale
m_foot3 = 1.069 / scale

J_torso = 0.08944 / scale
J_thigh1 = 0.01231 / scale
J_calf1 = 0.00716 / scale
J_foot1 = 0.002899 / scale
J_thigh3 = 0.01525 / scale
J_calf3 = 0.01660 / scale
J_foot3 = 0.005444 / scale

l_torso1 = 0.5
l_torso2 = 0.5
l_thigh1 = 0.266
l_calf1 = 0.212
l_foot1 = 0.140
l_thigh3 = 0.290
l_calf3 = 0.290
l_foot3 = 0.188

lc_thigh1 = 0.5 * l_thigh1
lc_calf1 = 0.5 * l_calf1
lc_foot1 = 0.5 * l_foot1

lc_thigh3 = 0.5 * l_thigh3
lc_calf3 = 0.5 * l_calf3
lc_foot3 = 0.5 * l_foot3

halfcheetah4 = Halfcheetah4(nq, nu, nw, nc,
				l_torso1, l_torso2, m_torso, J_torso,

				l_thigh1, lc_thigh1, m_thigh1, J_thigh1,
				l_calf1, lc_calf1, m_calf1, J_calf1,
				l_foot1, lc_foot1, m_foot1, J_foot1,

				l_thigh3, lc_thigh3, m_thigh3, J_thigh3,
				l_calf3, lc_calf3, m_calf3, J_calf3,
				l_foot3, lc_foot3, m_foot3, J_foot3,

	            friction_joint,
				joint_spring_stiffness,
				joint_spring_offset,
                friction_body_world,
                friction_foot_world,
                gravity)

halfcheetah4_contact_kinematics = [
	q -> kinematics_foot(halfcheetah4, q, leg=:leg1, mode=:ee),
	q -> kinematics_foot(halfcheetah4, q, leg=:leg3, mode=:ee),
	q -> kinematics_calf(halfcheetah4, q, leg=:leg1, mode=:ee),
    q -> kinematics_calf(halfcheetah4, q, leg=:leg3, mode=:ee),
    # q -> kinematics_thigh(halfcheetah4, q, leg=:leg1, mode=:ee),
    # q -> kinematics_thigh(halfcheetah4, q, leg=:leg3, mode=:ee),
    # q -> kinematics_hip(halfcheetah4, q, hip=:hip1),
    # q -> kinematics_hip(halfcheetah4, q, hip=:hip2)
    ]

halfcheetah4_contact_kinematics_jacobians = [
	q -> kinematics_jacobian_foot(halfcheetah4, q, leg=:leg1, mode=:ee),
	q -> kinematics_jacobian_foot(halfcheetah4, q, leg=:leg3, mode=:ee),
 	q -> kinematics_jacobian_calf(halfcheetah4, q, leg=:leg1, mode=:ee),
    q -> kinematics_jacobian_calf(halfcheetah4, q, leg=:leg3, mode=:ee),
    # q -> kinematics_jacobian_thigh(halfcheetah4, q, leg=:leg1, mode=:ee),
    # q -> kinematics_jacobian_thigh(halfcheetah4, q, leg=:leg3, mode=:ee),
    # q -> kinematics_jacobian_hip(halfcheetah4, q, hip=:hip1),
    # q -> kinematics_jacobian_hip(halfcheetah4, q, hip=:hip2)
]

name(::Halfcheetah4) = :halfcheetah4


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
function f(L)
    a = []
    for i in 1:L
        push!(a, rand())
    end
    a
end

function g(L)
    a = Float64[]
    for i in 1:L
        push!(a, rand())
    end
    a
end

function h(L)
    a = Vector{Float64}(undef, L)
    for i in eachindex(a)
        a[i] = rand()
    end
    a
end

f(1); g(1); h(1); # compile

# %%
print("a = [] and push!:                           ")
@time a = f(10^6)
print("a = Float64[] and push!:                    ")
@time b = g(10^6)
print("a = Vector{Float64}(undef, L) and a[i] =...:")
@time c = h(10^6);

# %%
using BenchmarkTools

println("Test of constructors")
print("a = [] and push!:                          ")
@btime a = f(10^6)
print("a = Float64[] and push!:                    ")
@btime b = g(10^6)
print("a = Vector{Float64}(undef, L) and a[i] =...:")
@btime c = h(10^6);

# %%
println("Test of sum")
print("a = [] and push!:                        ")
@btime sum($a)
print("a = Float64[] and push!:                    ")
@btime sum($b)
print("a = Vector{Float64}(undef, L) and a[i] =...:")
@btime sum($c);

# %%
function h_inbounds(L)
    a = Vector{Float64}(undef, L)
    @inbounds for i in eachindex(a)
        a[i] = rand()
    end
    a
end

@btime h_inbounds(10^6);

# %%
@btime rand(10^6);

# %%

max_stars_repo_path: 0016/Don't use push!.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0016/Don't use push!.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0016/Don't use push!.jl

max_forks_repo_name: genkuroki/public


text:
function create_filter(y;
    model,
    constant_in_time::Bool = true,
    input_variables::Bool = false)

    input = Input(model)

    options = Options(constant_in_time)

    sizes = load_sizes(input, options)

    prb = Problem(input, options, sizes)

    return prb
end

function load_sizes(input::Input, options::Options)

    model = input.model

    p, n = size(model[:y])
    m, r = size(model[:R])

    sizes = Sizes(n, p, m, r)
    return sizes
end


text:
"""
    fit!(gcm::NBCopulaVCModel, solver=Ipopt.IpoptSolver)

Fit an `NBCopulaVCModel` object by block MLE using a nonlinear programming solver.
Start point should be provided in `gcm.β`, `gcm.θ`, `gcm.r`.
In our block updates, we fit 15 iterations of `gcm.β`, `gcm.θ` using IPOPT, followed by 10 iterations of
Newton on nuisance parameter `gcm.r`. Convergence is declared when difference of
successive loglikelihood is less than `tol`.

# Arguments
- `gcm`: A `NBCopulaVCModel` model object.
- `solver`: Specified solver to use. By default we use IPOPT with 100 quas-newton iterations with convergence tolerance 10^-6.
    (default `solver = Ipopt.IpoptSolver(print_level = 0, max_iter = 15, limited_memory_max_history = 20,
                            warm_start_init_point = "yes",  mu_strategy = "adaptive",
                            hessian_approximation = "limited-memory")`)

# Optional Arguments
- `tol`: Convergence tolerance for the max block iter updates (default `tol = 1e-6`).
- `maxBlockIter`: Number of maximum block iterations to update `gcm.β`, `gcm.θ` and  `gcm.r` (default `maxBlockIter = 10`).
"""
function fit!(
        gcm::NBCopulaVCModel,
        solver=Ipopt.IpoptSolver(print_level = 0, max_iter = 15, limited_memory_max_history = 20,
                                warm_start_init_point = "yes",  mu_strategy = "adaptive",
                                hessian_approximation = "limited-memory");
        tol::Float64 = 1e-6,
        maxBlockIter::Int=10
    )
    # initialize model
    initialize_model!(gcm)
    npar = gcm.p + gcm.m
    optm = MathProgBase.NonlinearModel(solver)
    # set lower bounds and upper bounds of parameters
    lb = fill(-Inf, npar)
    ub = fill( Inf, npar)
    offset = gcm.p + 1
    for k in 1:gcm.m
        lb[offset] = 0
        offset += 1
    end
    MathProgBase.loadproblem!(optm, npar, 0, lb, ub, Float64[], Float64[], :Max, gcm)
    # starting point
    par0 = zeros(npar)
    modelpar_to_optimpar!(par0, gcm)
    MathProgBase.setwarmstart!(optm, par0)
    logl0 = MathProgBase.getobjval(optm)
    println("Converging when tol ≤ $tol (max block iter = $maxBlockIter)")
    # optimize
    for i in 1:maxBlockIter
        MathProgBase.optimize!(optm)
        logl = MathProgBase.getobjval(optm)
        update_r!(gcm)
        if abs(logl - logl0) / (1 + abs(logl0)) ≤ tol # this is faster but has wider confidence intervals
        # if abs(logl - logl0) ≤ tol # this is slower but has very tight confidence intervals
            break
        else
            println("Block iter $i r = $(round(gcm.r[1], digits=2))," *
            " logl = $(round(logl, digits=2)), tol = $(abs(logl - logl0) / (1 + abs(logl0)))")
            logl0 = logl
        end
    end
    # update parameters and refresh gradient
    optimpar_to_modelpar!(gcm, MathProgBase.getsolution(optm))
    loglikelihood!(gcm, true, false)
    # gcm
end

"""
    fit!(gcm::NBCopulaARModel, solver=Ipopt.IpoptSolver)

Fit an `NBCopulaARModel` object by block MLE using a nonlinear programming solver.
Start point should be provided in `gcm.β`, `gcm.ρ`, `gcm.σ2`, `gcm.r`.
In our block updates, we fit 15 iterations of `gcm.β`, `gcm.ρ`, `gcm.σ2` using IPOPT, followed by 10 iterations of
Newton on nuisance parameter `gcm.r`. Convergence is declared when difference of
successive loglikelihood is less than `tol`.

# Arguments
- `gcm`: A `NBCopulaARModel` model object.
- `solver`: Specified solver to use. By default we use IPOPT with 100 quas-newton iterations with convergence tolerance 10^-6.
    (default `solver = Ipopt.IpoptSolver(print_level = 0, max_iter = 15, limited_memory_max_history = 20,
                            warm_start_init_point = "yes",  mu_strategy = "adaptive",
                            hessian_approximation = "limited-memory")`)

# Optional Arguments
- `tol`: Convergence tolerance for the max block iter updates (default `tol = 1e-6`).
- `maxBlockIter`: Number of maximum block iterations to update `gcm.β`, `gcm.θ` and  `gcm.r` (default `maxBlockIter = 10`).
"""
function fit!(
    gcm::NBCopulaARModel,
    solver=Ipopt.Ipopt.IpoptSolver(print_level = 0, max_iter = 15,
                            limited_memory_max_history = 20,
                            warm_start_init_point = "yes",  mu_strategy = "adaptive",
                            hessian_approximation = "limited-memory");
    tol::Float64 = 1e-6,
    maxBlockIter::Int=10
    )
    initialize_model!(gcm)
    npar = gcm.p + 2 # rho and sigma squared
    optm = MathProgBase.NonlinearModel(solver)
    # set lower bounds and upper bounds of parameters
    lb   = fill(-Inf, npar)
    ub   = fill(Inf, npar)
    offset = gcm.p + 1
    ub[offset] = 1
    for k in 1:2
        lb[offset] = 0
        offset += 1
    end
    MathProgBase.loadproblem!(optm, npar, 0, lb, ub, Float64[], Float64[], :Max, gcm)
    # starting point
    par0 = zeros(npar)
    modelpar_to_optimpar!(par0, gcm)
    MathProgBase.setwarmstart!(optm, par0)
    logl0 = MathProgBase.getobjval(optm)
    println("Converging when tol ≤ $tol (max block iter = $maxBlockIter)")
    # optimize
    for i in 1:maxBlockIter
        MathProgBase.optimize!(optm)
        logl = MathProgBase.getobjval(optm)
        update_r!(gcm)
        if abs(logl - logl0) / (1 + abs(logl0)) ≤ tol # this is faster but has wider confidence intervals
            # if abs(logl - logl0) ≤ tol # this is slower but has very tight confidence intervals
            break
        else
            println("Block iter $i r = $(round(gcm.r[1], digits=2))," *
            " logl = $(round(logl, digits=2)), tol = $(abs(logl - logl0) / (1 + abs(logl0)))")
            logl0 = logl
        end
    end
    # update parameters and refresh gradient
    optimpar_to_modelpar!(gcm, MathProgBase.getsolution(optm))
    loglikelihood!(gcm, true, false)
    # gcm
end

"""
    fit!(gcm::NBCopulaCSModel, solver=Ipopt.IpoptSolver)

Fit an `NBCopulaCSModel` object by block MLE using a nonlinear programming solver.
Start point should be provided in `gcm.β`, `gcm.ρ`, `gcm.σ2`, `gcm.r`.
In our block updates, we fit 15 iterations of `gcm.β`, `gcm.ρ`, `gcm.σ2` using IPOPT, followed by 10 iterations of
Newton on nuisance parameter `gcm.r`. Convergence is declared when difference of
successive loglikelihood is less than `tol`.

# Arguments
- `gcm`: A `NBCopulaCSModel` model object.
- `solver`: Specified solver to use. By default we use IPOPT with 100 quas-newton iterations with convergence tolerance 10^-6.
    (default `solver = Ipopt.IpoptSolver(print_level = 0, max_iter = 15, limited_memory_max_history = 20,
                            warm_start_init_point = "yes",  mu_strategy = "adaptive",
                            hessian_approximation = "limited-memory")`)

# Optional Arguments
- `tol`: Convergence tolerance for the max block iter updates (default `tol = 1e-6`).
- `maxBlockIter`: Number of maximum block iterations to update `gcm.β`, `gcm.θ` and  `gcm.r` (default `maxBlockIter = 10`).
"""
function fit!(
    gcm::NBCopulaCSModel,
    solver=Ipopt.IpoptSolver(print_level = 0, max_iter = 15, limited_memory_max_history = 20,
                            warm_start_init_point = "yes", mu_strategy = "adaptive",
                             hessian_approximation = "limited-memory");
    tol::Float64 = 1e-6,
    maxBlockIter::Int=10
    )
    initialize_model!(gcm)
    npar = gcm.p + 2 # rho and sigma squared
    optm = MathProgBase.NonlinearModel(solver)
    # set lower bounds and upper bounds of parameters
    lb   = fill(-Inf, npar)
    ub   = fill(Inf, npar)
    offset = gcm.p + 1
    # rho
    ub[offset] = 1
    # lb[offset] = 0
    lb[offset] = -inv(gcm.data[1].n - 1)
    offset += 1
    # sigma2
    lb[offset] = 0
    MathProgBase.loadproblem!(optm, npar, 0, lb, ub, Float64[], Float64[], :Max, gcm)
    # starting point
    par0 = zeros(npar)
    modelpar_to_optimpar!(par0, gcm)
    MathProgBase.setwarmstart!(optm, par0)
    logl0 = MathProgBase.getobjval(optm)
    println("Converging when tol ≤ $tol (max block iter = $maxBlockIter)")
    # optimize
    for i in 1:maxBlockIter
        MathProgBase.optimize!(optm)
        logl = MathProgBase.getobjval(optm)
        update_r!(gcm)
        if abs(logl - logl0) / (1 + abs(logl0)) ≤ tol # this is faster but has wider confidence intervals
            # if abs(logl - logl0) ≤ tol # this is slower but has very tight confidence intervals
            println("Block iter $i r = $(round(gcm.r[1], digits=2))," *
            " logl = $(round(logl, digits=2)), tol = $(abs(logl - logl0) / (1 + abs(logl0)))")
            break
        else
            println("Block iter $i r = $(round(gcm.r[1], digits=2))," *
            " logl = $(round(logl, digits=2)), tol = $(abs(logl - logl0) / (1 + abs(logl0)))")
            logl0 = logl
        end
    end
    # update parameters and refresh gradient
    optimpar_to_modelpar!(gcm, MathProgBase.getsolution(optm))
    loglikelihood!(gcm, true, false)
    # gcm
end

"""
    modelpar_to_optimpar!(par, gcm)

Translate model parameters in `gcm` to optimization variables in `par`.
"""
function modelpar_to_optimpar!(
        par :: Vector,
        gcm :: NBCopulaVCModel
    )
    # β
    copyto!(par, gcm.β)
    # L
    offset = gcm.p + 1
    @inbounds for k in 1:gcm.m
        par[offset] = gcm.θ[k]
        offset += 1
    end
    par
end

"""
    optimpar_to_modelpar!(gcm, par)

Translate optimization variables in `par` to the model parameters in `gcm`.
"""
function optimpar_to_modelpar!(
        gcm :: NBCopulaVCModel,
        par :: Vector
    )
    # β
    copyto!(gcm.β, 1, par, 1, gcm.p)
    # L
    offset = gcm.p + 1
    @inbounds for k in 1:gcm.m
        gcm.θ[k] = par[offset]
        offset   += 1
    end
    gcm
end

function MathProgBase.initialize(
    gcm::NBCopulaVCModel,
    requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in [:Grad, :Hess])
            error("Unsupported feature $feat")
        end
    end
end

MathProgBase.features_available(gcm::NBCopulaVCModel) = [:Grad, :Hess]

function MathProgBase.eval_f(
        gcm :: NBCopulaVCModel,
        par :: Vector
    )
    optimpar_to_modelpar!(gcm, par)
    loglikelihood!(gcm, false, false) # don't need gradient here
end

function MathProgBase.eval_grad_f(
        gcm    :: NBCopulaVCModel,
        grad :: Vector,
        par  :: Vector
    )
    optimpar_to_modelpar!(gcm, par)
    obj = loglikelihood!(gcm, true, false)
    # gradient wrt β
    copyto!(grad, gcm.∇β)
    # gradient wrt L
    offset = gcm.p + 1
    @inbounds for k in 1:gcm.m
        grad[offset] = gcm.∇θ[k]
        offset += 1
    end
    obj
end

MathProgBase.eval_g(gcm::NBCopulaVCModel, g, par) = nothing
MathProgBase.jac_structure(gcm::NBCopulaVCModel) = Int[], Int[]
MathProgBase.eval_jac_g(gcm::NBCopulaVCModel, J, par) = nothing

# """
#     ◺(n::Integer)
# Triangular number n * (n+1) / 2
# """
# @inline ◺(n::Integer) = (n * (n + 1)) >> 1

function MathProgBase.hesslag_structure(gcm::NBCopulaVCModel)
    m◺ = ◺(gcm.m)
    # we work on the upper triangular part of the Hessian
    arr1 = Vector{Int}(undef, ◺(gcm.p) + m◺)
    arr2 = Vector{Int}(undef, ◺(gcm.p) + m◺)
    # Hββ block
    idx  = 1
    for j in 1:gcm.p
        for i in j:gcm.p
            arr1[idx] = i
            arr2[idx] = j
            idx += 1
        end
    end
    # Haa block
    for j in 1:gcm.m
        for i in 1:j
            arr1[idx] = gcm.p + i
            arr2[idx] = gcm.p + j
            idx += 1
        end
    end
    return (arr1, arr2)
end

function MathProgBase.eval_hesslag(
        gcm   :: NBCopulaVCModel,
        H   :: Vector{T},
        par :: Vector{T},
        σ   :: T,
        μ   :: Vector{T}
    ) where {T}
    optimpar_to_modelpar!(gcm, par)
    loglikelihood!(gcm, true, true)
    # Hβ block
    idx = 1
    @inbounds for j in 1:gcm.p, i in 1:j
        H[idx] = gcm.Hβ[i, j]
        idx   += 1
    end
    # Haa block
    @inbounds for j in 1:gcm.m, i in 1:j
        H[idx] = gcm.Hθ[i, j]
        idx   += 1
    end
    # lmul!(σ, H)
    H .*= σ
end


text:
#! /usr/bin/julia

# Rosetta Code, Knapsack problem/0-1

using MathProgBase

immutable KPDSupply{S<:String, T<:Integer}
    item::S
    weight::T
    value::T
    quant::T
end
function KPDSupply{S<:String, T<:Integer}(item::S, weight::T, value::T)
    KPDSupply(item, weight, value, one(T))
end

function solve{S<:String, T<:Integer}(gear::Array{KPDSupply{S,T},1},
                                      capacity::T)
    w = map(x->x.weight, gear)
    v = map(x->x.value, gear)
    sol = mixintprog(-v, w', '<', capacity, :Bin, 0, 1)
    sol.status == :Optimal || error("This Problem could not be solved")
    gear[sol.sol .== 1.0]
end

gear = [KPDSupply("map", 9, 150),
        KPDSupply("compass", 13, 35),
        KPDSupply("water", 153, 200),
        KPDSupply("sandwich", 50, 160),
        KPDSupply("glucose", 15, 60),
        KPDSupply("tin", 68, 45),
        KPDSupply("banana", 27, 60),
        KPDSupply("apple", 39, 40),
        KPDSupply("cheese", 23, 30),
        KPDSupply("beer", 52, 10),
        KPDSupply("suntan cream", 11, 70),
        KPDSupply("camera", 32, 30),
        KPDSupply("T-shirt", 24, 15),
        KPDSupply("trousers", 48, 10),
        KPDSupply("umbrella", 73, 40),
        KPDSupply("waterproof trousers", 42, 70),
        KPDSupply("waterproof overclothes", 43, 75),
        KPDSupply("note-case", 22, 80),
        KPDSupply("sunglasses", 7, 20),
        KPDSupply("towel", 18, 12),
        KPDSupply("socks", 4, 50),
        KPDSupply("book", 30, 10)]

pack = solve(gear, 400)

println("The hiker should pack:")
for s in pack
    println("    ", s.item)
end
println()
println("Packed Weight:  ", mapreduce(x->x.weight, +, pack))
println("Packed Value:  ", mapreduce(x->x.value, +, pack))

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
struct ByKeyTrait end
const ByKey = ByKeyTrait()

struct ByValueTrait end
const ByValue = ByValueTrait()

"""
    index_by(a::AbstractIndex, i)

Returns `ByKey` if the index of `i` should be identified by searching the keys
and returns `ByValue` if the index of `i` should be identified directly from the
values.
"""
index_by(a::AbstractIndex{name,K}, i::Function) where {name,K} = ByKey

index_by(a::AbstractIndex{name,K}, i::K) where {name,K} = ByKey
index_by(a::AbstractIndex{name,K}, i::AbstractVector{K}) where {name,K} = ByKey

index_by(a::AbstractIndex{name,K}, i::K) where {name,K<:Integer} = ByKey
index_by(a::AbstractIndex{name,K}, i::AbstractVector{K}) where {name,K<:Integer} = ByKey

index_by(a::AbstractIndex{name,K}, i::I) where {name,K,I<:Integer} = ByValue
index_by(a::AbstractIndex{name,K}, i::AbstractVector{I}) where {name,K,I<:Integer} = ByValue

@inline @propagate_inbounds function Base.to_index(a::AbstractIndex, i)
    return _to_index(index_by(a, i), a, i)
end

# _to_index
# TODO find_all should be filter where possible
@propagate_inbounds function _to_index(b::ByKeyTrait, a, i::Function)
    return __to_index(a, i, find_all(i, keys(a)))
end
@propagate_inbounds function _to_index(b::ByKeyTrait, a, i::AbstractVector)
    return __to_index(a, i, find_all(in(i), keys(a)))
end
@propagate_inbounds function _to_index(b::ByKeyTrait, a, i)
    return __to_index(a, i, find_first(==(i), keys(a)))
end
@propagate_inbounds function _to_index(b::ByValueTrait, a, i::Function)
    return __to_index(a, i, find_all(i, values(a)))
end
@propagate_inbounds function _to_index(b::ByValueTrait, a, i)
    @boundscheck if !checkindex(Bool, values(a), i)
        throw(BoundsError(a, i))
    end
    return @inbounds getindex(values(a), i)
end
@propagate_inbounds function _to_index(b::ByValueTrait, a, i::AbstractVector)
    @boundscheck if !checkindex(Bool, values(a), i)
        throw(BoundsError(a, i))
    end
    return @inbounds unsafe_reindex(a, i)
end

# __to_index
@propagate_inbounds function __to_index(a, i, idx::T) where {T<:Union{Integer,Nothing}}
    @boundscheck if T <: Nothing
        throw(BoundsError(a, i))
    end
    return @inbounds getindex(values(a), idx)
end
@propagate_inbounds function __to_index(a, i, idx::AbstractVector{T}) where {T<:Union{Integer,Nothing}}
    @boundscheck if !(T<:Integer)
        throw(BoundsError(a, i))
    end
    return unsafe_reindex(a, idx)
end


text:
module Test07

import Automa
import Automa.RegExp: @re_str
using Test

@testset "Test07" begin
    re1 = re"a.*b"
    machine = Automa.compile(re1)
    ctx = Automa.CodeGenContext()
    @eval function ismatch1(data)
        $(Automa.generate_init_code(ctx, machine))
        p_end = p_eof = lastindex(data)
        $(Automa.generate_exec_code(ctx, machine))
        return cs == 0
    end
    @test ismatch1(b"ab")
    @test ismatch1(b"azb")
    @test ismatch1(b"azzzb")
    @test !ismatch1(b"azzz")
    @test !ismatch1(b"zzzb")

    re2 = re"a\.*b"
    machine = Automa.compile(re2)
    ctx = Automa.CodeGenContext()
    @eval function ismatch2(data)
        $(Automa.generate_init_code(ctx, machine))
        p_end = p_eof = lastindex(data)
        $(Automa.generate_exec_code(ctx, machine))
        return cs == 0
    end
    @test ismatch2(b"ab")
    @test ismatch2(b"a.b")
    @test ismatch2(b"a...b")
    @test !ismatch2(b"azzzb")
    @test !ismatch2(b"a...")
    @test !ismatch2(b"...b")

    re3 = re"a\.\*b"
    machine = Automa.compile(re3)
    ctx = Automa.CodeGenContext()
    @eval function ismatch3(data)
        $(Automa.generate_init_code(ctx, machine))
        p_end = p_eof = lastindex(data)
        $(Automa.generate_exec_code(ctx, machine))
        return cs == 0
    end
    @test ismatch3(b"a.*b")
    @test !ismatch3(b"a...b")
end

end


text:
#=
Basic types for characters and strings

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md

Encodings inspired from collaborations on the following packages:
https://github.com/quinnj/Strings.jl with @quinnj (Jacob Quinn)
https://github.com/nalimilan/StringEncodings.jl with @nalimilan (Milan Bouchet-Valat)
=#
const BIG_ENDIAN    = (ENDIAN_BOM == 0x01020304)
const LITTLE_ENDIAN = !BIG_ENDIAN

const STR_KEEP_NUL    = true  # keep nul byte placed by String

struct CharSet{CS}   end
struct Encoding{Enc} end
struct CSE{CS, ENC}  end

CharSet(s)  = CharSet{Symbol(s)}()
Encoding(s) = Encoding{Symbol(s)}()
CSE(cs, e)  = CSE{CharSet(cs), Encoding(e)}()

macro cs_str(s)
    :(CharSet{$(quotesym(s))}())
end
macro enc_str(s)
    :(Encoding{$(quotesym(s))}())
end
macro cse(cs, e)
    :(CSE{CharSet{$(quotesym(cs)), $(quotesym(e))}()})
end

# Define symbols used for characters, codesets, codepoints

const _cpname1 =
    [:Text1,   # Unknown character set, 1 byte
     :ASCII,   # (7-bit subset of Unicode)
     :Latin]   # ISO-8859-1 (8-bit subset of Unicode)
const _cpname2 =
    [:Text2,   # Unknown character set, 2 byte
     :UCS2]    # BMP (16-bit subset of Unicode)
const _cpname4 =
    [:Text4,   # Unknown character set, 4 byte
     :UTF32]   # corresponding to codepoints (0-0xd7ff, 0xe000-0x10fff)
const _subsetnam =
    [:_Latin,
     :_UCS2,
     :_UTF32]
const _mbwname =
    [:UTF8,
     :UTF16]
const _binname = vcat(_cpname1, :Binary)

# List of basic character sets
const charsets = vcat(_cpname1, _cpname2, _cpname4)

for nam in vcat(charsets,
                :Binary,   # really, no character set at all, not text
                :UniPlus)  # valid Unicode, plus unknown characters (for String)
    charset = symstr(nam, "CharSet")
    @eval const $charset = CharSet{$(quotesym(nam))}
    @eval export $charset
end

# These are to indicate string types that must have at least one character of the type,
# for the internal types to make up the UniStr union type

const LatinSubSet  = CharSet{:LatinSubSet} # Has at least 1 character > 0x7f, all <= 0xff
const UCS2SubSet   = CharSet{:UCS2SubSet}  # Has at least 1 character > 0xff, all <= 0xffff
const UTF32SubSet  = CharSet{:UTF32SubSet} # Has at least 1 non-BMP character in string

export Native1Byte, UTF8Encoding
const Native1Byte  = Encoding(:Byte)
const UTF8Encoding = Encoding(:UTF8)
@eval show(io::IO, ::Type{UTF8Encoding})  = print(io, "UTF-8")
@eval show(io::IO, ::Type{Native1Byte}) = print(io, "8-bit")

# Allow handling different endian encodings

for (n, l, b, s) in (("2Byte", :LE2, :BE2, "16-bit"),
                     ("4Byte", :LE4, :BE4, "32-bit"),
                     ("UTF16", :UTF16LE, :UTF16BE, "UTF-16"))
    nat, swp = BIG_ENDIAN ? (b, l) : (l, b)
    natnam = symstr("Native",  n)
    swpnam = symstr("Swapped", n)
    @eval export $nat, $swp
    @eval const $natnam = Encoding($(quotesym("N", n)))
    @eval const $swpnam = Encoding($(quotesym("S", n)))
    @eval const $nat = $natnam
    @eval const $swp = $swpnam
    @eval show(io::IO, ::Type{$natnam}) = print(io, $s)
    @eval show(io::IO, ::Type{$swpnam}) = print(io, $(string(s, " ", BIG_ENDIAN ? "LE" : "BE")))
end

const _CSE{U} = Union{CharSet{U}, Encoding{U}} where {U}

print(io::IO, ::S) where {S<:_CSE{U}} where {U} =
    print(io, U)

show(io::IO, ::Type{CharSet{S}}) where {S}   = print(io, "CharSet{", string(S), "}")
show(io::IO, ::Type{Encoding{S}}) where {S}  = print(io, "Encoding{", string(S), "}")
show(io::IO, ::Type{CSE{CS,E}}) where {S,T,CS<:CharSet{S},E<:Encoding{T}} =
    print(io, "CSE{", string(S), ", ", string(T), "}")
print(io::IO, ::T) where {S,U,CS<:CharSet{S},E<:Encoding{U},T<:CSE{CS,E}} =
    (show(io, T); print(io, "()"))

const CodeUnitTypes = Union{UInt8, UInt16, UInt32}

# Note: this is still in transition to expressing character set, encoding
# and optional cached info for hashes, UTF-8/UTF-16 encodings, subsets, etc.
# via more type parameters

struct Str{T,SubStr,Cache,Hash} <: AbstractString
    data::String
    substr::SubStr
    cache::Cache
    hash::Hash

    (::Type{Str})(::Type{T}, v::String, s::S, c::C, h::H) where {T<:CSE,S,C,H} =
        new{T,S,C,H}(v,s,c,h)
end

(::Type{Str})(::Type{C}, v::String) where {C<:CSE} = Str(C, v, nothing, nothing, nothing)
(::Type{Str})(::Type{C}, v::Str) where {C<:CSE} = Str(C, v.data, nothing, nothing, nothing)

const MaybeSub{T} = Union{T, SubString{T}} where {T<:AbstractString}

struct Chr{CS<:CharSet,T<:CodeUnitTypes} <: AbstractChar
    v::T
    (::Type{Chr})(::Type{CS}, v::T) where {CS<:CharSet,T<:CodeUnitTypes} = new{CS,T}(v)
end

(::Type{<:Chr{CS,T}})(v::Number) where {CS<:CharSet,T<:CodeUnitTypes} = Chr(CS, T(v))

# Handle change from endof -> lastindex
@static if !isdefined(Base, :lastindex)
    lastindex(str::AbstractString) = Base.endof(str)
    lastindex(arr::AbstractArray) = Base.endof(arr)
    Base.endof(str::Str) = lastindex(str)
end
@static if !isdefined(Base, :firstindex)
    firstindex(str::AbstractString) = 1
    # AbstractVector might be an OffsetArray
    firstindex(str::Vector) = 1
end

# This needs to be redone, with character sets and the code unit as part of the type

for (names, typ) in ((_cpname1, UInt8), (_cpname2, UInt16), (_cpname4, UInt32)), nam in names
    chrnam = symstr(nam, "Chr")
    @eval const $chrnam = Chr{$(symstr(nam, "CharSet")), $typ}
    @eval export $chrnam
    @eval show(io::IO, ::Type{$chrnam}) = print(io, $(quotesym(chrnam)))
end

const _LatinChr = Chr{LatinSubSet, UInt8}
show(io::IO, ::Type{_LatinChr}) = print(io, :_LatinChr)

codepoint(ch::Chr) = ch.v
basetype(::Type{<:Chr{CS,B}}) where {CS,B} = B
charset(::Type{<:Chr{CS,B}}) where {CS,B} = CS
typemin(::Type{T}) where {CS,B,T<:Chr{CS,B}} = Chr(CS, typemin(B))
typemax(::Type{T}) where {CS,B,T<:Chr{CS,B}} = Chr(CS, typemax(B))

# Handle a few quirks
charset(::Type{<:AbstractChar}) = UTF32CharSet
charset(::Type{UInt8})          = BinaryCharSet  # UInt8 instead of "BinaryChr"
charset(::Type{Char})           = UniPlusCharSet # Char instead of "UniPlusChr"

const LatinChars   = Union{LatinChr, _LatinChr}
const ByteChars    = Union{ASCIIChr, LatinChr, _LatinChr, Text1Chr}
const WideChars    = Union{UCS2Chr, UTF32Chr}

# Definition of built-in CSEs (Character Set Encodings)

for (cs, enc) in ((Native1Byte, _binname), (Native2Byte, _cpname2), (Native4Byte, _cpname4)),
    nam in enc
    @eval const $(symstr(nam, "CSE")) = CSE{$(symstr(nam, "CharSet")), $cs}
end
const UTF8CSE   = CSE{UTF32CharSet, UTF8Encoding}
const UTF16CSE  = CSE{UTF32CharSet, NativeUTF16}

const _LatinCSE = CSE{LatinSubSet,  Native1Byte}
const _UCS2CSE  = CSE{UCS2SubSet,   Native2Byte}
const _UTF32CSE = CSE{UTF32SubSet,  Native4Byte}

export UniPlusCSE
const UniPlusCSE = CSE{UniPlusCharSet, UTF8Encoding}

for nam in vcat(charsets, :_Latin)
    @eval codepoint_cse(::Type{$(symstr(nam,"Chr"))}) = $(symstr(nam,"CSE"))
end
codepoint_cse(::Type{Char}) = UniPlusCSE

# Definition of built-in Str types

const empty_string = ""

for nam in vcat(_binname, _cpname2, _cpname4, _subsetnam, _mbwname)
    str = String(nam)
    sym = symstr(nam, "Str")
    cse = symstr(nam, "CSE")
    @eval const $sym = Str{$cse, Nothing, Nothing, Nothing}
    @eval (::Type{$sym})(v::Vector{UInt8}) = convert($sym, v)
    @eval show(io::IO, ::Type{$sym}) = print(io, $(quotesym(sym)))
    @eval show(io::IO, ::Type{$cse}) = print(io, $(quotesym(cse)))
    low = lowercase(str)
    if str[1] == '_'
        @eval empty_str(::Type{$cse}) = $(symstr("empty", low))
    else
        emp = symstr("empty_", low)
        @eval const $emp = Str($cse, empty_string)
        @eval empty_str(::Type{$cse}) = $emp
        @eval export $sym, $cse
    end
    @eval convert(::Type{$sym}, str::$sym) = str
end
empty_str(::Type{<:Str{C}}) where {C<:CSE} = empty_str(C)
empty_str(::Type{String}) = empty_string

typemin(::Type{T}) where {T<:Str} = empty_str(T)
typemin(::T) where {T<:Str} = empty_str(T)

"""Union type for fast dispatching"""
const UniStr = Union{ASCIIStr, _LatinStr, _UCS2Str, _UTF32Str}
show(io::IO, ::Type{UniStr}) = print(io, :UniStr)

# Display BinaryCSE as if String
show(io::IO, str::T) where {T<:Str{BinaryCSE}} = show(io, str.data)
show(io::IO, str::SubString{T}) where {T<:Str{BinaryCSE}} =
    @inbounds show(io, SubString(str.string.data, str.offset+1, str.offset+lastindex(str)))

_allocate(len) = Base._string_n((len+STR_KEEP_NUL-1)%Csize_t)

function _allocate(::Type{T}, len) where {T <: CodeUnitTypes}
    buf = _allocate((len+STR_KEEP_NUL-1) * sizeof(T))
    buf, reinterpret(Ptr{T}, pointer(buf))
end

# Various useful groups of character set types

# These should be done via traits
const Binary_CSEs   = Union{Text1CSE, BinaryCSE}
const Latin_CSEs    = Union{LatinCSE, _LatinCSE}
const UTF8_CSEs     = Union{UTF8CSE,  UniPlusCSE}
const UCS2_CSEs     = Union{UCS2CSE,  _UCS2CSE}
const UTF32_CSEs    = Union{UTF32CSE, _UTF32CSE}
const SubSet_CSEs   = Union{_LatinCSE, _UCS2CSE, _UTF32CSE}

const Byte_CSEs     = Union{Text1CSE, BinaryCSE, ASCIICSE, LatinCSE, _LatinCSE, UTF8CSE}
const Word_CSEs     = Union{Text2CSE, UCS2CSE, _UCS2CSE, UTF16CSE} # 16-bit characters
const Quad_CSEs     = Union{Text4CSE, UTF32CSE, _UTF32CSE}         # 32-bit code units

basecse(::Type{C}) where {C<:CSE} = C
basecse(::Type{_LatinCSE}) = LatinCSE
basecse(::Type{_UCS2CSE})  = UCS2CSE
basecse(::Type{_UTF32CSE}) = UTF32CSE

const MS_UTF8     = MaybeSub{<:Str{UTF8CSE}}
const MS_UTF16    = MaybeSub{<:Str{UTF16CSE}}
const MS_UTF32    = MaybeSub{<:Str{UTF32CSE}}
const MS_SubUTF32 = MaybeSub{<:Str{_UTF32CSE}}
const MS_ByteStr  = MaybeSub{<:Union{String,Str{BinaryCSE},Str{Text1CSE}}}

const AbsChar = @static isdefined(Base, :AbstractChar) ? AbstractChar : Union{Char, Chr}

## Get the character set / encoding used by a string type
cse(::Type{<:AbstractString}) = UniPlusCSE     # allows invalid sequences
cse(::Type{<:SubString{T}}) where {T} = cse(T)
cse(::Type{<:SubString{<:Str{C}}}) where {C<:SubSet_CSEs} = basecse(C)
cse(::Type{<:Str{C}}) where {C<:CSE} = C
cse(str::AbstractString) = cse(typeof(str))

charset(::Type{T}) where {T<:AbstractString} = charset(cse(T)) # Default unless overridden
charset(::Type{C}) where {CS,C<:CSE{CS}} = CS
charset(str::AbstractString) = charset(cse(str))

encoding(::Type{T}) where {T<:AbstractString} = encoding(cse(T)) # Default unless overridden
encoding(::Type{C}) where {CS,E,C<:CSE{CS,E}} = E
encoding(str::AbstractString) = encoding(cse(str))

promote_rule(::Type{T}, ::Type{T}) where {T<:Chr} = T
promote_rule(::Type{Text2Chr}, ::Type{Text1Chr}) = Text2Chr
promote_rule(::Type{Text4Chr}, ::Type{Text1Chr}) = Text4Chr
promote_rule(::Type{Text4Chr}, ::Type{Text2Chr}) = Text4Chr

promote_rule(::Type{T}, ::Type{ASCIIChr}) where {T} = T
promote_rule(::Type{LatinChr}, ::Type{_LatinChr}) = LatinChr
promote_rule(::Type{UTF32Chr}, ::Type{UCS2Chr}) = UTF32Chr
promote_rule(::Type{T}, ::Type{<:ByteChars}) where {T<:WideChars} = T

promote_rule(::Type{T}, ::Type{T}) where {T<:Str} = T

promote_rule(::Type{Text2CSE}, ::Type{Text1CSE}) = Text2CSE
promote_rule(::Type{Text4CSE}, ::Type{Text1CSE}) = Text4CSE
promote_rule(::Type{Text4CSE}, ::Type{Text2CSE}) = Text4CSE

promote_rule(::Type{T}, ::Type{ASCIICSE}) where {T<:CSE} = T
promote_rule(::Type{T},
             ::Type{<:Latin_CSEs}) where {T<:Union{UTF8CSE,UTF16CSE,UCS2_CSEs,UTF32_CSEs}} = T

promote_rule(::Type{T}, ::Type{_LatinCSE}) where {T<:Union{ASCIICSE,LatinCSE}} = LatinCSE

promote_rule(::Type{T}, ::Type{_UCS2CSE}) where {T<:Union{ASCIICSE,Latin_CSEs,UCS2CSE}} = UCS2CSE
promote_rule(::Type{T}, ::Type{_UTF32CSE}) where {T<:CSE} = UTF32CSE
promote_rule(::Type{T}, ::Type{UTF32CSE}) where {T<:UCS2_CSEs} = UTF32CSE

promote_rule(::Type{String}, ::Type{<:Str}) = String

promote_rule(::Type{<:Str{S}}, ::Type{<:Str{T}}) where {S,T} =
    (P = promote_rule(S,T)) === Union{} ? Union{} : Str{P}

sizeof(s::Str) = sizeof(s.data) + 1 - STR_KEEP_NUL

"""Codeunits of string as a Vector"""
_data(s::Vector{UInt8}) = s
_data(s::String)        = s
_data(s::Str{<:Byte_CSEs}) =
    @static V6_COMPAT ? Vector{UInt8}(s.data) : unsafe_wrap(Vector{UInt8}, s.data)

"""Pointer to codeunits of string"""
pointer(s::Str{<:Byte_CSEs}) = pointer(s.data)
pointer(s::Str{<:Word_CSEs}) = reinterpret(Ptr{UInt16}, pointer(s.data))
pointer(s::Str{<:Quad_CSEs}) = reinterpret(Ptr{UInt32}, pointer(s.data))

const CHUNKSZ = sizeof(UInt64) # used for fast processing of strings
const CHUNKMSK = (CHUNKSZ-1)%UInt64

_pnt64(s::Union{String,Vector{UInt8}}) = reinterpret(Ptr{UInt64}, pointer(s))
_pnt64(s::Str) = reinterpret(Ptr{UInt64}, pointer(s.data))

"""Length of string in codeunits"""
ncodeunits(s::Str)              = sizeof(s)
ncodeunits(s::Str{<:Word_CSEs}) = sizeof(s) >>> 1
ncodeunits(s::Str{<:Quad_CSEs}) = sizeof(s) >>> 2

# For convenience
@inline _calcpnt(str, siz) = (pnt = _pnt64(str) - CHUNKSZ;  (pnt, pnt + siz))

@inline _mask_bytes(n) = (1%UInt << ((n & CHUNKMSK) << 3)) - 0x1

Base.need_full_hex(c::Chr) = is_hex_digit(c)
Base.escape_nul(c::Chr) = ('0' <= c <= '7') ? "\\x00" : "\\0"

# Support for SubString of Str

Base.SubString(str::Str{C}) where {C<:SubSet_CSEs} =
    SubString(Str(basecse(C), str))
Base.SubString(str::Str{C}, off::Int) where {C<:SubSet_CSEs} =
    SubString(Str(basecse(C), str), off)
Base.SubString(str::Str{C}, off::Int, fin::Int) where {C<:SubSet_CSEs} =
    SubString(Str(basecse(C), str), off, fin)

# pointer conversions of ASCII/UTF8/UTF16/UTF32 strings:
pointer(str::Str, pos::Integer) = bytoff(pointer(str), pos - 1)

# pointer conversions of SubString of ASCII/UTF8/UTF16/UTF32:
pointer(x::SubString{<:Str}) = bytoff(pointer(x.string), x.offset)
pointer(x::SubString{<:Str}, pos::Integer) = bytoff(pointer(x.string), x.offset + pos - 1)

# CodePoints iterator

struct CodePoints{T}
    xs::T
end

"""
    codepoints(str)

An iterator that generates the code points of a string

# Examples
```jldoctest
julia> a = Str("abc\U1f596")

julia> collect(a)

julia> collect(codepoints(a))
```
"""
codepoints(xs) = CodePoints(xs)
eltype(::Type{<:CodePoints{S}}) where {S} = eltype(S)
length(it::CodePoints) = length(it.xs)
start(it::CodePoints) = start(it.xs)
done(it::CodePoints, state) = done(it.xs, state)
@propagate_inbounds next(it::CodePoints{T}, state) where {T<:Str} =
    _next(CodePointStyle(T), eltype(T), it.xs, state)


text:
# import Base:
#     VersionNumber
# VERSION = VersionNumber(0,0,1,("git",),(0005,))
# module version
# end

## julia version info

# Include build number if we've got at least some distance from a tag (e.g. a release)
try
    build_number = GIT_VERSION_INFO.build_number != 0 ? "+$(GIT_VERSION_INFO.build_number)" : ""
    global const VERSION = convert(VersionNumber, "$(VERSION_STRING)$(build_number)")
catch e
    println("while creating Base.VERSION, ignoring error $e")
    global const VERSION = VersionNumber(0)
end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using Printf
using JSON
using DataStructures
using GZip
import Base: getindex, time

const INSTANCES_URL = "https://axavier.org/UnitCommitment.jl/0.3/instances"

"""
    read_benchmark(name::AbstractString)::UnitCommitmentInstance

Read one of the benchmark unit commitment instances included in the package.
See "Instances" section of the documentation for the entire list of benchmark
instances available.

Example
-------

    import UnitCommitment
    instance = UnitCommitment.read_benchmark("matpower/case3375wp/2017-02-01")
"""
function read_benchmark(name::AbstractString)::UnitCommitmentInstance
    basedir = dirname(@__FILE__)
    filename = "$basedir/../../instances/$name.json.gz"
    url = "$INSTANCES_URL/$name.json.gz"
    if !isfile(filename)
        @info "Downloading: $(url)"
        dpath = download(url)
        mkpath(dirname(filename))
        cp(dpath, filename)
    end
    return UnitCommitment.read(filename)
end

"""
    read(path::AbstractString)::UnitCommitmentInstance

Read a unit commitment instance from a file. The file may be gzipped.

Example
-------

    import UnitCommitment
    instance = UnitCommitment.read("/path/to/input.json.gz")
"""
function read(path::AbstractString)::UnitCommitmentInstance
    if endswith(path, ".gz")
        return _read(gzopen(path))
    else
        return _read(open(path))
    end
end

function _read(file::IO)::UnitCommitmentInstance
    return _from_json(
        JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing)),
    )
end

function _read_json(path::String)::OrderedDict
    if endswith(path, ".gz")
        file = GZip.gzopen(path)
    else
        file = open(path)
    end
    return JSON.parse(file, dicttype = () -> DefaultOrderedDict(nothing))
end

function _from_json(json; repair = true)
    units = Unit[]
    buses = Bus[]
    contingencies = Contingency[]
    lines = TransmissionLine[]
    loads = PriceSensitiveLoad[]

    function scalar(x; default = nothing)
        x !== nothing || return default
        return x
    end

    time_horizon = json["Parameters"]["Time (h)"]
    if time_horizon === nothing
        time_horizon = json["Parameters"]["Time horizon (h)"]
    end
    time_horizon !== nothing || error("Missing parameter: Time horizon (h)")
    time_step = scalar(json["Parameters"]["Time step (min)"], default = 60)
    (60 % time_step == 0) ||
        error("Time step $time_step is not a divisor of 60")
    time_multiplier = 60 ÷ time_step
    T = time_horizon * time_multiplier

    name_to_bus = Dict{String,Bus}()
    name_to_line = Dict{String,TransmissionLine}()
    name_to_unit = Dict{String,Unit}()

    function timeseries(x; default = nothing)
        x !== nothing || return default
        x isa Array || return [x for t in 1:T]
        return x
    end

    # Read parameters
    power_balance_penalty = timeseries(
        json["Parameters"]["Power balance penalty (\$/MW)"],
        default = [1000.0 for t in 1:T],
    )
    shortfall_penalty = timeseries(
        json["Parameters"]["Reserve shortfall penalty (\$/MW)"],
        default = [-1.0 for t in 1:T],
    )

    # Read buses
    for (bus_name, dict) in json["Buses"]
        bus = Bus(
            bus_name,
            length(buses),
            timeseries(dict["Load (MW)"]),
            Unit[],
            PriceSensitiveLoad[],
        )
        name_to_bus[bus_name] = bus
        push!(buses, bus)
    end

    # Read units
    for (unit_name, dict) in json["Generators"]
        bus = name_to_bus[dict["Bus"]]

        # Read production cost curve
        K = length(dict["Production cost curve (MW)"])
        curve_mw = hcat(
            [timeseries(dict["Production cost curve (MW)"][k]) for k in 1:K]...,
        )
        curve_cost = hcat(
            [timeseries(dict["Production cost curve (\$)"][k]) for k in 1:K]...,
        )
        min_power = curve_mw[:, 1]
        max_power = curve_mw[:, K]
        min_power_cost = curve_cost[:, 1]
        segments = CostSegment[]
        for k in 2:K
            amount = curve_mw[:, k] - curve_mw[:, k-1]
            cost = (curve_cost[:, k] - curve_cost[:, k-1]) ./ amount
            replace!(cost, NaN => 0.0)
            push!(segments, CostSegment(amount, cost))
        end

        # Read startup costs
        startup_delays = scalar(dict["Startup delays (h)"], default = [1])
        startup_costs = scalar(dict["Startup costs (\$)"], default = [0.0])
        startup_categories = StartupCategory[]
        for k in 1:length(startup_delays)
            push!(
                startup_categories,
                StartupCategory(
                    startup_delays[k] .* time_multiplier,
                    startup_costs[k],
                ),
            )
        end

        # Read and validate initial conditions
        initial_power = scalar(dict["Initial power (MW)"], default = nothing)
        initial_status = scalar(dict["Initial status (h)"], default = nothing)
        if initial_power === nothing
            initial_status === nothing ||
                error("unit $unit_name has initial status but no initial power")
        else
            initial_status !== nothing ||
                error("unit $unit_name has initial power but no initial status")
            initial_status != 0 ||
                error("unit $unit_name has invalid initial status")
            if initial_status < 0 && initial_power > 1e-3
                error("unit $unit_name has invalid initial power")
            end
            initial_status *= time_multiplier
        end

        unit = Unit(
            unit_name,
            bus,
            max_power,
            min_power,
            timeseries(dict["Must run?"], default = [false for t in 1:T]),
            min_power_cost,
            segments,
            scalar(dict["Minimum uptime (h)"], default = 1) * time_multiplier,
            scalar(dict["Minimum downtime (h)"], default = 1) * time_multiplier,
            scalar(dict["Ramp up limit (MW)"], default = 1e6),
            scalar(dict["Ramp down limit (MW)"], default = 1e6),
            scalar(dict["Startup limit (MW)"], default = 1e6),
            scalar(dict["Shutdown limit (MW)"], default = 1e6),
            initial_status,
            initial_power,
            timeseries(
                dict["Provides spinning reserves?"],
                default = [true for t in 1:T],
            ),
            startup_categories,
        )
        push!(bus.units, unit)
        name_to_unit[unit_name] = unit
        push!(units, unit)
    end

    # Read reserves
    reserves = Reserves(zeros(T))
    if "Reserves" in keys(json)
        reserves.spinning =
            timeseries(json["Reserves"]["Spinning (MW)"], default = zeros(T))
    end

    # Read transmission lines
    if "Transmission lines" in keys(json)
        for (line_name, dict) in json["Transmission lines"]
            line = TransmissionLine(
                line_name,
                length(lines) + 1,
                name_to_bus[dict["Source bus"]],
                name_to_bus[dict["Target bus"]],
                scalar(dict["Reactance (ohms)"]),
                scalar(dict["Susceptance (S)"]),
                timeseries(
                    dict["Normal flow limit (MW)"],
                    default = [1e8 for t in 1:T],
                ),
                timeseries(
                    dict["Emergency flow limit (MW)"],
                    default = [1e8 for t in 1:T],
                ),
                timeseries(
                    dict["Flow limit penalty (\$/MW)"],
                    default = [5000.0 for t in 1:T],
                ),
            )
            name_to_line[line_name] = line
            push!(lines, line)
        end
    end

    # Read contingencies
    if "Contingencies" in keys(json)
        for (cont_name, dict) in json["Contingencies"]
            affected_units = Unit[]
            affected_lines = TransmissionLine[]
            if "Affected lines" in keys(dict)
                affected_lines =
                    [name_to_line[l] for l in dict["Affected lines"]]
            end
            if "Affected units" in keys(dict)
                affected_units =
                    [name_to_unit[u] for u in dict["Affected units"]]
            end
            cont = Contingency(cont_name, affected_lines, affected_units)
            push!(contingencies, cont)
        end
    end

    # Read price-sensitive loads
    if "Price-sensitive loads" in keys(json)
        for (load_name, dict) in json["Price-sensitive loads"]
            bus = name_to_bus[dict["Bus"]]
            load = PriceSensitiveLoad(
                load_name,
                bus,
                timeseries(dict["Demand (MW)"]),
                timeseries(dict["Revenue (\$/MW)"]),
            )
            push!(bus.price_sensitive_loads, load)
            push!(loads, load)
        end
    end

    instance = UnitCommitmentInstance(
        buses_by_name = Dict(b.name => b for b in buses),
        buses = buses,
        contingencies_by_name = Dict(c.name => c for c in contingencies),
        contingencies = contingencies,
        lines_by_name = Dict(l.name => l for l in lines),
        lines = lines,
        power_balance_penalty = power_balance_penalty,
        price_sensitive_loads_by_name = Dict(ps.name => ps for ps in loads),
        price_sensitive_loads = loads,
        reserves = reserves,
        shortfall_penalty = shortfall_penalty,
        time = T,
        units_by_name = Dict(g.name => g for g in units),
        units = units,
    )
    if repair
        UnitCommitment.repair!(instance)
    end
    return instance
end


text:
# This file was generated, do not modify it. # hide
fit!(clf)
ŷ = predict(clf, X2)
cross_entropy(ŷ, y) |> mean

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
# 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 = "SDPABuilder"
version = v"7.3.8"

# Collection of sources required to build SDPABuilder
sources = [
    "https://sourceforge.net/projects/sdpa/files/sdpa/sdpa_7.3.8.tar.gz" =>
    "c7541333da2f0bb2d18e90dbf758ac7cc099f3f7da3f256b284b0725f96d4117",
    "./bundled"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd sdpa-7.3.8/
update_configure_scripts

for path in ${LD_LIBRARY_PATH//:/ }; do
    for file in $(ls $path/*.la); do
        echo "$file"
        baddir=$(sed -n "s|libdir=||p" $file)
        sed -i~ -e "s|$baddir|'$path'|g" $file
    done
done
if [ $target = "x86_64-apple-darwin14" ]; then
  # seems static linking requires apple's ar
  export AR=/opt/x86_64-apple-darwin14/bin/x86_64-apple-darwin14-ar
fi

## First build SDPA

patch -p1 < $WORKSPACE/srcdir/patches/shared.diff
mv configure.in configure.ac
patch -p1 < $WORKSPACE/srcdir/patches/lt_init.diff
autoreconf -i

./configure --prefix=$prefix --host=${target}  lt_cv_deplibs_check_method=pass_all \
--with-blas="-L${prefix}/lib -lcoinblas -lgfortran -lcoinmumps -lcoinmetis" \
--with-lapack="-L${prefix}/lib -lcoinlapack -lcoinmumps -lcoinmetis" \
--with-mumps-libs="-L${prefix}/lib -lcoinmumps -lcoinmetis" --with-mumps-include="-I$prefix/include/coin/ThirdParty"

make
make install

## Then build the libcxxwrap-julia wrapper

cd $WORKSPACE/srcdir
cd sdpawrap

mkdir build
cd build
if [[ $target == i686-* ]] || [[ $target == arm-* ]]; then
    export processor=pentium4
else
    export processor=x86-64
fi

cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DSDPA_DIR=$prefix -DMUMPS_INCLUDE_DIR="$prefix/include/coin/ThirdParty" \
-DCMAKE_FIND_ROOT_PATH=$prefix -DJulia_PREFIX=$prefix  -DSDPA_LIBRARY="-lsdpa" -DCMAKE_CXX_FLAGS="-march=$processor" \
-D_GLIBCXX_USE_CXX11_ABI=1 ..
cmake --build . --config Release --target install

if [[ $target == *w64-mingw32* ]] ; then
    cp $WORKSPACE/destdir/lib/libsdpawrap.dll $WORKSPACE/destdir/bin
fi
"""

# platforms are restricted by libcxxwrap-julia, which requires gcc7 or gcc8
# and hence will not work with the official binaries for windows (which uses gcc4)

platforms = Platform[
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc7)),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc8)),
    Linux(:x86_64, compiler_abi=CompilerABI(:gcc7, :cxx11)),
    Linux(:x86_64, compiler_abi=CompilerABI(:gcc8, :cxx11)),
]

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "sdpa", :sdpa),
    LibraryProduct(prefix, "libsdpa", :libsdpa),
    LibraryProduct(prefix, "libsdpawrap", :libsdpawrap)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaOpt/COINBLASBuilder/releases/download/v1.4.6-1-static/build_COINBLASBuilder.v1.4.6.jl",
    "https://github.com/JuliaOpt/COINLapackBuilder/releases/download/v1.5.6-1-static/build_COINLapackBuilder.v1.5.6.jl",
    "https://github.com/JuliaOpt/COINMetisBuilder/releases/download/v1.3.5-1-static/build_COINMetisBuilder.v1.3.5.jl",
    "https://github.com/JuliaOpt/COINMumpsBuilder/releases/download/v1.6.0-1-static-nm/build_COINMumpsBuilder.v1.6.0.jl",
    "https://github.com/JuliaInterop/libcxxwrap-julia/releases/download/v0.5.1/build_libcxxwrap-julia-1.0.v0.5.1.jl",
    "https://github.com/JuliaPackaging/JuliaBuilder/releases/download/v1.0.0-2/build_Julia.v1.0.0.jl"
]

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


text:
module Xerxes

using Logging

mutable struct Program
    A::Matrix
    b::Vector
    c::Vector
    nonbasic_vars::Vector
    basic_vars::Vector
    x::Vector # current primal solution
    z::Vector # current dual solution
end

@enum Status Optimal Nonoptimal Unbounded Infeasible

function primal_iterate!(prog::Program)
    B = view(prog.A, :, prog.basic_vars)
    N = view(prog.A, :, prog.nonbasic_vars)
    @debug "Submatrices" B N

    # Check for optimality
    if all(prog.z[prog.nonbasic_vars] .>= 0)
        @info "Program is optimal"
        return Optimal
    end
    @debug "z has negative values, therefore the program is not optimal"

    # Select entering variable
    j = argmin(prog.z)
    if prog.z[j] >= 0
        @debug "No entering variable"
        return Infeasible
    end
    @debug "Entering variable" j

    # Compute primal step direction
    ej = zeros(length(prog.x))
    ej[j] += 1
    primal_step = B \ (N * ej[prog.nonbasic_vars])
    @debug "Primal step direction" primal_step

    # Compute primal step length
    inv_primal_step_length = -Inf
    i = 0
    for k in eachindex(prog.basic_vars)
        if primal_step[k] == 0 && prog.x[prog.basic_vars[k]] == 0
            t = 0
        else
            t = primal_step[k] / prog.x[prog.basic_vars[k]]
        end
        if inv_primal_step_length < t
            inv_primal_step_length = t
            i = prog.basic_vars[k]
        end
    end
    if inv_primal_step_length <= 0
        return Unbounded
    end
    primal_step_length = 1 / inv_primal_step_length
    @debug "Primal step length" primal_step_length
    primal_step *= primal_step_length

    # Select leaving variable (already done, i above)
    @debug "Leaving variable" i

    # Compute dual step direction
    ei = zeros(length(prog.z))
    ei[i] += 1
    dual_step = - N' * (B' \ ei[prog.basic_vars])
    @debug "Dual step direction" dual_step

    # Compute dual step length
    idx = indexin(j, prog.nonbasic_vars)[1]
    dual_step_length = prog.z[j] / dual_step[idx]
    @debug "Dual step length" dual_step_length
    dual_step *= dual_step_length

    # Update current dual and primal solutions
    prog.x[j] = primal_step_length
    prog.x[prog.basic_vars] = prog.x[prog.basic_vars] - primal_step
    prog.z[i] = dual_step_length
    prog.z[prog.nonbasic_vars] = prog.z[prog.nonbasic_vars] - dual_step

    # Update basis
    prog.basic_vars = union(setdiff(prog.basic_vars, [i]), [j])
    prog.nonbasic_vars = union(setdiff(prog.nonbasic_vars, [j]), [i])

    return Nonoptimal
end

end

max_stars_repo_name: dlozeve/Xerxes

max_issues_repo_name: dlozeve/Xerxes

max_forks_repo_name: dlozeve/Xerxes