text:
VERSION < v"0.7.0-beta2.199" && __precompile__()

module NLopt

export Opt, NLOPT_VERSION, algorithm, algorithm_name, ForcedStop,
       lower_bounds!, lower_bounds, upper_bounds!, upper_bounds, stopval!, stopval, ftol_rel!, ftol_rel, ftol_abs!, ftol_abs, xtol_rel!, xtol_rel, xtol_abs!, xtol_abs, maxeval!, maxeval, maxtime!, maxtime, force_stop!, force_stop, force_stop!, population!, population, vector_storage!, vector_storage, initial_step!, initial_step, default_initial_step!, local_optimizer!,
       min_objective!, max_objective!, equality_constraint!, inequality_constraint!, remove_constraints!,
       optimize!, optimize, Algorithm, Result

import MathProgBase.SolverInterface
import MathProgBase.SolverInterface.optimize!

import Libdl
const depsjl_path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("NLopt not installed properly; run Pkg.build(\"NLopt\"), restart Julia, and try again.")
end
include(depsjl_path)

function __init__()
    check_deps()
    v = version()
    v >= v"2.5" || error("NLopt $v < 2.5 is too old")
end

############################################################################
# Mirrors of NLopt's C enum constants:

@enum Algorithm::Cint begin
    GN_DIRECT=0
    GN_DIRECT_L=1
    GN_DIRECT_L_RAND=2
    GN_DIRECT_NOSCAL=3
    GN_DIRECT_L_NOSCAL=4
    GN_DIRECT_L_RAND_NOSCAL=5
    GN_ORIG_DIRECT=6
    GN_ORIG_DIRECT_L=7
    GD_STOGO=8
    GD_STOGO_RAND=9
    LD_LBFGS_NOCEDAL=10
    LD_LBFGS=11
    LN_PRAXIS=12
    LD_VAR1=13
    LD_VAR2=14
    LD_TNEWTON=15
    LD_TNEWTON_RESTART=16
    LD_TNEWTON_PRECOND=17
    LD_TNEWTON_PRECOND_RESTART=18
    GN_CRS2_LM=19
    GN_MLSL=20
    GD_MLSL=21
    GN_MLSL_LDS=22
    GD_MLSL_LDS=23
    LD_MMA=24
    LN_COBYLA=25
    LN_NEWUOA=26
    LN_NEWUOA_BOUND=27
    LN_NELDERMEAD=28
    LN_SBPLX=29
    LN_AUGLAG=30
    LD_AUGLAG=31
    LN_AUGLAG_EQ=32
    LD_AUGLAG_EQ=33
    LN_BOBYQA=34
    GN_ISRES=35
    AUGLAG=36
    AUGLAG_EQ=37
    G_MLSL=38
    G_MLSL_LDS=39
    LD_SLSQP=40
    LD_CCSAQ=41
    GN_ESCH=42
    GN_AGS=43
end

const sym2alg = Dict(Symbol(i)=>i for i in instances(Algorithm))

# enum nlopt_result
@enum Result::Cint begin
    FORCED_STOP=-5
    ROUNDOFF_LIMITED=-4
    OUT_OF_MEMORY=-3
    INVALID_ARGS=-2
    FAILURE=-1
    SUCCESS=1
    STOPVAL_REACHED=2
    FTOL_REACHED=3
    XTOL_REACHED=4
    MAXEVAL_REACHED=5
    MAXTIME_REACHED=6
end

# so that result < 0 checks continue to work
Base.isless(x::Integer, r::Result) = isless(x, Cint(r))
Base.isless(r::Result, x::Integer) = isless(Cint(r), x)
# so that == :Foo checks continue to work
Base.:(==)(s::Symbol, r::Result) = s == Symbol(r)
Base.:(==)(r::Result, s::Symbol) = s == r

############################################################################
# wrapper around nlopt_opt type

const _Opt = Ptr{Cvoid} # nlopt_opt

# pass both f and o to the callback so that we can handle exceptions
struct Callback_Data
    f::Function
    o::Any # should be Opt, but see Julia issue #269
end

mutable struct Opt
    opt::_Opt

    # need to store callback data for objective and constraints in
    # Opt so that they aren't garbage-collected.  cb[1] is the objective.
    cb::Vector{Callback_Data}

    function Opt(p::_Opt)
        opt = new(p, Array{Callback_Data}(undef,1))
        finalizer(destroy,opt)
        opt
    end
    function Opt(algorithm::Algorithm, n::Integer)
        n >= 0 || throw(ArgumentError("invalid dimension $n < 0"))
        p = ccall((:nlopt_create,libnlopt), _Opt, (Algorithm, Cuint),
                  algorithm, n)
        if p == C_NULL
            error("Error in nlopt_create")
        end
        Opt(p)
    end
    Opt(alg::Integer, n::Integer) = Opt(Algorithm(alg), n)
    Opt(algorithm::Symbol, n::Integer) = Opt(try sym2alg[algorithm]
                                             catch
                         throw(ArgumentError("unknown algorithm $algorithm"))
                                             end, n)
end

Base.unsafe_convert(::Type{_Opt}, o::Opt) = getfield(o, :opt) # for passing to ccall

destroy(o::Opt) = ccall((:nlopt_destroy,libnlopt), Cvoid, (_Opt,), o)

Base.ndims(o::Opt) = Int(ccall((:nlopt_get_dimension,libnlopt), Cuint, (_Opt,), o))
algorithm(o::Opt) = ccall((:nlopt_get_algorithm,libnlopt), Algorithm, (_Opt,), o)

Base.show(io::IO, o::Opt) = print(io, "Opt($(algorithm(o)), $(ndims(o)))")

############################################################################
# copying is a little tricky because we have to tell NLopt to use
# new Callback_Data.

# callback wrapper for nlopt_munge_data in NLopt 2.4
function munge_callback(p::Ptr{Cvoid}, f_::Ptr{Cvoid})
    f = unsafe_pointer_to_objref(f_)::Function
    f(p)::Ptr{Cvoid}
end

function Base.copy(o::Opt)
    p = ccall((:nlopt_copy,libnlopt), _Opt, (_Opt,), o)
    if p == C_NULL
        error("Error in nlopt_copy")
    end
    n = Opt(p)

    cb = getfield(o.cb)
    n.cb = similar(cb)
    for i = 1:length(cb)
        try
            n.cb[i] = Callback_Data(cb[i].f, n)
        catch e
            # if objective has not been set, cb[1] will throw
            # an UndefRefError, which is okay.
            if i != 1 || !isa(e, UndefRefError)
                rethrow(e) # some not-okay exception
            end
        end
    end

    try
        # n.o, for each callback, stores a pointer to an element of o.cb,
        # and we need to convert this into a pointer to the corresponding
        # element of n.cb.  nlopt_munge_data allows us to call a function
        # to transform each stored pointer in n.o, and we use the cbi
        # dictionary to convert pointers to indices into o.cb, whence
        # we obtain the corresponding element of n.cb.
        cbi = Dict{Ptr{Cvoid},Int}()
        for i in 1:length(cb)
            try
                cbi[pointer_from_objref(cb[i])] = i
            catch
            end
        end
        munge_callback_ptr = @cfunction(munge_callback, Ptr{Cvoid},
                                        (Ptr{Cvoid}, Ptr{Cvoid}))
        ccall((:nlopt_munge_data,libnlopt), Cvoid, (_Opt, Ptr{Cvoid}, Any),
              n, munge_callback_ptr,
              p::Ptr{Cvoid} -> p==C_NULL ? C_NULL :
                              pointer_from_objref(n.cb[cbi[p]]))
    catch e0
        # nlopt_munge_data not available, punt unless there is
        # no callback data
        try
            cb[1]
        catch e
            if length(cb) == 1 && isa(e, UndefRefError)
                return n
            end
        end
        error("copy(o::Opt) not supported for NLopt version < 2.4")
    end

    return n
end

############################################################################
# converting error results into exceptions

struct ForcedStop <: Exception end

# cache current exception for forced stop
nlopt_exception = nothing

function errmsg(o::Opt)
    msg = ccall((:nlopt_get_errmsg,libnlopt), Ptr{UInt8}, (_Opt,), o)
    return msg == C_NULL ? nothing : unsafe_string(msg)
end

function _errmsg(o::Opt)
    s = errmsg(o)
    return s === nothing || isempty(s) ? "" : ": "*s
end

# check result and throw an exception if necessary
function chk(o::Opt, result::Result)
    if result < 0 && result != ROUNDOFF_LIMITED
        if result == INVALID_ARGS
            throw(ArgumentError("invalid NLopt arguments"*_errmsg(o)))
        elseif result == OUT_OF_MEMORY
            throw(OutOfMemoryError())
        elseif result == FORCED_STOP
            global nlopt_exception
            e = nlopt_exception
            if e !== nothing && !isa(e, ForcedStop)
                nlopt_exception = nothing
                rethrow(e)
            end
        else
            error("nlopt failure $result", _errmsg(o))
        end
    end
    return nothing
end

############################################################################
# getting and setting scalar and vector parameters

# make a quoted symbol expression out of the arguments
qsym(args...) = Expr(:quote, Symbol(string(args...)))

# scalar parameters p of type T
macro GETSET(T, p)
    Tg = T == :Cdouble ? :Real : (T == :Cint || T == :Cuint ? :Integer : :Any)
    ps = Symbol(string(p, "!"))
    quote
        $(esc(p))(o::Opt) = ccall(($(qsym("nlopt_get_", p)),libnlopt),
                                  $T, (_Opt,), o)
        $(esc(ps))(o::Opt, val::$Tg) =
          chk(o, ccall(($(qsym("nlopt_set_", p)),libnlopt),
                     Result, (_Opt, $T), o, val))
    end
end

# Vector{Cdouble} parameters p
macro GETSET_VEC(p)
    ps = Symbol(string(p, "!"))
    quote
        function $(esc(p))(o::Opt, v::Vector{Cdouble})
            if length(v) != ndims(o)
                throw(BoundsError())
            end
            chk(o, ccall(($(qsym("nlopt_get_", p)),libnlopt),
                      Result, (_Opt, Ptr{Cdouble}), o, v))
            v
        end
        $(esc(p))(o::Opt) = $(esc(p))(o, Array{Cdouble}(undef, ndims(o)))
        function $(esc(ps))(o::Opt, v::Vector{Cdouble})
            if length(v) != ndims(o)
                throw(BoundsError())
            end
            chk(o, ccall(($(qsym("nlopt_set_", p)),libnlopt),
                      Result, (_Opt, Ptr{Cdouble}), o, v))
        end
        $(esc(ps))(o::Opt, v::AbstractVector{<:Real}) =
          $(esc(ps))(o, Array{Cdouble}(v))
        $(esc(ps))(o::Opt, val::Real) =
          chk(o, ccall(($(qsym("nlopt_set_", p, "1")),libnlopt),
                     Result, (_Opt, Cdouble), o, val))
    end
end

############################################################################
# Optimizer parameters

@GETSET_VEC lower_bounds
@GETSET_VEC upper_bounds
@GETSET Cdouble stopval
@GETSET Cdouble ftol_rel
@GETSET Cdouble ftol_abs
@GETSET Cdouble xtol_rel
@GETSET_VEC xtol_abs
@GETSET Cint maxeval
@GETSET Cdouble maxtime
@GETSET Cint force_stop
@GETSET Cuint population
@GETSET Cuint vector_storage

force_stop!(o::Opt) = force_stop!(o, 1)

local_optimizer!(o::Opt, lo::Opt) =
  chk(o, ccall((:nlopt_set_local_optimizer,libnlopt),
             Result, (_Opt, _Opt), o, lo))

# the initial-stepsize stuff is a bit different than GETSET_VEC,
# since the heuristics depend on the position x.

function default_initial_step!(o::Opt, x::Vector{Cdouble})
    if length(x) != ndims(o)
        throw(BoundsError())
    end
    chk(o, ccall((:nlopt_set_default_initial_step,libnlopt),
               Result, (_Opt, Ptr{Cdouble}), o, x))
end
default_initial_step!(o::Opt, x::AbstractVector{<:Real}) =
  default_initial_step!(o, Array{Cdouble}(x))

function initial_step!(o::Opt, dx::Vector{Cdouble})
    if length(dx) != ndims(o)
        throw(BoundsError())
    end
    chk(o, ccall((:nlopt_set_initial_step,libnlopt),
               Result, (_Opt, Ptr{Cdouble}), o, dx))
end
initial_step!(o::Opt, dx::AbstractVector{<:Real}) =
  initial_step!(o, Array{Cdouble}(dx))
initial_step!(o::Opt, dx::Real) =
  chk(o, ccall((:nlopt_set_initial_step1,libnlopt),
             Result, (_Opt, Cdouble), o, dx))

function initial_step(o::Opt, x::Vector{Cdouble}, dx::Vector{Cdouble})
    if length(x) != ndims(o) || length(dx) != ndims(o)
        throw(BoundsError())
    end
    chk(o, ccall((:nlopt_get_initial_step,libnlopt),
               Result, (_Opt, Ptr{Cdouble}, Ptr{Cdouble}), o, x, dx))
    dx
end
initial_step(o::Opt, x::AbstractVector{<:Real}) =
    initial_step(o, Array{Cdouble}(x),
                 Array{Cdouble}(undef, ndims(o)))

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

function algorithm_name(a::Algorithm)
    s = ccall((:nlopt_algorithm_name,libnlopt), Ptr{UInt8}, (Algorithm,), a)
    if s == C_NULL
        throw(ArgumentError("invalid algorithm $a"))
    end
    return unsafe_string(s)
end

algorithm_name(a::Integer) = algorithm_name(Algorithm(a))
algorithm_name(a::Symbol) = algorithm_name(try sym2alg[a]
                                           catch
                             throw(ArgumentError("unknown algorithm $a"))
                                           end)
algorithm_name(o::Opt) = algorithm_name(algorithm(o))

function Base.show(io::IO, ::MIME"text/plain", a::Algorithm)
    show(io, a)
    print(io, ": ", algorithm_name(a))
end

numevals(o::Opt) = ccall((:nlopt_get_numevals,libnlopt), Cint, (_Opt,), o)

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

function version()
    major = Ref{Cint}()
    minor = Ref{Cint}()
    patch = Ref{Cint}()
    ccall((:nlopt_version,libnlopt), Cvoid, (Ref{Cint},Ref{Cint},Ref{Cint}),
          major, minor, patch)
    return VersionNumber(major[], minor[], patch[])
end

const NLOPT_VERSION = version()

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

srand(seed::Integer) = ccall((:nlopt_srand,libnlopt),
                             Cvoid, (Culong,), seed)
srand_time() = ccall((:nlopt_srand_time,libnlopt), Cvoid, ())

############################################################################
# Objective function:

const empty_grad = Cdouble[] # for passing when grad == C_NULL

function nlopt_callback_wrapper(n::Cuint, x::Ptr{Cdouble},
                                grad::Ptr{Cdouble}, d_::Ptr{Cvoid})
    d = unsafe_pointer_to_objref(d_)::Callback_Data
    try
        res = convert(Cdouble,
                      d.f(unsafe_wrap(Array, x, (convert(Int, n),)),
                          grad == C_NULL ? empty_grad
                          : unsafe_wrap(Array, grad, (convert(Int, n),))))
        return res::Cdouble
    catch e
        global nlopt_exception
        nlopt_exception = e
        force_stop!(d.o::Opt)
        return 0.0 # ignored by nlopt
    end
end

for m in (:min, :max)
    mf = Symbol(string(m,"_objective!"))
    @eval function $mf(o::Opt, f::Function)
        getfield(o, :cb)[1] = Callback_Data(f, o)
        nlopt_callback_wrapper_ptr = @cfunction(nlopt_callback_wrapper,
            Cdouble, (Cuint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
        chk(o, ccall(($(qsym("nlopt_set_", m, "_objective")),libnlopt),
                   Result, (_Opt, Ptr{Cvoid}, Any),
                   o, nlopt_callback_wrapper_ptr,
                   getfield(o, :cb)[1]))
    end
end

############################################################################
# Nonlinear constraints:

for c in (:inequality, :equality)
    cf = Symbol(string(c, "_constraint!"))
    @eval function $cf(o::Opt, f::Function, tol::Real=0.0)
        push!(getfield(o, :cb), Callback_Data(f, o))
        nlopt_callback_wrapper_ptr = @cfunction(nlopt_callback_wrapper,
            Cdouble, (Cuint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
        chk(o, ccall(($(qsym("nlopt_add_", c, "_constraint")),libnlopt),
                   Result, (_Opt, Ptr{Cvoid}, Any, Cdouble),
                   o, nlopt_callback_wrapper_ptr,
                   getfield(o, :cb)[end], tol))
    end
end

function remove_constraints!(o::Opt)
    resize!(getfield(o, :cb), 1)
    chkn(ccall((:nlopt_remove_inequality_constraints,libnlopt),
               Result, (_Opt,), o))
    chkn(ccall((:nlopt_remove_equality_constraints,libnlopt),
               Result, (_Opt,), o))
end

############################################################################
# Vector-valued constraints

const empty_jac = Array{Cdouble}(undef,0,0) # for passing when grad == C_NULL

function nlopt_vcallback_wrapper(m::Cuint, res::Ptr{Cdouble},
                                 n::Cuint, x::Ptr{Cdouble},
                                 grad::Ptr{Cdouble}, d_::Ptr{Cvoid})
    d = unsafe_pointer_to_objref(d_)::Callback_Data
    try
        d.f(unsafe_wrap(Array, res, (convert(Int, m),)),
            unsafe_wrap(Array, x, (convert(Int, n),)),
            grad == C_NULL ? empty_jac
            : unsafe_wrap(Array, grad, (convert(Int, n),convert(Int, m))))
    catch e
        global nlopt_exception
        nlopt_exception = e
        force_stop!(d.o::Opt)
    end
    nothing
end

for c in (:inequality, :equality)
    cf = Symbol(string(c, "_constraint!"))
    @eval begin
        function $cf(o::Opt, f::Function, tol::Vector{Cdouble})
            push!(getfield(o, :cb), Callback_Data(f, o))
            nlopt_vcallback_wrapper_ptr = @cfunction(nlopt_vcallback_wrapper, Cvoid,
                  (Cuint, Ptr{Cdouble}, Cuint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
            chk(o, ccall(($(qsym("nlopt_add_", c, "_mconstraint")),
                        libnlopt),
                       Result, (_Opt, Cuint, Ptr{Cvoid}, Any, Ptr{Cdouble}),
                       o, length(tol), nlopt_vcallback_wrapper_ptr,
                       getfield(o, :cb)[end], tol))
        end
        $cf(o::Opt, f::Function, tol::AbstractVector{<:Real}) =
           $cf(o, f, Array{Float64}(tol))
        $cf(o::Opt, m::Integer, f::Function, tol::Real=0.0) =
           $cf(o, f, fill!(Cdouble(tol), m))
    end
end

############################################################################
# property-based getters setters opt.foo for Julia 0.7
# … at some point we will deprecate the old interface.

function Base.getproperty(o::Opt, p::Symbol)
    if p === :lower_bounds
        return lower_bounds(o)
    elseif p === :upper_bounds
        return upper_bounds(o)
    elseif p === :stopval
        return stopval(o)
    elseif p === :ftol_rel
        return ftol_rel(o)
    elseif p === :ftol_abs
        return ftol_abs(o)
    elseif p === :xtol_rel
        return xtol_rel(o)
    elseif p === :xtol_abs
        return xtol_abs(o)
    elseif p === :maxeval
        return maxeval(o)
    elseif p === :maxtime
        return maxtime(o)
    elseif p === :force_stop
        return force_stop(o)
    elseif p === :population
        return population(o)
    elseif p === :vector_storage
        return vector_storage(o)
    elseif p === :initial_step
        return initial_step(o)
    elseif p === :algorithm
        return algorithm(o)
    elseif p === :numevals
        return numevals(o)
    elseif p === :errmsg
        return errmsg(o)
    else
        error("type Opt has no readable property $p")
    end
end

function Base.setproperty!(o::Opt, p::Symbol, x)
    if p === :lower_bounds
        lower_bounds!(o, x)
    elseif p === :upper_bounds
        upper_bounds!(o, x)
    elseif p === :stopval
        stopval!(o, x)
    elseif p === :ftol_rel
        ftol_rel!(o, x)
    elseif p === :ftol_abs
        ftol_abs!(o, x)
    elseif p === :xtol_rel
        xtol_rel!(o, x)
    elseif p === :xtol_abs
        xtol_abs!(o, x)
    elseif p === :maxeval
        maxeval!(o, x)
    elseif p === :maxtime
        maxtime!(o, x)
    elseif p === :force_stop
        force_stop!(o, x)
    elseif p === :population
        population!(o, x)
    elseif p === :vector_storage
        vector_storage!(o, x)
    elseif p === :local_optimizer
        local_optimizer!(o, x)
    elseif p === :default_initial_step
        default_initial_step!(o, x)
    elseif p === :initial_step
        initial_step!(o, x)
    elseif p === :min_objective
        min_objective!(o, x)
    elseif p === :max_objective
        max_objective!(o, x)
    elseif p === :inequality_constraint
        inequality_constraint!(o, x)
    elseif p === :equality_constraint
        equality_constraint!(o, x)
    else
        error("type Opt has no writable property $p")
    end
    return x
end

Base.propertynames(o::Opt) =
   (:lower_bounds, :upper_bounds, :stopval, :ftol_rel, :ftol_abs, :xtol_rel, :xtol_abs, :maxeval, :maxtime, :force_stop, :population, :vector_storage, :initial_step, :algorithm, :local_optimizer, :default_initial_step, :initial_step, :min_objective, :max_objective, :inequality_constraint, :equality_constraint, :numevals, :errmsg)

############################################################################
# Perform the optimization:

function optimize!(o::Opt, x::Vector{Cdouble})
    if length(x) != ndims(o)
        throw(BoundsError())
    end
    opt_f = Array{Cdouble}(undef,1)
    ret = ccall((:nlopt_optimize,libnlopt), Result, (_Opt, Ptr{Cdouble},
                                                     Ptr{Cdouble}),
                o, x, opt_f)
    ret == INVALID_ARGS && chk(o, ret)
    return (opt_f[1], x, Symbol(ret))
end

optimize(o::Opt, x::AbstractVector{<:Real}) =
  optimize!(o, copyto!(Array{Cdouble}(undef,length(x)), x))

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

include("NLoptSolverInterface.jl")

end # module


text:
export LPIPLinearProblem

mutable struct LPIPVariables{T}
    x::Vector{T}
    s::Vector{T}
    p::Vector{T}

    function LPIPVariables{T}(m::Int, n::Int) where T
        return new{T}(
                    ones(T, n),
                    ones(T, n),
                    ones(T, m)
                )
    end
end

struct LPIPLinearProblem{T}
    A::AbstractMatrix{T}
    b::AbstractVector{T}
    c::AbstractVector{T}
    d::T
    variables::LPIPVariables{T}
    m::Int # Number of constraints
    n::Int # Number of variables
    f::Int # Number of equality constraints
    l::Int # Number of nonnegative constraints

    function LPIPLinearProblem{T}(A::AbstractMatrix{T}, b::AbstractVector{T}, c::AbstractVector{T},
                                  d::T, l::Int) where T

        # The number of equality constraints must be
        # number of constraints - number of nonnegative
        f = l > size(A, 1) ? error("The number of nonnegative constraints must be smaller than the number of rows") : size(A, 1) - l
        extended_A = extend_A(A, f, l)
        extended_c = extend_c(c, l)
        m, n = size(extended_A)
        return new{T}(
                    extended_A, b,
                    extended_c, d,
                    LPIPVariables{T}(m, n),
                    m, n, f, l
               )
    end
end

function extend_A(A::AbstractMatrix{T}, f::Int, l::Int) where T
    if l > 0 && f != 0
        return [A [zeros(f, l); Matrix{T}(I, l, l)]]
    elseif l > 0 && f == 0
        return [A Matrix{T}(I, l, l)]
    end
    return A
end

function extend_c(c::AbstractVector{T}, l::Int) where T
    if l > 0
        return [c; zeros(T, l)]
    end
    return c
end

mutable struct Result{T}
    variables::LPIPVariables{T}
    obj_val::T
    dual_obj_val::T
    status::Int
    iter::Int
    time::Float64

    function Result{T}(vars::LPIPVariables{T}, obj_val::T, dual_obj_val::T,
                       status::Int, iter::Int, time::Float64) where T
        return new(vars, obj_val, dual_obj_val, status, iter, time)
    end
end

function Result(lpip_pb::LPIPLinearProblem{T}, status::Int, iter::Int, time::Float64) where T
    obj_val = dot(lpip_pb.c, lpip_pb.variables.x) + lpip_pb.d
    dual_obj_val = dot(lpip_pb.b, lpip_pb.variables.p) + lpip_pb.d
    return Result{T}(lpip_pb.variables, obj_val, dual_obj_val, status, iter, time)
end

mutable struct Params
    verbose::Bool
    max_iter::Int
    rho::Float64
    alpha::Float64
    tol::Float64
    time_limit::Float64
    infeasible_unbounded_tol::Float64

    function Params(;verbose::Bool = false, max_iter::Int = 100, rho::Float64 = 0.1,
                    alpha::Float64 = 0.9, tol::Float64 = 1e-6, time_limit::Float64 = 1e10,
                    infeasible_unbounded_tol::Float64 = 1e20)
        (rho <= 0 && rho > 1) && error("rho must be between in [0, 1)")
        (alpha < 0 && alpha > 1) && error("alpha must be in (0, 1)")
        return new(verbose, max_iter, rho, alpha, tol, time_limit, infeasible_unbounded_tol)
    end
end
function Params(args)
    params = Params()
    parse_args!(params, args)
    return params
end
function parse_args!(params::Params, args)
    for i in args
        parse_arg!(params::Params, i)
    end
    return nothing
end
function parse_arg!(params::Params, arg)
    fields = fieldnames(Params)
    name = arg[1]
    value = arg[2]
    if name in fields
        setfield!(params, name, value)
    end
    return nothing
end

ext: jl


text:
using InteractiveFormatter
using Test

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


text:
# External packages.
using .Convex, .SCS

import .WiringDiagramLayouts: has_port_layout_method, solve_isotonic

has_port_layout_method(::Val{:isotonic}) = true

function solve_isotonic(y::Vector; solver=()->SCS.Optimizer(verbose=false),
                        loss=sumsquares, lower::Number=-Inf, upper::Number=Inf,
                        pad::Number=0)
  if isempty(y)
    return empty(y)
  end
  x = Variable(length(y))
  constraints = Constraint[ x[i] - x[i-1] >= pad for i in 2:length(x) ]
  if lower != -Inf
    push!(constraints, x[1] >= lower)
  end
  if upper != Inf
    push!(constraints, x[end] <= upper)
  end
  solve!(minimize(loss(x-y), constraints), solver)
  length(x) == 1 ? [ x.value ] : dropdims(x.value, dims=2) # XXX: MATLAB-ism.
end


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

# Dates.DateTime arithmetic
a = Dates.DateTime(2013,1,1,0,0,0,1)
b = Dates.DateTime(2013,1,1,0,0,0,0)
@test a - b == Dates.Millisecond(1)
@test Dates.DateTime(2013,1,2) - b == Dates.Millisecond(86400000)

# Dates.DateTime-Year arithmetic
dt = Dates.DateTime(1999,12,27)
@test dt + Dates.Year(1) == Dates.DateTime(2000,12,27)
@test dt + Dates.Year(100) == Dates.DateTime(2099,12,27)
@test dt + Dates.Year(1000) == Dates.DateTime(2999,12,27)
@test dt - Dates.Year(1) == Dates.DateTime(1998,12,27)
@test dt - Dates.Year(100) == Dates.DateTime(1899,12,27)
@test dt - Dates.Year(1000) == Dates.DateTime(999,12,27)
dt = Dates.DateTime(2000,2,29)
@test dt + Dates.Year(1) == Dates.DateTime(2001,2,28)
@test dt - Dates.Year(1) == Dates.DateTime(1999,2,28)
@test dt + Dates.Year(4) == Dates.DateTime(2004,2,29)
@test dt - Dates.Year(4) == Dates.DateTime(1996,2,29)
dt = Dates.DateTime(1972,6,30,23,59,59)
@test dt + Dates.Year(1) == Dates.DateTime(1973,6,30,23,59,59)
@test dt - Dates.Year(1) == Dates.DateTime(1971,6,30,23,59,59)
@test dt + Dates.Year(-1) == Dates.DateTime(1971,6,30,23,59,59)
@test dt - Dates.Year(-1) == Dates.DateTime(1973,6,30,23,59,59)

# Wrapping arithemtic for Months
# This ends up being trickier than expected because
# the user might do 2014-01-01 + Month(-14)
# monthwrap figures out the resulting month
# when adding/subtracting months from a date
@test Dates.monthwrap(1,-14) == 11
@test Dates.monthwrap(1,-13) == 12
@test Dates.monthwrap(1,-12) == 1
@test Dates.monthwrap(1,-11) == 2
@test Dates.monthwrap(1,-10) == 3
@test Dates.monthwrap(1,-9) == 4
@test Dates.monthwrap(1,-8) == 5
@test Dates.monthwrap(1,-7) == 6
@test Dates.monthwrap(1,-6) == 7
@test Dates.monthwrap(1,-5) == 8
@test Dates.monthwrap(1,-4) == 9
@test Dates.monthwrap(1,-3) == 10
@test Dates.monthwrap(1,-2) == 11
@test Dates.monthwrap(1,-1) == 12
@test Dates.monthwrap(1,0) == 1
@test Dates.monthwrap(1,1) == 2
@test Dates.monthwrap(1,2) == 3
@test Dates.monthwrap(1,3) == 4
@test Dates.monthwrap(1,4) == 5
@test Dates.monthwrap(1,5) == 6
@test Dates.monthwrap(1,6) == 7
@test Dates.monthwrap(1,7) == 8
@test Dates.monthwrap(1,8) == 9
@test Dates.monthwrap(1,9) == 10
@test Dates.monthwrap(1,10) == 11
@test Dates.monthwrap(1,11) == 12
@test Dates.monthwrap(1,12) == 1
@test Dates.monthwrap(1,13) == 2
@test Dates.monthwrap(1,24) == 1
@test Dates.monthwrap(12,-14) == 10
@test Dates.monthwrap(12,-13) == 11
@test Dates.monthwrap(12,-12) == 12
@test Dates.monthwrap(12,-11) == 1
@test Dates.monthwrap(12,-2) == 10
@test Dates.monthwrap(12,-1) == 11
@test Dates.monthwrap(12,0) == 12
@test Dates.monthwrap(12,1) == 1
@test Dates.monthwrap(12,2) == 2
@test Dates.monthwrap(12,11) == 11
@test Dates.monthwrap(12,12) == 12
@test Dates.monthwrap(12,13) == 1

# yearwrap figures out the resulting year
# when adding/subtracting months from a date
@test Dates.yearwrap(2000,1,-3600) == 1700
@test Dates.yearwrap(2000,1,-37) == 1996
@test Dates.yearwrap(2000,1,-36) == 1997
@test Dates.yearwrap(2000,1,-35) == 1997
@test Dates.yearwrap(2000,1,-25) == 1997
@test Dates.yearwrap(2000,1,-24) == 1998
@test Dates.yearwrap(2000,1,-23) == 1998
@test Dates.yearwrap(2000,1,-14) == 1998
@test Dates.yearwrap(2000,1,-13) == 1998
@test Dates.yearwrap(2000,1,-12) == 1999
@test Dates.yearwrap(2000,1,-11) == 1999
@test Dates.yearwrap(2000,1,-2) == 1999
@test Dates.yearwrap(2000,1,-1) == 1999
@test Dates.yearwrap(2000,1,0) == 2000
@test Dates.yearwrap(2000,1,1) == 2000
@test Dates.yearwrap(2000,1,11) == 2000
@test Dates.yearwrap(2000,1,12) == 2001
@test Dates.yearwrap(2000,1,13) == 2001
@test Dates.yearwrap(2000,1,23) == 2001
@test Dates.yearwrap(2000,1,24) == 2002
@test Dates.yearwrap(2000,1,25) == 2002
@test Dates.yearwrap(2000,1,36) == 2003
@test Dates.yearwrap(2000,1,3600) == 2300
@test Dates.yearwrap(2000,2,-2) == 1999
@test Dates.yearwrap(2000,3,10) == 2001
@test Dates.yearwrap(2000,4,-4) == 1999
@test Dates.yearwrap(2000,5,8) == 2001
@test Dates.yearwrap(2000,6,-18) == 1998
@test Dates.yearwrap(2000,6,-6) == 1999
@test Dates.yearwrap(2000,6,6) == 2000
@test Dates.yearwrap(2000,6,7) == 2001
@test Dates.yearwrap(2000,6,19) == 2002
@test Dates.yearwrap(2000,12,-3600) == 1700
@test Dates.yearwrap(2000,12,-36) == 1997
@test Dates.yearwrap(2000,12,-35) == 1998
@test Dates.yearwrap(2000,12,-24) == 1998
@test Dates.yearwrap(2000,12,-23) == 1999
@test Dates.yearwrap(2000,12,-14) == 1999
@test Dates.yearwrap(2000,12,-13) == 1999
@test Dates.yearwrap(2000,12,-12) == 1999
@test Dates.yearwrap(2000,12,-11) == 2000
@test Dates.yearwrap(2000,12,-2) == 2000
@test Dates.yearwrap(2000,12,-1) == 2000
@test Dates.yearwrap(2000,12,0) == 2000
@test Dates.yearwrap(2000,12,1) == 2001
@test Dates.yearwrap(2000,12,11) == 2001
@test Dates.yearwrap(2000,12,12) == 2001
@test Dates.yearwrap(2000,12,13) == 2002
@test Dates.yearwrap(2000,12,24) == 2002
@test Dates.yearwrap(2000,12,25) == 2003
@test Dates.yearwrap(2000,12,36) == 2003
@test Dates.yearwrap(2000,12,37) == 2004
@test Dates.yearwrap(2000,12,3600) == 2300

dt = Dates.DateTime(1999,12,27)
@test dt + Dates.Month(1) == Dates.DateTime(2000,1,27)
@test dt + Dates.Month(-1) == Dates.DateTime(1999,11,27)
@test dt + Dates.Month(-11) == Dates.DateTime(1999,1,27)
@test dt + Dates.Month(11) == Dates.DateTime(2000,11,27)
@test dt + Dates.Month(-12) == Dates.DateTime(1998,12,27)
@test dt + Dates.Month(12) == Dates.DateTime(2000,12,27)
@test dt + Dates.Month(13) == Dates.DateTime(2001,1,27)
@test dt + Dates.Month(100) == Dates.DateTime(2008,4,27)
@test dt + Dates.Month(1000) == Dates.DateTime(2083,4,27)
@test dt - Dates.Month(1) == Dates.DateTime(1999,11,27)
@test dt - Dates.Month(-1) == Dates.DateTime(2000,1,27)
@test dt - Dates.Month(100) == Dates.DateTime(1991,8,27)
@test dt - Dates.Month(1000) == Dates.DateTime(1916,8,27)
dt = Dates.DateTime(2000,2,29)
@test dt + Dates.Month(1) == Dates.DateTime(2000,3,29)
@test dt - Dates.Month(1) == Dates.DateTime(2000,1,29)
dt = Dates.DateTime(1972,6,30,23,59,59)
@test dt + Dates.Month(1) == Dates.DateTime(1972,7,30,23,59,59)
@test dt - Dates.Month(1) == Dates.DateTime(1972,5,30,23,59,59)
@test dt + Dates.Month(-1) == Dates.DateTime(1972,5,30,23,59,59)

dt = Dates.DateTime(1999,12,27)
@test dt + Dates.Week(1) == Dates.DateTime(2000,1,3)
@test dt + Dates.Week(100) == Dates.DateTime(2001,11,26)
@test dt + Dates.Week(1000) == Dates.DateTime(2019,2,25)
@test dt - Dates.Week(1) == Dates.DateTime(1999,12,20)
@test dt - Dates.Week(100) == Dates.DateTime(1998,1,26)
@test dt - Dates.Week(1000) == Dates.DateTime(1980,10,27)
dt = Dates.DateTime(2000,2,29)
@test dt + Dates.Week(1) == Dates.DateTime(2000,3,7)
@test dt - Dates.Week(1) == Dates.DateTime(2000,2,22)
dt = Dates.DateTime(1972,6,30,23,59,59)
@test dt + Dates.Week(1) == Dates.DateTime(1972,7,7,23,59,59)
@test dt - Dates.Week(1) == Dates.DateTime(1972,6,23,23,59,59)
@test dt + Dates.Week(-1) == Dates.DateTime(1972,6,23,23,59,59)

dt = Dates.DateTime(1999,12,27)
@test dt + Dates.Day(1) == Dates.DateTime(1999,12,28)
@test dt + Dates.Day(100) == Dates.DateTime(2000,4,5)
@test dt + Dates.Day(1000) == Dates.DateTime(2002,9,22)
@test dt - Dates.Day(1) == Dates.DateTime(1999,12,26)
@test dt - Dates.Day(100) == Dates.DateTime(1999,9,18)
@test dt - Dates.Day(1000) == Dates.DateTime(1997,4,1)
dt = Dates.DateTime(1972,6,30,23,59,59)
@test dt + Dates.Day(1) == Dates.DateTime(1972,7,1,23,59,59)
@test dt - Dates.Day(1) == Dates.DateTime(1972,6,29,23,59,59)
@test dt + Dates.Day(-1) == Dates.DateTime(1972,6,29,23,59,59)

dt = Dates.DateTime(1999,12,27)
@test dt + Dates.Hour(1) == Dates.DateTime(1999,12,27,1)
@test dt + Dates.Hour(100) == Dates.DateTime(1999,12,31,4)
@test dt + Dates.Hour(1000) == Dates.DateTime(2000,2,6,16)
@test dt - Dates.Hour(1) == Dates.DateTime(1999,12,26,23)
@test dt - Dates.Hour(100) == Dates.DateTime(1999,12,22,20)
@test dt - Dates.Hour(1000) == Dates.DateTime(1999,11,15,8)
dt = Dates.DateTime(1972,6,30,23,59,59)
@test dt + Dates.Hour(1) == Dates.DateTime(1972,7,1,0,59,59)
@test dt - Dates.Hour(1) == Dates.DateTime(1972,6,30,22,59,59)
@test dt + Dates.Hour(-1) == Dates.DateTime(1972,6,30,22,59,59)

dt = Dates.DateTime(1999,12,27)
@test dt + Dates.Minute(1) == Dates.DateTime(1999,12,27,0,1)
@test dt + Dates.Minute(100) == Dates.DateTime(1999,12,27,1,40)
@test dt + Dates.Minute(1000) == Dates.DateTime(1999,12,27,16,40)
@test dt - Dates.Minute(1) == Dates.DateTime(1999,12,26,23,59)
@test dt - Dates.Minute(100) == Dates.DateTime(1999,12,26,22,20)
@test dt - Dates.Minute(1000) == Dates.DateTime(1999,12,26,7,20)
dt = Dates.DateTime(1972,6,30,23,59,59)
@test dt + Dates.Minute(1) == Dates.DateTime(1972,7,1,0,0,59)
@test dt - Dates.Minute(1) == Dates.DateTime(1972,6,30,23,58,59)
@test dt + Dates.Minute(-1) == Dates.DateTime(1972,6,30,23,58,59)

dt = Dates.DateTime(1999,12,27)
@test dt + Dates.Second(1) == Dates.DateTime(1999,12,27,0,0,1)
@test dt + Dates.Second(100) == Dates.DateTime(1999,12,27,0,1,40)
@test dt + Dates.Second(1000) == Dates.DateTime(1999,12,27,0,16,40)
@test dt - Dates.Second(1) == Dates.DateTime(1999,12,26,23,59,59)
@test dt - Dates.Second(100) == Dates.DateTime(1999,12,26,23,58,20)
@test dt - Dates.Second(1000) == Dates.DateTime(1999,12,26,23,43,20)

dt = Dates.DateTime(1999,12,27)
@test dt + Dates.Millisecond(1) == Dates.DateTime(1999,12,27,0,0,0,1)
@test dt + Dates.Millisecond(100) == Dates.DateTime(1999,12,27,0,0,0,100)
@test dt + Dates.Millisecond(1000) == Dates.DateTime(1999,12,27,0,0,1)
@test dt - Dates.Millisecond(1) == Dates.DateTime(1999,12,26,23,59,59,999)
@test dt - Dates.Millisecond(100) == Dates.DateTime(1999,12,26,23,59,59,900)
@test dt - Dates.Millisecond(1000) == Dates.DateTime(1999,12,26,23,59,59)
dt = Dates.DateTime(1972,6,30,23,59,59)
@test dt + Dates.Millisecond(1) == Dates.DateTime(1972,6,30,23,59,59,1)
@test dt - Dates.Millisecond(1) == Dates.DateTime(1972,6,30,23,59,58,999)
@test dt + Dates.Millisecond(-1) == Dates.DateTime(1972,6,30,23,59,58,999)

dt = Dates.Date(1999,12,27)
@test dt + Dates.Year(1) == Dates.Date(2000,12,27)
@test dt + Dates.Year(100) == Dates.Date(2099,12,27)
@test dt + Dates.Year(1000) == Dates.Date(2999,12,27)
@test dt - Dates.Year(1) == Dates.Date(1998,12,27)
@test dt - Dates.Year(100) == Dates.Date(1899,12,27)
@test dt - Dates.Year(1000) == Dates.Date(999,12,27)
dt = Dates.Date(2000,2,29)
@test dt + Dates.Year(1) == Dates.Date(2001,2,28)
@test dt - Dates.Year(1) == Dates.Date(1999,2,28)
@test dt + Dates.Year(4) == Dates.Date(2004,2,29)
@test dt - Dates.Year(4) == Dates.Date(1996,2,29)

dt = Dates.Date(1999,12,27)
@test dt + Dates.Month(1) == Dates.Date(2000,1,27)
@test dt + Dates.Month(100) == Dates.Date(2008,4,27)
@test dt + Dates.Month(1000) == Dates.Date(2083,4,27)
@test dt - Dates.Month(1) == Dates.Date(1999,11,27)
@test dt - Dates.Month(100) == Dates.Date(1991,8,27)
@test dt - Dates.Month(1000) == Dates.Date(1916,8,27)
dt = Dates.Date(2000,2,29)
@test dt + Dates.Month(1) == Dates.Date(2000,3,29)
@test dt - Dates.Month(1) == Dates.Date(2000,1,29)

dt = Dates.Date(1999,12,27)
@test dt + Dates.Week(1) == Dates.Date(2000,1,3)
@test dt + Dates.Week(100) == Dates.Date(2001,11,26)
@test dt + Dates.Week(1000) == Dates.Date(2019,2,25)
@test dt - Dates.Week(1) == Dates.Date(1999,12,20)
@test dt - Dates.Week(100) == Dates.Date(1998,1,26)
@test dt - Dates.Week(1000) == Dates.Date(1980,10,27)
dt = Dates.Date(2000,2,29)
@test dt + Dates.Week(1) == Dates.Date(2000,3,7)
@test dt - Dates.Week(1) == Dates.Date(2000,2,22)

dt = Dates.Date(1999,12,27)
@test dt + Dates.Day(1) == Dates.Date(1999,12,28)
@test dt + Dates.Day(100) == Dates.Date(2000,4,5)
@test dt + Dates.Day(1000) == Dates.Date(2002,9,22)
@test dt - Dates.Day(1) == Dates.Date(1999,12,26)
@test dt - Dates.Day(100) == Dates.Date(1999,9,18)
@test dt - Dates.Day(1000) == Dates.Date(1997,4,1)

# Vectorized arithmetic
a = Dates.Date(2014,1,1)
dr = [a,a,a,a,a,a,a,a,a,a]
b = a + Dates.Year(1)
@test dr .+ Dates.Year(1) == repmat([b],10)
b = a + Dates.Month(1)
@test dr .+ Dates.Month(1) == repmat([b],10)
b = a + Dates.Day(1)
@test dr .+ Dates.Day(1) == repmat([b],10)
b = a - Dates.Year(1)
@test dr .- Dates.Year(1) == repmat([b],10)
b = a - Dates.Month(1)
@test dr .- Dates.Month(1) == repmat([b],10)
b = a - Dates.Day(1)
@test dr .- Dates.Day(1) == repmat([b],10)

# Vectorized arithmetic
b = a + Dates.Year(1)
@test dr .+ Dates.Year(1) == repmat([b],10)
b = a + Dates.Month(1)
@test dr .+ Dates.Month(1) == repmat([b],10)
b = a + Dates.Day(1)
@test dr .+ Dates.Day(1) == repmat([b],10)
b = a - Dates.Year(1)
@test dr .- Dates.Year(1) == repmat([b],10)
b = a - Dates.Month(1)
@test dr .- Dates.Month(1) == repmat([b],10)
b = a - Dates.Day(1)
@test dr .- Dates.Day(1) == repmat([b],10)

# Month arithmetic minimizes "edit distance", or number of changes
# needed to get a correct answer
# This approach results in a few cases of non-associativity
a = Dates.Date(2012,1,29)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = Dates.Date(2012,1,30)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = Dates.Date(2012,2,29)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = Dates.Date(2012,3,30)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = Dates.Date(2012,4,30)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = Dates.Date(2012,5,30)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = Dates.Date(2012,6,30)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = Dates.Date(2012,8,30)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = Dates.Date(2012,9,30)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = Dates.Date(2012,10,30)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)
a = Dates.Date(2012,11,30)
@test (a+Dates.Day(1))+Dates.Month(1) != (a+Dates.Month(1))+Dates.Day(1)

dt = Dates.DateTime(2000,1,1,12,30,45,500)
dt2 = dt + Dates.Year(1)
@test Dates.year(dt2) == 2001
@test Dates.month(dt2) == 1
@test Dates.day(dt2) == 1
@test Dates.hour(dt2) == 12
@test Dates.minute(dt2) == 30
@test Dates.second(dt2) == 45
@test Dates.millisecond(dt2) == 500

t1 = [Dates.Date(2009,1,1) Dates.Date(2009,1,2) Dates.Date(2009,1,3); Dates.Date(2009,2,1) Dates.Date(2009,2,2) Dates.Date(2009,2,3)]
t2 = [Dates.Date(2009,1,2) Dates.Date(2009,2,2) Dates.Date(2010,1,3); Dates.Date(2010,2,1) Dates.Date(2009,3,2) Dates.Date(2009,2,4)]
t3 = [Dates.DateTime(2009,1,1), Dates.DateTime(2009,1,2), Dates.DateTime(2009,1,3)]
t4 = [Dates.DateTime(2009,1,1,0,0,1), Dates.DateTime(2009,1,2,0,1), Dates.DateTime(2009,1,3,1)]

# TimeType, Array{TimeType}
@test Dates.Date(2010,1,1) .- t1 == [Dates.Day(365) Dates.Day(364) Dates.Day(363); Dates.Day(334) Dates.Day(333) Dates.Day(332)]
@test t1 .- Dates.Date(2010,1,1) == [Dates.Day(-365) Dates.Day(-364) Dates.Day(-363); Dates.Day(-334) Dates.Day(-333) Dates.Day(-332)]
@test Dates.DateTime(2009,1,1) .- t3 == [Dates.Millisecond(0), Dates.Millisecond(-86400000), Dates.Millisecond(-172800000)]
@test t3 .- Dates.DateTime(2009,1,1) == [Dates.Millisecond(0), Dates.Millisecond(86400000), Dates.Millisecond(172800000)]
@test Dates.Date(2010,1,1) - t1 == [Dates.Day(365) Dates.Day(364) Dates.Day(363); Dates.Day(334) Dates.Day(333) Dates.Day(332)]
@test t1 - Dates.Date(2010,1,1) == [Dates.Day(-365) Dates.Day(-364) Dates.Day(-363); Dates.Day(-334) Dates.Day(-333) Dates.Day(-332)]
@test Dates.DateTime(2009,1,1) - t3 == [Dates.Millisecond(0), Dates.Millisecond(-86400000), Dates.Millisecond(-172800000)]
@test t3 - Dates.DateTime(2009,1,1) == [Dates.Millisecond(0), Dates.Millisecond(86400000), Dates.Millisecond(172800000)]

# GeneralPeriod, Array{TimeType}
@test Dates.Day(1) .+ t1 == [Dates.Date(2009,1,2) Dates.Date(2009,1,3) Dates.Date(2009,1,4); Dates.Date(2009,2,2) Dates.Date(2009,2,3) Dates.Date(2009,2,4)]
@test Dates.Hour(1) .+ t3 == [Dates.DateTime(2009,1,1,1), Dates.DateTime(2009,1,2,1), Dates.DateTime(2009,1,3,1)]
@test t1 .+ Dates.Day(1) == [Dates.Date(2009,1,2) Dates.Date(2009,1,3) Dates.Date(2009,1,4); Dates.Date(2009,2,2) Dates.Date(2009,2,3) Dates.Date(2009,2,4)]
@test t3 .+ Dates.Hour(1) == [Dates.DateTime(2009,1,1,1), Dates.DateTime(2009,1,2,1), Dates.DateTime(2009,1,3,1)]
@test Dates.Day(1) + t1 == [Dates.Date(2009,1,2) Dates.Date(2009,1,3) Dates.Date(2009,1,4); Dates.Date(2009,2,2) Dates.Date(2009,2,3) Dates.Date(2009,2,4)]
@test Dates.Hour(1) + t3 == [Dates.DateTime(2009,1,1,1), Dates.DateTime(2009,1,2,1), Dates.DateTime(2009,1,3,1)]
@test t1 + Dates.Day(1) == [Dates.Date(2009,1,2) Dates.Date(2009,1,3) Dates.Date(2009,1,4); Dates.Date(2009,2,2) Dates.Date(2009,2,3) Dates.Date(2009,2,4)]
@test t3 + Dates.Hour(1) == [Dates.DateTime(2009,1,1,1), Dates.DateTime(2009,1,2,1), Dates.DateTime(2009,1,3,1)]

@test (Dates.Month(1) + Dates.Day(1)) .+ t1 == [Dates.Date(2009,2,2) Dates.Date(2009,2,3) Dates.Date(2009,2,4); Dates.Date(2009,3,2) Dates.Date(2009,3,3) Dates.Date(2009,3,4)]
@test (Dates.Hour(1) + Dates.Minute(1)) .+ t3 == [Dates.DateTime(2009,1,1,1,1), Dates.DateTime(2009,1,2,1,1), Dates.DateTime(2009,1,3,1,1)]
@test t1 .+ (Dates.Month(1) + Dates.Day(1)) == [Dates.Date(2009,2,2) Dates.Date(2009,2,3) Dates.Date(2009,2,4); Dates.Date(2009,3,2) Dates.Date(2009,3,3) Dates.Date(2009,3,4)]
@test t3 .+ (Dates.Hour(1) + Dates.Minute(1)) == [Dates.DateTime(2009,1,1,1,1), Dates.DateTime(2009,1,2,1,1), Dates.DateTime(2009,1,3,1,1)]
@test (Dates.Month(1) + Dates.Day(1)) + t1 == [Dates.Date(2009,2,2) Dates.Date(2009,2,3) Dates.Date(2009,2,4); Dates.Date(2009,3,2) Dates.Date(2009,3,3) Dates.Date(2009,3,4)]
@test (Dates.Hour(1) + Dates.Minute(1)) + t3 == [Dates.DateTime(2009,1,1,1,1), Dates.DateTime(2009,1,2,1,1), Dates.DateTime(2009,1,3,1,1)]
@test t1 + (Dates.Month(1) + Dates.Day(1)) == [Dates.Date(2009,2,2) Dates.Date(2009,2,3) Dates.Date(2009,2,4); Dates.Date(2009,3,2) Dates.Date(2009,3,3) Dates.Date(2009,3,4)]
@test t3 + (Dates.Hour(1) + Dates.Minute(1)) == [Dates.DateTime(2009,1,1,1,1), Dates.DateTime(2009,1,2,1,1), Dates.DateTime(2009,1,3,1,1)]

@test Dates.Day(1) .- t1 == [Dates.Date(2008,12,31) Dates.Date(2009,1,1) Dates.Date(2009,1,2); Dates.Date(2009,1,31) Dates.Date(2009,2,1) Dates.Date(2009,2,2)]
@test Dates.Hour(1) .- t3 == [Dates.DateTime(2008,12,31,23), Dates.DateTime(2009,1,1,23), Dates.DateTime(2009,1,2,23)]
@test t1 .- Dates.Day(1) == [Dates.Date(2008,12,31) Dates.Date(2009,1,1) Dates.Date(2009,1,2); Dates.Date(2009,1,31) Dates.Date(2009,2,1) Dates.Date(2009,2,2)]
@test t3 .- Dates.Hour(1) == [Dates.DateTime(2008,12,31,23), Dates.DateTime(2009,1,1,23), Dates.DateTime(2009,1,2,23)]
@test Dates.Day(1) - t1 == [Dates.Date(2008,12,31) Dates.Date(2009,1,1) Dates.Date(2009,1,2); Dates.Date(2009,1,31) Dates.Date(2009,2,1) Dates.Date(2009,2,2)]
@test Dates.Hour(1) - t3 == [Dates.DateTime(2008,12,31,23), Dates.DateTime(2009,1,1,23), Dates.DateTime(2009,1,2,23)]
@test t1 - Dates.Day(1) == [Dates.Date(2008,12,31) Dates.Date(2009,1,1) Dates.Date(2009,1,2); Dates.Date(2009,1,31) Dates.Date(2009,2,1) Dates.Date(2009,2,2)]
@test t3 - Dates.Hour(1) == [Dates.DateTime(2008,12,31,23), Dates.DateTime(2009,1,1,23), Dates.DateTime(2009,1,2,23)]

@test (Dates.Month(1) + Dates.Day(1)) .- t1 == [Dates.Date(2008,11,30) Dates.Date(2008,12,1) Dates.Date(2008,12,2); Dates.Date(2008,12,31) Dates.Date(2009,1,1) Dates.Date(2009,1,2)]
@test (Dates.Hour(1) + Dates.Minute(1)) .- t3 == [Dates.DateTime(2008,12,31,22,59), Dates.DateTime(2009,1,1,22,59), Dates.DateTime(2009,1,2,22,59)]
@test t1 .- (Dates.Month(1) + Dates.Day(1)) == [Dates.Date(2008,11,30) Dates.Date(2008,12,1) Dates.Date(2008,12,2); Dates.Date(2008,12,31) Dates.Date(2009,1,1) Dates.Date(2009,1,2)]
@test t3 .- (Dates.Hour(1) + Dates.Minute(1)) == [Dates.DateTime(2008,12,31,22,59), Dates.DateTime(2009,1,1,22,59), Dates.DateTime(2009,1,2,22,59)]
@test (Dates.Month(1) + Dates.Day(1)) - t1 == [Dates.Date(2008,11,30) Dates.Date(2008,12,1) Dates.Date(2008,12,2); Dates.Date(2008,12,31) Dates.Date(2009,1,1) Dates.Date(2009,1,2)]
@test (Dates.Hour(1) + Dates.Minute(1)) - t3 == [Dates.DateTime(2008,12,31,22,59), Dates.DateTime(2009,1,1,22,59), Dates.DateTime(2009,1,2,22,59)]
@test t1 - (Dates.Month(1) + Dates.Day(1)) == [Dates.Date(2008,11,30) Dates.Date(2008,12,1) Dates.Date(2008,12,2); Dates.Date(2008,12,31) Dates.Date(2009,1,1) Dates.Date(2009,1,2)]
@test t3 - (Dates.Hour(1) + Dates.Minute(1)) == [Dates.DateTime(2008,12,31,22,59), Dates.DateTime(2009,1,1,22,59), Dates.DateTime(2009,1,2,22,59)]

# Array{TimeType}, Array{TimeType}
@test t2 - t1 == [Dates.Day(1) Dates.Day(31) Dates.Day(365); Dates.Day(365) Dates.Day(28) Dates.Day(1)]
@test t4 - t3 == [Dates.Millisecond(1000), Dates.Millisecond(60000), Dates.Millisecond(3600000)]
@test (Dates.Date(2009,1,1):Dates.Week(1):Dates.Date(2009,1,21)) - (Dates.Date(2009,1,1):Dates.Day(1):Dates.Date(2009,1,3)) == [Dates.Day(0), Dates.Day(6), Dates.Day(12)]
@test (Dates.DateTime(2009,1,1,1,1,1):Dates.Second(1):Dates.DateTime(2009,1,1,1,1,3)) - (Dates.DateTime(2009,1,1,1,1):Dates.Second(1):Dates.DateTime(2009,1,1,1,1,2)) == [Dates.Second(1), Dates.Second(1), Dates.Second(1)]

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
#
# Julia System Desciption Language
#
module Causal

using UUIDs
using DifferentialEquations
using Sundials
using LightGraphs
using DataStructures
using JLD2
using Plots
using ProgressMeter
using Logging
using LinearAlgebra
using Dates
using NLsolve
using Interpolations
using LibGit2
using DocStringExtensions
import GraphPlot.gplot
import FileIO: load
import Base: show, display, write, read, close, setproperty!, mv, cp, open,  istaskdone, istaskfailed,
    getindex, setindex!, size, isempty

include("utilities/utils.jl")
export equip

include("utilities/callback.jl")
export Callback,
    enable!,
    disable!,
    isenabled,
    applycallbacks

include("utilities/buffer.jl")
export BufferMode,
    LinearMode,
    CyclicMode,
    Buffer,
    Normal,
    Cyclic,
    Fifo,
    Lifo,
    write!,
    isfull,
    ishit,
    content,
    mode,
    snapshot,
    datalength,
    inbuf,
    outbuf

include("connections/link.jl")
export Link,
    launch,
    refresh!

include("connections/pin.jl")
export AbstractPin,
    Outpin,
    Inpin,
    connect!,
    disconnect!,
    isconnected,
    isbound

include("connections/port.jl")
export AbstractPort,
    Inport,
    Outport,
    datatype

include("components/componentsbase/hierarchy.jl")
export AbstractComponent,
    AbstractSource,
    AbstractSystem,
    AbstractSink,
    AbstractStaticSystem,
    AbstractDynamicSystem,
    AbstractSubSystem,
    AbstractMemory,
    AbstractDiscreteSystem,
    AbstractODESystem,
    AbstractRODESystem,
    AbstractDAESystem,
    AbstractSDESystem,
    AbstractDDESystem

include("components/componentsbase/macros.jl")
include("components/componentsbase/interpolant.jl")
export Interpolant

include("components/componentsbase/takestep.jl")
export readtime!,
    readstate,
    readinput!,
    writeoutput!,
    computeoutput,
    evolve!,
    takestep!,
    drive!,
    approve!

include("components/constants.jl")

include("components/sources/clock.jl")
export Clock,
    isrunning,
    ispaused,
    isoutoftime,
    set!,
    stop!,
    pause!

include("components/sources/generators.jl")
export @def_source,
    FunctionGenerator,
    SinewaveGenerator,
    DampedSinewaveGenerator,
    SquarewaveGenerator,
    TriangularwaveGenerator,
    ConstantGenerator,
    RampGenerator,
    StepGenerator,
    ExponentialGenerator,
    DampedExponentialGenerator

include("components/systems/staticsystems/staticsystems.jl")
export @def_static_system,
    StaticSystem,
    Adder,
    Multiplier,
    Gain,
    Terminator,
    Memory,
    Coupler,
    Differentiator

include("components/systems/dynamicalsystems/init.jl")

include("components/systems/dynamicalsystems/odesystems.jl")
export @def_ode_system,
    ODESystem,
    ContinuousLinearSystem,
    LorenzSystem,
    ChenSystem,
    ChuaSystem,
    RosslerSystem,
    VanderpolSystem,
    ForcedLorenzSystem,
    ForcedChenSystem,
    ForcedChuaSystem,
    ForcedRosslerSystem,
    ForcedVanderpolSystem,
    Integrator

include("components/systems/dynamicalsystems/discretesystems.jl")
export @def_discrete_system,
    DiscreteSystem,
    DiscreteLinearSystem,
    HenonSystem,
    LoziSystem,
    BogdanovSystem,
    GingerbreadmanSystem,
    LogisticSystem

include("components/systems/dynamicalsystems/sdesystems.jl")
export @def_sde_system,
    SDESystem,
    NoisyLorenzSystem,
    ForcedNoisyLorenzSystem

include("components/systems/dynamicalsystems/daesystems.jl")
export @def_dae_system,
    DAESystem,
    RobertsonSystem,
    PendulumSystem,
    RLCSystem

include("components/systems/dynamicalsystems/rodesystems.jl")
export @def_rode_system,
    RODESystem,
    MultiplicativeNoiseLinearSystem

include("components/systems/dynamicalsystems/ddesystems.jl")
export @def_dde_system,
    DDESystem,
    DelayFeedbackSystem

include("components/sinks/sinks.jl")
export @def_sink,
    Writer,
    Printer,
    Scope,
    write!,
    fwrite!,
    fread,
    update!

include("models/taskmanager.jl")
export TaskManager,
    checktaskmanager

include("models/simulation.jl")
export Simulation,
    setlogger,
    closelogger,
    report

include("models/model.jl")
export @defmodel,
    Model,
    inspect!,
    initialize!,
    run!,
    terminate!,
    simulate!,
    getloops,
    breakloop!,
    Node,
    Branch,
    addnode!,
    getnode,
    addbranch!,
    getbranch,
    deletebranch!,
    signalflow,
    troubleshoot

include("plugins/loadplugins.jl")
export AbstractPlugin,
    process,
    add,
    remove,
    enable,
    disable,
    check,
    @def_plugin

end  # module


text:
using SuperResolution
using Test
using DataFrames

@testset "importSR.jl" begin
    @test typeof(SuperResolution.DataSim)==DataFrame
    @test typeof(SuperResolution.RealLocs)==DataFrame
    # Load Data with colnames:
    # mol_ID|loc_ID|framenumber|pos_x|pos_y
    # [Int64,Int64,Int64,Float64,Float64,Float64,Float64,Float64,Float64,Float64]
    # Write your tests here.
end

max_stars_repo_name: PGSch/SuperResolution

max_issues_repo_name: PGSch/SuperResolution

max_forks_repo_name: PGSch/SuperResolution


text:
@doc """
A SparseRollCall object is a variant of the ZeroOneRollCall encoding that
does not explicitly represent any of the missing or inactive votes. This
can be more efficient for model fitting since the missing and inactive notes
are ignored by many models.
""" ->
immutable SparseRollCall
    n_legislators::Int
    n_bills::Int
    legislators::Vector{Int}
    bills::Vector{Int}
    votes::Vector{Float64}
end

@doc """
Convert a RollCall object into a SparseRollCall object.
""" ->
function Base.convert(::Type{SparseRollCall}, roll_call::RollCall)
    encoding, old_votes = roll_call.encoding, roll_call.votes
    n_legislators, n_bills = size(old_votes)
    count_nans = 0
    for vote in old_votes
        count_nans += ifelse(isnan(transcode(vote, encoding)), 1, 0)
    end
    count_not_nans = length(old_votes) - count_nans
    legislators = Array(Int, count_not_nans)
    bills = Array(Int, count_not_nans)
    votes = Array(Float64, count_not_nans)
    i = 0
    for l in 1:n_legislators
        for b in 1:n_bills
            vote = transcode(old_votes[l, b], encoding)
            if !isnan(vote)
                i += 1
                legislators[i] = l
                bills[i] = b
                votes[i] = vote
            end
        end
    end
    SparseRollCall(n_legislators, n_bills, legislators, bills, votes)
end

max_stars_repo_name: johnmyleswhite/RollCallDataIO

max_issues_repo_name: johnmyleswhite/RollCallDataIO

max_forks_repo_name: johnmyleswhite/RollCallDataIO


text:
# Example from Improved relaxations for the parametric solutionsof ODEs using differential inequalities
id = :ScottA2013a

pL = [-6.5; 0.01]
pU = [1.0; 0.05]
ScottA2013a_x0(p::Vector{T}) where T = T[one(T); 0.5*one(T)]
function ScottA2013a_f!(du::Vector{T}, u::Vector{T}, p::Vector{T}, t) where T
    du[1] = -(2.0 + sin(p[1]/3.0))*u[1]^2 + p[2]*u[1]*u[2]
    du[2] = sin(p[1]/3.0)*u[1]^2 - p[2]*u[1]*u[2]
    return
end
tspan = (0.0, 2.0)

prob = ODERelaxProb(ScottA2013a_f!, tspan, ScottA2013a_x0, pL, pU)
url = "https://link.springer.com/article/10.1007%2Fs10898-012-9909-0"
source = "ScottA2013"
desc = "Sinusodial Example 1"


text:
using Statistics

@testset "Testing core Scheduler creation" begin
    es = EmptyState()
    f = (a,b)->nothing
    sReal = Scheduler{EmptyState, Real}(1.2, Vector{EventSimulation.Action{Real}}(),
                                        es, f)
    @test isa(sReal, Scheduler{EmptyState, Real})
    @test sReal.now == 1.2
    @test isempty(sReal.event_queue)
    @test isa(sReal.state, EmptyState)
    @test isa(sReal.monitor(0, sReal), Nothing)
end

@testset "Testing if register! and go! produce correct order" begin
    s = Scheduler()

    mlen = 100
    ts = (1.0:mlen)[randperm(mlen)]

    for x in ts
        register!(s, z -> push!(v, z.now), x)
    end

    v = Float64[]
    go!(s, mlen)
    @test v == [1.0:mlen;]
    @test length(s.event_queue) == 0
end

@testset "Testing register! and terminate! for correct number of elements" begin
    s = Scheduler()

    mlen = 100
    ts = (1.0:mlen)[randperm(mlen)]
    for x in ts
        register!(s, z -> push!(v, z.now), x)
    end

    @test length(s.event_queue) == mlen
    terminate!(s)
    @test length(s.event_queue) == 0
end

@testset "Testing repeat_register!" begin
    s = Scheduler()
    rrv = []
    f1(x) = push!(rrv, x.now)
    ifun(x) = x.now + 1
    repeat_register!(s, f1, ifun)
    go!(s, 50)
    @test rrv == [1,3,7,15,31]

    s = Scheduler()
    rrv = []
    f2(x) = push!(rrv, x.now)
    ifun2(x) = length(rrv) < 5 ? x.now + 1 : NaN
    repeat_register!(s, f2, ifun2)
    go!(s)
    @test length(rrv) == 5
end

@testset "Testing ordered bulk_register!" begin
    s = Scheduler()
    v2 = []
    who = [1:5;]
    for x in randperm(10)
        bulk_register!(s, who, (a,b) -> push!(v2, (s.now, b)), Float64(x), false)
    end
    go!(s, 11)
    v2t = []
    for i in 1.0:10.0, j in 1:5
        push!(v2t, (i, j))
    end

    @test v2 == v2t
end

@testset "Testing random bulk_register!" begin
    tmp_store = []
    s = Scheduler()
    v2 = []
    who = [1:100;]
    for x in randperm(10)
        bulk_register!(s, who, (a,b) -> push!(tmp_store, (s.now, b)), Float64(x), true)
    end
    go!(s, 11)
    v = getindex.(tmp_store, 1)
    @test v == sort(v)
    l = unique(v)
    for k in l
        r = [tmp_store[i][2] for i in 1:length(tmp_store) if tmp_store[i][1] == k]
        @test r != sort(r)
    end
end

@testset "Testing ordered repeat_bulk_register!" begin
    s = Scheduler()
    v2 = []
    who = [1:5;]
    repeat_bulk_register!(s, who, (a,b) -> push!(v2, (s.now, b)), x -> 1.0, false)
    go!(s, 10)
    v2t = []
    for i in 1.0:10.0, j in 1:5
        push!(v2t, (i, j))
    end

    @test v2 == v2t

    s = Scheduler()
    v2 = []
    who = [1:5;]
    repeat_bulk_register!(s, who, (a,b) -> push!(v2, (s.now, b)),
                          x -> length(v2) < length(v2t) ? 1.0 : nothing, false)
    go!(s)

    @test v2 == v2t
end

@testset "Testing random repeat_bulk_register!. Visual inspection required." begin
    tmp_store = []
    s = Scheduler()
    v2 = []
    who = [1:100;]
    repeat_bulk_register!(s, who, (a,b) -> push!(tmp_store, (s.now, b)), x -> 1.0, true)
    go!(s, 10)

    v = getindex.(tmp_store, 1)
    @test v == sort(v)
    l = unique(v)
    for k in l
        r = [tmp_store[i][2] for i in 1:length(tmp_store) if tmp_store[i][1] == k]
        @test r != sort(r)
    end
end

@testset "Testing interrupt." begin
    s = Scheduler()
    @test interrupt!(s, EventSimulation.Action(identity, 1.0)) == false
    for i in 1:10
        register!(s, identity, rand())
    end

    OK = []
    a = register!(s, x -> push!(OK, false), rand())

    for i in 1:10
        register!(s, identity, rand())
    end

    @test interrupt!(s, a) == true
    go!(s, 1.0)
    @test isempty(OK)
end

@testset "Testing monitor." begin
    s = Scheduler()
    deltas = Float64[]
    s.monitor = (s,Δ) -> push!(deltas, Δ)
    repeat_register!(s, x -> nothing, x -> rand())
    go!(s, 100_000)
    m = mean(deltas)
    v = var(deltas)
    @test max(abs(m-0.5), abs(v-1/12)) < 0.005
end


text:
################################################################
# used inside benchmarks.jl

include("distributions.jl")

################################################################
# used inside benchmarks.jl

"""
    numArgsDims(in)
Finding number of arguments and number of dimension sets

# Examples
```julia
numArgs, numDimsSets = numArgsDims(in)
```
"""
function numArgsDims(in)
    dimsDims = ndims(in)
    dimsSize  = size(in)
    if dimsDims == 1
        numArgs = dimsSize[1]
        numDimsSets = 1
    elseif dimsDims == 2
        numArgs = dimsSize[1]
        numDimsSets = dimsSize[2]
    end
    return numArgs, numDimsSets
end

################################################################
# used inside plots.jl

flatten(A) =reduce(hcat, A)
uniqueflatten(A) =  unique(flatten(A))

function stringMatrix(A)
    sprint(Base.print_matrix, A)
end


text:
using Turing

const n_samples = 10^4

#Sample from transformed distributions, check moments.
within_4_stdev(dist::Distribution, samples) = begin
    s = sqrt(var(dist))
    m = mean(dist)
    n = length(samples)
    stdevs = 4
    @assert 0.95*s < sqrt(var(samples)) < 1.05*s
    @assert m-stdevs*s/sqrt(n) < mean(samples) < m+stdevs*s/sqrt(n)
end

macro test_dist(Dist)
    dist = symbol(lowercase(string(Dist.args[1])))
    args = Dist.args[2:end]
    quote
        within_4_stdev($Dist, $dist(State(), $(args...); dims=($n_samples,)))
    end
    #dist = deepcopy(Dist)
    #dist.args[1] = symbol(lowercase(string(dist.args[1])))
    #insert!(dist.args, 2, :(State()))
    #push!(dist.args, :(($n_samples,)))
    #:(within_4_stdev($Dist, $dist))
end

@test_dist Normal(2, 3)
@test_dist LogNormal(1., 0.5)
#@test_dist Chisq(2)
#MvNormal
m = [1., 2.]
v = [5. -1; -1 3.]
samples = [mvnormal(State(), m, v) for i =1:n_samples]
println((m, mean(samples)))
println((v, cov(hcat(samples...)')))

@test_dist Exponential(3.)
@test_dist Bernoulli(0.7)
#@test_dist Uniform()
@test_dist Laplace(1., 3.)
@test_dist Categorical([0.1, 0.3, 0.6])
@test_dist Gamma(1., 3.)
@test_dist Beta(1., 3.)
#
#Dirichlet
println(); println(); println();
alpha = [1., 1.5, 0.5]
dist = Dirichlet(alpha)
samples = [dirichlet(State(), alpha) for i =1:n_samples]
println((mean(dist), mean(samples)))
println((cov(dist), cov(hcat(samples...)')))

#Multidimensional sampling fallbacks work
laplace(State(), ones(3), 1.)
laplace(State(), ones(3), ones(3))
gamma(State(), ones(3), 1.)
gamma(State(), ones(3), ones(3))


text:
# # add_format(format"Pulseq", "# Pulseq sequence file", ".seq", [:Koma=>UUID("a9056882-1c2f-47b4-848d-dcb4a04f1994")])
# module Pulseq
# This file is a copy of the MATLAB Pulseq v1.4.0 read functions

"""
read_version Read the [VERSION] section of a sequence file.
   defs=read_version(fid) Read Pulseq version from file
   identifier of an open MR sequence file and return it
"""
function read_version(io)
    major =    @scanf(readline(io), "major %i", Int)[end]
    minor =    @scanf(readline(io), "minor %i", Int)[end]
    revision = @scanf(readline(io), "revision %i", Int)[end]

    version_combined = 1_000_000*major+1_000*minor+revision

    @assert major == 1 "Unsupported version_major $major"
    if     version_combined < 1002000
        @error "Unsupported version $major.$minor.$revision, only file format revision 1.2.0 and above are supported"
    elseif version_combined < 1003001
        @warn "Loading older Pulseq format file (version $major.$minor.$revision) some code may not function as expected"
    end

    major, minor, revision, version_combined
end

"""
read_definitions Read the [DEFINITIONS] section of a sequence file.
   defs=read_definitions(fid) Read user definitions from file
   identifier of an open MR sequence file and return a map of
   key/value entries.
"""
function read_definitions(io)
    def = Dict{String,Any}()
    while true
        r,key,value_string1,value_string2,value_string3 = @scanf(readline(io), "%s %s %s %s", String, String, String, String)
        if r == 2
            value_numeric = tryparse.(Float64, value_string1)
            def[key] = (value_numeric === nothing) ? value_string1 : value_numeric
        end
        if r == 4
            value_string_array = [value_string1,value_string2,value_string3]
            value_numeric = tryparse.(Float64, value_string_array)
            def[key] = value_numeric[value_numeric .!== nothing]
        end
        (r == 2 || r == 4) || break #Break on white space
    end
    #Default values
    if !haskey(def,"BlockDurationRaster")       def["BlockDurationRaster"] = 1e-5       end
    if !haskey(def,"GradientRasterTime")        def["GradientRasterTime"] = 1e-5        end
    if !haskey(def,"RadiofrequencyRasterTime")  def["RadiofrequencyRasterTime"] = 1e-6  end
    if !haskey(def,"AdcRasterTime")             def["AdcRasterTime"] = 1e-7             end
    def
end

"""
read_blocks Read the [BLOCKS] section of a sequence file.
   library=read_blocks(fid) Read blocks from file identifier of an
   open MR sequence file and return the event table.
"""
function read_blocks(io, blockDurationRaster, version_combined)
    eventTable = Dict()
    blockDurations = Dict()
    delayIDs_tmp = Dict()
    while true
        if version_combined <= 1002001
            NumberBlockEvents = 7
        else
            NumberBlockEvents = 8
        end

        fmt = Scanf.Format("%i "^NumberBlockEvents)
        r, blockEvents... = scanf(readline(io), fmt, zeros(Int,NumberBlockEvents)...)

        if blockEvents[1] != 0
            if version_combined <= 1002001
                eventTable[blockEvents[1]] = [0 blockEvents[3:end]... 0]
            else
                eventTable[blockEvents[1]] = [0 blockEvents[3:end]...]
            end

            if version_combined >= 1004000
                blockDurations[blockEvents[1]] = blockEvents[2]*blockDurationRaster
            else
                delayIDs_tmp[blockEvents[1]] = blockEvents[2]
            end
        end

        r == NumberBlockEvents || break #Break on white space
    end
    sort(eventTable), sort(blockDurations), sort(delayIDs_tmp)
end

"""
read_events Read an event section of a sequence file.
   library=read_events(fid) Read event data from file identifier of
   an open MR sequence file and return a library of events.

   library=read_events(fid,scale) Read event data and scale
   elements according to column vector scale.

   library=read_events(fid,scale,type) Attach the type string to
   elements of the library.

   library=read_events(...,library) Append new events to the given
   library.
"""
function read_events(io, scale; type=-1, eventLibrary=Dict())
    while true
        EventLength = length(scale) + 1
        fmt = Scanf.Format("%i"*"%f "^(EventLength-1))
        r, data... = scanf(readline(io), fmt, Int, zeros(Float64,EventLength-1)...)
        id = floor(Int, data[1])
        data = scale .* [data[2:end]...]'
        if type != -1
            eventLibrary[id] = Dict("data"=>data, "type"=>type)
        else
            eventLibrary[id] = Dict("data"=>data)
        end
        r == EventLength || break #Break on white space
    end
    sort(eventLibrary)
end

"""
read_shapes Read the [SHAPES] section of a sequence file.
   library=read_shapes(fid) Read shapes from file identifier of an
   open MR sequence file and return a library of shapes.
"""
function read_shapes(io, forceConvertUncompressed)
    shapeLibrary = Dict()
    readline(io)
    while true #Reading shapes
        r, id = @scanf(readline(io), "shape_id %i", Int)
        r == 1 || break #Break if no "shape_id id" is identified
        _, num_samples = @scanf(readline(io), "num_samples %i", Int)
        shape = Float64[]
        while true #Reading shape data
            r, data_point = @scanf(readline(io), "%f", Float64)
            r == 1 || break #Break if no sample
            append!(shape, data_point)
        end
        # check if conversion is needed: in v1.4.x we use length(data)==num_samples
        # as a marker for the uncompressed (stored) data. In older versions this condition could occur by chance
        if forceConvertUncompressed && length(shape)==num_samples
            shape = compress_shape(decompress_shape(shape, num_samples; forceDecompression=true))
            data = (num_samples, shape)
        else
            data = (num_samples, shape)
        end
        shapeLibrary[id] = data
    end

    sort(shapeLibrary)
end

"""
compress_shape Compress a gradient or pulse shape.
   s=compress_shape(w) Compress the waveform using a run-length compression
   scheme on the derivative. This strategy encodes constant and linear
   waveforms with very few samples. Returns:
     num_samples - the number of samples in the uncompressed waveform
     data - containing the compressed waveform

   See also decompress_shape
"""
function compress_shape(w; forceCompression=false)
    num_samples = length(w)

    if !forceCompression && num_samples <= 4
        data = w[:]
    else
        quant_fac = 1e-7
        ws = w ./ quant_fac;
        datq = round.([ws[1]; diff(ws[:])])
        qerr = ws[:] .- cumsum(datq)
        qcor = [0; diff(round.(qerr))]
        datd = datq .+ qcor
        maskChanges = [true; diff(datd) .!= 0]
        vals = datd[maskChanges] .* quant_fac # Elements without repetitions

        k = findall([maskChanges; true]) # Indices of changes
        n = diff(k) # Number of repetitions

        # Encode in Pulseq format
        nExtra = convert(Vector{Any}, n .- 2)
        vals2 = convert(Vector{Any}, vals)
        vals2[nExtra .< 0] .= NaN;
        nExtra[nExtra .< 0] .= NaN;
        v = [vals'; vals2'; nExtra']
        v = convert(Vector{Float64}, v[isfinite.(v)])
        v[abs.(v) .<= 1e-10] .= 0

        # decide whether compression makes sense, otherwise store the original
        if forceCompression || num_samples > length(v)
            data = v
        else
            data = w
        end
    end

    num_samples, data
end

"""
decompress_shape Decompress a gradient or pulse shape.
   w=decompress_shape(shape) Decompress the shape compressed with a run-length
   compression scheme on the derivative. Returns:
     num_samples - the number of samples in the uncompressed waveform
     data - containing the compressed waveform

   See also compress_shape
"""
function decompress_shape(num_samples, data; forceDecompression = false)
    dataPack = data
    dataPackLen = length(dataPack)
    numSamples = num_samples

    if !forceDecompression && numSamples == dataPackLen
        w = dataPack
    else
        w = zeros(numSamples) #pre-allocate the results
        #Decompression starts
        dataPackDiff = dataPack[2:end] .- dataPack[1:end-1]
        #when dataPackDiff == 0 the subsequent samples are equal ==> marker for
        #repeats (run-length encoding)
        dataPackMarkers = findall(dataPackDiff .== 0)

        countPack = 1       # counter 1: points to the current compressed sample
        countUnpack = 1     # counter 2: points to the current uncompressed sample
        for i = 1:length(dataPackMarkers)
            nextPack = dataPackMarkers[i] # careful, this index may have "false positives" , e.g. if the value 3 repeats 3 times, then we will have 3 3 3
            currUnpackSamples = nextPack - countPack
            if currUnpackSamples < 0 # this rejects false positives
                continue
            elseif currUnpackSamples > 0 #do we have an unpacked block to copy?
                w[countUnpack:(countUnpack+currUnpackSamples-1)] .= dataPack[countPack:(nextPack-1)]
                countPack += currUnpackSamples
                countUnpack += currUnpackSamples
            end
            # now comes the packed/repeated section
            rep = floor(Int, dataPack[countPack+2] + 2)
            w[countUnpack:(countUnpack+rep-1)] .= dataPack[countPack]
            countPack += 3
            countUnpack += rep
        end
        # samples left?
        if countPack <= dataPackLen
            # println("$dataPackLen $countPack $numSamples $countUnpack")
            @assert dataPackLen-countPack == numSamples-countUnpack "Unsuccessful unpacking of samples"
            # copy the rest of the shape, it is unpacked
            w[countUnpack:end] .= dataPack[countPack:end]
        end

        w = cumsum(w)
    end

    w
end

"""
READ Load sequence from file.
   READ(seqObj, filename, ...) Read the given filename and load sequence
   data into sequence object.

   optional parwameter 'detectRFuse' can be given to let the function
   infer the currently missing flags concerning the intended use of the RF
   pulses (excitation, refocusing, etc). These are important for the
   k-space trajectory calculation

   Examples:
   Load the sequence defined in gre.seq in my_sequences directory

       read(seqObj,'my_sequences/gre.seq')

 See also  write
"""
function read_seq(filename)
    println("")
    @info "Loading sequence $(basename(filename)) ..."
    version_combined = 0
    version_major = 0
    version_minor = 0
    gradLibrary = Dict()
    def = Dict()
    blockEvents = Dict()
    blockDurations = Dict()
    delayInd_tmp = Dict()
    rfLibrary = Dict()
    adcLibrary = Dict()
    tmp_delayLibrary = Dict()
    shapeLibrary = Dict()
    extensionLibrary = Dict()
    #Reading file and storing data
    open(filename) do io
        while !eof(io)
            section = readline(io)
            if      section == "[DEFINITIONS]"
                def = read_definitions(io)
            elseif  section == "[VERSION]"
                version_major, version_minor, _, version_combined = read_version(io)
            elseif  section == "[BLOCKS]"
                if version_combined == 0
                    @error "Pulseq file MUST include [VERSION] section prior to [BLOCKS] section"
                end
                blockEvents, blockDurations, delayInd_tmp = read_blocks(io, def["BlockDurationRaster"], version_combined)
            elseif  section == "[RF]"
                if version_combined >= 1004000
                    rfLibrary = read_events(io, [1/γ 1 1 1 1e-6 1 1]) # this is 1.4.x format
                else
                    rfLibrary = read_events(io, [1/γ 1 1 1e-6 1 1]) # this is 1.3.x and below
                    # we will have to scan through the library later after all the shapes have been loaded
                end
            elseif  section == "[GRADIENTS]"
                if version_combined >= 1004000
                    gradLibrary = read_events(io, [1/γ 1 1 1e-6]; type='g', eventLibrary=gradLibrary) # this is 1.4.x format
                else
                    gradLibrary = read_events(io, [1/γ 1 1e-6];   type='g', eventLibrary=gradLibrary) # this is 1.3.x and below
                end
            elseif  section == "[TRAP]"
                gradLibrary = read_events(io, [1/γ 1e-6 1e-6 1e-6 1e-6]; type='t', eventLibrary=gradLibrary);
            elseif  section == "[ADC]"
                adcLibrary = read_events(io, [1 1e-9 1e-6 1 1])
            elseif  section == "[DELAYS]"
                if version_combined >= 1004000
                    @error "Pulseq file revision 1.4.0 and above MUST NOT contain [DELAYS] section"
                end
                tmp_delayLibrary = read_events(io, 1e-6);
            elseif  section == "[SHAPES]"
                shapeLibrary = read_shapes(io, (version_major==1 && version_minor<4))
            elseif  section == "[EXTENSIONS]"
                extensionLibrary = read_events(io,1) #For now, it reads the extensions but it does not take it them into account
            elseif  section == "[SIGNATURE]"
                #Not implemented yet
            end

        end
    end
    # fix blocks, gradients and RF objects imported from older versions
    if version_combined < 1004000
        # scan through the RF objects
        for i = 0:length(rfLibrary)-1
            rfLibrary[i]["data"] = [rfLibrary[i]["data"][1:3]' 0 rfLibrary[i]["data"][4:end]']
        end
        # scan through the gradient objects and update 't'-s (trapezoids) und 'g'-s (free-shape gradients)
        for i = 0:length(gradLibrary)-1
            if gradLibrary[i]["type"] == 't'
                #(1)amplitude (2)rise (2)flat (3)fall (4)delay
                if gradLibrary[i]["data"][2] == 0 #rise
                    if abs(gradLibrary[i]["data"][1]) == 0 && gradLibrary[i]["data"][3] > 0
                        gradLibrary[i]["data"][3] -= def["gradRasterTime"]
                        gradLibrary[i]["data"][2]  = def["gradRasterTime"]
                    end
                end
                if gradLibrary[i]["data"][4] == 0 #delay
                    if abs(gradLibrary[i]["data"][1]) == 0 && gradLibrary[i]["data"][3] > 0
                        gradLibrary[i]["data"][3] -= def["gradRasterTime"]
                        gradLibrary[i]["data"][4]  = def["gradRasterTime"]
                    end
                end
            end
            if gradLibrary[i]["type"] == 'g'
                #(1)amplitude (2)amp_shape_id (3)time_shape_id (4)delay
                gradLibrary[i]["data"] = [gradLibrary[i]["data"][1:2]; 0; gradLibrary[i]["data"][3:end]]
            end
        end
        # for versions prior to 1.4.0 blockDurations have not been initialized
        blockDurations = zeros(length(blockEvents))
        for i = 1:length(blockEvents)
        idelay = delayInd_tmp[i]
            if idelay > 0
                delay = tmp_delayLibrary[idelay]["data"][1]
                blockDurations[i] = delay
            end
        end
    end
    #Sequence
    obj = Dict(
        "blockEvents"=>blockEvents,
        "blockDurations"=>blockDurations,
        "delayInd_tmp"=>delayInd_tmp,
        "gradLibrary"=>gradLibrary,
        "rfLibrary"=>rfLibrary,
        "adcLibrary"=>adcLibrary,
        "tmp_delayLibrary"=>tmp_delayLibrary,
        "shapeLibrary"=>shapeLibrary,
        "extensionLibrary"=>extensionLibrary,
        "definitions"=>def)
    #Transforming Dictionary to Sequence object
    #This should only work for Pulseq files >=1.4.0
    seq = Sequence()
    for i = 1:length(blockEvents)
        seq += get_block(obj,i)
    end
    #Final details
    seq = seq[2:end] #Removed dummy seq block at the start
    seq.DEF = obj["definitions"]
    seq.DEF["FileName"] = basename(filename)
    seq.DEF["PulseqVersion"] = version_combined
    if !haskey(seq.DEF,"Nx")
        Nx = maximum(seq.ADC.N)
        RF_ex = (get_flip_angles(seq) .<= 90.01) .* is_RF_on.(seq)
        Nz = max(length(unique(seq.RF[RF_ex].Δf)), 1)
        Ny = sum(is_ADC_on.(seq)) / Nz |> x->floor(Int,x)

        seq.DEF["Nx"] = Nx  #Number of samples per ADC
        seq.DEF["Ny"] = Ny  #Number of ADC events
        seq.DEF["Nz"] = Nz  #Number of unique RF frequencies, in a 3D acquisition this should not work
        @warn "Pulseq file did not contain [DEFINITIONS] for Nx, Ny, and Nz. We added Nx=$Nx, Ny=$Ny, and Nz=$Nz to the Sequence object for the reconstruction."
    end
    #Koma sequence
    println("Successfully loaded $(basename(filename))!")
    seq
end

#To Sequence
function read_Grad(gradLibrary, shapeLibrary, Δt_gr, i)
    G = Grad(0,0)
    if gradLibrary[i]["type"] == 't' #if trapezoidal gradient
        #(1)amplitude (2)rise (3)flat (4)fall (5)delay
        g_A, g_rise, g_T, g_fall, g_delay = gradLibrary[i]["data"]
        G = Grad(g_A,g_T,g_rise,g_fall,g_delay)
    elseif gradLibrary[i]["type"] == 'g' #Arbitrary gradient waveform
        #(1)amplitude (2)amp_shape_id (3)time_shape_id (4)delay
        g = gradLibrary[i]["data"]
        amplitude =     g[1]
        amp_shape_id =  g[2] |> x->floor(Int64,x)
        time_shape_id = g[3] |> x->floor(Int64,x)
        delay =         g[4]
        #Amplitude
        gA = amplitude * decompress_shape(shapeLibrary[amp_shape_id]...)
        Nrf = length(gA) - 1
        #Creating timings
        if time_shape_id == 0 #no time waveform
            gT = Nrf * Δt_gr
        else
            gt = decompress_shape(shapeLibrary[time_shape_id]...)
            gT = (gt[2:end] .- gt[1:end-1]) * Δt_gr
        end
        G = Grad(gA,gT,(time_shape_id==0)*Δt_gr/2,0,delay)
    end
    G
end

function read_RF(rfLibrary, shapeLibrary, Δt_rf, i)
    #Unpacking
    #(1)amplitude (2)mag_id (3)phase_id (4)time_shape_id (5)delay (6)freq (7)phase
    r = rfLibrary[i]["data"]
    amplitude =     r[1]
    mag_id =        r[2] |> x->floor(Int64,x)
    phase_id =      r[3] |> x->floor(Int64,x)
    time_shape_id = r[4] |> x->floor(Int64,x)
    delay =         r[5] + (time_shape_id==0)*Δt_rf/2
    freq =          r[6]
    phase =         r[7]
    #Amplitude and phase waveforms
    if amplitude != 0 && mag_id != 0
        rfA = decompress_shape(shapeLibrary[mag_id]...)[1:end-1]
        rfϕ = decompress_shape(shapeLibrary[phase_id]...)[1:end-1]
        Nrf = shapeLibrary[mag_id][1] - 1
        rfAϕ = amplitude .* rfA .* exp.(-1im*(2π*rfϕ .+ phase))
    else
        rfA = 0
        rfϕ = 0
        Nrf = 0
        rfAϕ = 0
    end
    #Creating timings
    if time_shape_id == 0 #no time waveform
        rfT = Nrf * Δt_rf
    else
        rft = decompress_shape(shapeLibrary[time_shape_id]...)
        rfT = (rft[2:end] .- rft[1:end-1]) * Δt_rf
    end
    R = [RF(rfAϕ,rfT,freq,delay);;]
    R
end

function read_ADC(adcLibrary, i)
    #Unpacking
    #(1)num (2)dwell (3)delay (4)freq (5)phase
    a = adcLibrary[i]["data"]
    num =   a[1] |> x->floor(Int64,x)
    dwell = a[2]
    delay = a[3] + dwell/2
    freq =  a[4]
    phase = a[5]
    #Definition
    T = (num-1) * dwell
    A = [ADC(num,T,delay,freq,phase)]
    A
end

function get_block(obj, i)
    #Unpacking
    idelay, irf, ix, iy, iz, iadc, iext = obj["blockEvents"][i]

    #Gradient definition
    Δt_gr = obj["definitions"]["GradientRasterTime"]
    Gx = read_Grad(obj["gradLibrary"], obj["shapeLibrary"], Δt_gr, ix)
    Gy = read_Grad(obj["gradLibrary"], obj["shapeLibrary"], Δt_gr, iy)
    Gz = read_Grad(obj["gradLibrary"], obj["shapeLibrary"], Δt_gr, iz)
    G = [Gx;Gy;Gz;;]

    #RF definition
    Δt_rf = obj["definitions"]["RadiofrequencyRasterTime"]
    R = read_RF(obj["rfLibrary"], obj["shapeLibrary"], Δt_rf, irf)

    #ADC definition
    A = read_ADC(obj["adcLibrary"], iadc)

    #DUR
    D = [obj["blockDurations"][i]]

    #Sequence block definition
    s = Sequence(G,R,A,D)
    s
end


text:
#=
module Gaugefields_4D_module
    using LinearAlgebra
    import ..AbstractGaugefields_module:AbstractGaugefields,Shifted_Gaugefields,shift_U,
                        Adjoint_Gaugefields,set_wing_U!,Abstractfields,construct_staple!,clear_U!,
                        calculate_Plaquette
    import Base
    =#

    abstract type Gaugefields_4D{NC} <: AbstractGaugefields{NC,4}
    end

    include("./gaugefields_4D_wing.jl")
    include("./gaugefields_4D_nowing.jl")

    function __init__()
        @require MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" begin
            include("./gaugefields_4D_mpi.jl")
            include("./gaugefields_4D_mpi_nowing.jl")
            include("../2D/gaugefields_2D_mpi_nowing.jl")
            include("./TA_gaugefields_4D_mpi.jl")
            include("../2D/TA_gaugefields_2D_mpi.jl")
        end
    end

    function Base.size(U::Gaugefields_4D{NC}) where NC
        return NC,NC,U.NX,U.NY,U.NZ,U.NT
    end

    function clear_U!(U::Array{T,1}) where T <: Gaugefields_4D
        for μ=1:4
            clear_U!(U[μ])
        end

    end

    #=
    function calculate_Plaquet(U::Array{T,1}) where T <: Gaugefields_4D
        error("calculate_Plaquet is not implemented in type $(typeof(U)) ")
    end
    =#

#end


text:
using ComponentArrays
using DifferentialEquations
using SimulationLogs
using Plots; gr()
using UnPack

## Simulation functions
# Simple car with velocity-square drag
function car!(D, vars, p, t; u=0.0)
    @unpack vel = vars
    @unpack c, m = p

    @log drag = c*vel^2

    D.pos = vel
    D.vel = (-drag*sign(vel) + u)/m
end

# Proportional-integral control
function PI!(D, vars, p, t; e)
    @unpack ∫e = vars
    @unpack kp, ki = p

    @log u = kp*e + ki*∫e

    D.∫e = e
    return u
end

# Car with cruise control
function feedback_car!(D, vars, p, t)
    @log r = p.ref(t)
    @log e = r - vars.car.vel

    u = PI!(D.control, vars.control, p.control, t; e=e)
    car!(D.car, vars.car, p.car, t; u)
end

## Inputs
# Parameters
p = (
    ref = t -> 10*(t>1),
    car = (
        m = 1000,
        c = 5,
    ),
    control = (
        kp = 800,
        ki = 40,
    )
)

# Initial conditions
ic = ComponentArray(
    car = (
        pos = 0.0,
        vel = 0.0,
    ),
    control = (
        ∫e = 0.0,
    ),
)

## Setup problem and solve
prob = ODEProblem(feedback_car!, ic, (0.0, 20.0), p)
sol = solve(prob)

## Plots
p1 = plot(sol, vars=2, title="Reference Tracking")
scope!(sol, :r)

p2 = scope(sol, [:u, :drag], title="Forces")

p3 = plot(sol, vars=3, title="Errors")
scope!(sol, :e)

plot(p1, p2, p3, layout=(3,1), size=(600,800))

max_stars_repo_name: jonniedie/SimulationLogs

max_issues_repo_name: jonniedie/SimulationLogs

max_forks_repo_name: jonniedie/SimulationLogs


text:
function compact_string(x)
    io = IOContext(IOBuffer(), :compact => true)
    print(io, x)
    return String(take!(io.io))
end

function compact_string_row(x::Integer)
    io = IOContext(IOBuffer(), :compact => true)
    print(io, x)
    return String(take!(io.io))
end

function compact_string_row(x::CartesianIndex{2})
    io = IOContext(IOBuffer(), :compact => true)
    print(io, x.I[1], ",", x.I[2])
    return String(take!(io.io))
end

function truncate_string(str, max_width)
    string_width = textwidth(str)
    if string_width <= max_width
        return str
    end

    newstring = ""
    string_width = 0
    for grapheme ∈ graphemes(str)
        grapheme_width = textwidth(grapheme)
        if string_width + grapheme_width >= max_width
            return newstring * "…"
        end
        newstring = newstring * grapheme
        string_width += grapheme_width
    end
    # Shouldn't get here, given that textwidth(s) = reduce(charwidth, 0, s)
end

# Two columns, for Series
function balance_widths(width1::Int, width2::Int, max_width::Int, min_width::Int = 20)
    @assert max_width >= 0
    biggest_width = max(width1, width2)

    while width1 + width2 > max_width
        biggest_width -= 1
        if biggest_width < min_width
            break
        end
        width1 = min(biggest_width, width1)
        width2 = min(biggest_width, width2)
    end

    return (width1, width2)
end

# Multiple columns, for Table
function balance_widths!(widths::Vector{Int}, ncols::Int, max_width::Int, buffer::Int = 2, min_width::Int = 20)
    @assert max_width >= 0
    biggest_width = maximum(widths)

    # First pass: shorten all the column widths
    while sum(widths) + buffer*length(widths) > max_width
        biggest_width -= 1
        if biggest_width < min_width
            break
        end
        for i = 1:length(widths)
            @inbounds widths[i] = min(biggest_width, widths[i])
        end
    end

    # Second pass: truncate columns, if necessary
    if sum(widths) + buffer*length(widths) > max_width
        total_width = 0
        for i = 1:length(widths)
            width = widths[i] + buffer
            if i == length(widths)
                ncols_shown = i-1
            else
                if total_width + width + buffer + 1 > max_width # we anticipate at least a "⋯" column to follow this one
                    ncols_shown = i-1
                    break
                end
            end
            total_width += width
        end
    else
        ncols_shown = length(widths)
    end

    return ncols_shown
end

show(io::IO, ::MIME"text/plain", t::Table) = showtable(io, t)
show(io::IO, t::Table) = showtable(io, t)
show(io::IO, ::MIME"text/plain", t::ProductTable) = showtable(io, t)
show(io::IO, t::ProductTable) = showtable(io, t)
function showtable(io::IO, @nospecialize t)
    row_inds = keys(t)
    col_inds = colnames(t)
    nrows = length(row_inds)::Int
    ncols = length(col_inds)::Int
    display_width = displaysize(io)[2]
    max_cols = display_width ÷ 3 # assuming one-width columns with two spaces in-between

    typename = typeof(t).name.name
    print(io, "$typename with $ncols column$(ncols == 1 ? "" : "s") and $nrows row$(nrows == 1 ? "" : "s")")

    n_indexes = 0
    if t isa Table
        n_indexes = length(accelerators(t))
    elseif t isa ProductTable
        n_indexes = length(accelerators(t.t1)) + length(accelerators(t.t2))
    end
    if n_indexes > 0
        print(io, " and $n_indexes accelerator$(n_indexes == 1 ? "" : "s")")
    end

    max_show_rows = displaysize(io)[1] - 7

    ncols_shown = min(ncols+1, max_cols)::Int # First "column" shown is the indices

    strings = Vector{Vector{String}}(undef, ncols_shown)
    for i ∈ 1:ncols_shown
        for j ∈ 1:(min(max_show_rows, nrows) + 1)
            if j == 1
                if i == 1
                    strings[i] = [""]
                else
                    strings[i] = [compact_string(col_inds[i-1])]
                end
            else
                if i == 1
                    push!(strings[i], compact_string_row(row_inds[j-1]))
                else
                    push!(strings[i], compact_string(t[row_inds[j-1]][col_inds[i-1]]))
                end
            end
        end
    end

    max_column_widths = [max(1, maximum(textwidth, str_vec)) for str_vec ∈ strings]

    if sum(max_column_widths) + 2*length(max_column_widths) > display_width
        # Shorten each column and reduce the number of columns shown
        # If they all still don't fit, the final shown column will be filled with "⋯"

        ncols_shown = balance_widths!(max_column_widths, ncols, display_width)
        strings = strings[1:ncols_shown]
        max_column_widths = max_column_widths[1:ncols_shown]
        for i = 1:length(strings)
            map!(str -> truncate_string(str, max_column_widths[i]), strings[i], strings[i])
        end
    end

    if ncols_shown < ncols+1
        push!(strings, fill("⋯", length(strings[1])))
        push!(max_column_widths, 1)
    end

    if nrows > max_show_rows
        for i ∈ 1:length(strings)
            if i == length(strings) && strings[i][1] == "⋯"
                push!(strings[i], (" " ^ div(max_column_widths[i]-1, 2)) *  "⋱")
            else
                push!(strings[i], (" " ^ div(max_column_widths[i]-1, 2)) *  "⋮")
            end
        end
    end

    # Now produce output

    # Header: " --  Column1  Column2"
    print(io, ":\n ")
    n_spaces = max_column_widths[1] + 3
    print(io, " " ^ n_spaces)
    for i ∈ 2:length(strings)
        @inbounds col_str = strings[i][1]
        print(io, col_str)
        if i != length(strings)
            n_spaces = max_column_widths[i] - textwidth(col_str) + 2
            print(io, " " ^ n_spaces)
        end
    end

    # Seperator: " ┌────────"
    print(io, "\n ")
    print(io, " " ^ max_column_widths[1])
    print(io, " ┌─")
    print(io, "─" ^ (sum(max_column_widths) - max_column_widths[1] + 2*(length(max_column_widths) - 2)))

    # Body " rowind │ val1  val2"
    for j = 2:length(strings[1])
        print(io, "\n ")
        @inbounds row_str = strings[1][j]
        print(io, row_str)
        n_spaces = max_column_widths[1] - textwidth(row_str)
        if n_spaces > 0
            print(io, " " ^ n_spaces)
        end
        print(io, " │ ")

        for i = 2:length(strings)
            @inbounds val_str = strings[i][j]
            print(io, val_str)

            if i != length(strings)
                n_spaces = max_column_widths[i] - textwidth(val_str) + 2
                print(io, " " ^ n_spaces)
            end
        end
    end
end