text:
@doc """`pole(sys)`

Compute the poles of system `sys`.""" ->
pole(sys::StateSpace) = eig(sys.A)[1]
pole(sys::TransferFunction) = [map(pole, sys.matrix)...;]
pole(sys::SisoTf) = error("pole is not implemented for type $(typeof(sys))")

@doc """`dcgain(sys)`

Compute the gain of SISO system `sys`.""" ->
function dcgain(sys::StateSpace, zs::Vector=tzero(sys))
    !issiso(sys) && error("Gain only defined for siso systems")
    return ( sys.C*(-sys.A\sys.B) + sys.D)[1]
end

@doc """`dcgain(sys)`

Compute the dcgain of system `sys`.

equal to G(0) for continuous-time systems and G(1) for discrete-time systems.""" ->
function dcgain(sys::TransferFunction)
    !issiso(sys) && error("Gain only defined for siso systems")
    if sys.Ts > 0
        return map(s->sum(s.num.a)/sum(s.den.a), sys.matrix)
    end

    return map(s -> s.num[end]/s.den[end], sys.matrix)
end

@doc """`markovparam(sys, n)`

Compute the `n`th markov parameter of state-space system `sys`. This is defined
as the following:

`h(0) = D`

`h(n) = C*A^(n-1)*B`""" ->
function markovparam(sys::StateSpace, n::Integer)
    n < 0 && error("n must be >= 0")
    return n == 0 ? sys.D : sys.C * sys.A^(n-1) * sys.B
end

@doc """`z, p, k = zpkdata(sys)`

Compute the zeros, poles, and gains of system `sys`.

### Returns
`z` : Matrix{Vector{Complex128}}, (ny x nu)

`p` : Matrix{Vector{Complex128}}, (ny x nu)

`k` : Matrix{Float64}, (ny x nu)""" ->
function zpkdata(sys::LTISystem)
    ny, nu = size(sys)
    zs = Array(Vector{Complex128}, ny, nu)
    ps = Array(Vector{Complex128}, ny, nu)
    ks = Array(Float64, ny, nu)
    for j = 1:nu
        for i = 1:ny
            zs[i, j], ps[i, j], ks[i, j] = _zpk_kern(sys, i, j)
        end
    end
    return zs, ps, ks
end

function _zpk_kern(sys::StateSpace, iy::Int, iu::Int)
    A, B, C = struct_ctrb_obsv(sys.A, sys.B[:, iu:iu], sys.C[iy:iy, :])
    D = sys.D[iy:iy, iu:iu]
    z = tzero(A, B, C, D)
    nx = size(A, 1)
    nz = length(z)
    k = nz == nx ? D[1] : (C*(A^(nx - nz - 1))*B)[1]
    return z, eigvals(A), k
end

function _zpk_kern(sys::TransferFunction, iy::Int, iu::Int)
    s = sys.matrix[iy, iu]
    return _zpk_kern(s)
end

function _zpk_kern(s::SisoRational)
  return roots(s.num), roots(s.den), s.num[1]/s.den[1]
end

function _zpk_kern(s::SisoZpk)
    return s.z, s.p, s.k
end

@doc """`Wn, zeta, ps = damp(sys)`

Compute the natural frequencies, `Wn`, and damping ratios, `zeta`, of the
poles, `ps`, of `sys`""" ->
function damp(sys::LTISystem)
    ps = pole(sys)
    if !iscontinuous(sys)
        Ts = sys.Ts == -1 ? 1 : sys.Ts
        ps = log(ps)/Ts
    end
    Wn = abs(ps)
    order = sortperm(Wn)
    Wn = Wn[order]
    ps = ps[order]
    ζ = -cos(angle(ps))
    return Wn, ζ, ps
end

@doc """`dampreport(sys)`

Display a report of the poles, damping ratio, natural frequency, and time
constant of the system `sys`""" ->
function dampreport(io::IO, sys::LTISystem)
    Wn, zeta, ps = damp(sys)
    t_const = 1./(Wn.*zeta)
    header =
    ("|     Pole      |   Damping     |   Frequency   | Time Constant |\n"*
    "|               |    Ratio      |   (rad/sec)   |     (sec)     |\n"*
    "+---------------+---------------+---------------+---------------+")
    println(io, header)
    for i=eachindex(ps)
        p, z, w, t = ps[i], zeta[i], Wn[i], t_const[i]
        @printf(io, "|  %-13.3e|  %-13.3e|  %-13.3e|  %-13.3e|\n", p, z, w, t)
    end
end
dampreport(sys::LTISystem) = dampreport(STDOUT, sys)

@doc """`tzero(sys)`

Compute the invariant zeros of the system `sys`. If `sys` is a minimal
realization, these are also the transmission zeros.""" ->
function tzero(sys::TransferFunction)
    if issiso(sys)
        return tzero(sys.matrix[1,1])
    else
        return tzero(ss(sys))
    end
end

# Implements the algorithm described in:
# Emami-Naeini, A. and P. Van Dooren, "Computation of Zeros of Linear
# Multivariable Systems," Automatica, 18 (1982), pp. 415–430.
#
# Note that this returns either Vector{Complex64} or Vector{Float64}
tzero(sys::StateSpace) = tzero(sys.A, sys.B, sys.C, sys.D)
function tzero(A::Matrix{Float64}, B::Matrix{Float64}, C::Matrix{Float64},
        D::Matrix{Float64})
    # Balance the system
    A, B, C = balance_statespace(A, B, C)

    # Compute a good tolerance
    meps = 10*eps()*norm([A B; C D])
    A, B, C, D = reduce_sys(A, B, C, D, meps)
    A, B, C, D = reduce_sys(A', C', B', D', meps)
    if isempty(A)   return Float64[]    end

    # Compress cols of [C D] to [0 Df]
    mat = [C D]
    # To ensure type-stability, we have to annote the type here, as qrfact
    # returns many different types.
    W = full(qrfact(mat')[:Q], thin=false)::Matrix{Float64}
    W = flipdim(W,2)
    mat = mat*W
    if fastrank(mat', meps) > 0
        nf = size(A, 1)
        m = size(D, 2)
        Af = ([A B] * W)[1:nf, 1:nf]
        Bf = ([eye(nf) zeros(nf, m)] * W)[1:nf, 1:nf]
        zs = eig(Af, Bf)[1]
    else
        zs = Float64[]
    end
    return zs
end

"""
Implements REDUCE in the Emami-Naeini & Van Dooren paper. Returns transformed
A, B, C, D matrices. These are empty if there are no zeros.
"""
function reduce_sys(A::Matrix{Float64}, B::Matrix{Float64}, C::Matrix{Float64}, D::Matrix{Float64}, meps::Float64)
    Cbar, Dbar = C, D
    if isempty(A)
        return A, B, C, D
    end
    while true
        # Compress rows of D
        U = full(qrfact(D)[:Q], thin=false)::Matrix{Float64}
        D = U'*D
        C = U'*C
        sigma = fastrank(D, meps)
        Cbar = C[1:sigma, :]
        Dbar = D[1:sigma, :]
        Ctilde = C[(1 + sigma):end, :]
        if sigma == size(D, 1)
            break
        end

        # Compress columns of Ctilde
        V = full(qrfact(Ctilde')[:Q], thin=false)::Matrix{Float64}
        V = flipdim(V,2)
        Sj = Ctilde*V
        rho = fastrank(Sj', meps)
        nu = size(Sj, 2) - rho

        if rho == 0
            break
        elseif nu == 0
            # System has no zeros, return empty matrices
            A = B = Cbar = Dbar = Float64[]
            break
        end
        # Update System
        n, m = size(B)
        Vm = [V zeros(n, m); zeros(m, n) eye(m)]
        if sigma > 0
            M = [A B; Cbar Dbar]
            Vs = [V' zeros(n, sigma) ; zeros(sigma, n) eye(sigma)]
        else
            M = [A B]
            Vs = V'
        end
        sigma, rho, nu
        M = Vs * M * Vm
        A = M[1:nu, 1:nu]
        B = M[1:nu, (nu + rho + 1):end]
        C = M[(nu + 1):end, 1:nu]
        D = M[(nu + 1):end,  (nu + rho + 1):end]
    end
    return A, B, Cbar, Dbar
end

# Determine the number of non-zero rows, with meps as a tolerance. For an
# upper-triangular matrix, this is a good proxy for determining the row-rank.
function fastrank(A::Matrix{Float64}, meps::Float64)
    n, m = size(A, 1, 2)
    if n*m == 0     return 0    end
    norms = Array(Float64, n)
    for i = 1:n
        norms[i] = norm(A[i, :])
    end
    mrank = sum(norms .> meps)
    return mrank
end

"""
`ωgₘ, gₘ, ωϕₘ, ϕₘ = margin{S<:Real}(sys::LTISystem, w::AbstractVector{S}; full=false, allMargins=false)`

returns frequencies for gain margins, gain margins, frequencies for phase margins, phase margins

If `!allMargins`, return only the smallest margin

If `full` return also `fullPhase`

"""
function margin{S<:Real}(sys::LTISystem, w::AbstractVector{S}; full=false, allMargins=false)
    ny, nu = size(sys)
    vals = (:wgm, :gm, :wpm, :pm, :fullPhase)
    if allMargins
        for val in vals
            eval(:($val = Array{Array{Float64,1}}($ny,$nu)))
        end
    else
        for val in vals
            eval(:($val = Array{Float64,2}($ny,$nu)))
        end
    end
    for j=1:nu
        for i=1:ny
            wgm[i,j], gm[i,j], wpm[i,j], pm[i,j], fullPhase[i,j] = sisomargin(sys[i,j], w, full=true, allMargins=allMargins)
        end
    end
    if full
        wgm, gm, wpm, pm, fullPhase
    else
        wgm, gm, wpm, pm
    end
end

"""
`ωgₘ, gₘ, ωϕₘ, ϕₘ = sisomargin{S<:Real}(sys::LTISystem, w::AbstractVector{S}; full=false, allMargins=false)`

returns frequencies for gain margins, gain margins, frequencies for phase margins, phase margins
"""
function sisomargin{S<:Real}(sys::LTISystem, w::AbstractVector{S}; full=false, allMargins=false)
    ny, nu = size(sys)
    if ny !=1 || nu != 1
        error("System must be SISO, use `margin` instead")
    end
    mag, phase, w = bode(sys, w)
    wgm, = _allPhaseCrossings(w, phase)
    gm = similar(wgm)
    for i = eachindex(wgm)
        gm[i] = 1./abs(freqresp(sys,[wgm[i]])[1][1])
    end
    wpm, fi = _allGainCrossings(w, mag)
    pm = similar(wpm)
    for i = eachindex(wpm)
        pm[i] = mod(rad2deg(angle(freqresp(sys,[wpm[i]])[1][1])),360)-180
    end
    if !allMargins #Only output the smallest margins
        gm, idx = findmin([gm;Inf])
        wgm = [wgm;NaN][idx]
        fi = [fi;NaN][idx]
        pm, idx = findmin([abs(pm);Inf])
        wpm = [wpm;NaN][idx]
        if full
            if !isnan(fi) #fi may be NaN, fullPhase is a scalar
                fullPhase = interpolate(fi, phase)
            else
                fullPhase = NaN
            end
        end
    else
        if full #We know that all values are defined and fullPhase is a vector
            fullPhase = interpolate(fi, phase)
        end
    end
    if full
        wgm, gm, wpm, pm, fullPhase
    else
        wgm, gm, wpm, pm
    end
end
margin(system::LTISystem; kwargs...) =
margin(system, _default_freq_vector(system, :bode); kwargs...)
#margin(sys::LTISystem, args...) = margin(LTISystem[sys], args...)

# Interpolate the values in "list" given the floating point "index" fi
function interpolate(fi, list)
    fif = floor(Integer, fi)
    fic = ceil(Integer, fi)
    list[fif]+mod(fi,1).*(list[fic]-list[fif])
end

function _allGainCrossings(w, mag)
    _findCrossings(w,mag.>1,mag-1)
end

function _allPhaseCrossings(w, phase)
    #Calculate numer of times real axis is crossed on negative side
    n =  Array{Float64,1}(length(w)) #Nbr of crossed
    ph = Array{Float64,1}(length(w)) #Residual
    for i = eachindex(w) #Found no easier way to do this
        n[i], ph[i] = fldmod(phase[i]+180,360)#+180
    end
    _findCrossings(w, n, ph)
end

function _findCrossings(w, n, res)
    wcross = Array{Float64,1}()
    tcross = Array{Float64,1}()
    for i in 1:(length(w)-1)
        if res[i] == 0
            wcross = [wcross; w[i]]
            tcross = [tcross; i]
        elseif n[i] != n[i+1]
            #Interpolate to approximate crossing
            t = res[i]/(res[i]-res[i+1])
            tcross = [tcross; i+t]
            wt = w[i]+t*(w[i+1]-w[i])
            wcross = [wcross; wt]
        end
    end
    if res[end] == 0 #Special case if multiple points
        wcross = [wcross; w[end]]
        tcross = [tcross; length(w)]
    end
    wcross, tcross
end

@doc """`dₘ = delaymargin(G::LTISystem)`

Only supports SISO systems""" ->
function delaymargin(G::LTISystem)
    # Phase margin in radians divided by cross-over frequency in rad/s.
    if G.nu + G.ny > 2
        error("delaymargin only supports SISO systems")
    end
    m   = margin(G,allMargins=true)
    ϕₘ,i= findmin(m[4])
    ϕₘ *= π/180
    ωϕₘ = m[3][i]
    dₘ  = ϕₘ/ωϕₘ
    if G.Ts > 0
        dₘ /= G.Ts # Give delay margin in number of sample times, as matlab does
    end
    dₘ
end

@doc """`S,D,N,T = gangoffour(P,C)`, `gangoffour(P::AbstractVector,C::AbstractVector)`

Given a transfer function describing the Plant `P` and a transferfunction describing the controller `C`, computes the four transfer functions in the Gang-of-Four.

`S = 1/(1+PC)` Sensitivity function

`D = P/(1+PC)`

`N = C/(1+PC)`

`T = PC/(1+PC)` Complementary sensitivity function

Only supports SISO systems""" ->
function gangoffour(P::TransferFunction,C::TransferFunction)
    if P.nu + P.ny + C.nu + C.ny > 4
        error("gangoffour only supports SISO systems")
    end
    S = minreal(1/(1+P*C))
    D = minreal(P*S)
    N = minreal(C*S)
    T = minreal(P*N)
    return S, D, N, T
end

function gangoffour(P::AbstractVector, C::AbstractVector)
    if P[1].nu + P[1].ny + C[1].nu + C[1].ny > 4
        error("gangoffour only supports SISO systems")
    end
    length(P) == length(C) || error("P has to be the same length as C")
    n = length(P)
    S = [minreal(1/(1+P[i]*C[i])) for i in 1:n]
    D = [minreal(P[i]*S[i]) for i in 1:n]
    N = [minreal(C[i]*S[i]) for i in 1:n]
    T = [minreal(P[i]*N[i]) for i in 1:n]
    return S, D, N, T
end

function gangoffour(P::TransferFunction, C::AbstractVector)
    gangoffour(fill(P,length(C)), C)
end

function gangoffour(P::AbstractVector, C::TransferFunction)
    gangoffour(P, fill(C,length(P)))
end

@doc """`S, D, N, T, RY, RU, RE = gangofseven(P,C,F)`

Given transfer functions describing the Plant `P`, the controller `C` and a feed forward block `F`,
computes the four transfer functions in the Gang-of-Four and the transferfunctions corresponding to the feed forward.

`S = 1/(1+PC)` Sensitivity function

`D = P/(1+PC)`

`N = C/(1+PC)`

`T = PC/(1+PC)` Complementary sensitivity function

`RY = PCF/(1+PC)`

`RU = CF/(1+P*C)`

`RE = F/(1+P*C)`

Only supports SISO systems""" ->
function gangofseven(P::TransferFunction,C::TransferFunction,F::TransferFunction)
    if P.nu + P.ny + C.nu + C.ny + F.nu + F.ny > 6
        error("gof only supports SISO systems")
    end
    S, D, N, T = gangoffour(P,C)
    RY = T*F
    RU = N*F
    RE = S*F
    return S, D, N, T, RY, RU, RE
end


text:
mutable struct ArrayPairWithLifetime
    turns_left::Int
    cached_array::CachedArray
    reference::Array
end

function epoch(
    manager::CachedArrays.CacheManager,
    arrays::Set{ArrayPairWithLifetime},
    generator::F;
    allocate = true,
    prefetch_percent = 0.20,
    evict_percent = 0.02,
    defrag_percent = 0.01,
) where {F}
    if allocate
        size, lifetime = generator()
        nelements = ceil(Int, size / sizeof(Float32))
        reference = randn(Float32, nelements)
        cached_array = CachedArrays.CachedArray{Float32}(undef, manager, nelements)
        CachedArrays.writable(cached_array) .= reference
        push!(arrays, ArrayPairWithLifetime(lifetime, cached_array, reference))
    end

    # Go through each element in the "arrays" set.
    # Decrement lifetime by 1.
    # If lifetime is at zero, check that the cached array is still the same as the reference
    # array and pop it from the `arrays` set - making it available for garbage collection.
    #
    # If it's not time to free the array, possibly call "evict" of "prefetch".
    arrays_to_delete = []
    for pair in arrays
        turns_left = pair.turns_left
        if iszero(turns_left)
            @test CachedArrays.readable(pair.cached_array) == pair.reference
            push!(arrays_to_delete, pair)
        else
            pair.turns_left -= 1

            # Try to prefetch or evict.
            if rand() < prefetch_percent
                CachedArrays.prefetch!(pair.cached_array)
            end

            if rand() < evict_percent
                CachedArrays.evict!(pair.cached_array)
            end
        end
    end

    # Maybe Defrag.
    if rand() < defrag_percent
        array_evicted = false
        for pair in arrays
            if rand() < evict_percent
                CachedArrays.evict!(pair.cached_array)
                array_evicted = true
            end
        end
        if array_evicted
            CachedArrays.defrag!(manager)
        end
    end

    for array in arrays_to_delete
        delete!(arrays, array)
    end
end

@testset "Performing Stress Test 1" begin
    target_allocations = 100_000
    target_mean_lifetime = 200
    target_mean_allocation = 8192 # bytes
    allocation_range = 8191 # bytes

    manager = CachedArrays.CacheManager(
        CachedArrays.AlignedAllocator(),
        CachedArrays.AlignedAllocator();
        localsize = 100 * target_mean_allocation,
        # Make the remote size pretty large so we don't need to worry about running out
        # of memory.
        remotesize = ceil(Int, target_allocations * target_mean_allocation),
        minallocation = 12,
    )

    # Create the generator function.
    Random.seed!(12345)
    exp_distribution = Exponential(1 / target_mean_lifetime)
    function generator()
        _size = target_mean_allocation + rand(-allocation_range:allocation_range)
        _lifetime = clamp(rand(exp_distribution), 2, 1000)
        return _size, _lifetime
    end

    @show CachedArrays.gettimer()
    CachedArrays.reset_timer!()
    arrays = Set{ArrayPairWithLifetime}()
    println("Beginning Stress Test")
    @time ProgressMeter.@showprogress 1 for _i in 1:target_allocations
        epoch(manager, arrays, generator; allocate = true)
        CachedArrays.check(manager) || error()
    end

    println("Cleaning Up Stress Test")
    @time while !isempty(arrays)
        epoch(manager, arrays, generator; allocate = false)
        CachedArrays.check(manager) || error()
    end

    @show manager
    @show CachedArrays.gettimer()
end

ext: jl


text:
function EuropeanDownAndOutCallMC(K, T, S, sig, r, div, H, N, M)

dt = T/N
nudt = (r - div - 0.5 * sig^2) * dt
sigsdt = sig * sqrt(dt)

## variables to store the sums of the option values

sum_CT = 0
sum_CT2 = 0

## outer loop is the number of simulations bigger is better

for j = 1:M

        St = S
        BARRIERCROSSED = false
    ## inner loop is simulating each price path of length N
        for i = 1:N

        St = St * exp( nudt + sigsdt * randn()) ## note, it only add the changes, does not store the path

        if St <= H

                BARRIERCROSSED = true
                break
        end

        end

    ## wraps up and stores all the values of each path

if BARRIERCROSSED == true
            CT = 0

        else

            CT = max(0, St - K)
        end

sum_CT = sum_CT + CT
sum_CT2 = sum_CT2 + CT^2

    end

  ## discounting the average call value
call_value = sum_CT / M * exp(-r*T)

SD = (sqrt(sum_CT2 - 1/M * (sum_CT^2)) * exp(-2*r*T)) / (M-1)
SE = SD / sqrt(M)

println("The Call Value is $call_value")
println("The SD is $SD")
println("The SE is $SE")

MC = [call_value, SD, SE]

end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    Techs(p::REoptInputs, s::BAUScenario)

Create a `Techs` struct for the BAUInputs
"""
function Techs(p::REoptInputs, s::BAUScenario)
    pvtechs = String[pv.name for pv in s.pvs]

    all_techs = copy(pvtechs)
    elec = copy(pvtechs)
    techs_no_turndown = copy(pvtechs)
    techs_no_curtail = String[]
    segmented_techs = String[]
    gentechs = String[]
    pbi_techs = String[]
    heating_techs = String[]
    boiler_techs = String[]
    chp_techs = String[]

    if p.s.generator.existing_kw > 0
        push!(all_techs, "Generator")
        push!(gentechs, "Generator")
        push!(elec, "Generator")
    end

    if p.s.existing_boiler.max_kw > 0
        push!(all_techs, "ExistingBoiler")
        push!(heating_techs, "ExistingBoiler")
        push!(boiler_techs, "ExistingBoiler")
    end

    fuel_burning_techs = union(gentechs, boiler_techs, chp_techs)
    thermal_techs = union(heating_techs, boiler_techs)

    Techs(
        all_techs,
        elec,
        pvtechs,
        gentechs,
        pbi_techs,
        techs_no_curtail,
        techs_no_turndown,
        segmented_techs,
        heating_techs,
        boiler_techs,
        fuel_burning_techs,
        thermal_techs,
        chp_techs
    )
end

"""
    Techs(s::Scenario)

Create a `Techs` struct for the REoptInputs.
"""
function Techs(s::Scenario)
    pvtechs = String[pv.name for pv in s.pvs]
    if length(Base.Set(pvtechs)) != length(pvtechs)
        error("PV names must be unique, got $(pvtechs)")
    end

    all_techs = copy(pvtechs)
    elec = copy(pvtechs)
    techs_no_turndown = copy(pvtechs)
    gentechs = String[]
    pbi_techs = String[]
    techs_no_curtail = String[]
    segmented_techs = String[]
    heating_techs = String[]
    boiler_techs = String[]
    chp_techs = String[]
    if s.wind.max_kw > 0
        push!(all_techs, "Wind")
        push!(elec, "Wind")
    end
    if s.generator.max_kw > 0
        push!(all_techs, "Generator")
        push!(gentechs, "Generator")
        push!(elec, "Generator")
    end

    if s.existing_boiler.max_kw > 0
        push!(all_techs, "ExistingBoiler")
        push!(heating_techs, "ExistingBoiler")
        push!(boiler_techs, "ExistingBoiler")
    end

    if "Wind" in all_techs
        append!(techs_no_turndown, ["Wind"])
    end

    if !isnothing(s.chp)
        push!(all_techs, "CHP")
        push!(elec, "CHP")
        push!(chp_techs, "CHP")
    end

    thermal_techs = union(heating_techs, boiler_techs, chp_techs)
    fuel_burning_techs = union(gentechs, boiler_techs, chp_techs)

    Techs(
        all_techs,
        elec,
        pvtechs,
        gentechs,
        pbi_techs,
        techs_no_curtail,
        techs_no_turndown,
        segmented_techs,
        heating_techs,
        boiler_techs,
        fuel_burning_techs,
        thermal_techs,
        chp_techs
    )
end

"""
    Techs(s::MPCScenario)

Create a `Techs` struct for the MPCInputs
"""
function Techs(s::MPCScenario)
    pvtechs = String[pv.name for pv in s.pvs]
    if length(Base.Set(pvtechs)) != length(pvtechs)
        error("PV names must be unique, got $(pvtechs)")
    end

    all_techs = copy(pvtechs)
    techs_no_turndown = copy(pvtechs)
    gentechs = String[]
    if s.generator.size_kw > 0
        push!(all_techs, "Generator")
        push!(gentechs, "Generator")
    end

    Techs(
        all_techs,
        all_techs,
        pvtechs,
        gentechs,
        String[],
        String[],
        techs_no_turndown,
        String[],
        String[],
        String[],
        String[],
        String[],
        String[]
    )
end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 84086dd0-cee2-11eb-0054-6daab703dc7f
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots"),
        Pkg.PackageSpec(name="PlutoUI"),
        Pkg.PackageSpec(name="Distributions"),
		Pkg.PackageSpec(name="StatsPlots"),

    ])
    using Plots, PlutoUI, Distributions, StatsPlots
end

# ╔═╡ 01c84816-938d-48f9-bd98-367cdb539baa
begin
	# Initialize the arrays with zeros
	x = zeros(100)
	y = zeros(100)

	# Initialize the starting position.
	x[1] = 5.0
	y[1] = 5.0

	# Loop over the walk and add the random noise.
	for ii=2:length(x)
		newx = x[ii-1]
		newy = y[ii-1]
		while true
			dx = rand(Normal())*0.5
			dy = rand(Normal())*0.5
			newx = x[ii-1] + dx
			newy = y[ii-1] + dy
			if (0<=newx<=10.0) && (0<=newy<=10.0)
				x[ii] = newx
				y[ii] = newy
				break
			end
		end

	end
end

# ╔═╡ 6ec0c117-9867-4669-9758-3795884f8d73
begin
	anim = @animate for i ∈ 1:length(x)
    	scatter([x[i]],[y[i]],
			xlims=(0,10),
			ylims=(0,10),
			legend=false,
			framestyle=:box)
	end

	gif(anim, "anim_fps15.gif", fps = 15)
end

# ╔═╡ d91a7845-dfc1-4daf-99d7-d4d16596d48c
begin
	N = 100
	p = 10000
	# Initialize the arrays with zeros
	xx = zeros(p,N)
	yy = zeros(p,N)

	for jj=1:N

		# Initialize the starting position.
		xx[1,jj] = rand(Uniform(0, 10))
		yy[1,jj] = rand(Uniform(0, 10))

		# Loop over the walk and add the random noise.
		for ii=2:p
			newx = xx[ii-1,jj]
			newy = yy[ii-1,jj]
			while true
				dx = rand(Normal())*0.1
				dy = rand(Normal())*0.1
				newx = xx[ii-1,jj] + dx
				newy = yy[ii-1,jj] + dy
				if (0<=newx<=10.0) && (0<=newy<=10.0)
					xx[ii,jj] = newx
					yy[ii,jj] = newy
					break
				end
			end

		end
	end
end

# ╔═╡ 1f0cdab6-e09f-482a-9447-6932ab87f5d1
begin
	diffusion = 0.5
	NF = 10
	ff = zeros(p,N)
	fF = zeros(p,NF)
	fF[1,:] .= 100   # foragers have food = 100
	intDist = 0.3 # distance at which interactions are initiated

	# Initialize the arrays with zeros
	xxF = zeros(p,NF)
	yyF = zeros(p,NF)

	for jj=1:NF

		# Initialize the starting position.
		xxF[1,jj] = 0
		yyF[1,jj] = 0

		# Loop over the walk and add the random noise.
		for ii=2:p
			newx = xxF[ii-1,jj]
			newy = yyF[ii-1,jj]
			while true
				dx = rand(Normal())*diffusion
				dy = rand(Normal())*diffusion
				newx = xxF[ii-1,jj] + dx
				newy = yyF[ii-1,jj] + dy
				if (0<=newx<=10.0) && (0<=newy<=10.0)
					xxF[ii,jj] = newx
					yyF[ii,jj] = newy
					fF[ii,jj] = fF[ii-1,jj]
					for j=1:N
						ff[ii,j] = ff[ii-1,j]
						dist = sqrt((xxF[ii-1,jj]-xx[ii-1,j])^2+(yyF[ii-1,jj]-yy[ii-1,j])^2)
						if dist<intDist
							tofill = 100-ff[ii-1,j]
							if tofill <=0
								continue
							end
							transfer = rand(Exponential(tofill*.2))
							println(transfer)
							if transfer>tofill
								transfer = tofill
							end
							if transfer>=fF[ii-1,jj]
								fF[ii,jj]=100
								#xxF[ii,jj]=0
								#yyF[ii,jj]=0
								ff[ii,j] = ff[ii-1,j] + fF[ii-1,jj]
							else
								fF[ii,jj] = fF[ii-1,jj] - transfer
								ff[ii,j] = ff[ii-1,j] + transfer
							end
						end
					end
					break
				end
			end

		end
	end
end

# ╔═╡ 2ecfb056-2486-4dc1-be45-825b38975066
xx[1,:]

# ╔═╡ 83151cec-2b31-41ca-a718-89a599494698
begin
	anim3 = @animate for i ∈ 1:p
    	scatter([xx[i,:]],[yy[i,:]],
			xlims=(0,10),
			ylims=(0,10),
			legend = false,
			framestyle=:box,
			color=:red)

			scatter!([xxF[i,:]],[yyF[i,:]],
			xlims=(0,10),
			ylims=(0,10),
			legend = false,
			framestyle=:box,
			color=:blue)

	end

	gif(anim3, "anim_fps15_3.gif", fps = 15)
end

# ╔═╡ fcb59fa7-8eb9-41a7-b794-134598abb560
plot(ff[1:1000,10])

# ╔═╡ c1c07da4-cdbd-4c78-88bb-81ea965f9e37
begin
	ff_avg = sum(ff,dims=2)/100
	plot(ff_avg)
end

# ╔═╡ fdfb1502-e071-4d89-8439-4f5130544edc
histogram(ff[8000,:],bins=20)

# ╔═╡ Cell order:
# ╠═84086dd0-cee2-11eb-0054-6daab703dc7f
# ╠═01c84816-938d-48f9-bd98-367cdb539baa
# ╠═6ec0c117-9867-4669-9758-3795884f8d73
# ╠═d91a7845-dfc1-4daf-99d7-d4d16596d48c
# ╠═1f0cdab6-e09f-482a-9447-6932ab87f5d1
# ╠═2ecfb056-2486-4dc1-be45-825b38975066
# ╠═83151cec-2b31-41ca-a718-89a599494698
# ╠═fcb59fa7-8eb9-41a7-b794-134598abb560
# ╠═c1c07da4-cdbd-4c78-88bb-81ea965f9e37
# ╠═fdfb1502-e071-4d89-8439-4f5130544edc

max_stars_repo_name: hstrey/ants

max_issues_repo_name: hstrey/ants

max_forks_repo_name: hstrey/ants


text:
__precompile__()

module BandedMatrices
using Base, Compat

import Base: getindex, setindex!, *, .*, +, .+, -, .-, ==, <, <=, >,
                >=, ./, /, .^, ^, \, transpose, showerror

import Base: convert, size

import Base.BLAS: libblas

import Base.LinAlg: BlasInt,
                    BlasReal,
                    BlasFloat,
                    BlasComplex,
                    A_ldiv_B!,
                    At_ldiv_B!,
                    Ac_ldiv_B!,
                    copy_oftype

import Base.LAPACK: gbtrs!,
                    gbtrf!

import Base: lufact
import Compat: view

export BandedMatrix,
       bandrange,
       bzeros,
       beye,
       brand,
       bones,
       bandwidth,
       BandError,
       band,
       BandRange

# BLAS/linear algebra overrides

@inline dot(x...) = Base.dot(x...)
@inline dot{T<:Union{Float64,Float32}}(M::Int,a::Ptr{T},incx::Int,b::Ptr{T},incy::Int) =
    BLAS.dot(M,a,incx,b,incy)
@inline dot{T<:Union{Complex128,Complex64}}(M::Int,a::Ptr{T},incx::Int,b::Ptr{T},incy::Int) =
    BLAS.dotc(M,a,incx,b,incy)

dotu{T<:Union{Complex64,Complex128}}(f::Vector{T},g::Vector{T}) =
    BLAS.dotu(f,g)
dotu{N<:Real}(f::Vector{Complex{Float64}},g::Vector{N}) = dot(conj(f),g)
dotu{N<:Real,T<:Number}(f::Vector{N},g::Vector{T}) = dot(f,g)

normalize!(w::AbstractVector) = scale!(w,inv(norm(w)))
normalize!{T<:BlasFloat}(w::Vector{T}) = normalize!(length(w),w)
normalize!{T<:Union{Float64,Float32}}(n,w::Union{Vector{T},Ptr{T}}) =
    BLAS.scal!(n,inv(BLAS.nrm2(n,w,1)),w,1)
normalize!{T<:Union{Complex128,Complex64}}(n,w::Union{Vector{T},Ptr{T}}) =
    BLAS.scal!(n,T(inv(BLAS.nrm2(n,w,1))),w,1)

# AbstractBandedMatrix must implement

abstract AbstractBandedMatrix{T} <: AbstractSparseMatrix{T,Int}

bandwidths(A::AbstractBandedMatrix) = bandwidth(A,1),bandwidth(A,2)
bandinds(A::AbstractBandedMatrix) = -bandwidth(A,1),bandwidth(A,2)
bandinds(A::AbstractBandedMatrix,k::Integer) = k==1?-bandwidth(A,1):bandwidth(A,2)
bandrange(A::AbstractBandedMatrix) = -bandwidth(A,1):bandwidth(A,2)
bandwidth(A::AbstractMatrix,k::Integer) = k==1?size(A,1)-1:size(A,2)-1

function getindex(A::AbstractBandedMatrix,k::Integer,j::Integer)
    if k>size(A,1) || j>size(A,2)
        throw(BoundsError(A,(k,j)))
    elseif (bandinds(A,1)≤j-k≤bandinds(A,2))
        unsafe_getindex(A,k,j)
    else
        zero(eltype(A))
    end
end

# override bandwidth(A,k) for each AbstractBandedMatrix
# override unsafe_getindex(A,k,j)

## Gives an iterator over the banded entries of a BandedMatrix

immutable BandedIterator
    m::Int
    n::Int
    l::Int
    u::Int
end

Base.start(B::BandedIterator)=(1,1)

Base.next(B::BandedIterator,state)=
    state,ifelse(state[1]==min(state[2]+B.l,B.m),
                (max(state[2]+1-B.u,1),state[2]+1),
                (state[1]+1,  state[2])
                 )

Base.done(B::BandedIterator,state)=state[2]>B.n || state[2]>B.m+B.u

Base.eltype(::Type{BandedIterator})=Tuple{Int,Int}

# Commented out since there's an error

# for m x n matrix, assuming m ≥ n
# see notes
function bandslength(m,n,l,u)
    if m ≥ n
        (1 + m - n)*n + (-l - l^2 + m + 2l*m - m^2 - n + n^2)÷2 +
            (-u + 2n*u - u^2)÷2
    else
        n*(l+1)+u*n-(u^2+u)÷2
    end
end

function Base.length(B::BandedIterator)
    if B.m ≥ B.n
        bandslength(B.m,B.n,B.l,B.u)
    else # transpose
        bandslength(B.n,B.m,B.u,B.l)
    end
end

# returns an iterator of each index in the banded part of a matrix
eachbandedindex(B) = BandedIterator(size(B,1),size(B,2),bandwidth(B,1),bandwidth(B,2))

##
# Represent a banded matrix
# [ a_11 a_12
#   a_21 a_22 a_23
#   a_31 a_32 a_33 a_34
#        a_42 a_43 a_44  ]
# ordering the data like  (columns first)
#       [ *      a_12   a_23    a_34
#         a_11   a_22   a_33    a_44
#         a_21   a_32   a_43    *
#         a_31   a_42   *       *       ]
###
type BandedMatrix{T} <: AbstractBandedMatrix{T}
    data::Matrix{T}  # l+u+1 x n (# of columns)
    m::Int #Number of rows
    l::Int # lower bandwidth ≥0
    u::Int # upper bandwidth ≥0
    function BandedMatrix(data::Matrix{T},m,l,u)
        if l < 0 || u < 0
            error("Bandwidths $l,$u must both be non-negative")
        elseif size(data,1)!=l+u+1
            error("Data matrix must have number rows equal to number of bands")
        else
            new(data,m,l,u)
        end
    end
end

# BandedSubMatrix are also banded
typealias BandedSubMatrix{T} Union{
                SubArray{T,2,BandedMatrix{T},Tuple{UnitRange{Int},UnitRange{Int}}},
                SubArray{T,2,BandedMatrix{T},Tuple{Colon,UnitRange{Int}}},
                SubArray{T,2,BandedMatrix{T},Tuple{Colon,Colon}}
            }

include("BandedLU.jl")
include("BandedQR.jl")
include("blas.jl")

## Constructors

BandedMatrix(data::Matrix,m::Integer,a::Integer,b::Integer) = BandedMatrix{eltype(data)}(data,m,a,b)

# Use zeros to avoid unallocated entries for bigfloat
BandedMatrix{T<:BlasFloat}(::Type{T},n::Integer,m::Integer,a::Integer,b::Integer) =
    BandedMatrix{T}(Array(T,b+a+1,m),n,a,b)
BandedMatrix{T<:Number}(::Type{T},n::Integer,m::Integer,a::Integer,b::Integer) =
    BandedMatrix{T}(zeros(T,b+a+1,m),n,a,b)
BandedMatrix{T}(::Type{T},n::Integer,m::Integer,a::Integer,b::Integer) =
    BandedMatrix{T}(Array(T,b+a+1,m),n,a,b)

BandedMatrix{T}(::Type{T},n::Integer,a::Integer,b::Integer) = BandedMatrix(T,n,n,a,b)
BandedMatrix{T}(::Type{T},n::Integer,::Colon,a::Integer,b::Integer) = BandedMatrix(T,n,n+b,a,b)

BandedMatrix(data::Matrix,m::Integer,a) = BandedMatrix(data,m,-a[1],a[end])
BandedMatrix{T}(::Type{T},n::Integer,m::Integer,a) = BandedMatrix(T,n,m,-a[1],a[end])
BandedMatrix{T}(::Type{T},n::Integer,::Colon,a) = BandedMatrix(T,n,:,-a[1],a[end])
BandedMatrix{T}(::Type{T},n::Integer,a) = BandedMatrix(T,n,-a[1],a[end])

Base.convert{V}(::Type{BandedMatrix{V}},M::BandedMatrix) = BandedMatrix{V}(convert(Matrix{V},M.data),M.m,M.l,M.u)
function Base.convert{BM<:BandedMatrix}(::Type{BM},M::Matrix)
    ret=BandedMatrix(eltype(BM)==Any?eltype(M):promote_type(eltype(BM),eltype(M)),size(M,1),size(M,2),size(M,1)-1,size(M,2)-1)
    for k=1:size(M,1),j=1:size(M,2)
        ret[k,j]=M[k,j]
    end
    ret
end

Base.promote_rule{T,V}(::Type{BandedMatrix{T}},::Type{BandedMatrix{V}})=BandedMatrix{promote_type(T,V)}

for (op,bop) in ((:(Base.rand),:brand),(:(Base.zeros),:bzeros),(:(Base.ones),:bones))
    @eval begin
        $bop{T}(::Type{T},n::Integer,m::Integer,a::Integer,b::Integer)=BandedMatrix($op(T,b+a+1,m),n,a,b)
        $bop{T}(::Type{T},n::Integer,a::Integer,b::Integer)=$bop(T,n,n,a,b)
        $bop{T}(::Type{T},n::Integer,::Colon,a::Integer,b::Integer)=$bop(T,n,n+b,a,b)
        $bop{T}(::Type{T},::Colon,m::Integer,a::Integer,b::Integer)=$bop(T,m+a,m,a,b)
        $bop(n::Integer,m::Integer,a::Integer,b::Integer)=$bop(Float64,n,m,a,b)
        $bop(n::Integer,a::Integer,b::Integer)=$bop(n,n,a,b)

        $bop{T}(::Type{T},n::Integer,m::Integer,a)=$bop(T,n,m,-a[1],a[end])
        $bop{T}(::Type{T},n::Number,::Colon,a)=$bop(T,n,:,-a[1],a[end])
        $bop{T}(::Type{T},::Colon,m::Integer,a)=$bop(T,:,m,-a[1],a[end])
        $bop{T}(::Type{T},n::Integer,a)=$bop(T,n,-a[1],a[end])
        $bop(n::Integer,m::Integer,a)=$bop(Float64,n,m,-a[1],a[end])
        $bop(n::Integer,a)=$bop(n,-a[1],a[end])

        $bop(B::AbstractBandedMatrix)=$bop(eltype(B),size(B,1),size(B,2),bandwidth(B,1),bandwidth(B,2))
    end
end

Base.similar(B::AbstractBandedMatrix)=BandedMatrix(eltype(B),size(B,1),size(B,2),bandwidth(B,1),bandwidth(B,2))

function beye{T}(::Type{T},n::Integer,a...)
    ret=bzeros(T,n,a...)
    for k=1:n
         ret[k,k]=one(T)
    end
    ret
end
beye{T}(::Type{T},n::Integer) = beye(T,n,0,0)
beye(n::Integer) = beye(n,0,0)
beye(n::Integer,a...) = beye(Float64,n,a...)

## Abstract Array Interface

Base.size(A::BandedMatrix,k) = ifelse(k==2,size(A.data,2),A.m)
Base.size(A::BandedMatrix) = A.m,size(A.data,2)

Base.linearindexing{T}(::Type{BandedMatrix{T}})=Base.LinearSlow()

unsafe_getindex(A::BandedMatrix,k::Integer,j::Integer)=A.data[k-j+A.u+1,j]

function getindex(A::BandedMatrix,kr::UnitRange,jr::UnitRange)
    shft=first(kr)-first(jr)
    if A.l-shft ≥ 0 && A.u+shft ≥ 0
        BandedMatrix(A.data[:,jr],length(kr),A.l-shft,A.u+shft)
    else
        #TODO: Make faster
        ret=bzeros(eltype(A),length(kr),length(jr),max(0,A.l-shft),max(0,A.u+shft))
        for (k,j) in eachbandedindex(ret)
            ret[k,j]=A[kr[k],jr[j]]
        end
        ret
    end
end

getindex(A::BandedMatrix,::Colon,jr::UnitRange)=A[1:size(A,1),jr]
getindex(A::BandedMatrix,kr::UnitRange,jr::Colon)=A[kr,1:size(A,2)]
getindex(A::BandedMatrix,kr::Colon,jr::Colon)=copy(A)

# Additional get index overrides
unsafe_getindex(A::BandedMatrix,kr::Range,jr::Integer)=vec(A.data[kr-j+A.u+1,j])

getindex(A::BandedMatrix,kr::Range,j::Integer)=-A.l≤j-kr[1]≤j-kr[end]≤A.u?unsafe_getindex(A,kr,j):[A[k,j] for k=kr]

function getindex(A::BandedMatrix,kr::UnitRange,j::Integer)
    if -A.l≤j-kr[1]≤j-kr[end]≤A.u
        unsafe_getindex(A,kr,j)
    else
        eltype(A)[A[k,j] for k=kr]
    end
end

getindex(A::BandedMatrix,kr::Range,jr::Range)=[A[k,j] for k=kr,j=jr]

# ~~ setindex! ~~

# TODO
# ~ implement indexing with vectors of indices
# ~ implement scalar/vector - band - integer
# ~ implement scalar/vector - band - range

# ~ Utilities ~

# prepare for v0.5 new bound checking facilities
if VERSION < v"0.5"
    macro boundscheck(ex)
        ex
    end
end

# ~~ Type to set\get data along a band
immutable Band
    i::Int
end
band(i::Int) = Band(i)

# ~~ Indexing on the i-th row/column within band range
immutable BandRange end

# ~~ Out of band error
immutable BandError <: Exception
    A::BandedMatrix
    i::Int
end

function showerror(io::IO, e::BandError)
    A, i, u, l = e.A, e.i, e.A.u, e.A.l
    print(io, "attempt to access $(typeof(A)) with bandwidths " *
              "($(-l), $u) at band $i")
end

# start/stop indices of the i-th column/row, bounded by actual matrix size
@inline colstart(A, i::Integer) = min(max(i-bandwidth(A,2), 1), size(A, 2))
@inline  colstop(A, i::Integer) = min(i+bandwidth(A,1), size(A, 1))
@inline rowstart(A, i::Integer) = min(max(i-bandwidth(A,1), 1), size(A, 1))
@inline  rowstop(A, i::Integer) = min(i+bandwidth(A,2), size(A, 2))

@inline colstart(A::BandedMatrix, i::Integer) = min(max(i-A.u, 1), size(A, 2))
@inline  colstop(A::BandedMatrix, i::Integer) = min(i+A.l, size(A, 1))
@inline rowstart(A::BandedMatrix, i::Integer) = min(max(i-A.l, 1), size(A, 1))
@inline  rowstop(A::BandedMatrix, i::Integer) = min(i+A.u, size(A, 2))

@inline colrange(A, i::Integer) = colstart(A,i):colstop(A,i)
@inline rowrange(A, i::Integer) = rowstart(A,i):rowstop(A,i)

# length of i-the column/row
@inline collength(A, i::Integer) = max(colstop(A, i) - colstart(A, i) + 1, 0)
@inline rowlength(A, i::Integer) = max(rowstop(A, i) - rowstart(A, i) + 1, 0)

# length of diagonal
@inline diaglength(A::BandedMatrix, b::Band) = diaglength(A, b.i)
@inline function diaglength(A::BandedMatrix, i::Integer)
    min(size(A, 2), size(A, 1)+i) - max(0, i)
end

# return id of first empty diagonal intersected along row k
function _firstdiagrow(A, k)
    a, b = rowstart(A, k), rowstop(A, k)
    c = a == 1 ? b+1 : a-1
    c-k
end

# return id of first empty diagonal intersected along column j
function _firstdiagcol(A, j)
    a, b = colstart(A, j), colstop(A, j)
    r = a == 1 ? b+1 : a-1
    j-r
end

# ~ bound checking functions ~

checkbounds(A::BandedMatrix, k::Integer, j::Integer) =
    (0 < k ≤ size(A, 1) && 0 < j ≤ size(A, 2) || throw(BoundsError(A, (k,j))))

checkbounds(A::BandedMatrix, kr::Range, j::Integer) =
    (checkbounds(A, first(kr), j); checkbounds(A,  last(kr), j))

checkbounds(A::BandedMatrix, k::Integer, jr::Range) =
    (checkbounds(A, k, first(jr)); checkbounds(A, k,  last(jr)))

checkbounds(A::BandedMatrix, kr::Range, jr::Range) =
    (checkbounds(A, kr, first(jr)); checkbounds(A, kr,  last(jr)))

# check indices fall in the band
checkband(A::BandedMatrix, i::Integer) =
    (bandinds(A, 1) ≤ i ≤ bandinds(A, 2) || throw(BandError(A, i)))

checkband(A::BandedMatrix, b::Band) = checkband(A, b.i)

checkband(A::BandedMatrix, k::Integer, j::Integer) = checkband(A, j-k)

checkband(A::BandedMatrix, kr::Range, j::Integer) =
    (checkband(A, first(kr), j); checkband(A,  last(kr), j))

checkband(A::BandedMatrix, k::Integer, jr::Range) =
    (checkband(A, k, first(jr)); checkband(A, k,  last(jr)))

checkband(A::BandedMatrix, kr::Range, jr::Range) =
    (checkband(A, kr, first(jr)); checkband(A, kr,  last(jr)))

# check dimensions of inputs
checkdimensions(sizedest::Tuple{Int, Vararg{Int}}, sizesrc::Tuple{Int, Vararg{Int}}) =
    (sizedest == sizesrc ||
        throw(DimensionMismatch("tried to assign $(sizesrc) sized " *
                                "array to $(sizedest) destination")) )

checkdimensions(dest::AbstractVector, src::AbstractVector) =
    checkdimensions(size(dest), size(src))

checkdimensions(ldest::Int, src::AbstractVector) =
    checkdimensions((ldest, ), size(src))

checkdimensions(kr::Range, jr::Range, src::AbstractMatrix) =
    checkdimensions((length(kr), length(jr)), size(src))

# ~ Special setindex methods ~

# slow fall back method
@inline unsafe_setindex!(A::BandedMatrix, v, k::Integer, j::Integer) =
    unsafe_setindex!(A.data, A.u, v, k, j)

# fast method used below
@inline unsafe_setindex!{T}(data::Matrix{T}, u::Integer, v, k::Integer, j::Integer) =
    @inbounds data[u + k - j + 1, j] = convert(T, v)::T

# scalar - integer - integer
function setindex!(A::BandedMatrix, v, k::Integer, j::Integer)
    @boundscheck  checkbounds(A, k, j)
    if bandinds(A, 1) ≤ j-k ≤ bandinds(A, 2)
        unsafe_setindex!(A.data, A.u, v, k, j)
    elseif v ≠ 0  # allow setting outside bands to zero
        throw(BandError(A,j-k))
    end
end

# scalar - colon - colon
setindex!{T}(A::BandedMatrix{T}, v, ::Colon, ::Colon) =
    @inbounds A.data[:] = convert(T, v)::T

# scalar - colon
setindex!{T}(A::BandedMatrix{T}, v, ::Colon) =
    @inbounds A.data[:] = convert(T, v)::T

# ~ indexing along a band

# scalar - band - colon
function setindex!{T}(A::BandedMatrix{T}, v, b::Band)
    @boundscheck checkband(A, b)
    @inbounds A.data[A.u - b.i + 1, :] = convert(T, v)::T
end

# vector - band - colon
function setindex!{T}(A::BandedMatrix{T}, V::AbstractVector, b::Band)
    @boundscheck checkband(A, b)
    @boundscheck checkdimensions(diaglength(A, b), V)
    row = A.u - b.i + 1
    data, i = A.data, max(b.i - A.u + 2, 1)
    for v in V
        @inbounds data[row, i] = convert(T, v)::T
        i += 1
    end
    V
end

# ~ indexing along columns

# scalar - colon - integer -- A[:, 1] = 2 - not allowed
setindex!(A::BandedMatrix, v, kr::Colon, j::Integer) =
    throw(BandError(A, _firstdiagcol(A, j)))

# vector - colon - integer -- A[:, 1] = [1, 2, 3] - not allowed
setindex!(A::BandedMatrix, V::AbstractVector, kr::Colon, j::Integer) =
    throw(BandError(A, _firstdiagcol(A, j)))

# scalar - BandRange - integer -- A[1, BandRange] = 2
setindex!{T}(A::BandedMatrix{T}, v, ::Type{BandRange}, j::Integer) =
    (A[colstart(A, j):colstop(A, j), j] = convert(T, v)::T) # call range method

# vector - BandRange - integer -- A[1, BandRange] = 2
setindex!(A::BandedMatrix, V::AbstractVector, ::Type{BandRange}, j::Integer) =
    (A[colstart(A, j):colstop(A, j), j] = V) # call range method

# scalar - range - integer -- A[1:2, 1] = 2
function setindex!(A::BandedMatrix, v, kr::Range, j::Integer)
    @boundscheck checkbounds(A, kr, j)
    @boundscheck  checkband(A, kr, j)
    data, u = A.data, A.u
    for k in kr
        unsafe_setindex!(data, u, v, k, j)
    end
    v
end

# vector - range - integer -- A[1:3, 1] = [1, 2, 3]
function setindex!(A::BandedMatrix, V::AbstractVector, kr::Range, j::Integer)
    @boundscheck checkbounds(A, kr, j)
    @boundscheck checkdimensions(kr, V)
    @boundscheck  checkband(A, kr, j)
    data, u, i = A.data, A.u, 0
    for v in V
        unsafe_setindex!(data, u, v, kr[i+=1], j)
    end
    V
end

# ~ indexing along a row

# scalar - integer - colon -- A[1, :] = 2 - not allowed
setindex!(A::BandedMatrix, v, k::Integer, jr::Colon) =
    throw(BandError(A, _firstdiagrow(A, k)))

# vector - integer - colon -- A[1, :] = [1, 2, 3] - not allowed
setindex!(A::BandedMatrix, V::AbstractVector, k::Integer, ::Colon) =
    throw(BandError(A, _firstdiagrow(A, k)))

# scalar - integer - BandRange -- A[1, BandRange] = 2
setindex!{T}(A::BandedMatrix{T}, v, k::Integer, ::Type{BandRange}) =
    (A[k, rowstart(A, k):rowstop(A, k)] = convert(T, v)::T) # call range method

# vector - integer - BandRange -- A[1, BandRange] = [1, 2, 3]
setindex!(A::BandedMatrix, V::AbstractVector, k::Integer, ::Type{BandRange}) =
    (A[k, rowstart(A, k):rowstop(A, k)] = V) # call range method

# scalar - integer - range -- A[1, 2:3] = 3
function setindex!(A::BandedMatrix, v, k::Integer, jr::Range)
    @boundscheck checkbounds(A, k, jr)
    @boundscheck checkband(A, k, jr)
    data, u = A.data, A.u
    for j in jr
        unsafe_setindex!(data, u, v, k, j)
    end
    v
end

# vector - integer - range -- A[1, 2:3] = [3, 4]
function setindex!(A::BandedMatrix, V::AbstractVector, k::Integer, jr::Range)
    @boundscheck checkbounds(A, k, jr)
    @boundscheck checkband(A, k, jr)
    @boundscheck checkdimensions(jr, V)
    data, u, i = A.data, A.u, 0
    for v in V
        unsafe_setindex!(data, u, v, k, jr[i+=1])
    end
    V
end

# ~ indexing over a rectangular block

# scalar - range - range
function setindex!(A::BandedMatrix, v, kr::Range, jr::Range)
    @boundscheck checkbounds(A, kr, jr)
    @boundscheck checkband(A, kr, jr)
    data, u = A.data, A.u
    for j in jr, k in kr
        unsafe_setindex!(data, u, v, k, j)
    end
    v
end

# matrix - range - range
function setindex!(A::BandedMatrix, V::AbstractMatrix, kr::Range, jr::Range)
    @boundscheck checkbounds(A, kr, jr)
    @boundscheck checkband(A, kr, jr)
    @boundscheck checkdimensions(kr, jr, V)
    data, u = A.data, A.u
    jj = 1
    for j in jr
        kk = 1
        for k in kr
            # we index V manually in column-major order
            @inbounds unsafe_setindex!(data, u, V[kk, jj], k, j)
            kk += 1
        end
        jj += 1
    end
    V
end

# ~~ end setindex! ~~

function Base.convert(::Type{Matrix},A::BandedMatrix)
    ret=zeros(eltype(A),size(A,1),size(A,2))
    for (k,j) in eachbandedindex(A)
        @inbounds ret[k,j]=A[k,j]
    end
    ret
end

Base.full(A::BandedMatrix) = convert(Matrix,A)

## Band range

bandwidth(A::BandedMatrix,k::Integer) = k==1?A.l:A.u

function Base.sparse(B::BandedMatrix)
    i=Array(Int,length(B.data));j=Array(Int,length(B.data))
    n,m=size(B.data)
    Bn=size(B,1)
    vb=copy(vec(B.data))
    for κ=1:n,ℓ=1:m
        j[κ+n*(ℓ-1)]=ℓ
        ii=κ+ℓ-B.u-1
        if ii <1 || ii > Bn
            vb[κ+n*(ℓ-1)] = 0
        end
        i[κ+n*(ℓ-1)]=min(max(ii,1),Bn)
    end

    sparse(i,j,vb,Bn,m)
end

# pass standard routines to full matrix

Base.norm(B::BandedMatrix,opts...)=norm(full(B),opts...)

# We turn off bound checking to allow nicer syntax without branching
#setindex!(A::BandedMatrix,v,k::Integer,j::Integer)=((A.l≤j-k≤A.u)&&k≤A.n)?ussetindex!(A,v,k,j):throw(BoundsError())
#setindex!(A::BandedMatrix,v,kr::Range,j::Integer)=(A.l≤j-kr[end]≤j-kr[1]≤A.u&&kr[end]≤A.n)?ussetindex!(A,v,kr,j):throw(BoundsError())

## ALgebra and other functions

function Base.maximum(B::BandedMatrix)
    m=zero(eltype(B))
    for (k,j) in eachbandedindex(B)
        m=max(B[k,j],m)
    end
    m
end

for OP in (:*,:.*,:+,:.+,:-,:.-)
    @eval begin
        $OP(B::BandedMatrix{Bool},x::Bool) = BandedMatrix($OP(B.data,x),B.m,B.l,B.u)
        $OP(x::Bool,B::BandedMatrix{Bool}) = BandedMatrix($OP(x,B.data),B.m,B.l,B.u)
        $OP(B::BandedMatrix,x::Number) = BandedMatrix($OP(B.data,x),B.m,B.l,B.u)
        $OP(x::Number,B::BandedMatrix) = BandedMatrix($OP(x,B.data),B.m,B.l,B.u)
    end
end

function +{T,V}(A::BandedMatrix{T},B::BandedMatrix{V})
    if size(A) != size(B)
        throw(DimensionMismatch("+"))
    end
    n,m=size(A,1),size(A,2)

    ret = bzeros(promote_type(T,V),n,m,max(A.l,B.l),max(A.u,B.u))
    BLAS.axpy!(1.,A,ret)
    BLAS.axpy!(1.,B,ret)

    ret
end

function -{T,V}(A::BandedMatrix{T},B::BandedMatrix{V})
    if size(A) != size(B)
        throw(DimensionMismatch("+"))
    end
    n,m=size(A,1),size(A,2)

    ret = bzeros(promote_type(T,V),n,m,max(A.l,B.l),max(A.u,B.u))
    BLAS.axpy!(1.,A,ret)
    BLAS.axpy!(-1.,B,ret)

    ret
end

function *(A::BandedMatrix,B::BandedMatrix)
    if size(A,2)!=size(B,1)
        throw(DimensionMismatch("*"))
    end
    n,m=size(A,1),size(B,2)

    ret=BandedMatrix(promote_type(eltype(A),eltype(B)),n,m,A.l+B.l,A.u+B.u)
    for (k,j) in eachbandedindex(ret)
        νmin=max(1,k-bandwidth(A,1),j-bandwidth(B,2))
        νmax=min(size(A,2),k+bandwidth(A,2),j+bandwidth(B,1))

        ret[k,j]=A[k,νmin]*B[νmin,j]
        for ν=νmin+1:νmax
            ret[k,j]+=A[k,ν]*B[ν,j]
        end
    end

    ret
end

function *{T<:Number,V<:Number}(A::BandedMatrix{T},B::BandedMatrix{V})
    if size(A,2)!=size(B,1)
        throw(DimensionMismatch("*"))
    end
    n,m=size(A,1),size(B,2)

    A_mul_B!(bzeros(promote_type(T,V),n,m,A.l+B.l,A.u+B.u),A,B)
end

function *{T<:Number,V<:Number}(A::BandedMatrix{T},B::Matrix{V})
    if size(A,2)!=size(B,1)
        throw(DimensionMismatch("*"))
    end
    n,m=size(A,1),size(B,2)

    A_mul_B!(Array(promote_type(T,V),n,m),A,B)
end

*{T<:Number,V<:Number}(A::Matrix{T},B::BandedMatrix{V}) =
    A*full(B)

*{T<:BlasFloat}(A::BandedMatrix{T},b::Vector{T}) =
    BLAS.gbmv('N',A.m,A.l,A.u,one(T),A.data,b)

function *{T}(A::BandedMatrix{T},b::Vector{T})
    ret = zeros(T,size(A,1))
    for (k,j) in eachbandedindex(A)
        @inbounds ret[k]+=A[k,j]*b[j]
    end
    ret
end

function *(A::BandedMatrix,b::Vector)
    T=promote_type(eltype(A),eltype(b))
    convert(BandedMatrix{T},A)*convert(Vector{T},b)
end

function Base.transpose(B::BandedMatrix)
    Bt=bzeros(eltype(B),size(B,2),size(B,1),B.u,B.l)
    for (k,j) in eachbandedindex(B)
       Bt[j,k]=B[k,j]
    end
    Bt
end

function Base.ctranspose(B::BandedMatrix)
    Bt=bzeros(eltype(B),size(B,2),size(B,1),B.u,B.l)
    for (k,j) in eachbandedindex(B)
       Bt[j,k]=conj(B[k,j])
    end
    Bt
end

function Base.diag{T}(A::BandedMatrix{T})
    n=size(A,1)
    @assert n==size(A,2)

    vec(A.data[A.u+1,1:n])
end

# basic interface
(\){T<:BlasFloat}(A::Union{BandedLU{T}, BandedMatrix{T}}, B::StridedVecOrMat{T}) =
    A_ldiv_B!(A, copy(B)) # makes a copy

## Matrix.*Matrix

function .*(A::BandedMatrix,B::BandedMatrix)
    @assert size(A,1)==size(B,1)&&size(A,2)==size(B,2)

    l=min(A.l,B.l);u=min(A.u,B.u)
    T=promote_type(eltype(A),eltype(B))
    ret=BandedMatrix(T,size(A,1),size(A,2),l,u)

    for (k,j) in eachbandedindex(ret)
        @inbounds ret[k,j]=A[k,j]*B[k,j]
    end
    ret
end

#implements fliplr(flipud(A))
function fliplrud(A::BandedMatrix)
    n,m=size(A)
    l=A.u+n-m
    u=A.l+m-n
    ret=BandedMatrix(eltype(A),n,m,l,u)
    for (k,j) in eachbandedindex(ret)
        @inbounds ret[k,j]=A[n-k+1,m-j+1]
    end
    ret
end

for OP in (:(Base.real),:(Base.imag))
    @eval $OP(A::BandedMatrix) =
        BandedMatrix($OP(A.data),A.m,A.l,A.u)
end

## Show

type PrintShow
    str
end
Base.show(io::IO,N::PrintShow) = print(io,N.str)

if VERSION < v"0.5.0-rc4"
    showarray(io,M;opts...) = Base.showarray(io,M;opts...)
    function Base.showarray(io::IO,B::AbstractBandedMatrix;
                   header::Bool=true, limit::Bool=Base._limit_output,
                   sz = (s = Base.tty_size(); (s[1]-4, s[2])), repr=false)
       header && print(io,summary(B))

       if !isempty(B) && size(B,1) ≤ 1000 && size(B,2) ≤ 1000
           header && println(io,":")
           M=Array(Any,size(B)...)
           fill!(M,PrintShow(""))
           for (k,j) in eachbandedindex(B)
               M[k,j]=B[k,j]
           end

           showarray(io,M;header=false)
       end
   end

else
    showarray(io,M;opts...) = Base.showarray(io,M,false;opts...)
    function Base.showarray(io::IO,B::AbstractBandedMatrix,repr::Bool = true; header = true)
        header && print(io,summary(B))

        if !isempty(B) && size(B,1) ≤ 1000 && size(B,2) ≤ 1000
            header && println(io,":")
            M=Array(Any,size(B)...)
            fill!(M,PrintShow(""))
            for (k,j) in eachbandedindex(B)
                M[k,j]=B[k,j]
            end

            showarray(io,M;header=false)
        end
    end
end

## SubArray routines

bandshift(S) = parentindexes(S)[1][1]-parentindexes(S)[2][1]

bandwidth{T}(S::BandedSubMatrix{T},k::Integer) = bandwidth(parent(S),k) + (k==1?-1:1)*bandshift(S)

unsafe_getindex{T}(S::BandedSubMatrix{T},k,j) = unsafe_getindex(parent(S),parentindexes(S)[1][k],parentindexes(S)[2][j])

end #module


text:
?println

println("¡Estoy emocionado por aprender Julia!")

123
456

my_answer = 42
typeof(my_answer)

my_pi = 3.14159
typeof(my_pi)

my_name = "Jane"
typeof(my_name)

my_answer = my_name

typeof(my_answer)

# Puedes dejar un comentario en una sola línea usando la tecla de gato

#=

Para comentarios de varias lineas,
usa la secuencia de '#= =#' .

=#

suma = 3 + 7

resta = 10 - 3

producto = 20 * 5

cociente = 100 / 10

potencia = 10 ^ 2

módulo = 101 % 2

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
#=
lm-aa
2019-08-07 Jeff Fessler, University of Michigan
=#

export lmaa_hcat, lmaa_vcat, lmaa_hvcat

using LinearMaps: LinearMap
using LinearAlgebra: UniformScaling, I
import LinearAlgebra: issymmetric, ishermitian, isposdef
import LinearAlgebra: mul!, lmul!, rmul!
import SparseArrays: sparse

Indexer = AbstractVector{Int}

"""
`mutable struct LinearMapAA{T} <: AbstractMatrix{T}`

old way may not properly allow `setindex!` to work as desired
because it may change the type of the lmap and of the prop:
`struct LinearMapAA{T, M <: LinearMap, P <: NamedTuple} <: AbstractMatrix{T}`
"""
mutable struct LinearMapAA{T} <: AbstractMatrix{T}
#{T, M <: LinearMap, P <: NamedTuple}
#	_lmap::M
	_lmap::LinearMap
#	_prop::P
	_prop::NamedTuple
#=
	function LinearMapAA{T}(L::M, p::P) where {T, M <: LinearMap, P <: NamedTuple}
	function LinearMapAA(L::LinearMap, p::NamedTuple) # where {T, M <: LinearMap, P <: NamedTuple}
	#	new{T,M,P}(L, p)
		new(L, p)
	end
=#
end

include("setindex.jl")

# constructors

"""
`A = LinearMapAA(L::LinearMap [, prop::NamedTuple ; T = eltype(L)])`
constructor

`prop` cannot include the field `_lmap` or `_prop`
"""
function LinearMapAA(L::LinearMap, prop::NamedTuple ; T = eltype(L))
	:_lmap in propertynames(prop) && throw("cannot use _lmap")
	:_prop in propertynames(prop) && throw("cannot use _prop")
	return LinearMapAA{T}(L, prop)
end
LinearMapAA(L::LinearMap) = LinearMapAA(L, (none=nothing,))

"""
`A = LinearMapAA(L::AbstractMatrix [, prop::NamedTuple])`
constructor
"""
LinearMapAA(L::AbstractMatrix, prop::NamedTuple) =
	LinearMapAA(LinearMap(L), prop)
LinearMapAA(L::AbstractMatrix) = LinearMapAA(L, (none=nothing,))

"""
`A = LinearMapAA(f::Function, fc::Function, D::Dims{2} [, prop::NamedTuple)] ; T::DataType`
constructor
"""
LinearMapAA(f::Function, fc::Function, D::Dims{2}, prop::NamedTuple ;
	T::DataType = Float32) =
	LinearMapAA(LinearMap{T}(f, fc, D[1], D[2]), prop)
LinearMapAA(f::Function, fc::Function, D::Dims{2}, T::DataType = Float32) =
	LinearMapAA(f, fc, D, (none=nothing,) ; T=T)

"""
`A = LinearMapAA(f::Function, D::Dims{2} [, prop::NamedTuple)]`
constructor
"""
LinearMapAA(f::Function, D::Dims{2}, prop::NamedTuple ; T::DataType = Float32) =
	LinearMapAA(LinearMap{T}(f, D[1], D[2]), prop)
LinearMapAA(f::Function, D::Dims{2} ; T::DataType = Float32) =
	LinearMapAA(f, D, (none=nothing,), T=T)

# copy
Base.copy(A::LinearMapAA) = LinearMapAA(A._lmap, A._prop)

# Matrix
Base.Matrix(A::LinearMapAA) = Matrix(A._lmap)

# ndims
# Base.ndims(A::LinearMapAA) = ndims(A._lmap) # 2 for AbstractMatrix

# display
Base.display(A::LinearMapAA) =
	begin
		println("LinearMapAA: $(size(A,1))×$(size(A,2)) $(A._prop)")
		tmp = "$(A._lmap)"[1:77]
		println(" $tmp ..")
	#	display(A._lmap)
	#	display(A._prop)
	end

# size
Base.size(A::LinearMapAA) = size(A._lmap)
Base.size(A::LinearMapAA, d::Int) = size(A._lmap, d)

# adjoint
Base.adjoint(A::LinearMapAA) = LinearMapAA(A._lmap', A._prop)

# transpose
Base.transpose(A::LinearMapAA) = LinearMapAA(transpose(A._lmap), A._prop)

# eltype
Base.eltype(A::LinearMapAA) = eltype(A._lmap)

# LinearMap algebraic properties
issymmetric(A::LinearMapAA) = issymmetric(A._lmap)
#ishermitian(A::LinearMapAA{<:Real}) = issymmetric(A._lmap)
ishermitian(A::LinearMapAA) = ishermitian(A._lmap)
isposdef(A::LinearMapAA) = isposdef(A._lmap)

# comparison of LinearMapAA objects, sufficient but not necessary
Base.:(==)(A::LinearMapAA, B::LinearMapAA) =
	(eltype(A) == eltype(B) && A._lmap == B._lmap && A._prop == B._prop)

# convert to sparse
sparse(A::LinearMapAA) = sparse(A._lmap)

# cat (hcat, vcat, hvcat) are tricky for avoiding type piracy
# It is especially hard to handle AbstractMatrix,
# so typically force the user to wrap it in LinearMap(AA) first.
#LMcat{T} = Union{LinearMapAA{T}, LinearMap{T}, UniformScaling{T},AbstractMatrix{T}}
LMcat{T} = Union{LinearMapAA{T}, LinearMap{T}, UniformScaling{T}} # settle
#LMelse{T} = Union{LinearMap{T},UniformScaling{T},AbstractMatrix{T}} # non AA

# convert to something suitable for LinearMap.*cat
function lm_promote(A::LMcat)
	A isa LinearMapAA ? A._lmap :
	A isa AbstractMatrix ? LinearMap(A) :
	A isa UniformScaling ? A : # leave unchanged - ok for LinearMaps.*cat
	# A isa LinearMap ?
	A # otherwise it is this
end

# single-letter codes for cat objects, e.g., [A I A] becomes "AIA"
#= not so useful
function lm_code(A)
	isa(A, LinearMapAA) ? "A" :
	isa(A, AbstractMatrix) ? "M" :
	isa(A, UniformScaling) ? "I" :
	isa(A, LinearMap) ? "L" :
	"?"
end
=#

# concatenate the single-letter codes, e.g., [A I A] becomes "AIA"
# lm_name = As -> *(lm_code.(As)...)
lm_name = As -> nothing

# these rely on LinearMap.*cat methods
"`B = lmaa_hcat(A1, A2, ...)` `hcat` of multiple objects"
lmaa_hcat(As::LMcat...) =
	LinearMapAA(hcat(lm_promote.(As)...), (hcat=lm_name(As),))
"`B = lmaa_vcat(A1, A2, ...)` `vcat` of multiple objects"
lmaa_vcat(As::LMcat...) =
	LinearMapAA(vcat(lm_promote.(As)...), (vcat=lm_name(As),))
"`B = lmaa_hvcat(rows, A1, A2, ...)` `hvcat` of multiple objects"
lmaa_hvcat(rows::NTuple{nr,Int} where {nr}, As::LMcat...) =
	LinearMapAA(hvcat(rows, lm_promote.(As)...), (hvcat=lm_name(As),))

# a single leading LinearMapAA followed by others is clear
Base.hcat(A1::LinearMapAA, As::LMcat...) =
	lmaa_hcat(A1, As...)
Base.vcat(A1::LinearMapAA, As::LMcat...) =
	lmaa_vcat(A1, As...)
Base.hvcat(rows::NTuple{nr,Int} where {nr}, A1::LinearMapAA, As::LMcat...) =
	lmaa_hvcat(rows, A1, As...)
# or in 2nd position, beyond that, user can use lmaa_*
#Base.hcat(A1::LMelse, A2::LinearMapAA, As::LMcat...) = # fails!?
Base.hcat(A1::UniformScaling, A2::LinearMapAA, As::LMcat...) =
	lmaa_hcat(A1, A2, As...)
Base.vcat(A1::UniformScaling, A2::LinearMapAA, As::LMcat...) =
	lmaa_vcat(A1, A2, As...)
Base.hvcat(rows::NTuple{nr,Int} where nr,
		A1::UniformScaling, A2::LinearMapAA, As::LMcat...) =
	lmaa_hvcat(rows, A1, A2, As...)

# multiply with vectors

mul!(y::AbstractVector, A::LinearMapAA, x::AbstractVector) = mul!(y, A._lmap, x)

#= these seem pointless; see multiplication with scalars below
lmul!(s::Number, A::LinearMapAA) = lmul!(s, A._lmap)
rmul!(A::LinearMapAA, s::Number) = rmul!(A._lmap, s)
=#

#=
function A_mul_B!(y::AbstractVector, A::LinearMapAA, x::AbstractVector)
	A_mul_B!(y, A._lmap, x)
	return y
end

function At_mul_B!(x::AbstractVector, A::LinearMapAA, y::AbstractVector)
	At_mul_B!(x, A._lmap, y)
	return x
end

function Ac_mul_B!(x::AbstractVector, A::LinearMapAA, y::AbstractVector)
	Ac_mul_B!(x, A._lmap, y)
	return x
end
=#

# add or subtract objects
Base.:(+)(A::LinearMapAA, B::LinearMapAA) =
	LinearMapAA(A._lmap + B._lmap, (sum=nothing,))
Base.:(+)(A::LinearMapAA, B::AbstractMatrix) =
	LinearMapAA(A._lmap + LinearMap(B), A._prop)
Base.:(+)(A::AbstractMatrix, B::LinearMapAA) =
	LinearMapAA(LinearMap(A) + B._lmap, B._prop)
Base.:(-)(A::LinearMapAA, B::LinearMapAA) = A + (-1)*B
Base.:(-)(A::LinearMapAA, B::AbstractMatrix) = A + (-1)*B
Base.:(-)(A::AbstractMatrix, B::LinearMapAA) = A + (-1)*B

# multiply objects
Base.:(*)(A::LinearMapAA, B::LinearMapAA) =
	LinearMapAA(A._lmap * B._lmap, (prod=nothing,))
Base.:(*)(A::LinearMapAA, B::AbstractMatrix) =
	LinearMapAA(A._lmap * LinearMap(B), A._prop)
Base.:(*)(A::AbstractMatrix, B::LinearMapAA) =
	LinearMapAA(LinearMap(A) * B._lmap, B._prop)

# multiply with I or s*I (identity or scaled identity)
Base.:(*)(A::LinearMapAA, B::UniformScaling{Bool}) = B.λ ? A : (A * B.λ)
Base.:(*)(A::LinearMapAA, B::UniformScaling) = (B.λ == 1) ? A : (A * B.λ)
Base.:(*)(B::UniformScaling{Bool}, A::LinearMapAA) = B.λ ? A : (B.λ * A)
Base.:(*)(B::UniformScaling, A::LinearMapAA) = (B.λ == 1) ? A : (B.λ * A)

# multiply with vector
Base.:(*)(A::LinearMapAA, v::AbstractVector{<:Number}) = A._lmap * v

# multiply with scalars
Base.:(*)(s::Number, A::LinearMapAA) = LinearMapAA((s*I) * A._lmap, A._prop)
Base.:(*)(A::LinearMapAA, s::Number) = LinearMapAA(A._lmap * (s*I), A._prop)

# kronecker
#=
todo after LM is updated
Base.kron(A::LinearMapAA, M::AbstractMatrix) =
	LinearMapAA(kron(A._lmap, M), A._prop)
Base.kron(M::AbstractMatrix, A::LinearMapAA) =
	LinearMapAA(kron(A._lmap, M), A._prop)
Base.kron(A::LinearMapAA, B::LinearMapAA) =
	LinearMapAA(kron(A._lmap, B._lmap), (kron=nothing,))
=#

# A.?
Base.getproperty(A::LinearMapAA, s::Symbol) =
	s in (:_lmap, :_prop) ? getfield(A, s) :
#	s == :m ? size(A._lmap, 1) :
	haskey(A._prop, s) ? getfield(A._prop, s) :
		throw("unknown key $s")

Base.propertynames(A::LinearMapAA) = propertynames(A._prop)

# indexing

# [end]
function Base.lastindex(A::LinearMapAA)
	return prod(size(A._lmap))
end

# [?,end] and [end,?]
function Base.lastindex(A::LinearMapAA, d::Int)
	return size(A._lmap, d)
end

# A[i,j]
function Base.getindex(A::LinearMapAA, i::Int, j::Int)
	T = eltype(A)
	e = zeros(T, size(A._lmap,2)); e[j] = one(T)
	tmp = A._lmap * e
	return tmp[i]
end

# A[:,j]
# it is crucial to provide this function rather than to inherit from
# Base.getindex(A::AbstractArray, ::Colon, ::Int)
# because Base.getindex does this by iterating (I think).
function Base.getindex(A::LinearMapAA, ::Colon, j::Int)
	T = eltype(A)
	e = zeros(T, size(A,2)); e[j] = one(T)
	return A * e
end

# A[ii,j]
Base.getindex(A::LinearMapAA, ii::Indexer, j::Int) = A[:,j][ii]

# A[i,jj]
Base.getindex(A::LinearMapAA, i::Int, jj::Indexer) = A[i,:][jj]

# A[:,jj]
# this one is also important for efficiency
Base.getindex(A::LinearMapAA, ::Colon, jj::AbstractVector{Bool}) =
	A[:,findall(jj)]
Base.getindex(A::LinearMapAA, ::Colon, jj::Indexer) =
	hcat([A[:,j] for j in jj]...)

# A[ii,:]
# trick: cannot use A' for a FunctionMap with no fc
function Base.getindex(A::LinearMapAA, ii::Indexer, ::Colon)
	if (:fc in propertynames(A._lmap)) && isnothing(A._lmap.fc)
		return hcat([A[ii,j] for j in 1:size(A,2)]...)
	else
		return A'[:,ii]'
	end
end

# A[ii,jj]
Base.getindex(A::LinearMapAA, ii::Indexer, jj::Indexer) = A[:,jj][ii,:]

# A[k]
function Base.getindex(A::LinearMapAA, k::Int)
	c = CartesianIndices(size(A._lmap))[k] # is there a more elegant way?
	return A[c[1], c[2]]
end

# A[kk]
Base.getindex(A::LinearMapAA, kk::AbstractVector{Bool}) = A[findall(kk)]
Base.getindex(A::LinearMapAA, kk::Indexer) = [A[k] for k in kk]

# A[i,:]
# trick: one row slice returns a 1D ("column") vector
Base.getindex(A::LinearMapAA, i::Int, ::Colon) = A[[i],:][:]

# A[:,:] = Matrix(A)
Base.getindex(A::LinearMapAA, ::Colon, ::Colon) = Matrix(A._lmap)

# A[:]
Base.getindex(A::LinearMapAA, ::Colon) = A[:,:][:]

# test
using Test: @test, @test_throws

"""
`LinearMapAA_test_getindex(A::LinearMapAA)`
tests for `getindex`
"""
function LinearMapAA_test_getindex(A::LinearMapAA)
	B = Matrix(A)
	@test all(size(A) .>= (4,4)) # required by tests

	tf1 = [false; trues(size(A,1)-1)]
	tf2 = [false; trues(size(A,2)-2); false]
	ii1 = (3, 2:4, [2,4], :, tf1)
	ii2 = (2, 3:4, [1,4], :, tf2)
	for i2 in ii2
		for i1 in ii1
		#	@show i1,i2
			@test B[i1,i2] == A[i1,i2]
		end
	end

	L = A._lmap
	test_adj = !((:fc in propertynames(L)) && isnothing(L.fc))
	if test_adj
		for i1 in ii2
			for i2 in ii1
				@test B'[i1,i2] == A'[i1,i2]
			end
		end
	end

	# "end"
	@test B[3,end-1] == A[3,end-1]
	@test B[end-2,3] == A[end-2,3]
	if test_adj
		@test B'[3,end-1] == A'[3,end-1]
	end

	# [?]
	@test B[1] == A[1]
	@test B[7] == A[7]
	if test_adj
		@test B'[3] == A'[3]
	end
	@test B[end] == A[end] # lastindex

	kk = [k in [3,5] for k = 1:length(A)] # Bool
	@test B[kk] == A[kk]

	# Some tests could rely on the fact that LinearMapAA <:i AbstractMatrix,
	# by inheriting from general Base.getindex, but all are provided here.
	@test B[[1, 3, 4]] == A[[1, 3, 4]]
	@test B[4:7] == A[4:7]

	true
end

"""
`LinearMapAA_test_vmul(A::LinearMapAA)`
tests for multiply with vector and `lmul!` and `rmul!` for scalars too
"""
function LinearMapAA_test_vmul(A::LinearMapAA)
	B = Matrix(A)

	u = rand(size(A,1))
	v = rand(size(A,2))

	y = A * v
	x = A' * u
	@test isapprox(B * v, y)
	@test isapprox(B' * u, x)

	mul!(y, A, v)
	mul!(x, A', u)
	@test isapprox(B * v, y)
	@test isapprox(B' * u, x)

	s = 5.1
	C = s * A
	@test isapprox(Matrix(C), s * B)
	C = A * s
	@test isapprox(Matrix(C), B * s)

#=
	s = 5.1
	C = copy(A)
	lmul!(s, C)
	@test isapprox(s * B * v, C * v)

	C = copy(A)
	rmul!(C, s)
	@test isapprox(s * B * v, C * v)
=#

	true
end

#=
	this approach using eval() works only in the global scope!

	N = 6
	forw = x -> [cumsum(x); 0] # non-square to stress test
	back = y -> reverse(cumsum(reverse(y[1:N])))
	prop = (name="cumsum", extra=1)
	A = LinearMapAA(forw, back, (N+1, N), prop)

	list1 = [
		:([A I]), :([I A]), :([2*A 3*A]),
		:([A; I]), :([I; A]), :([2*A; 3*A]),
		:([A A I]), :([A I A]), :([2*A 3*A 4*A]),
		:([A; A; I]), :([A; I; A]), :([2*A; 3*A; 4*A]),
		:([I A I]), :([I A A]),
		:([I; A; I]), :([I; A; A]),
	#	:([I I A]), :([I; I; A]), # unsupported
		:([A A; A A]), :([A I; 2A 3I]), :([I A; 2I 3A]),
	#	:([I A; A I]), :([A I; I A]), # weird sizes
		:([A I A; 2A 3I 4A]), :([I A I; 2I 3A 4I]),
	#	:([A I A; I A A]), :([I A 2A; 3A 4I 5A]), # weird
	#	:([I I A; I A I]), # unsupported
		:([A A I; 2A 3A 4I]),
		:([A I I; 2A 3I 4I]),
	]

	M = Matrix(A)
	list2 = [
		:([M I]), :([I M]), :([2*M 3*M]),
		:([M; I]), :([I; M]), :([2*M; 3*M]),
		:([M M I]), :([M I M]), :([2*M 3*M 4*M]),
		:([M; M; I]), :([M; I; M]), :([2*M; 3*M; 4*M]),
		:([I M I]), :([I M M]),
		:([I; M; I]), :([I; M; M]),
	#	:([I I M]), :([I; I; M]), # unsupported
		:([M M; M M]), :([M I; 2M 3I]), :([I M; 2I 3M]),
	#	:([I M; M I]), :([M I; I M]), # weird sizes
		:([M I M; 2M 3I 4M]), :([I M I; 2I 3M 4I]),
	#	:([M I M; I M M]), :([I M 2M; 3M 4I 5M]), # weird
	#	:([I I M; I M I]), # unsupported
		:([M M I; 2M 3M 4I]),
		:([M I I; 2M 3I 4I]),
	]

	length(list1) != length(list2) && throw("bug")

	for ii in 1:length(list1)
		e1 = list1[ii]
		b1 = eval(e1)
		e2 = list2[ii]
		b2 = eval(e2)
		@test b1 isa LinearMapAA
	end
=#

"""
`LinearMapAA_test_cat(A::LinearMapAA)`
test hcat vcat hvcat
"""
function LinearMapAA_test_cat(A::LinearMapAA)
#	L = LinearMap(x -> A*x, y -> A'*y, size(A,1), size(A,2))
	M = Matrix(A)

#=
	# cannot get cat with AbstractMatrix to work
	M1 = reshape(1:35, N+1, N-1)
	H2 = [A M1]
	@test H2 isa LinearMapAA
	@test Matrix(H2) == [Matrix(A) H2]
	H1 = [M1 A]
	@test H1 isa LinearMapAA
	@test Matrix(H1) == [M1 Matrix(A)]

	M2 = reshape(1:(3*N), 3, N)
	V1 = [M2; A]
	@test V1 isa LinearMapAA
	@test Matrix(V1) == [M2; Matrix(A)]
	V2 = [A; M2]
	@test V2 isa LinearMapAA
	@test Matrix(V2) == [Matrix(A); M2]
=#

	fun0 = A -> [
		[A I], [I A], [2A 3A],
		[A; I], [I; A], [2A; 3A],
		[A A I], [A I A], [2A 3A 4A],
		[A; A; I], [A; I; A], [2A; 3A; 4A],
		[I A I], [I A A],
		[I; A; I], [I; A; A],
	#	[I I A], [I; I; A], # unsupported
		[A A; A A], [A I; 2A 3I], [I A; 2I 3A],
	#	[I A; A I], [A I; I A], # weird sizes
		[A I A; 2A 3I 4A], [I A I; 2I 3A 4I],
	#	[A I A; I A A], [I A 2A; 3A 4I 5A], # weird
	#	[I I A; I A I], # unsupported
		[A A I; 2A 3A 4I],
		[A I I; 2A 3I 4I],
	]

	list1 = fun0(A)
	list2 = fun0(M)

	for ii in 1:length(list1)
		b1 = list1[ii]
		b2 = list2[ii]
		@test b1 isa LinearMapAA
		@test b2 isa AbstractMatrix
		@test Matrix(b1) == b2
		@test Matrix(b1') == Matrix(b1)'
	end

	true
end

"""
`LinearMapAA(:test)`
self test
"""
function LinearMapAA(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))

	B = 1:6
	L = LinearMap(x -> B*x, y -> B'*y, 6, 1)

	B = reshape(1:6, 6, 1)
	@test Matrix(LinearMapAA(B)) == B

	N = 6; M = N+1
	forw = x -> [cumsum(x); 0] # non-square to stress test
	back = y -> reverse(cumsum(reverse(y[1:N])))

	prop = (name="cumsum", extra=1)
	@test LinearMapAA(forw, (M, N)) isa LinearMapAA
	@test LinearMapAA(forw, (M, N), prop, T=Float64) isa LinearMapAA

	L = LinearMap{Float32}(forw, back, M, N)
	A = LinearMapAA(L, prop)
	Lm = Matrix(L)

#	@test lm_name((Lm, A, I, L, "")) == "MAIL?"

	display(A)

	@test A._lmap == LinearMapAA(L)._lmap
#	@test A == LinearMapAA(forw, back, M, N, prop)
	@test A._prop == LinearMapAA(forw, back, (M, N), prop)._prop
	@test A._lmap == LinearMapAA(forw, back, (M, N), prop)._lmap
	@test A == LinearMapAA(forw, back, (M, N), prop)
	@test A._lmap == LinearMapAA(forw, back, (M, N))._lmap
	@test LinearMapAA(forw, back, (M, N)) isa LinearMapAA
	@test propertynames(A) == (:name, :extra)

	@test issymmetric(A) == false
	@test ishermitian(A) == false
	@test ishermitian(im * A) == false
	@test isposdef(A) == false
	@test issymmetric(A' * A) == true

	@test Matrix(LinearMapAA(L, prop)) == Lm
	@test Matrix(LinearMapAA(L)) == Lm
	@test Matrix(sparse(A)) == Lm

	@test eltype(A) == eltype(L)
	@test Base.eltype(A) == eltype(L) # codecov
	@test ndims(A) == 2
	@test size(A) == size(L)

	B = copy(A)
	@test B == A
	@test !(B === A)

	@test A._prop == prop
	@test A.name == prop.name

	@test_throws String A.bug

	@test Matrix(A)' == Matrix(A')
	@test Matrix(A)' == Matrix(transpose(A))
	@test LinearMapAA_test_getindex(A)
	@test LinearMapAA_test_vmul(A)

	@test LinearMapAA_test_cat(A)
	@test LinearMapAA_test_setindex(A)

	# add / subtract
	@test 2A + 6A isa LinearMapAA
	@test 7A - 2A isa LinearMapAA
	@test Matrix(2A + 6A) == 8 * Matrix(A)
	@test Matrix(7A - 2A) == 5 * Matrix(A)
	@test Matrix(7A - 2*ones(size(A))) == 7 * Matrix(A) - 2*ones(size(A))
	@test Matrix(3*ones(size(A)) - 5A) == 3*ones(size(A)) - 5 * Matrix(A)

	# multiply with identity
	@test Matrix(A * 6I) == 6 * Matrix(A)
	@test Matrix(7I * A) == 7 * Matrix(A)
	@test Matrix((false*I) * A) == zeros(size(A))
	@test Matrix(A * (false*I)) == zeros(size(A))
	@test 1.0I * A === A
	@test A * 1.0I === A
	@test I * A === A
	@test A * I === A

	# multiply
	D = A * A'
	@test D isa LinearMapAA
	@test Matrix(D) == Lm * Lm'
	@test issymmetric(D) == true
	E = A * Lm'
	@test E isa LinearMapAA
	@test Matrix(E) == Lm * Lm'
	F = Lm' * A
	@test F isa LinearMapAA
	@test Matrix(F) == Lm' * Lm
	@test LinearMapAA_test_getindex(F)

	# non-adjoint version
	Af = LinearMapAA(forw, (M, N))
	@test LinearMapAA_test_getindex(Af)
	@test LinearMapAA_test_setindex(Af)

	# kron
#= todo
	M1 = rand(ComplexF64, 3, 3)
	M2 = rand(ComplexF64, 2, 2)
	M3 = kron(M1, M2)
	A1 = LinearMap(M1)
	A2 = LinearMap(M2)
	AK = kron(A1, A2)
	@test AK isa LinearMapAA
	@test Matrix(AK) == M3
	@test Matrix(AK)' == Matrix(AK')
=#

#= todo: cut
	a = ones(3,2)
	b = zeros(3,4)

	A = LinearMapAA(x -> a*x, y -> a'*y, (3, 2))
	B = LinearMapAA(x -> b*x, y -> b'*y, (3, 4))

	Al = [a b I]
	C = [A B I] # now returns a new LinearMap :)
	@test typeof(C) <: LinearMapAA
#	display(C)
#	@show C
	@test Matrix(C) == Matrix(Al)
	@test size(C) == (3,2+4+3)
	@test C * ones(size(C,2)) == 3. * ones(3)
	@test C' * ones(size(C,1)) == [3.0, 3.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0]
	@test Matrix(C)' == Matrix(C')
=#

	true
end

#= mwe for LinearMaps bug
A = LinearMap(ones(3,2))
B = [I A; A I]
Matrix(B) # works fine
Matrix(B') # fails with error DimensionMismatch("mul!")
=#


text:
α = 0.3
δ = 0.1
β = 0.9

kupper = 2
klower = 0.001
n = 10
kgrid = collect(range(klower, stop = kupper, length = n))
v = zeros(length(kgrid))

c = zeros(n, n)
for i in 1:n
    for j in 1:n
        c[i,j] = kgrid[i]^α + (1-δ)*kgrid[i] - kgrid[j]
    end
    # vnew = maximum(log.(c) .+ β*v, dims = 2)
end

c

using Plots
heatmap(kgrid, kgrid, c, dims = 2,
        xlab = "k prime",
        ylab = "k",
        xmirror = true,
        yflip = true)

## Fixing negative consumption problem
vnew = Array{Float64, n}
c = zeros(n, n)
for i in 1:n
    for j in 1:n
        c[i,j] = kgrid[i]^α + (1-δ)*kgrid[i] - kgrid[j]
        if c[i,j] <= 0
            c[i,j] = 10^-10
        end
    end
    vnew = maximum(log.(c) .+ β*v', dims = 2)
end

heatmap(kgrid, kgrid, c, dims = 2,
        xlab = "k prime",
        ylab = "k",
        xmirror = true,
        yflip = true)

scatter(kgrid, vnew)

## Iterate
tolerance = 0.001
imax = 1000
vnew = zeros(length(kgrid))
v = vnew .+ 2*tolerance
global cartesianindex = Array{CartesianIndex{2}, n}
i = 1
while maximum(abs.(v - vnew)) > tolerance && i<=imax
    v = vnew;
    c = zeros(n, n);
    for i in 1:n
        for j in 1:n
            c[i,j] = kgrid[i]^α + (1-δ)*kgrid[i] - kgrid[j];
            if c[i,j] <= 0
                c[i,j] = 10^-10;
            end
        end
        # vnew = maximum(log.(c) .+ β*v', dims = 2);
        (vnew, cartesianindex) = findmax(log.(c) .+ β*v', dims = 2);
    end
    i += 1
end

scatter(kgrid, vnew, title = "v(k)")

kprimeindex = getindex.(cartesianindex, 2)
kprime = kgrid[kprimeindex]
scatter(kgrid, kprime, title = "k'(k)")


text:
#using Javis

function ground(args...)
    background("black") # canvas background
    sethue("white") # pen color
end

function object(p=O, color="black")
    sethue(color)
    circle(p, 25, :fill)
    return p
end

function axes(p1,p2,action,color)
    sethue(color)
    return line(p1,p2,action)
end

function pointer(p)
    sethue("red")
    return circle(p,2,:fill)
end

function connector(con,i)
    map(x->line(x[1],x[2],:stroke),con[1:i])
end

function make_animation()

    connection = []

    frames = 200
    npoints = 30

    points_x = [Point(-200,x) for x in range(-200,stop = 200,length = npoints)][2:end]
    points_y = [Point(x,200) for x in range(-200,stop = 200,length = npoints)][2:end]

    myvideo = Video(500,500)
    Background(1:frames,ground)

    # make the axes
    yaxis = Object(1:frames,(args...)->axes(Point(-200,-200),Point(-200,200),:stroke,"#03eeff"))
    xaxis = Object(1:frames,(args...)->axes(Point(-200,200),Point(200,200),:stroke,"#03eeff"))

    # add some pointers
    Object(1:frames, (args...)->pointer.(points_x))
    Object(1:frames, (args...)->pointer.(points_y))

    # make the connections

    joints = collect(zip(points_x,points_y))
    loop = range(1,stop = frames,length = npoints-1) |> collect .|> x -> round(Int,x)

    for (idx,frame) in enumerate(loop)
        Object(frame:frames,(args...) -> connector(joints,idx))
    end

    render(myvideo; pathname = "parabola_with_lines.gif",framerate = 60)
end

make_animation()

max_stars_repo_name: Ved-Mahajan/graphic-dev

max_issues_repo_name: Ved-Mahajan/JavisIsAwsome

max_forks_repo_name: Ved-Mahajan/JavisIsAwsome


text:
try success(`nvcc --version`)
    cd("../src") do
        run(`make libknet8.so`)
    end
catch
    warn("CUDA not installed, GPU support will not be available.")
end
Base.compilecache("Knet")


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIB = MOI.Bridges
const MOIT = MOI.Test
const MOIU = MOI.Utilities

MOIU.@model(ModelData,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
            (MOI.SecondOrderCone,),
            (MOI.SOS1, MOI.SOS2),
            (),
            (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction,))
const CACHE = MOIU.UniversalFallback(ModelData{Float64}())
const CACHED = MOIU.CachingOptimizer(CACHE, SCIP.Optimizer(display_verblevel=0))
const BRIDGED2 = MOIB.full_bridge_optimizer(CACHED, Float64)

const CONFIG_CACHED = MOIT.Config(
    atol=5e-3, rtol=1e-4,
    exclude=Any[
        MOI.ConstraintDual, MOI.ConstraintName, MOI.DualObjectiveValue, MOI.VariableBasisStatus, MOI.ConstraintBasisStatus,
    ],
)

@testset "MOI Modification tests" begin
    exclude_list = copy(MOI_BASE_EXCLUDED)
    append!(
        exclude_list,
        [
            # Can not delete variable while model contains constraints
            "test_basic_ScalarQuadraticFunction_ZeroOne",
            "test_basic_VectorAffineFunction_NormOneCone",
            "test_basic_VectorAffineFunction_SOS1",
            "test_basic_VectorAffineFunction_SOS2",
            "test_basic_VectorQuadraticFunction_NormOneCone",
            "test_basic_VectorQuadraticFunction_SOS1",
            "test_basic_VectorQuadraticFunction_SOS2",
            "test_linear_integration",
            "test_quadratic_duplicate_terms",
            "test_quadratic_nonhomogeneous",
        ]
    )
    MOIT.runtests(BRIDGED2, CONFIG_CACHED, exclude = exclude_list)
end


text:
#=
This file is auto-generated. Do not edit.
=#

mutable struct TapTransformer <: ACBranch
    name::String
    available::Bool
    activepower_flow::Float64
    reactivepower_flow::Float64
    arc::Arc
    r::Float64  # System per-unit value
    x::Float64  # System per-unit value
    primaryshunt::Float64  # System per-unit value
    tap::Float64
    rate::Union{Nothing, Float64}
    _forecasts::InfrastructureSystems.Forecasts
    internal::InfrastructureSystemsInternal
end

function TapTransformer(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, rate, _forecasts=InfrastructureSystems.Forecasts(), )
    TapTransformer(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, rate, _forecasts, InfrastructureSystemsInternal())
end

function TapTransformer(; name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, rate, _forecasts=InfrastructureSystems.Forecasts(), )
    TapTransformer(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, rate, _forecasts, )
end

# Constructor for demo purposes; non-functional.

function TapTransformer(::Nothing)
    TapTransformer(;
        name="init",
        available=false,
        activepower_flow=0.0,
        reactivepower_flow=0.0,
        arc=Arc(Bus(nothing), Bus(nothing)),
        r=0.0,
        x=0.0,
        primaryshunt=0.0,
        tap=1.0,
        rate=0.0,
        _forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get TapTransformer name."""
get_name(value::TapTransformer) = value.name
"""Get TapTransformer available."""
get_available(value::TapTransformer) = value.available
"""Get TapTransformer activepower_flow."""
get_activepower_flow(value::TapTransformer) = value.activepower_flow
"""Get TapTransformer reactivepower_flow."""
get_reactivepower_flow(value::TapTransformer) = value.reactivepower_flow
"""Get TapTransformer arc."""
get_arc(value::TapTransformer) = value.arc
"""Get TapTransformer r."""
get_r(value::TapTransformer) = value.r
"""Get TapTransformer x."""
get_x(value::TapTransformer) = value.x
"""Get TapTransformer primaryshunt."""
get_primaryshunt(value::TapTransformer) = value.primaryshunt
"""Get TapTransformer tap."""
get_tap(value::TapTransformer) = value.tap
"""Get TapTransformer rate."""
get_rate(value::TapTransformer) = value.rate
"""Get TapTransformer _forecasts."""
get__forecasts(value::TapTransformer) = value._forecasts
"""Get TapTransformer internal."""
get_internal(value::TapTransformer) = value.internal