text:
########################################################
####                  cehck syntax                  ####
########################################################

#------ truncated-normal ---------

function checksyn(::Type{Trun})
    for k in (:depvar, :frontier, :μ, :σᵤ², :σᵥ²)
        if  (_dicM[k] === nothing)
            throw("For the truncated-normal model, the $k equation is missing in sfmodel_spec().")
        end
    end

    for k in (:hscale, :timevar, :σₐ², :idvar, :gamma)
        if  (_dicM[k] !== nothing)
            throw("For the truncated-normal model, the $k equation is not needed in sfmodel_spec().")
        end
    end
end

 #------- half-normal -------------

function checksyn(::Type{Half})
    for k in (:depvar, :frontier, :σᵤ², :σᵥ²)
        if  (_dicM[k] === nothing)
            throw("For the half-normal model, the $k equation is missing in sfmodel_spec().")
        end
    end

    for k in (:μ, :hscale, :timevar, :σₐ², :idvar, :gamma)
        if  (_dicM[k] !== nothing)
            throw("For the half-normal model, the $k equation is not needed in sfmodel_spec().")
        end
    end
end

# ------- Exponential ----------------

function checksyn(::Type{Expo})
    for k in (:depvar, :frontier, :σᵤ², :σᵥ²)
        if  (_dicM[k] === nothing)
            throw("For the exponential-normal model, the $k equation is missing in sfmodel_spec().")
        end
    end

    for k in (:μ, :hscale, :timevar, :σₐ², :idvar, :gamma )
        if  (_dicM[k] !== nothing)
            throw("For the exponential-normal model, the $k equation is not needed in sfmodel_spec().")
        end
    end
end

# ------- scaling property -------------

function checksyn(::Type{Trun_Scale})
    for k in (:depvar, :frontier, :hscale, :μ, :σᵤ², :σᵥ²)
        if  (_dicM[k] === nothing)
            throw("For the scaling-property model, the $k equation is required but is currently missing in your sfmodel_spec().")
        end
    end
    for k in ( :timevar, :idvar, :σₐ², :gamma)
        if  (_dicM[k] !== nothing)
            throw("For the scaling-property model, the $k equation is not needed in sfmodel_spec().")
        end
    end
end

# ------- panel Wang and Ho 2020, with Truncated Dist -------------

function checksyn(::Type{PFEWHT})
    for k in (:depvar, :frontier, :hscale, :μ, :σᵤ², :σᵥ², :timevar, :idvar)
        if  (_dicM[k] === nothing)
            throw("For the panel fixed effect model of Wang and Ho (2010, JE) with the truncatd-normal assumption, the $k equation is required but is currently missing in your sfmodel_spec().")
        end
    end
    for k in (:gamma, :σₐ²)
        if  (_dicM[k] !== nothing)
            throw("For the panel fixed effect model of Wang and Ho (2010, JE) with the truncated-normal assumption, the $k equation is not needed in sfmodel_spec().")
        end
    end
end

# ------- panel Wang and Ho 2020, with Halfnormal Dist -------------

function checksyn(::Type{PFEWHH})
    for k in (:depvar, :frontier, :hscale, :σᵤ², :σᵥ², :timevar, :idvar)
        if  (_dicM[k] === nothing)
            throw("For the panel fixed effect model of Wang and Ho (2010, JE) with the half-normal assumption, the $k equation is required but is currently missing in your sfmodel_spec().")
        end
    end
    for k in (:μ, :σₐ², :gamma)
        if  (_dicM[k] !== nothing)
            throw("For the panel fixed effect model of Wang and Ho (2010, JE) with the half-normal assumption, the $k equation is not needed in sfmodel_spec().")
        end
    end
end

# ------- panel Time Decay of BC1992 -------------

function checksyn(::Type{PanDecay})
    for k in (:depvar, :frontier, :μ, :gamma, :σᵤ², :σᵥ², :idvar)
        if  (_dicM[k] === nothing)
            throw("For the panel time-decay model of Battese and Coelli (1992), the $k equation is required but is currently missing in your sfmodel_spec().")
        end
    end
    for k in (:σₐ², :hscale, )
        if  (_dicM[k] !== nothing)
            throw("For the panel time-decay model of Battese and Coelli (1992), the $k equation is not needed in sfmodel_spec().")
        end
    end
end

# ------- panel Kumbhakar 1990 model  -------------

function checksyn(::Type{PanKumb90})
    for k in (:depvar, :frontier, :μ, :gamma, :σᵤ², :σᵥ², :idvar)
        if  (_dicM[k] === nothing)
            throw("For the panel Kumbhakar (1990) model, the $k equation is required but is currently missing in your sfmodel_spec().")
        end
    end
    for k in (:σₐ², :hscale, )
        if  (_dicM[k] !== nothing)
            throw("For the panel Kumbhakar (1990) model, the $k equation is not needed in sfmodel_spec().")
        end
    end
end

# ------- panel CSW (2014 JoE) CSN model, with Halfnormal Dist ----------

function checksyn(::Type{PFECSWH})
    for k in (:depvar, :frontier, :σᵤ², :σᵥ², :timevar, :idvar)
        if  (_dicM[k] === nothing)
            throw("For the panel fixed effect model of CSW (2014, JE) with the half-normal assumption, the $k equation is required but is currently missing in your sfmodel_spec().")
        end
    end
    for k in (:μ, :gamma, :σₐ², :hscale)
        if  (_dicM[k] !== nothing)
            throw("For the panel fixed effect model of CSW (2014, JE) with the half-normal assumption, the $k equation is not needed in sfmodel_spec().")
        end
    end
end

# ------- panel TRE with Halfnormal Dist ----------

function checksyn(::Type{PTREH})
    for k in (:depvar, :frontier, :σₐ², :σᵤ², :σᵥ², :timevar, :idvar)
        if  (_dicM[k] === nothing)
            throw("For the panel true random effect model of Greene (2005, JoE) with the half-normal assumption, the $k equation is required but is currently missing in your sfmodel_spec().")
        end
    end
    for k in (:μ, :gamma, :hscale)
        if  (_dicM[k] !== nothing)
            throw("For the panel true random effect model of Greene (2005, JoE) with the half-normal assumption, the $k equation is not needed in sfmodel_spec().")
        end
    end
end

# ------- panel TRE with truncated normal dist ----------

function checksyn(::Type{PTRET})
    for k in (:depvar, :frontier, :μ, :σₐ², :σᵤ², :σᵥ², :timevar, :idvar)
        if  (_dicM[k] === nothing)
            throw("For the panel true random effect model of Greene (2005, JoE) with the truncated-normal assumption, the $k equation is required but is currently missing in your sfmodel_spec().")
        end
    end
    for k in ( :gamma, :hscale)
        if  (_dicM[k] !== nothing)
            throw("For the panel true random effect model of Greene (2005, JoE) with the truncated-normal assumption, the $k equation is not needed in sfmodel_spec().")
        end
    end
end

########################################################
####               check constant                   ####
########################################################

macro requireConst(arg)
    return arg
end

    # ---- check overall constant --------------

function checkConst(mat::Array, S::Symbol, requireConst::Int)
    for i in 1:size(mat,2)
        aa = mat[:,i]
        # varstd = sqrt(sum((aa .- sum(aa)/length(aa)).^2)/length(aa)) ; then varstd <= 1e-8
        bb = unique(aa)
        if requireConst == 1
            (length(bb) == 1) || throw("The variable in the $(S) function has to be a constant for this model.")
        else
            (length(bb)  > 1) || throw("The $(S) function cannot have a constant for this model.")
        end
    end
end

     # ---- check constant panel-by-panel -----------

function checkConst(mat::Array, S::Symbol, requireConst::Int, rowid::Vector)
    for i in 1:size(mat,2)
        for j in 1:size(rowid,1)
            aa = mat[rowid[j], i]
            # varstd = sqrt(sum((aa .- sum(aa)/length(aa)).^2)/length(aa))
            bb = unique(aa)
            if requireConst == 1
               (length(bb) == 1) || throw("The variable in the $(S) function has to be a constant within a panel for this model.")
            else
               (length(bb)  > 1)  || throw("The $(S) function cannot have a constant within a panel for this model.")
            end
        end
    end
end

########################################################
####          Chi-equare related tables             ####
########################################################

#* --------- mixed Chi-square table ---------
"""
    sfmodel_MixTable(<keyword arguments>)

Display and return critical values of the mixed `χ²` (ch-square) distribution.
The values are taken from Table 1, Kodde and Palm (1986, Econometrica).

## Argument
- `dof::Integer`: the degree of freedom. Currently support `dof` between 1 and 40.

## Examples
```julia
julia> sfmodel_MixTable(3)

  * Significance levels and critical values of the mixed χ² distribution
┌─────┬───────┬───────┬───────┬────────┐
│ dof │  0.10 │  0.05 │ 0.025 │   0.01 │
├─────┼───────┼───────┼───────┼────────┤
│ 3.0 │ 5.528 │ 7.045 │ 8.542 │ 10.501 │
└─────┴───────┴───────┴───────┴────────┘

source: Table 1, Kodde and Palm (1986, Econometrica).
```

"""
function sfmodel_MixTable(dof::Real=9999) # not using ::Int64 in order to thow informative error

    (dof isa(Integer)) || throw("The degree of freedom (dof) has to be an integer.")

    @isdefined(_dicOPT) || sfmodel_opt()
    sf_table = _dicOPT[:table_format]

    table = Array{Float32}(undef, 40, 5)

    table[1 ,1] = 1;   table[1 ,2] =   1.642; table[1 ,3] =   2.705;  table[1 ,4] =   3.841;   table[1 ,5]=    5.412 ;
    table[2 ,1] = 2;   table[2 ,2] =   3.808; table[2 ,3] =   5.138;  table[2 ,4] =   6.483;   table[2 ,5]=    8.273 ;
    table[3 ,1] = 3;   table[3 ,2] =   5.528; table[3 ,3] =   7.045;  table[3 ,4] =   8.542;   table[3 ,5]=   10.501 ;
    table[4 ,1] = 4;   table[4 ,2] =   7.094; table[4 ,3] =   8.761;  table[4 ,4] =  10.383;   table[4 ,5]=   12.483 ;
    table[5 ,1] = 5;   table[5 ,2] =   8.574; table[5 ,3] =  10.371;  table[5 ,4] =   12.103;  table[5 ,5] =  14.325 ;
    table[6 ,1] = 6;   table[6 ,2] =   9.998; table[6 ,3] =  11.911;  table[6 ,4] =   13.742;  table[6 ,5] =  16.704 ;
    table[7 ,1] = 7;   table[7 ,2] =  11.383; table[7 ,3] =  13.401;  table[7 ,4] =   15.321;  table[7 ,5] =  17.755 ;
    table[8 ,1] = 8;   table[8 ,2] =  12.737; table[8 ,3] =  14.853;  table[8 ,4] =   16.856;  table[8 ,5] =  19.384 ;
    table[9 ,1] = 9;   table[9 ,2] =  14.067; table[9 ,3] =  16.274;  table[9 ,4] =   18.354;  table[9 ,5] =  20.972 ;
    table[10,1] = 10;  table[10,2] =  15.377; table[10,3] =  17.670;  table[10,4] =   19.824;  table[10,5] =  22.525 ;
    table[11,1] = 11;  table[11,2] =  16.670; table[11,3] =  19.045;  table[11,4] =   21.268;  table[11,5] =  24.049 ;
    table[12,1] = 12;  table[12,2] =  17.949; table[12,3] =  20.410;  table[12,4] =   22.691;  table[12,5] =  25.549 ;
    table[13,1] = 13;  table[13,2] =  19.216; table[13,3] =  21.742;  table[13,4] =   24.096;  table[13,5] =  27.026 ;
    table[14,1] = 14;  table[14,2] =  20.472; table[14,3] =  23.069;  table[14,4] =   25.484;  table[14,5] =  28.485 ;
    table[15,1] = 15;  table[15,2] =  21.718; table[15,3] =  24.384;  table[15,4] =   26.856;  table[15,5] =  29.927 ;
    table[16,1] = 16;  table[16,2] =  22.956; table[16,3] =  25.689;  table[16,4] =   28.219;  table[16,5] =  31.353 ;
    table[17,1] = 17;  table[17,2] =  24.186; table[17,3] =  26.983;  table[17,4] =   29.569;  table[17,5] =  32.766 ;
    table[18,1] = 18;  table[18,2] =  25.409; table[18,3] =  28.268;  table[18,4] =   30.908;  table[18,5] =  34.167 ;
    table[19,1] = 19;  table[19,2] =  26.625; table[19,3] =  29.545;  table[19,4] =   32.237;  table[19,5] =  35.556 ;
    table[20,1] = 20;  table[20,2] =  27.835; table[20,3] =  30.814;  table[20,4] =   33.557;  table[20,5] =  36.935 ;
    table[21,1] = 21;  table[21,2] =  29.040; table[21,3] =  32.077;  table[21,4] =   34.869;  table[21,5] =  38.304 ;
    table[22,1] = 22;  table[22,2] =  30.240; table[22,3] =  33.333;  table[22,4] =   36.173;  table[22,5] =  39.664 ;
    table[23,1] = 23;  table[23,2] =  31.436; table[23,3] =  34.583;  table[23,4] =   37.470;  table[23,5] =  41.016 ;
    table[24,1] = 24;  table[24,2] =  32.627; table[24,3] =  35.827;  table[24,4] =   38.761;  table[24,5] =  42.360 ;
    table[25,1] = 25;  table[25,2] =  33.813; table[25,3] =  37.066;  table[25,4] =   40.045;  table[25,5] =  43.696 ;
    table[26,1] = 26;  table[26,2] =  34.996; table[26,3] =  38.301;  table[26,4] =   41.324;  table[26,5] =  45.026 ;
    table[27,1] = 27;  table[27,2] =  36.176; table[27,3] =  39.531;  table[27,4] =   42.597;  table[27,5] =  46.349 ;
    table[28,1] = 28;  table[28,2] =  37.352; table[28,3] =  40.756;  table[28,4] =   43.865;  table[28,5] =  47.667 ;
    table[29,1] = 29;  table[29,2] =  38.524; table[29,3] =  41.977;  table[29,4] =   45.128;  table[29,5] =  48.978 ;
    table[30,1] = 30;  table[30,2] =  39.694; table[30,3] =  43.194;  table[30,4] =   46.387;  table[30,5] =  50.284 ;
    table[31,1] = 31;  table[31,2] =  40.861; table[31,3] =  44.408;  table[31,4] =   47.641;  table[31,5] =  51.585 ;
    table[32,1] = 32;  table[32,2] =  42.025; table[32,3] =  45.618;  table[32,4] =   48.891;  table[32,5] =  52.881 ;
    table[33,1] = 33;  table[33,2] =  43.186; table[33,3] =  46.825;  table[33,4] =   50.137;  table[33,5] =  54.172 ;
    table[34,1] = 34;  table[34,2] =  44.345; table[34,3] =  48.029;  table[34,4] =   51.379;  table[34,5] =  55.459 ;
    table[35,1] = 35;  table[35,2] =  45.501; table[35,3] =  49.229;  table[35,4] =   52.618;  table[35,5] =  56.742 ;
    table[36,1] = 36;  table[36,2] =  46.655; table[36,3] =  50.427;  table[36,4] =   53.853;  table[36,5] =  58.020 ;
    table[37,1] = 37;  table[37,2] =  47.808; table[37,3] =  51.622;  table[37,4] =   55.085;  table[37,5] =  59.295 ;
    table[38,1] = 38;  table[38,2] =  48.957; table[38,3] =  52.814;  table[38,4] =   56.313;  table[38,5] =  60.566 ;
    table[39,1] = 39;  table[39,2] =  50.105; table[39,3] =  54.003;  table[39,4] =   57.539;  table[39,5] =  61.833 ;
    table[40,1] = 40;  table[40,2] =  51.251; table[40,3] =  55.190;  table[40,4] =   58.762;  table[40,5] =  63.097 ;

    println()
    printstyled("  * Significance levels and critical values of the mixed χ² distribution\n"; color = :yellow)

    pretty_table(if 1<=dof<=40       # correct instance
                    table[dof,:]'
                 elseif dof == 9999  # user did not specify dof
                    printstyled("  * Use `sfmodel_MixTable(d)` to show the specific critical values for the d.o.f. equal to `d`.\n"; color = :red)
                    table            # show all
                 elseif (dof < 1) || (dof > 40)
                    throw("The table only allows the degree of freedom between 1 and 40.")
                 else
                    throw("There is a problem in your specification of the degree of freedom.")
                 end,
                 header=["dof", "0.10", "0.05", "0.025", "0.01"],
                 formatters = ft_printf("%2.3f", 2:5),
                 compact_printing = true,
                 backend = Val(sf_table))
    println()
    printstyled("source: Table 1, Kodde and Palm (1986, Econometrica)."; color = :yellow)
    println()

    if 1<=dof<=40
        return table[dof, :]'
    elseif dof == 9999
        return table
    end
end

#* ---------  Chi-square table ---------

function sfmodel_ChiSquareTable(dof::Real=1) # not using ::Int64 in order to throw informative error

    (dof isa(Integer)) || throw("The degree of freedom (dof) has to be an integer.")

    table = Array{Float64}(undef, 1, 5)

    table[1, 1] = dof
    table[1, 2] = chisqinvcdf(dof, 1-0.10)
    table[1, 3] = chisqinvcdf(dof, 1-0.05)
    table[1, 4] = chisqinvcdf(dof, 1-0.025)
    table[1, 5] = chisqinvcdf(dof, 1-0.01)

    println()
    printstyled("  * significance levels and critical values of the χ² distribution\n"; color = :yellow)

    pretty_table(table,
                 ["dof", "0.10", "0.05", "0.025", "0.01"],
                 formatters = ft_printf( "%2.3f", 2:5),
                 compact_printing = true,
                 backend = Val(sf_table))

end

#* ---- check multi-collinearity in the data --------

function checkCollinear(modelid, xvar, zvar, qvar, wvar, vvar)

   # check multicollinearity in the dataset (matrix)

    for i in 1:5

        if i == 1
            themat = xvar
            eqname = " `frontier` "
            eqsym  = :frontier
        elseif i == 2
            themat = zvar
            eqname = " `μ` "
            eqsym  = :μ
        elseif i == 3
            themat = qvar
            if modelid == PanDecay
                eqname = " `gamma` "
                eqsym  = :gamma
            elseif modelid == TRE
                eqname = " `σₐ²` "
                eqsym  = :σₐ²
            else
                eqname = " `hscale` "
                eqsym  = :hscale
            end
        elseif i == 4
            themat = wvar
            eqname = " `σᵤ²` "
            eqsym  =  :σᵤ²
        elseif i == 5
            themat = vvar
            eqname = " `σᵥ²` "
            eqsym  = :σᵥ²
        end

        if length(themat) > 0 && size(themat, 2) > 1  # check only if have more than 1 var

            _, pivots = rref_with_pivots(themat)  # pivots has unique,
                                                  # non-collinear columns;
                                                  # require `using RowEchelon`

            if length(pivots) != size(themat, 2) # number of columns are different, meaning collinear is dropped
                theColl = filter(x -> x ∉ pivots, 1:size(themat, 2)) # positions of problematic columns
                bb = ""
                for i in theColl  # get problematic var's name from the equation
                    bb = bb*String(_dicM[eqsym][i])*String(", ")
                end
                printstyled(bb, "appear to have perfect collinearity with other variables in the equation", eqname, ". Try dropping this/these variable(s) from the equation.\n"; color = :red)
                throw("Multi-collinearity in the variables.")
            end  # if length(pivots)
        end # if size(themat, 2) > 1
    end # for i in 1:5
  end


text:
include("euler/euler.jl")
using .Numbers: count_divisors
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    i, triangle = 1, 1
    while count_divisors(triangle) < n
        i += 1
        triangle += i
    end
    return triangle
end

compute(5)

compute(500)

@benchmark compute(500)

ext: jl

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using SampledSignals

#####
##### `SynthesizerSource`
#####

struct SynthesizerSource{T<:AbstractSynthComponent} <: SampleSource
    output::T
end

function SampledSignals.unsafe_read!(source::SynthesizerSource, buf::Array, frameoffset, framecount)
    for i in (1:framecount) .+ frameoffset
        buf[i] = next!(source.output)
    end
    return framecount
end

SampledSignals.samplerate(::SynthesizerSource) = SAMPLE_RATE[]
SampledSignals.nchannels(::SynthesizerSource) = 1
SampledSignals.eltype(::SynthesizerSource{T}) where {T} = output_type(T)


text:
mutable struct LinprogSolution
    status
    objval
    sol
    attrs
end

InputVector{T<:Union{Real,Char}} = Union{Vector{T},Real,Char}

function expandvec(x,len::Integer)
    if isa(x,Vector)
        if length(x) != len
            error("Input size mismatch. Expected vector of length $len but got $(length(x))")
        end
        return x
    else
        return fill(x,len)
    end
end

function no_lp_solver()
    error("""
          A linear programming solver must be specified as the final argument.

          You can use any solver listed on the solvers table of http://www.juliaopt.org/ that has a checkmark in the Linear/Quadratic column.

          A (free) example is Clp.jl. Once Clp is installed and imported via "using Clp", you can specify ClpSolver() as the solver.
          Solver options are specified by using keyword arguments to ClpSolver().
          """)
end

function buildlp(c::InputVector, A::AbstractMatrix, rowlb::InputVector, rowub::InputVector, lb::InputVector, ub::InputVector, solver::AbstractMathProgSolver)
    m = LinearQuadraticModel(solver)
    nrow,ncol = size(A)

    c = expandvec(c, ncol)
    rowlbtmp = expandvec(rowlb, nrow)
    rowubtmp = expandvec(rowub, nrow)
    lb = expandvec(lb, ncol)
    ub = expandvec(ub, ncol)

    # rowlb is allowed to be vector of senses
    if eltype(rowlbtmp) == Char
        realtype = eltype(rowubtmp)
        warn_no_inf(realtype)
        sense = rowlbtmp
        rhs = rowubtmp
        @assert realtype <: Real

        rowlb = Array{realtype}(undef, nrow)
        rowub = Array{realtype}(undef, nrow)
        for i in 1:nrow
            if sense[i] == '<'
                rowlb[i] = typemin(realtype)
                rowub[i] = rhs[i]
            elseif sense[i] == '>'
                rowlb[i] = rhs[i]
                rowub[i] = typemax(realtype)
            elseif sense[i] == '='
                rowlb[i] = rhs[i]
                rowub[i] = rhs[i]
            else
                error("Unrecognized sense '$(sense[i])'")
            end
        end
    else
        rowlb = rowlbtmp
        rowub = rowubtmp
    end

    loadproblem!(m, A, lb, ub, c, rowlb, rowub, :Min)

    return m
end

function solvelp(m)
    optimize!(m)
    stat = status(m)
    if stat == :Optimal
        attrs = Dict()
        attrs[:redcost] = getreducedcosts(m)
        attrs[:lambda] = getconstrduals(m)
        return LinprogSolution(stat, getobjval(m), getsolution(m), attrs)
    elseif stat == :Infeasible
        attrs = Dict()
        try
            attrs[:infeasibilityray] = getinfeasibilityray(m)
        catch
            @warn("Problem is infeasible, but infeasibility ray (\"Farkas proof\") is unavailable; check that the proper solver options are set.")
        end
        return LinprogSolution(stat, nothing, [], attrs)
    elseif stat == :Unbounded
        attrs = Dict()
        try
            attrs[:unboundedray] = getunboundedray(m)
        catch
            @warn("Problem is unbounded, but unbounded ray is unavailable; check that the proper solver options are set.")
        end
        return LinprogSolution(stat, nothing, [], attrs)
    else
        return LinprogSolution(stat, nothing, [], Dict())
    end
end

function linprog(c::InputVector, A::AbstractMatrix, rowlb::InputVector, rowub::InputVector, lb::InputVector, ub::InputVector, solver::AbstractMathProgSolver)
    m = buildlp(c, A, rowlb, rowub, lb, ub, solver)
    return solvelp(m)
end

linprog(c,A,rowlb,rowub, solver::AbstractMathProgSolver) = linprog(c,A,rowlb,rowub,0,Inf, solver)

buildlp(c,A,rowlb,rowub, solver::AbstractMathProgSolver) = buildlp(c,A,rowlb,rowub,0,Inf, solver)

# Old versions
linprog(c::InputVector, A::AbstractMatrix, rowlb::InputVector, rowub::InputVector, lb::InputVector, ub::InputVector) = no_lp_solver()
linprog(c,A,rowlb,rowub) = no_lp_solver()

export linprog, buildlp, solvelp


text:
"""
    MOInputIsotopicByFeatures(x::AbstractVector, out_dim::Integer)

`MOInputIsotopicByFeatures(x, out_dim)` has length `out_dim * length(x)`.

```jldoctest
julia> x = [1, 2, 3];

julia> KernelFunctions.MOInputIsotopicByFeatures(x, 2)
6-element KernelFunctions.MOInputIsotopicByFeatures{Int64, Vector{Int64}}:
 (1, 1)
 (1, 2)
 (2, 1)
 (2, 2)
 (3, 1)
 (3, 2)
```

Accommodates modelling multi-dimensional output data where all outputs are always observed.

As shown above, an `MOInputIsotopicByFeatures` represents a vector of tuples.
The first `out_dim` elements represent all outputs for the first input, the second
`out_dim` elements represent the outputs for the second input, etc.

See [Inputs for Multiple Outputs](@ref) in the docs for more info.
"""
struct MOInputIsotopicByFeatures{S,T<:AbstractVector{S}} <: AbstractVector{Tuple{S,Int}}
    x::T
    out_dim::Integer
end

"""
    MOInputIsotopicByOutputs(x::AbstractVector, out_dim::Integer)

`MOInputIsotopicByOutputs(x, out_dim)` has length `length(x) * out_dim`.

```jldoctest
julia> x = [1, 2, 3];

julia> KernelFunctions.MOInputIsotopicByOutputs(x, 2)
6-element KernelFunctions.MOInputIsotopicByOutputs{Int64, Vector{Int64}}:
 (1, 1)
 (2, 1)
 (3, 1)
 (1, 2)
 (2, 2)
 (3, 2)
```

Accommodates modelling multi-dimensional output data where all outputs are always observed.

As shown above, an `MOInputIsotopicByOutputs` represents a vector of tuples.
The first `length(x)` elements represent the inputs for the first output, the second
`length(x)` elements represent the inputs for the second output, etc.
"""
struct MOInputIsotopicByOutputs{S,T<:AbstractVector{S}} <: AbstractVector{Tuple{S,Int}}
    x::T
    out_dim::Integer
end

const IsotopicMOInputsUnion = Union{MOInputIsotopicByFeatures,MOInputIsotopicByOutputs}

function Base.getindex(inp::MOInputIsotopicByOutputs, ind::Integer)
    @boundscheck checkbounds(inp, ind)
    output_index, feature_index = fldmod1(ind, length(inp.x))
    feature = @inbounds inp.x[feature_index]
    return feature, output_index
end

function Base.getindex(inp::MOInputIsotopicByFeatures, ind::Integer)
    @boundscheck checkbounds(inp, ind)
    feature_index, output_index = fldmod1(ind, inp.out_dim)
    feature = @inbounds inp.x[feature_index]
    return feature, output_index
end

Base.size(inp::IsotopicMOInputsUnion) = (inp.out_dim * length(inp.x),)

function Base.vcat(x::MOInputIsotopicByFeatures, y::MOInputIsotopicByFeatures)
    x.out_dim == y.out_dim || throw(DimensionMismatch("out_dim mismatch"))
    return MOInputIsotopicByFeatures(vcat(x.x, y.x), x.out_dim)
end

function Base.vcat(x::MOInputIsotopicByOutputs, y::MOInputIsotopicByOutputs)
    x.out_dim == y.out_dim || throw(DimensionMismatch("out_dim mismatch"))
    return MOInputIsotopicByOutputs(vcat(x.x, y.x), x.out_dim)
end

"""
    MOInput(x::AbstractVector, out_dim::Integer)

A data type to accommodate modelling multi-dimensional output data.
`MOInput(x, out_dim)` has length `length(x) * out_dim`.

```jldoctest
julia> x = [1, 2, 3];

julia> MOInput(x, 2)
6-element KernelFunctions.MOInputIsotopicByOutputs{Int64, Vector{Int64}}:
 (1, 1)
 (2, 1)
 (3, 1)
 (1, 2)
 (2, 2)
 (3, 2)
```
As shown above, an `MOInput` represents a vector of tuples.
The first `length(x)` elements represent the inputs for the first output, the second
`length(x)` elements represent the inputs for the second output, etc.
See [Inputs for Multiple Outputs](@ref) in the docs for more info.

`MOInput` will be deprecated in version 0.11 in favour of `MOInputIsotopicByOutputs`,
and removed in version 0.12.
"""
const MOInput = MOInputIsotopicByOutputs

"""
    prepare_isotopic_multi_output_data(x::AbstractVector, y::ColVecs)

Utility functionality to convert a collection of `N = length(x)` inputs `x`, and a
vector-of-vectors `y` (efficiently represented by a `ColVecs`) into a format suitable for
use with multi-output kernels.

`y[n]` is the vector-valued output corresponding to the input `x[n]`.
Consequently, it is necessary that `length(x) == length(y)`.

For example, if outputs are initially stored in a `num_outputs × N` matrix:
```julia
julia> x = [1.0, 2.0, 3.0];

julia> Y = [1.1 2.1 3.1; 1.2 2.2 3.2]
2×3 Matrix{Float64}:
 1.1  2.1  3.1
 1.2  2.2  3.2

julia> inputs, outputs = prepare_isotopic_multi_output_data(x, ColVecs(Y));

julia> inputs
6-element KernelFunctions.MOInputIsotopicByFeatures{Float64, Vector{Float64}}:
 (1.0, 1)
 (1.0, 2)
 (2.0, 1)
 (2.0, 2)
 (3.0, 1)
 (3.0, 2)

julia> outputs
6-element Vector{Float64}:
 1.1
 1.2
 2.1
 2.2
 3.1
 3.2
```

See also [`prepare_heterotopic_multi_output_data`](@ref).
"""
function prepare_isotopic_multi_output_data(x::AbstractVector, y::ColVecs)
    length(x) == length(y) || throw(ArgumentError("length(x) not equal to length(y)."))
    return MOInputIsotopicByFeatures(x, size(y.X, 1)), vec(y.X)
end

"""
    prepare_isotopic_multi_output_data(x::AbstractVector, y::RowVecs)

Utility functionality to convert a collection of `N = length(x)` inputs `x` and output
vectors `y` (efficiently represented by a `RowVecs`) into a format suitable for
use with multi-output kernels.

`y[n]` is the vector-valued output corresponding to the input `x[n]`.
Consequently, it is necessary that `length(x) == length(y)`.

For example, if outputs are initial stored in an `N × num_outputs` matrix:
```jldoctest
julia> x = [1.0, 2.0, 3.0];

julia> Y = [1.1 1.2; 2.1 2.2; 3.1 3.2]
3×2 Matrix{Float64}:
 1.1  1.2
 2.1  2.2
 3.1  3.2

julia> inputs, outputs = prepare_isotopic_multi_output_data(x, RowVecs(Y));

julia> inputs
6-element KernelFunctions.MOInputIsotopicByOutputs{Float64, Vector{Float64}}:
 (1.0, 1)
 (2.0, 1)
 (3.0, 1)
 (1.0, 2)
 (2.0, 2)
 (3.0, 2)

julia> outputs
6-element Vector{Float64}:
 1.1
 2.1
 3.1
 1.2
 2.2
 3.2
```

See also [`prepare_heterotopic_multi_output_data`](@ref).
"""
function prepare_isotopic_multi_output_data(x::AbstractVector, y::RowVecs)
    length(x) == length(y) || throw(ArgumentError("length(x) not equal to length(y)."))
    return MOInputIsotopicByOutputs(x, size(y.X, 2)), vec(y.X)
end

"""
    prepare_heterotopic_multi_output_data(
        x::AbstractVector, y::AbstractVector{<:Real}, output_indices::AbstractVector{Int},
    )

Utility functionality to convert a collection of inputs `x`, observations `y`, and
`output_indices` into a format suitable for use with multi-output kernels.
Handles the situation in which only one (or a subset) of outputs are observed at each
feature.
Ensures that all arguments are compatible with one another, and returns a vector of inputs
and a vector of outputs.

`y[n]` should be the observed value associated with output `output_indices[n]` at feature
`x[n]`.

```jldoctest
julia> x = [1.0, 2.0, 3.0];

julia> y = [-1.0, 0.0, 1.0];

julia> output_indices = [3, 2, 1];

julia> inputs, outputs = prepare_heterotopic_multi_output_data(x, y, output_indices);

julia> inputs
3-element Vector{Tuple{Float64, Int64}}:
 (1.0, 3)
 (2.0, 2)
 (3.0, 1)

julia> outputs
3-element Vector{Float64}:
 -1.0
  0.0
  1.0
```

See also [`prepare_isotopic_multi_output_data`](@ref).
"""
function prepare_heterotopic_multi_output_data(
    x::AbstractVector, y::AbstractVector{<:Real}, output_indices::AbstractVector{Int}
)
    # Ensure validity of arguments.
    if length(x) != length(y)
        throw(ArgumentError("length(x) != length(y)"))
    end
    if length(x) != length(output_indices)
        throw(ArgumentError("length(x) != length(output_indices)"))
    end

    # Construct inputs and outputs for multi-output kernel.
    x_mogp = map(tuple, x, output_indices)
    return x_mogp, y
end


text:
module TestObjectiveSlack

using Test

using VecMathOptInterface
const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_SlackBridge_ObjectiveSense_error()
    inner = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = MOI.Bridges.Objective.Slack{Float64}(inner)
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.1, x)], 1.2)
    @test_throws(
        ErrorException(
            "Set `MOI.ObjectiveSense` before `MOI.ObjectiveFunction` when" *
            " using `MOI.Bridges.Objective.SlackBridge`.",
        ),
        MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f),
    )
    return
end

function test_SlackBridge_ObjectiveSense_modify()
    inner = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = MOI.Bridges.Objective.Slack{Float64}(inner)
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.1, x)], -1.2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    @test_throws(
        ArgumentError,
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    @test MOI.get(model, MOI.ObjectiveFunction{typeof(f)}()) ≈ f
    return
end

function test_SlackBridge_ObjectiveFunction_modify()
    inner = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = MOI.Bridges.Objective.Slack{Float64}(inner)
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.1, x)], -1.2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    @test MOI.get(model, MOI.ObjectiveFunction{typeof(f)}()) ≈ f
    g = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-1.1, x)], 1.2)
    MOI.set(model, MOI.ObjectiveFunction{typeof(g)}(), g)
    @test MOI.get(model, MOI.ObjectiveFunction{typeof(g)}()) ≈ g
    return
end

function test_SlackBridge_get_ObjectiveFunction_MIN()
    inner = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = MOI.Bridges.Objective.Slack{Float64}(inner)
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.1, x)], -1.2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    @test MOI.get(model, MOI.ObjectiveFunction{typeof(f)}()) ≈ f
    return
end

function test_SlackBridge_get_ObjectiveFunction_MAX()
    inner = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = MOI.Bridges.Objective.Slack{Float64}(inner)
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.1, x)], -1.2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    @test MOI.get(model, MOI.ObjectiveFunction{typeof(f)}()) ≈ f
    return
end

function test_SlackBridge_get_ObjectiveSense()
    inner = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = MOI.Bridges.Objective.Slack{Float64}(inner)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    return
end

function test_SlackBridge_NumberOfVariables()
    inner = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = MOI.Bridges.Objective.Slack{Float64}(inner)
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.1, x)], 1.2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1
    @test MOI.get(inner, MOI.NumberOfVariables()) == 2
    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1
    @test MOI.get(inner, MOI.NumberOfVariables()) == 1
    return
end

function test_SlackBridge_ListOfModelAttributesSet()
    inner = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = MOI.Bridges.Objective.Slack{Float64}(inner)
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.1, x)], 1.2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    attr = MOI.get(model, MOI.ListOfModelAttributesSet())
    @test length(attr) == 2
    @test MOI.ObjectiveSense() in attr
    @test MOI.ObjectiveFunction{typeof(f)}() in attr
    attr = MOI.get(inner, MOI.ListOfModelAttributesSet())
    @test length(attr) == 2
    @test MOI.ObjectiveSense() in attr
    @test MOI.ObjectiveFunction{MOI.VariableIndex}() in attr
    return
end

function test_SlackBridge_ListOfVariableIndices()
    inner = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = MOI.Bridges.Objective.Slack{Float64}(inner)
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.1, x)], 1.2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    @test MOI.get(model, MOI.ListOfVariableIndices()) == [x]
    @test length(MOI.get(inner, MOI.ListOfVariableIndices())) == 2
    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    @test MOI.get(model, MOI.ListOfVariableIndices()) == [x]
    @test MOI.get(inner, MOI.ListOfVariableIndices()) == [x]
    return
end

function test_SlackBridge_ListOfConstraintTypesPresent()
    inner = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = MOI.Bridges.Objective.Slack{Float64}(inner)
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.1, x)], 1.2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    @test MOI.get(model, MOI.ListOfConstraintTypesPresent()) == []
    @test MOI.get(inner, MOI.ListOfConstraintTypesPresent()) ==
          [(MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64})]
    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    @test MOI.get(model, MOI.ListOfConstraintTypesPresent()) == []
    @test MOI.get(inner, MOI.ListOfConstraintTypesPresent()) == []
    return
end

function test_SlackBridge_ListOfConstraintIndices()
    inner = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = MOI.Bridges.Objective.Slack{Float64}(inner)
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.1, x)], 1.2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    attr = MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{Float64},
        MOI.GreaterThan{Float64},
    }()
    @test MOI.get(model, attr) == []
    @test length(MOI.get(inner, attr)) == 1
    attr = MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{Float64},
        MOI.LessThan{Float64},
    }()
    @test MOI.get(model, attr) == []
    @test length(MOI.get(inner, attr)) == 0
    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    @test MOI.get(model, attr) == []
    @test length(MOI.get(inner, attr)) == 0
    return
end

function test_SlackBridge_ObjectiveFunctionValue()
    inner = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = MOI.Bridges.Objective.Slack{Float64}(inner)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.GreaterThan(2.0))
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.1, x)], -1.2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.Utilities.set_mock_optimize!(
        inner,
        mock -> MOI.Utilities.mock_optimize!(mock, [2.0, 1.0]),
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1.0
    return
end

function test_original()
    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    bridged_mock = MOI.Bridges.Objective.Slack{Float64}(mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0, 2.0, 5.0]),
        ),
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0, 2.0, 7.0]),
        ),
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0, 2.0, 2.0]),
        ),
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [1.0, 2.0, 7.0]),
        ),
    )
    MOI.Test.test_objective_qp_ObjectiveFunction_edge_cases(
        bridged_mock,
        MOI.Test.Config(),
    )

    @test MOI.Bridges.is_objective_bridged(bridged_mock)
    @test MOI.get(bridged_mock, MOI.ObjectiveFunctionType()) ==
          MOI.ScalarQuadraticFunction{Float64}
    @test MOI.get(bridged_mock, MOI.ListOfModelAttributesSet()) == [
        MOI.ObjectiveSense(),
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
    ]

    var_names = ["x", "y"]
    xy = MOI.get(bridged_mock, MOI.ListOfVariableIndices())
    MOI.set(bridged_mock, MOI.VariableName(), xy, var_names)
    abs = MOI.get(mock, MOI.ListOfVariableIndices())
    @test length(abs) == 3
    MOI.set(mock, MOI.VariableName(), abs[3], "s")
    cquad = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarQuadraticFunction{Float64},
            MOI.LessThan{Float64},
        }(),
    )
    MOI.set(bridged_mock, MOI.ConstraintName(), cquad[1], "quad")

    s = """
    variables: x, y, s
    x >= 1.0
    y >= 2.0
    quad: 1.0 * x * x + 1.0 * x * y + 1.0 * y * y + -1.0 * s <= 0.0
    minobjective: s
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        [var_names; "s"],
        ["quad"],
        [
            ("x", MOI.GreaterThan{Float64}(1.0)),
            ("y", MOI.GreaterThan{Float64}(2.0)),
        ],
    )
    bridged_var_names = ["x", "y"]
    s = """
    variables: x, y
    x >= 1.0
    y >= 2.0
    minobjective: 1.0 * x * x + 1.0 * x * y + 1.0 * y * y
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        var_names,
        String[],
        [
            ("x", MOI.GreaterThan{Float64}(1.0)),
            ("y", MOI.GreaterThan{Float64}(2.0)),
        ],
    )
    err = ArgumentError(
        "Objective bridge of type `$(VecMathOptInterface.Bridges.Objective.SlackBridge{Float64,VecMathOptInterface.ScalarQuadraticFunction{Float64},VecMathOptInterface.ScalarQuadraticFunction{Float64}})`" *
        " does not support modifying the objective sense. As a workaround, set" *
        " the sense to `MOI.FEASIBILITY_SENSE` to clear the objective function" *
        " and bridges.",
    )
    @test_throws err MOI.set(bridged_mock, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    obj = MOI.get(
        bridged_mock,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
    )
    MOI.set(bridged_mock, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    @test !MOI.Bridges.is_objective_bridged(bridged_mock)
    @test MOI.get(bridged_mock, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE
    @test MOI.get(bridged_mock, MOI.ListOfModelAttributesSet()) ==
          [MOI.ObjectiveSense()]
    MOI.set(bridged_mock, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(
        bridged_mock,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
        obj,
    )
    abs = MOI.get(mock, MOI.ListOfVariableIndices())
    @test length(abs) == 3
    MOI.set(mock, MOI.VariableName(), abs[3], "s")
    cquad = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarQuadraticFunction{Float64},
            MOI.GreaterThan{Float64},
        }(),
    )
    MOI.set(bridged_mock, MOI.ConstraintName(), cquad[1], "quad")

    s = """
    variables: x, y, s
    x >= 1.0
    y >= 2.0
    quad: 1.0 * x * x + 1.0 * x * y + 1.0 * y * y + -1.0 * s >= 0.0
    maxobjective: s
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        [var_names; "s"],
        ["quad"],
        [
            ("x", MOI.GreaterThan{Float64}(1.0)),
            ("y", MOI.GreaterThan{Float64}(2.0)),
        ],
    )
    s = """
    variables: x, y
    x >= 1.0
    y >= 2.0
    maxobjective: 1.0 * x * x + 1.0 * x * y + 1.0 * y * y
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        var_names,
        String[],
        [
            ("x", MOI.GreaterThan{Float64}(1.0)),
            ("y", MOI.GreaterThan{Float64}(2.0)),
        ],
    )
    _test_delete_objective(
        bridged_mock,
        2,
        (
            (MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}, 0),
            (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}, 0),
        ),
    )
    return
end

function test_quadratic_integration()
    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    bridged_mock = MOI.Bridges.Objective.Slack{Float64}(mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [4 / 7, 3 / 7, 6 / 7, 13 / 7],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [5 / 7, 6 / 7],
        ),
    )
    MOI.Test.test_quadratic_integration(bridged_mock, MOI.Test.Config())
    return
end

function test_quadratic_duplicate_terms()
    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    bridged_mock = MOI.Bridges.Objective.Slack{Float64}(mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [4 / 7, 3 / 7, 6 / 7, 13 / 7],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [5 / 7, 6 / 7],
        ),
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [4 / 7, 3 / 7, 6 / 7, -2 * 13 / 7],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [10 / 7, 12 / 7],
        ),
    )
    MOI.Test.test_quadratic_duplicate_terms(bridged_mock, MOI.Test.Config())
    return
end

function test_quadratic_nonhomogeneous()
    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    bridged_mock = MOI.Bridges.Objective.Slack{Float64}(mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1 / 4, 3 / 4, 2.875],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                [11 / 4],
            (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1.0],
        ),
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1.0, 0.0, 3.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                [-2.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [1.0],
        ),
    )
    MOI.Test.test_quadratic_nonhomogeneous(bridged_mock, MOI.Test.Config())
    return
end

end  # module

TestObjectiveSlack.runtests()


text:
is_l1(p::Penalty) = isa(p, L1R)
function is_elnet(cp::CompositePenalty)
    length(cp.penalties) == 2 || return false
    isa(cp.penalties[1], L2R) && isa(cp.penalties[2], L1R) && return true
    isa(cp.penalties[1], L1R) && isa(cp.penalties[2], L2R) && return true
    return false
end

get_l1(p::ScaledPenalty{L1Penalty}) = p
get_l1(cp::CompositePenalty) = cp.penalties[findfirst(e->isa(e, L1R), cp.penalties)]
get_l2(p::ScaledPenalty{L1Penalty}) = NoPenalty()
get_l2(cp::CompositePenalty) = cp.penalties[findfirst(e->isa(e, L2R), cp.penalties)]

getscale_l1(p1::Union{NoPenalty,L1R}) = p1 |> getscale
getscale_l1(cp::CompositePenalty) = is_elnet(cp) ? cp |> get_l1 |> getscale :
                                                   @error "Case not implemented."
getscale_l2(p1::Union{NoPenalty,L1R}) = 0.0
getscale_l2(p2::L2R) = p2 |> getscale
getscale_l2(cp::CompositePenalty) = is_elnet(cp) ? cp |> get_l2 |> getscale :
                                                   @error "Case not implemented."


text:
module MainModule
using Gtk
using Gtk.ShortNames
using eve
using DataFrames

include("GtkUtils/GObjectImpls.jl")
include("GtkUtils/GtkTreeModelIface.jl")
include("MarginDataView.jl")
include("DataFrameStore.jl")
include("DataFrameView.jl")
include("JuliaView.jl")
include("MainWindow.jl")
end

max_stars_repo_name: barcharcraz/eve-julia-gui

max_issues_repo_name: barcharcraz/eve-julia-gui

max_forks_repo_name: barcharcraz/eve-julia-gui


text:
# Autogenerated wrapper script for PCRE2_jll for aarch64-apple-darwin
export libpcre2_16, libpcre2_32, libpcre2_8

JLLWrappers.@generate_wrapper_header("PCRE2")
JLLWrappers.@declare_library_product(libpcre2_16, "@rpath/libpcre2-16.0.dylib")
JLLWrappers.@declare_library_product(libpcre2_32, "@rpath/libpcre2-32.0.dylib")
JLLWrappers.@declare_library_product(libpcre2_8, "@rpath/libpcre2-8.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libpcre2_16,
        "lib/libpcre2-16.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcre2_32,
        "lib/libpcre2-32.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcre2_8,
        "lib/libpcre2-8.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Plots
using CPUTime

function makeplot()
    anim = @animate for i in 1:200
        scatter(
            t->sinpi(t+i/100),
            t->cospi(t+i/100),
            [0],
            size=(400,400),
            xlims = (-1.5, 1.5),
            ylims = (-1.5, 1.5),
            aspect_ratio =1,
        )
    end
    gif(anim, "example_3_by_jl.gif", fps = 6)
end

@time makeplot()


text:
#
# NaiveLocalFilters.jl -
#
# Naive (non-optimized) implementation of local filters.  This module provides
# reference implementations of methods for checking results and benchmarking.
#
#------------------------------------------------------------------------------
#
# This file is part of the `LocalFilters.jl` package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2017-2018, Éric Thiébaut.
#

module NaiveLocalFilters

using Compat
using LocalFilters
using LocalFilters:
    Kernel,
    Neighborhood,
    RectangularBox,
    _store!,
    _typeofsum,
    axes,
    cartesianregion,
    coefs,
    limits,
    offset

import LocalFilters:
    bottom_hat,
    cartesianregion,
    closing,
    convolve!,
    convolve,
    dilate!,
    dilate,
    erode!,
    erode,
    localextrema!,
    localextrema,
    localfilter!,
    localmean,
    localmean!,
    opening,
    top_hat

localmean(variant::Val, A::AbstractArray{T,N}, B=3) where {T,N} =
    localmean(variant, A, Neighborhood{N}(B))

localmean(variant::Val, A::AbstractArray{T}, B::RectangularBox) where {T} =
    localmean!(variant, similar(Array{float(T)}, axes(A)), A, B)

localmean(variant::Val, A::AbstractArray{T}, B::Kernel{Bool}) where {T} =
    localmean!(variant, similar(Array{float(T)}, axes(A)), A, B)

localmean(variant::Val, A::AbstractArray{T}, B::Kernel{K}) where {T,K} =
    localmean!(variant, similar(Array{float(promote_type(T,K))},
                                axes(A)), A, B)

convolve(variant::Val, A::AbstractArray{T,N}, B=3) where {T,N} =
    convolve(variant, A, Neighborhood{N}(B))

convolve(variant::Val, A::AbstractArray{T}, B::RectangularBox) where {T} =
    convolve!(variant, similar(Array{_typeofsum(T)}, axes(A)), A, B)

convolve(variant::Val, A::AbstractArray{T}, B::Kernel{Bool}) where {T} =
    convolve!(variant, similar(Array{_typeofsum(T)}, axes(A)), A, B)

convolve(variant::Val, A::AbstractArray{T}, B::Kernel{K}) where {T,K} =
    convolve!(variant, similar(Array{_typeofsum(promote_type(T,K))},
                               axes(A)), A, B)

dilate(variant::Val, A::AbstractArray, args...) =
    dilate!(variant, similar(A), A, args...)

erode(variant::Val, A::AbstractArray, args...) =
    erode!(variant, similar(A), A, args...)

closing(variant::Val, A::AbstractArray{T,N}, B=3) where {T,N} =
    closing(variant, A, Neighborhood{N}(B))

closing(variant::Val, A::AbstractArray, B::Neighborhood) =
    erode(variant, dilate(variant, A, B), B)

opening(variant::Val, A::AbstractArray{T,N}, B=3) where {T,N} =
    opening(variant, A, Neighborhood{N}(B))

opening(variant::Val, A::AbstractArray, B::Neighborhood) =
    dilate(variant, erode(variant, A, B), B)

top_hat(variant::Val, a, r=3) = a .- opening(variant, a, r)
top_hat(variant::Val, a, r, s) =
    top_hat(variant, closing(variant, a, s), r)

bottom_hat(variant::Val, a, r=3) = closing(variant, a, r) .- a
bottom_hat(variant::Val, a, r, s) =
    bottom_hat(variant, opening(variant, a, s), r)

localextrema(variant::Val, A::AbstractArray, args...) =
    localextrema!(variant, similar(A), similar(A), A, args...)

function erode!(::Val{:Naive},
                dst::AbstractArray{T,N},
                A::AbstractArray{T,N},
                B::RectangularBox{N}) where {T,N}
    @assert axes(dst) == axes(A)
    localfilter!(dst, A, B,
                 (a)     -> typemax(T),
                 (v,a,b) -> min(v, a),
                 (d,i,v) -> d[i] = v)
end

function dilate!(::Val{:Naive},
                 dst::AbstractArray{T,N},
                 A::AbstractArray{T,N},
                 B::RectangularBox{N}) where {T,N}
    @assert axes(dst) == axes(A)
    localfilter!(dst, A, B,
                 (a)     -> typemin(T),
                 (v,a,b) -> max(v, a),
                 (d,i,v) -> d[i] = v)
end

#------------------------------------------------------------------------------
# Variants for computing interesecting regions.  This is the most critical
# part for indexing a neighborhood (apart from using a non-naive algorithm).

# "Base" variant: use constructors and methods provided by the Base package.
@inline function _cartesianregion(::Val{:Base},
                                  imin::CartesianIndex{N},
                                  imax::CartesianIndex{N},
                                  i::CartesianIndex{N},
                                  off::CartesianIndex{N}) where N
    cartesianregion(max(imin, i - off), min(imax, i + off))
end
@inline function _cartesianregion(::Val{:Base},
                                  imin::CartesianIndex{N},
                                  imax::CartesianIndex{N},
                                  i::CartesianIndex{N},
                                  kmin::CartesianIndex{N},
                                  kmax::CartesianIndex{N}) where N
    cartesianregion(max(imin, i - kmax), min(imax, i - kmin))
end

# "NTuple" variant: use `ntuple()` to expand expressions.
@inline function _cartesianregion(::Val{:NTuple},
                                  imin::CartesianIndex{N},
                                  imax::CartesianIndex{N},
                                  i::CartesianIndex{N},
                                  off::CartesianIndex{N}) where N
    cartesianregion(ntuple(k -> (max(imin[k], i[k] - off[k]) :
                                 min(imax[k], i[k] + off[k])), N))
end
@inline function _cartesianregion(::Val{:NTuple},
                                  imin::CartesianIndex{N},
                                  imax::CartesianIndex{N},
                                  i::CartesianIndex{N},
                                  kmin::CartesianIndex{N},
                                  kmax::CartesianIndex{N}) where N
    cartesianregion(ntuple(k -> (max(imin[k], i[k] - kmax[k]) :
                                 min(imax[k], i[k] - kmin[k])), N))
end

# "NTupleVal" variant: use `ntuple()` to expand expressions and `Val(...)` to
# force specialized versions.
for N in (1,2,3,4)
    @eval begin
        @inline function _cartesianregion(::Val{:NTupleVar},
                                          imin::CartesianIndex{$N},
                                          imax::CartesianIndex{$N},
                                          i::CartesianIndex{$N},
                                          off::CartesianIndex{$N})
            cartesianregion(ntuple(k -> (max(imin[k], i[k] - off[k]) :
                                         min(imax[k], i[k] + off[k])),
                                   Val($N)))
        end
        @inline function _cartesianregion(::Val{:NTupleVar},
                                          imin::CartesianIndex{$N},
                                          imax::CartesianIndex{$N},
                                          i::CartesianIndex{$N},
                                          kmin::CartesianIndex{$N},
                                          kmax::CartesianIndex{$N})
            cartesianregion(ntuple(k -> (max(imin[k], i[k] - kmax[k]) :
                                         min(imax[k], i[k] - kmin[k])),
                                   Val($N)))
        end
    end
end

# "Map" variant: Use `map()` to expand expressions.
@inline _range(imin::Int, imax::Int, i::Int, off::Int) =
    max(imin, i - off) : min(imax, i + off)
@inline _range(imin::Int, imax::Int, i::Int, kmin::Int, kmax::Int) =
    max(imin, i - kmax) : min(imax, i - kmin)
@inline function _cartesianregion(::Val{:Map},
                                  imin::CartesianIndex{N},
                                  imax::CartesianIndex{N},
                                  i::CartesianIndex{N},
                                  off::CartesianIndex{N}) where N
    cartesianregion(map(_range, imin.I, imax.I, i.I, off.I))
end
@inline function _cartesianregion(::Val{:Map},
                                  imin::CartesianIndex{N},
                                  imax::CartesianIndex{N},
                                  i::CartesianIndex{N},
                                  kmin::CartesianIndex{N},
                                  kmax::CartesianIndex{N}) where N
    cartesianregion(map(_range, imin.I, imax.I, i.I, kmin.I, kmax.I))
end

#------------------------------------------------------------------------------
# Methods for rectangular boxes (with optimization for centered boxes).

function localmean!(variant::Val,
                    dst::AbstractArray{Td,N},
                    A::AbstractArray{Ts,N},
                    B::RectangularBox{N}) where {Td,Ts,N}
    @assert axes(dst) == axes(A)
    R = cartesianregion(A)
    imin, imax = limits(R)
    kmin, kmax = limits(B)
    T = _typeofsum(Ts)
    if kmin === -kmax
        off = kmax
        @inbounds for i in R
            n, s = 0, zero(T)
            for j in _cartesianregion(variant, imin, imax, i, off)
                s += A[j]
                n += 1
            end
            _store!(dst, i, s/n)
        end
    else
        @inbounds for i in R
            n, s = 0, zero(T)
            for j in _cartesianregion(variant, imin, imax, i, kmin, kmax)
                s += A[j]
                n += 1
            end
            _store!(dst, i, s/n)
        end
    end
    return dst
end

function erode!(variant::Val,
                Amin::AbstractArray{T,N},
                A::AbstractArray{T,N},
                B::RectangularBox{N}) where {T,N}
    @assert axes(Amin) == axes(A)
    R = cartesianregion(A)
    imin, imax = limits(R)
    kmin, kmax = limits(B)
    tmax = typemax(T)
    if kmin === -kmax
        off = kmax
        @inbounds for i in R
            vmin = tmax
            for j in _cartesianregion(variant, imin, imax, i, off)
                vmin = min(vmin, A[j])
            end
            Amin[i] = vmin
        end
    else
        @inbounds for i in R
            vmin = tmax
            for j in _cartesianregion(variant, imin, imax, i, kmin, kmax)
                vmin = min(vmin, A[j])
            end
            Amin[i] = vmin
        end
    end
    return Amin
end

function dilate!(variant::Val,
                 Amax::AbstractArray{T,N},
                 A::AbstractArray{T,N},
                 B::RectangularBox{N}) where {T,N}
    @assert axes(Amax) == axes(A)
    R = cartesianregion(A)
    imin, imax = limits(R)
    kmin, kmax = limits(B)
    tmin = typemin(T)
    if kmin === -kmax
        off = kmax
        @inbounds for i in R
            vmax = tmin
            for j in _cartesianregion(variant, imin, imax, i, off)
                vmax = max(vmax, A[j])
            end
            Amax[i] = vmax
        end
    else
        @inbounds for i in R
            vmax = tmin
            for j in _cartesianregion(variant, imin, imax, i, kmin, kmax)
                vmax = max(vmax, A[j])
            end
            Amax[i] = vmax
        end
    end
    return Amax
end

function localextrema!(variant::Val,
                       Amin::AbstractArray{T,N},
                       Amax::AbstractArray{T,N},
                       A::AbstractArray{T,N},
                       B::RectangularBox{N}) where {T,N}
    @assert axes(Amin) == axes(Amax) == axes(A)
    R = cartesianregion(A)
    imin, imax = limits(R)
    kmin, kmax = limits(B)
    tmin, tmax = limits(T)
    if kmin === -kmax
        off = kmax
        @inbounds for i in R
            vmin, vmax = tmax, tmin
            for j in _cartesianregion(variant, imin, imax, i, off)
                vmin = min(vmin, A[j])
                vmax = max(vmax, A[j])
            end
            Amin[i] = vmin
            Amax[i] = vmax
        end
    else
        @inbounds for i in R
            vmin, vmax = tmax, tmin
            for j in _cartesianregion(variant, imin, imax, i, kmin, kmax)
                vmin = min(vmin, A[j])
                vmax = max(vmax, A[j])
            end
            Amin[i] = vmin
            Amax[i] = vmax
        end
    end
    return Amin, Amax
end

#------------------------------------------------------------------------------
# Nethods for kernels of booleans.

function localmean!(variant::Val,
                    dst::AbstractArray{Td,N},
                    A::AbstractArray{Ts,N},
                    B::Kernel{Bool,N}) where {Td,Ts,N}
    @assert axes(dst) == axes(A)
    R = cartesianregion(A)
    imin, imax = limits(R)
    kmin, kmax = limits(B)
    ker, off = coefs(B), offset(B)
    T = _typeofsum(Ts)
    @inbounds for i in R
        n, s = 0, zero(T)
        k = i + off
        for j in _cartesianregion(variant, imin, imax, i, kmin, kmax)
            if ker[k-j]
                n += 1
                s += A[j]
            end
        end
        _store!(dst, i, s/n)
    end
    return dst
end

function erode!(variant::Val,
                Amin::AbstractArray{T,N},
                A::AbstractArray{T,N},
                B::Kernel{Bool,N}) where {T,N}
    @assert axes(Amin) == axes(A)
    R = cartesianregion(A)
    imin, imax = limits(R)
    kmin, kmax = limits(B)
    ker, off = coefs(B), offset(B)
    tmax = typemax(T)
    @inbounds for i in R
        vmin = tmax
        k = i + off
        for j in _cartesianregion(variant, imin, imax, i, kmin, kmax)
            #if ker[k-j] && A[j] < vmin
            #    vmin = A[j]
            #end
            vmin = ker[k-j] && A[j] < vmin ? A[j] : vmin
        end
        Amin[i] = vmin
    end
    return Amin
end

function dilate!(variant::Val,
                 Amax::AbstractArray{T,N},
                 A::AbstractArray{T,N},
                 B::Kernel{Bool,N}) where {T,N}
    @assert axes(Amax) == axes(A)
    R = cartesianregion(A)
    imin, imax = limits(R)
    kmin, kmax = limits(B)
    ker, off = coefs(B), offset(B)
    tmin = typemin(T)
    @inbounds for i in R
        vmax = tmin
        k = i + off
        for j in _cartesianregion(variant, imin, imax, i, kmin, kmax)
            #if ker[k-j] && A[j] > vmax
            #    vmax = A[j]
            #end
            vmax = ker[k-j] && A[j] > vmax ? A[j] : vmax
        end
        Amax[i] = vmax
    end
    return Amax
end

function localextrema!(variant::Val,
                       Amin::AbstractArray{T,N},
                       Amax::AbstractArray{T,N},
                       A::AbstractArray{T,N},
                       B::Kernel{Bool,N}) where {T,N}
    @assert axes(Amin) == axes(Amax) == axes(A)
    R = cartesianregion(A)
    imin, imax = limits(R)
    kmin, kmax = limits(B)
    tmin, tmax = limits(T)
    ker, off = coefs(B), offset(B)
    @inbounds for i in R
        vmin, vmax = tmax, tmin
        k = i + off
        for j in _cartesianregion(variant, imin, imax, i, kmin, kmax)
            #if ker[k-j]
            #    vmin = min(vmin, A[j])
            #    vmax = max(vmax, A[j])
            #end
            vmin = ker[k-j] && A[j] < vmin ? A[j] : vmin
            vmax = ker[k-j] && A[j] > vmax ? A[j] : vmax
        end
        Amin[i] = vmin
        Amax[i] = vmax
    end
    return Amin, Amax
end

#------------------------------------------------------------------------------
# Methods for other kernels.

function localmean!(variant::Val,
                    dst::AbstractArray{Td,N},
                    A::AbstractArray{Ts,N},
                    B::Kernel{Tk,N}) where {Td,Ts,Tk,N}
    @assert axes(dst) == axes(A)
    R = cartesianregion(A)
    imin, imax = limits(R)
    kmin, kmax = limits(B)
    ker, off = coefs(B), offset(B)
    T = _typeofsum(promote_type(Ts, Tk))
    @inbounds for i in R
        s1, s2 = zero(T), zero(T)
        k = i + off
        for j in _cartesianregion(variant, imin, imax, i, kmin, kmax)
            w = ker[k-j]
            s1 += w*A[j]
            s2 += w
        end
        _store!(dst, i, s1/s2)
    end
    return dst
end

function erode!(variant::Val,
                Amin::AbstractArray{T,N},
                A::AbstractArray{T,N},
                B::Kernel{T,N}) where {T<:AbstractFloat,N}
    @assert axes(Amin) == axes(A)
    R = cartesianregion(A)
    imin, imax = limits(R)
    kmin, kmax = limits(B)
    ker, off = coefs(B), offset(B)
    tmax = typemax(T)
    @inbounds for i in R
        vmin = tmax
        k = i + off
        for j in _cartesianregion(variant, imin, imax, i, kmin, kmax)
            vmin = min(vmin, A[j] - ker[k-j])
        end
        Amin[i] = vmin
    end
    return Amin
end

function dilate!(variant::Val,
                 Amax::AbstractArray{T,N},
                 A::AbstractArray{T,N},
                 B::Kernel{T,N}) where {T<:AbstractFloat,N}
    @assert axes(Amax) == axes(A)
    R = cartesianregion(A)
    imin, imax = limits(R)
    kmin, kmax = limits(B)
    ker, off = coefs(B), offset(B)
    tmin = typemin(T)
    @inbounds for i in R
        vmax = tmin
        k = i + off
        for j in _cartesianregion(variant, imin, imax, i, kmin, kmax)
            vmax = max(vmax, A[j] + ker[k-j])
        end
        Amax[i] = vmax
    end
    return Amax
end

function convolve!(variant::Val,
                   dst::AbstractArray{Td,N},
                   A::AbstractArray{Ts,N},
                   B::Kernel{Tk,N}) where {Td,Ts,Tk,N}
    @assert axes(dst) == axes(A)
    R = cartesianregion(A)
    imin, imax = limits(R)
    kmin, kmax = limits(B)
    ker, off = coefs(B), offset(B)
    T = _typeofsum(promote_type(Ts, Tk))
    @inbounds for i in R
        v = zero(T)
        k = i + off
        for j in _cartesianregion(variant, imin, imax, i, kmin, kmax)
            v += A[j]*ker[k-j]
        end
        _store!(dst, i, v)
    end
    return dst
end

function localextrema!(variant::Val,
                       Amin::AbstractArray{T,N},
                       Amax::AbstractArray{T,N},
                       A::AbstractArray{T,N},
                       B::Kernel{T,N}) where {T<:AbstractFloat,N}
    @assert axes(Amin) == axes(Amax) == axes(A)
    R = cartesianregion(A)
    imin, imax = limits(R)
    kmin, kmax = limits(B)
    tmin, tmax = limits(T)
    ker, off = coefs(B), offset(B)
    @inbounds for i in R
        vmin, vmax = tmax, tmin
        k = i + off
        for j in _cartesianregion(variant, imin, imax, i, kmin, kmax)
            vmin = min(vmin, A[j] - ker[k-j])
            vmax = max(vmax, A[j] + ker[k-j])
        end
        Amin[i] = vmin
        Amax[i] = vmax
    end
    return Amin, Amax
end

end # module


text:
using MultiScaleArrays, OrdinaryDiffEq, DiffEqBase, StochasticDiffEq, SafeTestsets
using Test

@time @testset "Tuple Nodes" begin include("tuple_nodes.jl") end
@time @testset "Bisect Search Tests" begin include("bisect_search_tests.jl") end
@time @testset "Indexing and Creation Tests" begin include("indexing_and_creation_tests.jl") end
@time @testset "Values Indexing" begin include("values_indexing.jl") end
@time @testset "Get Indices Tests" begin include("get_indices.jl") end
@time @testset "Additional Fields Test" begin include("additional_fields_test.jl") end
@time @testset "Dynamic DiffEq Tests" begin include("dynamic_diffeq.jl") end
@time @testset "Single Layer DiffEq Tests" begin include("single_layer_diffeq.jl") end
@time @safetestset "New Nodes Tests" begin include("new_nodes.jl") end


text:
using Test
using ShuffleProofs: tobig, int2bytes, bitlength, interpret

x = big(2)^100

@test tobig(int2bytes(x)) == x

x = 2^63 - 1

@test int2bytes(big(x)) == reinterpret(UInt8, [x])

y = UInt64(2)^63 + UInt64(1)

@test int2bytes(big(y)) == reinterpret(UInt8, [y])

z = UInt128(2)^128 - UInt128(1)

@test int2bytes(big(z)) == reinterpret(UInt8, [z])

u = 300
q = interpret(Vector{UInt8}, u)
#@test Int(frombytes(UInt64, copy(q))) == u

@test interpret(UInt64, copy(q)) == UInt64(u)

#@test frombytes(BigInt, copy(q)) == big(u)
@test interpret(BigInt, copy(q)) == big(u)

# Test cases can be generated on:
# https://compiler.javatpoint.com/opr/test.jsp?filename=BigIntegerBitLengthExample

x = 2323535352 ## Java says 32
@test bitlength(x) == bitlength(big(x)) == 32

x = 121232 # Java says 17
@test bitlength(x) == bitlength(big(x)) == 17

x = 23235323423415352
@test bitlength(x) == bitlength(big(x)) == 55


text:
function cubic_ab(model::ABCubicModel,V,T,z=SA[1.0],n=sum(z))
    invn2 = (one(n)/n)^2
    a = model.params.a.values
    b = model.params.b.values
    α = @f(α_function,model.alpha)
    if length(z)>1
        ā,b̄ = @f(mixing_rule,model.mixing,α,a,b)
    else
        ā = a[1,1]*α[1]
        b̄ = b[1,1]
    end
    return ā ,b̄
end

function second_virial_coefficient(model::ABCubicModel,T,z = SA[1.0])
    #@info "fast shortcut"
    a,b = cubic_ab(model,1e-4,T,z)
    return b-a/(R̄*T)
end

# function x0_volume_sc(model::ABCubicModel,p,T,z)
#     Zc = cubic_zc(model)
#     return Zc*R̄*T/p
# end

function volume(model::ABCubicModel,p,T,z=SA[1.0];phase=:unknown,threaded=false)
    lb_v   =lb_volume(model,z)
    xx = z/sum(z)
    RTp = R̄*T/p
    _poly = cubic_poly(model,p,T,z)

    #sols = Solvers.poly3(_poly)
    sols = Solvers.roots3(_poly)
    #xa = real.(sols) .* RTp
    #@show xa
    #with this, real sol is first
    function imagfilter(x)
        absx = abs(imag(x))
        return absx <8*eps(typeof(absx))
    end

    isreal = map(imagfilter,sols)
    _sols = sort(real.(sols)) .* RTp

    #@show _sols[isreal]
    #@show model.params.b.values
    #@show sum(isreal)
    if sum(isreal) == 3 #3 roots
        sols = sort(real.(sols))
        vg = last(sols)*RTp
        _vl = first(sols)*RTp
        vl = ifelse(_vl>lb_v,_vl,vg)
    elseif  sum(isreal) == 1
        i = findfirst(imagfilter,sols)
        vl = real(sols[i])*RTp
        vg = real(sols[i])*RTp
    elseif  sum(isreal) == 0
        @show sols
        error("weird")
    end

    #fp(_v) = pressure(model,_v,T,z) - p
    #@show vl
    #@show vg
    function gibbs(v)
        _df,f =  ∂f(model,v,T,z)
        dv,dt = _df
        if abs((p+dv)/p) > 0.03
            return Inf
        else
            return f  +p*v
        end
    end
    #this catches the supercritical phase as well
    if vl ≈ vg
        return vl
    end
    if phase == :unknown
        gg = gibbs(vg)
        gl = gibbs(vl)
        #@show vg,vl
        #@show gg,gl
        return ifelse(gg<gl,vg,vl)
    elseif is_liquid(phase)
        return vl
    elseif is_vapour(phase)
        return vg
    # else
    #     gg = gibbs(vg)
    #     gl = gibbs(vl)
    #     #@show vg,vl
    #     #@show gg,gl
    #     return ifelse(gg<gl,vg,vl)
    end
end

function ab_consts(model::CubicModel)
    return ab_consts(typeof(model))
end


text:
using Test, PolynomialBases
import SymPy

tol = 5.e-15

for p in 0:6
    basis_sympy = LobattoLegendre(p, SymPy.Sym)
    basis_float = LobattoLegendre(p, Float64)
    @test maximum(abs.( float.(basis_sympy.nodes) - basis_float.nodes )) < tol
    @test maximum(abs.( float.(basis_sympy.weights) - basis_float.weights )) < tol
    @test maximum(abs.( float.(basis_sympy.baryweights) - basis_float.baryweights )) < tol
    @test maximum(abs.( float.(basis_sympy.D) - basis_float.D, )) < 2tol
end
@test_throws ArgumentError LobattoLegendre(7, SymPy.Sym)

for p in 0:4
    basis_sympy = GaussLegendre(p, SymPy.Sym)
    basis_float = GaussLegendre(p, Float64)
    @test maximum(abs.( float.(basis_sympy.nodes) - basis_float.nodes )) < tol
    @test maximum(abs.( float.(basis_sympy.weights) - basis_float.weights )) < tol
    @test maximum(abs.( float.(basis_sympy.baryweights) - basis_float.baryweights )) < tol
    @test maximum(abs.( float.(basis_sympy.D) - basis_float.D, )) < 2tol
end
@test_throws ArgumentError GaussLegendre(5, SymPy.Sym)

interpolation_matrix([-1, 1], LobattoLegendre(4, SymPy.Sym))
interpolation_matrix([-1, 1], GaussLegendre(4, SymPy.Sym))


text:
using PlotlyJS

function two_hists()
    x0 = randn(500)
    x1 = x0+1

    trace1 = histogram(x=x0, opacity=0.75)
    trace2 = histogram(x=x1, opacity=0.75)
    data = [trace1, trace2]
    layout = Layout(barmode="overlay")
    plot(data, layout)
end


text:
#///////////////////////////////////////
#// File Name: runtests.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2020/11/03
#// Description: Test script for RATiLQR
#///////////////////////////////////////

using Test

using Distributed
if nprocs() < 2
    addprocs(1)
end

@everywhere using RATiLQR

@testset "RATiLQR Unit Tests" begin
@info "Executing iLEQG Test"
include("ileqg_test.jl")
@info "Executing Cross Entropy Bilevel Optimization Test"
include("cross_entropy_bilevel_optimization_test.jl")
@info "Executing Nelder-Mead Simplex Bilevel Optimization Test"
include("nelder_mead_bilevel_optimization_test.jl")
@info "Executing PETS Test"
include("pets_test.jl")
end


text:
include("../src/env_config.jl")

module EnvConfigTest
using ..EnvConfig

greet2()

end  # module

using .EnvConfig

greet1()


text:
# --------------------------------------------------
# Activation Functions
# --------------------------------------------------
immutable IdentityAct <: ActivationFunc
end

@inline function forward(f::IdentityAct, x::Float32)
    x
end

@inline function backward(f::IdentityAct, x::Float32)
    x
end

# --------------------------------------------------
function sigmoid(x)
    one(x) / (one(x) + exp(-x))
end

function sigmoid(X::Signal)
    Z = zeros(X)
    for (row, col) in zip(findn(X)...)
        Z[row, col] = sigmoid(X[row, col])
    end
    Z
end

function sigmoid(diff::Symbol, X::Signal)
    Zd = zeros(X)
    for (row, col) in zip(findn(X)...)
        s = sigmoid(X[row, col])
        Zd[row, col] = s * (one(eltype(X)) - s)
    end
    Zd
end

# --------------------------------------------------
function relu(X::Signal)
    bottom = zero(eltype(X))
    Z = zeros(X)
    for (row, col) in zip(findn(X)...)
        Z[row, col] = ifelse((X[row, col] < bottom), bottom, X[row, col])
    end
    Z
end

function relu(diff::Symbol, X::Signal)
    bottom = zero(eltype(X))
    Zd = zeros(X)
    for (row, col) in zip(findn(X)...)
        Zd[row, col] = ifelse((X[row, col] < bottom), bottom, one(eltype(X)))
    end
    Zd
end

# --------------------------------------------------
immutable LeakyRelu <: ActivationFunc
    α::Float32
    LeakyRelu() = new(0.05f0)
end

@inline function forward(f::LeakyRelu, x::Float32)
    ifelse(x < 0.0f0, x * f.α, x)
end

@inline function backward(f::LeakyRelu, x::Float32)
    ifelse(x < 0.0f0, f.α, 1.0f0)
end

# --------------------------------------------------
immutable Softmax <: ActivationFunc
end

function forward!(f::Softmax, X::Signal, Z::Signal)
    @inbounds for n in 1:size(X,2)
        x     = X[:,n]
        exp_x = exp(x .- maximum(x))
        total = sum(exp_x)
        @inbounds for k in 1:size(Z,1)
            Z[k,n] = exp_x[k] / total
        end
    end
    Z
end

function backward!(f::Softmax, X::Signal, Z::Signal)
    copy!(Z, X)
end

# --------------------------------------------------
function forward!(f::ActivationFunc, U::Signal, Z::Signal)
    @inbounds for col=1:size(Z,2), row=1:size(Z,1)
        Z[row, col] = forward(f, U[row,col])
    end
    Z
end

function backward!(f::ActivationFunc, U::Signal, Z::Signal)
    @inbounds for col=1:size(Z,2), row=1:size(Z,1)
        Z[row, col] = backward(f, U[row,col])
    end
    Z
end

function forward!(f::ActivationFunc, U::AbstractArray{Float32,1}, Z::AbstractArray{Float32,1})
    @inbounds for i=1:length(Z)
        Z[i] = forward(f, U[i])
    end
    Z
end

function backward!(f::ActivationFunc, U::AbstractArray{Float32,1}, Z::AbstractArray{Float32,1})
    @inbounds for i=1:length(Z)
        Z[i] = backward(f, U[i])
    end
    Z
end


text:
module ParametricInversionTestUtils

include("testfuncs.jl")

end # module


text:
# abstracttensormap.jl
#
# This file defines an AbstractTensorMap type for creating linear maps or transformations
# between tensors. Concrete subtypes should provide a method domain and codomain that describe
# the vector space in which the input and output tensor are living, and a method A_mul_B!
# (and optionally Ac_mul_B!) for describing the action of the map on a tensor. AbstractTensorMap
# provides functionality for acting multiplication with Julia's standard Vector type.

abstract AbstractTensorMap{S<:IndexSpace,P<:TensorSpace,T}

spacetype{S}(::AbstractTensorMap{S})=S
spacetype{S<:IndexSpace}(::Type{AbstractTensorMap{S}})=S
spacetype{S<:IndexSpace,P<:TensorSpace}(::Type{AbstractTensorMap{S,P}})=S
spacetype{S<:IndexSpace,P<:TensorSpace,T}(::Type{AbstractTensorMap{S,P,T}})=S
spacetype{TM<:AbstractTensorMap}(::Type{TM})=spacetype(super(TM))

tensortype{S,P}(::AbstractTensorMap{S,P})=P
tensortype{S<:IndexSpace,P<:TensorSpace}(::Type{AbstractTensorMap{S,P}})=P
tensortype{S<:IndexSpace,P<:TensorSpace,T}(::Type{AbstractTensorMap{S,P,T}})=P
tensortype{TM<:AbstractTensorMap}(::Type{TM})=spacetype(super(TM))

Base.eltype{S,P,T}(::AbstractTensorMap{S,P,T}) = T
Base.eltype{S<:IndexSpace,P<:TensorSpace,T}(::Type{AbstractTensorMap{S,P,T}})=T
Base.eltype{TM<:AbstractTensor}(::Type{TM}) = eltype(super(TM))

domain(A::AbstractTensorMap)=throw(MethodError(domain,(AbstractTensorMap))) # this should be implemented by subtypes; defined here to allow import
codomain(A::AbstractTensorMap)=throw(MethodError(codomain,(AbstractTensorMap))) # this should be implemented by subtypes; defined here to allow import

# Generate non-mutating multiplication methods for maps that only implement mutating methods
*{S,P}(A::AbstractTensorMap{S,P},x::AbstractTensor{S,P})=(x in domain(A) ? Base.A_mul_B!(similar(x,promote_type(eltype(A),eltype(x)),codomain(A)),A,x) : throw(SpaceError("Tensor not in domain of map")))
Base.At_mul_B{S,P}(A::AbstractTensorMap{S,P},x::AbstractTensor{S,P})=(x in dual(codomain(A)) ? Base.A_mul_B!(similar(x,promote_type(eltype(A),eltype(x)),dual(domain(A))),A,x) : throw(SpaceError("Tensor not in domain of map")))
Base.Ac_mul_B{S,P}(A::AbstractTensorMap{S,P},x::AbstractTensor{S,P})=(x in dual(conj(codomain(A)))s ? Base.A_mul_B!(similar(x,promote_type(eltype(A),eltype(x)),dual(conj(domain(A)))),A,x) : throw(SpaceError("Tensor not in domain of map")))

# Other properties
dual(A::AbstractTensorMap)=transpose(A)
adjoint{S<:HilbertSpace}(A::AbstractTensorMap{S})=(domain(A)==dual(conj(codomain(A))) ? ctranspose(A) : throw(SpaceError("Not an operator, i.e. domain != codomain")))

# The following methods allow to multiply AbstractTensorMap objects with standard Julia Vector objects, and to use it in
# general methods such as eigs, ...
Base.size(A::AbstractTensorMap)=(dim(codomain(A)),dim(domain(A)))
Base.size(A::AbstractTensorMap,n::Int)=(n==1 ? dim(codomain(A)) : (n==2 ? dim(domain(A)) : 1))

Base.isreal{S<:IndexSpace,P<:TensorSpace,T<:Real}(::AbstractTensorMap{S,P,T}) = true
Base.isreal(::AbstractTensorMap) = false # default assumption
Base.issym(::AbstractTensorMap)=false # default assumption
Base.ishermitian(A::AbstractTensorMap)=false # default assumption
Base.isposdef(::AbstractTensorMap)=false # default assumption

Base.A_mul_B!(y::Vector,A::AbstractTensorMap,x::Vector)=vec(Base.A_mul_B!(tensor(y,codomain(A)), A, tensor(x,domain(A))))
Base.At_mul_B!(y::Vector,A::AbstractTensorMap,x::Vector)=vec(Base.At_mul_B!(tensor(y,dual(domain(A))), A, tensor(x,dual(codomain(A)))))
Base.Ac_mul_B!(y::Vector,A::AbstractTensorMap,x::Vector)=vec(Base.Ac_mul_B!(tensor(y,conj(dual(domain(A)))), A, tensor(x,conj(dual(codomain(A))))))

*(A::AbstractTensorMap,x::Vector)=Base.A_mul_B!(similar(x,promote_type(eltype(A),eltype(x)),dim(codomain(A))), A, x)
Base.At_mul_B(A::AbstractTensorMap,x::Vector)=Base.At_mul_B!(similar(x,promote_type(eltype(A),eltype(x)),dim(domain(A))), A, x)
Base.Ac_mul_B(y::Vector,A::AbstractTensorMap,x::Vector)=Base.Ac_mul_B!(similar(x,promote_type(eltype(A),eltype(x)),dim(domain(A))), A, x)

# Compute full matrix representation
function Base.full(A::AbstractTensorMap)
    M,N=size(A)
    T=eltype(A)
    mat=zeros(T,(M,N))
    v=Array(T,N)
    w=Array(T,M)
    for j=1:N
        fill!(v,zero(T))
        v[j]=one(T)
        x=tensor(v,domain(A))
        y=tensor(w,codomain(A))
        A_mul_B!(y,A,x)
        mat[:,j]=w
    end
    return mat
end


text:
# ----------------------------------------------------------
# -------------- The real magic happens here ---------------
# ----------------------------------------------------------
# --- Allocate buffers for operator outputs ----------------
# ------------------- and create the aggragated function ---
# ----------------------------------------------------------

# This WeakKeyDict stores all already allocated buffers. It serves like a "L2-cache": when newBuffer
# function cannot find a proper buffer in storage ("L1-cache"), then it tries to find one in bufferPool.
# And as it is a WeakKeyDict, it allows Garbage Collector to release buffers not used any more.
const bufferPool = Array{WeakKeyDict{Buffer, Int}, 1}(undef, 0)

const bufferCounter = Counter(Threads.Atomic{Int}(1))
getNextBufferNum() = Threads.atomic_add!(bufferCounter.num, 1)

# Special constructor of Buffer type
# We get the list of all previously allocated buffers (storage) and first we try to find
# there a buffer with appropriate size which is also marked as available.
# If we find a good one, we don't allocate a new one.
function newBuffer(datatype::Type, new_size::Tuple{Vararg{Int}}, storage::Vector{Buffer})
    for i in range(length(storage), stop=1, step=-1)
        if length(storage[i].buffer) == prod(new_size) && eltype(storage[i].buffer) == datatype && storage[i].available
            return size(storage[i].buffer) == new_size ?
                storage[i] :
                Buffer(reshape(storage[i].buffer, new_size), storage[i].name, storage[i].number, true)
        end
    end
    new_buffer = nothing
    for (buffer, bufferLength) in bufferPool[Threads.threadid()]
        if bufferLength == prod(new_size) && eltype(buffer.buffer) == datatype && buffer.available
            new_buffer = size(buffer.buffer) == new_size ?
                buffer :
                Buffer(reshape(buffer.buffer, new_size), buffer.name, buffer.number, true)
            info("Re-using $(buffer.number) from bufferPool")
            break
        end
    end
    if new_buffer isa Nothing
        number = getNextBufferNum()
        name = "buffer"*string(number)
        info("Allocation of $name, size: $(new_size)")
        new_buffer = Buffer(Array{datatype}(undef, new_size), name, number, true)
        push!(bufferPool[Threads.threadid()], new_buffer => length(new_buffer.buffer))
    end
    push!(storage, new_buffer)
    new_buffer
end

#  --- getPlan functions ---
# Traverse the computational tree recursively, and tries to allocate as few buffer
# as it is (safely) possible. The goal is to get three values:
#   - A function that accepts a buffer and an input vector, it has identical functionality
#        as the operation defined by the node, and it stores the result in the given buffer
#   - A buffer which is capable of store the output of operation
#   - A string representation of the function created
#
# How is the function created that aggregates the functionality of combined operators?
# When processing a node in the computational tree, we create a new function that invokes the children nodes' aggregated functions. Trivially, it means "x -> OP1.forw(OP2.forw(x))" in case of "OP1 * OP2", and "x -> OP1.forw(x) + OP2.forw(x)" in case of "OP1 + OP2". If the node is marked as adjoint, then roles of right and left nodes are switched and it toggles their adjoint flag by xor operation (⊻); but adjoint of addition/substraction cannot by calculated that way, so an error message is thrown.
#
# How the buffers are handled?
# The buffers and passed up and down in the computational tree. Buffers descending in the tree are always the ones that are capable to store the value of output of the parent node. If it also fits the requirents of the children nodes, then the children nodes will use it to store their intermediate results. If they can't use it, they allocate their own buffers to store their intermediate results. There are three main cases:
#  - In case of multiplication of operators, the parent node passes down the buffer received from above to its right children node initiating the traversion of right node. When the traversion of right node is done, the parent receives a buffer from the right node, which will hold the output value of right node operation. Then the parent node passes the same buffer to the left node, which was previously passed to the right node, and initiates the traversion of left node. If any of the children node return buffer different from the one they received, then that buffer will be captured in the closure.
#  - In case of addition/substraction of operators, the parent allocates first a buffer for its input (to avoid calculating it two times by left and right children), and then initiates the traversal of left node. The buffer returned from the left node might be used by the right node, if exactly one of left and right is mutating, otherwise the result of left node would be overwritten by the right node before the addition. Then the closures are created accordingly. (One additional trick: if the left is mutating, but the right is not, then they are switched.)
#  - In case of leaf is reached, allocates a new buffer, if the size of the received one doesn't fit its requirements, otherwise use the provided one.

# "In case of multiplication of operators..."
function getPlanMul(FO::FunctionOperatorComposite, buffer::Buffer, adjoint::Bool, inside::String,
        storage::Vector{Buffer})
    left, right = adjoint ? (FO.right, FO.left) : (FO.left, FO.right)
    rFunc, rBuffer, rText = getPlan(right, buffer, adjoint ⊻ right.adjoint, inside, storage)
    lFunc, lBuffer, lText = getPlan(left, buffer, adjoint ⊻ left.adjoint, rText, storage)
    if rBuffer.name == buffer.name == lBuffer.name
        (buffer, x) -> lFunc(buffer, rFunc(buffer, x)), lBuffer, lText
    elseif lBuffer.name == buffer.name != rBuffer.name
        (buffer, x) -> lFunc(buffer, rFunc(rBuffer.buffer, x)), lBuffer, lText
    #elseif lBuffer.name != buffer.name == rBuffer.name
        # This case never happens
        #((buffer, x) -> lFunc(lBuffer.buffer, rFunc(buffer, x)), lBuffer, lText)
    else
        (buffer, x) -> lFunc(lBuffer.buffer, rFunc(rBuffer.buffer, x)), lBuffer, lText
    end
end

# Some code generation macro to avoid repetitive code in getPlanAddSub
macro createReturnValue(left, right, op)
    lBuf = left == :buffer ? :buffer : :(lBuffer.buffer)
    rBuf = right == :buffer ? :buffer : :(rBuffer.buffer)
    name = left == :buffer ? :(buffer.name) : :(lBuffer.name)
    esc(quote
        #if FO.left.ismutating && !FO.right.ismutating
        #   if $op == +
        #       ((buffer, x) -> begin
        #            inBuff.buffer .= x;
        #            broadcast!(+, $lBuf, rFunc($rBuf, inBuff.buffer), lFunc($lBuf, inBuff.buffer))
        #        end, lBuffer, "($(inBuff.name) .= $inside; broadcast!(+, $($name), $rText, $lText))")
        #    else
        #       ((buffer, x) -> begin
        #            inBuff.buffer .= x;
        #            broadcast!((x,y) -> y-x, $lBuf, rFunc($rBuf, inBuff.buffer), lFunc($lBuf, inBuff.buffer))
        #       end, lBuffer, "($(inBuff.name) .= $inside; broadcast!((x,y) -> y-x, , $($name), $rText, $lText))")
        #    end
        #else
           ((buffer, x) -> begin
                inBuff.buffer .= x;
                broadcast!($op, $lBuf, lFunc($lBuf, inBuff.buffer), rFunc($rBuf, inBuff.buffer))
            end, lBuffer, "($(inBuff.name) .= $inside; broadcast!($op, $($name), $lText, $rText))")
        #end
    end)
end

# "In case of addition/substraction of operators..."
function getPlanAddSub(FO::FunctionOperatorComposite, buffer::Buffer, adjoint::Bool, inside::String,
        op::Symbol, storage::Vector{Buffer})
    adjoint && error("Sorry, I don't know how to calculate the adjoint of $(FO.name)")
    buffer.available = false
    inBuff = newBuffer(eltype(FO), FO.inDims, storage)
    inBuff.available = false
    lFunc, lBuffer, lText = getPlan(FO.left, buffer, adjoint ⊻ FO.left.adjoint, inBuff.name, storage)
    lBuffer.available = false
    inBuff.available = true
    rBuffer = #FO.left.ismutating && FO.right.ismutating ?
        newBuffer(eltype(FO), FO.outDims, storage) #: lBuffer
    rFunc, rBuffer, rText = getPlan(FO.right, rBuffer, adjoint ⊻ FO.right.adjoint, inBuff.name, storage)
    lBuffer.available = buffer.available = true
    if rBuffer.name == buffer.name == lBuffer.name
        op == :+ ? @createReturnValue(buffer, buffer, +) : @createReturnValue(buffer, buffer, -)
    elseif lBuffer.name == buffer.name != rBuffer.name
        op == :+ ? @createReturnValue(buffer, rBuffer, +) : @createReturnValue(buffer, rBuffer, -)
    #elseif lBuffer.name != buffer.name == rBuffer.name
        # This case never happens
        #op == :+ ? @createReturnValue(lBuffer, buffer, +) : @createReturnValue(lBuffer, buffer, -)
    else
        op == :+ ? @createReturnValue(lBuffer, rBuffer, +) : @createReturnValue(lBuffer, rBuffer, -)
    end
end

# Pretty much self-explanatory...
function getPlan(FO::FunctionOperatorComposite, buffer::Buffer, adjoint::Bool, inside::String,
        storage::Vector{Buffer})
    if FO.operator == :*
        getPlanMul(FO, buffer, adjoint, inside, storage)
    elseif FO.operator in (:+, :-)
        getPlanAddSub(FO, buffer, adjoint, inside, FO.operator, storage)
    else
        error("Unknown operator: $(FO.operator)")
    end
end

# "In case of leaf is reached..."
function getPlan(FO::FunctionOperator, buffer::Buffer, adjoint::Bool, inside::String,
        storage::Array{Buffer,1})
    # We need to check "adjoint != FO.adjoint" because adjoint operation propagates downward in the computation tree.
    # Adjoint symbol switches inDims and outDims of the top-level
    # FunctionOperatorComposite object, but leaves intact its descendants.
    # Thus, sometimes inDims and outDims are already switched, but sometimes not,
    # that's why we need this awkward expression: adjoint != FO.adjoint ? FO.inDims : FO.outDims
    (inDims, outDims) = adjoint != FO.adjoint ? (FO.outDims, FO.inDims) : (FO.inDims, FO.outDims)
    size(buffer.buffer) != outDims && (buffer = newBuffer(eltype(FO), outDims, storage))
    if checkTwoInputs(FO.forw)
        text = FO.scaling ?
            "broadcast!(*, $(buffer.name), $(adjoint ? conj(FO.getScale()) : FO.getScale()), $inside)" :
            FO.name*(adjoint ? ".backw" : ".forw")*"($(buffer.name), $inside)"
        (adjoint ? (b,x) -> begin FO.backw(b,x); b; end : (b,x) -> begin FO.forw(b,x); b; end, buffer, text)
    else
        text = buffer.name*" .= "*FO.name*(adjoint ? ".backw" : ".forw")*"($inside)"
        (adjoint ?
            (b, x) -> b .= reshape(FO.backw(x), size(b)) :
            (b, x) -> b .= reshape(FO.forw(x), size(b)),
        buffer, text)
    end
end


text:
"""
    SweepOutput

Output generator for sweep information.
"""
abstract type SweepOutput end

"""
    init(::SweepOutput)

Initialize the `SweepOutput` before sweeping.
"""
init(::SweepOutput) = nothing

"""
    step(::SweepOutput, ::SweepState)

Update the `SweepOutput` with `SweepState` after a single step.
"""
step(::SweepOutput, ::SweepState) = nothing

"""
    sweep(::SweepOutput, ::SweepState)

Update the `SweepOutput` with `SweepState` after an entire sweep.
"""
sweep(::SweepOutput, ::SweepState) = nothing

"""
    done(::SweepOutput)

Finalize the `SweepOutput` before sweeping.
"""
done(::SweepOutput) = nothing

"""
    SweepOutputFile

`SweepOutput` that writes to a file.
"""
struct SweepOutputFile <: SweepOutput
    io::IOStream
end

"""
    SweepOutputFile(path::String)

Create a `SweepOutputFile` that will write to `path`.
"""
SweepOutputFile(path::String) = SweepOutputFile(open(path, "w"))

function init(so::SweepOutputFile)
    println(so.io, "# sweep_num energy dH2 SvN mps_max_bond_dim duration")
    flush(so.io)

    nothing
end

function sweep(so::SweepOutputFile, state::SweepState)
    items = Any[state.sweep_num, state.energy, state.dH2,
                S_vn(state.middle_eigvals), state.mps_max_bond_dim,
                state.duration]
    println(so.io, join(items, " "))
    flush(so.io)

    nothing
end

"""
    SweepOutputDynamic

`SweepOutput` that shows sweep information in real time.
"""
mutable struct SweepOutputDynamic <: SweepOutput
    "Minimum time (in seconds) between outputs."
    time_step::Float64
    "Time of last output (in seconds)."
    time_prev::Float64
end

"""
    SweepOutputDynamic(time_step::Float64=0.1)

Create a `SweepOutputDynamic` that will output no more than once every
`time_step` seconds.
"""
SweepOutputDynamic(time_step::Float64=0.1) = SweepOutputDynamic(time_step, 0.0)

function init(::SweepOutputDynamic)
    println()

    for _ in 1:7
        println()
    end

    nothing
end

function make_string(::SweepOutputDynamic, state::SweepState{L};
                     sweep=false) where {L}
    io = IOBuffer()

    for _ in 1:7
        # Go up, clear line.
        print(io, "\u1b[A\u1b[K")
    end

    # Don't exceed 80 columns.
    if L > 78
        scaling = Int(ceil((L-1)/77))
        L_eff = Int(ceil((L-1)/scaling))+1
        site_eff = div(state.site-1, scaling)+1
    else
        L_eff = L
        site_eff = state.site
    end

    print(io, "[")
    for j in 1:(site_eff-1)
        printstyled(io, "-"; color=:blue)
    end
    if sweep
        printstyled(io, "--"; color=:blue)
    else
        if state.dir == Right
            printstyled(io, ">>"; color=:red)
        elseif state.dir == Left
            printstyled(io, "<<"; color=:red)
        end
    end
    for j in (site_eff+2):L_eff
        printstyled(io, "-"; color=:blue)
    end
    print(io, "]\n")

    if sweep
        @printf(io, "done sweep %d\n", state.sweep_num)
    else
        @printf(io, "on sweep %d at site %d of %d\n",
                state.sweep_num, state.site, L)
    end
    @printf(io, "maximum MPS bond dimension is %d\n", state.mps_max_bond_dim)
    if isnan(state.duration)
        println(io, "")
    else
        @printf(io, "previous sweep took %.3f s\n", state.duration)
    end
    if isnan(state.energy)
        println(io, "")
    else
        @printf(io, "energy is %.5g\n", state.energy)
    end
    if isempty(state.middle_eigvals)
        println(io, "")
    else
        @printf(io, "von Neumann entropy across middle bond is %.5g\n",
                S_vn(state.middle_eigvals))
    end
    if isnan(state.dH2)
        println(io, "")
    else
        @printf(io, "Hamiltonian fluctuation (dH2) is %.5g\n", state.dH2)
    end

    io |> take! |> String
end

function step(so::SweepOutputDynamic, state::SweepState)
    time_now = time()
    time_now >= so.time_prev + so.time_step || return

    print(make_string(so, state; sweep=false))

    so.time_prev = time_now

    nothing
end

function sweep(so::SweepOutputDynamic, state::SweepState)
    time_now = time()

    print(make_string(so, state; sweep=true))

    # Don't rate limit the end-of-sweep updates, but still update the time.
    so.time_prev = time_now

    nothing
end

function done(::SweepOutputDynamic)
    println()

    nothing
end