text:
include("readfile.jl")
include("combat.jl")
include("recursivecombat.jl")

function playscore(player)
    score = 0
    iterator = 1
    for value = Iterators.reverse(player)
       score += value * iterator
       iterator += 1
    end
    return score
end

(player1, player2) = deckplayers("data\\final")

cloneplayer1 = copy(player1)
cloneplayer2 = copy(player2)

combat!(cloneplayer1, cloneplayer2)
result = playscore(cloneplayer1) + playscore(cloneplayer2)
println("Result combat: $result")

winner = recursivecombat!(player1, player2, false)

println("recursive combat winner: $winner")

println("recursive combat player 1: $(playscore(player1))")
println("recursive combat player 2: $(playscore(player2))")


text:
struct Events
    stamps::Vector{Float64} # contains the time stamps
    event::Vector{Symbol}  # name of the function that is being recorded
    startstop::Vector{Symbol} # if the time stamp corresponds to start or to stop
    i::Ref{Int}
end

function Events(n::Int)
    Events(Vector{Float64}(undef, n+1), Vector{Symbol}(undef, n+1), Vector{Symbol}(undef, n+1), Ref{Int}(0))
end

function Base.show(io::IO, calls::Events)
    offset = 0
    if calls.i[] >= length(calls.stamps)
        @warn "The recording buffer was too small, consider increasing it"
    end
    for i in 1:min(calls.i[], length(calls.stamps))
        offset -= calls.startstop[i] == :stop
        foreach(_ -> print(io, " "), 1:max(offset, 0))
        rel_time = calls.stamps[i] - calls.stamps[1]
        println(io, calls.event[i], ": ", rel_time)
        offset += calls.startstop[i] == :start
    end
end

global const to = [Events(1000) for i in 1:1]

function initbuffer!(n)
    while length(to) < Threads.nthreads()
        push!(to, Events(n))
    end
    clear!()
    to
end

"""
    record_start(ev::Symbol)

    record the start of the event, the time stamp is recorded after all counters are
    appropriately increased
"""
record_start(ev::Symbol) = record_start(to[Threads.threadid()], ev)
# function record_start(ev::Symbol)
#     @show (Threads.threadid(), ev)
#     record_start(to[Threads.threadid()], ev)
# end
function record_start(calls, ev::Symbol)
    n = calls.i[] = calls.i[] + 1
    n > length(calls.stamps) && return
    calls.event[n] = ev
    calls.startstop[n] = :start
    calls.stamps[n] = time_ns()
end

"""
    record_end(ev::Symbol)

    record the end of the event, the time stamp is recorded before all counters are
    appropriately increased
"""
record_end(ev::Symbol) = record_end(to[Threads.threadid()], ev::Symbol)
function record_end(calls, ev::Symbol)
    t = time_ns()
    n = calls.i[] = calls.i[] + 1
    n > length(calls.stamps) && return
    calls.event[n] = ev
    calls.startstop[n] = :stop
    calls.stamps[n] = t
end

clear!() = foreach(t -> t.i[] = 0, to)

function Base.resize!(calls::Events, n::Integer)
  resize!(calls.stamps, n)
  resize!(calls.event, n)
  resize!(calls.startstop, n)
end

resizebuffer!(n::Integer) = foreach(t -> resize!(t, n), to)
recorded() = maximum(t.i[] for t in to)
function adjustbuffer!()
    resizebuffer!(2*recorded())
    clear!()
end

ext: jl

max_stars_repo_name: pevnak/LoggingProfiler

max_issues_repo_name: pevnak/LoggingProfiler

max_forks_repo_name: pevnak/LoggingProfiler


text:
#=
Implements the model described in Fogli and Veldkamp, "Germs, Social Networks, and Growth."
@author: Arnav Sood
@date: 2017-09-27

Runs the necessary simulations.
=#

# Import parameters.
include("networksBackend.jl")
include("networksReporting.jl")
using Distributions, networksBackend, LightGraphs, Gadfly, GraphPlot

# Set parameters and create holder objects.
ξ = Bernoulli(0.2)
T = 500 # Paper page 21.
N = 400 # Paper page 22.
F1 = true
F2 = true
λ = 0.008
δ = 80
A₀ = 1 # Paper page 20.
π = 0.12
ϕ = 0.35
θ = Bernoulli(0.14)
nf = [4, 6]
m = [0, 0.2]
Δ = 7
Λ = 0.0008
Π = 0.12

# Bundle the three node parameter distributions into one multivariate node type distribution.
α = NodeDist(θ, nf, m)

# Create model.
m = Model(ξ, T, N, F1, F2, λ, δ, A₀, π, ϕ, α, Δ, Λ, Π)

# Run model 100 times, storing output in memory.
results = Array{ModelState}(100)
@time map!(entry -> runModel(m), results, results) # Can parallelize this eventually.

# Inspect an arbitrary element.
element = results[50]
g = element.g

print(element.nodeList)
adjacency_matrix(g)
gplot(g, layout=circular_layout)

# Run our specialized reporting functions.
rep = report(element)
print(rep)

max_stars_repo_name: ArnavSood/networks-paper

max_issues_repo_name: arnavs/networks-paper

max_forks_repo_name: arnavs/networks-paper


text:
## Matrix factorizations and decompositions

abstract Factorization{T}

macro assertposdef(A, info)
   :(($info)==0 ? $A : throw(PosDefException($info)))
end

macro assertnonsingular(A, info)
   :(($info)==0 ? $A : throw(SingularException($info)))
end

##########################
# Cholesky Factorization #
##########################
immutable Cholesky{T} <: Factorization{T}
    UL::Matrix{T}
    uplo::Char
end
immutable CholeskyPivoted{T} <: Factorization{T}
    UL::Matrix{T}
    uplo::Char
    piv::Vector{BlasInt}
    rank::BlasInt
    tol::Real
    info::BlasInt
end

function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0)
    uplochar = string(uplo)[1]
    if pivot
        A, piv, rank, info = LAPACK.pstrf!(uplochar, A, tol)
        return CholeskyPivoted{T}(A, uplochar, piv, rank, tol, info)
    else
        C, info = LAPACK.potrf!(uplochar, A)
        return @assertposdef Cholesky(C, uplochar) info
    end
end
cholfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0) = cholfact!(copy(A), uplo, pivot=pivot, tol=tol)
cholfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0) = (S = promote_type(typeof(sqrt(one(T))),Float32); S != T ? cholfact!(convert(AbstractMatrix{S},A), uplo, pivot=pivot, tol=tol) : cholfact!(copy(A), uplo, pivot=pivot, tol=tol)) # When julia Cholesky has been implemented, the promotion should be changed.
cholfact(x::Number) = @assertposdef Cholesky(fill(sqrt(x), 1, 1), :U) !(imag(x) == 0 && real(x) > 0)

chol(A::Union(Number, AbstractMatrix), uplo::Symbol) = cholfact(A, uplo)[uplo]
chol(A::Union(Number, AbstractMatrix)) = triu!(cholfact(A, :U).UL)

convert{T}(::Type{Cholesky{T}},C::Cholesky) = Cholesky(convert(AbstractMatrix{T},C.UL),C.uplo)
convert{T}(::Type{Factorization{T}}, C::Cholesky) = convert(Cholesky{T}, C)
convert{T}(::Type{CholeskyPivoted{T}},C::CholeskyPivoted) = CholeskyPivoted(convert(AbstractMatrix{T},C.UL),C.uplo,C.piv,C.rank,C.tol,C.info)
convert{T}(::Type{Factorization{T}}, C::CholeskyPivoted) = convert(CholeskyPivoted{T}, C)

function full{T<:BlasFloat}(C::Cholesky{T})
    if C.uplo == 'U'
        BLAS.trmm!('R', C.uplo, 'N', 'N', one(T), C.UL, tril!(C.UL'))
    else
        BLAS.trmm!('L', C.uplo, 'N', 'N', one(T), C.UL, triu!(C.UL'))
    end
end

size(C::Union(Cholesky, CholeskyPivoted)) = size(C.UL)
size(C::Union(Cholesky, CholeskyPivoted), d::Integer) = size(C.UL,d)

function getindex(C::Cholesky, d::Symbol)
    d == :U && return Triangular(triu!(symbol(C.uplo) == d ? C.UL : C.UL'),:U)
    d == :L && return Triangular(tril!(symbol(C.uplo) == d ? C.UL : C.UL'),:L)
    d == :UL && return Triangular(C.UL, symbol(C.uplo))
    throw(KeyError(d))
end
function getindex{T<:BlasFloat}(C::CholeskyPivoted{T}, d::Symbol)
    d == :U && return triu!(symbol(C.uplo) == d ? C.UL : C.UL')
    d == :L && return tril!(symbol(C.uplo) == d ? C.UL : C.UL')
    d == :p && return C.piv
    if d == :P
        n = size(C, 1)
        P = zeros(T, n, n)
        for i=1:n
            P[C.piv[i],i] = one(T)
        end
        return P
    end
    throw(KeyError(d))
end

show(io::IO, C::Cholesky) = (println("$(typeof(C)) with factor:");show(io,C[symbol(C.uplo)]))

A_ldiv_B!{T<:BlasFloat}(C::Cholesky{T}, B::StridedVecOrMat{T}) = LAPACK.potrs!(C.uplo, C.UL, B)
A_ldiv_B!(C::Cholesky, B::StridedVecOrMat) = C.uplo=='L' ? Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B)) : A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B))

function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedVector{T})
    chkfullrank(C)
    ipermute!(LAPACK.potrs!(C.uplo, C.UL, permute!(B, C.piv)), C.piv)
end
function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedMatrix{T})
    chkfullrank(C)
    n = size(C, 1)
    for i=1:size(B, 2)
        permute!(sub(B, 1:n, i), C.piv)
    end
    LAPACK.potrs!(C.uplo, C.UL, B)
    for i=1:size(B, 2)
        ipermute!(sub(B, 1:n, i), C.piv)
    end
    B
end
A_ldiv_B!(C::CholeskyPivoted, B::StridedVector) = C.uplo=='L' ? Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B[C.piv]))[invperm(C.piv)] : A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B[C.piv]))[invperm(C.piv)]
A_ldiv_B!(C::CholeskyPivoted, B::StridedMatrix) = C.uplo=='L' ? Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B[C.piv,:]))[invperm(C.piv),:] : A_ldiv_B!(Triangular(C.UL,C.uplo,'N'), Ac_ldiv_B!(Triangular(C.UL,C.uplo,'N'), B[C.piv,:]))[invperm(C.piv),:]

function det{T}(C::Cholesky{T})
    dd = one(T)
    for i in 1:size(C.UL,1) dd *= abs2(C.UL[i,i]) end
    dd
end

det{T}(C::CholeskyPivoted{T}) = C.rank<size(C.UL,1) ? real(zero(T)) : prod(abs2(diag(C.UL)))

function logdet{T}(C::Cholesky{T})
    dd = zero(T)
    for i in 1:size(C.UL,1) dd += log(C.UL[i,i]) end
    dd + dd # instead of 2.0dd which can change the type
end

inv(C::Cholesky) = copytri!(LAPACK.potri!(C.uplo, copy(C.UL)), C.uplo, true)

function inv(C::CholeskyPivoted)
    chkfullrank(C)
    ipiv = invperm(C.piv)
    copytri!(LAPACK.potri!(C.uplo, copy(C.UL)), C.uplo, true)[ipiv, ipiv]
end

chkfullrank(C::CholeskyPivoted) = C.rank<size(C.UL, 1) && throw(RankDeficientException(C.info))

rank(C::CholeskyPivoted) = C.rank

####################
# QR Factorization #
####################

immutable QR{T} <: Factorization{T}
    factors::Matrix{T}
    τ::Vector{T}
end
# Note. For QRCompactWY factorization without pivoting, the WY representation based method introduced in LAPACK 3.4
immutable QRCompactWY{S} <: Factorization{S}
    factors::Matrix{S}
    T::Matrix{S}
end

immutable QRPivoted{T} <: Factorization{T}
    factors::Matrix{T}
    τ::Vector{T}
    jpvt::Vector{BlasInt}
end

qrfact!{T<:BlasFloat}(A::StridedMatrix{T}; pivot=false) = pivot ? QRPivoted{T}(LAPACK.geqp3!(A)...) : QRCompactWY(LAPACK.geqrt!(A, min(minimum(size(A)), 36))...)
function qrfact!{T}(A::AbstractMatrix{T}; pivot=false)
    pivot && warn("pivoting only implemented for Float32, Float64, Complex64 and Complex128")
    m, n = size(A)
    τ = zeros(T, min(m,n))
    @inbounds begin
        for k = 1:min(m-1+!(T<:Real),n)
            τk = elementaryLeft!(A, k, k)
            τ[k] = τk
            for j = k+1:n
                vAj = A[k,j]
                for i = k+1:m
                    vAj += conj(A[i,k])*A[i,j]
                end
                vAj = conj(τk)*vAj
                A[k,j] -= vAj
                for i = k+1:m
                    A[i,j] -= A[i,k]*vAj
                end
            end
        end
    end
    QR(A, τ)
end
qrfact{T<:BlasFloat}(A::StridedMatrix{T}; pivot=false) = qrfact!(copy(A),pivot=pivot)
qrfact{T}(A::StridedMatrix{T}; pivot=false) = (S = typeof(one(T)/norm(one(T)));S != T ? qrfact!(convert(AbstractMatrix{S},A), pivot=pivot) : qrfact!(copy(A),pivot=pivot))
qrfact(x::Number) = qrfact(fill(x,1,1))

function qr(A::Union(Number, AbstractMatrix); pivot=false, thin::Bool=true)
    F = qrfact(A, pivot=pivot)
    full(F[:Q], thin=thin), F[:R]
end

convert{T}(::Type{QR{T}},A::QR) = QR(convert(AbstractMatrix{T}, A.factors), convert(Vector{T}, A.τ))
convert{T}(::Type{Factorization{T}}, A::QR) = convert(QR{T}, A)
convert{T}(::Type{QRCompactWY{T}},A::QRCompactWY) = QRCompactWY(convert(AbstractMatrix{T}, A.factors), convert(AbstractMatrix{T}, A.T))
convert{T}(::Type{Factorization{T}}, A::QRCompactWY) = convert(QRCompactWY{T}, A)
convert{T}(::Type{QRPivoted{T}},A::QRPivoted) = QRPivoted(convert(AbstractMatrix{T}, A.factors), convert(Vector{T}, A.τ), A.jpvt)
convert{T}(::Type{Factorization{T}}, A::QRPivoted) = convert(QRPivoted{T}, A)

function getindex(A::QR, d::Symbol)
    d == :R && return triu(A.factors[1:minimum(size(A)),:])
    d == :Q && return QRPackedQ(A.factors,A.τ)
    throw(KeyError(d))
end
function getindex(A::QRCompactWY, d::Symbol)
    d == :R && return triu(A.factors[1:minimum(size(A)),:])
    d == :Q && return QRCompactWYQ(A.factors,A.T)
    throw(KeyError(d))
end
function getindex{T}(A::QRPivoted{T}, d::Symbol)
    d == :R && return triu(A.factors[1:minimum(size(A)),:])
    d == :Q && return QRPackedQ(A.factors,A.τ)
    d == :p && return A.jpvt
    if d == :P
        p = A[:p]
        n = length(p)
        P = zeros(T, n, n)
        for i in 1:n
            P[p[i],i] = one(T)
        end
        return P
    end
    throw(KeyError(d))
end

immutable QRPackedQ{T} <: AbstractMatrix{T}
    factors::Matrix{T}
    τ::Vector{T}
end
immutable QRCompactWYQ{S} <: AbstractMatrix{S}
    factors::Matrix{S}
    T::Matrix{S}
end

convert{T}(::Type{QRPackedQ{T}}, Q::QRPackedQ) = QRPackedQ(convert(AbstractMatrix{T}, Q.factors), convert(Vector{T}, Q.τ))
convert{T}(::Type{AbstractMatrix{T}}, Q::QRPackedQ) = convert(QRPackedQ{T}, Q)
convert{S}(::Type{QRCompactWYQ{S}}, Q::QRCompactWYQ) = QRCompactWYQ(convert(AbstractMatrix{S}, Q.factors), convert(AbstractMatrix{S}, Q.T))
convert{S}(::Type{AbstractMatrix{S}}, Q::QRCompactWYQ) = convert(QRCompactWYQ{S}, Q)

size(A::Union(QR,QRCompactWY,QRPivoted), dim::Integer) = size(A.factors, dim)
size(A::Union(QR,QRCompactWY,QRPivoted)) = size(A.factors)
size(A::Union(QRPackedQ,QRCompactWYQ), dim::Integer) = 0 < dim ? (dim <= 2 ? size(A.factors, 1) : 1) : throw(BoundsError())
size(A::Union(QRPackedQ,QRCompactWYQ)) = size(A, 1), size(A, 2)

full{T}(A::Union(QRPackedQ{T},QRCompactWYQ{T}); thin::Bool=true) = A_mul_B!(A, thin ? eye(T, size(A.factors,1), minimum(size(A.factors))) : eye(T, size(A.factors,1)))

## Multiplication by Q
### QB
A_mul_B!{T<:BlasFloat}(A::QRCompactWYQ{T}, B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','N',A.factors,A.T,B)
A_mul_B!{T<:BlasFloat}(A::QRPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr!('L','N',A.factors,A.τ,B)
function A_mul_B!{T}(A::QRPackedQ{T}, B::AbstractVecOrMat{T})
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    mA == mB || throw(DimensionMismatch(""))
    Afactors = A.factors
    @inbounds begin
        for k = min(mA,nA):-1:1
            for j = 1:nB
                vBj = B[k,j]
                for i = k+1:mB
                    vBj += conj(Afactors[i,k])*B[i,j]
                end
                vBj = A.τ[k]*vBj
                B[k,j] -= vBj
                for i = k+1:mB
                    B[i,j] -= Afactors[i,k]*vBj
                end
            end
        end
    end
    B
end
function *{TA,Tb}(A::Union(QRPackedQ{TA},QRCompactWYQ{TA}), b::StridedVector{Tb})
    TAb = promote_type(TA,Tb)
    Anew = convert(AbstractMatrix{TAb},A)
    bnew = size(A.factors,1) == length(b) ? (Tb == TAb ? copy(b) : convert(Vector{TAb}, b)) : (size(A.factors,2) == length(b) ? [b,zeros(TAb, size(A.factors,1)-length(b))] : throw(DimensionMismatch("")))
    A_mul_B!(Anew,bnew)
end
function *{TA,TB}(A::Union(QRPackedQ{TA},QRCompactWYQ{TA}), B::StridedMatrix{TB})
    TAB = promote_type(TA,TB)
    Anew = convert(AbstractMatrix{TAB},A)
    Bnew = size(A.factors,1) == size(B,1) ? (TB == TAB ? copy(B) : convert(AbstractMatrix{TAB}, B)) : (size(A.factors,2) == size(B,1) ? [B;zeros(TAB, size(A.factors,1)-size(B,1),size(B,2))] : throw(DimensionMismatch("")))
    A_mul_B!(Anew,Bnew)
end
### QcB
Ac_mul_B!{T<:BlasReal}(A::QRCompactWYQ{T}, B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','T',A.factors,A.T,B)
Ac_mul_B!{T<:BlasComplex}(A::QRCompactWYQ{T}, B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','C',A.factors,A.T,B)
Ac_mul_B!{T<:BlasReal}(A::QRPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr!('L','T',A.factors,A.τ,B)
Ac_mul_B!{T<:BlasComplex}(A::QRPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr!('L','C',A.factors,A.τ,B)
function Ac_mul_B!{T}(A::QRPackedQ{T}, B::AbstractVecOrMat{T})
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    mA == mB || throw(DimensionMismatch(""))
    Afactors = A.factors
    @inbounds begin
        for k = 1:min(mA,nA)
            for j = 1:nB
                vBj = B[k,j]
                for i = k+1:mB
                    vBj += conj(Afactors[i,k])*B[i,j]
                end
                vBj = conj(A.τ[k])*vBj
                B[k,j] -= vBj
                for i = k+1:mB
                    B[i,j] -= Afactors[i,k]*vBj
                end
            end
        end
    end
    B
end
function Ac_mul_B{TQ<:Number,TB<:Number,N}(Q::Union(QRPackedQ{TQ},QRCompactWYQ{TQ}), B::StridedArray{TB,N})
    TQB = promote_type(TQ,TB)
    Ac_mul_B!(convert(AbstractMatrix{TQB}, Q), TB == TQB ? copy(B) : convert(AbstractArray{TQB,N}, B))
end
### AQ
A_mul_B!{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) = LAPACK.gemqrt!('R','N', B.factors, B.T, A)
A_mul_B!(A::StridedVecOrMat, B::QRPackedQ) = LAPACK.ormqr!('R', 'N', B.factors, B.τ, A)
function A_mul_B!{T}(A::StridedMatrix{T},Q::QRPackedQ{T})
    mQ, nQ = size(Q.factors)
    mA, nA = size(A,1), size(A,2)
    nA == mQ || throw(DimensionMismatch(""))
    Qfactors = Q.factors
    @inbounds begin
        for k = 1:min(mQ,nQ)
            for i = 1:mA
                vAi = A[i,k]
                for j = k+1:mQ
                    vAi += A[i,j]*Qfactors[j,k]
                end
                vAi = vAi*Q.τ[k]
                A[i,k] -= vAi
                for j = k+1:nA
                    A[i,j] -= vAi*conj(Qfactors[j,k])
                end
            end
        end
    end
    A
end
function *{TA,TQ,N}(A::StridedArray{TA,N}, Q::Union(QRPackedQ{TQ},QRCompactWYQ{TQ}))
    TAQ = promote_type(TA, TQ)
    A_mul_B!(TA==TAQ ? copy(A) : convert(AbstractArray{TAQ,N}, A), convert(AbstractMatrix{TAQ}, Q))
end
### AQc
A_mul_Bc!{T<:BlasReal}(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) = LAPACK.gemqrt!('R','T',B.factors,B.T,A)
A_mul_Bc!{T<:BlasComplex}(A::StridedVecOrMat{T}, B::QRCompactWYQ{T}) = LAPACK.gemqrt!('R','C',B.factors,B.T,A)
A_mul_Bc!{T<:BlasReal}(A::StridedVecOrMat{T}, B::QRPackedQ{T}) = LAPACK.ormqr!('R','T',B.factors,B.τ,A)
A_mul_Bc!{T<:BlasComplex}(A::StridedVecOrMat{T}, B::QRPackedQ{T}) = LAPACK.ormqr!('R','C',B.factors,B.τ,A)
function A_mul_Bc!{T}(A::AbstractMatrix{T},Q::QRPackedQ{T})
    mQ, nQ = size(Q.factors)
    mA, nA = size(A,1), size(A,2)
    nA == mQ || throw(DimensionMismatch(""))
    Qfactors = Q.factors
    @inbounds begin
        for k = min(mQ,nQ):-1:1
            for i = 1:mA
                vAi = A[i,k]
                for j = k+1:mQ
                    vAi += A[i,j]*Qfactors[j,k]
                end
                vAi = vAi*conj(Q.τ[k])
                A[i,k] -= vAi
                for j = k+1:nA
                    A[i,j] -= vAi*conj(Qfactors[j,k])
                end
            end
        end
    end
    A
end
function A_mul_Bc{TA,TB,N}(A::AbstractArray{TA,N}, B::Union(QRCompactWYQ{TB},QRPackedQ{TB}))
    TAB = promote_type(TA,TB)
    A_mul_Bc!(size(A,2)==size(B.factors,1) ? (TA == TAB ? copy(A) : convert(AbstractMatrix{TAB,N}, A)) : (size(A,2)==size(B.factors,2) ? [A zeros(TAB, size(A, 1), size(B.factors, 1) - size(B.factors, 2))] : throw(DimensionMismatch(""))),convert(AbstractMatrix{TAB}, B))
end

# Julia implementation similarly to xgelsy
function A_ldiv_B!{T<:BlasFloat}(A::Union(QRCompactWY{T},QRPivoted{T}), B::StridedMatrix{T}, rcond::Real)
    mA, nA = size(A.factors)
    nr = min(mA,nA)
    nrhs = size(B, 2)
    if nr == 0 return zeros(0, nrhs), 0 end
    ar = abs(A.factors[1])
    if ar == 0 return zeros(nr, nrhs), 0 end
    rnk = 1
    xmin = ones(T, nr)
    xmax = ones(T, nr)
    tmin = tmax = ar
    while rnk < nr
        tmin, smin, cmin = LAPACK.laic1!(2, sub(xmin, 1:rnk), tmin, sub(A.factors, 1:rnk, rnk + 1), A.factors[rnk + 1, rnk + 1])
        tmax, smax, cmax = LAPACK.laic1!(1, sub(xmax, 1:rnk), tmax, sub(A.factors, 1:rnk, rnk + 1), A.factors[rnk + 1, rnk + 1])
        tmax*rcond > tmin && break
        xmin[1:rnk + 1] = [smin*sub(xmin, 1:rnk), cmin]
        xmax[1:rnk + 1] = [smax*sub(xmin, 1:rnk), cmax]
        rnk += 1
        # if cond(r[1:rnk, 1:rnk])*rcond < 1 break end
    end
    C, τ = LAPACK.tzrzf!(A.factors[1:rnk,:])
    A_ldiv_B!(Triangular(C[1:rnk,1:rnk],:U),sub(Ac_mul_B!(A[:Q],sub(B, 1:mA, 1:nrhs)),1:rnk,1:nrhs))
    B[rnk+1:end,:] = zero(T)
    LAPACK.ormrz!('L', iseltype(B, Complex) ? 'C' : 'T', C, τ, sub(B,1:nA,1:nrhs))
    return isa(A,QRPivoted) ? B[invperm(A[:p]),:] : B[1:nA,:], rnk
end
A_ldiv_B!{T<:BlasFloat}(A::Union(QRCompactWY{T},QRPivoted{T}), B::StridedVector{T}) = A_ldiv_B!(A,reshape(B,length(B),1))[:]
A_ldiv_B!{T<:BlasFloat}(A::Union(QRCompactWY{T},QRPivoted{T}), B::StridedVecOrMat{T}) = A_ldiv_B!(A, B, sqrt(eps(real(float(one(eltype(B)))))))[1]
function A_ldiv_B!{T}(A::QR{T},B::StridedMatrix{T})
    m, n = size(A)
    minmn = min(m,n)
    mB, nB = size(B)
    Ac_mul_B!(A[:Q],sub(B,1:m,1:nB)) # Reconsider when arrayviews are merged.
    R = A[:R]
    @inbounds begin
        if n > m # minimum norm solution
            τ = zeros(T,m)
            for k = m:-1:1 # Trapezoid to triangular by elementary operation
                τ[k] = elementaryRightTrapezoid!(R,k)
                for i = 1:k-1
                    vRi = R[i,k]
                    for j = m+1:n
                        vRi += R[i,j]*R[k,j]
                    end
                    vRi *= τ[k]
                    R[i,k] -= vRi
                    for j = m+1:n
                        R[i,j] -= vRi*R[k,j]
                    end
                end
            end
        end
        for k = 1:nB # solve triangular system. When array views are implemented, consider exporting    to function.
            for i = minmn:-1:1
                for j = i+1:minmn
                    B[i,k] -= R[i,j]*B[j,k]
                end
                B[i,k] /= R[i,i]
            end
        end
        if n > m # Apply elemenary transformation to solution
            B[m+1:mB,1:nB] = zero(T)
            for j = 1:nB
                for k = 1:m
                    vBj = B[k,j]
                    for i = m+1:n
                        vBj += B[i,j]*conj(R[k,i])
                    end
                    vBj *= τ[k]
                    B[k,j] -= vBj
                    for i = m+1:n
                        B[i,j] -= R[k,i]*vBj
                    end
                end
            end
        end
    end
    return B[1:n,:]
end
A_ldiv_B!(A::QR, B::StridedVector) = A_ldiv_B!(A, reshape(B, length(B), 1))[:]
A_ldiv_B!(A::QRPivoted, B::StridedVector) = A_ldiv_B!(QR(A.factors,A.τ),B)[invperm(A.jpvt)]
A_ldiv_B!(A::QRPivoted, B::StridedMatrix) = A_ldiv_B!(QR(A.factors,A.τ),B)[invperm(A.jpvt),:]
function \{TA,Tb}(A::Union(QR{TA},QRCompactWY{TA},QRPivoted{TA}),b::StridedVector{Tb})
    S = promote_type(TA,Tb)
    m,n = size(A)
    m == length(b) || throw(DimensionMismatch("left hand side has $(m) rows, but right hand side has length $(length(b))"))
    n > m ? A_ldiv_B!(convert(Factorization{S},A),[b,zeros(S,n-m)]) : A_ldiv_B!(convert(Factorization{S},A), S == Tb ? copy(b) : convert(AbstractVector{S}, b))
end
function \{TA,TB}(A::Union(QR{TA},QRCompactWY{TA},QRPivoted{TA}),B::StridedMatrix{TB})
    S = promote_type(TA,TB)
    m,n = size(A)
    m == size(B,1) || throw(DimensionMismatch("left hand side has $(m) rows, but right hand side has $(size(B,1)) rows"))
    n > m ? A_ldiv_B!(convert(Factorization{S},A),[B;zeros(S,n-m,size(B,2))]) : A_ldiv_B!(convert(Factorization{S},A), S == TB ? copy(B) : convert(AbstractMatrix{S}, B))
end

##TODO:  Add methods for rank(A::QRP{T}) and adjust the (\) method accordingly
##       Add rcond methods for Cholesky, LU, QR and QRP types
## Lower priority: Add LQ, QL and RQ factorizations

# FIXME! Should add balancing option through xgebal
immutable Hessenberg{T} <: Factorization{T}
    factors::Matrix{T}
    τ::Vector{T}
end
Hessenberg(A::StridedMatrix) = Hessenberg(LAPACK.gehrd!(A)...)

hessfact!{T<:BlasFloat}(A::StridedMatrix{T}) = Hessenberg(A)
hessfact{T<:BlasFloat}(A::StridedMatrix{T}) = hessfact!(copy(A))
hessfact{T}(A::StridedMatrix{T}) = (S = promote_type(Float32,typeof(one(T)/norm(one(T)))); S != T ? hessfact!(convert(AbstractMatrix{S},A)) : hessfact!(copy(A)))

immutable HessenbergQ{T} <: AbstractMatrix{T}
    factors::Matrix{T}
    τ::Vector{T}
end
HessenbergQ(A::Hessenberg) = HessenbergQ(A.factors, A.τ)
size(A::HessenbergQ, args...) = size(A.factors, args...)

function getindex(A::Hessenberg, d::Symbol)
    d == :Q && return HessenbergQ(A)
    d == :H && return triu(A.factors, -1)
    throw(KeyError(d))
end

full(A::HessenbergQ) = LAPACK.orghr!(1, size(A.factors, 1), copy(A.factors), A.τ)

# Also printing of QRQs
print_matrix(io::IO, A::Union(QRPackedQ,QRCompactWYQ,HessenbergQ), rows::Integer, cols::Integer, punct...) = print_matrix(io, full(A), rows, cols, punct...)

#######################
# Eigendecompositions #
#######################

# Eigenvalues
immutable Eigen{T,V} <: Factorization{T}
    values::Vector{V}
    vectors::Matrix{T}
end

# Generalized eigenvalue problem.
immutable GeneralizedEigen{T,V} <: Factorization{T}
    values::Vector{V}
    vectors::Matrix{T}
end

function getindex(A::Union(Eigen,GeneralizedEigen), d::Symbol)
    d == :values && return A.values
    d == :vectors && return A.vectors
    throw(KeyError(d))
end

isposdef(A::Union(Eigen,GeneralizedEigen)) = all(A.values .> 0)

function eigfact!{T<:BlasReal}(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true)
    n = size(A, 2)
    n==0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
    issym(A) && return eigfact!(Symmetric(A))
    A, WR, WI, VL, VR, _ = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)
    all(WI .== 0.) && return Eigen(WR, VR)
    evec = zeros(Complex{T}, n, n)
    j = 1
    while j <= n
        if WI[j] == 0.0
            evec[:,j] = VR[:,j]
        else
            evec[:,j]   = VR[:,j] + im*VR[:,j+1]
            evec[:,j+1] = VR[:,j] - im*VR[:,j+1]
            j += 1
        end
        j += 1
    end
    return Eigen(complex(WR, WI), evec)
end

function eigfact!{T<:BlasComplex}(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true)
    n = size(A, 2)
    n == 0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
    ishermitian(A) && return eigfact!(Hermitian(A))
    return Eigen(LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)[[2,4]]...)
end
eigfact{T}(A::AbstractMatrix{T}, args...; kwargs...) = (S = promote_type(Float32,typeof(one(T)/norm(one(T)))); S != T ? eigfact!(convert(AbstractMatrix{S}, A), args...; kwargs...) : eigfact!(copy(A), args...; kwargs...))
eigfact(x::Number) = Eigen([x], fill(one(x), 1, 1))

# function eig(A::Union(Number, AbstractMatrix); permute::Bool=true, scale::Bool=true)
#     F = eigfact(A, permute=permute, scale=scale)
#     F[:values], F[:vectors]
# end
function eig(A::Union(Number, AbstractMatrix), args...; kwargs...)
    F = eigfact(A, args..., kwargs...)
    F[:values], F[:vectors]
end
#Calculates eigenvectors
eigvecs(A::Union(Number, AbstractMatrix), args...; kwargs...) = eigfact(A, args...; kwargs...)[:vectors]

function eigvals!{T<:BlasReal}(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true)
    issym(A) && return eigvals!(Symmetric(A))
    _, valsre, valsim, _ = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'N', 'N', A)
    return all(valsim .== 0) ? valsre : complex(valsre, valsim)
end
function eigvals!{T<:BlasComplex}(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true)
    ishermitian(A) && return eigvals(Hermitian(A))
    return LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'N', 'N', A)[2]
end
eigvals{T}(A::AbstractMatrix{T}, args...; kwargs...) = (S = promote_type(Float32,typeof(one(T)/norm(one(T)))); S != T ? eigvals!(convert(AbstractMatrix{S}, A), args...; kwargs...) : eigvals!(copy(A), args...; kwargs...))
eigvals{T<:Number}(x::T; kwargs...) = (val = convert(promote_type(Float32,typeof(one(T)/norm(one(T)))),x); imag(val) == 0 ? [real(val)] : [val])

#Computes maximum and minimum eigenvalue
function eigmax(A::Union(Number, AbstractMatrix); kwargs...)
    v = eigvals(A; kwargs...)
    iseltype(v,Complex) ? error("DomainError: complex eigenvalues cannot be ordered") : maximum(v)
end
function eigmin(A::Union(Number, AbstractMatrix); kwargs...)
    v = eigvals(A; kwargs...)
    iseltype(v,Complex) ? error("DomainError: complex eigenvalues cannot be ordered") : minimum(v)
end

inv(A::Eigen) = A.vectors/Diagonal(A.values)*A.vectors'
det(A::Eigen) = prod(A.values)

# Generalized eigenproblem
function eigfact!{T<:BlasReal}(A::StridedMatrix{T}, B::StridedMatrix{T})
    issym(A) && isposdef(B) && return eigfact!(Symmetric(A), Symmetric(B))
    n = size(A, 1)
    alphar, alphai, beta, _, vr = LAPACK.ggev!('N', 'V', A, B)
    all(alphai .== 0) && return GeneralizedEigen(alphar ./ beta, vr)

    vecs = zeros(Complex{T}, n, n)
    j = 1
    while j <= n
        if alphai[j] == 0.0
            vecs[:,j] = vr[:,j]
        else
            vecs[:,j  ] = vr[:,j] + im*vr[:,j+1]
            vecs[:,j+1] = vr[:,j] - im*vr[:,j+1]
            j += 1
        end
        j += 1
    end
    return GeneralizedEigen(complex(alphar, alphai)./beta, vecs)
end

function eigfact!{T<:BlasComplex}(A::StridedMatrix{T}, B::StridedMatrix{T})
    ishermitian(A) && isposdef(B) && return eigfact!(Hermitian(A), Hermitian(B))
    alpha, beta, _, vr = LAPACK.ggev!('N', 'V', A, B)
    return GeneralizedEigen(alpha./beta, vr)
end
eigfact{TA,TB}(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}) = (S = promote_type(Float32,typeof(one(TA)/norm(one(TA))),TB); eigfact!(S != TA ? convert(AbstractMatrix{S},A) : copy(A), S != TB ? convert(AbstractMatrix{S},B) : copy(B)))

function eigvals!{T<:BlasReal}(A::StridedMatrix{T}, B::StridedMatrix{T})
    issym(A) && isposdef(B) && return eigvals!(Symmetric(A), Symmetric(B))
    alphar, alphai, beta, vl, vr = LAPACK.ggev!('N', 'N', A, B)
    (all(alphai .== 0) ? alphar : complex(alphar, alphai))./beta
end
function eigvals!{T<:BlasComplex}(A::StridedMatrix{T}, B::StridedMatrix{T})
    ishermitian(A) && isposdef(B) && return eigvals!(Hermitian(A), Hermitian(B))
    alpha, beta, vl, vr = LAPACK.ggev!('N', 'N', A, B)
    alpha./beta
end
eigvals{TA,TB}(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}) = (S = promote_type(Float32,typeof(one(TA)/norm(one(TA))),TB); eigvals!(S != TA ? convert(AbstractMatrix{S},A) : copy(A), S != TB ? convert(AbstractMatrix{S},B) : copy(B)))

# SVD
immutable SVD{T<:BlasFloat,Tr} <: Factorization{T}
    U::Matrix{T}
    S::Vector{Tr}
    Vt::Matrix{T}
end
function svdfact!{T<:BlasFloat}(A::StridedMatrix{T}; thin::Bool=true)
    m,n = size(A)
    if m == 0 || n == 0
        u,s,vt = (eye(T, m, thin ? n : m), real(zeros(T,0)), eye(T,n,n))
    else
        u,s,vt = LAPACK.gesdd!(thin ? 'S' : 'A', A)
    end
    SVD(u,s,vt)
end
svdfact{T<:BlasFloat}(A::StridedMatrix{T};thin=true) = svdfact!(copy(A),thin=thin)
svdfact{T}(A::StridedVecOrMat{T};thin=true) = (S = promote_type(Float32,typeof(one(T)/norm(one(T)))); S != T ? svdfact!(convert(AbstractMatrix{S},A),thin=thin) : svdfact!(copy(A),thin=thin))
svdfact(x::Number; thin::Bool=true) = SVD(x == 0 ? fill(one(x), 1, 1) : fill(x/abs(x), 1, 1), [abs(x)], fill(one(x), 1, 1))
svdfact(x::Integer; thin::Bool=true) = svdfact(float(x), thin=thin)

function svd(A::Union(Number, AbstractArray); thin::Bool=true)
    F = svdfact(A, thin=thin)
    F.U, F.S, F.Vt'
end

function getindex(F::SVD, d::Symbol)
    d == :U && return F.U
    d == :S && return F.S
    d == :Vt && return F.Vt
    d == :V && return F.Vt'
    throw(KeyError(d))
end

svdvals!{T<:BlasFloat}(A::StridedMatrix{T}) = any([size(A)...].==0) ? zeros(T, 0) : LAPACK.gesdd!('N', A)[2]
svdvals{T<:BlasFloat}(A::StridedMatrix{T}) = svdvals!(copy(A))
svdvals{T}(A::StridedMatrix{T}) = (S = promote_type(Float32,typeof(one(T)/norm(one(T)))); S != T ? svdvals!(convert(AbstractMatrix{S}, A)) : svdvals!(copy(A)))
svdvals(x::Number) = [abs(x)]

# SVD least squares
function \{T<:BlasFloat}(A::SVD{T}, B::StridedVecOrMat{T})
    n = length(A.S)
    Sinv = zeros(T, n)
    Sinv[A.S .> sqrt(eps())] = 1.0 ./ A.S
    scale(A.Vt', Sinv) * A.U[:,1:n]'B
end

# Generalized svd
immutable GeneralizedSVD{T} <: Factorization{T}
    U::Matrix{T}
    V::Matrix{T}
    Q::Matrix{T}
    a::Vector
    b::Vector
    k::Int
    l::Int
    R::Matrix{T}
end

function svdfact!{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T})
    U, V, Q, a, b, k, l, R = LAPACK.ggsvd!('U', 'V', 'Q', A, B)
    GeneralizedSVD(U, V, Q, a, b, int(k), int(l), R)
end
svdfact{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}) = svdfact!(copy(A),copy(B))
svdfact{TA,TB}(A::StridedMatrix{TA}, B::StridedMatrix{TB}) = (S = promote_type(Float32,typeof(one(TA)/norm(one(TA))),TB); svdfact!(S != TA ? convert(AbstractMatrix{S},A) : copy(A), S != TB ? convert(AbstractMatrix{S},B) : copy(B)))

function svd(A::AbstractMatrix, B::AbstractMatrix)
    F = svdfact(A, B)
    F[:U], F[:V], F[:Q], F[:D1], F[:D2], F[:R0]
end

function getindex{T}(obj::GeneralizedSVD{T}, d::Symbol)
    d == :U && return obj.U
    d == :V && return obj.V
    d == :Q && return obj.Q
    (d == :alpha || d == :a) && return obj.a
    (d == :beta || d == :b) && return obj.b
    (d == :vals || d == :S) && return obj.a[1:obj.k + obj.l] ./ obj.b[1:obj.k + obj.l]
    if d == :D1
        m = size(obj.U, 1)
        if m - obj.k - obj.l >= 0
            return [eye(T, obj.k) zeros(T, obj.k, obj.l); zeros(T, obj.l, obj.k) diagm(obj.a[obj.k + 1:obj.k + obj.l]); zeros(T, m - obj.k - obj.l, obj.k + obj.l)]
        else
            return [eye(T, m, obj.k) [zeros(T, obj.k, m - obj.k); diagm(obj.a[obj.k + 1:m])] zeros(T, m, obj.k + obj.l - m)]
        end
    end
    if d == :D2
        m = size(obj.U, 1)
        p = size(obj.V, 1)
        if m - obj.k - obj.l >= 0
            return [zeros(T, obj.l, obj.k) diagm(obj.b[obj.k + 1:obj.k + obj.l]); zeros(T, p - obj.l, obj.k + obj.l)]
        else
            return [zeros(T, p, obj.k) [diagm(obj.b[obj.k + 1:m]); zeros(T, obj.k + p - m, m - obj.k)] [zeros(T, m - obj.k, obj.k + obj.l - m); eye(T, obj.k + p - m, obj.k + obj.l - m)]]
        end
    end
    d == :R && return obj.R
    if d == :R0
        n = size(obj.Q, 1)
        return [zeros(T, obj.k + obj.l, n - obj.k - obj.l) obj.R]
    end
    throw(KeyError(d))
end

function svdvals!{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T})
    _, _, _, a, b, k, l, _ = LAPACK.ggsvd!('N', 'N', 'N', A, B)
    a[1:k + l] ./ b[1:k + l]
end
svdvals{T<:BlasFloat}(A::StridedMatrix{T},B::StridedMatrix{T}) = svdvals!(copy(A),copy(B))
svdvals{TA,TB}(A::StridedMatrix{TA}, B::StridedMatrix{TB}) = (S = promote_type(Float32,typeof(one(T)/norm(one(TA))),TB); svdvals!(S != TA ? convert(AbstractMatrix{S}, A) : copy(A), S != TB ? convert(AbstractMatrix{S}, B) : copy(B)))

immutable Schur{Ty<:BlasFloat} <: Factorization{Ty}
    T::Matrix{Ty}
    Z::Matrix{Ty}
    values::Vector
end

schurfact!{T<:BlasFloat}(A::StridedMatrix{T}) = Schur(LinAlg.LAPACK.gees!('V', A)...)
schurfact{T<:BlasFloat}(A::StridedMatrix{T}) = schurfact!(copy(A))
schurfact{T}(A::StridedMatrix{T}) = (S = promote_type(Float32,typeof(one(T)/norm(one(T)))); S != T ? schurfact!(convert(AbstractMatrix{S},A)) : schurfact!(copy(A)))

function getindex(F::Schur, d::Symbol)
    (d == :T || d == :Schur) && return F.T
    (d == :Z || d == :vectors) && return F.Z
    d == :values && return F.values
    throw(KeyError(d))
end

function schur(A::AbstractMatrix)
    SchurF = schurfact(A)
    SchurF[:T], SchurF[:Z], SchurF[:values]
end

immutable GeneralizedSchur{Ty<:BlasFloat} <: Factorization{Ty}
    S::Matrix{Ty}
    T::Matrix{Ty}
    alpha::Vector
    beta::Vector{Ty}
    Q::Matrix{Ty}
    Z::Matrix{Ty}
end

schurfact!{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}) = GeneralizedSchur(LinAlg.LAPACK.gges!('V', 'V', A, B)...)
schurfact{T<:BlasFloat}(A::StridedMatrix{T},B::StridedMatrix{T}) = schurfact!(copy(A),copy(B))
schurfact{TA,TB}(A::StridedMatrix{TA}, B::StridedMatrix{TB}) = (S = promote_type(Float32,typeof(one(TA)/norm(one(TA))),TB); schurfact!(S != TA ? convert(AbstractMatrix{S},A) : copy(A), S != TB ? convert(AbstractMatrix{S},B) : copy(B)))

function getindex(F::GeneralizedSchur, d::Symbol)
    d == :S && return F.S
    d == :T && return F.T
    d == :alpha && return F.alpha
    d == :beta && return F.beta
    d == :values && return F.alpha./F.beta
    (d == :Q || d == :left) && return F.Q
    (d == :Z || d == :right) && return F.Z
    throw(KeyError(d))
end

function schur(A::AbstractMatrix, B::AbstractMatrix)
    SchurF = schurfact(A, B)
    SchurF[:S], SchurF[:T], SchurF[:Q], SchurF[:Z]
end

### General promotion rules
convert{T}(::Type{Factorization{T}}, F::Factorization{T}) = F
inv{T}(F::Factorization{T}) = A_ldiv_B!(F, eye(T, size(F,1)))
function \{TF<:Number,TB<:Number,N}(F::Factorization{TF}, B::AbstractArray{TB,N})
    TFB = typeof(one(TF)/one(TB))
    A_ldiv_B!(convert(Factorization{TFB}, F), TB == TFB ? copy(B) : convert(AbstractArray{TFB,N}, B))
end

function Ac_ldiv_B{TF<:Number,TB<:Number,N}(F::Factorization{TF}, B::AbstractArray{TB,N})
    TFB = typeof(one(TF)/one(TB))
    Ac_ldiv_B!(convert(Factorization{TFB}, F), TB == TFB ? copy(B) : convert(AbstractArray{TFB,N}, B))
end

function At_ldiv_B{TF<:Number,TB<:Number,N}(F::Factorization{TF}, B::AbstractArray{TB,N})
    TFB = typeof(one(TF)/one(TB))
    At_ldiv_B!(convert(Factorization{TFB}, F), TB == TFB ? copy(B) : convert(AbstractArray{TFB,N}, B))
end

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


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

using Markdown
using InteractiveUtils

# ╔═╡ c7faa3e6-72fd-11eb-1b1c-cff28ec840c7
#using Pkg; Pkg.activate(@__DIR__())

# ╔═╡ 2365b56c-72fd-11eb-2fb3-fb247a857b01
using CairoMakie

# ╔═╡ 01044994-72fe-11eb-1e12-c1cf87f4216c
using DataFrames

# ╔═╡ 5cfe5194-72ff-11eb-329a-9934d15a20a9
md"""
# A Dummy Notebook for PackageCompiler.jl
"""

# ╔═╡ 415a2582-72fd-11eb-207e-a7459239116c
CairoMakie.activate!(type = "png")

# ╔═╡ 47bbca66-72fd-11eb-2485-ffe43809bbc3
begin
	fig = Figure()
	ax = Axis(fig[1,1])

	scatter!(ax, rand(100), rand(100), label = "bla")

	Legend(fig[1,2], ax)

	fig |> display
end

# ╔═╡ 033b6f1c-72fe-11eb-1cba-7da643d154c9
DataFrame(x = rand(100))

# ╔═╡ Cell order:
# ╟─5cfe5194-72ff-11eb-329a-9934d15a20a9
# ╠═c7faa3e6-72fd-11eb-1b1c-cff28ec840c7
# ╠═2365b56c-72fd-11eb-2fb3-fb247a857b01
# ╠═415a2582-72fd-11eb-207e-a7459239116c
# ╠═47bbca66-72fd-11eb-2485-ffe43809bbc3
# ╠═01044994-72fe-11eb-1e12-c1cf87f4216c
# ╠═033b6f1c-72fe-11eb-1cba-7da643d154c9

max_stars_repo_name: greimel/pluto-binder

max_issues_repo_name: greimel/pluto-binder

max_forks_repo_name: greimel/pluto-binder


text:
@testset verbose=true "datapoint.jl " begin
    ## Check _to_tuple functions
    @testset "_to_tuple    " begin
        # _to_tuple should distribute real numbers and promote mixed tuples
        @test LatSpec._to_tuple(1)     == (1, 1)        # int   -> tuple{int}
        @test LatSpec._to_tuple(1.5)   == (1.5, 1.5)    # float -> tuple{float}
        @test LatSpec._to_tuple((1,1)) == (1, 1)        # tuple{int} -> tuple{int}
        @test typeof(LatSpec._to_tuple((1.5,1))) != typeof((1.5, 1))    # tuple{float, int} -> tuple{float, float}
        @test_throws MethodError LatSpec._to_tuple("test")  # Should not exist
    end

    # Check Constructors
    # Right hand side should always call the base cosntructor
    @testset verbose=true "DataPoint    " begin
        @test (local D = DataPoint(); all(isnan.([D.val D.stat_err... D.sys_err...])))
        # Check handling of tuples and numbers as input
        @testset "dispatch    " begin
            # uncertainties are numbers
            @test DataPoint(1.0, 2.0, 3.0) == DataPoint(1.0,(2.0,2.0),(3.0,3.0))
            # uncertainties are tuple+number
            @test DataPoint(1.0, (2.0,3.0), 4.0) == DataPoint(1.0,(2.0,3.0),(4.0,4.0))
            # uncertainties are tuples
            @test DataPoint(1.0, (2.0,3.0), (4.0,5.0)) == DataPoint(1.0,(2.0,3.0),(4.0,5.0))
            # uncertainties are number+tuple
            @test DataPoint(1.0, 2.0, (3.0,4.0)) == DataPoint(1.0,(2.0,2.0),(3.0,4.0))
            # unertainties passed by keywords
            @test DataPoint(1.0, sys = 3.0) == DataPoint(1.0,(0.0,0.0),(3.0,3.0))
            @test DataPoint(1.0, sys = (3.0,4.0)) == DataPoint(1.0,(0.0,0.0),(3.0,4.0))
            @test DataPoint(1.0, stat = (3.0,4.0)) == DataPoint(1.0,(3.0,4.0),(0.0,0.0))
            @test DataPoint(1.0, stat = 3.0) == DataPoint(1.0,(3.0,3.0),(0.0,0.0))
        end
        # Check promotion of datatypes
        @testset "promotes    " begin
            # promote to highest (float)
            @test DataPoint(1.0, 2, 3.0) == DataPoint(1.0,(2.0,2.0),(3.0,3.0))
            @test DataPoint(1.0, (2,3.0), 4.0) == DataPoint(1.0,(2.0,3.0),(4.0,4.0))
            @test DataPoint(1.0, 2.0, 3) == DataPoint(1.0,(2.0,2.0),(3.0,3.0))
            @test DataPoint(1.0, 2.0, (3.0,4)) == DataPoint(1.0,(2.0,2.0),(3.0,4.0))
            # promote to highest (Im)
            @test DataPoint(1.0 - im, 2, 3.0) == DataPoint(1.0 - im,(2.0 + im * 0.0,2.0 + im * 0.0),(3.0 + im * 0.0,3.0 + im * 0.0))
            @test DataPoint(1.0 - im, (2,3.0), 4.0) == DataPoint(1.0 - im,(2.0 + im * 0.0,3.0 + im * 0.0),(4.0 + im * 0.0,4.0 + im * 0.0))
            @test DataPoint(1.0 - im, 2.0, 3) == DataPoint(1.0 - im,(2.0 + im * 0.0,2.0 + im * 0.0),(3.0 + im * 0.0,3.0 + im * 0.0))
            @test DataPoint(1.0 - im, 2.0, (3.0,4)) == DataPoint(1.0 - im,(2.0 + im * 0.0,2.0 + im * 0.0),(3.0 + im * 0.0,4.0 + im * 0.0))
        end
        # Check correct handling when value is of type integer
        # cast all numbers to integer and ceil errors
        @testset "integers    " begin
            # unchanged for all integers
            @test DataPoint(1,2,3) == DataPoint(1,(2,2),(3,3))
            # mixtures (old promotion)
            # @test DataPoint(1,2.5,3) == DataPoint(1,(3,3),(3,3))
            # @test DataPoint(1,(2.5,3.5),3) == DataPoint(1,(3,4),(3,3))
            # @test DataPoint(1,2.5,(3,4)) == DataPoint(1,(3,3),(3,4))
            # @test DataPoint(1,2.5,(3.1,4)) == DataPoint(1,(3,3),(4,4))
            # mixtures (new promotion)
            @test DataPoint(1,2.5,3) == DataPoint(1.0,(2.5,2.5),(3.0,3.0))
            @test DataPoint(1,(2.5,3.5),3) == DataPoint(1.0,(2.5,3.5),(3.0,3.0))
            @test DataPoint(1,2.5,(3,4)) == DataPoint(1.0,(2.5,2.5),(3.0,4.0))
            @test DataPoint(1,2.5,(3.1,4)) == DataPoint(1,(2.5,2.5),(3.1,4.0))
            # unertainties passed by keywords (old promotion)
            # @test DataPoint(1, sys = 2.5) == DataPoint(1,(0,0),(3,3))
            # @test DataPoint(1, stat = 2.5) == DataPoint(1,(3,3),(0,0))
            # unertainties passed by keywords (new promotion)
            @test DataPoint(1, sys = 2.5) == DataPoint(1.0,(0.0,0.0),(2.5,2.5))
            @test DataPoint(1, stat = 2.5) == DataPoint(1.0,(2.5,2.5),(0.0,0.0))
        end
        @testset "conversion  " begin
            # Usual conversions
            @test convert(DataPoint{Float32},DataPoint(1.0,2.0,3.0)) == DataPoint(1.0f0,(2.0f0,2.0f0),(3.0f0,3.0f0))
            @test convert(DataPoint{Float64},DataPoint(1,2,3)) == DataPoint(1.0,(2.0,2.0),(3.0,3.0))
            # to integer conversion
            @test convert(DataPoint{Int64},DataPoint(1.0,2.0,3.0)) == DataPoint(1,(2,2),(3,3))
            # to integer conversion with rounding
            @test convert(DataPoint{Int64},DataPoint(1.0,2.5,3.5)) == DataPoint(1,(3,3),(4,4))
            # to integer conversion only works if val is convertible to int
            @test_throws InexactError convert(DataPoint{Int64},DataPoint(1.5,2.5,3.5))
        end
        @testset "errors      " begin
            @test_throws AssertionError DataPoint(1,-2.5,3)
            @test_throws MethodError DataPoint("test",-2.5,3)
        end
    end
    # Check Construction via operator
    @testset "Operator: ±  " begin
        # single operator
        @test 1.0 ± 2.0 == DataPoint(1.0, (2.0,2.0), (0.0,0.0))
        @test 1.0 ± (sys = 2.0,) == DataPoint(1.0, (0.0,0.0), (2.0,2.0))
        @test 1.0 ± (stat = 2.0,) == DataPoint(1.0, (2.0,2.0), (0.0,0.0))
        @test 1.0 ± (stat = 2.0, sys = 3.0) == DataPoint(1.0, (2.0,2.0), (3.0,3.0))
        @test 1.0 ± (other = 2.0,) == DataPoint(1.0, (0.0,0.0), (0.0,0.0))
        # double operator
        @test 1.0 ± (2.0,3.0) ± (4.0,5.0) == DataPoint(1.0,(2.0,3.0),(4.0,5.0))
        # old integer cast
        # @test 15 ± (1, 2.5) ± pi == DataPoint(15, (1,3), (4, 4))
        # new integer cast
        @test 15 ± (1, 2.5) ± pi == DataPoint(15.0, (1.0,2.5), (pi, pi))
    end
    @testset "Helper functions" begin
        x = DataPoint(1.0, 0.1, 0.2)
        @test value(x) == 1.0
        @test staterr(x) == (0.1,0.1)
        @test syserr(x) == (0.2,0.2)
        x = DataPoint(1.0, (0.1,0.12), (0.2,0.23))
        @test value(x) == 1.0
        @test staterr(x) == (0.1,0.12)
        @test syserr(x) == (0.2,0.23)
    end
    @testset "Printing     " begin
        io = IOBuffer()
        D1 = 1.0 ± (2, 3)
        D2 = 1.0 ± (2, 3) ± (4, 5)
        D3 = DataPoint(1.0, sys = (4,5))

        print(io, D1)
        @test String(take!(io)) == "DataPoint{Float64}:\nValue → 1.0\tStatistic unc. → (2.0, -3.0)\tSystematic unc. → (0.0, -0.0)"
        print(io, D2)
        @test String(take!(io)) == "DataPoint{Float64}:\nValue → 1.0\tStatistic unc. → (2.0, -3.0)\tSystematic unc. → (4.0, -5.0)"
        print(io, D3)
        @test String(take!(io)) == "DataPoint{Float64}:\nValue → 1.0\tStatistic unc. → (0.0, -0.0)\tSystematic unc. → (4.0, -5.0)"

        print(IOContext(io, :compact => true), D1)
        @test String(take!(io)) == "1.0 ± (2.0, 3.0) ± (0.0, 0.0)"
        print(IOContext(io, :compact => true), D2)
        @test String(take!(io)) == "1.0 ± (2.0, 3.0) ± (4.0, 5.0)"
        print(IOContext(io, :compact => true), D3)
        @test String(take!(io)) == "1.0 ± (0.0, 0.0) ± (4.0, 5.0)"

        show(io, MIME("text/plain"), D1)
        @test String(take!(io)) == "DataPoint{Float64}:\nValue → 1.0\tStatistic unc. → (2.0, -3.0)"
        show(io, MIME("text/plain"), D2)
        @test String(take!(io)) == "DataPoint{Float64}:\nValue → 1.0\tStatistic unc. → (2.0, -3.0)\tSystematic unc. → (4.0, -5.0)"
        show(io, MIME("text/plain"), D3)
        @test String(take!(io)) == "DataPoint{Float64}:\nValue → 1.0\tSystematic unc. → (4.0, -5.0)"
    end
end


text:
logistic(x) = 1 ./ (1 + exp.(-x))
@diffrule logistic(x::Number) x (logistic(x) .* (1 .- logistic(x)) .* dy)

softplus(x) = log(exp(x) + 1)
@diffrule softplus(x::Number) x logistic(x) .* dy

function xavier_init(dim_in, dim_out; c=1)
    low = -c * sqrt(6.0 / (dim_in + dim_out))
    high = c * sqrt(6.0 / (dim_in + dim_out))
    return rand(Uniform(low, high), dim_in, dim_out)
end


text:
using IntervalWavelets
using Base.Test

println("Testing boundary scaling functions...")

#=
The boundary wavelets are mutually orthonormal and orthogonal to the interior wavelets.
=#

EPS = 1e-3
R = 11

# ------------------------------------------------------------
# Boundary to boundary

for side in ['L'; 'R'], p in 2:8
	#@show side, p
	x, Y = DaubScaling(p, side, R)

	F = bfilter(p, side)
	S = support(F)

	for k in 1:p
		#@show IntervalWavelets.l2norm(x, Y[:,k])
		@test_approx_eq_eps IntervalWavelets.l2norm(x, Y[:,k]) 1.0 10*EPS

		for l in k+1:p
			intersect_supp = intersect( support(F,k-1), support(F,l-1) )
			xx = dyadic_rationals(intersect_supp, R)
			startx = x2index( left(intersect_supp), S, R )
			stopx = x2index( right(intersect_supp), S, R )

			#@show IntervalWavelets.inner(xx, Y[startx:stopx,k], Y[startx:stopx,l])
			@test_approx_eq_eps IntervalWavelets.inner(xx, Y[startx:stopx,k], Y[startx:stopx,l]) 0.0 EPS
		end
	end
end

# ------------------------------------------------------------
# Boundary to interior
# The interior function are translated such that there is a non-empty
# overlap with support of the boundary function.

EPS = 1e-5

for side in ['L'; 'R'], p in 2:8
	C = ifilter(p,true)
	interior = DaubScaling(C,R)
	suppC = support(C)
	int_index = [0:length(interior)-1;]

	B = bfilter(p, side)
	bound = DaubScaling(B, C, R)'

	# Boundary functions
	for k in 0:p-1
		# Translations of interior functions
		suppB = support(B,k)
		translation_index = (side == 'L' ? (p:2*p+k-2) : -2*p-k+1:-p-1)

		for l in translation_index
			suppI = intersect( suppC+l, suppB )
			x = dyadic_rationals(suppI, R)

			startB = x2index( left(suppI), suppB, R )
			stopB = x2index( right(suppI), suppB, R )

			startC = x2index( left(suppI), suppC+l, R )
			stopC = x2index( right(suppI), suppC+l, R )

			#@show IntervalWavelets.inner(x, interior[startC:stopC], bound[startB:stopB,k+1])
			@test_approx_eq_eps IntervalWavelets.inner(x, interior[startC:stopC], bound[startB:stopB,k+1]) 0.0 EPS
		end
	end
end


text:
module SubVectorsTests
using Gridap
using Gridap.Arrays

a = collect(1:10)
pini = 1
pend = 5
b = SubVector(a,pini,pend)
test_array(b,collect(pini:pend))

pini = 7
pend = 9
b = SubVector(a,pini,pend)
test_array(b,collect(pini:pend))

inds = Int32[1,4,3]
a = rand(10)
sa = Gridap.Arrays.SubVector(a,1,5)
vsa = view(sa,inds)

b = rand(10)
vb = view(b,inds)

x = Base.broadcasted(+,vb,vsa)
Base.materialize!(vb,x)

end # module


text:
using Test
using SparseArrays
# using Revise
using CutCellDG
include("../useful_routines.jl")

@test CutCellDG.node_to_dof_id(9,1,2) == 17
@test CutCellDG.node_to_dof_id(5,2,2) == 10

nodeids = [14,15,16,19,20,21,24,25,26]
edofs = CutCellDG.element_dofs(nodeids,2)
testedofs = [27,28,29,30,31,32,37,38,39,40,41,42,47,48,49,50,51,52]
@test allequal(edofs,testedofs)

rows,cols = CutCellDG.element_dofs_to_operator_dofs(1:4,[5,6,7])
testrows = vcat(1:4,1:4,1:4)
@test allequal(testrows,rows)
testcols = vcat([5,5,5,5],[6,6,6,6],[7,7,7,7])
@test allequal(testcols,cols)


text:
using AbstractAlgebra

@testset "ff f4 noreduce degrevlex" begin

R, (x, y) = PolynomialRing(GF(2^31 - 1), ["x", "y"], ordering=:degrevlex)

fs = [
    x + y^2,
    x*y - y^2
]
gb = Groebner.groebner(fs, reduced=false)
# println(gb)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))
#=
    ans is
    y^2 + x,
    x^2 + x,
    x*y + x
=#

fs = [
    x + y,
    x^2 + y
]
gb = Groebner.groebner(fs, reduced=false)
# println(gb)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

fs = [
    x + y,
    x^2 + y
]
gb = Groebner.groebner(fs, reduced=false)
# println(gb)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

fs = [
    x^2 + 5,
    2y^2 + 3
]
gb = Groebner.groebner(fs, reduced=false)
# println(gb)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

fs = [
    y,
    x*y + x
]
gb = Groebner.groebner(fs, reduced=false)
# println(gb)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

fs = [
    5y,
    4x*y + 8x
]
gb = Groebner.groebner(fs, reduced=false)
# println(gb)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

fs = [
    y^2 + x,
    x^2*y + y,
    y + 1
]
gb = Groebner.groebner(fs, reduced=false)
# println(gb)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))
@test gb == [R(1)]

fs = [
    x*y + y,
    y - 5
]
gb = Groebner.groebner(fs, reduced=false)
# println(gb)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

root = Groebner.change_ordering(Groebner.rootn(3, ground=GF(2^31 - 1)), :degrevlex)
gb = Groebner.groebner(root, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

root = Groebner.change_ordering(Groebner.rootn(4, ground=GF(2^31 - 1)), :degrevlex)
gb = Groebner.groebner(root, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

for proot in Groebner.Combinatorics.permutations(root)
    gb = Groebner.groebner(proot, reduced=false)
    @test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))
end

root = Groebner.change_ordering(Groebner.rootn(8, ground=GF(2^31 - 1)), :degrevlex)
gb = Groebner.groebner(root, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

noon = Groebner.change_ordering(Groebner.noonn(3, ground=GF(2^31 - 1)), :degrevlex)
gb = Groebner.groebner(noon, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

for pnoon in Groebner.Combinatorics.permutations(noon)
    gb = Groebner.groebner(pnoon, reduced=false)
    @test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))
end

noon = Groebner.change_ordering(Groebner.noonn(4, ground=GF(2^31 - 1)), :degrevlex)
gb = Groebner.groebner(noon, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

eco = Groebner.change_ordering(Groebner.eco5(), :degrevlex)
gb = Groebner.groebner(eco, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

ku = Groebner.change_ordering(Groebner.ku10(ground=GF(2^31 - 1)), :degrevlex)
gb = Groebner.groebner(ku, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

end

@testset "ff f4 noreduce lex" begin

R, (x, y) = PolynomialRing(GF(2^31 - 1), ["x", "y"], ordering=:lex)

fs = [
    x + y^2,
    x*y - y^2
]
gb = Groebner.groebner(fs, reduced=false)
@test gb ≂ [x + y^2, x*y + 2147483646*y^2, y^3 + y^2]

fs = [
    x + y,
    x^2 + y
]
gb = Groebner.groebner(fs, reduced=false)
@test gb ≂ [x + y, x^2 + y, y^2 + y]

fs = [
    y,
    x*y + x
]
gb = Groebner.groebner(fs, reduced=false)
@test gb ≂ [x, y]

fs = [
    x^2 + 5,
    2y^2 + 3
]
gb = Groebner.groebner(fs, reduced=false)
@test gb ≂ [ y^2 + 1073741825, x^2 + 5]

fs = [
    y^2 + x,
    x^2*y + y,
    y + 1
]
gb = Groebner.groebner(fs, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))
@test gb == [R(1)]

fs = [
    x^2*y^2,
    2*x*y^2 + 3*x*y
]
gb = Groebner.groebner(fs, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))
@test gb ≂ [ x*y^2 + 1073741825*x*y, x^2*y]

root = Groebner.rootn(3, ground=GF(2^31 - 1))
gb = Groebner.groebner(root, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

root = Groebner.rootn(4, ground=GF(2^31 - 1))
gb = Groebner.groebner(root, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

#=
noon = Groebner.noon3(ground=GF(2^31 - 1))
gb = Groebner.groebner(noon, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))
=#

end

@testset "ff f4 noreduce deglex" begin

R, (x, y) = PolynomialRing(GF(2^31 - 1), ["x", "y"], ordering=:deglex)

fs = [
    x + y^2,
    x*y - y^2
]
gb = Groebner.groebner(fs, reduced=true)
@test Groebner._isgroebner_reference(gb, initial_gens=fs)

fs = [
    x + y,
    x^2 + y
]
gb = Groebner.groebner(fs, reduced=true)
@test Groebner._isgroebner_reference(gb, initial_gens=fs)

fs = [
    y,
    x*y + x
]
gb = Groebner.groebner(fs, reduced=false)
@test Groebner._isgroebner_reference(gb, initial_gens=fs)

fs = [
    x^2 + 5,
    2y^2 + 3
]
gb = Groebner.groebner(fs, reduced=false)
@test Groebner._isgroebner_reference(gb, initial_gens=fs)

fs = [
    y^2 + x,
    x^2*y + y,
    y + 1
]
gb = Groebner.groebner(fs, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))
@test Groebner._isgroebner_reference(gb, initial_gens=fs)

fs = [
    x^2*y^2,
    2*x*y^2 + 3*x*y
]
gb = Groebner.groebner(fs, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))
@test Groebner._isgroebner_reference(gb, initial_gens=fs)

root = Groebner.change_ordering(Groebner.rootn(3, ground=GF(2^31 - 1)), :deglex)
gb = Groebner.groebner(root, reduced=true)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

root = Groebner.change_ordering(Groebner.rootn(4, ground=GF(2^31 - 1)), :deglex)
gb = Groebner.groebner(root, reduced=true)
@test Groebner._isgroebner_reference(gb, initial_gens=root)

noon = Groebner.change_ordering(Groebner.noonn(3, ground=GF(2^31 - 1)), :deglex)
gb = Groebner.groebner(noon, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

noon = Groebner.change_ordering(Groebner.noonn(4, ground=GF(2^31 - 1)), :deglex)
gb = Groebner.groebner(noon, reduced=false)
@test Groebner._isgroebner_reference(Groebner._reducegb_reference(gb))

end

@testset "ff f4 corner cases" begin
    R, (x, y) = PolynomialRing(GF(2^31 - 1), ["x", "y"], ordering=:degrevlex)

    fs = [x, x, x]
    @test Groebner.groebner(fs) ≂ [x]

    fs = [x, x, y, y, y, x, x, y]
    @test Groebner.groebner(fs) ≂ [x, y]

    fs = [R(1)]
    @test Groebner.groebner(fs) ≂ [R(1)]

    fs = [x^2 + y, y*x - 1, R(1), y^4]
    @test Groebner.groebner(fs) ≂ [R(1)]

    fs = [2x, 3y, 4x + 5y]
    @test Groebner.groebner(fs) ≂ [x, y]

    R, (x1, x2) = PolynomialRing(GF(2^31 - 1), ["x1", "x2"], ordering=:degrevlex)

    fs = [1*x1^2*x2^2 + 2*x1^2*x2,
        1*x1^2*x2^2 + 3*x1^2*x2 + 5*x1*x2^2,
        1*x1^2*x2^2]
    gb = Groebner.groebner(fs)
    @test gb ≂ [x1*x2^2, x1^2*x2]

    fs = [x1*x2^2 + x1*x2, x1^2*x2^2 + 2*x1*x2, 2*x1^2*x2^2 + x1*x2]
    gb = Groebner.groebner(fs)
    @test gb ≂ [x1*x2]

end


text:
type Hedge <: BanditAlgorithm
  eta::Float64
  counts::Vector{Int64}
  values::Vector{Float64}
end

function Hedge(n_arms::Int64, eta::Float64)
  Hedge(eta, zeros(Int64, n_arms), zeros(n_arms))
end

function initialize(algo::Hedge, n_arms::Int64)
  algo.counts = zeros(Int64, n_arms)
  algo.values = zeros(n_arms)
end

function initialize(algo::Hedge, arms::Vector)
  initialize(algo, length(arms))
end

function select_arm(algo::Hedge)
  z = sum(exp(algo.values * algo.eta))
  probs = exp(algo.values * algo.eta) ./ z
  return rand(Categorical(probs))
end

function update(algo::Hedge, chosen_arm::Int64, reward::Real)
  algo.counts[chosen_arm] += 1

  value = algo.values[chosen_arm]
  algo.values[chosen_arm] = value + reward
end

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
"""
Add variables to the OptimizationContainer for any component.
"""
function add_variables!(
    optimization_container::OptimizationContainer,
    ::Type{T},
    devices::Union{Vector{U}, IS.FlattenIteratorWrapper{U}},
    formulation::Union{AbstractDeviceFormulation, AbstractServiceFormulation},
) where {T <: VariableType, U <: PSY.Component}
    add_variable!(optimization_container, T(), devices, formulation)
end

"""
Add variables to the OptimizationContainer for a service.
"""
function add_variables!(
    optimization_container::OptimizationContainer,
    ::Type{T},
    service::U,
    devices::Vector{V},
    formulation::AbstractReservesFormulation,
) where {T <: VariableType, U <: PSY.Reserve, V <: PSY.Device}
    add_variable!(optimization_container, T(), devices, service, formulation)
end

@doc raw"""
Adds a variable to the optimization model and to the affine expressions contained
in the optimization_container model according to the specified sign. Based on the inputs, the variable can
be specified as binary.

# Bounds

``` lb_value_function <= varstart[name, t] <= ub_value_function ```

If binary = true:

``` varstart[name, t] in {0,1} ```

# LaTeX

``  lb \ge x^{device}_t \le ub \forall t ``

``  x^{device}_t \in {0,1} \forall t iff \text{binary = true}``

# Arguments
* optimization_container::OptimizationContainer : the optimization_container model built in PowerSimulations
* devices : Vector or Iterator with the devices
* var_name::Symbol : Base Name for the variable
* binary::Bool : Select if the variable is binary
* expression_name::Symbol : Expression_name name stored in optimization_container.expressions to add the variable
* sign::Float64 : sign of the addition of the variable to the expression_name. Default Value is 1.0

# Accepted Keyword Arguments
* ub_value : Provides the function over device to obtain the value for a upper_bound
* lb_value : Provides the function over device to obtain the value for a lower_bound. If the variable is meant to be positive define lb = x -> 0.0
* initial_value : Provides the function over device to obtain the warm start value

"""
function add_variable!(
    optimization_container::OptimizationContainer,
    variable_type::VariableType,
    devices::U,
    formulation,
) where {U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}}} where {D <: PSY.Component}
    @assert !isempty(devices)
    time_steps = model_time_steps(optimization_container)
    settings = get_settings(optimization_container)
    var_name = make_variable_name(typeof(variable_type), D)
    binary = get_variable_binary(variable_type, D, formulation)
    expression_name = get_variable_expression_name(variable_type, D)
    sign = get_variable_sign(variable_type, D, formulation)

    variable = add_var_container!(
        optimization_container,
        var_name,
        [PSY.get_name(d) for d in devices],
        time_steps,
    )

    for t in time_steps, d in devices
        name = PSY.get_name(d)
        variable[name, t] = JuMP.@variable(
            optimization_container.JuMPmodel,
            base_name = "$(var_name)_{$(name), $(t)}",
            binary = binary
        )

        ub = get_variable_upper_bound(variable_type, d, formulation)
        !(ub === nothing) && JuMP.set_upper_bound(variable[name, t], ub)

        lb = get_variable_lower_bound(variable_type, d, formulation)
        !(lb === nothing) && !binary && JuMP.set_lower_bound(variable[name, t], lb)

        if get_warm_start(settings)
            init = get_variable_initial_value(variable_type, d, formulation)
            !(init === nothing) && JuMP.set_start_value(variable[name, t], init)
        end

        if !((expression_name === nothing))
            bus_number = PSY.get_number(PSY.get_bus(d))
            add_to_expression!(
                get_expression(optimization_container, expression_name),
                bus_number,
                t,
                variable[name, t],
                get_variable_sign(variable_type, eltype(devices), formulation),
            )
        end
    end

    return
end

# TODO: refactor this function when ServiceModel is updated to include service name
function add_variable!(
    optimization_container::OptimizationContainer,
    variable_type::VariableType,
    devices::U,
    service::T,
    formulation::AbstractReservesFormulation,
) where {
    T <: PSY.Service,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
} where {D <: PSY.Component}
    @assert !isempty(devices)
    time_steps = model_time_steps(optimization_container)

    var_name = make_variable_name(PSY.get_name(service), T)
    binary = get_variable_binary(variable_type, T, formulation)
    expression_name = get_variable_expression_name(variable_type, T)
    sign = get_variable_sign(variable_type, T, formulation)

    variable = add_var_container!(
        optimization_container,
        var_name,
        [PSY.get_name(d) for d in devices],
        time_steps,
    )

    for t in time_steps, d in devices
        name = PSY.get_name(d)
        variable[name, t] = JuMP.@variable(
            optimization_container.JuMPmodel,
            base_name = "$(var_name)_{$(name), $(t)}",
            binary = binary
        )

        ub = get_variable_upper_bound(
            variable_type,
            service,
            d,
            optimization_container.settings,
        )
        !(ub === nothing) && JuMP.set_upper_bound(variable[name, t], ub)

        lb = get_variable_lower_bound(
            variable_type,
            service,
            d,
            optimization_container.settings,
        )
        !(lb === nothing) && !binary && JuMP.set_lower_bound(variable[name, t], lb)

        init = get_variable_initial_value(variable_type, d, optimization_container.settings)
        !(init === nothing) && JuMP.set_start_value(variable[name, t], init)

        if !((expression_name === nothing))
            bus_number = PSY.get_number(PSY.get_bus(d))
            add_to_expression!(
                get_expression(optimization_container, expression_name),
                bus_number,
                t,
                variable[name, t],
                get_variable_sign(variable_type, eltype(devices), formulation),
            )
        end
    end

    return
end

@doc raw"""
Adds a bounds to a variable in the optimization model.

# Bounds

``` bounds.min <= varstart[name, t] <= bounds.max  ```

# LaTeX

``  x^{device}_t >= bound^{min;} \forall t ``

``  x^{device}_t <= bound^{max} \forall t ``

# Arguments
* optimization_container::OptimizationContainer : the optimization_container model built in PowerSimulations
* bounds::DeviceRangeConstraintInfo : contains names and vector of min / max
* var_type::AbstractString : type of the variable
* T: type of the device

"""
function set_variable_bounds!(
    optimization_container::OptimizationContainer,
    bounds::Vector{DeviceRangeConstraintInfo},
    var_type::AbstractString,
    ::Type{T},
) where {T <: PSY.Component}
    var = get_variable(optimization_container, var_type, T)
    for t in model_time_steps(optimization_container), bound in bounds
        _var = var[get_component_name(bound), t]
        JuMP.set_upper_bound(_var, bound.limits.max)
        JuMP.set_lower_bound(_var, bound.limits.min)
    end
end


text:
loadleaf!(dst, src, err) = dst
loadleaf!(dst::AbstractArray, src, err) =
  error("Tried to copy $src into an array destination; this is not allowed.")
loadleaf!(dst, src::AbstractArray, err) =
  error("Tried to copy an array to $dst; this is not allowed.")
function loadleaf!(dst::AbstractArray, src::Bool, err)
  if iszero(src)
    dst .= src
  else
    error("Cannot copy boolean parameter == true to non-zero parameter.")
  end
  return dst
end
loadleaf!(dst::Bool, src::AbstractArray, err) = iszero(dst) ? dst :
  error("Cannot copy non-zero parameter to boolean parameter == true.")
function loadleaf!(dst::AbstractArray, src::AbstractArray, err)
  (size(dst) == size(src)) || throw(err)
  copyto!(dst, src)
end

_tie_check(dst::Bool, src::AbstractArray) = iszero(dst) ||
  error("Encountered tied parameter with boolean source at some nodes and non-boolean sources at others.")
_tie_check(dst::AbstractArray, src::Bool) = (iszero(dst) && iszero(src)) ||
  error("Encountered tied parameter with boolean source at some nodes and non-boolean sources at others.")
_tie_check(dst::AbstractArray, src::AbstractArray) = (dst == src) ||
  error("Encountered tied destination parameters with untied and mismatched sources.")
_tie_check(dst, src) = true

_bool_tie_check(dst, src) = true

"""
    loadmodel!(dst, src)

Copy all the parameters (trainable and non-trainable) from `src` into `dst`.

Recursively walks `dst` and `src` together using [`Functors.children`](@ref),
and calling `copyto!` on parameter arrays or throwing an error when there is a mismatch.
Non-array elements (such as activation functions) are not copied and need not match.
Zero bias vectors and `bias=false` are considered equivalent
(see extended help for more details).

# Examples
```julia
julia> dst = Chain(Dense(Flux.ones32(2, 5, tanh)), Dense(2 => 1; bias = [1f0]))
Chain(
  Dense(5 => 2, tanh),                  # 12 parameters
  Dense(2 => 1),                        # 3 parameters
)                   # Total: 4 arrays, 15 parameters, 316 bytes.

julia> dst[1].weight ≈ ones(2, 5)  # by construction
true

julia> src = Chain(Dense(5 => 2, relu), Dense(2 => 1, bias=false));

julia> Flux.loadmodel!(dst, src);

julia> dst[1].weight ≈ ones(2, 5)  # values changed
false

julia> iszero(dst[2].bias)
true
```

# Extended help

Throws an error when:
- `dst` and `src` do not share the same fields (at any level)
- the sizes of leaf nodes are mismatched between `dst` and `src`
- copying non-array values to/from an array parameter
  (except inactive parameters described below)
- `dst` is a "tied" parameter (i.e. refers to another parameter) and
  loaded into multiple times with mismatched source values

Inactive parameters can be encoded by using the boolean value `false` instead of an array.
If `dst == false` and `src` is an all-zero array, no error will be raised (and no values copied);
however, attempting to copy a non-zero array to an inactive parameter will throw an error.
Likewise, copying a `src` value of `false` to any `dst` array is valid,
but copying a `src` value of `true` will error.
"""
function loadmodel!(dst, src; cache = Base.IdSet())
  ldsts, _ = functor(dst)
  lsrcs, _ = functor(src)
  (keys(ldsts) == keys(lsrcs)) ||
    throw(ArgumentError("Tried to load $src into $dst but the structures do not match."))

  err = DimensionMismatch("Tried to load $src into $dst but the parameter sizes do not match.")
  foreach(ldsts, lsrcs) do ldst, lsrc
    if ldst in cache # we already loaded this parameter before
      _tie_check(ldst, lsrc) && return ldst
    elseif Functors.isleaf(ldst) # our first time loading this leaf
      push!(cache, ldst)
      loadleaf!(ldst, lsrc, err)
    else # this isn't a leaf
      loadmodel!(ldst, lsrc; cache = cache)
    end
  end

  return dst
end


text:
module Day7
export day7

function parseLine(line)
    m = match(r"^([\w\s]*) bags contain ([\w\s,]*)\.$", line)
    if occursin("no other bags", m[2])
        return m[1], ()
    end
    m[1], map(split(m[2], ',')) do s
        let matches = match(r"(\d*) ([a-z\s]*) bags?", strip(s))
            parse(Int, matches[1]), matches[2]
        end
    end
end

function readInput(inputfile="inputs/input07.txt")
    data = Dict{String, Any}()
    for line in eachline(inputfile)
        newBag, contains = parseLine(line)
        data[newBag] = contains
    end
    data
end

function day7(inputfile="inputs/input07.txt")
    data = readInput(inputfile)
    (day7Part1(data), day7Part2(data))
end

function day7Part1(data)
    ## invert relationship data[bag] -> contains
    ## need map bag -> contained in
    containedIn = Dict{String, Any}()
    for (bag, contains) in data
        for (amount, containedBag) in contains
            if haskey(containedIn, containedBag)
                push!(containedIn[containedBag], bag)
            else
                containedIn[containedBag] = [bag]
            end
        end
    end

    ## find all connected nodes with "shiny gold"
    found = Set{String}()
    todo = containedIn["shiny gold"]
    while length(todo) > 0
        current = pop!(todo)
        if !(current in found) # have we looked at current already?
            push!(found, current)
            if haskey(containedIn, current) # current might be a root
                append!(todo, containedIn[current])
            end
        end
    end
    length(found)

end

function day7Part2(data)
    countBagsContained(data, "shiny gold")
end

## Little benchmark showed: memoization offers in fact NO speedup :(
countBagsContained(data, bag) = countBagsContained(data, bag, Dict{String, Int}())
function countBagsContained(data, bag, alreadyComputed)
    if !haskey(alreadyComputed, bag)
        amountContained = 0
        for (amount, subBag) in get(data, bag, ())
            amountContained += amount * (1 + countBagsContained(data, subBag, alreadyComputed))
        end
        alreadyComputed[bag] = amountContained
    end
    alreadyComputed[bag]
end

end #module


text:
@withkw @trait immutable PoissonParameter
    rate::Float = 1Hz
end

@withkw @trait type Poisson
    param::PoissonParameter = PoissonParameter()
    N::Int = 100
    randcache::Vector{Float} = rand(N)
    fire::Vector{Bool} = zeros(Bool, N)
    records::Dict = Dict()
end

@replace function integrate!(p::Poisson, param::PoissonParameter, dt::Float)
    prob = rate * dt
    rand!(randcache)
    @inbounds for i = 1:N
        fire[i] = randcache[i] < prob
    end
end


text:
# For systems of rigid bodies, this constructs the right-hand side of the equations
# with the velocity prescribed by the motions in ml
# The state vector is a concatenation of the 3 x 1 vectors for each body
function TimeMarching.r₁(u::Vector{Float64},t::Real,ml::Vector{RigidBodyMotion})
    du = deepcopy(u)
    cnt = 0
    for ib = 1:length(ml)
        _,ċ,_,_,α̇,_ = ml[ib](t)
        du[cnt+1:cnt+3] = [real(ċ),imag(ċ),α̇]
        cnt += 3
    end
    return du
end


text:
using InformationMeasures
using Base.Test

arr1 = rand(1000)
arr2 = rand(1000)
arr3 = rand(1000)

include("testDiscretization.jl")
include("testEntropy.jl")
include("testConditionalEntropy.jl")
include("testMutualInformation.jl")
include("testConditionalMutualInformation.jl")
include("testInteractionInformation.jl")
include("testTotalCorrelation.jl")
include("testPartialInformationDecomposition.jl")


text:
@reexport module MaintenanceScheduling

# Standard library
using Statistics
using Printf
using Dates

using RecipesBase
using Distributions
using Parameters
using StochasticPrograms
using Flux
using HydroModels
using HydroModels: AbstractHydroModel, StochasticHydroModel, River, Plant, Area, Scenario, Q̄, μ, %

import HydroModels: modelindices

using Plots: font, text, Shape

export
    MaintenanceSchedulingData,
    MaintenanceSchedulingScenario,
    RecurrentMaintenanceSchedulingSampler,
    MaintenanceSchedulingModel

include("data.jl")
include("scenarios.jl")
include("indices.jl")
include("model.jl")

end


text:
instructions = string.(split(strip(read("input.txt", String)), "\n"))

struct Bitmask
    andmask::UInt64
    ormask::UInt64
end
Bitmask(raw::String) = Bitmask([parse(Int, replace(raw[8:end], "X" => l), base=2) for l in "10"]...)
mask(number::UInt64, bitmask::Bitmask) = (number & bitmask.andmask) | bitmask.ormask

struct Mutator
    address::UInt64
    value::UInt64
end
Mutator(raw::String) = Mutator(parse.(Int, match(r"\[(\d+)\] = (\d+)", raw).captures)...)

function run(instructions::Array{String})::UInt64
    program = [(startswith(instruction, "mask") ? Bitmask : Mutator)(instruction) for instruction in instructions]
    memory = zeros(UInt64, maximum(m -> m.value, filter(x -> isa(x, Mutator), program)))
    bitmask = nothing

    for instruction in program
        if instruction isa Bitmask
            bitmask = instruction
        elseif instruction isa Mutator
            memory[instruction.address] = mask(instruction.value, bitmask)
        end
    end

    sum(memory)
end

run(program)

function maskedaddresses(address::String, mask::String)
    result = []
    for maskno in 0:(1 << count("X", mask))
        b = string(maskno, base=2)
        push!(result, parse(Int64, address) | b)
    end
    result
end

using SparseArrays

function otherrun(instructions::Array{String})::UInt64
    memory = spzeros(UInt64, 2^36)
    bitmask = nothing

    for instruction in instructions
        if startswith("mask", instruction)
            bitmask = instruction[8:end]
        else
            address, value = parse.(Int, match(r"\[(\d+)\] = (\d+)", raw).captures)
            for maskedaddress in maskedaddresses(address, bitmask)
                memory[maskedaddress] = value
            end
        end
    end

    sum(memory)
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
@testset "1320.minimum-distance-to-type-a-word-using-two-fingers.jl" begin
    @test minimum_distance("HAPPY") == 6
    @test minimum_distance("CAKE") == 3
    @test minimum_distance("NEW") == 3
    @test minimum_distance("YEAR") == 7
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

using FEMBasis
using Test

@testset "FEMBasis.jl" begin
    @testset "test substitution" begin include("test_substitute.jl") end
    @testset "create Vandermonde matrix" begin include("test_vandermonde.jl") end
    @testset "create interpolation polynomials" begin include("test_generate_polynomials.jl") end
    #@testset "test_create_basis" begin include("test_create_basis.jl") end
    @testset "test_lagrange_elements" begin include("test_lagrange_elements.jl") end
    @testset "test_nurbs_elements" begin include("test_nurbs_elements.jl") end
    @testset "test_math" begin include("test_math.jl") end
end


text:
@compat abstract type ParentLinks end

"""
  Indicates that this tree stores parent links explicitly. The implementation
  is responsible for defining the parentind function to expose this
  information.
"""
immutable StoredParents <: ParentLinks; end
immutable ImplicitParents <: ParentLinks; end

parentlinks(::Type) = ImplicitParents()
parentlinks(tree) = parentlinks(typeof(tree))

@compat abstract type SiblingLinks end

"""
  Indicates that this tree stores sibling links explicitly, or can compute them
  quickly (e.g. because the tree has a (small) fixed branching ratio, so the
  current index of a node can be determined by quick linear search). The
  implementation is responsible for defining the relative_state function
  to expose this information.
"""
immutable StoredSiblings <: SiblingLinks; end
immutable ImplicitSiblings <: SiblingLinks; end

siblinglinks(::Type) = ImplicitSiblings()
siblinglinks(tree) = siblinglinks(typeof(tree))

@compat abstract type TreeKind end
immutable RegularTree <: TreeKind; end
immutable IndexedTree <: TreeKind; end

treekind(tree::Type) = RegularTree()
treekind(tree) = treekind(typeof(tree))
children(tree, node, ::RegularTree) = children(node)
children(tree, node, ::IndexedTree) = (tree[y] for y in childindices(tree, node))
children(tree, node) = children(tree, node, treekind(tree))

function rootstate()
end

childindices(tree, node) =
  tree == node ? childindices(tree, rootstate(tree)) :
  error("Must implement childindices(tree, node)")
function childindices()
end

function parentind
end

nodetype(tree) = Any
idxtype(tree) = Int

max_stars_repo_name: femtocleaner/DeprecationBotTest

max_issues_repo_name: femtocleaner/DeprecationBotTest

max_forks_repo_name: femtocleaner/DeprecationBotTest


text:
# Here are all the deprecations

# deprecations from 0.4 on
@deprecate ccm89(x::AbstractArray, r_v::Real = 3.1) ccm89.(x, r_v)
@deprecate od94(x::AbstractArray, r_v::Real = 3.1) od94.(x, r_v)

# deprecations from 0.7 on
@deprecate redden(f::Real, λ::Real, Av::Real; Rv = 3.1, law = CCM89) redden(law, λ, f; Av = Av, Rv = 3.1)
@deprecate redden(f::Quantity, λ::Quantity, Av::Real; Rv = 3.1, law = CCM89) redden(law, λ, f; Av = Av, Rv = 3.1)

@deprecate deredden(f::Real, λ::Real, Av::Real; Rv = 3.1, law = CCM89) deredden(law, λ, f; Av = Av, Rv = 3.1)
@deprecate deredden(f::Quantity, λ::Quantity, Av::Real; Rv = 3.1, law = CCM89) deredden(law, λ, f; Av = Av, Rv = 3.1)

@deprecate ccm89(wave, Rv = 3.1) CCM89(Rv = Rv)(wave)
@deprecate od94(wave, Rv = 3.1) OD94(Rv = Rv)(wave)
@deprecate cal00(wave, Rv = 4.05) CAL00(Rv = Rv)(wave)
@deprecate gcc09(wave, Rv = 3.1) GCC09(Rv = Rv)(wave)
@deprecate vcg04(wave, Rv = 3.1) VCG04(Rv = Rv)(wave)


text:
using FactCheck, Images, Color, FixedPointNumbers, TestImages

facts("OS X reader") do
    context("Autumn leaves") do
        img = testimage("autumn_leaves")
        @fact colorspace(img) --> "RGBA"
        @fact ndims(img) --> 2
        @fact colordim(img) --> 0
        @fact eltype(img) --> AlphaColorValue{RGB{Ufixed16}, Ufixed16}
    end
    context("Camerman") do
        img = testimage("cameraman")
        @fact colorspace(img) --> "Gray"
        @fact ndims(img) --> 2
        @fact colordim(img) --> 0
        @fact eltype(img) --> Gray{Ufixed8}
    end
    context("Earth Apollo") do
        img = testimage("earth_apollo17")
        @fact colorspace(img) --> "RGB4"
        @fact ndims(img) --> 2
        @fact colordim(img) --> 0
        @fact eltype(img) --> RGB4{Ufixed8}
    end
    context("Fabio") do
    img = testimage("fabio")
        @fact colorspace(img) --> "Gray"
        @fact ndims(img) --> 2
        @fact colordim(img) --> 0
        @fact eltype(img) --> Gray{Ufixed8}
    end
    context("House") do
        img = testimage("house")
        @fact colorspace(img) --> "GrayAlpha"
        @fact ndims(img) --> 2
        @fact colordim(img) --> 0
        @fact eltype(img) --> AlphaColorValue{Gray{Ufixed8}, Ufixed8}
    end
    context("Jetplane") do
        img = testimage("jetplane")
        @fact colorspace(img) --> "GrayAlpha"
        @fact ndims(img) --> 2
        @fact colordim(img) --> 0
        @fact eltype(img) --> AlphaColorValue{Gray{Ufixed8}, Ufixed8}
    end
    context("Lighthouse") do
        img = testimage("lighthouse")
        @fact colorspace(img) --> "RGB4"
        @fact ndims(img) --> 2
        @fact colordim(img) --> 0
        @fact eltype(img) --> RGB4{Ufixed8}
    end
    context("Mandrill") do
        img = testimage("mandrill")
        @fact colorspace(img) --> "RGB"
        @fact ndims(img) --> 2
        @fact colordim(img) --> 0
        @fact eltype(img) --> RGB{Ufixed8}
    end
    context("Moonsurface") do
        img = testimage("moonsurface")
        @fact colorspace(img) --> "Gray"
        @fact ndims(img) --> 2
        @fact colordim(img) --> 0
        @fact eltype(img) --> Gray{Ufixed8}
    end
    context("Mountainstream") do
        img = testimage("mountainstream")
        @fact colorspace(img) --> "RGB4"
        @fact ndims(img) --> 2
        @fact colordim(img) --> 0
        @fact eltype(img) --> RGB4{Ufixed8}
    end
end


text:
using DataFrames, Query, Turing, StatsFuns

d = DataFrame(species = ["afarensis", "africanus", "habilis", "boisei", "rudolfensis", "ergaster", "sapiens"],
         brain = [438, 452, 612, 521, 752, 871, 1350],
         mass = [37.0, 35.5, 34.5, 41.5, 55.5, 61.0, 53.5])

# d = d |> @mutate(mass_std = (_.mass - mean(d.mass)),
#          brain_std = _.brain / maximum(d.brain)) |> DataFrame

# TODO: use population standard deviation to be consistent with Rethinking
scale(x) = (x .- mean(x)) ./ std(x, corrected = true)
d.mass_std = scale(d.mass)
d.brain_std = d.brain ./ maximum(d.brain)

@model model7_1(brain_std, mass_std) = begin
    log_σ ~ Normal(0, 1)
    σ = exp(log_σ)
    α ~ Normal(0.5, 1)
    β ~ Normal(0, 10)
    for i ∈ eachindex(brain_std)
        brain_std[i] ~ Normal(α + β * mass_std[i], σ)
    end
end

brain_std = d.brain_std
mass_std = d.mass_std

model = model7_1(brain_std, mass_std)

m7_1 = sample(model, NUTS(), 2000, save_state = true)

describe(m7_1)

loglik =
[logprob"brain_std = brain_std[i] | mass_std = mass_std[i], chain = m7_1" for i ∈ eachindex(brain_std)]

ll = DataFrame(loglik)

function WAIC(ll::DataFrame)
    n_samples, n_obs = size(ll)
    logsums = [logsumexp(col) for col ∈ eachcol(ll)]
    lppd = logsums .- log(n_samples)
    pWAIC = [var(col) for col ∈ eachcol(ll)]
    sum_lppd = sum(lppd)
    sum_pWAIC =sum(pWAIC)
    waic = -2 * (sum_lppd - sum_pWAIC)
    waic_vec = -2 .* (lppd .- pWAIC)
    se = sqrt(n_obs * var(waic_vec))

    return (waic = waic, se = se, lppd = sum_lppd, pwaic = sum_pWAIC)
end

waic, se, lppd, pWAIC = WAIC(ll)


text:
function f(n::Int)
    steps = []
    while !isborder(robot, Sud)
        move!(robot, Sud)
        push!(steps, Nord)
    end
    while !isborder(robot, West)
        move!(robot, West)
        push!(steps, Ost)
    end
    k = 0
    done = false
    while !done
        k = (k + 1) % 2
        for i = 1:n
            j = 0
            while !isborder(robot, Ost)
                if trunc(j / n) % 2 != k
                    putmarker!(robot)
                end
                j += 1
                move!(robot, Ost)
            end
            if trunc(j / n) % 2 != k
                putmarker!(robot)
            end
            while !isborder(robot, West)
                move!(robot, West)
            end
            if !isborder(robot, Nord)
                move!(robot, Nord)
            else
                done = true
                break
            end
        end
    end
    while !isborder(robot, West)
        move!(robot, West)
    end
    while !isborder(robot, Sud)
        move!(robot, Sud)
    end
    for i in steps
        move!(robot, i)
    end
end

f(2)

max_stars_repo_name: TheXlebushek/mirea-progs

max_issues_repo_name: TheXlebushek/mirea-progs

max_forks_repo_name: TheXlebushek/mirea-progs


text:
# syntax: proto2
# COV_EXCL_START
using ProtoBuf
import ProtoBuf.meta

const RpcKindProto = (;[
    Symbol("RPC_BUILTIN") => Int32(0),
    Symbol("RPC_WRITABLE") => Int32(1),
    Symbol("RPC_PROTOCOL_BUFFER") => Int32(2),
]...)

mutable struct RPCTraceInfoProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RPCTraceInfoProto(; kwargs...)
        obj = new(meta(RPCTraceInfoProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RPCTraceInfoProto
const __meta_RPCTraceInfoProto = Ref{ProtoMeta}()
function meta(::Type{RPCTraceInfoProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RPCTraceInfoProto)
            __meta_RPCTraceInfoProto[] = target = ProtoMeta(RPCTraceInfoProto)
            allflds = Pair{Symbol,Union{Type,String}}[:traceId => Int64, :parentId => Int64]
            meta(target, RPCTraceInfoProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RPCTraceInfoProto[]
    end
end
function Base.getproperty(obj::RPCTraceInfoProto, name::Symbol)
    if name === :traceId
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :parentId
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

mutable struct RPCCallerContextProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RPCCallerContextProto(; kwargs...)
        obj = new(meta(RPCCallerContextProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RPCCallerContextProto
const __meta_RPCCallerContextProto = Ref{ProtoMeta}()
function meta(::Type{RPCCallerContextProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RPCCallerContextProto)
            __meta_RPCCallerContextProto[] = target = ProtoMeta(RPCCallerContextProto)
            req = Symbol[:context]
            allflds = Pair{Symbol,Union{Type,String}}[:context => AbstractString, :signature => Array{UInt8,1}]
            meta(target, RPCCallerContextProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RPCCallerContextProto[]
    end
end
function Base.getproperty(obj::RPCCallerContextProto, name::Symbol)
    if name === :context
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :signature
        return (obj.__protobuf_jl_internal_values[name])::Array{UInt8,1}
    else
        getfield(obj, name)
    end
end

const RpcRequestHeaderProto_OperationProto = (;[
    Symbol("RPC_FINAL_PACKET") => Int32(0),
    Symbol("RPC_CONTINUATION_PACKET") => Int32(1),
    Symbol("RPC_CLOSE_CONNECTION") => Int32(2),
]...)

mutable struct RpcRequestHeaderProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RpcRequestHeaderProto(; kwargs...)
        obj = new(meta(RpcRequestHeaderProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RpcRequestHeaderProto
const __meta_RpcRequestHeaderProto = Ref{ProtoMeta}()
function meta(::Type{RpcRequestHeaderProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RpcRequestHeaderProto)
            __meta_RpcRequestHeaderProto[] = target = ProtoMeta(RpcRequestHeaderProto)
            req = Symbol[:callId,:clientId]
            val = Dict{Symbol,Any}(:retryCount => -1)
            wtype = Dict(:callId => :sint32, :retryCount => :sint32)
            allflds = Pair{Symbol,Union{Type,String}}[:rpcKind => Int32, :rpcOp => Int32, :callId => Int32, :clientId => Array{UInt8,1}, :retryCount => Int32, :traceInfo => RPCTraceInfoProto, :callerContext => RPCCallerContextProto, :stateId => Int64]
            meta(target, RpcRequestHeaderProto, allflds, req, ProtoBuf.DEF_FNUM, val, ProtoBuf.DEF_PACK, wtype, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RpcRequestHeaderProto[]
    end
end
function Base.getproperty(obj::RpcRequestHeaderProto, name::Symbol)
    if name === :rpcKind
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :rpcOp
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :callId
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :clientId
        return (obj.__protobuf_jl_internal_values[name])::Array{UInt8,1}
    elseif name === :retryCount
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :traceInfo
        return (obj.__protobuf_jl_internal_values[name])::RPCTraceInfoProto
    elseif name === :callerContext
        return (obj.__protobuf_jl_internal_values[name])::RPCCallerContextProto
    elseif name === :stateId
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

const RpcResponseHeaderProto_RpcStatusProto = (;[
    Symbol("SUCCESS") => Int32(0),
    Symbol("ERROR") => Int32(1),
    Symbol("FATAL") => Int32(2),
]...)

const RpcResponseHeaderProto_RpcErrorCodeProto = (;[
    Symbol("ERROR_APPLICATION") => Int32(1),
    Symbol("ERROR_NO_SUCH_METHOD") => Int32(2),
    Symbol("ERROR_NO_SUCH_PROTOCOL") => Int32(3),
    Symbol("ERROR_RPC_SERVER") => Int32(4),
    Symbol("ERROR_SERIALIZING_RESPONSE") => Int32(5),
    Symbol("ERROR_RPC_VERSION_MISMATCH") => Int32(6),
    Symbol("FATAL_UNKNOWN") => Int32(10),
    Symbol("FATAL_UNSUPPORTED_SERIALIZATION") => Int32(11),
    Symbol("FATAL_INVALID_RPC_HEADER") => Int32(12),
    Symbol("FATAL_DESERIALIZING_REQUEST") => Int32(13),
    Symbol("FATAL_VERSION_MISMATCH") => Int32(14),
    Symbol("FATAL_UNAUTHORIZED") => Int32(15),
]...)

mutable struct RpcResponseHeaderProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RpcResponseHeaderProto(; kwargs...)
        obj = new(meta(RpcResponseHeaderProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RpcResponseHeaderProto
const __meta_RpcResponseHeaderProto = Ref{ProtoMeta}()
function meta(::Type{RpcResponseHeaderProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RpcResponseHeaderProto)
            __meta_RpcResponseHeaderProto[] = target = ProtoMeta(RpcResponseHeaderProto)
            req = Symbol[:callId,:status]
            val = Dict{Symbol,Any}(:retryCount => -1)
            wtype = Dict(:retryCount => :sint32)
            allflds = Pair{Symbol,Union{Type,String}}[:callId => UInt32, :status => Int32, :serverIpcVersionNum => UInt32, :exceptionClassName => AbstractString, :errorMsg => AbstractString, :errorDetail => Int32, :clientId => Array{UInt8,1}, :retryCount => Int32, :stateId => Int64]
            meta(target, RpcResponseHeaderProto, allflds, req, ProtoBuf.DEF_FNUM, val, ProtoBuf.DEF_PACK, wtype, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RpcResponseHeaderProto[]
    end
end
function Base.getproperty(obj::RpcResponseHeaderProto, name::Symbol)
    if name === :callId
        return (obj.__protobuf_jl_internal_values[name])::UInt32
    elseif name === :status
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :serverIpcVersionNum
        return (obj.__protobuf_jl_internal_values[name])::UInt32
    elseif name === :exceptionClassName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :errorMsg
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :errorDetail
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :clientId
        return (obj.__protobuf_jl_internal_values[name])::Array{UInt8,1}
    elseif name === :retryCount
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :stateId
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

const RpcSaslProto_SaslState = (;[
    Symbol("SUCCESS") => Int32(0),
    Symbol("NEGOTIATE") => Int32(1),
    Symbol("INITIATE") => Int32(2),
    Symbol("CHALLENGE") => Int32(3),
    Symbol("RESPONSE") => Int32(4),
    Symbol("WRAP") => Int32(5),
]...)

mutable struct RpcSaslProto_SaslAuth <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RpcSaslProto_SaslAuth(; kwargs...)
        obj = new(meta(RpcSaslProto_SaslAuth), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RpcSaslProto_SaslAuth
const __meta_RpcSaslProto_SaslAuth = Ref{ProtoMeta}()
function meta(::Type{RpcSaslProto_SaslAuth})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RpcSaslProto_SaslAuth)
            __meta_RpcSaslProto_SaslAuth[] = target = ProtoMeta(RpcSaslProto_SaslAuth)
            req = Symbol[:method,:mechanism]
            allflds = Pair{Symbol,Union{Type,String}}[:method => AbstractString, :mechanism => AbstractString, :protocol => AbstractString, :serverId => AbstractString, :challenge => Array{UInt8,1}]
            meta(target, RpcSaslProto_SaslAuth, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RpcSaslProto_SaslAuth[]
    end
end
function Base.getproperty(obj::RpcSaslProto_SaslAuth, name::Symbol)
    if name === :method
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :mechanism
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :protocol
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :serverId
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :challenge
        return (obj.__protobuf_jl_internal_values[name])::Array{UInt8,1}
    else
        getfield(obj, name)
    end
end

mutable struct RpcSaslProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RpcSaslProto(; kwargs...)
        obj = new(meta(RpcSaslProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct RpcSaslProto
const __meta_RpcSaslProto = Ref{ProtoMeta}()
function meta(::Type{RpcSaslProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RpcSaslProto)
            __meta_RpcSaslProto[] = target = ProtoMeta(RpcSaslProto)
            req = Symbol[:state]
            allflds = Pair{Symbol,Union{Type,String}}[:version => UInt32, :state => Int32, :token => Array{UInt8,1}, :auths => Base.Vector{RpcSaslProto_SaslAuth}]
            meta(target, RpcSaslProto, allflds, req, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RpcSaslProto[]
    end
end
function Base.getproperty(obj::RpcSaslProto, name::Symbol)
    if name === :version
        return (obj.__protobuf_jl_internal_values[name])::UInt32
    elseif name === :state
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :token
        return (obj.__protobuf_jl_internal_values[name])::Array{UInt8,1}
    elseif name === :auths
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{RpcSaslProto_SaslAuth}
    else
        getfield(obj, name)
    end
end

export RpcKindProto, RPCTraceInfoProto, RPCCallerContextProto, RpcRequestHeaderProto_OperationProto, RpcRequestHeaderProto, RpcResponseHeaderProto_RpcStatusProto, RpcResponseHeaderProto_RpcErrorCodeProto, RpcResponseHeaderProto, RpcSaslProto_SaslState, RpcSaslProto_SaslAuth, RpcSaslProto
# COV_EXCL_STOP


text:
# !! auto-generated code, do not edit !!
module shape_msgs
using RobotOSData.Messages
using RobotOSData.StdMsgs
using RobotOSData.CommonMsgs
struct MeshTriangle <: Readable
    vertex_indices::SVector{3, UInt32}
end
struct SolidPrimitive <: Readable
    type_::UInt8
    dimensions::Vector{Float64}
end
struct Plane <: Readable
    coef::SVector{4, Float64}
end
struct Mesh <: Readable
    triangles::Vector{MeshTriangle}
    vertices::Vector{geometry_msgs.Point}
end
end