text:
"""
   Lanczos bidiagonalization with partial reorthogonalization, A can be matrix or
forward operator. k is the steps. r0 is the starting vector. the output U is a m*k
matrix, Bk is a k*k bidiagonal matrix and V is a n*k matrix and the residue vector
r0 make A*V_{k} = U_{k} * B_{k} + u_{k+1} * beta_{k+1} * e_{k}' and r0 = u_{k+1} * beta[k+1]

the call of the functon is
(U, Bk, V, r0, ierr, workload) = lanbpro(A, k, r0)

workload =[numReoU, numIptU; numReoV, numIptV]

numReoU: number of orthgonalization for U
numIptU: number of inner product for U
numReoV: number of orthgonalization for V
numIptV: number of inner product for V

"""
function lanbpro(A, k::Ti; m=[], n=[], params=[],
                 delta=[], Anorm=[], Uold=[], Bold=[], Vold=[], r0=[],
                 examine_flag=false, elr_flag=true) where {Ti<:Int64}

    # decide A is a function or operator
    matrix_flag = isa(A, Matrix)

    # A is a matrix
    if matrix_flag
       Tv = eltype(A)
       (m, n) = size(A)

    # A is function
    else
       if m == [] || n==[]
          error("must specify the size of operator A")
       end
       etype = Float64
    end

    # initialize some parameters
    smleps = eps(etype)

    # relaxiation coefficient for the inner product after reorthogonalization
    m2 = etype(3./2.)
    n2 = etype(3./2.)

    # desired level of orthogonaloty (the defult value)
    delta = delta == [] ? sqrt(smleps/k) : delta

    # critia for selecting subset of vectors
    eta   = etype(smleps^(3./4.) / sqrt(k))

    # the estimated error of matrix multiplication
    eps1  = etype(sqrt(maximum([m n])) * eps(etype) / 2.0)

    # force orthogonalization at every step (still select subset of vectors to orthogonal aginst)
    forceRO_flag = delta == 0 ? true : false

    # parameter for Gram-schmit method
    gamma   = etype(1.0 / sqrt(2.0))
    gs_iflag = 2 # modified iterated Gram-schmidt

    # initialize the estimate of the maximum of the eigenvalue of A
    if Anorm == []
       Anorm         = zero(etype)
       estAnorm_flag = true

       # relaxiation coefficient for the estimation of ||A||_2
       fudge         = etype(1.01)

       if Bold != []
          error("must provide Anorm when Bold is not empty")
       end
    end

    # the initial input
    if r0 != []
       p = copy(r0)
    else
       p = rand(etype, m) - etype(0.5)
    end

    # initialize output
    U = zeros(etype, m, k)
    B = zeros(etype, k, k)
    V = zeros(etype, n, k)
    beta  = zeros(etype, k+1)
    alpha = zeros(etype, k  )
    indexSet = []
    ierr     = 0

    # start with the first vector, no existed decompositions
    if Uold == [] && Bold == [] && Vold == []
       beta[1] = norm(p)

       iptU = zeros(etype, k+1); iptU[1] = one(etype);
       iptV = zeros(etype, k  ); iptV[1] = one(etype);

       iptUmax = zeros(etype, k)
       iptVmax = zeros(etype, k)

       # count the number of inner product
       numReoU = 0; numIptU = 0;
       numReoV = 0; numIptV = 0;

       # the flag to remove coupling
       coupling_flag = false

       j0 = 1

    elseif Uold != [] && Bold != [] && Vold != []

       # input with the existing factorizations
       j0 = size(Uold, 2)
       U[:, 1:j0] = Uold
       V[:, 1:j0] = Vold

       alpha[1:j0]    = diag(Bold)
       if j0 > 1
          beta[2:j0]  = diag(Bold, -1)
       end
       beta[j0+1] = norm(p)

       # check for convergence
       if j0 < k && beta[j0+1] * delta < Anorm * smleps
          forceRO_flag = true
          ierr         = j0
       end
       indexSet = collect(1:j0)
       (p, beta[j0+1], rr) = reorthogonalization(U, p, beta[j0+1], indexSet,
                                                 gamma=gamma, iflag=gs_flag)

       # statistic
       numReoU = 1; numIptU = rr * j0;
       numReoV = 0; numIptV = 0;

       # apply the same orthogonalization to V
       coupling_flag = true

       if estAnorm_flag
          Anorm = fudge * sqrt(norm(B'*B, 1) )
       end

       iptU = m2 * smleps * ones(k+1)
       iptV = zeros(etype, k)

       iptUmax = zeros(etype, k)
       iptVmax = zeros(etype, k)

       j0 = j0 + 1
    end

    # examine the accuracy of the estimated inner product
    # allocate some space for the true inner product and keep the history
    if examine_flag
       iptU_his   = zeros(etype, k, k); iptV_his   = zeros(etype, k, k);
       iptU_true  = zeros(etype, k, k); iptV_true  = zeros(etype, k, k);
       iptU_after = zeros(etype, k, k); iptV_after = zeros(etype, k, k);
    end

    # ==========================================================================
    #               start the main loop of lanczos iteration
    # ==========================================================================
    j = j0
    for j = j0 : k

        # each iteration start with assign U vectors
        if beta[j] != 0.0
           U[:,j]   = p / beta[j]
        else
           # restarted with a new vector
           U[:,j]   = p
        end

        # replace the estimated Anorm with the largest Ritz value
        # and the l2 norm of A is converged
        if j == 6
           Btmp = vcat(diagm(alpha[1:j-1])+diagm(beta[2:j-1],-1), zeros(1, j-1))
           Btmp[j,j-1] = beta[j]

           # last time update the estimate of the norm of A
           Anorm = fudge * norm(Btmp, 2)
           estAnorm_flag = false
        end

        # lanczos step to generate v_j
        if j == 1

           # for the first vector
           if matrix_flag
              r = A' * U[:,1]
           else
              r = A(U[:,1], 2; params...)
           end

           alpha[1] = norm(r)
           if estAnorm_flag
              Anorm = fudge * alpha[1]
           end

        else

           # from second V vector to the last one
           if matrix_flag
              r = A' * U[:,j] - beta[j] * V[:,j-1]
           else
              r = A(U[:,j], 2; params...) - beta[j] * V[:,j-1]
           end

           alpha[j] = norm(r)

           # extended local reorthogonalize
           if alpha[j] < gamma * beta[j] && elr_flag && !forceRO_flag

              normold = alpha[j]
              condition = true
              while condition

                    t = dot(V[:,j-1], r)
                    r = r - V[:,j-1] * t
                    beta[j] = beta[j] + t

                    alpha[j] = norm(r)

                    if alpha[j] < gamma * normold
                       normold = alpha[j]
                    else
                       condition = false
                    end
                    # do I need to count this inner product?
                    # numIptV = numIptV + 1
                    # numReoV = numReoV + 1
              end
           end

           #  update the estimate of the inner product of V
           if estAnorm_flag
              if j == 2
                 tmp = fudge * sqrt(alpha[  1]^2 + beta[2]^2 + alpha[  2]*beta[  2])
              else
                 tmp = fudge * sqrt(alpha[j-1]^2 + beta[j]^2 + alpha[j-1]*beta[j-1] + alpha[j]*beta[j])
              end
              Anorm = maximum([Anorm, tmp])
           end

           # update the estimate of the inner product of V
           # only estimate when full RO is off
           if !forceRO_flag && alpha[j] != 0.0
              iptV       = updateIptV(iptV, iptU, j, alpha, beta, Anorm, eps1)
              iptVmax[j] = maximum(abs.(iptV[1:j-1]))
           end

           if examine_flag && j > 1
              iptV_his[ 1:j-1, j-1] = iptV[1:j-1]
              iptV_true[1:j-1, j-1] = V[: ,1:j-1]' * (r/alpha[j])
           end

           if elr_flag
              iptV[j-1] = n2 * smleps
           end

           # if the estimated orthogonality is worse than delta
           # select the set of vectors to Reorthogonalize against
           if (forceRO_flag || iptVmax[j] > delta || coupling_flag) && alpha[j] > 0.0

              #  force full orthogonalize or η == 0
              if forceRO_flag || eta == 0.
                 indexSet = collect(1:j-1)
              elseif !coupling_flag
                 indexSet = selectVectors(iptV, j-1, delta, eta, 0, 0, 0)
              end

              # Reorthogonalize against
              (r, alpha[j], rr) = reorthogonalization(V, r, alpha[j], indexSet, gamma, gs_iflag)

              numIptV = numIptV + rr * length(indexSet)
              numReoV = numReoV + 1

              iptV[indexSet] = n2 * eps(etype)

              coupling_flag = coupling_flag ? false : true

           end

        end

        # Check for convergence or failure to maintain semiorthogonality
        if alpha[j] < maximum([m,n]) * Anorm * smleps && j < k
           # we deflate by setting it to 0 and attempt to restart with a basis by
           # replacing r with a random starting vector
           alpha[j] = 0
           bailout  = true

           # try three times
           for attempt = 1 : 3

               r = rand(etype, m) - etype(0.5)

               if matrix_flag
                  r = A' * r
               else
                  r = A(r, 2; params...)
               end

               normold = norm(r)

               indexSet = collect(1:j-1)
               (r, normnew, rr) = reorthogonalization(V, r, normold, indexSet, gamma, gs_iflag)

               numIptV = numIptV + rr * length(indexSet)
               numReoV = numReoV + 1

               iptV[indexSet] = n2 * smleps

               if normnew > 0
                  # a new vector is find
                  bailout = false
                  break
               end

            end

            # converged
            if bailout
               j = j-1
               ierr = -j
               break
            else
               r = r / normnew
               coupling_flag = true

               # turn-off forced full reorthogonalize if delta > 0.0
               if delta > 0.0
                  forceRO_flag = false
               end
            end

        elseif j < k && !forceRO_flag && Anorm * smleps > delta*alpha[j]

            ierr = j

        end  # end of the bailout

        if alpha[j] != 0.0
           V[:,j] = r / alpha[j]
        else
           V[:,j] = r
        end

        # kepp the history
        if j > 1 && examine_flag
           iptV_after[1:j-1, j-1] = iptV[1:j-1]
           iptV_true[ 1:j-1, j-1] = V[:,1:j-1]' * r / alpha[j]
        end

        # ======================================================================
        #                                 U[:,j+1]
        # ======================================================================
        if matrix_flag
           p = A * V[:,j] - alpha[j] * U[:,j]
        else
           p = A(V[:,j], 1; params...) - alpha[j] * U[:,j]
        end
        beta[j+1] = norm(p)

        # extended local orthogonalization
        if beta[j+1] < gamma*alpha[j] && elr_flag && !forceRO_flag
           normold = beta[j+1]
           condition = true
           while condition
                 t = dot(U[:,j], p)
                 p = p - U[:,j]*t
                 beta[j+1] = norm(p)
                 alpha[j]  = alpha[j] + t

                 if beta[j+1] < gamma * normold
                    normold = beta[j+1]
                 else
                    condition = false
                 end
           end
        end

        # update the estimate of largest singular value of A
        if estAnorm_flag

           if j == 1
              tmp = fudge * sqrt(alpha[1]^2 + beta[2  ]^2)
           else
              tmp = fudge * sqrt(alpha[j]^2 + beta[j+1]^2 + alpha[j]*beta[j])
           end
           Anorm = maximum([Anorm, tmp])
        end

        # update the estimate of the orthogonality for the columns of U
        if !forceRO_flag && beta[j+1] != 0
           iptU = updateIptU(iptU, iptV, j, alpha, beta, Anorm, eps1)
           iptUmax[j] = maximum(abs.(iptU[1:j]))
        end

        # keep all the history
        if examine_flag
           # keep the history of estimated inner product
           iptU_his[ 1:j,j] = iptU[1:j]

           # keep the history of true inner product
           tmp = p / beta[j+1]
           iptU_true[1:j,j] = U[:,1:j]' * tmp
        end

        if elr_flag
           iptU[j] = m2 * smleps
        end

        # if the maximum of iptU is worse than delta, then orthogonalize U[:,]+1]
        # against subset of U[:,i] with 1 <= i <=j
        # priority forceRO > coupling > iptUmax
        if (forceRO_flag || iptUmax[j] > delta || coupling_flag) && beta[j+1] > 0.0

           # full orthogonalize
           if forceRO_flag || eta == 0
              indexSet = collect(1:j)

           # select subset which iptU[i] >= eta
           elseif !coupling_flag
              indexSet = selectVectors(iptU, j, delta, eta, 0, 0, 0)

           # aviod propagation from V to U
           else
              indexSet = vcat(indexSet, maximum(indexSet)+1)
           end

           (p, beta[j+1], rr) = reorthogonalization(U, p, beta[j+1], indexSet, gamma, gs_iflag)

           numIptU = numIptU + rr * length(indexSet)
           numReoU = numReoU + 1

           iptU[indexSet] = m2 * smleps

           # if this orthogonalization is to aviod coupling,
           # then next time doesn't need orth
           coupling_flag = coupling_flag ? false : true

        end

        # check for convergence or failure to maintain,
        # converged before j ==k
        if beta[j+1] < maximum([m,n]) * Anorm * smleps && j < k

           #  restart it with new random vector
           beta[j+1] = zero(etype)
           bailout   = true

           for attempt = 1 : 3
               p = rand(etype, n) - etype(0.5)
               if matrix_flag
                  p = A * p
               else
                  p = A(p, 1; params...)
               end

               normold = norm(p)
               indexSet = collect(1:j)
               (p, normnew, rr) = reorthogonalization(U, p, normold, indexSet, gamma, gs_iflag)

               numIptU = numIptU + rr * length(indexSet)
               numReoU = numReoU + 1

               iptU[indexSet] = m2 * smleps

               if normnew > 0.0
                  bailout = false
                  break
               end

           end

           # indicate the step error happens
           if bailout
              ierr = -j
              break
           else

              p = p / normnew
              coupling_flag = true

              if delta > 0.0
                 forceRO_flag = false
              end
           end

        elseif j<k && !forceRO_flag && Anorm * smleps > delta * beta[j+1]

           ierr = j

        end # end of convergence check

        # keep the history of the iteration
        if examine_flag
           iptU_after[1:j, j] = iptU[1:j]
           iptU_true[ 1:j, j] = U[:,1:j]' * (p/beta[j+1])
        end

    end # end of lanczos iteration

    # finish the iterations
    # incase the iteration finished before kth one
    k = j<k ? j : k
    Bk = spdiagm((alpha[1:k], beta[2:k+1]), (0,-1), k+1, k)
    if beta[k+1] > 0
       U = hcat(U[:,1:k], p/beta[k+1])
    else
       U = hcat(U[:,1:k], p)
    end
    V  = V[:,1:k]

    if examine_flag
       optout= Dict(:numReoU=>numReoU, :numReoV=>numReoV, :numIptU=>numIptU, :numIptV=>numIptV,
                    :iptU_his=>iptU_his, :iptU_true=>iptU_true, :iptU_after=>iptU_after,
                    :iptV_his=>iptV_his, :iptV_true=>iptV_true, :iptV_after=>iptV_after,
                    :iptUmax=>iptUmax, :iptVmax=>iptVmax)
       return U, Bk, V, optout
    else
       return U, Bk, V
    end

end

"""
   Gram-Schmidt reorthogonalization of r with respect to the subset of the columns
of Q, the subset is specified by indexSet. The columns of the matrix Q must be orthonormal.
"""
function reorthogonalization(Q::Matrix{Tv}, r::Vector{Tv}, normold::Tv,
                             index::Vector{Ti}, gamma::Tv, iflag::Ti) where {Tv<:Number, Ti<:Int64}

    (m, n)   = size(Q)
    normnew  = copy(normold)

    numReo  = 0
    re = copy(r)

    while normnew < gamma * normold || numReo == 0

       if iflag == 1

          t  = Q[:, index]' * re
          re = re - Q[:, index] * t

       else

          for i in index
              t = dot(Q[:,i], re)
              re= re - Q[:,i] * t
          end
       end

       normold = normnew
       normnew = norm(re)
       numReo  = numReo + 1

       if numReo > 4
          re = zeros(r)
          normnew = 0
          return re, normnew, numReo
       end

    end

    return re, normnew, numReo

end

"""
   update the estimate of the inner product of the left lanczos vectors recursively
   μ[i,i]    = 1 for 1<= i <= j+1  and  ν[j,0] = 0

   for 1 <= i <= j
   μ[j+1,i]' = α[i]ν[j,i] + β[i]ν[j,i-1] - α[j]μ[j,i]
   μ[j+1,i]  = ( μ[j+1,i]' + sign(μ[j+1,i]')ϵ ) / β[j+1]

   μ[j+1, j+1] = 1
"""
# function updateIptU(iptU_old::Vector{Tv}, iptV::Vector{Tv}, j::Ti,
#                     alpha::Vector{Tv}, beta::Vector{Tv}, Anorm::Tv, eps1::Tv) where {Tv<:AbstractFloat, Ti<:Int64}
#
#     etype = eltype(alpha)
#     iptU_new = zeros(iptU_old)
#
#     T = Anorm * eps1
#
#     # for the first iteration
#     if j == 1
#        # j=1, μ[2,1] = α[1]ν[1,1] + β[1]ν[1,0] - α[1]μ[1,1] = α[1] + 0 - α[1] = 0
#        iptU_new[1] = T / beta[2]
#
#     # later iterations
#     else
#
#        # j != 1 and i = 1, μ[j+1,1]' = α[1]ν[j,1] + β[i]ν[j,0] - α[j]μ[j,1] = α[1]ν[j,1] + 0 - α[j]μ[j,1]
#        iptU_new[1]     = alpha[1]*iptV[1]                     - alpha[j] * iptU_old[1]
#        iptU_new[1]     = (iptU_new[1] + sign(iptU_new[1]) * T) / beta[j+1]
#
#        # i = 2 : j-1,  μ[j+1,i]' = α[i]ν[j,i] + β[i]ν[j,i-1] - α[j]μ[j,i]
#        for k = 2 : j-1
#            iptU_new[k] = alpha[k]*iptV[k] + beta[k]*iptV[k-1] - alpha[j] * iptU_old[k]
#            iptU_new[k] = (iptU_new[k] + sign(iptU_new[k]) * T) / beta[j+1]
#        end
#
#        # i = j, μ[j+1,j]' = α[i]ν[j,j] + β[j]ν[j,j-1] - α[j]μ[j,j] = α[i] + β[i]ν[j,j-1] - α[j] = β[j]ν[j,j-1]
#        iptU_new[j]     =                    beta[j]*iptV[j-1]
#        iptU_new[j]     = (iptU_new[j] + sign(iptU_new[j]) * T) / beta[j+1]
#
#     end
#
#     # j+1
#     iptU_new[j+1] = one(etype)
#
#     return iptU_new
#
# end
function updateIptU(iptU_old::Vector{Tv}, iptV::Vector{Tv}, j::Ti,
                    alpha::Vector{Tv}, beta::Vector{Tv}, Anorm::Tv, eps1::Tv) where {Tv<:AbstractFloat, Ti<:Int64}

    etype = eltype(alpha)
    iptU_new = zeros(iptU_old)

    # for the first iteration
    if j == 1
       # j=1, μ[2,1] = α[1]ν[1,1] + β[1]ν[1,0] - α[1]μ[1,1] = α[1] + 0 - α[1] = 0
       T = eps1 * (sqrt(alpha[1]^2 + beta[2]^2) + sqrt(alpha[1]^2+beta[1]^2))
       T = T + eps1 * Anorm
       iptU_new[1] = T / beta[2]

    # later iterations
    else

       # j != 1 and i = 1, μ[j+1,1]' = α[1]ν[j,1] + β[i]ν[j,0] - α[j]μ[j,1] = α[1]ν[j,1] + 0 - α[j]μ[j,1]
       iptU_new[1]     = alpha[1]*iptV[1]                     - alpha[j] * iptU_old[1]
       T = eps1 * (sqrt(alpha[j]^2 + beta[j+1]^2) + sqrt(alpha[1]^2+beta[1]^2))
       T = T + eps1 * Anorm
       iptU_new[1]     = (iptU_new[1] + sign(iptU_new[1]) * T) / beta[j+1]

       # i = 2 : j-1,  μ[j+1,i]' = α[i]ν[j,i] + β[i]ν[j,i-1] - α[j]μ[j,i]
       for k = 2 : j-1
           iptU_new[k] = alpha[k]*iptV[k] + beta[k]*iptV[k-1] - alpha[j] * iptU_old[k]
           T = eps1 * (sqrt(alpha[j]^2 + beta[j+1]^2) + sqrt(alpha[k]^2+beta[k]^2))
           T = T + eps1 * Anorm
           iptU_new[k] = (iptU_new[k] + sign(iptU_new[k]) * T) / beta[j+1]
       end

       # i = j, μ[j+1,j]' = α[i]ν[j,j] + β[j]ν[j,j-1] - α[j]μ[j,j] = α[i] + β[i]ν[j,j-1] - α[j] = β[j]ν[j,j-1]
       T = eps1 * (sqrt(alpha[j]^2 + beta[j+1]^2) + sqrt(alpha[j]^2+beta[j]^2))
       T = T + eps1 * Anorm
       iptU_new[j]     =                    beta[j]*iptV[j-1]
       iptU_new[j]     = (iptU_new[j] + sign(iptU_new[j]) * T) / beta[j+1]

    end

    # j+1
    iptU_new[j+1] = one(etype)

    return iptU_new

end

"""
   update the estimate of the inner product of the right lanczos vectors recursively
   ν[i,i]    = 1 for 1<= i <= j  and  μ[j,0] = 0

   for 1 <= i <= j-1
   ν[j,i]' = β[i+1]μ[j,i+1] + α[i]μ[j,i] - β[j]ν[j-1,i]
   ν[j,i]  = ( ν[j,i]' + sign(ν[j,i]')ϵ ) / α[j]

   ν[j,j] = 1
"""
# function updateIptV(iptV_old::Vector{Tv}, iptU::Vector{Tv}, j::Ti,
#                      alpha::Vector{Tv}, beta::Vector{Tv}, Anorm::Tv, eps1::Tv) where {Tv<:AbstractFloat, Ti<:Int64}
#
#     etype = eltype(alpha)
#     iptV_new = zeros(iptV_old)
#
#     T     = Anorm * eps1
#
#     for k = 1 : j-1
#
#         iptV_new[k] = beta[k+1] * iptU[k+1] + alpha[k] * iptU[k] - beta[j] * iptV_old[k]
#         iptV_new[k] = (iptV_new[k] + sign(iptV_new[k]) * T) / alpha[j]
#
#     end
#
#     iptV_new[j] = one(etype)
#
#     return iptV_new
#
# end
function updateIptV(iptV_old::Vector{Tv}, iptU::Vector{Tv}, j::Ti,
                     alpha::Vector{Tv}, beta::Vector{Tv}, Anorm::Tv, eps1::Tv) where {Tv<:AbstractFloat, Ti<:Int64}

    etype = eltype(alpha)

    iptV_new = zeros(iptV_old)

    for k = 1 : j-1

        T = eps1 * (sqrt(alpha[k]^2 + beta[k+1]^2) + sqrt(alpha[j]^2 + beta[j]^2))
        T = T + eps1 * Anorm
        iptV_new[k] = beta[k+1] * iptU[k+1] + alpha[k] * iptU[k] - beta[j] * iptV_old[k]
        iptV_new[k] = (iptV_new[k] + sign(iptV_new[k]) * T) / alpha[j]

    end

    iptV_new[j] = one(etype)

    return iptV_new

end

"""
   select the sets of vectors to orthogonalize against it. Three strategy are provided
"""
function selectVectors(ipt::Vector{Tv}, j::Ti, delta::Tv, eta::Tv,
                       LL::Ti, strategy::Ti, extra::Ti) where {Tv<:AbstractFloat, Ti<:Int64}

    if delta < eta
       error("delta must be larger than eta")
    end

    # orthogonalize against the neighbours of ipt[i] > delta and the neighbours > eta
    if strategy == 0

       I0 = find(x->(abs(x) >= delta), ipt[1:j])
       if length(I0) == 0
          (maxipt, I0) = findmax(ipt)
       end

       index = zeros(Int64, j)

       for i = 1 : length(I0)

           # left neighbour
           r = I0[i]
           for r = I0[i] : -1 : 1
               if abs(ipt[r]) < eta || index[r] == 1  # belong to the neighbour of other vector
                  break
               else
                  index[r] = 1
               end
           end
           if extra != 0
              lower = maximum([1, r-extra+1])
              index[lower:r] = 1
           end

           # right neighbour
           s = I0[i]+1
           for s = I0[i]+1 : j
               if abs(ipt[s]) < eta || index[s] == 1
                  break
               else
                  index[s] = 1
               end
           end
           if extra != 0
              upper = minimum([j, s+extra-1])
              index[s:upper] = 1
           end

       end

       if LL > 0
          index[1:LL] = 0
       end
       index = find(x->(x>0), index)

    elseif strategy == 1

       I0 = find(x->(abs(x) > eta), ipt[1:j])
       lower = maximum([LL+1, minimum(I0)-extra])
       upper = minimum([j   , maximum(I0)+extra])
       if lower <= upper
          index = collect(lower:upper)
       end

    else

       index = find(x->(abs(x) > eta), ipt[1:j])

    end

    return index
end


text:
#Macros to make it easier to write drivers for SCPI instruments

macro scpifloat(name, instrType, scpiStr, scale, units)
	setterFunc = symbol(string(name)*"!")
	setterProto = :($(setterFunc)(instr::$instrType, val::Real))
	setterBody = :(write(instr, $scpiStr * " $val"))
	setter = Expr(:function, esc(setterProto), esc(setterBody))
	getCmd = scpiStr*"?"
	getterProto = :($name(instr::$instrType))
	getterBody = :(query(instr, $getCmd))
	getter = Expr(:function, esc(getterProto), esc(getterBody))
	Expr(:block, setter, getter)
end

macro scpibool(name, instrType, scpiStr)
	setterFunc = symbol(string(name)*"!")
	setterProto = :($(setterFunc)(instr::$instrType, val::Bool))
	setterBody = :(write(instr, $scpiStr * " " * string(int(val)) ))
	setter = Expr(:function, esc(setterProto), esc(setterBody))
	getCmd = scpiStr*"?"
	getterProto = :($name(instr::$instrType))
	getterBody = :(bool(int(query(instr, $getCmd))))
	getter = Expr(:function, esc(getterProto), esc(getterBody))
	Expr(:block, setter, getter)
end

ext: jl


text:
const a = rand(10000,10000);
@time svd(a);


text:
module Clipper
    export PolyType, PolyTypeSubject, PolyTypeClip,
           ClipType, ClipTypeIntersection, ClipTypeUnion, ClipTypeDifference, ClipTypeXor,
           PolyFillType, PolyFillTypeEvenOdd, PolyFillTypeNonZero, PolyFillTypePositive, PolyFillTypeNegative,
           JoinType, JoinTypeSquare, JoinTypeRound, JoinTypeMiter,
           EndType, EndTypeClosedPolygon, EndTypeClosedLine, EndTypeOpenSquare, EndTypeOpenRound, EndTypeOpenButt,
           Clip, add_path!, add_paths!, execute, clear!, get_bounds,
           IntPoint, IntRect, orientation, area,
           ClipperOffset

    @enum PolyType PolyTypeSubject=0 PolyTypeClip=1

    @enum ClipType ClipTypeIntersection=0 ClipTypeUnion=1 ClipTypeDifference=2 ClipTypeXor=3

    @enum PolyFillType PolyFillTypeEvenOdd=0 PolyFillTypeNonZero=1 PolyFillTypePositive=2 PolyFillTypeNegative=3

    @enum JoinType JoinTypeSquare=0 JoinTypeRound=1 JoinTypeMiter=2

    @enum EndType EndTypeClosedPolygon=0 EndTypeClosedLine=1 EndTypeOpenSquare=2 EndTypeOpenRound=3 EndTypeOpenButt=4

    @static if is_windows()
      const library_path = joinpath(dirname(@__FILE__), "cclipper.dll")
    end

    @static if is_unix()
        const library_path = joinpath(dirname(@__FILE__), "cclipper.so")
    end

    immutable IntPoint
        X::Int64
        Y::Int64
    end

    function Base.show(io::IO, point::IntPoint)
      print(io, "[$(point.X),$(point.Y)]")
    end

    function append!(outputArray::Ptr{Void}, polyIndex::Csize_t, point::IntPoint)
        ourArray = unsafe_pointer_to_objref(outputArray)::Vector{Vector{IntPoint}}

        while (polyIndex + 1) > length(ourArray)
          push!(ourArray, Vector{IntPoint}())
        end

        push!(ourArray[polyIndex + 1], point)
    end

    #==============================================================#
  	# Static functions
  	#==============================================================#
    function orientation(path::Vector{IntPoint})
        ccall((:orientation, library_path), Cuchar, (Ptr{IntPoint}, Csize_t),
              path,
              length(path)) == 1 ? true : false
    end

    function area(path::Vector{IntPoint})
      ccall((:area, library_path), Float64, (Ptr{IntPoint}, Csize_t),
            path,
            length(path))
    end

    #==============================================================#
  	# Clipper object
  	#==============================================================#
    type Clip
        clipper_ptr::Ptr{Void}

        function Clip()
            clipper = new(ccall((:get_clipper, library_path), Ptr{Void}, ()))
            finalizer(clipper, c -> ccall((:delete_clipper, library_path), Void, (Ptr{Void},), c.clipper_ptr))

            clipper
        end
    end

    function add_path!(c::Clip, path::Vector{IntPoint}, polyType::PolyType, closed::Bool)
        ccall((:add_path, library_path), Cuchar, (Ptr{Void}, Ptr{IntPoint}, Csize_t, Cint, Cuchar),
              c.clipper_ptr,
              path,
              length(path),
              Int(polyType),
              closed) == 1 ? true : false
    end

    function add_paths!(c::Clip, paths::Vector{Vector{IntPoint}}, polyType::PolyType, closed::Bool)
        lengths = Vector{UInt64}()
        for path in paths
            push!(lengths, length(path))
        end

        ccall((:add_paths, library_path), Cuchar, (Ptr{Void}, Ptr{Ptr{IntPoint}}, Ptr{Csize_t}, Csize_t, Cint, Cuchar),
              c.clipper_ptr,
              paths,
              lengths,
              length(paths),
              Int(polyType),
              closed) == 1 ? true : false
    end

    function execute(c::Clip, clipType::ClipType, subjFillType::PolyFillType, clipFillType::PolyFillType)
        polys = Vector{Vector{IntPoint}}()

        result = ccall((:execute, library_path), Cuchar, (Ptr{Void}, Cint, Cint, Cint, Any, Ptr{Void}),
                        c.clipper_ptr,
                        Int(clipType),
                        Int(subjFillType),
                        Int(clipFillType),
                        polys,
                        cfunction(append!, Any, (Ptr{Void}, Csize_t, IntPoint)))

        return result == 1 ? true : false, polys
    end

    function clear!(c::Clip)
        ccall((:clear, library_path), Void, (Ptr{Void},), c.clipper_ptr)
    end

    type IntRect
      left::Int64
      top::Int64
      right::Int64
      bottom::Int64
    end

    function get_bounds(c::Clip)
        ccall((:get_bounds, library_path), IntRect, (Ptr{Void}, ), c.clipper_ptr)
    end

    #==============================================================#
  	# ClipperOffset object
  	#==============================================================#
    type ClipperOffset
        clipper_ptr::Ptr{Void}

        function ClipperOffset(miterLimit::Float64 = 2.0, roundPrecision::Float64 = 0.25)
            clipper = new(ccall((:get_clipper_offset, library_path), Ptr{Void}, (Cdouble, Cdouble), miterLimit, roundPrecision))
            finalizer(clipper, c -> ccall((:delete_clipper_offset, library_path), Void, (Ptr{Void},), c.clipper_ptr))

            clipper
        end
    end

    function add_path!(c::ClipperOffset, path::Vector{IntPoint}, joinType::JoinType, endType::EndType)
        ccall((:add_offset_path, library_path), Void, (Ptr{Void}, Ptr{IntPoint}, Csize_t, Cint, Cint),
              c.clipper_ptr,
              path,
              length(path),
              Int(joinType),
              Int(endType))
    end

    function add_paths!(c::ClipperOffset, paths::Vector{Vector{IntPoint}}, joinType::JoinType, endType::EndType)
        lengths = Vector{UInt64}()
        for path in paths
            push!(lengths, length(path))
        end

        ccall((:add_offset_paths, library_path), Void, (Ptr{Void}, Ptr{Ptr{IntPoint}}, Ptr{Csize_t}, Csize_t, Cint, Cint),
              c.clipper_ptr,
              paths,
              lengths,
              length(paths),
              Int(joinType),
              Int(endType))
    end

    function clear!(c::ClipperOffset)
        ccall((:clear_offset, library_path), Void, (Ptr{Void},), c.clipper_ptr)
    end

    function execute(c::ClipperOffset, delta::Float64)
        polys = Vector{Vector{IntPoint}}()

        result = ccall((:execute_offset, library_path), Void, (Ptr{Void}, Cdouble, Any, Ptr{Void}),
                        c.clipper_ptr,
                        delta,
                        polys,
                        cfunction(append!, Any, (Ptr{Void}, Csize_t, IntPoint)))

        return polys
    end
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

abstract type AbstractRay{T<:Real,N} end
export AbstractRay

"""
    Ray{T,N} <: AbstractRay{T,N}

Purely geometric ray, defined as `origin + alpha * direction`.

```julia
Ray(origin::SVector{N,T}, direction::SVector{N,T})
```

Has the following accessor methods:
```julia
direction(ray::Ray{T,N}) -> SVector{N,T}
origin(ray::Ray{T,N}) -> SVector{N,T}
```
"""
struct Ray{T,N} <: AbstractRay{T,N}
    origin::SVector{N,T}
    direction::SVector{N,T}

    function Ray(origin::SVector{N,T}, direction::SVector{N,T}) where {T<:Real,N}
        return new{T,N}(origin, normalize(direction))
    end

    function Ray(origin::AbstractArray{T,1}, direction::AbstractArray{T,1}) where {T<:Real}
        @assert length(origin) == length(direction)
        N = length(origin)
        return new{T,N}(SVector{N,T}(origin), normalize(SVector{N,T}(direction)))
    end
end
export Ray

direction(ray::Ray{T,N}) where {T<:Real,N} = ray.direction
origin(ray::Ray{T,N}) where {T<:Real,N} = ray.origin
export origin, direction

function Base.print(io::IO, a::Ray{T,N}) where {T,N}
    println(io, "$(rpad("Origin:", 20)) $(origin(a))")
    println(io, "$(rpad("Direction:", 20)) $(direction(a))")
end

"""
    point(ray::AbstractRay{T,N}, alpha::T) -> SVector{T, N}

Returns a point on the ray at origin + alpha * direction. Alpha must be >= 0.
"""
function point(ray::AbstractRay{T,N}, alpha::T) where {N,T<:Real}
    @assert alpha >= zero(T) "Alpha must be nonnegative. alpha value: $alpha"
    return origin(ray) + alpha * direction(ray)
end

"""
    closestpointonray(r::Ray{T,N}, point::SVector{N,T}) -> SVector{T,N

Returns the point on the ray closest to point.
"""
function closestpointonray(r::AbstractRay{T,N}, point::SVector{N,T}) where {T,N}
    #find t that gives the smallest distance between the ray and the point
    #ray = p0 + tr̂, distance = ||line - point||
    o = origin(r)
    r̂ = direction(r)
    t = dot(point .- o, r̂)

    if t >= 0.0
        return o .+ t .* r̂
    else
        return o
    end
end

"""
    distance(r::Ray{T,N}, point::SVector{N,T}) -> Union{T,Nothing}

Returns distance to the position on the ray closest to point. If t < 0 returns nothing.
"""
function distance(r::AbstractRay{T,N}, point::SVector{N,T}) where {T,N}
    #find t that gives the smallest distance between the ray and the point
    #ray = p0 + tr̂, distance = ||line - point||
    # this has redundant code with closestpointonray but it is faster to do this way than to call closestpointonray and then comput norm(result .- point)

    o = origin(r)
    r̂ = direction(r)
    t = dot(point .- o, r̂)

    if t >= 0.0
        temp = o .+ t .* r̂
        return norm(temp .- point)
    else
        return nothing
    end
end

"""
    α(ray::AbstractRay{T,N}, point::SVector{N,T}) -> T

Computes the alpha corresponding to the closest position on the ray to point
"""
α(ray::AbstractRay{T,N}, point::SVector{N,T}) where {T<:Real,N} = dot(direction(ray), (point .- origin(ray)))

"""
Apply a Transform to a Ray object
"""
function Base.:*(a::Transform{T}, r::Ray{T,3})::Ray{T,3} where {T}
    return Ray(a * origin(r), Geometry.rotate(a, direction(r)))
end


text:
function L1(x, y)
    return sum(abs.(x .- y))
end

function L2(x, y)
    return sum((x .- y).^2)
end

function χ²(x, y)
    S = Float64.((x .- y).^2)
    S[x .!= 0] .= S[x .!= 0]./x[x .!= 0]
    return sum(S)
end

function L1l(x, y, limit)
    sum(abs.((x .≥ limit) .- (y .≥ limit)))
end

function L2l(x, y, limit)
    sum(((x .≥ limit) .- (y .≥ limit)).^2)
end

function splitmetric(x, y; metric=L1, threshold=0.5)
    projection = sum(x, dims=2) |> vec
    sp = argthreshold(projection, threshold)
    # The regions should weigh `threshold` and `(1-threshold)`
    # Add penalization equal to their difference
    l1 = metric(x[:, 1:sp],   y[:, 1:sp])
    l2 = metric(x[:, sp:end], y[:, sp:end])
    diff = abs(threshold*l1 - (1-threshold)*l2)
    return l1 + l2 + diff
end

#=
Return the index where the cumulative sum of `x` is
`threshold` of the total sum.
=#
function argthreshold(x, threshold=0.5)
    csum = cumsum(x)
    limit = threshold*csum[end]
    for (i, x_) in enumerate(csum)
        if x_ ≥ limit
            return i
        end
    end
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
module MultiArmBandit

using Statistics, StatsBase, Distributions
using DataFrames
using RCall

export Bandit, CreateBandit, Train!
export GetHistory, TrackHistory!, ClearHistory!
export GetHistoryDF

"""Bandit

struct holding multi arm bandit state

args:

QN: Dict{NamedTuple,Dict{Symbol,Float64}} Holds reward Q and sample size N for each group

Y: Contains data for each group in a Dict

keylist: vector of keys for hashing the bandit

history: Contains history of QN values if track_hist == true in CreateBandit(). User should
use GetHistory() or GetHistoryDF() to access this in a cleaner data structure.

"""
mutable struct Bandit
    QN::Dict{NamedTuple,Dict{Symbol,Float64}}
    Y::Dict{NamedTuple,Array{Float64,1}}
    keylist::Vector{NamedTuple}
    history:: Union{Vector{Dict{NamedTuple,Dict{Symbol,Float64}}},Nothing}
end

"""CreateBandit(data,groups,response;Q₀=0,track_hist=false)

This function instantiates the Multi Arm Bandit object where each arm is started with
a Q-value drawn from Unif(Q₀,Q₀+1e-10). The slight perturbation of Q from Q₀ with
uniform noise breaks ties on step 1 when the bandit is trained in the ϵ-greedy manner
to maximize Q.

args:

data: DataFrame containing data

groups: Vector of Symbols containing Groups (arms) to group by

response: Symbol corresponding to response Y

Q₀: Initial Q value to start each group at (before), 0 is default. Starting higher than the expected
reward for any group will lead to more exploration by the bandit

track_hist: boolean to track history of the bandit

"""
function CreateBandit(data,groups,response;Q₀=0,track_hist=false)
    gdf = groupby(data,groups)
    QN = Dict{NamedTuple,Dict{Symbol,Float64}}()
    Y = Dict{NamedTuple,Array{Float64,1}}()
    for key in keys(gdf)
        QN[NamedTuple(key)] = Dict(:Q=>rand(Uniform(Q₀,Q₀+1e-10)),:n=>0.0)
        Y[NamedTuple(key)] = gdf[key][:,response]
    end
    if track_hist==false
        return(Bandit(QN,Y,collect(keys(QN)),nothing))
    elseif track_hist==true
        return(Bandit(QN,Y,collect(keys(QN)),[deepcopy(QN)]))
    end
end

#Allows bandit to be hashed like a dict

Base.getindex(bandit::Bandit,key) = Base.getindex(bandit.QN,key)

Base.keys(bandit::Bandit) = Base.keys(bandit.QN)
Base.values(bandit::Bandit) = Base.values(bandit.QN)

#
"""
    argmax(bandit::Bandit)

Returns the group which has the current maximum Q value

"""
function Base.argmax(bandit::Bandit)
    keylist = bandit.keylist
    idx = argmax(getindex.(values(bandit.QN),:Q))
    return(keylist[idx])
end

""" Train!(bandit::Bandit,steps;ϵ=0.05,η=nothing)

This function trains the bandit for a given number of steps in an ϵ-greedy manner. η is an
optional learning rate to modify the update rule:

Default: Q = Q + (R-Q)/N

With Learning Rate: Q = Q + (R-Q)*η

At each step, the bandit selects a Y (response) from the group with the maximum Q at the moment as
the reward R. ϵ of the time, the bandit will select a group at random for which to select the reward.

The default update rule is derived from computing the sample mean in an online manner.

"""
function Train!(bandit::Bandit,steps;ϵ=0.05,η=nothing)
    lrn = ifelse(isnothing(η),0,1)
    track_hist = ifelse(isnothing(bandit.history),0,1)
    for i=1:steps
        u = rand(Uniform(0,1))
        if u <= ϵ
            action = sample(bandit.keylist,1)[1]
        else
            action = argmax(bandit)
        end
        Qcurrent = bandit[action][:Q]
        Reward = sample(bandit.Y[action],1)[1]
        diff = Reward - Qcurrent
        curr = bandit[action]
        curr[:n] += 1
        if lrn == 0
            curr[:Q] += 1/curr[:n] .* diff
        else
            curr[:Q] += η .* diff
        end
        if track_hist == 1
            push!(bandit.history,deepcopy(bandit.QN))
        end
    end
    return(bandit)
end

"""GetHistory(band::Bandit)

Returns the bandit Q history per group in a Dict if history was tracked

"""
function GetHistory(band::Bandit)
    isnothing(band.history) && throw(DomainError("No History was tracked"))
    intermediate_dict = Dict()
    for x in band.history, (k, v) in x
        push!(get!(intermediate_dict, k, Vector{typeof(v)}()), v)
    end
    histresult = Dict{NamedTuple,Dict{Symbol,Vector{Float64}}}()
    for key in keys(intermediate_dict)
        temp = Dict{Symbol,Vector{Float64}}()
        for x in intermediate_dict[key], (k, v) in x
            push!(get!(temp, k, Vector{Float64}()), v)
        end
        histresult[key] = temp
    end
    return(histresult)
end

"""TrackHistory!(band::Bandit)

Turns on bandit history tracking if history was not being tracked before

"""
function TrackHistory!(band::Bandit)
    isnothing(band.history)==false && throw(DomainError("Bandit History is already being tracked"))
    band.history = [deepcopy(band.QN)]
    print("\n Bandit training history will now be tracked \n")
end

""" ClearHistory(band::Bandit;off=false)

Clears the bandit history. Can also turn off history tracking if off==true

"""
function ClearHistory!(band::Bandit;off=false)
    if off==true
        band.history = nothing
        print("\n Bandit training history cleared and turned off \n")
        return(nothing)
    else
        band.history = [deepcopy(band.QN)]
        print("\n Bandit training history cleared")
        return(nothing)
    end
end

""" GetHistoryDF(band::Bandit)

Returns history of the bandit per group in a DataFrame in long format. This is recommended
if one wants to plot the history with other packages using Grammar of Graphics, such as
Gadfly in Julia or ggplot2 in R.

"""
function GetHistoryDF(band::Bandit)
    hist = GetHistory(band)
    histdf = DataFrame()
    for k in keys(hist)
        Q = hist[k][:Q]
        n_Arm = hist[k][:n]
        n=length(Q)
        d = DataFrame(steps=1:n,Group=fill(k,n),Q=Q,n_Arm=n_Arm)
        histdf = vcat(histdf,d)
    end
    histdf.Group = categorical(string.(histdf.Group))
    return(histdf)
end

""" summary(band::Bandit)

Returns a summary of the current bandit status (Q and N) as a DataFrame

"""
function Base.summary(band::Bandit)
    summarydf = DataFrame()
    for (k,v) in band.QN
        df = DataFrame(Group=k,n=v[:n],Q=v[:Q])
        summarydf = vcat(summarydf,df)
    end
    return(summarydf)
end

end

max_stars_repo_name: rokapre/MultiArmBandit

max_issues_repo_name: rokapre/MultiArmBandit

max_forks_repo_name: rokapre/MultiArmBandit


text:
module TestFinCats
using Test

using Catlab, Catlab.Theories, Catlab.CategoricalAlgebra, Catlab.Graphs
using Catlab.Graphs.BasicGraphs: TheoryGraph, TheoryReflexiveGraph

# Discrete categories
#####################

C = FinCat(FinSet(3))
@test C isa FinCat{Int,Int}
@test collect(ob_generators(C)) == 1:3
@test isempty(hom_generators(C))
@test (dom(C, 1), codom(C, 1)) == (1, 1)
@test (id(C, 2), compose(C, 2, 2)) == (2, 2)

# Categories on graphs
######################

# Free categories
#----------------

g = parallel_arrows(Graph, 3)
C = FinCat(g)
@test graph(C) == g
@test Ob(C) == FinSet(2)
@test is_free(C)
@test hom(C, 1) == Path(g, 1)
@test ob_generators(C) == 1:2
@test hom_generators(C) == 1:3

h = Graph(4)
add_edges!(h, [1,1,2,3], [2,3,4,4])
D = FinCat(h)
f = id(D, 2)
@test (src(f), tgt(f)) == (2, 2)
@test isempty(edges(f))
f = compose(D, 1, 3)
@test edges(f) == [1,3]

# Functors between free categories.
C = FinCat(parallel_arrows(Graph, 2))
F = FinFunctor((V=[1,4], E=[[1,3], [2,4]]), C, D)
@test dom(F) == C
@test codom(F) == D
@test is_functorial(F)
@test Ob(F) == FinFunction([1,4], FinSet(4))

@test ob_map(F, 2) == 4
@test hom_map(F, 1) == Path(h, [1,3])
@test collect_ob(F) == [1,4]
@test collect_hom(F) == [Path(h, [1,3]), Path(h, [2,4])]

# Composition of functors.
g, h, k = path_graph(Graph, 2), path_graph(Graph, 3), path_graph(Graph, 5)
C, D, E = FinCat(g), FinCat(h), FinCat(k)
F = FinFunctor([1,3], [[1,2]], C, D)
G = FinFunctor([1,3,5], [[1,2],[3,4]], D, E)
@test is_functorial(G)
@test hom_map(G, Path(h, [1,2])) == Path(k, [1,2,3,4])
@test F⋅G == FinFunctor([1,5], [[1,2,3,4]], C, E)
@test id(C)⋅F == F
@test F⋅id(D) == F

# Free diagrams
#--------------

C = FinCat(parallel_arrows(Graph, 2))
f, g = FinFunction([1,3], 3), FinFunction([2,3], 3)
F = FinDomFunctor([FinSet(2), FinSet(3)], [f,g], C)
@test is_functorial(F)
@test dom(F) == C
@test codom(F) isa TypeCat{<:FinSet{Int},<:FinFunction{Int}}
@test ob_map(F, 1) == FinSet(2)
@test hom_map(F, 2) == g

# `FreeDiagrams` interop.
diagram = FreeDiagram(ParallelPair(f, g))
@test FreeDiagram(F) == diagram
@test FinDomFunctor(diagram) == F

# Diagram interface.
@test diagram_type(F) <: Tuple{FinSet{Int},FinFunction{Int}}
@test cone_objects(F) == [FinSet(2), FinSet(3)]
@test cocone_objects(F) == [FinSet(2), FinSet(3)]
@test ob(limit(F)) == FinSet(1)
@test ob(colimit(F)) == FinSet(2)

# Commutative square as natural transformation.
C = FinCat(path_graph(Graph, 2))
F = FinDomFunctor([FinSet(4), FinSet(2)], [FinFunction([1,1,2,2])], C)
α₀, α₁ = FinFunction([3,4,1,2]), FinFunction([2,1])
α = FinTransformation([α₀, α₁], F, F)
@test is_natural(α)
@test (α[1], α[2]) == (α₀, α₁)
@test components(α) == [α₀, α₁]
@test α⋅α == FinTransformation([FinFunction(1:4), FinFunction(1:2)], F, F)

# Path equations
#---------------

# Simplex category truncated to one dimension.
Δ¹_graph = Graph(2)
add_edges!(Δ¹_graph, [1,1,2], [2,2,1])
Δ¹ = FinCat(Δ¹_graph, [ [1,3] => empty(Path, Δ¹_graph, 1),
                        [2,3] => empty(Path, Δ¹_graph, 1) ])
@test graph(Δ¹) == Δ¹_graph
@test length(equations(Δ¹)) == 2
@test !is_free(Δ¹)

# Symbolic categories
#####################

@present Simplex1D(FreeCategory) begin
  (V, E)::Ob
  (δ₀, δ₁)::Hom(V, E)
  σ₀::Hom(E, V)

  δ₀ ⋅ σ₀ == id(V)
  δ₁ ⋅ σ₀ == id(V)
end

Δ¹ = FinCat(Simplex1D)
@test Δ¹ isa FinCat{FreeCategory.Ob,FreeCategory.Hom}
@test ob(Δ¹, :V) isa FreeCategory.Ob
@test hom(Δ¹, :δ₀) isa FreeCategory.Hom
@test first.(ob_generators(Δ¹)) == [:V, :E]
@test first.(hom_generators(Δ¹)) == [:δ₀, :δ₁, :σ₀]
@test length(equations(Δ¹)) == 2
@test !is_free(Δ¹)

# Graph as set-valued functor on a free category.
F = FinDomFunctor(TheoryGraph, path_graph(Graph, 3))
C = dom(F)
@test is_functorial(F)
@test ob_map(F, :V) == FinSet(3)
@test hom_map(F, :src) == FinFunction([1,2], 3)
@test F(ob(C, :E)) == FinSet(2)
@test F(hom(C, :tgt)) == FinFunction([2,3], 3)
@test F(id(ob(C, :E))) == id(FinSet(2))

# Reflexive graph as set-valued functor on a category with equations.
G = FinDomFunctor(TheoryReflexiveGraph, path_graph(ReflexiveGraph, 3))
@test is_functorial(G)

# Graph homomorphisms as natural transformations.
g = parallel_arrows(Graph, 2)
add_edges!(g, [2,2], [2,2])
G = FinDomFunctor(TheoryGraph, g)
α = FinTransformation(F, G, V=FinFunction([1,2,2]), E=FinFunction([1,3],4))
@test dom_ob(α) == C
@test codom_ob(α) isa TypeCat{<:FinSet{Int},<:FinFunction{Int}}
@test is_natural(α)
@test α[:V](3) == 2
σ = FinTransformation(G, G, V=id(FinSet(2)), E=FinFunction([2,1,4,3]))
@test σ⋅σ == FinTransformation(G, G, V=id(FinSet(2)), E=FinFunction(1:4))
@test α⋅σ == FinTransformation(F, G, V=FinFunction([1,2,2]), E=FinFunction([2,4]))

# Pullback data migration by pre-whiskering.
ιV = FinFunctor([:V], [], FinCat(1), FinCat(TheoryGraph))
αV = ιV * α
@test ob_map(dom(αV), 1) == ob_map(F, :V)
@test ob_map(codom(αV), 1) == ob_map(G, :V)
@test component(αV, 1) == component(α, :V)

# Post-whiskering and horizontal composition.
ιE = FinFunctor([:E], [], FinCat(1), FinCat(TheoryGraph))
ϕ = FinTransformation([:src], ιE, ιV)
@test is_natural(ϕ)
@test component(ϕ*F, 1) == hom_map(F, :src)
@test component(ϕ*α, 1) == hom_map(F, :src) ⋅ α[:V]

end


text:
# ====================  Find Classes  ====================
"""
    find_classes(adj::BitArray{2})

returns an array of length `size(adj,1)` with each entry with index `i` being
a unique identifier for the equivalency class of the node `i`.
Uses adjacency matrix as input. I.e. graph needs to be constructed before.

# Examples
```
julia> find_classes(convert(BitArray, [0 1 0; 1 0 0; 0 0 0]))
[1, 1, 2]
```
"""
function find_classes(adj::BitArray{2})
  vlR = UnitRange{UInt32}(1:size(adj,1))
  classes = collect(vlR)                        # every node has its own class at the start
  openL = trues(size(adj,1))                    # no vertices have been visited yet
  for i in vlR                                  # check all vertices at least once
    !openL[i] && continue                       # if node already checked, continue
    classes[i] = i                              # current node has class i
    indL = [i]                                  # add current vertex to search list
    while !isempty(indL)                        # search through all reachable vertices
      j = pop!(indL)                            # first vertex in index List is next candidate
      !openL[j] && continue
      openL[j] = false                          # set current node to visited
      neighbors = findall(adj[i,:] .& openL)    # add all adjacent AND open vertices to class
      indL = union(neighbors,indL)              # without double entries
      classes[indL] .= i                        # set all neighbors to current class
    end
    sum(openL) == 0 && break                    # no more open vertices left, return
  end
  return classes
end

"""
    find_classes(m::Mapping, vl; vl_len=length(vl), sorted=false)

returns an array of length `size(adj,1)` with each entry with index `i` being
a unique identifier for the equivalency class of the node `i`.
Uses a mapping function `m(x) = y` and an vertex list `vl` to construct graph i
n place.
`m` has to be symmetric (i.e. `m(x) = y => m(y) = x`)!
In case the vertex list is closed under the mapping (i.e. there exist no `x` in
`vl` such that `m(x)` is not in `vl`), `closed` can be set to `false` in order
to improve performace.
for that class. By specifying `continuous = true`,
`vl_len` can be specified explicitly in cases where `length(vl)` does not work
(e.g. nested generator objects).
For performace reasons this function can also be called with a plain function
instead of a Mapping.

# Examples
```
```
"""
find_classes(m::Mapping, vl::AbstractArray{T,1}; vl_len=length(vl)) where T = find_classes(m.f, vl, vl_len=vl_len)

function find_classes(m::Function, vl::AbstractArray{T,1}; vl_len=length(vl)) where T
  vlR = UnitRange{UInt32}(1:vl_len)
  classes = OrderedDict(zip(vl,vlR))
  openL = Dict(zip(vl,trues(vl_len)))           # mark all vertices as open (not visited)
  searchL = Stack{eltype(vl)}()
  @time @inbounds for i in vlR                  # visit each entry at least once
    vi = vl[i]
    if openL[vi]                              # if entry already checked, continue
      push!(searchL, vi)                    # add current vertex to search list
      while !isempty(searchL)               # search through all reachable entries
        vj = pop!(searchL)                # next vertex in index List is next candidate
        classes[vj] = i
        if openL[vj]
          openL[vj] = false             # set current node to visited
          for (j,el) in enumerate(m(vj))           # add all adjacent AND open vertices to class, skip entries outside vl
            if !(el in searchL) && haskey(openL,el) && openL[el]
                push!(searchL, el)
            end
          end
        end
      end
    end
  end
  return @inbounds ReduceMap(classes)
end

"""
    find_classes(pred::Predicate, vl; isSymmetric=false)

returns an array of length `length(vl)` with each entry with index `i` being
a unique identifier for the equivalency class of the node `i`.
See `find_classes(m::Mapping, vl)` for more information.
"""
function find_classes(pred::Predicate, vl::AbstractArray; isSymmetric=false)
    adj = isSymmetric ? build_adj_matrix(pred, vl) :
        build_adj_matrix(Predicate((x,y) ->(pred(x,y) || pred(y,x))), vl)
    @inbounds res = ReduceMap(OrderedDict(zip(vl, find_classes(adj))))
    return res
end

# ====================  Find Paths  ====================
function find_classes(m::Function, vl::AbstractArray{T,1}, operations::Array{UInt32,1}; vl_len=length(vl)) where T
  vlR = UnitRange{UInt32}(1:vl_len)
  parent = Dict(zip(vl,vl))
  parent_edge = Dict(zip(vl,zeros(UInt32,vl_len)))
  openL = Dict(zip(vl,trues(vl_len)))       # mark all vertices as open (not visited)
  searchL = Stack{eltype(vl)}()
  for i in vlR              # visit each entry at least once
    vi = vl[i]
    if openL[vi]                            # if entry already checked, continue
      push!(searchL, vi)                    # add current vertex to search list
      parent[vi] = vi
      parent_edge[vi] = 0x00000
      while !isempty(searchL)               # search through all reachable entries
        vj = pop!(searchL)                  # next vertex in index List is next candidate
        if openL[vj]
          openL[vj] = false                 # set current node to visited
          for (j,el) in enumerate(m(vj))    # add all adjacent AND open vertices to class, skip entries outside vl
              if !(el in searchL) && haskey(openL,el) && openL[el]
                push!(searchL, el)
                parent[el] = vj
                parent_edge[el] = operations[j]
            end
          end
        end
      end
    end
  end
  return parent, parent_edge
end

# ====================  Auxiliary functions ====================

"""
    build_adj_matrix(pred, vl)

Constructs adjacency matrix with vertices from `vl` and edges
from `pred(i,j)` where `i,j ∈ vl`
"""
function build_adj_matrix(pred::Predicate, vl)
  adj = BitArray(undef, length(vl),length(vl))
  for (i,ind) in enumerate(vl)
    adj[i,i] = true
    for (j,ind2) in enumerate(Iterators.drop(vl,i))
      adj[i,j+i] = pred(ind,ind2) || pred(ind2,ind)
      adj[j+i,i] = adj[i,j+i]
    end
  end
  return adj
end

#merge_classes()

# ====================  EquivalenceClasses Wrappers  ====================
"""
    EquivalenceClasses(pred::Predicate, vl::AbstractArray)

Constructs an EquivalenceClasses struct from a predicate `pred`
over `n` indices in the list `vl`.
This inputs requires `n^2` operations in order to create the
adjecency matrix. For large inputs a mapping is therefore
preferabel.

# Examples
```
julia> EquivalenceClasses(Predicate((x,y)->all(x .== -y)),
                          [(i,j) for i in -2:2 for j in 4:7])
```
"""
function EquivalenceClasses(pred::Predicate, vl::AbstractArray)
    return find_classes(pred,vl)
end

"""
    EquivalenceClasses(m::Mapping, vl::AbstractArray)

Constructs an EquivalenceClasses struct from a mapping `m`
over the indices in the list `vl`.

# Examples
```
julia> EquivalenceClasses(Mapping((x)->[-x]),
                          [(i,j) for i in -2:2 for j in 4:7])
```
"""
function EquivalenceClasses(m::Mapping, vl::AbstractArray)
  return find_classes(m,vl)
end

#function parent_to_classes(parent_map::Dict{T,T})
#end


text:
using Revise, Pkg

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

Pkg.activate(pwd())
using LeverageChanges

#LeverageChanges.recompressfiles()
#LeverageChanges.recompressfile("CRSP-D")
#LeverageChanges.recompressfile("CRSP-M")

#refreshcomp (PreprocessCOMP): refreshes the compustat data from the csv.lz4 file
#refreshcrsp (PreprocessCRSP): refreshes the crsp data from the csv.lz4 file
#refreshmerge[1] (MergeCRSPCOMP): refreshes the compustat-crsp merge.
#refreshfactors (FormFactors): refreshes the compustat-crsp merge.
#[1]True if refreshcomp, refreshcrsp, or validatemerged are true
#[2]True if refreshcomp, refreshcrsp, refreshmerge, or validatemerged are true
#hello
#usequarterlydates
function leveragechangesscript()
  #LeverageChanges.recompressfile("CCM-A")
  #LeverageChanges.recompressfile("CRSP-D-I")
  #LeverageChanges.recompressfile("CRSP-D-II")
  #LeverageChanges.recompressfile("CRSP-M")
  #LeverageChanges.recompressfile("CCM")
  #LeverageChanges.recompressfile("COMP-A")

  #WHERE DOES THIS HANG? PROFILE? MAYBE GO LONG INSTEAD OF WIDE?
  LeverageChanges.replicate(
    refreshcomp = false,
    refreshtvolmcapcrsp = false,
    refreshtretcrsp = false,
    refreshdataseries=false,
    refreshportfolios=true,
    parallel = true,
    trialrunonly = true,
    table1 = false)

end

#LeverageChanges.archivefile("lev-daily.jls.lz4")
#LeverageChanges.unarchivefile("lev-daily.jls.lz4")
@time leveragechangesscript()

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# Ground state potential energy surface.

"""
    ground_pes(sys::System{S,2}, extent::NTuple{4,Float64}; lengths::NTuple{2,Float64}=(501, 501))

Compute the ground state PES for the 2-mode system `sys`.

The extent for each mode is given in `extent`, and the number of points is
given in `lengths`.
"""
function ground_pes(sys::System{S,2}, extent::NTuple{4,Float64}; lengths::NTuple{2,Int}=(501, 501)) where {S}
    M = 2
    q1_min, q1_max, q2_min, q2_max = extent
    q1_length, q2_length = lengths
    pes = zeros(q2_length, q1_length)
    for (i, q1) in enumerate(range(q1_min; stop=q1_max, length=q1_length))
        for (j, q2) in enumerate(range(q2_min; stop=q2_max, length=q2_length))
            # Lowest eigenvalue of potential matrix.
            pes[j, i] = eigvals(potential(sys, [q1, q2]))[1]
        end
    end
    pes
end

"""
    path_mean_std(sys::DiagonalSystem, beta::Float64, P::Int, s::Int, m::Int)

Compute the standard deviation of the mean of the path for surface `s` and mode
`m` in `sp`.
"""
function path_mean_std(sys::DiagonalSystem{S,M}, beta::Float64, P::Int, s::Int, m::Int) where {S,M}
    1 <= s <= S || throw(DomainError(s, "Invalid surface index."))
    1 <= m <= M || throw(DomainError(m, "Invalid mode index."))
    sp = SamplingParameters(sys, beta, P)
    centroid_std = sqrt(1 / eigvals(sp.precs[m, s])[1])
    centroid_std / sqrt(P)
end


text:
"""
    eigsolve(A::AbstractMatrix, [x₀, howmany = 1, which = :LM, T = eltype(A)]; kwargs...)
    eigsolve(f, n::Int, [howmany = 1, which = :LM, T = Float64]; kwargs...)
    eigsolve(f, x₀, [howmany = 1, which = :LM]; kwargs...)
    # expert version:
    eigsolve(f, x₀, howmany, which, algorithm)

Compute at least `howmany` eigenvalues from the linear map encoded in the matrix `A` or by
the function `f`. Return eigenvalues, eigenvectors and a `ConvergenceInfo` structure.

### Arguments:

The linear map can be an `AbstractMatrix` (dense or sparse) or a general function or
callable object. If an `AbstractMatrix` is used, a starting vector `x₀` does not need to be
provided, it is then chosen as `rand(T, size(A,1))`. If the linear map is encoded more
generally as a a callable function or method, the best approach is to provide an explicit
starting guess `x₀`. Note that `x₀` does not need to be of type `AbstractVector`; any type
that behaves as a vector and supports the required methods (see KrylovKit docs) is accepted.
If instead of `x₀` an integer `n` is specified, it is assumed that `x₀` is a regular vector
and it is initialized to `rand(T,n)`, where the default value of `T` is `Float64`, unless
specified differently.

The next arguments are optional, but should typically be specified. `howmany` specifies how
many eigenvalues should be computed; `which` specifies which eigenvalues should be
targeted. Valid specifications of `which` are given by

  - `:LM`: eigenvalues of largest magnitude
  - `:LR`: eigenvalues with largest (most positive) real part
  - `:SR`: eigenvalues with smallest (most negative) real part
  - `:LI`: eigenvalues with largest (most positive) imaginary part, only if `T <: Complex`
  - `:SI`: eigenvalues with smallest (most negative) imaginary part, only if `T <: Complex`
  - [`EigSorter(f; rev = false)`](@ref): eigenvalues `λ` that appear first (or last if
    `rev == true`) when sorted by `f(λ)`

!!! note "Note about selecting `which` eigenvalues"

    Krylov methods work well for extremal eigenvalues, i.e. eigenvalues on the periphery of
    the spectrum of the linear map. All of they valid `Symbol`s for `which` have this
    property, but could also be specified using `EigSorter`, e.g. `:LM` is equivalent to
    `Eigsorter(abs; rev = true)`. Note that smallest magnitude sorting is obtained using
    e.g. `EigSorter(abs; rev = false)`, but since no (shift-and)-invert is used, this will
    only be successful if you somehow know that eigenvalues close to zero are also close
    to the periphery of the spectrum.

!!! warning "Degenerate eigenvalues"

    From a theoretical point of view, Krylov methods can at most find a single eigenvector
    associated with a targetted eigenvalue, even if the latter is degenerate. In the case of
    a degenerate eigenvalue, the specific eigenvector that is returned is determined by the
    starting vector `x₀`. For large problems, this turns out to be less of an issue in
    practice, as often a second linearly independent eigenvector is generated out of the
    numerical noise resulting from the orthogonalisation steps in the Lanczos or Arnoldi
    iteration. Nonetheless, it is important to take this into account and to try not to
    depend on this potentially fragile behaviour, especially for smaller problems.

The argument `T` acts as a hint in which `Number` type the computation should be performed,
but is not restrictive. If the linear map automatically produces complex values, complex
arithmetic will be used even though `T<:Real` was specified. However, if the linear map and
initial guess are real, approximate eigenvalues will be searched for using a partial Schur
factorization, which implies that complex conjugate eigenvalues come in pairs and cannot
be split. It is then illegal to choose `which` in a way that would treat `λ` and `conj(λ)`
differently, i.e. `:LI` and `:SI` are invalid, as well as any `EigSorter` that would lead
to `by(λ) != by(conj(λ))`.

### Return values:

The return value is always of the form `vals, vecs, info = eigsolve(...)` with

  - `vals`: a `Vector` containing the eigenvalues, of length at least `howmany`, but could
    be longer if more eigenvalues were converged at the same cost. Eigenvalues will be real
    if [`Lanczos`](@ref) was used and complex if [`Arnoldi`](@ref) was used (see below).
  - `vecs`: a `Vector` of corresponding eigenvectors, of the same length as `vals`. Note
    that eigenvectors are not returned as a matrix, as the linear map could act on any
    custom Julia type with vector like behavior, i.e. the elements of the list `vecs` are
    objects that are typically similar to the starting guess `x₀`, up to a possibly
    different `eltype`. In particular  for a general matrix (i.e. with `Arnoldi`) the
    eigenvectors are generally complex and are therefore always returned in a complex
    number format. When the linear map is a simple `AbstractMatrix`, `vecs` will be
    `Vector{Vector{<:Number}}`.
  - `info`: an object of type [`ConvergenceInfo`], which has the following fields

      + `info.converged::Int`: indicates how many eigenvalues and eigenvectors were actually
        converged to the specified tolerance `tol` (see below under keyword arguments)
      + `info.residual::Vector`: a list of the same length as `vals` containing the
        residuals `info.residual[i] = f(vecs[i]) - vals[i] * vecs[i]`
      + `info.normres::Vector{<:Real}`: list of the same length as `vals` containing the
        norm of the residual `info.normres[i] = norm(info.residual[i])`
      + `info.numops::Int`: number of times the linear map was applied, i.e. number of times
        `f` was called, or a vector was multiplied with `A`
      + `info.numiter::Int`: number of times the Krylov subspace was restarted (see below)

!!! warning "Check for convergence"

    No warning is printed if not all requested eigenvalues were converged, so always check
    if `info.converged >= howmany`.

### Keyword arguments:

Keyword arguments and their default values are given by:

  - `verbosity::Int = 0`: verbosity level, i.e. 0 (no messages), 1 (single message
    at the end), 2 (information after every iteration), 3 (information per Krylov step)
  - `tol::Real`: the requested accuracy (corresponding to the 2-norm of the residual for
    Schur vectors, not the eigenvectors). If you work in e.g. single precision (`Float32`),
    you should definitely change the default value.
  - `krylovdim::Integer`: the maximum dimension of the Krylov subspace that will be
    constructed. Note that the dimension of the vector space is not known or checked, e.g.
    `x₀` should not necessarily support the `Base.length` function. If you know the actual
    problem dimension is smaller than the default value, it is useful to reduce the value of
    `krylovdim`, though in principle this should be detected.
  - `maxiter::Integer`: the number of times the Krylov subspace can be rebuilt; see below
    for further details on the algorithms.
  - `orth::Orthogonalizer`: the orthogonalization method to be used, see
    [`Orthogonalizer`](@ref)
  - `issymmetric::Bool`: if the linear map is symmetric, only meaningful if `T<:Real`
  - `ishermitian::Bool`: if the linear map is hermitian
  - `eager::Bool = false`: if true, eagerly compute the eigenvalue or Schur decomposition
    after every expansion of the Krylov subspace to test for convergence, otherwise wait
    until the Krylov subspace has dimension `krylovdim`. This can result in a faster return,
    for example if the initial guess is very good, but also has some overhead, as many more
    dense Schur factorizations need to be computed.

The default values are given by `tol = KrylovDefaults.tol`,
`krylovdim = KrylovDefaults.krylovdim`, `maxiter = KrylovDefaults.maxiter`,
`orth = KrylovDefaults.orth`; see [`KrylovDefaults`](@ref) for details.

The default value for the last two parameters depends on the method. If an `AbstractMatrix`
is used, `issymmetric` and `ishermitian` are checked for that matrix, otherwise the default
values are `issymmetric = false` and `ishermitian = T <: Real && issymmetric`. When values
for the keyword arguments are provided, no checks will be performed even in the matrix case.

### Algorithm

The final (expert) method, without default values and keyword arguments, is the one that is
finally called, and can also be used directly. Here, one specifies the algorithm explicitly
as either [`Lanczos`](@ref), for real symmetric or complex hermitian problems, or
[`Arnoldi`](@ref), for general problems. Note that these names refer to the process for
building the Krylov subspace, but the actual algorithm is an implementation of the
Krylov-Schur algorithm, which can dynamically shrink and grow the Krylov subspace, i.e. the
restarts are so-called thick restarts where a part of the current Krylov subspace is kept.

!!! note "Note about convergence"

    In case of a general problem, where the `Arnoldi` method is used, convergence of an
    eigenvalue is not based on the norm of the residual `norm(f(vecs[i]) - vals[i]*vecs[i])`
    for the eigenvector but rather on the norm of the residual for the corresponding Schur
    vectors.

    See also [`schursolve`](@ref) if you want to use the partial Schur decomposition
    directly, or if you are not interested in computing the eigenvectors, and want to work
    in real arithmetic all the way true (if the linear map and starting guess are real).
"""
function eigsolve end

"""
    EigSorter(by; rev = false)

A simple `struct` to be used in combination with [`eigsolve`](@ref) or [`schursolve`](@ref)
to indicate which eigenvalues need to be targeted, namely those that appear first when
sorted by `by` and possibly in reverse order if `rev == true`.
"""
struct EigSorter{F}
    by::F
    rev::Bool
end
EigSorter(f::F; rev = false) where {F} = EigSorter{F}(f, rev)

const Selector = Union{Symbol,EigSorter}

eigsolve(
    A::AbstractMatrix,
    howmany::Int = 1,
    which::Selector = :LM,
    T::Type = eltype(A);
    kwargs...
) = eigsolve(A, rand(T, size(A, 1)), howmany, which; kwargs...)

eigsolve(f, n::Int, howmany::Int = 1, which::Selector = :LM, T::Type = Float64; kwargs...) =
    eigsolve(f, rand(T, n), howmany, which; kwargs...)
function eigsolve(f, x₀, howmany::Int = 1, which::Selector = :LM; kwargs...)
    Tx = typeof(x₀)
    Tfx = Core.Compiler.return_type(apply, Tuple{typeof(f),Tx})
    T = Core.Compiler.return_type(dot, Tuple{Tx,Tfx})
    alg = eigselector(f, T; kwargs...)
    checkwhich(which) || error("Unknown eigenvalue selector: which = $which")
    if alg isa Lanczos
        if which == :LI || which == :SI
            error(
                "Eigenvalue selector which = $which invalid: real eigenvalues expected with Lanczos algorithm"
            )
        end
    elseif T <: Real
        if which == :LI ||
           which == :SI ||
           (which isa EigSorter && which.by(+im) != which.by(-im))
            error(
                "Eigenvalue selector which = $which invalid because it does not treat
          `λ` and `conj(λ)` equally: work in complex arithmetic by providing a complex starting vector `x₀`"
            )
        end
    end
    return eigsolve(f, x₀, howmany, which, alg)
end

function eigselector(
    f,
    T::Type;
    issymmetric::Bool = false,
    ishermitian::Bool = issymmetric && !(T <: Complex),
    krylovdim::Int = KrylovDefaults.krylovdim,
    maxiter::Int = KrylovDefaults.maxiter,
    tol::Real = KrylovDefaults.tol,
    orth::Orthogonalizer = KrylovDefaults.orth,
    eager::Bool = false,
    verbosity::Int = 0
)
    if (issymmetric && !(T <: Complex)) || ishermitian
        return Lanczos(
            krylovdim = krylovdim,
            maxiter = maxiter,
            tol = tol,
            orth = orth,
            eager = eager,
            verbosity = verbosity
        )
    else
        return Arnoldi(
            krylovdim = krylovdim,
            maxiter = maxiter,
            tol = tol,
            orth = orth,
            eager = eager,
            verbosity = verbosity
        )
    end
end
function eigselector(
    A::AbstractMatrix,
    T::Type;
    issymmetric::Bool = T <: Real && LinearAlgebra.issymmetric(A),
    ishermitian::Bool = issymmetric || LinearAlgebra.ishermitian(A),
    krylovdim::Int = KrylovDefaults.krylovdim,
    maxiter::Int = KrylovDefaults.maxiter,
    tol::Real = KrylovDefaults.tol,
    orth::Orthogonalizer = KrylovDefaults.orth,
    eager::Bool = false,
    verbosity::Int = 0
)
    if (T <: Real && issymmetric) || ishermitian
        return Lanczos(
            krylovdim = krylovdim,
            maxiter = maxiter,
            tol = tol,
            orth = orth,
            eager = eager,
            verbosity = verbosity
        )
    else
        return Arnoldi(
            krylovdim = krylovdim,
            maxiter = maxiter,
            tol = tol,
            orth = orth,
            eager = eager,
            verbosity = verbosity
        )
    end
end

checkwhich(::EigSorter) = true
checkwhich(s::Symbol) = s in (:LM, :LR, :SR, :LI, :SI)

eigsort(s::EigSorter) = s.by, s.rev
function eigsort(which::Symbol)
    if which == :LM
        by = abs
        rev = true
    elseif which == :LR
        by = real
        rev = true
    elseif which == :SR
        by = real
        rev = false
    elseif which == :LI
        by = imag
        rev = true
    elseif which == :SI
        by = imag
        rev = false
    else
        error("invalid specification of which eigenvalues to target: which = $which")
    end
    return by, rev
end


text:
module Solution
  export areaOrPerimeter
  function areaOrPerimeter(l::T, w::T)::T where T<:Number
      l == w ? l*w : 2*(l+w)
  end
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
struct HomogeneousArray{T,N,A<:Union{<:AbstractVector,<:NTuple},P} <: AbstractArray{T,N}
    data::A
    size::NTuple{N,Int}
    HomogeneousArray{T,N,A,P}(data::A,size::NTuple{N,Int}) where {T,N,A,P} = (@assert 1 <= P <= N; new{T,N,A,P}(data,size))
end

@inline HomogeneousArray{P}(data,size) where P = HomogeneousArray{eltype(data),length(size),typeof(data),P}(data,size)

Base.IndexStyle(::Type{<:HomogeneousArray}) = IndexCartesian()

@inline Base.size(a::HomogeneousArray) = a.size

@inline Base.@propagate_inbounds function Base.getindex(v::HomogeneousArray{T,N,A,P},I::Vararg{Int,N}) where {T,N,A,P}
    d = v.data
    i = @inbounds I[P]
    return d[i]
end


text:
export entropy, entropy!, entanglement_entropy

"""
    entropy(pvector)
    entropy(ρ)

Compute the shannon entropy: S(p) = -∑ᵢ pᵢ log pᵢ
Where pᵢ are the eigenvalues of ρ.
"""
function entropy(evals::AbstractVector)
	T = typeof(log(one(eltype(evals))))
	cutoff = convert(T, 1e-16)
	- sum(eval*log2(eval) for eval in Iterators.filter(>(cutoff), evals); init=zero(T))
end

entropy(ρ::AbstractMatrix) = entropy(eigvals(ρ))

"""
    entropy!(ρ)

Destructive version of [`entropy`](@ref).
"""
entropy!(ρ::AbstractMatrix) = entropy(eigvals!(ρ))

"""
    entanglement_entropy(ψ[, N₁])

Compute the entanglement entropy of the left N₁ spins of the chain. Default value for N₁ is
half of the chain.
"""
function entanglement_entropy(state, N1=div(nspins(state), 2))
    return entropy(abs2.(svdvals(reshape(state, 2^N1, :))))
end


text:
#=
Задача 10
ДАНО: Робот - в юго-западном углу поля, на котором расставлено некоторое количество маркеров
РЕЗУЛЬТАТ: Функция верула значение средней температуры всех замаркированных клеток
=#
function temperature!(r::Robot)::Real
    side = Ost
    count = 0
    total = 0
    res = 0
    total = search(r,side,res,count)
    return total
end

function search(r::Robot,side::HorizonSide, count::Int, res::Int)::Real
    while (isborder(r,side)==false)
        if (ismarker(r) == true)
            res += temperature(r)
            count += 1
        end
        move!(r,side)
    end
    if (isborder(r,Nord)==false)
        move!(r,Nord)
        if (ismarker(r) == true)
            res += temperature(r)
            count += 1
        end
        side = invers(side)
        search(r,side,count,res)
    end
    return (res/ count)
end


text:
# @testset "Test optimization with iterative solver" begin
#     # Test 1
#     Ne = 20
#     C = SparseRadialMapComponent(2, [0,0])
#     X = randn(2, Ne).*randn(2, Ne)
#
#     xricardo = optimize(C, X, 0.0, 0.0)
#     xiter = iterative(C, X, 0.0, 0.0)
#
#     @test norm(xricardo-xiter)<1e-3
#
#     # Test  2
#     Ne = 200
#     C = SparseRadialMapComponent(3, [2,2,0])
#     X = randn(3, Ne).*randn(3, Ne)
#
#     λ = 0.1
#
#
#     xricardo = optimize(C, X, λ, 0.0)
#     xiter = iterative(C, X, 0.1, 0.0)
#
#     @test norm(xricardo-xiter)<1e-3
#
#     # Test 3
#     Ne = 200
#     C = SparseRadialMapComponent(6, [2, -1,2,-1, 2,0])
#     X = randn(6, Ne).*randn(6, Ne)
#
#     λ = 0.1
#
#
#     xricardo = optimize(C, X, λ, 0.0)
#     xiter = iterative(C, X, 0.1, 0.0)
#     if norm(xricardo-xiter)>1e-3
#         xiter = iterative(C, X, 0.1, 0.0)
#     end
#     @test norm(xricardo-xiter)<1e-3
#
# end


text:
using RDatasets, DataFrames
iris_dataframe = dataset("datasets", "iris")

iris_dataframe[:id] = 1:size(iris_dataframe, 1)
iris_stack = stack(iris, [1:4])

iris_stack_petal = stack(iris_dataframe, [:PetalLength, :PetalWidth], :Species)

iris_melt = melt(iris_dataframe, :Species)

iris_default_stack = stack(iris_dataframe)

max_stars_repo_path: Chapter 2/7.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 2/7.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 2/7.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
using DelimitedFiles

elv = readdlm("./temp/elevation.txt")
ssp = readdlm("./temp/ssp.txt")

# Normalise ground level
elv[:,2] = elv[:,2] .- 1000.0
ssp[:,1] = ssp[:,1] .- 1000.0

#### Geometry
# Point source 100 Hz at 80 m above ground level
source = Source(
    frequency = 100,
    height = -80 + (-169) # 169 is ground level
)

# Multiple recievers on grid
receiver = Receiver(
    depth_point = 1001,
    range_point = 10001,
    depth = Vec2(-1000f0,0f0),
    range = Vec2(0f0,10f0)
)

# Terrain geometry
terrain = Terrain(
    interp_type = "C",
    profile = (Vec2(elv[:,1],elv[:,2]))
)

### Boundary condition
Zc = 12.81 + 11.62im
Theta, Rmag, Rphase = R_coeff(Zc;len=100)

brc = Vec3(Theta,Rmag,Rphase)
trc = Vec3([0f0,45f0,90f0],[0f0,0f0,0f0],[0f0,0f0,0f0])

reflection = Reflection_Coeff(
    top_coeff = trc,
    bottom_coeff = brc
)

### Sound speed profile

sspl = Vec2(ssp[:,1],ssp[:,2]) # case 1
ssp = SSP(sound_speed_profile= sspl)

### Analysis
opt = Analysis(
    filename = "RealCase_Bellhop_f100",
    analyse = "CG",
    option1 = "CFW",
    option2 = "F*",
    num_ray = 1601,
    alpha = Vec2(-80.0f0,80.0f0),
    box = Vec2(10f0,1000.0f0),
    step= 0
)

Environment(opt,source, receiver,ssp,terrain,reflection)

# Run Bellhop
fn= opt.filename
filename = "temp\\$fn"
run_bellhop = `bellhop $filename`
@time run(run_bellhop)

p1 = PlotRay("$filename.ray",
        xlabs = "Range, m",
        ylabs = "Height, m")

        plot!(p1,elv[:,1]*1000,elv[:,2],
        lw=0,fill = 0, color = "#8c510a",legend = false)
        yflip!(true)

        scatter!(p1,[0], [source.height],
        markersize = 6, color = "#4daf4a")

#savefig(p1,"ray_real.png")

p2 = PlotShd("$filename.shd";
        xlabs = "Range, m",
        ylabs = "Transmission loss, dB",
        cblabs = "dB",
        climb = (40,120))

        plot!(p2,elv[:,1]*1000,elv[:,2],
        lw=0,fill = 0, color = "#8c510a", legend = false)
        yflip!(true)

        scatter!(p2,[0], [source.height],
        markersize = 6, color = "#4daf4a",
        label = "Source")

#savefig(p2,"trans_real.png")
rmfile(filename)


text:
module AsyncFinalizersTests

include("utils.jl")
include("test_bags.jl")
include("test_finalizers.jl")
include("test_chaos.jl")

end  # module AsyncFinalizersTests


text:
include("./euclideanGeometry.jl")
include("./sparseEuclideanFunctions.jl")


text:
module VideoIO

using Libdl
using FixedPointNumbers, ColorTypes, ImageCore, Requires, Dates

libpath = joinpath(@__DIR__, "..", "deps", "usr", "bin")

include("init.jl")
include("util.jl")
include(joinpath(av_load_path, "AVUtil", "src", "AVUtil.jl"))
include(joinpath(av_load_path, "AVCodecs", "src", "AVCodecs.jl"))
include(joinpath(av_load_path, "AVFormat", "src", "AVFormat.jl"))
include(joinpath(av_load_path, "AVDevice", "src", "AVDevice.jl"))
include(joinpath(av_load_path, "SWScale", "src", "SWScale.jl"))

using .AVUtil
using .AVCodecs
using .AVFormat
using .SWScale

if have_avdevice()
    import .AVDevice
end

include("info.jl")
include("avdictionary.jl")
include("avio.jl")
include("encoding.jl")
include("testvideos.jl")
using .TestVideos

if Sys.islinux()
    import Glob
    function init_camera_devices()
        append!(CAMERA_DEVICES, Glob.glob("video*", "/dev"))
        DEFAULT_CAMERA_FORMAT[] = AVFormat.av_find_input_format("video4linux2")
    end
    function init_camera_settings()
        DEFAULT_CAMERA_OPTIONS["framerate"] = 30
        DEFAULT_CAMERA_DEVICE[] = isempty(CAMERA_DEVICES) ? "" : CAMERA_DEVICES[1]
    end
end

if Sys.iswindows()
    function init_camera_devices()
        append!(CAMERA_DEVICES, get_camera_devices(ffmpeg, "dshow", "dummy"))
        DEFAULT_CAMERA_FORMAT[] = AVFormat.av_find_input_format("dshow")
    end
    function init_camera_settings()
        DEFAULT_CAMERA_OPTIONS["framerate"] = 30
        DEFAULT_CAMERA_DEVICE[] = string(
            "video=",
            isempty(CAMERA_DEVICES) ? "0" : CAMERA_DEVICES[1]
        )
    end
end

if Sys.isapple()
    function init_camera_devices()
        try
            append!(CAMERA_DEVICES, get_camera_devices(ffmpeg, "avfoundation", "\"\""))
            DEFAULT_CAMERA_FORMAT[] = AVFormat.av_find_input_format("avfoundation")
        catch
            try
                append!(CAMERA_DEVICES, get_camera_devices(ffmpeg, "qtkit", "\"\""))
                DEFAULT_CAMERA_FORMAT[] = AVFormat.av_find_input_format("qtkit")
            catch
            end
        end
    end
    function init_camera_settings()
        DEFAULT_CAMERA_OPTIONS["framerate"] = 30
        # Note: "Integrated" is another possible default value
        DEFAULT_CAMERA_OPTIONS["pixel_format"] = "uyvy422"
        DEFAULT_CAMERA_DEVICE[] = isempty(CAMERA_DEVICES) ? "0" : CAMERA_DEVICES[1]
    end
end

#Helper functions to explain about Makie load order requirement
function play(f; flipx=false, flipy=false)
    error("Makie must be loaded before VideoIO to provide video playback functionality. Try a new session with `using Makie, VideoIO`")
end
function playvideo(video;flipx=false,flipy=false)
    error("Makie must be loaded before VideoIO to provide video playback functionality. Try a new session with `using Makie, VideoIO`")
end
function viewcam(device=DEFAULT_CAMERA_DEVICE, format=DEFAULT_CAMERA_FORMAT)
    if have_avdevice()
        error("Makie must be loaded before VideoIO to provide camera playback functionality. Try a new session with `using Makie, VideoIO`")
    else
        error("No AV device detected")
    end
end

function __init__()
    # Always check your dependencies from `deps.jl`
    # TODO remove uncessary ENV["LD_LIBRARY_PATH"] from check_deps, so that
    # it doesn't mess with LD_LIBRARY_PATH, which was causing CI download issues due to issues with julia's curl
    # since check_deps is optional, I hope this is ok for now

    #check_deps()

    read_packet[] = @cfunction(_read_packet, Cint, (Ptr{AVInput}, Ptr{UInt8}, Cint))

    av_register_all()

    if have_avdevice()
        AVDevice.avdevice_register_all()
        init_camera_devices()
        init_camera_settings()
    end

    @require Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" begin
        # Define read and retrieve for Images
        function play(f; flipx=false, flipy=false)
            scene = Makie.Scene(resolution = (f.width, f.height))
            buf = read(f)
            makieimg = Makie.image!(scene,buf, show_axis = false, scale_plot = false)[end]
            Makie.rotate!(scene, -0.5pi)
            if flipx && flipy
                Makie.scale!(scene, -1, -1, 1)
            else
                flipx && Makie.scale!(scene, -1, 1, 1)
                flipy && Makie.scale!(scene, 1, -1, 1)
            end
            display(scene)
            while !eof(f) && isopen(scene)
                read!(f, buf)
                makieimg[1] = buf
                sleep(1 / f.framerate)
            end
        end

        function playvideo(video;flipx=false,flipy=false)
            f = VideoIO.openvideo(video)
            play(f,flipx=flipx,flipy=flipy)
        end

        if have_avdevice()
            function viewcam(device=DEFAULT_CAMERA_DEVICE, format=DEFAULT_CAMERA_FORMAT)
                init_camera_settings()
                camera = opencamera(device[], format[])
                play(camera, flipx=true)
                close(camera)
            end
        end
    end
end

end # VideoIO


text:
function build_wall!(vis::Visualizer, walls::AbstractVector{<:AbstractWall}; name::Symbol=:Env, α=0.2)
    for i in eachindex(walls)
        build_wall!(vis, walls[i], name=name, id=i, α=α)
    end
    return nothing
end

function build_cylinder!(vis::Visualizer, walls::AbstractVector{<:AbstractWall}; radius=0.0, name::Symbol=:Env, α=0.2)
    for i in eachindex(walls)
        build_cylinder!(vis, walls[i], radius=radius, name=name, id=i, α=α)
    end
    return nothing
end

function build_wall!(vis::Visualizer, wall::Wall3D; name::Symbol=:Env, id::Int, α=0.2)
    p1 = wall.p1
    p2 = wall.p2
    p3 = wall.p3
    p4 = p3 + (p1 - p2)
    l1 = norm(p1 - p2)
    l2 = norm(p3 - p2)
    v1 = (p1 - p2)/l1
    v2 = (p3 - p2)/l2
    orange_mat, blue_mat, black_mat = get_material(;α=α)

    setobject!(vis[name]["wall$id"][:slab],
        MeshCat.HyperRectangle(MeshCat.Vec(0.0, 0.0, 0.0), MeshCat.Vec(l1, l2, 0.0)),
        black_mat)
    for j = 1:3
        setobject!(vis[name]["wall$id"]["p$j"],
            GeometryBasics.Sphere(MeshCat.Point3f0(0.0), 0.05),
            black_mat)
    end

    rt = MeshCat.LinearMap([v1 v2 cross(v1, v2)])
    settransform!(vis[name]["wall$id"][:slab], MeshCat.compose(MeshCat.Translation(p2...), rt))
    settransform!(vis[name]["wall$id"]["p1"], MeshCat.Translation(p1...))
    settransform!(vis[name]["wall$id"]["p2"], MeshCat.Translation(p2...))
    settransform!(vis[name]["wall$id"]["p3"], MeshCat.Translation(p3...))
    return nothing
end

function build_cylinder!(vis::Visualizer, wall::CylinderWall; radius=0.0, name::Symbol=:Env, id::Int, α=0.2)
    p = wall.p
    v = wall.v
    v_ = Int.([v == :x, v == :y, v == :z])
    l = convert(Float32, wall.l)
    r = convert(Float32, wall.r-radius)
    orange_mat, blue_mat, black_mat = get_material(;α=α)

    setobject!(vis[name]["cylinder$id"][:body],
        MeshCat.Cylinder(MeshCat.Point3f0(0.0), MeshCat.Point3f0((l*v_)...), r),
        black_mat)

    # MeshCat.HyperRectangle(MeshCat.Vec(0.0, 0.0, 0.0),
    #     MeshCat.Vec(l1, l2, 0.0)), MeshPhongMaterial(color = MeshCat.RGBA(0.5, 0.5, 0.5, α)))
    setobject!(vis[name]["cylinder$id"]["p"], GeometryBasics.Sphere(MeshCat.Point3f0(0.0), 0.05), black_mat)

    # rt = MeshCat.LinearMap([v1 v2 cross(v1, v2)])
    # settransform!(vis[name]["cylinder$id"][:body], MeshCat.compose(MeshCat.Translation(p...), rt))
    settransform!(vis[name]["cylinder$id"][:body], MeshCat.Translation(p...))
    settransform!(vis[name]["cylinder$id"]["p"], MeshCat.Translation(p...))
    return nothing
end