text:
####
# Associated
####

"""
    AssociatedWeighted(P)

We normalise so that `orthogonalityweight(::Associated)` is a probability measure.
"""
struct AssociatedWeight{T,OPs<:AbstractQuasiMatrix{T}} <: Weight{T}
    P::OPs
end
axes(w::AssociatedWeight) = (axes(w.P,1),)

sum(::AssociatedWeight{T}) where T = one(T)

"""
    Associated(P)

constructs the associated orthogonal polynomials for P, which have the Jacobi matrix

    jacobimatrix(P)[2:end,2:end]

and constant first term. Or alternatively

    w = orthogonalityweight(P)
    A = recurrencecoefficients(P)[1]
    Associated(P) == (w/(sum(w)*A[1]))'*((P[:,2:end]' - P[:,2:end]) ./ (x' - x))

where `x = axes(P,1)`.
"""

struct Associated{T, OPs<:AbstractQuasiMatrix{T}} <: OrthogonalPolynomial{T}
    P::OPs
end

associated(P) = Associated(P)

axes(Q::Associated) = axes(Q.P)
==(A::Associated, B::Associated) = A.P == B.P

orthogonalityweight(Q::Associated) = AssociatedWeight(Q.P)

function associated_jacobimatrix(X::Tridiagonal)
    c,a,b = subdiagonaldata(X),diagonaldata(X),supdiagonaldata(X)
    Tridiagonal(c[2:end], a[2:end], b[2:end])
end

function associated_jacobimatrix(X::SymTridiagonal)
    a,b = diagonaldata(X),supdiagonaldata(X)
    SymTridiagonal(a[2:end], b[2:end])
end
jacobimatrix(a::Associated) = associated_jacobimatrix(jacobimatrix(a.P))

associated(::ChebyshevT{T}) where T = ChebyshevU{T}()
associated(::ChebyshevU{T}) where T = ChebyshevU{T}()

const StieltjesPoint{T,W<:Number,V,D} = BroadcastQuasiMatrix{T,typeof(inv),Tuple{BroadcastQuasiMatrix{T,typeof(-),Tuple{W,QuasiAdjoint{V,Inclusion{V,D}}}}}}
const LogKernelPoint{T<:Real,C,W<:Number,V,D} = BroadcastQuasiMatrix{T,typeof(log),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{BroadcastQuasiMatrix{C,typeof(-),Tuple{W,QuasiAdjoint{V,Inclusion{V,D}}}}}}}}
const ConvKernel{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(-),Tuple{D1,QuasiAdjoint{T,D2}}}
const Hilbert{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(inv),Tuple{ConvKernel{T,Inclusion{T,D1},Inclusion{T,D2}}}}
const LogKernel{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(log),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{ConvKernel{T,Inclusion{T,D1},Inclusion{T,D2}}}}}}
const PowKernel{T,D1,D2,F<:Real} = BroadcastQuasiMatrix{T,typeof(^),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{ConvKernel{T,Inclusion{T,D1},Inclusion{T,D2}}}},F}}

# recognize structure of W = ((t .- x).^a
const PowKernelPoint{T,V,D,F} =  BroadcastQuasiVector{T, typeof(^), Tuple{ContinuumArrays.AffineQuasiVector{T, V, Inclusion{V, D}, T}, F}}

@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, w::ChebyshevTWeight)
    T = promote_type(eltype(H), eltype(w))
    zeros(T, axes(w,1))
end

@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, w::ChebyshevUWeight)
    T = promote_type(eltype(H), eltype(w))
    convert(T,π) * axes(w,1)
end

@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, w::LegendreWeight)
    T = promote_type(eltype(H), eltype(w))
    x = axes(w,1)
    log.(x .+ one(T)) .- log.(one(T) .- x)
end

@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wT::Weighted{<:Any,<:ChebyshevT})
    T = promote_type(eltype(H), eltype(wT))
    ChebyshevU{T}() * _BandedMatrix(Fill(-convert(T,π),1,∞), ℵ₀, -1, 1)
end

@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wU::Weighted{<:Any,<:ChebyshevU})
    T = promote_type(eltype(H), eltype(wU))
    ChebyshevT{T}() * _BandedMatrix(Fill(convert(T,π),1,∞), ℵ₀, 1, -1)
end

@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wP::Weighted{<:Any,<:OrthogonalPolynomial})
    P = wP.P
    w = orthogonalityweight(P)
    A = recurrencecoefficients(P)[1]
    Q = associated(P)
    (-A[1]*sum(w))*[zero(axes(P,1)) Q] + (H*w) .* P
end

@simplify *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, P::Legendre) = H * Weighted(P)

##
# OffHilbert
##

@simplify function *(H::Hilbert{<:Any,<:Any,<:ChebyshevInterval}, W::Weighted{<:Any,<:ChebyshevU})
    tol = eps()
    x = axes(H,1)
    T̃ = chebyshevt(x)
    ψ_1 = T̃ \ inv.(x .+ sqrtx2.(x)) # same ψ_1 = x .- sqrt(x^2 - 1) but with relative accuracy as x -> ∞
    M = Clenshaw(T̃ * ψ_1, T̃)
    data = zeros(eltype(ψ_1), ∞, ∞)
    # Operator has columns π * ψ_1^k
    copyto!(view(data,:,1), convert(eltype(data),π)*ψ_1)
    for j = 2:∞
        mul!(view(data,:,j),M,view(data,:,j-1))
        norm(view(data,:,j)) ≤ tol && break
    end
    # we wrap in a Padded to avoid increasing cache size
    T̃ * PaddedArray(chop(paddeddata(data), tol), size(data)...)
end

# Hilbert transform of weighted ChebyshevU is well-defined for x∈ℝ via extended ChebyshevT.
# Domain of integration can be whole real line....

@simplify function *(H::Hilbert{<:Any,<:DomainSets.RealNumbers,<:DomainSets.RealNumbers}, W::Weighted{<:Any,<:ChebyshevU})
    T = promote_type(eltype(H), eltype(W))
    ExtendedChebyshevT{T}() * _BandedMatrix(Fill(convert(T,π),1,∞), ℵ₀, 1, -1)
end

# ... or just [-1,1] and we have the same answer.
@simplify function *(H::Hilbert{<:Any,<:DomainSets.RealNumbers,<:ChebyshevInterval}, W::Weighted{<:Any,<:ChebyshevU})
    T = promote_type(eltype(H), eltype(W))
    ExtendedChebyshevT{T}() * _BandedMatrix(Fill(convert(T,π),1,∞), ℵ₀, 1, -1)
end

###
# LogKernel
###

@simplify function *(L::LogKernel{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wT::Weighted{<:Any,<:ChebyshevT})
    T = promote_type(eltype(L), eltype(wT))
    ChebyshevT{T}() * Diagonal(Vcat(-convert(T,π)*log(2*one(T)),-convert(T,π)./(1:∞)))
end

@simplify function *(H::LogKernel, wT::Weighted{<:Any,<:SubQuasiArray{<:Any,2,<:ChebyshevT,<:Tuple{AbstractAffineQuasiVector,Slice}}})
    V = promote_type(eltype(H), eltype(wT))
    kr,jr = parentindices(wT.P)
    @assert axes(H,1) == axes(H,2) == axes(wT,1)
    T = parent(wT.P)
    x = axes(T,1)
    W = Weighted(T)
    A = kr.A
    T[kr,:] * Diagonal(Vcat(-convert(V,π)*(log(2*one(V))+log(abs(A)))/A,-convert(V,π)./(A * (1:∞))))
end

###
# PowKernel
###

@simplify function *(K::PowKernel, wT::Weighted{<:Any,<:Jacobi})
    T = promote_type(eltype(K), eltype(wT))
    cnv,α = K.args
    x,y = K.args[1].args[1].args
    @assert x' == y
    β = (-α-1)/2
    error("Not implemented")
    # ChebyshevT{T}() * Diagonal(1:∞)
end

####
# StieltjesPoint
####

stieltjesmoment_jacobi_normalization(n::Int,α::Real,β::Real) = 2^(α+β)*gamma(n+α+1)*gamma(n+β+1)/gamma(2n+α+β+2)

@simplify function *(S::StieltjesPoint, w::AbstractJacobiWeight)
    α,β = w.a,w.b
    z,_ = parent(S).args[1].args
    (x = 2/(1-z);stieltjesmoment_jacobi_normalization(0,α,β)*HypergeometricFunctions.mxa_₂F₁(1,α+1,α+β+2,x))
end

@simplify function *(S::StieltjesPoint, wP::Weighted)
    P = wP.P
    w = orthogonalityweight(P)
    X = jacobimatrix(P)
    z, xc = parent(S).args[1].args
    z in axes(P,1) && return transpose(view(inv.(xc' .- xc) * wP,z,:))
    transpose((X'-z*I) \ [-sum(w)*_p0(P); zeros(∞)])
end

sqrtx2(z::Number) = sqrt(z-1)*sqrt(z+1)
sqrtx2(x::Real) = sign(x)*sqrt(x^2-1)

@simplify function *(S::StieltjesPoint, wP::Weighted{<:Any,<:ChebyshevU})
    T = promote_type(eltype(S), eltype(wP))
    z, xc = parent(S).args[1].args
    z in axes(wP,1) && return  (convert(T,π)*ChebyshevT()[z,2:end])'
    ξ = inv(z + sqrtx2(z))
    transpose(convert(T,π) * ξ.^oneto(∞))
end

####
# LogKernelPoint
####

@simplify function *(L::LogKernelPoint, wP::Weighted{<:Any,<:ChebyshevU})
    T = promote_type(eltype(L), eltype(wP))
    z, xc = parent(L).args[1].args[1].args
    if z in axes(wP,1)
        Tn = Vcat(convert(T,π)*log(2*one(T)), convert(T,π)*ChebyshevT{T}()[z,2:end]./oneto(∞))
        return transpose((Tn[3:end]-Tn[1:end])/2)
    else
        # for U_k where k>=1
        ξ = inv(z + sqrtx2(z))
        ζ = (convert(T,π)*ξ.^oneto(∞))./oneto(∞)
        ζ = (ζ[3:end]- ζ[1:end])/2

        # for U_0
        ζ = Vcat(convert(T,π)*(ξ^2/4 - (log.(abs.(ξ)) + log(2*one(T)))/2), ζ)
        return transpose(ζ)
    end

end

##
# mapped
###

@simplify function *(H::Hilbert, w::SubQuasiArray{<:Any,1})
    T = promote_type(eltype(H), eltype(w))
    m = parentindices(w)[1]
    # TODO: mapping other geometries
    @assert axes(H,1) == axes(H,2) == axes(w,1)
    P = parent(w)
    x = axes(P,1)
    (inv.(x .- x') * P)[m]
end

@simplify function *(H::Hilbert, wP::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
    T = promote_type(eltype(H), eltype(wP))
    kr,jr = parentindices(wP)
    W = parent(wP)
    x = axes(H,1)
    t = axes(H,2)
    t̃ = axes(W,1)
    if x == t
        (inv.(t̃ .- t̃') * W)[kr,jr]
    else
        M = affine(t,t̃)
        @assert x isa Inclusion
        a,b = first(x),last(x)
        x̃ = Inclusion((M.A * a .+ M.b)..(M.A * b .+ M.b)) # map interval to new interval
        Q̃,M = arguments(*, inv.(x̃ .- t̃') * W)
        parent(Q̃)[affine(x,axes(parent(Q̃),1)),:] * M
    end
end

@simplify function *(S::StieltjesPoint, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
    P = parent(wT)
    z, x = parent(S).args[1].args
    z̃ = inbounds_getindex(parentindices(wT)[1], z)
    x̃ = axes(P,1)
    (inv.(z̃ .- x̃') * P)[:,parentindices(wT)[2]]
end

@simplify function *(L::LogKernelPoint, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
    P = parent(wT)
    z, xc = parent(L).args[1].args[1].args
    kr, jr = parentindices(wT)
    z̃ = inbounds_getindex(kr, z)
    x̃ = axes(P,1)
    c = inv(kr.A)
    LP = log.(abs.(z̃ .- x̃')) * P
    Σ = sum(P; dims=1)
    transpose((c*transpose(LP) + c*log(c)*vec(Σ))[jr])
end

@simplify function *(L::LogKernel, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}})
    V = promote_type(eltype(L), eltype(wT))
    wP = parent(wT)
    kr, jr = parentindices(wT)
    x = axes(wP,1)
    T = wP.P
    @assert T isa ChebyshevT
    D = T \ (log.(abs.(x .- x')) * wP)
    c = inv(2*kr.A)
    T[kr,:] * Diagonal(Vcat(2*convert(V,π)*c*log(c), 2c*D.diag.args[2]))
end

### generic fallback
for Op in (:Hilbert, :StieltjesPoint, :LogKernelPoint, :PowKernelPoint, :LogKernel, :PowKernel)
    @eval begin
        @simplify function *(H::$Op, wP::WeightedBasis{<:Any,<:Weight,<:Any})
            w,P = wP.args
            Q = OrthogonalPolynomial(w)
            (H * Weighted(Q)) * (Q \ P)
        end
        @simplify *(H::$Op, wP::Weighted{<:Any,<:SubQuasiArray{<:Any,2}}) = H * view(Weighted(parent(wP.P)), parentindices(wP.P)...)
    end
end

###
# Interlace
###

@simplify function *(H::Hilbert, W::PiecewiseInterlace)
    axes(H,2) == axes(W,1) || throw(DimensionMismatch())
    Hs = broadcast(function(a,b)
                x,t = axes(a,1),axes(b,1)
                H = inv.(x .- t') * b
                H
            end, [W.args...], permutedims([W.args...]))
    N = length(W.args)
    Ts = [broadcastbasis(+, broadcast(H -> H.args[1], Hs[k,:])...) for k=1:N]
    Ms = broadcast((T,H) -> unitblocks(T\H), Ts, Hs)
    PiecewiseInterlace(Ts...) * BlockBroadcastArray{promote_type(eltype(H),eltype(W))}(hvcat, N, permutedims(Ms)...)
end

@simplify function *(H::StieltjesPoint, S::PiecewiseInterlace)
    z, xc = parent(H).args[1].args
    axes(H,2) == axes(S,1) || throw(DimensionMismatch())
    @assert length(S.args) == 2
    a,b = S.args
    xa,xb = axes(a,1),axes(b,1)
    Sa = inv.(z .- xa') * a
    Sb = inv.(z .- xb') * b
    transpose(BlockBroadcastArray(vcat, unitblocks(transpose(Sa)), unitblocks(transpose(Sb))))
end

@simplify function *(L::LogKernelPoint, S::PiecewiseInterlace)
    z, xc = parent(L).args[1].args[1].args
    axes(L,2) == axes(S,1) || throw(DimensionMismatch())
    @assert length(S.args) == 2
    a,b = S.args
    xa,xb = axes(a,1),axes(b,1)
    Sa = log.(abs.(z .- xa')) * a
    Sb = log.(abs.(z .- xb')) * b
    transpose(BlockBroadcastArray(vcat, unitblocks(transpose(Sa)), unitblocks(transpose(Sb))))
end

#################################################
# ∫f(x)g(x)(t-x)^a dx evaluation where f and g given in coefficients
#################################################

####
# cached operator implementation
####
# Constructors support BigFloat and it's recommended to use them for high orders.
mutable struct PowerLawMatrix{T, PP<:Normalized{<:Any,<:Legendre{<:Any}}} <: AbstractCachedMatrix{T}
    P::PP # OPs - only normalized Legendre supported for now
    a::T  # naming scheme follows (t-x)^a
    t::T
    data::Matrix{T}
    datasize::Tuple{Int,Int}
    function PowerLawMatrix{T, PP}(P::PP, a::T, t::T) where {T, PP<:AbstractQuasiMatrix}
        new{T, PP}(P,a,t, gennormalizedpower(a,t,50),(50,50))
    end
end
PowerLawMatrix(P::AbstractQuasiMatrix, a::T, t::T) where T = PowerLawMatrix{T,typeof(P)}(P,a,t)
size(K::PowerLawMatrix) = (ℵ₀,ℵ₀) # potential to add maximum size of operator
copy(K::PowerLawMatrix{T,PP}) where {T,PP} = K # Immutable entries

# data filling
#TODO: fix the weird inds
cache_filldata!(K::PowerLawMatrix, inds, _) = fillcoeffmatrix!(K, inds)

# because it really only makes sense to compute this symmetric operator in square blocks, we have to slightly rework some of LazyArrays caching and resizing
function getindex(K::PowerLawMatrix, k::Int, j::Int)
    resizedata!(K, k, j)
    K.data[k, j]
end
function getindex(K::PowerLawMatrix, kr::AbstractUnitRange, jr::AbstractUnitRange)
    resizedata!(K, maximum(kr),maximum(jr))
    K.data[kr, jr]
end
function resizedata!(K::PowerLawMatrix, n::Integer, m::Integer)
    olddata = K.data
    νμ = size(olddata)
    nm = max.(νμ,max(n,m))
    if νμ ≠ nm
        K.data = similar(K.data, nm...)
        K.data[axes(olddata)...] = olddata
    end
    if maximum(nm) > maximum(νμ)
        inds = maximum(νμ):maximum(nm)
        cache_filldata!(K, inds, inds)
        K.datasize = nm
    end
    K
end

####
# methods
####
function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), K::PowKernelPoint, Q::Normalized{<:Any,<:Legendre})
    tx,a = K.args
    t = tx[zero(typeof(a))]
    return Q*PowerLawMatrix(Q,a,t)
end

####
# operator generation
####

# evaluates ∫(t-x)^a Pn(x)P_m(x) dx for the case m=0, i.e. first row of operator, via backwards recursion
function powerleg_backwardsfirstrow(a::T, t::T, ℓ::Int) where T <: Real
    ℓ = ℓ+200
    coeff = zeros(BigFloat,ℓ)
    coeff[end-1] =  one(BigFloat)
    for m = reverse(1:ℓ-2)
        coeff[m] = (coeff[m+2]-t/((a+m+2)/(2*m+1))*normconst_Pnadd1(m,t)*coeff[m+1])/(((a-m+1)/(2*m+1))/((a+m+2)/(2*m+1))*normconst_Pnsub1(m,t))
    end
    coeff = PLnorminitial00(t,a)/coeff[1].*coeff
    return T.(coeff[1:ℓ-200])
end
# modify recurrence coefficients to work for normalized Legendre
normconst_Pnadd1(m::Int, settype::T) where T<:Real = sqrt(2*m+3*one(T))/sqrt(2*m+one(T))
normconst_Pnsub1(m::Int, settype::T) where T<:Real = sqrt(2*m+3*one(T))/sqrt(2*m-one(T))
normconst_Pmnmix(n::Int, m::Int, settype::T) where T<:Real = sqrt(2*m+3*one(T))*sqrt(2*n+one(T))/(sqrt(2*m+one(T))*sqrt(2*n-one(T)))
# useful explicit initial case
function PLnorminitial00(t::Real, a::Real)
    return ((t+1)^(a+1)-(t-1)^(a+1))/(2*(a+1))
end

# compute r-th coefficient of product expansion of order p and order q normalized Legendre polynomials
productseriescfs(p::T, q::T, r::T) where T = sqrt((2*p+1)*(2*q+1)/((2*(p+q-2*r)+1)))*(2*(p+q-2*r)+1)/(2*(p+q-r)+1)*exp(loggamma(r+one(T)/2)+loggamma(p-r+one(T)/2)+loggamma(q-r+one(T)/2)-loggamma(q-r+one(T))-loggamma(p-r+one(T))-loggamma(r+one(T))-loggamma(q+p-r+one(T)/2)+loggamma(q+p-r+one(T)))/π

# # this generates the entire operator via normalized product Legendre decomposition
# # this is very stable but scales rather poorly with high orders, so we only use it for testing
# function productoperator(a::T, t::T, ℓ::Int) where T
#     op::Matrix{T} = zeros(T,ℓ,ℓ)
#     # first row where n arbitrary and m==0
#     first = powerleg_backwardsfirstrow(a,t,2*ℓ+1)
#     op[1,:] = first[1:ℓ]
#     # generate remaining rows
#     for p = 1:ℓ-1
#         for q = p:ℓ-1
#             productcfs = zeros(T,2*ℓ+1)
#             for i = 0:min(p,q)
#                 productcfs[1+q+p-2*i] = productseriescfs(p,q,i)
#             end
#             op[p+1,q+1] = dot(first,productcfs)
#         end
#     end
#     # matrix is symmetric
#     for m = 1:ℓ
#         for n = m+1:ℓ
#             op[n,m] = op[m,n]
#         end
#     end
#     return op
# end

# This function returns the full ℓ×ℓ dot product operator, relying on several different methods for first row, second row, diagonal and remaining elements. We don't use this outside of the initial block.
function gennormalizedpower(a::T, t::T, ℓ::Int) where T <: Real
    # initialization
    ℓ = ℓ+3
    coeff = zeros(T,ℓ,ℓ)
    # construct first row via stable backwards recurrence
    first = powerleg_backwardsfirstrow(a,t,2*ℓ+1)
    coeff[1,:] = first[1:ℓ]
    # contruct second row via normalized product Legendre decomposition
    @inbounds for q = 1:ℓ-1
        productcfs = zeros(T,2*ℓ+1)
        productcfs[q+2] = productseriescfs(1,q,0)
        productcfs[q] = productseriescfs(1,q,1)
        coeff[2,q+1] = dot(first,productcfs)
    end
    # contruct the diagonal via normalized product Legendre decomposition
    @inbounds for q = 2:ℓ-1
        productcfs = zeros(T,2*ℓ+1)
        @inbounds for i = 0:q
            productcfs[1+2*q-2*i] = productseriescfs(q,q,i)
        end
        coeff[q+1,q+1] = dot(first,productcfs)
    end
    #the remaining cases can be constructed iteratively by means of a T-shaped recurrence
    @inbounds for m = 2:ℓ-2
        # build remaining row elements
        @inbounds for j = 1:m-2
            coeff[j+2,m+1] = (t/((j+1)*(a+m+j+2)/((2*j+1)*(m+j+1)))*normconst_Pnadd1(j,t)*coeff[j+1,m+1]+((a+1)*m/(m*(m+1)-j*(j+1)))/((j+1)*(a+m+j+2)/((2*j+1)*(m+j+1)))*normconst_Pmnmix(m,j,t)*coeff[j+1,m]-(j*(a+m-j+1)/((2*j+1)*(m-j)))*1/((j+1)*(a+m+j+2)/((2*j+1)*(m+j+1)))*normconst_Pnsub1(j,t)*coeff[j,m+1])
        end
    end
    #matrix is symmetric
    @inbounds for m = 1:ℓ
        @inbounds for n = m+1:ℓ
            coeff[n,m] = coeff[m,n]
        end
    end
    return coeff[1:ℓ-3,1:ℓ-3]
end

# the following version takes a previously computed block that has been resized and fills in the missing data guided by indices in inds
function fillcoeffmatrix!(K::PowerLawMatrix, inds::AbstractUnitRange)
    # the remaining cases can be constructed iteratively
    a = K.a; t = K.t; T = eltype(promote(a,t));
    ℓ = maximum(inds)
    # fill in first row via stable backwards recurrence
    first = powerleg_backwardsfirstrow(a,t,2*ℓ+1)
    K.data[1,inds] = first[inds]
    # fill in second row via normalized product Legendre decomposition
    @inbounds for q = minimum(inds):ℓ-1
        productcfs = zeros(T,2*ℓ+1)
        productcfs[q+2] = productseriescfs(1,q,0)
        productcfs[q] = productseriescfs(1,q,1)
        K.data[2,q+1] = dot(first,productcfs)
    end
    # fill in the diagonal via normalized product Legendre decomposition
    @inbounds for q = minimum(inds):ℓ-1
        productcfs = zeros(T,2*ℓ+1)
        @inbounds for i = 0:q
            productcfs[1+2*q-2*i] = productseriescfs(q,q,i)
        end
        K.data[q+1,q+1] = dot(first,productcfs)
    end
    @inbounds for m in inds
        m = m-2
        # build remaining row elements
        @inbounds for j = 1:m-1
            K.data[j+2,m+2] = (t/((j+1)*(a+m+j+3)/((2*j+1)*(m+j+2)))*normconst_Pnadd1(j,t)*K.data[j+1,m+2]+((a+1)*(m+1)/((m+1)*(m+2)-j*(j+1)))/((j+1)*(a+m+j+3)/((2*j+1)*(m+j+2)))*normconst_Pmnmix(m+1,j,t)*K.data[j+1,m+1]-(j*(a+m-j+2)/((2*j+1)*(m+1-j)))*1/((j+1)*(a+m+j+3)/((2*j+1)*(m+j+2)))*normconst_Pnsub1(j,t)*K.data[j,m+2])
        end
    end
    # matrix is symmetric
    @inbounds for m in reverse(inds)
        K.data[m,1:end] = K.data[1:end,m]
    end
end

function dot(v::AbstractVector{T}, W::PowerLawMatrix, q::AbstractVector{T}) where T
    vpad, qpad = paddeddata(v), paddeddata(q)
    vl, ql = length(vpad), length(qpad)
    return dot(vpad,W[1:vl,1:ql]*qpad)
end


text:
"""
    surf(data::Array{<:Real,2}, target::Array{<:Integer,1}, m::Signed=-1,
         dist_func::Any=(e1, e2) -> sum(abs.(e1 .- e2), dims=2);
         f_type::String="continuous")::Array{Float64,1}

Compute feature weights using SURF algorithm. The f_type argument specifies whether the features are continuous or discrete
and can either have the value of "continuous" or "discrete".

---
# Reference:
- Casey S. Greene, Nadia M. Penrod, Jeff Kiralis, and Jason H. Moore.
Spatially uniform ReliefF (SURF) for computationally-efficient filtering
of gene-gene interactions. BioData mining, 2(1):5–5, Sep 2009.
"""
function surf(data::Array{<:Real,2}, target::Array{<:Integer,1}, m::Signed=-1,
              dist_func::Any=(e1, e2) -> sum(abs.(e1 .- e2), dims=2);
              f_type::String="continuous")::Array{Float64,1}

    # Initialize feature weights vector.
    weights = zeros(Float64, 1, size(data, 2))

    # Compute vectors of maximum and minimum feature values.
    max_f_vals = maximum(data, dims=1)
    min_f_vals = minimum(data, dims=1)

    # Sample m samples without replacement. If m has signal value -1, use all samples.
    sample_idxs = StatsBase.sample(1:size(data, 1), if (m==-1) size(data,1) else m end, replace=false)
    if (m == -1) m = size(data, 1) end # If m has signal value -1, set m to total number of examples.

    # Compute pairwise distances between samples (vector form).
    dists = Array{Float64}(undef, Int64((size(data, 1)^2 - size(data, 1))/2 + 1))

    dists[1] = 0  # Set first value of distances vector to 0 - accessed when i == j in square form indices.

    # Construct pairwise distances vector using vectorized distance function.
    top_ptr = 2
    @inbounds for idx = 1:size(data,1)-1
        upper_lim = top_ptr + size(data, 1) - idx - 1
        dists[top_ptr:upper_lim] = dist_func(data[idx:idx, :], data[idx+1:end, :])
        top_ptr = upper_lim + 1
    end

    # Get mean distance between all samples.
    mean_dist = Statistics.mean(dists[2:end])

    # Go over sampled indices.
    @inbounds for idx = sample_idxs

        # Row and column indices for querying pairwise distance vector.
        row_idxs = repeat([idx - 1], size(data, 1))
        col_idxs = collect(0:size(data, 1)-1)

        # Get indices in distance vector (from square form indices).
        neigh_idx = square_to_vec(row_idxs, col_idxs, size(data, 1)) .+ 2

        # Query distances to neighbors and get mask for neighbors that fall within
        # hypersphere with radius mean_dist.
        neigh_mask = dists[neigh_idx[neigh_idx .!= 1]] .<= mean_dist
        insert!(neigh_mask, idx, 0)

        # Get masks for nearest hits and nearest misses.
        hit_neigh_mask = neigh_mask .& (target .== target[idx])
        miss_neigh_mask = neigh_mask .& (target .!= target[idx])

        # Ignore samples that have no hits or misses in radius.
        if sum(hit_neigh_mask) == 0 || sum(miss_neigh_mask) == 0
            continue
        end

        # Compute probability weights for misses in considered region.
        miss_classes = target[miss_neigh_mask]

        # Compute weights for nearest misses and compute weighting vector.
        weights_mult = Array{Float64}(undef, length(miss_classes))  # Allocate weights multiplier vector.
        cm = countmap(miss_classes)                                 # Count unique values.
        u = collect(keys(cm))
        c = collect(values(cm))
        neighbor_weights = c ./ length(miss_classes)       # Compute misses' weights.
        @inbounds for (w, val) = zip(neighbor_weights, u)  # Build multiplier vector.
            find_res = findall(miss_classes .== val)
            weights_mult[find_res] .= w
        end

        ### Weights Update ###

        if f_type == "continuous"
            # If features continuous.

            # Penalty term
            penalty = sum(abs.(data[idx:idx, :] .- data[hit_neigh_mask, :])./(max_f_vals .- min_f_vals .+ eps(Float64)), dims=1)

            # Reward term
            reward = sum(weights_mult .* abs.(data[idx:idx, :] .- data[miss_neigh_mask, :])./(max_f_vals .- min_f_vals .+ eps(Float64)), dims=1)

            # Weights update
            weights = weights .- penalty ./ (size(data, 1)*sum(hit_neigh_mask) + eps(Float64)) .+
            reward ./ (size(data, 1)*sum(miss_neigh_mask) + eps(Float64))

        elseif f_type == "discrete"
            # If features discrete.

            # Penalty term
            penalty = sum(Int64.(data[idx:idx, :] .!= data[hit_neigh_mask, :]), dims=1)

            # Reward term
            reward = sum(weights_mult .* Int64.(data[idx:idx, :] .!= data[miss_neigh_mask, :]), dims=1)

            # Weights update
            weights = weights .- penalty ./ (size(data, 1)*sum(hit_neigh_mask) + eps(Float64)) .+
            reward ./ (size(data, 1)*sum(miss_neigh_mask) + eps(Float64))

        else
            throw(DomainError(f_type, "f_type can only be equal to \"continuous\" or \"discrete\"."))
        end

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

    end

    # Return computed feature weights.
    return vec(weights)
end

ext: jl


text:
# Homology: type and methods

abstract type AbstractHomology end
grouptype(::Type{AbstractHomology}) = Nothing
group(h::AbstractHomology, dim::Int; kw...) = throw(MethodError(group,(typeof(h),Int)))

"""
Homology group iterator for an abstract complex
"""
struct Homology{C<:AbstractComplex, PID} <: AbstractHomology
    complex::C
end
homology(::Type{PID}, c::C) where {C<:AbstractComplex, PID} = Homology{C, PID}(c)

show(io::IO, h::Homology{C, PID}) where {C<:AbstractComplex, PID} = print(io, "Homology($PID)[$(h.complex)]")

#
# Interface methods
#

grouptype(::Type{Homology{C, PID}}) where {C<:AbstractComplex, PID} = PID

function group(h::Homology{C, PID}, p::Int; Dₚ::Int=0) where {C<:AbstractComplex, PID}
    cdim = dim(h.complex)
    @assert cdim >= p "Cannot define $p-th homology group for $cdim-dimensional complex"

    M = boundary(h.complex, p+1, PID)
    F = smith(M)
    D = diagm(F)

    nₚ, nₚ₊₁ = size(D)
    Dₚ₊₁ = trivial = 0

    # empty boundary
    nₚ₊₁ == 0 && return (0, 0, (F.S, F.Sinv, F.T, F.Tinv, D, Dₚ₊₁))

    # Calculate rank B_p = rank D_{p+1}
    for i in 1:nₚ₊₁
        D[i,i] == zero(PID) && break
        Dₚ₊₁ += 1
    end

    # rank of torsion-free subgroup
    for i in 1:min(nₚ, nₚ₊₁)
        if D[i,i] == one(PID) || D[i,i] == -one(PID)
           trivial += 1
        end
    end

    # β_p = rank H p = rank Z_p - rank B_p = n_p − rank D_p - rank D_{p+1}
    βₚ = nₚ - Dₚ - Dₚ₊₁
    τₚ = Dₚ₊₁ - trivial

    return (βₚ, τₚ, (F.S, F.Sinv, F.T, F.Tinv, D, Dₚ₊₁))
end

#
# Iterator methods
#
"""Return homology group type: dimension, Betti & torsion numbers."""
eltype(h::Homology) = Tuple{Int, Int, Int}

length(h::Homology) = dim(h.complex)

function iterate(h::Homology{C, PID}, state=nothing) where {C<:AbstractComplex, PID}
    if state === nothing
        Z = zeros(PID,0,0)
        snfstate = (Z,Z,Z,Z,Z,0)
        return iterate(h, (0, snfstate))
    end

    p = state[1]
    p >= dim(h.complex) && return nothing

    βₚ, τₚ, snfstate = group(h, p, Dₚ = state[2][end])
    return ((p, βₚ, τₚ), (p+1, snfstate))
end

"""
Homology generator iterator
"""
struct WithGenerators{H <: AbstractHomology}
    homology::H
end

"""
    withgenerators(hom)

An iterator that yields homology generators from a homology iterator `hom`.

Returns homology group parameters from iterator `hom` and generators as pair of `Chain` **x** and coefficient **k** that compose boundary **kx**.
When **k** is zero, then **x** is a boundary without any coefficient.
"""
withgenerators(h::H) where {H <: AbstractHomology} = WithGenerators{H}(h)

#
# Iterator methods
#
length(g::WithGenerators) = length(g.homology)
eltype(g::WithGenerators) = Tuple{Int,Int,Int,Dict{Chain, grouptype(typeof(g.homology))}}

function iterate(g::WithGenerators, state=nothing)
    h = g.homology

    # calculate betti, torsion & SNF
    result = iterate(h, state)
    result === nothing && return nothing

    GT = grouptype(typeof(h))
    chains = Dict{AbstractChain, GT}()

    p, βₚ, τₚ = result[1]
    U, Uinv, V, Vinv, D, t = result[2][2]
    cplx = g.homology.complex
    trivial = t - τₚ
    cls = cells(h.complex, p)

    A = view(U, 1:size(U,1), t+1:size(U,2))
    B = view(Uinv, t+1:size(Uinv,1), 1:size(Uinv,2))
    C = if p == 0
        n = size(cplx, p)
        SparseMatrixCSC{GT}(I, n, n)
    else
        VinvPrev = state[2][4]
        tPrev = state[2][6]
        view(VinvPrev, 1:size(VinvPrev,1), tPrev+1:size(VinvPrev,2))
    end

    P = A * (B * C)
    F = smith(P)
    G = F.S * diagm(F)

    # betti generator
    for j in 1:size(G,2)
        ii, V = findnz(G[:,j])
        length(ii) == 0 && continue
        ch = Chain(p, GT)
        for (i,v) in zip(ii,V)
            push!(ch, hash(cls[i])=>v)
        end
        chains[ch] = 0
    end

    # torsion generator
    for j in 1:τₚ
        ch = Chain(p, GT)
        I, V = findnz(D[:,trivial+j])
        for (i,v) in zip(I,V)
            push!(ch, hash(cls[i])=>v)
        end
        chains[ch] = D[trivial+j, trivial+j]
    end

    return (p, βₚ, τₚ, chains), result[2]
end

#
# Auxiliary  methods
#

"""Calculate Betti number"""
betti(g::AbstractHomology) = tuple((gst[2] for gst in g)...)

"""Calculate Euler characteristic"""
euler(btn::Tuple) = sum(isodd(i) ? v : -v for (i,v) in enumerate(btn))

"""Return linearly independent generators"""
function generators(g::WithGenerators{H}) where {H <: AbstractHomology}
    chains = Dict{Int,Vector{Chain}}()
    for (d,b,t,gens) in g
        chains[d] = collect(keys(gens))
    end
    return chains
end


text:
module SuffixArrays

export suffixsort, lcp

const CodeUnits = Union{UInt8,UInt16}
const IndexTypes = Union{Int8,Int16,Int32,Int64}
const IndexVector = AbstractVector{<:IndexTypes}

include("sais.jl")

function suffixsort(V::AbstractVector{U}, base::Integer=1) where {U<:CodeUnits}
    0 ≤ base || throw(ArgumentError("unsupported negative indexing base: $base"))
    n = length(V)
    # unsigned index type to return
    T = n+base-1 ≤ typemax(UInt8)  ? UInt8  :
        n+base-1 ≤ typemax(UInt16) ? UInt16 :
        n+base-1 ≤ typemax(UInt32) ? UInt32 : UInt64
    n ≤ 1 && return fill(T(base), n)
    # signed index type for algorithm
    S = n ≤ typemax(Int8)  ? Int8  :
        n ≤ typemax(Int16) ? Int16 :
        n ≤ typemax(Int32) ? Int32 : Int64
    if sizeof(T) == sizeof(S)
        I = zeros(T, n)
        sais(V, reinterpret(S, I), 0, n, Int(typemax(U))+1, false)
        base ≠ 0 && (I .+= base)
        return I
    else
        I = zeros(S, n)
        sais(V, I, 0, n, Int(typemax(U))+1, false)
        I′ = Vector{T}(undef, n)
        @inbounds for (i, x) in enumerate(I)
            I′[i] = (x + base) % T
        end
        return I′
    end
end

function suffixsort(s::AbstractString, base::Integer=1)
    return suffixsort(codeunits(s), base)
end

"""
    lcp(sa, s[, base])

Compute the longest common prefix (LCP) array from the suffix array `sa`
associated with sequence `s`.

reference:
Linear-Time Longest-Common-Prefix Computation in Suffix Arrays and Its Applications
Kasai et. al.
  http://web.cs.iastate.edu/~cs548/references/linear_lcp.pdf
"""
function lcp(sa, V::AbstractVector{U}, base::Integer=1) where {U<:CodeUnits}
    T = eltype(sa)
    pos = sa .+ T(1-base)
    n = length(pos)
    lcparr = similar(pos)
    rank = invperm(pos)
    h = 0
    for i in 1:n
        if rank[i] == 1
            continue
        end
        j = pos[rank[i]-1]
        maxh = n - max(i, j)
        while h <= maxh && V[i+h] == V[j+h]
            h += 1
        end
        lcparr[rank[i]] = h
        h = max(h-1, 0)
    end
    lcparr[1] = 0
    lcparr
end

end # module


text:
macro define_struct(classname, fields...)
    quote
        mutable struct $classname
            $(fields...)
            function $classname($(fields...))
                new($(fields...))
            end
        end
    end
end
#@macroexpand @define_struct foo bar baz
#x = foo(1,2)
#x.bar
#x.baz


text:
@testset "contains" begin
    html = "foo"
    block = PlutoStaticHTML.code_block(html)
    @test contains(block, "<b>foo</b>")

    notebook = Notebook([
        Cell("x = 1 + 1"),
        Cell("using Images: load"),
        Cell("PKGDIR = \"$PKGDIR\""),
        Cell("""im_file(ext) = joinpath(PKGDIR, "test", "im", "im.\$ext")"""),
        Cell("""load(im_file("png"))""")
    ])
    html, nb = notebook2html_helper(notebook)
    lines = split(html, '\n')

    @test contains(lines[1], "1 + 1")
    @test contains(lines[2], "2")
    @test contains(lines[2], "class='documenter-example-output'")

    @test contains(lines[end-1], """<img src=\"data:image/png;base64,""")

    notebook = Notebook([
        Cell("struct A end"),
        Cell("""
            struct B
                x::Int
                a::A
            end
            """
            ),
        Cell("B(1, A())")
    ])
    use_distributed = false
    html, nb = notebook2html_helper(notebook; use_distributed)
    lines = split(html, '\n')
    @test contains(lines[end-1], "B(1, A())")

    notebook = Notebook([
        Cell("md\"my text\"")
    ])
    html, nb = notebook2html_helper(notebook, HTMLOptions(; hide_md_code=true); use_distributed)
    lines = split(html, '\n')
    @test lines[1] == ""

    html, nb = notebook2html_helper(notebook, HTMLOptions(; hide_md_code=false); use_distributed)
    lines = split(html, '\n')
    @test lines[1] != ""

    opts = HTMLOptions(; hide_md_code=false, hide_code=true)
    html, nb = notebook2html_helper(notebook, opts; use_distributed)
    lines = split(html, '\n')
    @test lines[1] == ""
end

@testset "use_distributed=false and pwd" begin
    dir = pwd()
    nb = Notebook([Cell("1 + 1")])
    _, _ = notebook2html_helper(nb, HTMLOptions(; hide_md_code=true); use_distributed=false)
    @test pwd() == dir
end

@testset "hide" begin
    nb = Notebook([
        Cell("x = 1 + 1 # hide"),
        Cell("""
            # hideall
            y = 3 + 3
            """)
        ])
    html, nb = notebook2html_helper(nb; use_distributed=false)
    @test contains(html, ">2<")
    @test !contains(html, "1 + 1")
    @test contains(html, ">6<")
    @test !contains(html, "3 + 3")
end

@testset "from_file" begin
    mktempdir() do dir
        file = joinpath(dir, "tmp.jl")
        content = pluto_notebook_content("x = 1 + 2")
        write(file, content)
        html = PlutoStaticHTML.notebook2html(file)
        @test contains(html, "3")
    end
end

@testset "run_notebook!_errors" begin
    # To avoid changing the current working directory; I don't know what causes it to change
    # exactly.
    cd(pwd()) do
        mktempdir() do dir
            text = pluto_notebook_content("sum(1, :b)")
            path = joinpath(dir, "notebook.jl")
            write(path, text)
            session = ServerSession()
            session.options.evaluation.workspace_use_distributed = false
            err = nothing
            try
                nb = PlutoStaticHTML.run_notebook!(path, session)
            catch err
            end

            @test err isa Exception
            msg = sprint(showerror, err)
            @test contains(msg, "notebook failed")
            @test contains(msg, "notebook.jl")
            @test contains(msg, "sum(1, :b)")
            @test contains(msg, "Closest candidates are")
            @test contains(msg, "_foldl_impl")
        end
    end
end

@testset "pluto-docs-binding" begin
    text = """
        "This is a docstring"
        foo(x) = x
        """
    nb = Notebook([
        Cell(text),
    ])
    html, nb = notebook2html_helper(nb; use_distributed=false)

    @test !contains(html, "pluto-docs-binding")
end

@testset "benchmark-hack" begin
    # Related to https://github.com/fonsp/Pluto.jl/issues/1664
    nb = Notebook([
        Cell("using BenchmarkTools"),
        Cell("@benchmark sum(x)"),
        Cell("x = [1, 2]")
    ])
    html, nb = notebook2html_helper(nb)
    @test contains(html, "BenchmarkTools.Trial")
end

@testset "show_output_above_code" begin
    nb = Notebook([
        Cell("x = 1 + 1020"),
    ])
    hopts = HTMLOptions(; show_output_above_code=true)
    html, _ = notebook2html_helper(nb, hopts; use_distributed=false)
    lines = split(html, '\n')

    @test contains(lines[1], "1021")
    @test contains(lines[2], "1 + 1020")
end

@testset "with_terminal" begin
    nb = Notebook([
        Cell("using PlutoUI"),
        Cell("f(x) = Base.inferencebarrier(x);"),
        Cell("""
            with_terminal() do
                @code_warntype f(1)
            end
            """)
    ])
    html, _ = notebook2html_helper(nb)
    # Basically, this only tests whether `_patch_with_terminal` is applied.
    @test contains(html, """
""")
end


text:
function update!(ctx::T, data::U) where {T <: RIPEMD160_CTX,
                                         U <: Union{DenseArray{UInt8, 1},
                                                    NTuple{N, UInt8} where N}}
    UIntXXX = typeof(ctx.count)

    len = convert(UIntXXX, length(data))
    data_idx = convert(UIntXXX, 0)
    usedspace = ctx.count % bytes_per_block(T)

    while len - data_idx + usedspace >= bytes_per_block(T)
        copyto!(ctx.buffer, usedspace + 1,
                data,       data_idx + 1,
                bytes_per_block(T))

        transform!(ctx)

        ctx.count += bytes_per_block(T) - usedspace
        data_idx += bytes_per_block(T) - usedspace
        usedspace = convert(UIntXXX, 0)
    end

    if len > data_idx
        copyto!(ctx.buffer, usedspace + 1,
                data,       data_idx + 1,
                len - data_idx)
        ctx.count += len - data_idx
    end
    return nothing
end

function pad_remainder!(ctx::T) where {T <: RIPEMD160_CTX}

    @inbounds begin

        usedspace = ctx.count % bytes_per_block(T)

        if usedspace > 0
            usedspace += 1
            ctx.buffer[usedspace] = 0x80

            # do we have space for a UInt64?
            if usedspace <= bytes_per_block(T) - sizeof(ctx.count)
                # space for UInt64 so fill with 0x0 except the last UInt64
                for i = (usedspace + 1):(length(ctx.buffer) - sizeof(ctx.count))
                    ctx.buffer[i] = 0x0
                end
            else
                # no space for UInt64 fill out everything, transform and fill with
                # 0x0
                for i = (usedspace + 1):length(ctx.buffer)
                    ctx.buffer[i] = 0x0
                end
                transform!(ctx)
                for i = 1:bytes_per_block(T)
                    ctx.buffer[i] = 0x0
                end
            end
        else # usedspace == 0
            ctx.buffer[1] = 0x80
            for i = 2:bytes_per_block(T)
                ctx.buffer[i] = 0x0
            end
        end
    end
    return nothing
end

function digest!(ctx::T) where {T <: RIPEMD160_CTX}
    pad_remainder!(ctx)

    # bitcount_idx = div(words_per_block(T), sizeof(ctx.count)) + 1
    pbuf = Ptr{typeof(ctx.count)}(pointer(ctx.buffer))
    unsafe_store!(pbuf, UInt64(ctx.count * 8), 8)

    transform!(ctx)

    # TODO: is this safe for 0.7? the [1:end] should make a copy and get around
    # the lazy reinterpret.
    return reinterpret(UInt8, ctx.state)[1:digest_length(T)]
end


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

module UUIDs

using Random

import SHA

export UUID, uuid1, uuid4, uuid5, uuid_version

import Base: UUID

"""
    uuid_version(u::UUID) -> Int

Inspects the given UUID and returns its version
(see [RFC 4122](https://www.ietf.org/rfc/rfc4122)).

# Examples
```jldoctest
julia> uuid_version(uuid4())
4
```
"""
uuid_version(u::UUID) = Int((u.value >> 76) & 0xf)

# Some UUID namespaces provided in the appendix of RFC 4122
# https://tools.ietf.org/html/rfc4122.html#appendix-C
const namespace_dns  = UUID(0x6ba7b8109dad11d180b400c04fd430c8) # 6ba7b810-9dad-11d1-80b4-00c04fd430c8
const namespace_url  = UUID(0x6ba7b8119dad11d180b400c04fd430c8) # 6ba7b811-9dad-11d1-80b4-00c04fd430c8
const namespace_oid  = UUID(0x6ba7b8129dad11d180b400c04fd430c8) # 6ba7b812-9dad-11d1-80b4-00c04fd430c8
const namespace_x500 = UUID(0x6ba7b8149dad11d180b400c04fd430c8) # 6ba7b814-9dad-11d1-80b4-00c04fd430c8

"""
    uuid1([rng::AbstractRNG=GLOBAL_RNG]) -> UUID

Generates a version 1 (time-based) universally unique identifier (UUID), as specified
by RFC 4122. Note that the Node ID is randomly generated (does not identify the host)
according to section 4.5 of the RFC.

# Examples
```jldoctest; filter = r"[a-z0-9]{8}-([a-z0-9]{4}-){3}[a-z0-9]{12}"
julia> rng = MersenneTwister(1234);

julia> uuid1(rng)
UUID("cfc395e8-590f-11e8-1f13-43a2532b2fa8")
```
"""
function uuid1(rng::AbstractRNG=Random.GLOBAL_RNG)
    u = rand(rng, UInt128)

    # mask off clock sequence and node
    u &= 0x00000000000000003fffffffffffffff

    # set the unicast/multicast bit and version
    u |= 0x00000000000010000000010000000000

    # 0x01b21dd213814000 is the number of 100 nanosecond intervals
    # between the UUID epoch and Unix epoch
    timestamp = round(UInt64, time() * 1e7) + 0x01b21dd213814000
    ts_low = timestamp & typemax(UInt32)
    ts_mid = (timestamp >> 32) & typemax(UInt16)
    ts_hi = (timestamp >> 48) & 0x0fff

    u |= UInt128(ts_low) << 96
    u |= UInt128(ts_mid) << 80
    u |= UInt128(ts_hi) << 64

    UUID(u)
end

"""
    uuid4([rng::AbstractRNG=GLOBAL_RNG]) -> UUID

Generates a version 4 (random or pseudo-random) universally unique identifier (UUID),
as specified by RFC 4122.

# Examples
```jldoctest
julia> rng = MersenneTwister(1234);

julia> uuid4(rng)
UUID("196f2941-2d58-45ba-9f13-43a2532b2fa8")
```
"""
function uuid4(rng::AbstractRNG=Random.GLOBAL_RNG)
    u = rand(rng, UInt128)
    u &= 0xffffffffffff0fff3fffffffffffffff
    u |= 0x00000000000040008000000000000000
    UUID(u)
end

"""
    uuid5(ns::UUID, name::String) -> UUID

Generates a version 5 (namespace and domain-based) universally unique identifier (UUID),
as specified by RFC 4122.

!!! compat "Julia 1.1"
    This function requires at least Julia 1.1.

# Examples
```jldoctest
julia> rng = MersenneTwister(1234);

julia> u4 = uuid4(rng)
UUID("196f2941-2d58-45ba-9f13-43a2532b2fa8")

julia> u5 = uuid5(u4, "julia")
UUID("b37756f8-b0c0-54cd-a466-19b3d25683bc")
```
"""
function uuid5(ns::UUID, name::String)
    nsbytes = zeros(UInt8, 16)
    nsv = ns.value
    for idx in Base.OneTo(16)
        nsbytes[idx] = nsv >> 120
        nsv = nsv << 8
    end
    hash_result = SHA.sha1(append!(nsbytes, convert(Vector{UInt8}, codeunits(unescape_string(name)))))
    # set version number to 5
    hash_result[7] = (hash_result[7] & 0x0F) | (0x50)
    hash_result[9] = (hash_result[9] & 0x3F) | (0x80)
    v = zero(UInt128)
    #use only the first 16 bytes of the SHA1 hash
    for idx in Base.OneTo(16)
        v = (v << 0x08) | hash_result[idx]
    end
    return UUID(v)
end

end

max_stars_repo_name: musm/julia

max_issues_repo_name: Tominous/julia

max_forks_repo_name: Tominous/julia


text:
"""
	filter1d(cmd0::String="", arg1=[], kwargs...)

Time domain filtering of 1-D data tables.

Full option list at [`filter1d`](http://gmt.soest.hawaii.edu/doc/latest/filter1d.html)

Parameters
----------

- **F** : **filter_type** : -- Str --   Flags = type width[modifiers]

    Sets the filter type. Choose among convolution and non-convolution filters. Append the
    filter code followed by the full filter width in same units as time column.
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/filter1d.html#f)
- **D** : **inc** : -- Number --        Flags = increment

    ``increment`` is used when series is NOT equidistantly sampled. Then increment will be the abscissae resolution.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/filter1d.html#d)
- **E** : **ends** : -- Bool or [] --

    Include Ends of time series in output. Default loses half the filter-width of data at each end.
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/filter1d.html#e)
- **N** : **time_col** : -- Int --      Flags = t_col

    Indicates which column contains the independent variable (time). The left-most column
    is # 0, the right-most is # (n_cols - 1). [Default is 0].
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/filter1d.html#n)
- **Q** : **quality** : -- Number --    Flags = q_factor

    Assess Quality of output value by checking mean weight in convolution. Enter q_factor between 0 and 1.
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/filter1d.html#q)
- **S** : **symetry** : -- Number --    Flags = symmetry_factor

    Checks symmetry of data about window center. Enter a factor between 0 and 1.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/filter1d.html#s)
- **T** : **equi_space** : -- List or Str --     Flags = [min/max/]inc[+a|n]

    Make evenly spaced time-steps from min to max by inc [Default uses input times].
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/filter1d.html#t)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_swap_xy)
"""
function filter1d(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("filter1d", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd = parse_V_params("", d)
	cmd, = parse_b(cmd, d)
	cmd, = parse_d(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_g(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_o(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'F', d, [:F :filter_type])
	cmd = add_opt(cmd, 'D', d, [:D :inc])
	cmd = add_opt(cmd, 'E', d, [:E :ends])
	cmd = add_opt(cmd, 'N', d, [:N :time_col])
	cmd = add_opt(cmd, 'Q', d, [:Q :quality])
	cmd = add_opt(cmd, 'S', d, [:S :symetry])
	cmd = add_opt(cmd, 'T', d, [:T :equi_space])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "filter1d", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
filter1d(arg1=[], cmd0::String=""; kw...) = filter1d(cmd0, arg1; kw...)


text:
"""
    HullReach(resolution::Float64, tight::Bool)

HullReach performs over-approximated reachability analysis to compute the over-approximated output reachable set for a network.

# Problem requirement
1. Network: any depth, any activation that is monotone
2. Input: `Hyperrectangle`
3. Output: `HPolytope`

# Return
`BasicResult`

# Property
Sound but not complete.
"""
@with_kw struct HullReach
    resolution::Float64 = 1.0
    tight::Bool         = false
end

# This is the main function
function solve(solver::HullReach, problem::Problem) #original
    result = true
    delta = solver.resolution
    lower, upper = low(problem.input), high(problem.input)
    n_hypers_per_dim = BigInt.(max.(ceil.(Int, (upper-lower) / delta), 1))

    # preallocate work arrays
    local_lower, local_upper, CI = similar(lower), similar(lower), similar(lower)
    for i in 1:prod(n_hypers_per_dim)
        n = i
        hull = false
        for j in firstindex(CI):lastindex(CI)
            n, CI[j] = fldmod1(n, n_hypers_per_dim[j])
            if CI[j] == 1 || CI[j] == n_hypers_per_dim[j]
                hull = true
            end
        end
        if hull
            @. local_lower = lower + delta * (CI - 1)
            @. local_upper = min(local_lower + delta, upper)
            hyper = Hyperrectangle(low = local_lower, high = local_upper)
            reach = forward_network(solver, problem.network, hyper)
            if !issubset(reach, problem.output)
                result = false
            end
        end
    end
    if result
        return BasicResult(:holds)
    end
    return BasicResult(:violated)
end

# This function is called by forward_network
function forward_layer(solver::HullReach, L::Layer, input::Hyperrectangle)
    (W, b, act) = (L.weights, L.bias, L.activation)
    center = zeros(size(W, 1))
    gamma  = zeros(size(W, 1))
    for j in 1:size(W, 1)
        node = Node(W[j,:], b[j], act)
        center[j], gamma[j] = forward_node(solver, node, input)
    end
    return Hyperrectangle(center, gamma)
end

function forward_node(solver::HullReach, node::Node, input::Hyperrectangle)
    output    = node.w' * input.center + node.b
    deviation = sum(abs.(node.w) .* input.radius)
    β    = node.act(output)  # TODO expert suggestion for variable name. beta? β? O? x?
    βmax = node.act(output + deviation)
    βmin = node.act(output - deviation)
    if solver.tight
        return ((βmax + βmin)/2, (βmax - βmin)/2)
    else
        return (β, max(abs(βmax - β), abs(βmin - β)))
    end
end


text:
export xNES, xnes_init, xnes_populate, xnes_generation, populate, evaluate, generation

# adapted from
# https://github.com/francescoalemanno/NaturalES.jl/blob/master/src/exponential_nes.jl

import Cambrian: populate, evaluate, generation

mutable struct xNESState <: ESState
    μ::Array{Float64}
    σ::Float64
    B::Array{Float64}
    u::Array{Float64}
    s::Array{Float64}
end

function xNESState(d::Int, n::Int;
                   μ::Array{Float64}=randn(d),
                   σ::Float64=1.0,
                   A=diagm(fill(σ, d)))
    σ = abs(det(A))^(1/d)
    B = A ./ σ
    u = zeros(n)
    u = max.(0, log(n/2+1) .- log.(1:n))
    u .= u ./ sum(u) .- 1/n
    s = randn(d, n)
    xNESState(μ, σ, B, u, s)
end

function xnes_config(d::Int; n_population = 4+ceil(Int, log(3*d)),
                     ημ = 1.0, ησ = (9 + 3 * log(d)) / (5 * d * sqrt(d)),
                     ηB = ησ)
    (n_population = n_population, ημ = ημ, ησ = ησ, ηB = ηB)
end

mutable struct xNES <: Cambrian.AbstractEvolution
    config::NamedTuple
    logger::CambrianLogger
    population::Array{AbstractESIndividual}
    elites::Array{AbstractESIndividual}
    state::xNESState
    fitness::Function
    gen::Int
end

function xnes_init(cfg::NamedTuple, state::ESState; T::Type=ESIndividual)
    population = Array{AbstractESIndividual}(undef, cfg.n_population)
    for i in 1:cfg.n_population
        genes = state.μ .+ state.σ .* (state.B * view(state.s, :, i))
        population[i] = T(genes, -Inf*ones(cfg.d_fitness))
    end
    population
end

function xNES(cfg::NamedTuple, fitness::Function, state::xNESState; T::Type=ESIndividual ,logfile=string("logs/", cfg.id, ".csv"))
    logger = CambrianLogger(logfile)
    population = xnes_init(cfg, state, T=T)
    elites = deepcopy([population[i] for i in 1:cfg.n_elite])
    xNES(cfg, logger, population, elites, state, fitness, 0)
end

"""
Exponential Natural Evolutionary Strategies default constructor
Will use a random starting point using N(0, 1)
If cfg contains keys from snes_config, these will be used
To provide initial state, see xNES(cfg, fitness, state)
"""
function xNES(cfg::NamedTuple, fitness::Function;  T::Type=ESIndividual, logfile=string("logs/", cfg.id, ".csv"))
    logger = CambrianLogger(logfile)
    cfg = merge(xnes_config(cfg.n_genes), cfg)
    state = xNESState(cfg.n_genes, cfg.n_population)
    xNES(cfg, fitness, state; T=T, logfile=logfile)
end

"generate next population"
function xnes_populate(e::xNES)
    for i in eachindex(e.population)
        e.population[i].genes .= e.state.μ .+ e.state.σ .* (e.state.B * view(e.state.s, :, i))
        e.population[i].fitness .= -Inf
    end
end

"update NES state, called after populate and evaluate"
function xnes_generation(e::xNES)
    d = e.config.n_genes
    n = e.config.n_population

    # copy population information
    F = zeros(n)
    for i in eachindex(e.population)
        F[i] = -e.population[i].fitness[1]
    end
    idx = sortperm(F)

    # compute gradients
    ∇δ = zeros(d)
    ∇M = zeros(d,d)
    for i in 1:n
        j = idx[i]
        for k2 in 1:d
            ∇δ[k2] += e.state.u[i] * e.state.s[k2,j]
            for k1 in 1:d
                ∇M[k1,k2] += e.state.u[i] * (e.state.s[k1,j] * e.state.s[k2,j] - (k1==k2))
            end
        end
    end
    ∇σ = tr(∇M) / d
    for i in 1:d
        ∇M[i,i] -=  ∇σ
    end

    # update state variables
    e.state.μ .+= e.config.ημ .* e.state.σ .* (e.state.B * ∇δ)
    e.state.σ *= exp(e.config.ησ/2 * ∇σ)
    e.state.B = e.state.B * exp(e.config.ηB/2 .* ∇M)
    randn!(e.state.s)
    Cambrian.elites_generation(e)
end

populate(e::xNES) = xnes_populate(e)
evaluate(e::xNES) = fitness_evaluate(e, e.fitness)
generation(e::xNES) = xnes_generation(e)


text:
function get_coefficients(obj::FieldShortPeriodTerms, arg0::FieldAbsoluteDate, arg1::Set)
    return jcall(obj, "getCoefficients", Map, (FieldAbsoluteDate, Set), arg0, arg1)
end

function get_coefficients_key_prefix(obj::FieldShortPeriodTerms)
    return jcall(obj, "getCoefficientsKeyPrefix", JString, ())
end

function value(obj::FieldShortPeriodTerms, arg0::FieldOrbit)
    return jcall(obj, "value", Vector{RealFieldElement}, (FieldOrbit,), arg0)
end


text:
using CovarianceMatrices
using Random
using GLM
using DataFrames
using BenchmarkTools

Random.seed!(1)

u = zeros(6000*50)
for j in 1:2999
    u[j+1] = 0.97*u[j] + randn()
end

df = DataFrame(y = randn(300*50) .+ u[15001:30000])

df2 = DataFrame(y = sqrt(2).*randn(250))
for j in Symbol.("x".*string.(collect(1:5)))
    df2[j] = randn(250)
end
df2[:cluster] = repeat(1:50, inner = [5])

frm = @formula(y ~ x1 + x2 + x3 + x4 + x5)
lm1 = glm(frm, df, Normal(), IdentityLink())
lm2 = glm(frm, df2, Normal(), IdentityLink())

SUITE = BenchmarkGroup()
SUITE["HAC Andrews"] = BenchmarkGroup()
SUITE["HAC Newey"] = BenchmarkGroup()
SUITE["HAC Fixed(10)"] = BenchmarkGroup()
SUITE["HAC Fixed(30)"] = BenchmarkGroup()

SUITE["CRHC (large)"] = BenchmarkGroup()
SUITE["CRHC (small)"] = BenchmarkGroup()

SUITE["HAC Andrews"]["Parzen"] = @benchmarkable vcov($ParzenKernel(), lm1)
SUITE["HAC Andrews"]["Truncated"] = @benchmarkable vcov($TruncatedKernel(), lm1)
SUITE["HAC Andrews"]["Bartlett"] = @benchmarkable vcov($BartlettKernel(), lm1)

SUITE["HAC Newey"]["Parzen"] = @benchmarkable vcov($ParzenKernel(NeweyWest), lm1)
SUITE["HAC Newey"]["Bartlett"] = @benchmarkable vcov($BartlettKernel(NeweyWest), lm1)

SUITE["HAC Fixed(10)"]["Parzen"] = @benchmarkable vcov($ParzenKernel(10), lm1)
SUITE["HAC Fixed(10)"]["Truncated"] = @benchmarkable vcov($TruncatedKernel(10), lm1)
SUITE["HAC Fixed(10)"]["Bartlett"] = @benchmarkable vcov($BartlettKernel(10), lm1)

SUITE["HAC Fixed(30)"]["Parzen"] = @benchmarkable vcov($ParzenKernel(30), lm1)
SUITE["HAC Fixed(30)"]["Truncated"] = @benchmarkable vcov($TruncatedKernel(30), lm1)
SUITE["HAC Fixed(30)"]["Bartlett"] = @benchmarkable vcov($BartlettKernel(30), lm1)

k0 = CRHC0(df[!,:cluster])
k2 = CRHC2(df[!,:cluster])
k3 = CRHC3(df[!,:cluster])

k20 = CRHC0(df2[!,:cluster])
k22 = CRHC2(df2[!,:cluster])
k23 = CRHC3(df2[!,:cluster])

SUITE["CRHC (large)"]["CRHC0"] = @benchmarkable vcov($k0, lm1)
SUITE["CRHC (large)"]["CRHC2"] = @benchmarkable vcov($k2, lm1)
SUITE["CRHC (large)"]["CRHC3"] = @benchmarkable vcov($k3, lm1)

SUITE["CRHC (small)"]["CRHC0"] = @benchmarkable vcov($k20, lm2)
SUITE["CRHC (small)"]["CRHC2"] = @benchmarkable vcov($k22, lm2)
SUITE["CRHC (small)"]["CRHC3"] = @benchmarkable vcov($k23, lm2)


text:
using Test, Clines, LinearFractionalTransformations

@testset "Construction" begin
    C = Circle(1 + 0im, im, -im)
    @test radius(C) == 1
    @test center(C) == 0

    D = Circle(0, 0, 1)
    @test C == D

    E = Circle(0im, 1)
    @test C == E

    P = three_points(E)
    @test C == Cline(P...)

    L = Line(1, 2)
    @test L == Line(-1, 3)
    @test L == Cline(1, 5, Inf)

    LL = dilate(L)
    @test LL == L
end

@testset "Properties" begin
    C = Circle(2, 3, 2)
    @test area(C) == 4π
    @test circumference(C) == 4π

    L = Cline(0, 1 + im, 2 + 2im)
    @test slope(L) == 1
end

@testset "Set Ops" begin
    C = Circle(0im, 6)
    D = Circle(4 + 0im, 2)
    @test D ⊆ C
    S = C ∩ D
    @test length(S) == 1

    D = Circle(4 + 0im, 3)
    S = C ∩ D
    @test length(S) == 2

    D = Circle(4 + 0im, 1)
    S = C ∩ D
    @test isempty(S)
end

@testset "Intersection" begin
    C = Circle(0, 0, 1)
    L = Line(0, 3 + 4im)
    @test length(C ∩ L) == 2
    D = Circle(2, 0, 1)
    @test length(C ∩ D) == 1

    D = Circle(5, 5, 1)
    @test C ∩ D == Set{Complex{Float64}}()

    X = Line(0, 1)
    Y = Line(im, 2 + im)
    @test length(X ∩ Y) == 0

    Y = Cline(-5, 2, Inf)
    @test X == Y
end

@testset "LFT" begin
    C = Circle(1 + im, 3)
    F = LFT()
    D = F(C)
    @test abs(center(C) - center(D)) <= get_tolerance()
    @test abs(radius(C) - radius(D)) <= get_tolerance()

    D = Circle(2 - im, 5)
    F = LFT(C, D)
    @test F(C) == D
    @test (inv(F))(D) == C

    F = LFT(C)
    X = Line(0, 1)
    @test (inv(F))(X) == C

end

@testset "Inversion" begin
    X = Line(0, 1)
    z = 3 - 2im
    @test inv(X, z) == conj(z)

    C = Circle(2, 3, 5)
    L = Line(0, 1 - im)

    LL = inv(C, L)
    @test inv(C, LL) == L

    F = inv(L)
    @test F(F(C)) == C

end

@testset "Collinearity" begin
    pts = 2 * collect(1:3) - im * collect(1:3)
    @test collinear(pts...)
    @test isa(Cline(pts...), Line)

    C = Circle(4, 6, 10)
    PP = three_points(C)
    @test !collinear(PP...)
end

@testset "Kiss" begin
    pts = [0, 2im, 3 + im]
    CC = kiss(pts...)
    @test length(CC[1] ∩ CC[2]) == 1
    S = soddy(CC...)
    @test length(CC[3] ∩ S) == 1
end


text:
using Serialization

struct RunnerProcResult
    exitcode::Int64
    result::Union{RunnerOutput, Nothing}
    stdout::AbstractString
    stderr::AbstractString
end

RunnerProcResult(s::Int64, r::Any, e::Nothing) = RunnerProcResult(s, r, "")

function readpipe!(p::Pipe)
    close(p.in)
    return read(p, String)
end

"""
    launch_runner_proc(testsuite, submissionfile, functionname[; debug])

Launch a runner subprocess and get the resulting `RunnerProcResult`.

**Note:** This function uses `deserialize` and so bugs may cause the Julia
process to crash without warning (no error, no backtrace, etc.) if the
serialized input is invalid.
"""
function launch_runner_proc(
    testsuite::TestSuite,
    submissionfile::String,
    functionname::Union{AbstractString, Symbol},
    ;
    debug::String=get(ENV, "JULIA_DEBUG", ""),
)
    functionname = Symbol(functionname)
    input_path, output_path = tempname(), tempname()
    runner_input = RunnerInput(testsuite, submissionfile, functionname, output_path)

    # Julia makes it *very* hard to actually get the id of a subprocess, so
    # we generate a random string here to be able to track what logs are comming
    # from what process. Note that we use an alphabetic rather than a number to
    # prevent confusion with the process id.
    log_id = String(rand('a':'z', 6))
    try
        open(io -> serialize(io, runner_input), input_path, "w")

        @debug "Launching runner subprocess." submissionfile JULIA_EXE
        submissionfile = realpath(submissionfile)
        # proc_in, proc_out, proc_err = Pipe(), Pipe(), Pipe()
        stdout_pipe, stderr_pipe = Pipe(), Pipe()
        # jl_cmd = "using ComparativeAutograder: runner_main; runner_main()"
        jl_env = copy(ENV)
        jl_env["JULIA_DEBUG"] = debug
        proc = run(pipeline(
            Cmd(`$(Base.julia_cmd()) -e $RUNNER_CMD $input_path`, env=jl_env),
            stdout=stdout_pipe,
            stderr=stderr_pipe,
        ); wait=false)
        @debug "Spawned runner process."

        stdout_task = @async readpipe!(stdout_pipe)
        stderr_task = @async readpipe!(stderr_pipe)

        # Wait for the subprocess to complete.
        wait(proc)
        status = proc.exitcode
        @debug "Runner process ($log_id) exited." status stdout_task stderr_task
        stdout_str, stderr_str = fetch.([stdout_task, stderr_task])

        # Fetch the RunnerOutput.
        runner_output = nothing
        try
            runner_output = open(io -> deserialize(io), output_path, "r")
            typeassert(runner_output, RunnerOutput)
        catch exc
            @error "An error occurred while trying to read the runner output file." exc
        end

        return RunnerProcResult(status, runner_output, stdout_str, stderr_str)
    finally
        rm.([input_path, output_path], force=true)
    end
end


text:
using CUDA
using CUDA.CUSPARSE
using ExaPF
using KernelAbstractions
using LinearAlgebra
using Random
using SparseArrays
using Test
using TimerOutputs
using FiniteDiff

import ExaPF: AutoDiff
Random.seed!(2713)

const INSTANCES_DIR = joinpath(dirname(@__FILE__), "..", "data")

@testset "ExaPF.PowerSystem" begin
    # Test PowerNetwork object and parsers
    include("powersystem.jl")
end

@testset "ExaPF.LinearSolvers" begin
    include("iterative_solvers.jl")
end

@testset "Polar formulation" begin
    # Test that behavior matches Matpower
    include("Polar/matpower.jl")
    # Test polar formulation
    include("Polar/polar_form.jl")
    # Autodiff's Jacobians
    include("Polar/autodiff.jl")
    # Reduced gradient
    include("Polar/gradient.jl")
    # Reduced Hessian
    include("Polar/hessian.jl")
end

@testset "Optimization evaluators" begin
    # Resolution of powerflow equations with NLPEvaluators
    include("Evaluators/powerflow.jl")
    # Test generic API
    include("Evaluators/interface.jl")
    # Test more in-depth each evaluator
    include("Evaluators/reduced_evaluator.jl")
    include("Evaluators/proxal_evaluator.jl")
    include("Evaluators/auglag.jl")
    include("Evaluators/MOI_wrapper.jl")
    # Test basic reduced gradient algorithm
    include("Evaluators/test_rgm.jl")
end

@testset "Documentation" begin
    include("quickstart.jl")
end

@testset "Benchmark script" begin
    empty!(ARGS)
    push!(ARGS, "KrylovBICGSTAB")
    push!(ARGS, "CPU")
    push!(ARGS, "case300.m")
    include("../benchmark/benchmarks.jl")
    @test convergence.has_converged
end


text:
# identity operators: joEye

export joEye
"""
    julia> op = joEye(m;[DDT=...,][RDT=...,][name=...])

Identity matrix - square

# Signature

    joEye(m::Integer;DDT::DataType=joFloat,RDT::DataType=DDT,name::String="joEye")

# Arguments

- `m`: sizes
- keywords
    - `DDT`: domain data type
    - `RDT`: range data type
    - `name`: custom name

# Examples

    joEye(m)

examples with DDT/RDT

    joEye(m; DDT=Float32)
    joEye(m; DDT=Float32,RDT=Float64)

"""
joEye(m::Integer;DDT::DataType=joFloat,RDT::DataType=DDT,name::String="joEye") =
    joMatrix{DDT,DDT}("joEye",m,m,
        v1->jo_convert(RDT,v1,false),
        v2->jo_convert(DDT,v2,false),
        v3->jo_convert(DDT,v3,false),
        v4->jo_convert(RDT,v4,false),
        v5->jo_convert(DDT,v5,false),
        v6->jo_convert(RDT,v6,false),
        v7->jo_convert(RDT,v7,false),
        v8->jo_convert(DDT,v8,false)
        )
"""
    julia> op = joEye(m,n;[DDT=...,][RDT=...,][name=...])

Identity matrix - rectangular

# Signature

    joEye(m::Integer,n::Integer;DDT::DataType=joFloat,RDT::DataType=DDT,name::String="joEye")

# Arguments

- `m,n`: sizes
- keywords
    - `DDT`: domain data type
    - `RDT`: range data type
    - `name`: custom name

# Examples

    joEye(m,n)

examples with DDT/RDT

    joEye(m,n; DDT=Float32)
    joEye(m,n; DDT=Float32,RDT=Float64)

"""
joEye(m::Integer,n::Integer;DDT::DataType=joFloat,RDT::DataType=DDT,name::String="joEye") =
    joMatrix{DDT,DDT}("joEye",m,n,
        v1->jo_convert(RDT,jo_speye(eltype(v1),m,n)*v1,false),
        v2->jo_convert(DDT,jo_speye(eltype(v2),n,m)*v2,false),
        v3->jo_convert(DDT,jo_speye(eltype(v3),n,m)*v3,false),
        v4->jo_convert(RDT,jo_speye(eltype(v4),m,n)*v4,false),
        @joNF, @joNF, @joNF, @joNF
        )
        #v5->jo_speye(EDT,m,n)\v5, v6->jo_speye(EDT,n,m)\v6, v7->jo_speye(EDT,n,m)\v7, v8->jo_speye(EDT,m,n)\v8
        #v5->jo_speye(EDT,n,m)*v5, v6->jo_speye(EDT,m,n)*v6, v7->jo_speye(EDT,m,n)*v7, v8->jo_speye(EDT,n,m)*v8


text:
# This file was generated, do not modify it. # hide
using MLJ
import RDatasets: dataset
import DataFrames: DataFrame, select, Not, describe
using Random
MLJ.color_off() # hide

data = dataset("datasets", "USArrests")
names(data)


text:
function repl_runcode_request(conn, params::ReplRunCodeRequestParams)
    source_filename = params.filename
    code_line = params.line
    code_column = params.column
    source_code = params.code
    mod = params.mod

    resolved_mod = try
        module_from_string(mod)
    catch err
        # maybe trigger error reporting here
        Main
    end

    show_code = params.showCodeInREPL
    show_result = params.showResultInREPL

    rendered_result = nothing

    hideprompt() do
        if isdefined(Main, :Revise) && isdefined(Main.Revise, :revise) && Main.Revise.revise isa Function
            let mode = get(ENV, "JULIA_REVISE", "auto")
                mode == "auto" && Main.Revise.revise()
            end
        end
        if show_code
            for (i,line) in enumerate(eachline(IOBuffer(source_code)))
                if i==1
                    printstyled("julia> ", color=:green)
                    print(' '^code_column)
                else
                    # Indent by 7 so that it aligns with the julia> prompt
                    print(' '^7)
                end

                println(line)
            end
        end

        withpath(source_filename) do
            res = try
                ans = Base.invokelatest(include_string, resolved_mod, '\n'^code_line * ' '^code_column *  source_code, source_filename)
                @eval Main ans = $(QuoteNode(ans))
            catch err
                EvalError(err, catch_backtrace())
            end

            if show_result
                if res isa EvalError
                    Base.display_error(stderr, res.err, res.bt)
                elseif res !== nothing && !ends_with_semicolon(source_code)
                    Base.invokelatest(display, res)
                end
            else
                try
                    Base.invokelatest(display, InlineDisplay(), res)
                catch err
                    if !(err isa MethodError)
                        printstyled(stderr, "Display Error: ", color = Base.error_color(), bold = true)
                        Base.display_error(stderr, err, catch_backtrace())
                    end
                end
            end

            rendered_result = safe_render(res)
        end
    end
    return rendered_result
end

function repl_getvariables_request(conn, params::Nothing)
    M = Main
    variables = []
    clear_lazy()

    for n in names(M, all=true, imported=true)
        !isdefined(M, n) && continue
        Base.isdeprecated(M, n) && continue

        x = getfield(M, n)
        x === vscodedisplay && continue
        x === VSCodeServer && continue
        x === Main && continue

        s = string(n)
        startswith(s, "#") && continue
        try
            push!(variables, treerender(SubTree(s, wsicon(x), x)))
        catch err
            printstyled("Internal Error: ", bold = true, color = Base.error_color())
            Base.display_error(err, catch_backtrace())
        end
    end

    return variables
end

function repl_getlazy_request(conn, params::Int)
    res = get_lazy(params)

    return res
end

function repl_showingrid_notification(conn, params::String)
    try
        var = Core.eval(Main, Meta.parse(params))

        Base.invokelatest(internal_vscodedisplay, var)
    catch err
        Base.display_error(err, catch_backtrace())
    end
end

function repl_loadedModules_request(conn, params::Nothing)
    res = string.(collect(get_modules()))

    return res
end

function repl_isModuleLoaded_request(conn, params::String)
    is_loaded = is_module_loaded(params)

    return is_loaded
end

function repl_startdebugger_request(conn, params::String, crashreporting_pipename)
    hideprompt() do
        debug_pipename = params
        try
            DebugAdapter.startdebug(debug_pipename)
        catch err
            DebugAdapter.global_err_handler(err, catch_backtrace(), crashreporting_pipename, "Debugger")
        end
    end
end

max_stars_repo_name: fonsp/julia-vscode

max_issues_repo_name: fonsp/julia-vscode

max_forks_repo_name: fonsp/julia-vscode


text:
module CrystalMod

using LinearAlgebra
using Printf
using ..Atomdata:atoms
using GZip

using ..ThreeBodyTB:global_length_units
using ..ThreeBodyTB:convert_length
using ..ThreeBodyTB:convert_force
using ..ThreeBodyTB:convert_stress
using ..ThreeBodyTB:Ang

export crystal
export makecrys
export generate_supercell
export generate_random_distortion
export write_poscar
export write_efs

#holds data
"""
    mutable struct crystal{T}

Holds basic crystal structure information, type T. Use `makecrys` to easily construct.

Note: you can create supercells like

```julia-repl
julia> c = makecrys([5.0 0 0; 0 5.0 0; 0 0 5.0], [0.0 0.0 0.0], ["H"])
A1=     5.00000  0.00000  0.00000
A2=     0.00000  5.00000  0.00000
A3=     0.00000  0.00000  5.00000

H    0.00000  0.00000  0.00000

julia> c*[2,2,2]
A1=     10.00000  0.00000  0.00000
A2=     0.00000  10.00000  0.00000
A3=     0.00000  0.00000  10.00000

H    0.00000  0.00000  0.00000
H    0.00000  0.00000  0.50000
H    0.00000  0.50000  0.00000
H    0.00000  0.50000  0.50000
H    0.50000  0.00000  0.00000
H    0.50000  0.00000  0.50000
H    0.50000  0.50000  0.00000
H    0.50000  0.50000  0.50000
```

# Holds
- `A::Array{T,2}` 3 × 3 lattice vectors, Bohr (atomic units) internally.
- `coords::Array{T,2}` num_atoms × 3  atomic positions, fractional units.
- `types::Array{String,1}` atomic names, like `"H"` or `"Zn"`.
- `types::Array{Symbol,1}` atomic names, but julia Symbols like `:H` or `:Zn`, for nominally faster internal evaluation.
- `nat::Int64` number of atoms.
"""
mutable struct crystal{T}

    A::Array{T,2}
    coords::Array{T,2}
    types::Array{String,1}
    stypes::Array{Symbol,1}
    nat::Int64

end

#printing
Base.show(io::IO, c::crystal) = begin

    println(io, "Units: "* global_length_units)
    println(io)
    if global_length_units == "Å"
        At = c.A * Ang
    else
        At = c.A
    end

    for i in 1:3
            @printf(io, "A%.1i=     %.5f  %.5f  %.5f\n", i, At[i,1], At[i,2], At[i,3])
    end
    @printf(io, "\n")
    for i in 1:c.nat
        @printf(io, "%-3s  %.5f  %.5f  %.5f\n", c.types[i], c.coords[i,1], c.coords[i,2], c.coords[i,3])
    end

end

"""
    function print_with_force_stress(c::crystal, force, stress)

pretty io function
"""
function print_with_force_stress(c::crystal, force, stress)

#    println(io, "Units: "* global_length_units)
    #    println(io)

    force = convert_force(force)
    stress = convert_stress(stress)

    if global_length_units == "Å"
        At = c.A * Ang
    else
        At = c.A
    end

    for i in 1:3
            @printf( "A%.1i=     %+.5f  %+.5f  %+.5f  |  %+.5f  %+.5f  %+.5f\n", i, At[i,1], At[i,2], At[i,3], stress[i,1],stress[i,2], stress[i,3] )
    end
    @printf( "\n")
    for i in 1:c.nat
        @printf( "%-3s  %+.5f  %+.5f  %+.5f |  %+.5f  %+.5f  %+.5f\n", c.types[i], c.coords[i,1], c.coords[i,2], c.coords[i,3], force[i,1],force[i,2],force[i,3])
    end

end

Base.:+(c1::crystal, c2::crystal) = begin

    if c1.nat != c2.nat
        println("nat not same ", c1.nat, c2.nat)
        return c1
    end
    A = c1.A .+ c2.A
    coords = c1.coords .+ c2.coords

    types = deepcopy(c1.types)

    return makecrys(A, coords, types, units="Bohr")
end

Base.:(==)(c1::crystal, c2::crystal) = begin

    if c1.nat != c2.nat
        return false
    end

    if c1.types != c2.types
        return false
    end

    if sum(abs.(c1.A - c2.A)) > 1e-11
        return false
    end

    if sum(abs.(c1.coords - c2.coords)) > 1e-11
        return false
    end

    return true

end

Base.:*(a::Real, c::crystal) = begin

    A = a * c.A
    coords = deepcopy( c.coords)

    types = deepcopy(c.types)

    return makecrys(A, coords, types, units="Bohr")
end

Base.:*(c::crystal,a::Real) = begin
    return a*c
end

Base.:*(a::Array{Int}, c::crystal) = begin

    A_new = zeros(3,3)
    A_new[1,:] = c.A[1,:] *a[1]
    A_new[2,:] = c.A[2,:] *a[2]
    A_new[3,:] = c.A[3,:] *a[3]

    n = prod(a)
    coords_new = zeros(n*c.nat, 3)
    types_new = String[]

    counter = 0
    for x = 1:a[1]
        for y = 1:a[2]
            for z = 1:a[3]
                for i = 1:c.nat
                    counter +=1
                    coords_new[counter,:] = (c.coords[i,:] + [x-1, y-1, z-1]) ./ [a[1], a[2], a[3]]
                    push!(types_new, c.types[i])
                end
            end
        end
    end

    return makecrys(A_new, coords_new, types_new, units="Bohr")
end

Base.:*(c::crystal,  a::Array{Int} ) = begin
    return a*c
end

#setup crystal
"""
    makecrys(A,coords,types; units=missing)

Return a crystal object from 3×3 lattice, nat × 3 coords in crystal units, and nat element strings.

Will use units set by `set_units`, with default to Ang. Can override with `units`.

Note: also export-ed directly from ThreeBodyTB for convenience

```julia-repl
julia> makecrys([10.0 0 0; 0 10.0 0; 0 0 10.0], [0.0 0.0 0.0], ["H"])
A1=     10.00000  0.00000  0.00000
A2=     0.00000  10.00000  0.00000
A3=     0.00000  0.00000  10.00000

H    0.00000  0.00000  0.00000
```
"""
function makecrys(A,coords,types; units=missing)
    T = typeof(A[1,1])

    if !ismissing(units)
        if units == "A" || units == "Ang" || units == "Ang." || units == "Angstrom" || units == "a" || units == "ang" || units == "ang." || units == "angstrom" || units == "Å"
            factor = 1.0 / Ang
        else
            factor = 1.0
        end
    else
        if global_length_units == "Å"
            factor = 1.0 / Ang
        else
            factor = 1.0
        end
    end

    #entering integer coords or A messes everything up later
    if T == Int64
        coords = Float64.(coords)
        T = typeof(coords[1,1])
    end
    if typeof(A[1,1]) == Int64
        A = Float64.(A)
    end

    A = A * factor

    if length(size(types)) == 2  # if are accidently given a 1 x nat types array
        types = types[:]
    end

    if typeof(types[1]) == Symbol
        types = String.(types)
    end

    nat = length(types)

    if nat != size(coords,1) #try to fix
        coords = coords'
    end

    if nat != size(coords,1)

        println(types)
        println(coords)
        error("Error initilzing crystal, nat doesn't match, nat from types: ", nat," size(coords,1): ", size(coords,1)  )

        return
    end

    if (3,3) != size(A)
        println(A)
        error("Error initilzing crystal, A wrong ", size(A))

        return
    end

    if 3 != size(coords,2)
        println(coords)
        error("Error initilzing crystal, coords wrong ", size(coords))

        return
    end

    A = convert(Array{T,2}, A)
    coords = convert(Array{T,2}, coords)

    if size(coords)[1] != size(types)[1]
        error("Error making crys, types and coords sizes don't match")
        return
    end

    nat = size(coords)[1]
    stypes = Symbol.(types)
    return crystal{T}(A, coords, types, stypes, nat)
end

"""
    makecrys(filename::String)

Read filename, return crystal object.
File can be POSCAR, or simple quantum espresso inputfile

The entire file can be in the string instead of the filename. The code decides which is which by looking for newlines"
"""
function makecrys(filename::String)

    if !occursin("\n", filename)
        if !isfile(filename)
            println("tried to open $filename to load crystal, file does not exist")
            throw(error("noFile"))
        end

        f = gzopen(filename, "r")
        lines = readlines(f)
        close(f)

        c = makecrys(lines)
        return c
    else
        println("newlines found, interpret as string containing file instead of filename")
        lines = String.(split(filename, "\n"))
        c = makecrys(lines)
    end

end

"""
    makecrys(lines::Array{String,1})

Read string array, return crystal object.
string array can be POSCAR, or simple quantum espresso inputfile
"""
function makecrys(lines::Array{String,1})

    intype="POSCAR"
    for line in lines
        sp = split(line)
        if size(sp)[1] >= 1 && sp[1] == "&control"
            intype="QE"
        end
        if  size(sp)[1] >= 1 && sp[1] == "A1="
            intype="MYFORMAT"
        end
    end

    if intype == "POSCAR"

        A, coords, types = parsePOSCAR(lines)

    elseif intype == "QE"

        A, coords, types = parseQEinput(lines)

    elseif intype == "MYFORMAT"

        A, coords, types = parseMYFORMATinput(lines)
        if global_length_units=="Å"
            A = A / Ang
        end

    end

    c = makecrys(A, coords,types, units="Bohr")
    return c

    end

function parseARRfloat(sp)

    return map(x->parse(Float64,x),sp)

end

function parseARRint(sp)

    return map(x->parse(Int64,x),sp)

end

function parseMYFORMATinput(lines)

    A = zeros(3,3)

    coords = zeros(1000, 3)
    types = []
    nat = 0
    for line in lines
        sp = split(line)
        if length(sp) >= 1
            if sp[1] == "A1="
                A[1,:] = parseARRfloat(sp[2:4])
            elseif sp[1] == "A2="
                A[2,:] = parseARRfloat(sp[2:4])
            elseif sp[1] == "A3="
                A[3,:] = parseARRfloat(sp[2:4])
            elseif length(sp) >= 4
                push!(types, String(sp[1]))
                nat += 1
                coords[nat,:] = parseARRfloat(sp[2:4])
            end
        end
    end

    coords = coords[1:nat, :]

    return A, coords, types

end

"""
    function parsePOSCAR(lines)

Parse a POSCAR from VASP

Called by makecrys, doesn't need to be called directly.
"""
function parsePOSCAR(lines)

    title = lines[1]

    a = parse(Float64, split(lines[2])[1])

    A = zeros(3,3)
    A[1,:] = parseARRfloat(split(lines[3]))
    A[2,:] = parseARRfloat(split(lines[4]))
    A[3,:] = parseARRfloat(split(lines[5]))

    A = A * a / 0.529177

    thetypes = split(lines[6])

    thenumbers = parseARRint(split(lines[7]))

    if size(thetypes)[1] != size(thenumbers)[1]

        return -1,-1,-1
    end

    nat = sum(thenumbers)

    dc = split(lines[8])[1][1]
    cart = false
    if dc == 'c' || dc == 'C'
        cart = true
    end

    types = String[]
    for (t,n) in zip(thetypes,thenumbers)

        for i = 1:n
            push!(types, t)
        end
    end

    coords = zeros(nat,3)
    for i = 1:nat

        coords[i,:] = parseARRfloat(split(lines[8+i])[1:3])
    end

    if cart
        coords = (coords / 0.529177 )* inv(A)
    end

    return A, coords, types
end

"""
    function parseQEinput(lines)

Parse a quantum espresso inputfile. Can only handle simple cases with explict CELL_PARAMETERS
Cannot handle nonzero ibrav. or celldm

Called by makecrys, doesn't need to be called directly.
"""
function parseQEinput(lines)

#    println("parseQE")

    posind = -1
    Aind = -1
    types = String[]
    A=zeros(3,3)
    coords = Float64[]
    nat = -1
    units = 1.0
    for line in lines
        liner = replace(replace(line, "," => " "), "=" => " ")

        sp = split(liner)

        if size(sp)[1] >= 1

            if posind >= 0
                posind += 1
                coords[posind,:] = parseARRfloat(sp[2:4])
                push!(types, sp[1])
                if posind >= nat
                    posind = -1
                end
            end

            if Aind >= 0
                Aind += 1
                A[Aind,:] = parseARRfloat(sp[1:3])
                if Aind >= 3
                    Aind = -1
                end
            end

            if sp[1] == "nat"
                nat = parse(Int64, sp[2])
                coords = zeros(nat,3)

            elseif  sp[1] == "ATOMIC_POSITIONS"
                posind = 0
                if length(sp) > 1
                    if sp[2] != "crystal" && sp[2] != "(crystal)"
                        println("warning, only crystal coords supported !!!!!!!!!!!!!!!!!!")
                    end
                end
            elseif sp[1] == "CELL_PARAMETERS"
                Aind = 0
                if length(sp) > 1
                    if sp[2] == "angstrom" || sp[2] == "Angstrom" || sp[2] == "(angstrom)"
                        units = Ang
                    else
                        println("warning alat or other CELL_PARAMETERS not supported !!!!!!!!!!!!!!!!!!!!!!!")
                    end
                end
            end

        end
    end

    A = A / units

    return A, coords, types

end

"""
    function generate_supercell(crys, cell)

Generate supercell. cell is `[1,1,2]`, etc.

Note, perfered notation is to use syntax:
`c * [1,1,2]` , where `c` is a crystal, thus
using overloading of the `*` operator,
rather than calling directly.
"""
function generate_supercell(crys, cell)

    if typeof(cell) == Int
        cell = [cell,cell,cell]
    end

    A = copy(crys.A)
    A[1,:] *= cell[1]
    A[2,:] *= cell[2]
    A[3,:] *= cell[3]

    cells = prod(cell)

#    println("cells ", cell, " " , cells)

    coords = zeros(crys.nat*cells, 3)

    types=String[]

    c = 0
    for i in 0.0:cell[1]-1
        for j in 0.0:cell[2]-1
            for k in 0.0:cell[3]-1
                for at in 1:crys.nat
                    c+= 1

                    coords[c,:] = (crys.coords[at,:] + [i, j, k]) ./ cell'
                    push!(types, crys.types[at])
                end
            end
        end
    end

    csuper = makecrys(A, coords, types, units="Bohr")
    return csuper
end

"""
    function generate_random_distortion(crys, amag, strain_mag)

Randomly distort a crystal. `amag` is the atom distance, `strain_mag` is the strain magnitude
"""
function generate_random_distortion(crys, amag, strain_mag)

    amag = convert_length(amag)

    st = (rand(3,3) .- 0.5) * strain_mag
    st = (st + st')/2.0

    A = crys.A * (I + st)
    coords_real = (crys.coords * A) + (rand(crys.nat, 3) .- 0.5)*amag
    coords = coords_real * inv(A)

    return makecrys(A, coords, copy(crys.types), units="Bohr")

end

"""
    function write_poscar(crys, filename)

Write a `crystal` to a POSCAR.
"""
function write_poscar(crys, filename)

    fil = open(filename, "w")
    write(fil, "title kfg\n")
    write(fil, "1.0000000\n")

    Aang = crys.A * Ang

    write(fil, "$(Aang[1,1])  $(Aang[1,2])  $(Aang[1,3]) \n")
    write(fil, "$(Aang[2,1])  $(Aang[2,2])  $(Aang[2,3]) \n")
    write(fil, "$(Aang[3,1])  $(Aang[3,2])  $(Aang[3,3]) \n")

    tstr = ""
    numstr = ""

    ut = Set(crys.types)
    for t in ut
        tstr = tstr * t * " "
        count = sum( t .== crys.types)
        numstr = numstr * "$count "
    end

#    tstr = ""
#    numstr = ""
 #   for t in crys.types
 #       tstr = tstr * t * " "
 #       numstr = numstr * "1 "
 #   end

    write(fil,tstr*'\n')
    write(fil,numstr*'\n')

    write(fil,"Direct\n")
    for t in ut
        for at in 1:crys.nat
            if crys.types[at] == t
                write(fil, "$(crys.coords[at,1])  $(crys.coords[at,2]) $(crys.coords[at,3])\n")
            end
        end
    end
    close(fil)

    return 0

end

"""
    function write_xsf(c::crystal, filename="t.xsf", force=missing)

Write file for xcrysden. Forces optional
"""
function write_xsf(c::crystal; filename="t.xsf", force=missing)

    outstr = " INFO
nunit      1    1    1
unit   cell
celltype   primcell
shape   parapipedal
 END_INFO
 DIM-GROUP
           3           1
 PRIMVEC
"
    As = string.(c.A * Ang)
    for i = 1:3
        outstr *= As[i,1] * "   " *  As[i,2] * "   " * As[i,3]* "\n"
    end

    outstr *= " CONVVEC
"
    for i = 1:3
        outstr *= As[i,1] * "   " *  As[i,2] * "   " * As[i,3]* "\n"
    end

    outstr *= " PRIMCOORD
"
    nat = c.nat
    outstr *= "          $nat       1\n"

    pos = string.(c.coords * c.A * Ang)
    if !ismissing(force)
        fstring = string.(forces)
    else
        fstring = string.(zeros(size(c.coords)))
    end

    for i = 1:nat
#        println(atoms[c.types[i]].Z)
        outstr *= string(atoms[c.types[i]].Z)*"   "* pos[i,1] * "   " *  pos[i,2] * "   " * pos[i,3]*"   " * fstring[i,1]*"   " * fstring[i,2]*"   " * fstring[i,3]*"\n"
    end

    write(filename, outstr)

end

"""
    function write_axsf(CRYSTAL; filename="t.axsf", FORCES=missing)

write animated axsf file for xcrysden
Takes in array of crystal and optionally forces
"""
function write_axsf(CRYSTAL; filename="t.axsf", FORCES=missing)

    if typeof(CRYSTAL) == crystal
        write_xsf(CRYSTAL, filename=filename, force=FORCES)
    end

    nstep = length(CRYSTAL)
    outstr= "ANIMSTEPS  $nstep
CRYSTAL
"
    for i in 1:nstep
        c = CRYSTAL[i]
        nat = c.nat
        if !ismissing(FORCES)
            f = FORCES[i]
        else
            f = zeros(size(c.coords))
        end
        outstr *= "PRIMVEC $i\n"

        As = string.(c.A * Ang)
        for i = 1:3
            outstr *= As[i,1] * "   " *  As[i,2] * "   " * As[i,3]* "\n"
        end
        outstr *= "PRIMCOORD $i\n"
        outstr *= "$nat        1\n"
        pos = string.(c.coords * c.A * Ang)
        fstring = string.(f)

        for i = 1:nat
            outstr *= string(atoms[c.types[i]].Z)*"   "* pos[i,1] * "   " *  pos[i,2] * "   " * pos[i,3]*"   " * fstring[i,1]*"   " * fstring[i,2]*"   " * fstring[i,3]*"\n"
        end
    end

    write(filename, outstr)

end

"""
    function function write_efs(crys, energy, forces, stress, filename)

Write `crystal`, energy, force, stress to fake quantum espresso output file.
This is for testing purposes only.
"""
function write_efs(crys, energy, forces, stress, filename)

    fil = open(filename, "w")

    write(fil,"Program PWSCF fake\n")
    write(fil,"number of atoms/cell = $(crys.nat) \n")
#    ntypes = pos.shape[0]
    ntypes = 1
    write(fil,"number of types = 1\n")
    write(fil,"celldm(1)= 1.00\n")

    write(fil,"a(1) = ( $(crys.A[1,1]) $(crys.A[1,2]) $(crys.A[1,3]) ) \n")
    write(fil,"a(2) = ( $(crys.A[2,1]) $(crys.A[2,2]) $(crys.A[2,3]) ) \n")
    write(fil,"a(3) = ( $(crys.A[3,1]) $(crys.A[3,2]) $(crys.A[3,3]) ) \n")

    write(fil,"     site n.     atom                  positions (cryst. coord.)\n")
    for na in 1:crys.nat
        write(fil,"         1       $(crys.types[na])   tau(   $na) = (  $(crys.coords[na,1]) $(crys.coords[na,2]) $(crys.coords[na,3])  )\n")
    end
    write(fil,"\n")

    write(fil,"!    total energy              =     $energy Ry\n")

    write(fil,"     Forces acting on atoms (Ry/au):\n")
    for na in 1:crys.nat
        write(fil,"     atom    $na type  1   force =   $(forces[na,1]) $(forces[na,2]) $(forces[na,3]) \n")
    end

    write(fil,"The non-local\n")

    write(fil,"          total   stress  (Ry/bohr**3)                   (kbar)     P=  ???\n")
    write(fil,"$(stress[1,1]) \t $(stress[1,2]) \t  $(stress[1,3])   0 0 0 \n")
    write(fil,"$(stress[2,1]) \t $(stress[2,2]) \t  $(stress[2,3])   0 0 0 \n")
    write(fil,"$(stress[3,1]) \t $(stress[3,2]) \t  $(stress[3,3])   0 0 0 \n")

    write(fil,"JOB DONE.\n")

    close(fil)

end

"""
    function get_grid(c, kden=55.0)

Get a default k-point grid size with `kden` density.
"""
function get_grid(c, kden=55.0)

    B = transpose(inv(c.A))
#    kden = 55.0

#    b1 = 1.0/norm(c.A[1,:])
#    b2 = 1.0/norm(c.A[2,:])
#    b3 = 1.0/norm(c.A[3,:])

    b1 = norm(B[1,:])
    b2 = norm(B[2,:])
    b3 = norm(B[3,:])

    k1 = convert(Int, round(kden * b1))
    k2 = convert(Int, round(kden * b2))
    k3 = convert(Int, round(kden * b3))

    if k1%2 == 1
        k1 = k1 + 1
    end
    if k2%2 == 1
        k2 = k2 + 1
    end
    if k3%2 == 1
        k3 = k3 + 1
    end
    k1 = max(k1, 2)
    k2 = max(k2, 2)
    k3 = max(k3, 2)

    k1 = min(k1, 14)
    k2 = min(k2, 14)
    k3 = min(k3, 14)

#    k1 = min(k1, 30)
#    k2 = min(k2, 30)
#    k3 = min(k3, 30)

    kpoints = [k1, k2, k3]
#    println("get grid ", kpoints)
    return kpoints

#    R = [0,0,0]
#    for i = 1:3
#        R[i] = Int64(round(20.0/sum(c.A[1,:].^2)^0.5))
#        R[i] = max(R[i], 2)
#    end
#    println("get grid ", R)
#    return R

end

"""
    function orbital_index(c::crystal)

Get correspondence between `crystal` and the TB orbital numbers.

`return ind2orb, orb2ind, etotal, nval`

- `ind2orb` dictionary which gives `[atom_number,atom_type, :orbital_symbol]` from TB index integer.
- `orb2ind` dictionary which gives TB index integer from `[atom_number,atom_type, :orbital_symbol]`
- `etotal` total DFT energy of atoms, for calculation atomization energy.
- `nval` number of valence orbitals.
"""
function orbital_index(c::crystal)

    ind2orb = Dict()
    orb2ind = Dict()

    atomtypes = []
    ntypes=0
    for t in c.stypes
        if !(t in atomtypes)
            ntypes += 1
            push!(atomtypes, t)
        end
    end

    wan_counter = 1

    nsemi = 0
    nval = 0
    nwan = 0
    ind = 0

    etotal = 0.0

#    for at in atomtypes
#        for (i, t) in enumerate(c.types)
#            if t == at
    for (i,t) in enumerate(c.stypes)
        atom = atoms[t]
        nsemi += atom.nsemicore
        nval += atom.nval
        nwan += atom.nwan

        etotal += atom.total_energy

        for o in atom.orbitals
            if o == :s
                ind += 1
                ind2orb[ind] = [i,t, :s]
            elseif o == :p
                ind += 1
                ind2orb[ind] = [i,t, :pz]
                ind += 1
                ind2orb[ind] = [i,t, :px]
                ind += 1
                ind2orb[ind] = [i,t, :py]
            elseif o == :d
                ind += 1
                ind2orb[ind] = [i,t, :dz2]
                ind += 1
                ind2orb[ind] = [i,t, :dxz]
                ind += 1
                ind2orb[ind] = [i,t, :dyz]
                ind += 1
                ind2orb[ind] = [i,t, :dx2_y2]
                ind += 1
                ind2orb[ind] = [i,t, :dxy]
            elseif o == :f
                ind += 1
                ind2orb[ind] = [i,t, :fz3]
                ind += 1
                ind2orb[ind] = [i,t, :fxz2]
                ind += 1
                ind2orb[ind] = [i,t, :fyz2]
                ind += 1
                ind2orb[ind] = [i,t, :fz_x2y2]
                ind += 1
                ind2orb[ind] = [i,t, :fxyz]
                ind += 1
                ind2orb[ind] = [i,t, :fx_x2_3y2]
                ind += 1
                ind2orb[ind] = [i,t, :fy_3x2_y2]
            else
                error("orbitals ????")
            end

        end
        if ind != round(nwan/2)
            error("counting error", ind, " " , nwan)
        end

    end

    nwan = 0
    for (i,t) in enumerate(c.types)
        atom = atoms[t]
        orb2ind[i] = 1+nwan : nwan + Int(round(atom.nwan/2))
        nwan += Int(round(atom.nwan/2))
    end

    #reverse dictionary
#    for ind = 1:convert(Int, round(nwan/2))
#        orb2ind[ind2orb[ind]] = ind
#        println(ind, " ", ind2orb[ind][1], " ", ind2orb[ind][2])
#
#    end

    return ind2orb, orb2ind, etotal, nval

end

end #end module

#using .CrystalMod
#Base.:(==)(x::crystal, y::crystal) = (x.A == y.A && x.coords == y.coords && x.types == y.types && x.nat==y.nat)


text:
"""
`x, y = noisy_function(fun, x; noise = 0.01, f_rand = randn)`

Generates a noisy response `y` for the given function `fun`
by adding `noise .* f_randn(length(x))` to the result of `fun(x)`.
"""
function noisy_function(fun::Function, x::AbstractVector{T}; noise::Real = 0.01, f_rand::Function = randn) where T<:Real
    x_vec = collect(x)
    y = fun.(x_vec) .+ noise .* f_rand.()
    x_vec, y
end

"""
`x, y = noisy_sin(n, start, stop; noise = 0.3, f_rand = randn)`

Generates `n` noisy equally spaces samples of a sinus from `start` to `stop`
by adding `noise .* f_randn(length(x))` to the result of `fun(x)`.
"""
function noisy_sin(n::Int = 50, start::Real = 0, stop::Real = 2π; noise::Real = 0.3, f_rand::Function = randn)
    noisy_function(sin, range(start, stop=stop, length=n); noise = noise, f_rand = f_rand)
end

"""
`x, y = noisy_poly(coef, x; noise = 0.01, f_rand = randn)`

Generates a noisy response for a polynomial of degree `length(coef)`
using the vector `x` as input and adding `noise .* f_randn(length(x))` to the result.
The vector `coef` contains the coefficients for the terms of the polynome.
The first element of `coef` denotes the coefficient for the term with
the highest degree, while the last element of `coef` denotes the intercept.
"""
function noisy_poly(coef::AbstractVector{R}, x::AbstractVector{T}; noise::Real = 0.1, f_rand::Function = randn) where {T<:Real,R<:Real}
    n = length(x)
    m = length(coef)
    x_vec = collect(x)
    y = zeros(n)
    @inbounds for i = 1:n
        for k = 1:m
            y[i] += coef[k] * x_vec[i]^(m-k)
        end
    end
    y .+= noise .* f_rand.()
    x_vec, y
end

"""
`x, y = noisy_spiral(n, a, theta, b; noise = 0.01, f_rand = randn)`

Generates `n` noisy responses for a spiral with two labels. Uses the radius, angle
and scaling arguments to space the points in 2D space and adding `noise .* f_randn(n)`
to the response.
"""
function noisy_spiral(n::Int = 97, a::Real = 6.5, theta::Real = 16.0, b::Real=104.0; noise::Real = 0.1, f_rand::Function = randn)
    x = zeros(Float64, (2, 2*n))
    y = zeros(Int, 2*n)
    index = 0:1.0:(n-1)
    for i = 1:n
        _angle = index[i]*pi/theta
    	_radius = a * (b-index[i]) / b
    	x_coord = _radius * sin(_angle)
    	y_coord = _radius * cos(_angle)
    	x[1, i] = x_coord
    	x[2, i] = y_coord
        x[1, n+i] = -(x_coord)
    	x[2, n+i] = -(y_coord)
        y[i] = 1
        y[n+i] = 0
    end
    x[1, :] .+= noise .* f_rand.()
    x[2, :] .+= noise .* f_rand.()
    x, y
end


text:
using FFTW,Plots;N=128;P=64N;dt=1/4N;NT=1024;W=200;w=W/P*N;n=zeros(N);
v=collect(1:P.>P/2).*2 .-1.;u()=(x.=mod.(x.+v/2*dt,1));F=fft;IF=ifft;x=rand(P);
k=im.*2π*vcat(1,1:N/2,-N/2+1:-1);sc(a,b)=scatter(a,b,ylims=(-3,3),legend=0==1);
g(y)=(o=Int(mod1(ceil(y),N));([o,mod1(o-1,N)],[1-o+y,o-y]));O(z)=(z[1]*=0;z);
@gif for t in 1:NT;
  u();E=real.(IF((E=F((n.*=0;for i∈x;j=g(i);n[j[1]].+=j[2];end;n*w))./k;O(E))));
  u();i=0;for y∈x;i+=1;j=g(y);v[i]+=E[j[1]]'*j[2]*dt;end;sc(x,v);
end every 8


text:
@testset "SpeakerListenerEnv" begin
    env = SpeakerListenerEnv()

    RLBase.test_interfaces!(env)
    RLBase.test_runnable!(env)
end


text:
using LinearAlgebra
import LinearAlgebra.Eigen
using SparseArrays
#this lib will be used to estimate the executon time
using BenchmarkTools
using Test
using Distributions

struct SVDCM
	U::Matrix{Float64}
	S::SparseMatrixCSC{Float64, Int64}
	V::Matrix{Float64}

	function SVDCM(A::Matrix{Float64})
		U, S, V = svdCm(A)
		return new( U, S, V)
	end
	function svdCm(A::Matrix{Float64})

		r,c = size(A)
		shorter=min(r,c)
		longer=max(r,c)

		Asq = zeros(shorter,shorter)
		mul!(Asq, A', A)
		a, V = eigen!(Asq)

		S = sparse(zeros(longer,shorter))
		Sin = sparse(zeros(shorter,longer))

		S[diagind(S)] .= (sqrt.(a))
		Sin[diagind(Sin)] .= 1 ./ S[diagind(S)]

		U1 = zeros(r,c)
		U = zeros(longer,longer)

		mul!(U1, A, V)
		mul!(U,U1,Sin)
		return U, S, V
	end
end

R = rand(Uniform(1., 100.), 10000,1000)
MM = [1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0; 1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0; 1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0; 1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0;1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0;1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0 ]

C = R
println("size of matrix: ", size(C) )

println("Julia SVD time:")
SVD = @time svd(C)

println("SVD time:")
mySvd = @time SVDCM(C)

F = mySvd.U * mySvd.S * mySvd.V'
#display(F)
@test F ≈ C atol=1e-5