text:
#------------------------
#--------  Functions

#--------  random generation
# using DoubleDouble

function doubledot(x::Vector{Float64}, y::Vector{Float64})
    dsum=Double(0.0)
    for k=1:length(x)
        dsum+=Double(x[k])*Double(y[k])
    end
    dsum
end

function GenHalfArrow(n::Integer, p::Int)
    # generates n x n half arrowhad matrix with/out arrow top
    # This is one without arrowhead point
    p==0 ?
    HalfArrow(rand(n-1).-0.5,rand(n-1).-0.5) :
    # This is one with the point
    HalfArrow(rand(n-1).-0.5,rand(n).-0.5)
end

#-------- Inverses

function inv(A::HalfArrow{T},i::Integer,τ::Vector{Float64}=[1e3, 10*length(A.D)]) where T
    # COMPUTES: inverse of a SymArrow matrix B, where B=A'*A, and A =[A.D, A.z]
    # is the HalfArrow. Here inv(B-B.D[i]*I) is again SymArrow.
    # Uses higher precision to compute top of the arrow element accurately, if
    # needed.
    # τ=[tolb,tolz] are tolerances, usually [1e3, 10*n]
    # [0.0,0.0] forces DoubleDouble, [1e50,1e50] would never use it
    # RETURNS:  SymArrow(D,z,b,i), Kb, Kz, Qout
    # Kb - condition Kb, Kz - condition Kz, Qout = 1 / 0 - double was / was not used

    n=length(A.D)
    D=Array{T}(undef,n)
    z=Array{T}(undef,n)
    z1=Array{T}(undef,n)
    z1.=A.z[1:n].*A.D
    wz=one(T)/(z1[i])
    σ=A.D[i]
    for k=1:i-1
        D[k]=one(T)/((A.D[k]-σ)*(A.D[k]+σ))
        z[k]=-z1[k]*D[k]*wz
    end
    for k=i+1:n
        D[k-1]=one(T)/((A.D[k]-σ)*(A.D[k]+σ))
        z[k-1]=-z1[k]*D[k-1]*wz
    end

    # D=[1./(A.D[1:i-1]-σ),1./(A.D[i+1:end]-σ)]
    # wz=1/A.z[i]
    # z=-[A.z[1:i-1], A.z[i+1:end]].*D1*wz1
    # Maybe this needs more accuracy??
    # if A.z[k]^2>A.D[i]^2 for some k, then we should use (A.z[k]+A.D[i])*(A.z[k]-A.d[i]) +
    # the rest, which is sum of squares
    # a=dot(A.z,A.z)-A.D[i]^2 # this was original formula
    ad=doubledot([A.z;A.D[i]],[A.z;-A.D[i]])
    a=ad.hi+ad.lo
    # compute the sum in a plain loop
    P=zero(T)
    Q=zero(T)
    Qout=0
    nn=n-1
    for k=1:i-1
        D[k]>0.0 ? P+=z1[k]^2*D[k] : Q+=z1[k]^2*D[k]
    end
    for k=i:nn
        D[k]>0.0 ? P+=z1[k+1]^2*D[k] : Q+=z1[k+1]^2*D[k]
    end
    a<0 ? P-=a : Q-=a

    Kb=(P-Q)/abs(P+Q)
    Kz=maximum(abs,z1)*abs(wz)
    if Kb<τ[1] ||  Kz<τ[2]
        b=(P+Q)*wz*wz
    else  # recompute in Double
        Qout=1
        σ₁=map(Double,A.D[i])
        Dd=[Double{Float64}[(Double(A.D[k])+σ₁)*(Double(A.D[k])-σ₁) for k=1:i-1];
            Double{Float64}[(Double(A.D[k])+σ₁)*(Double(A.D[k])-σ₁) for
                            k=i+1:length(A.D)]]
        wzd=Double(A.z[i])*σ₁
        # ad=doubledot(A.z,A.z)-σ₁^2 # we already have this if
        # we use doubledot above
        Pd,Qd=map(Double,(0.0,0.0))

        for k=1:i-1
            Dd[k].hi>0.0 ? Pd+=(Double(A.z[k])*Double(A.D[k]))^2/Dd[k] :
            Qd+=(Double(A.z[k])*Double(A.D[k]))^2/Dd[k]
        end
        for k=i:nn
            Dd[k].hi>0.0 ? Pd+=(Double(A.z[k+1])*Double(A.D[k+1]))^2/Dd[k] :
            Qd+=(Double(A.z[k+1])*Double(A.D[k+1]))^2/Dd[k]
            # @show P,Q
        end
        ad.hi<0 ?   Pd-=ad : Qd-=ad
        bd=(Pd+Qd)/(wzd*wzd)
        b=bd.hi+bd.lo
    end
    # SymArrow([D[1:A.i-2],zero(T),D[A.i-1:end]],[z[1:A.i-2],wz,z[A.i-1:end]],b,i),Kb,Kz,Qout
    # This reduces memory allocation (4 vectors) and is 30-40% faster
    D[n]=zero(T)
    z[n]=wz

    # return this
    SymArrow(D,z,b,i),Kb,Kz,Qout
end # inv

function inv!(B::SymArrow{T},A::HalfArrow{T},i::Integer,τ::Vector{Float64}=[1e3, 10*length(A.D)]) where T
    # In-place version!
    # COMPUTES: inverse of a SymArrow matrix B, where B=A'*A, and A =[A.D, A.z]
    # is the HalfArrow. Here inv(B-B.D[i]*I) is again SymArrow.
    # Uses higher precision to compute top of the arrow element accurately, if
    # needed.
    # τ=[tolb,tolz] are tolerances, usually [1e3, 10*n]
    # [0.0,0.0] forces DoubleDouble, [1e50,1e50] would never use it
    # RETURNS:  SymArrow(D,z,b,i), Kb, Kz, Qout
    # Kb - condition Kb, Kz - condition Kz, Qout = 1 / 0 - double was / was not used
    n=length(A.D)
    # D=Array{T}(undef,n)
    # z=Array{T}(undef,n)
    z1=Array{T}(undef,n)
    z1.=A.z[1:n].*A.D
    wz=one(T)/(z1[i])
    σ=A.D[i]
    for k=1:i-1
        B.D[k]=one(T)/((A.D[k]-σ)*(A.D[k]+σ))
        B.z[k]=-z1[k]*B.D[k]*wz
    end
    for k=i+1:n
        B.D[k-1]=one(T)/((A.D[k]-σ)*(A.D[k]+σ))
        B.z[k-1]=-z1[k]*B.D[k-1]*wz
    end

    ad=doubledot([A.z;A.D[i]],[A.z;-A.D[i]])
    a=ad.hi+ad.lo
    # compute the sum in a plain loop
    P=zero(T)
    Q=zero(T)
    Qout=0
    nn=n-1
    for k=1:i-1
        B.D[k]>0.0 ? P+=z1[k]^2*B.D[k] : Q+=z1[k]^2*B.D[k]
    end
    for k=i:nn
        B.D[k]>0.0 ? P+=z1[k+1]^2*B.D[k] : Q+=z1[k+1]^2*B.D[k]
    end
    a<0 ? P-=a : Q-=a

    Kb=(P-Q)/abs(P+Q)
    Kz=maximum(abs,z1)*abs(wz)
    if Kb<τ[1] ||  Kz<τ[2]
        b=(P+Q)*wz*wz
    else  # recompute in Double
        Qout=1
        σ₁=map(Double,A.D[i])
        Dd=[Double{Float64}[(Double(A.D[k])+σ₁)*(Double(A.D[k])-σ₁) for k=1:i-1];
            Double{Float64}[(Double(A.D[k])+σ₁)*(Double(A.D[k])-σ₁) for
                            k=i+1:length(A.D)]]
        wzd=Double(A.z[i])*σ₁
        # ad=doubledot(A.z,A.z)-σ₁^2 # we already have this if
        # we use doubledot above
        Pd,Qd=map(Double,(0.0,0.0))
        for k=1:i-1
            Dd[k].hi>0.0 ? Pd+=(Double(A.z[k])*Double(A.D[k]))^2/Dd[k] :
            Qd+=(Double(A.z[k])*Double(A.D[k]))^2/Dd[k]
        end
        for k=i:nn
            Dd[k].hi>0.0 ? Pd+=(Double(A.z[k+1])*Double(A.D[k+1]))^2/Dd[k] :
            Qd+=(Double(A.z[k+1])*Double(A.D[k+1]))^2/Dd[k]
            # @show P,Q
        end
        ad.hi<0 ?   Pd-=ad : Qd-=ad
        bd=(Pd+Qd)/(wzd*wzd)
        b=bd.hi+bd.lo
    end
    # SymArrow([D[1:A.i-2],zero(T),D[A.i-1:end]],[z[1:A.i-2],wz,z[A.i-1:end]],b,i),Kb,Kz,Qout
    # This reduces memory allocation (4 vectors) and is 30-40% faster
    B.D[n]=zero(T)
    B.z[n]=wz
    B.a=b
    B.i=i
    # return this
    # SymArrow(D,z,b,i),Kb,Kz,Qout
    return Kb, Kz, Qout
end # inv!

function inv(A::HalfArrow{T}, σ::Float64, τ::Float64=1.0e3) where T
    # COMPUTES: inverse of a SymArrow matrix B, where B=A'*A, and A =[A.D, A.z]
    # is the HalfArrow. Here inv(B-σ^2*I) is a SymDPR1.
    # uses DoubleDouble to compute ρ accurately, if needed.
    # τ is tolerance, usually 1e3,  0.0 forces Double, 1e50 would never use it
    # RETURNS: SymDPR1(D,u,ρ), Kρ, Qout
    # Kρ - condition Kρ, Qout = 1 / 0 - Double was / was not used

    n=length(A.D)
    D=Array{T}(undef,n+1)
    u=Array{T}(undef,n+1)
    z1=Array{T}(undef,n)
    z1.=A.z[1:n].*A.D
    for k=1:n
        D[k]=one(T)/((A.D[k]-σ)*(A.D[k]+σ))
        u[k]=z1[k]*D[k]
    end
    D[n+1]=zero(T)
    u[n+1]=-one(T)

    # compute ρ and Kρ
    #--- compute the sum in a plain loop
    P=zero(T)
    Q=zero(T)
    Qout=0
    # a=dot(A.z,A.z)-σ^2 # try with more accuracy
    ad=doubledot([A.z;σ],[A.z;-σ])
    a=ad.hi+ad.lo
    for k=1:n
        D[k]>0.0 ? P+=z1[k]^2*D[k] : Q+=z1[k]^2*D[k]
    end
    a<0 ? P-=a : Q-=a

    # Condition of ρ
    Kρ=(P-Q)/abs(P+Q)

    if Kρ<τ
        ρ=-1.0/(P+Q)
    else  # recompute in Double
        Qout=1
        Pd,Qd=map(Double,(0.0,0.0))
        σ₁=Double(σ)
        Dd=Double{Float64}[(Double(A.D[k])+σ₁)*(Double(A.D[k])-σ₁) for k=1:n]
        for k=1:n
            Dd[k].hi>0.0 ? Pd+=(Double(A.z[k])*Double(A.D[k]))^2/Dd[k] :
            Qd+=(Double(A.z[k])*Double(A.D[k]))^2/Dd[k]
        end
        ad.hi+ad.lo<0 ? Pd-=ad : Qd-=ad
        r=Double(1.0)/(Pd+Qd)
        ρ=-(r.hi+r.lo)
    end

    # returns the following
    SymDPR1(D,u,ρ), Kρ, Qout
end # inv

function inv(A::HalfArrow{T}, σ::Double) where T
    # COMPUTES: inverse of a SymArrow matrix B, where B=A'*A, and A =[A.D, A.z]
    # is the HalfArrow. Here inv(B-σ^2*I) is SymDPR1.
    # Here shift is Double so it uses Double to compute everything
    # RETURNS: SymDPR1(D1,u1,ρ), Qout
    # Qout = 1 on exit meaning Double was used

    n=length(A.D)
    D=Array{Double}(undef,n+1)
    u=Array{Double}(undef,n+1)
    oned=Double(1.0)
    zerod=Double(0.0)
    for k=1:n
        D[k]=(Double(A.D[k])-σ)*(Double(A.D[k])+σ)
        u[k]=Double(A.z[k])*Double(A.D[k])/D[k]
    end
    a=doubledot(A.z,A.z)-σ^2
    D[n+1]=zerod
    u[n+1]=Double(-1.0)

    # compute ρ and Kρ
    #--- compute the sum in a plain loop
    P,Q=zerod,zerod
    Qout=1
    for k=1:n
        D[k].hi>0.0 ? P+=(Double(A.z[k])*Double(A.D[k]))^2/D[k] :
        Q+=(Double(A.z[k])*Double(A.D[k]))^2/D[k]
        D[k]=oned/D[k]
    end
    a.hi+a.lo<0.0  ? P-=a : Q-=a
    r=oned/(P+Q)
    ρ=-(r.hi+r.lo)

    # returns the following
    D₁=Array{T}(undef,n+1)
    u₁=Array{T}(undef,n+1)
    for k=1:n+1
        D₁[k]=D[k].hi+D[k].lo
    end
    for k=1:n+1
        u₁[k]=u[k].hi+u[k].lo
    end
    SymDPR1(D₁,u₁,ρ), Qout
end # inv

function  svd(A::HalfArrow{T},Ainv::SymArrow{T},k::Integer,
    τ::Vector{Float64}=[1e3,10.0*length(A.D),1e3,1e3,1e3]) where T
    # COMPUTES: k-th singular value triple of an ordered irreducible HalfArrow
    # A = [Diagonal(A.D) A.z] with A.D > 0
    # τ=[tolb,tolz,tolnu,tolrho,tollambda] = [1e3,10.0*n,1e3,1e3,1e3]
    # RETURNS: λ, u, v, Info
    # where
    # λ - k-th singular value in descending order
    # u - λ's normalized left singular vector
    # v - λ's normalized right singular vector
    # Info = Sind, Kb, Kz, Kν, Kρ, Qout
    # Sind = shift index i for the k-th singular value
    # Kb, Kz, Kν, Kρ - condition numbers
    # Qout = 1 / 0 - Double was / was not used

    # Set the dimension
    n = length(A.D) + 1    # Set all conditions initially to zero
    Kb,Kz,Kν,Kρ=0.0,0.0,0.0,0.0
    Qout=0
    u=zeros(length(A.z))
    v=zeros(n)
    # z1=Array{Float64}(undef,n)
    z1=A.z[1:n-1].*A.D

    # Kz is former kappa_nu
    # Determine the shift σ, the shift index i, and whether λ
    # is on the left or the right side of the nearest pole
    # Exterior eigenvalues (k = 1 or k = n):

    if k == 1
        σ,i,side = A.D[1],1,'R'
    elseif k==n
        σ,i,side = A.D[n-1],n-1,'L'
    else
        # Interior eigenvalues (k in (2,...n-1) ):
        Dtemp = (A.D.-A.D[k]).*(A.D.+A.D[k])
        atemp = dot(A.z,A.z)-A.D[k]^2
        middle = Dtemp[k-1]/2.0
        Fmiddle = (atemp-middle)-sum(z1.^2 ./(Dtemp.-middle))
        σ,i,side = Fmiddle < 0.0 ? (A.D[k],k,'R') : (A.D[k-1],k-1,'L')
    end

    # Compute the inverse of the shifted matrix, (A'*A)_i^(-1), Kb and Kz
    # Ainv,Kb,Kz,Qout = inv(A,i)
    Kb,Kz,Qout = inv!(Ainv,A,i)

    # Compute the eigenvalue of the inverse shifted matrix
    ν = bisect( Ainv,side )

    #  ν=fastzero([invD1; 0; invD2], [w1;wz;w2], b, side); # implement later
    #  [ν-ν₁]/ν, [ν-nueva]/ν, pause  # and compare
    if abs(ν)==Inf
        # this is nonstandard
        # Deflation in aheig (ν=Inf)
        u[i]=1.0
        v[i]=1.0
        λ=σ
    else
        # standard case, full computation
        # ν₁ is the F- or 1-norm of the inverse of the shifted matrix
        # ν₁0=maximum([sum(abs,Ainv.z)+abs(Ainv.a), maximum(abs,Ainv.D)+abs.(Ainv.z))])
        ν₁=0.0
        for k=1:n-1
            ν₁=max(ν₁,abs(Ainv.D[k])+abs(Ainv.z[k]))
        end
        ν₁=max(sum(abs,Ainv.z)+abs(Ainv.a), ν₁)
        Kν=ν₁/abs(ν)
        if Kν<τ[3]
            # Accuracy is fine, compute the eigenvector
            μ = 1.0/ν
            #         v=[ z1./((A.D-σ).*(A.D+σ)-μ);-1.0]
            for k=1:n-1
	        v[k] = z1[k]/((A.D[k]-σ)*(A.D[k]+σ)-μ)
            end
            v[n]=-1.0
            normalize!(v)
            λ=sqrt(μ+σ^2) # this may have errors
            u[1:n-1]=λ*v[1:n-1]./A.D
            if length(A.z)==n
                u[n]=A.z[n]*v[n]/λ
            end
        else
            # Remedies according to Remark 3 - we shift between original
            # eigenvalues and compute DPR1 matrix
            # 1/ν₁+σ, 1/ν+σ
            ν = side=='R' ? abs(ν) : -abs(ν)
            ν₁=-sign(ν)*ν₁
            σ₁=(ν₁+ν)/(2.0*ν*ν₁)+σ^2

            if count(isequal(sqrt(σ₁)),abs.(A.D))>0 # we came back with a pole
                # recompute σᵥ more accurately according with dekker
                σᵥ=(Double(ν₁)+Double(ν))/(Double(2.0)*Double(ν)*Double(ν₁))+Double(σ)^2
                # Compute the inverse of the shifted arrowhead (DPR1)
                Ainv,Qout₁=inv(A,sqrt(σᵥ)) # Ainv is Float64, here it was sqrt(σᵥ)
                ν₁=bisect(Ainv,side)
                μ₁ = 1.0/ν₁
                D0=map(A.D,Double)
                D0=A.D.*A.D-σᵥ
                D1=D0.hi+D0.lo
                # Shift the eigenvalue back in Double
                lam = Double(1.0)/Double(ν₁)+σᵥ
                λ=sqrt(lam.hi+lam.lo)

                if count(isequal(μ₁),D1)>0
                    ind=findall(isequal(μ₁),D1)
                    u[ind]=1.0
                    v[ind]=1.0
                else
                    for k=1:n-1
	                v[k] = z1[k]/(D1[k]-μ₁)
                    end
                    v[n]=-1.0
                    normalize!(v)
                    u[1:n-1]=λ*v[1:n-1]./A.D
                    if length(A.z)==n
                        u[n]=A.z[n]*v[n]/λ
                    end
                    # v=[ A.z./(D1-μ₁);-1.0]
                end
            else
                # Compute the inverse of the shifted arrowhead (DPR1)
                σₛ=sqrt(σ₁)
                Ainv, Kρ,Qout₁=inv(A,σₛ,τ[4]) # Ainv is Float64
                # Compute the eigenvalue by bisect for DPR1
	        # Note: instead of bisection could use dlaed4 (need a wrapper) but
	        # it is not faster. There norm(u)==1
                ν₁= bisect(Ainv,side)
                μ₁=1.0/ν₁
                for k=1:n-1
	            v[k] = z1[k]/((A.D[k]-σₛ)*(A.D[k]+σₛ)-μ₁)
                end
                v[n]=-1.0
                normalize!(v)
                λ=sqrt(μ₁+σ₁)
                u[1:n-1]=λ*v[1:n-1]./A.D
                if length(A.z)==n
                    u[n]=A.z[n]*v[n]/λ
                end
            end
            Qout==1 && (Qout=Qout+2)
        end

        # Remedy according to Remark 1 - we recompute the the eigenvalue
        # near zero from the inverse of the original matrix (a DPR1 matrix).
        # Here this can only happen for k==n
        if k==n && (abs(A.D[i])^2+abs(1.0/ν))/abs(λ^2)>τ[5]
            # if k==1 && A.D[1]<0.0 || k==n && A.D[n-1]>0.0 || side=='L' && sign(A.D[i])+sign(A.D[i+1])==0 || side=='R' && sign(A.D[i])+sign(A.D[i-1])==0
            # Compute the inverse of the original arrowhead (DPR1)
            Ainv,Kρ,Qout₁ = inv(A,0.0,τ[4]) # Ainv is Float64
            Qout==1 && (Qout=Qout+4)
            if abs(Ainv.r)==Inf
                λ=0.0
            else
                # Here we do not need bisection. We compute the Rayleigh
                # quotient by using already computed vectors which are
                # componentwise accurate
                ν₁=sum(v.^2 .*Ainv.D)+Ainv.r*sum(v.*Ainv.u)^2;
                λ=sqrt(1.0/ν₁)
            end
        end
    end
    # Return this
    κ=Info(i,Kb,Kz,Kν,Kρ,Qout)
    return λ,u,v,κ
end # svd (k)

function svd(A::HalfArrow{T}, τ::Vector{Float64}=[1e3,10.0*length(A.D),1e3,1e3,1e3]) where T
    # COMPUTES: all singular values and singular vectors of a real HalfArrow
    # A = [diagm(A.D) A.z]
    # τ=[tolb,tolz,tolnu,tolrho,tollambda] = [1e3,10.0*n,1e3,1e3,1e3]
    # RETURNS: SVD(U, Σ, V), κ
    # where
    # U = left singular vectors
    # Σ = singular values in decreasing order
    # V = right singular vectors
    # κ[k]=Info(Sind[k], Kb[k], Kz[k], Kν[k], Kρ[k], Qout[k])
    # Sind[k] - shift index i for the k-th eigenvalue
    # Kb, Kz, Kν, Kρ [k] - respective conditions for the k-th singular value
    # Qout[k] = 1 / 0 - Double was / was not used when computing k-th singular value

    nd=length(A.D)
    n=length(A.z)
    n0=n
    # Ordering the matrix
    signD=sign.(A.D)
    D=abs.(A.D)
    is=sortperm(D,rev=true)
    D=D[is]
    signD=signD[is]
    z=A.z[1:n]
    z=z[is]
    if n==nd
        U=Array{T}(I,n,n)
        V=Array{T}(I,n+1,n)
    else
        U=Array{T}(I,n,n)
        V=Array{T}(I,n,n)
    end
    Σ=zeros(n)
    κ=[Info(0, 0.0, 0.0, 0.0, 0.0, 0) for i=1:n]

    #  test for deflation in z
    z0=findall(iszero,z)
    zx=findall(!iszero,z)
    if isempty(zx)  # nothing to do
        if n==nd
            Σ=[A.D]
            isΣ=sortperm(Σ,rev=true)
            Σ=Σ[isΣ]
            V[:,1:n]=view(V,:,isΣ) # V[:,isΣ]
            U=view(U,invperm(isΣ),:) # U[invperm(isΣ),:]
        else
            Σ=[A.D;z[n]]
            isΣ=sortperm(Σ,rev=true)
            Σ=Σ[isΣ]
            V=view(V,:,isΣ) # V[:,isΣ]
            U=view(U,invperm(isΣ),:) # U[invperm(isΣ),:]
        end
        return SVD(U,Σ,V), κ
    end

    if !isempty(z0)
        Σ[z0]=D[z0]
        for k=1:length(z0)
            V[z0[k],z0[k]]=signD[z0[k]]
        end
        D=D[zx]
        # signD=signD[zx]
        z=z[zx]
        if isempty(z)
            if n>nd
                Σ[n]=abs(A.z[n])
                V[n,n]=sign(A.z[n])
            end
            # return U,Σ,V,Sind,Kb,Kz,Kν,Kρ,Qout
        end
    else
        if n0==nd
            zxv=[zx;n0+1]
        else
            zx=zxv=[zx;n0]
            z=[z;A.z[n0]]
        end
        n=length(z)

        #  Test for deflation in D
        g=D[1:n-2]-D[2:n-1]
        # Can play with inexact deflation
        # g0=find(abs(g)<eps)
        # gx=find(abs(g)>=eps)
        # Only exact deflation !!
        g0=findall(iszero,g)
        gx=findall(!iszero,g)
        if !isempty(g0)
            # Deflation
            Dgx=D[gx]; zgx=z[gx]
            lg0=length(g0)
            R=Array{Tuple{LinearAlgebra.Givens{Float64},Float64}}(undef,lg0)
            for l=lg0:-1:1
                R[l]=givens(z[g0[l]],z[g0[l]+1],zx[g0[l]],zx[g0[l]+1])
                z[g0[l]]=R[l][2]; z[g0[l]+1]=0.0
                # A_mul_Bc!(U,R) # multiply R'*U later
                Σ[zx[g0[l]+1]]=D[g0[l]+1]
            end
            # remains
            gx=[0;gx].+1
            nn=length(gx)
            zxx=zx[[gx;n]]
            if n0==nd
                zxxv=[zxx,n0+1]
            else
                zxx=zxxv=[zxx,n0]
                # z=[z;A.z[n0]]
            end
            Axx=HalfArrow(D[gx],z[gx])
            Bxx=[SymArrow(Vector{T}(undef,nn),Vector{T}(undef,nn),one(T),nn+1) for i=1:Threads.nthreads()]
            # Threads.@threads
            for k=1:nn+1
                tid=Threads.threadid()
                Σ[zxx[k]],U[zxx,zxx[k]],V[zxxv,zxx[k]],κ[zxx[k]]=svd(Axx,Bxx[tid],k)
            end
            for l=1:lg0
                # U=R[l][1]'*U
                lmul!(R[l][1],U)
            end
        else
            # No deflation in D
            Ax=HalfArrow(D,z)
            Bx=[SymArrow(Vector{T}(undef,n),Vector{T}(undef,n),one(T),n+1) for i=1:Threads.nthreads()]
            # Threads.@threads
            for k=1:n
                tid=Threads.threadid()
                Σ[zx[k]],U[zx,zx[k]],V[zxv,zx[k]],κ[zx[k]]=svd(Ax,Bx[tid],k)
                V[zxv[1:end-1],zx[k]]=V[zxv[1:end-1],zx[k]].*signD[zxv[1:end-1]]
            end
        end
    end
    # back premutation of vectors
    isi=sortperm(is)
    # must sort Σ once more
    es=sortperm(Σ,rev=true)
    Σ.=Σ[es]
    if n0==nd
        U=view(U,isi,es)
        V=view(V,[isi;n0+1],es)
    else
        isiv=[isi;n0]
        U=view(U,isiv,es)
        V=view(V,isiv,es)
    end
    κ.=κ[es]
    # Return this
    return SVD(U,Σ,V),κ
end # svd (all)


text:
#
# Validation utility functions
#

function get_data(file_name::String)
    text = collect.(readlines(file_name));
    alphabet = [union(unique.(text)...)..., '\n'];
    N = length(alphabet)

    return text, alphabet, N
end

function descr_data(dat::Array)
    alphabet = [unique(union.(collect.(dat))...)...];
    stop = Char(maximum(float.(alphabet))+1)
    push!(alphabet, stop)
    N = length(alphabet)

    return alphabet, stop, N
end

ext: jl


text:
# Numerical_Integration.jl
# Date: 28.06.2019. Created by Joshua Simon.

# Define some functions you want to integrate

function com_simpson(f, a, b, n)
    # Composite Simpson's rule
    # Integration of function f over interval [a,b] where n is
    # the number of subintervals. Accuracy only depends on the
    # number of subintervals.

    sum = 0

    if(a > b)
        println("ERROR: Lower interval bound a has to be smaller than upper bound b.")
        return 0
    end

    # Step length h
    h = (b - a) / n

    # Calculate integral value with composite simpson's rule
    for k = 1:n

        x_k = a + k * h
		x_k1 = a + (k - 1) * h

		simpson = h/6 * ( f(x_k1) + f(x_k) + 4* f((x_k1 + x_k)/2) )

        sum += simpson
    end

    return sum

end

integral_k_1 = com_simpson(sin, 0.0, pi*2.0, 100000)

println(integral_k_1)


text:
using BoltzmannMachines
import BoltzmannMachinesPlots # only needed for plotting capabilities

# ==============================================================================
# Binary data: BernoulliRBM and BasicDBM
# ------------------------------------------------------------------------------

# Generate data and split it in training and test data set
nsamples = 100;
nvariables = 36;
using Random
Random.seed!(1);
x = barsandstripes(nsamples, nvariables);
x, xtest = splitdata(x, 0.3);
datadict = DataDict("Training data" => x, "Test data" => xtest);

# Fit a BernoulliRBM and collect monitoring data during training.
# For small models, the model's loglikelihood can be computed exactly.
monitor = Monitor();
rbm = fitrbm(x; nhidden = 12,
      epochs = 80, learningrate = 0.006,
      monitoring = (rbm, epoch) -> begin
            monitorexactloglikelihood!(monitor, rbm, epoch, datadict)
            monitorreconstructionerror!(monitor, rbm, epoch, datadict)
      end);

# Plot the collected data (requires Gadfly)
BoltzmannMachinesPlots.plotevaluation(monitor, monitorexactloglikelihood)
BoltzmannMachinesPlots.plotevaluation(monitor, monitorreconstructionerror)

Random.seed!(12);

# DBM-Fitting: Pretraining and Fine-Tuning combined in one function
dbm = fitdbm(x, nhiddens = [6;2], epochs = 20, learningrate = 0.05);

# .. with extensive monitoring
monitor = Monitor(); monitor1 = Monitor(); monitor2 = Monitor();
dbm = fitdbm(x, epochs = 20, learningrate = 0.05,
      monitoringdatapretraining = datadict,
      pretraining = [
            TrainLayer(nhidden = 6, monitoring = (rbm, epoch, datadict) ->
                  monitorreconstructionerror!(monitor1, rbm, epoch, datadict));
            TrainLayer(nhidden = 2, monitoring = (rbm, epoch, datadict) ->
                  monitorreconstructionerror!(monitor2, rbm, epoch, datadict))
            ],
      monitoring = (dbm, epoch) ->
            monitorexactloglikelihood!(monitor, dbm, epoch, datadict));

# Monitoring plots
BoltzmannMachinesPlots.plotevaluation(monitor1, monitorreconstructionerror)
BoltzmannMachinesPlots.plotevaluation(monitor2, monitorreconstructionerror)
BoltzmannMachinesPlots.plotevaluation(monitor, monitorexactloglikelihood)

# Evaluate final result with exact computation of likelihood
exactloglikelihood(dbm, xtest)

# If the model has more parameters, in this case more hidden nodes,
# the exact calclation is not feasible any more:
# We need to calculate the likelihood using AIS.

Random.seed!(12);
monitor = Monitor();
rbm = fitrbm(x, nhidden = 36, epochs = 100,
      learningrates = [0.008*ones(80); 0.001*ones(20)],
      monitoring = (rbm, epoch) ->
            monitorloglikelihood!(monitor, rbm, epoch, datadict));

BoltzmannMachinesPlots.plotevaluation(monitor, monitorloglikelihood; sdrange = 3.0)

# In the DBM, the estimation of the log likelihood is much slower than in the
# RBM because it needs a AIS run for each sample.
# The lower bound of the log probability is faster to estimate:
# (But on the other hand, in some cases it might be too conservative to be
# useful as it is only a lower bound.)

Random.seed!(2);
monitor = Monitor();
dbm = stackrbms(x; nhiddens = [36;10;5], predbm = true, learningrate = 0.05);
traindbm!(dbm, x; epochs = 50, learningrate = 0.008,
      monitoring = (rbm, epoch) ->
            monitorlogproblowerbound!(monitor, rbm, epoch, datadict));

BoltzmannMachinesPlots.plotevaluation(monitor, monitorlogproblowerbound; sdrange = 3.0)

# Evaluate final result with AIS-estimated likelihood
loglikelihood(dbm, xtest)

# ==============================================================================
# Real valued data: GaussianBernoulliRBM
# ------------------------------------------------------------------------------

# Use "iris" dataset as example data to train a GaussianBernoulliRBM
using DelimitedFiles
x = convert(Matrix{Float64}, readdlm(
            joinpath(dirname(pathof(BoltzmannMachines)), "..", "test/data/iris.csv"), ',',
            header = true)[1][:,1:4]);
Random.seed!(12);
x, xtest = splitdata(x, 0.3);
datadict = DataDict("Training data" => x, "Test data" => xtest);
monitor = Monitor();
rbm = fitrbm(x, rbmtype = GaussianBernoulliRBM,
      nhidden = 3, epochs = 70, learningrate = 0.001,
      monitoring = (rbm, epoch) ->
            monitorexactloglikelihood!(monitor, rbm, epoch, datadict));

BoltzmannMachinesPlots.plotevaluation(monitor, monitorexactloglikelihood)

# ==============================================================================
# Data with binomially distributed values in {0,1,2}: Binomial2BernoulliRBM
# ------------------------------------------------------------------------------

# Generate data and split it in training and test data set
nsamples = 100;
nvariables = 25;
Random.seed!(12);
x = barsandstripes(nsamples, nvariables) + barsandstripes(nsamples, nvariables);
x, xtest = splitdata(x, 0.3);
datadict = DataDict("Training data" => x, "Test data" => xtest);

# Fit Binomial2BernoulliRBM.
monitor = Monitor();
rbm = fitrbm(x, rbmtype = Binomial2BernoulliRBM,
      nhidden = 10, epochs = 100, learningrates = [0.003*ones(60); 0.001*ones(40)],
      monitoring = (rbm, epoch) ->
            monitorexactloglikelihood!(monitor, rbm, epoch, datadict));

BoltzmannMachinesPlots.plotevaluation(monitor, monitorexactloglikelihood)

# ==============================================================================
# Examples for cross-validation
# ------------------------------------------------------------------------------

nsamples = 500;
nvariables = 16;
x = barsandstripes(nsamples, nvariables);

# Determine the optimal number of training epochs for a RBM
monitor = crossvalidation(x,
      (x, datadict) ->
            begin
               monitor = BoltzmannMachines.Monitor()
               BoltzmannMachines.fitrbm(x, epochs = 50,
                     learningrate = 0.005,
                     monitoring = (rbm, epoch) ->
                           BoltzmannMachines.monitorexactloglikelihood!(
                                 monitor, rbm, epoch, datadict))
               monitor
            end);
BoltzmannMachinesPlots.crossvalidationcurve(monitor)

# Determine the optimal number of pretraining epochs for a DBM
# given the other parameters
using Distributed
@everywhere using BoltzmannMachines
@everywhere function my_pretraining_monitoring(
      x::Matrix{Float64}, datadict::DataDict, epoch::Int)

   monitor = Monitor()
   monitorlogproblowerbound!(monitor,
         stackrbms(x,
               nhiddens = [16; 12],
               learningrate = 0.001,
               predbm = true,
               epochs = epoch),
         epoch,
         datadict,
         parallelized = false)
end

monitor = crossvalidation(x, my_pretraining_monitoring, 10:10:200);
BoltzmannMachinesPlots.crossvalidationcurve(monitor, monitorlogproblowerbound)

# ==============================================================================
# Examples for custom optimization algorithm
# ------------------------------------------------------------------------------

x = barsandstripes(100, 16);

struct MyRegularizedOptimizer{R<:AbstractRBM} <: AbstractOptimizer{R}
   llopt::LoglikelihoodOptimizer{R}
end

function MyRegularizedOptimizer()
   MyRegularizedOptimizer(loglikelihoodoptimizer(learningrate = 0.01))
end

import BoltzmannMachines: initialized, computegradient!, updateparameters!

function initialized(myoptimizer::MyRegularizedOptimizer,
      rbm::R) where {R<: AbstractRBM}
   MyRegularizedOptimizer(initialized(myoptimizer.llopt, rbm))
end

function computegradient!(myoptimizer::MyRegularizedOptimizer,
      v, vmodel, h, hmodel, rbm)
   computegradient!(myoptimizer.llopt, v, vmodel, h, hmodel, rbm)
   # add L2 regularization term
   myoptimizer.llopt.gradient.weights .-= 0.1 .* rbm.weights
   nothing
end

function updateparameters!(rbm::R, myoptimizer::MyRegularizedOptimizer{R}
      ) where {R<: AbstractRBM}
   updateparameters!(rbm, myoptimizer.llopt)
end

# The optimizer can be used for fitting RBMs ...
rbm = fitrbm(x; optimizer = MyRegularizedOptimizer())

# and also for DBMs
dbm = fitdbm(x; optimizer = MyRegularizedOptimizer())


text:
using RecipesBase

@recipe f(gp::AbstractGP, x::AbstractArray) = gp(x)
@recipe f(gp::AbstractGP, x::AbstractRange) = gp(x)
@recipe f(gp::AbstractGP, xmin::Real, xmax::Real) = gp(collect(range(xmin, xmax, length=1000)))
@recipe function f(gp::FiniteGP)
    x = gp.x
    f = gp.f
    ms = marginals(gp)

    @series begin
        μ = mean.(ms)
        σ = std.(ms)
        ribbon := σ
        fillalpha --> 0.3
        linewidth --> 2
        x, μ
    end
end

"""
    sampleplot(GP::FiniteGP, samples)

Plot samples from the given `FiniteGP`. Make sure to run `using Plots` before using this
function.

# Example
```julia
using Plots
f = GP(SqExponentialKernel())
sampleplot(f(rand(10)), 10; markersize=5)
```
The given example plots 10 samples from the given `FiniteGP`. The `markersize` is modified
from default of 0.5 to 5.
"""
@userplot SamplePlot
@recipe function f(sp::SamplePlot)
    x = sp.args[1].x
    f = sp.args[1].f
    num_samples = sp.args[2]
    @series begin
        samples = rand(f(x, 1e-9), num_samples)
        seriestype --> :line
        linealpha --> 0.2
        markershape --> :circle
        markerstrokewidth --> 0.0
        markersize --> 0.5
        markeralpha --> 0.3
        seriescolor --> "red"
        label --> ""
        x, samples
    end
end


text:
mutable struct Block # MCU
    data::Matrix{Float32}
end

function split_into_blocks(x::Matrix{Float32}; blocksize=8)::Matrix{Block}
    blocks::Matrix{Block} = map(x->Block(Matrix{Float32}(undef, (blocksize,blocksize))), zeros(Int.(size(x)./blocksize)...))
    for (idxi,i)=enumerate(1:blocksize:size(x)[1]), (idxj,j)=enumerate(1:blocksize:size(x)[2])
        block = blocks[idxi,idxj]
        for (idxii,ii)=enumerate(i:i+(blocksize-1)), (idxjj,jj)=enumerate(j:j+(blocksize-1))
            block.data[idxii,idxjj] = x[ii,jj]
        end
    end
    return blocks
end

function join_blocks(x::Matrix{Block}; blocksize=8)::Matrix{Float32}
    retimg = Matrix{Float32}(undef, Int.(size(x).*blocksize))
    for (idxi,i)=enumerate(1:blocksize:size(retimg)[1]), (idxj,j)=enumerate(1:blocksize:size(retimg)[2])
        block = x[idxi,idxj]
        for (idxii,ii)=enumerate(i:i+blocksize-1), (idxjj,jj)=enumerate(j:j+blocksize-1)
            retimg[ii,jj] = block.data[idxii,idxjj]
        end
    end
    return retimg
end


text:
# four1 test

Fs = 1000;            # Sampling frequency
T = 1/Fs;             # Sampling period

function prepare_input(input_size)
	t = (0:input_size-1)*T;        # Time vector
	S = zeros(2^ceil(log2(input_size)))
	S[1:input_size] = 0.7*sin.(2pi*50t) + sin.(2pi*120t) + 2*randn(size(t));
	return S
end

function test_function(input_data)
	isign = 1
	nn = length(input_data)
	nnlog = log2(nn)
	if ceil(nnlog) != nnlog
		error("Length of input_data vector must be a power of 2.")
	end
	j::Int64 = 1
	res::Array{Complex{Float64},1} = copy(input_data)

	for i = 1:nn
		if j > i
			res[j], res[i] = input_data[i], input_data[j]
		end
		m = nn/2
		while m >= 2 && j > m
			j -= m
			m = m/2
		end
		j += m
	end

	mmax = 1;
	while nn > mmax
		istep = 2*mmax
		theta = isign*(pi/mmax)
		wp = -2*sin(theta/2)^2 + sin(theta)*im
		w = 1 + 0im

		for m = 1:mmax
			for i = m:istep:nn
				j = i + mmax
				temp = w*res[j]
				res[j] = res[i] - temp
				res[i] += temp
			end
			w += w*wp
		end
		mmax = istep
	end

	return res
end


text:
# A logical type map can be provided during schema construction.
# It contains mapping of a column to a logical type and the converter function to be applied.
# Columns can be indentified either by their actual type or column name (the full path in the schema)
const TLogicalTypeMap = Dict{Union{Int32,Vector{String}},Tuple{DataType,Function}}

# schema and helper methods
mutable struct Schema
    schema::Vector{SchemaElement}
    map_logical_types::TLogicalTypeMap
    name_lookup::Dict{Vector{String},SchemaElement}
    type_lookup::Dict{Vector{String},Union{DataType,Union}}
    nttype_lookup::Dict{Vector{String},Union{DataType,Union}}

    function Schema(elems::Vector{SchemaElement}, map_logical_types::TLogicalTypeMap=TLogicalTypeMap())
        name_lookup = Dict{Vector{String},SchemaElement}()
        name_stack = String[]
        nchildren_stack = Int[]

        for idx in 1:length(elems)
            sch = elems[idx]
            nested_name = [name_stack; sch.name]
            name_lookup[nested_name] = sch

            if !haskey(map_logical_types, nested_name)
                if is_logical_string(sch)
                    map_logical_types[nested_name] = (String, logical_string)
                elseif is_logical_timestamp(sch)
                    map_logical_types[nested_name] = (DateTime, logical_timestamp)
                elseif is_logical_decimal(sch)
                    map_logical_types[nested_name] = map_logical_decimal(sch.precision, sch.scale)
                end
            end

            if (idx > 1) && (num_children(sch) > 0)
                push!(nchildren_stack, sch.num_children)
                push!(name_stack, sch.name)
            elseif !isempty(nchildren_stack)
                if nchildren_stack[end] == 1
                    pop!(nchildren_stack)
                    pop!(name_stack)
                else
                    nchildren_stack[end] -= 1
                end
            end
        end
        new(elems, map_logical_types, name_lookup, Dict{Vector{String},Union{DataType,Union}}(), Dict{Vector{String},Union{DataType,Union}}())
    end
end

leafname(schname::T) where {T <: AbstractVector{String}} = [schname[end]]

parentname(schname::T) where {T <: AbstractVector{String}} = istoplevel(schname) ? schname : schname[1:(end-1)]

istoplevel(schname::Vector) = !(length(schname) > 1)

elem(sch::Schema, schname::T) where {T <: AbstractVector{String}} = sch.name_lookup[schname]
function elemindex(sch::Schema, schname::T) where {T <: AbstractVector{String}}
    schema_element = elem(sch, schname)
    findfirst(x->x===schema_element, sch.schema)
end

isrepetitiontype(schelem::SchemaElement, repetition_type) = hasproperty(schelem, :repetition_type) && (schelem.repetition_type == repetition_type)

isrequired(sch::Schema, schname::T) where {T <: AbstractVector{String}} = isrequired(elem(sch, schname))
isrequired(schelem::SchemaElement) = isrepetitiontype(schelem, FieldRepetitionType.REQUIRED)

isoptional(sch::Schema, schname::T) where {T <: AbstractVector{String}} = isoptional(elem(sch, schname))
isoptional(schelem::SchemaElement) = isrepetitiontype(schelem, FieldRepetitionType.OPTIONAL)

isrepeated(sch::Schema, schname::T) where {T <: AbstractVector{String}} = isrepeated(elem(sch, schname))
isrepeated(schelem::SchemaElement) = isrepetitiontype(schelem, FieldRepetitionType.REPEATED)

is_logical_string(sch::SchemaElement) = hasproperty(sch, :_type) && (sch._type === _Type.BYTE_ARRAY) && ((hasproperty(sch, :converted_type) && (sch.converted_type === ConvertedType.UTF8)) || (hasproperty(sch, :logicalType) && hasproperty(sch.logicalType, :STRING)))

# converted_type is usually not set for INT96 types, but they are used exclusively used for timestamps only
is_logical_timestamp(sch::SchemaElement) = hasproperty(sch, :_type) && (sch._type === _Type.INT96)

function is_logical_decimal(sch::SchemaElement)
    if hasproperty(sch, :_type)
        if (sch._type === _Type.FIXED_LEN_BYTE_ARRAY) || (sch._type === _Type.INT64)
            if (hasproperty(sch, :converted_type) && (sch.converted_type === ConvertedType.DECIMAL)) || (hasproperty(sch, :logicalType) && hasproperty(sch.logicalType, :DECIMAL))
                return true
            end
        end
    end
    false
end

function path_in_schema(sch::Schema, schelem::SchemaElement)
    for (n,v) in sch.name_lookup
        (v === schelem) && return n
    end
    error("schema element not found in schema")
end

function logical_converter(sch::Schema, schname::T) where {T <: AbstractVector{String}}
    elem = sch.name_lookup[schname]

    if schname in keys(sch.map_logical_types)
        _logical_type, converter = sch.map_logical_types[schname]
        return converter
    elseif hasproperty(elem, :_type) && (elem._type in keys(sch.map_logical_types))
        _logical_type, converter = sch.map_logical_types[elem._type]
        return converter
    else
        return identity
    end
end

function logical_convert(sch::Schema, schname::T, val) where {T <: AbstractVector{String}}
    elem = sch.name_lookup[schname]

    if schname in keys(sch.map_logical_types)
        logical_type, converter = sch.map_logical_types[schname]
        converter(val)::logical_type
    elseif hasproperty(elem, :_type) && (elem._type in keys(sch.map_logical_types))
        logical_type, converter = sch.map_logical_types[elem._type]
        converter(val)::logical_type
    else
        val
    end
end

elemtype(sch::Schema, schname::T) where {T <: AbstractVector{String}} = get!(sch.type_lookup, schname) do
    elem = sch.name_lookup[schname]

    if schname in keys(sch.map_logical_types)
        logical_type, _converter = sch.map_logical_types[schname]
        logical_type
    elseif hasproperty(elem, :_type) && (elem._type in keys(sch.map_logical_types))
        logical_type, _converter = sch.map_logical_types[elem._type]
        logical_type
    else
        elemtype(elem)
    end
end
function elemtype(schelem::SchemaElement)
    jtype = Nothing

    if hasproperty(schelem, :_type)
        jtype = PLAIN_JTYPES[schelem._type+1]
    else
        jtype = Dict{Symbol,Any} # this is a nested type
    end

    if (hasproperty(schelem, :_type) && (schelem._type == _Type.BYTE_ARRAY || schelem._type == _Type.FIXED_LEN_BYTE_ARRAY)) ||
       (hasproperty(schelem, :repetition_type) && (schelem.repetition_type == FieldRepetitionType.REPEATED))  # array type
        jtype = Vector{jtype}
    end

    jtype
end

ntcolstype(sch::Schema, schname::T) where {T <: AbstractVector{String}} = get!(sch.nttype_lookup, schname) do
    ntcolstype(sch, sch.name_lookup[schname])
end
function ntcolstype(sch::Schema, schelem::SchemaElement)
    @assert num_children(schelem) > 0
    idx = findfirst(x->x===schelem, sch.schema)
    children_range = (idx+1):(idx+schelem.num_children)
    names = [Symbol(x.name) for x in sch.schema[children_range]]
    types = [(num_children(x) > 0) ? ntelemtype(sch, path_in_schema(sch, x)) : elemtype(sch, path_in_schema(sch, x)) for x in sch.schema[children_range]]
    optionals = [isoptional(x) for x in sch.schema[children_range]]
    types = [Vector{opt ? Union{t,Missing} : t} for (t,opt) in zip(types, optionals)]
    NamedTuple{(names...,),Tuple{types...}}
end

ntelemtype(sch::Schema, schname::T) where {T <: AbstractVector{String}} = get!(sch.nttype_lookup, schname) do
    ntelemtype(sch, sch.name_lookup[schname])
end
function ntelemtype(sch::Schema, schelem::SchemaElement)
    @assert num_children(schelem) > 0
    idx = findfirst(x->x===schelem, sch.schema)
    children_range = (idx+1):(idx+schelem.num_children)
    repeated = hasproperty(schelem, :repetition_type) && (schelem.repetition_type == FieldRepetitionType.REPEATED)
    names = [Symbol(x.name) for x in sch.schema[children_range]]
    types = [(num_children(x) > 0) ? ntelemtype(sch, path_in_schema(sch, x)) : elemtype(sch, path_in_schema(sch, x)) for x in sch.schema[children_range]]
    optionals = [isoptional(x) for x in sch.schema[children_range]]
    types = [opt ? Union{t,Missing} : t for (t,opt) in zip(types, optionals)]
    T = NamedTuple{(names...,),Tuple{types...}}
    repeated ? Vector{T} : T
end

bit_or_byte_length(sch::Schema, schname::Vector{String}) = bit_or_byte_length(elem(sch, schname))
bit_or_byte_length(schelem::SchemaElement) = hasproperty(schelem, :type_length) ? schelem.type_length : 0

num_children(schelem::SchemaElement) = hasproperty(schelem, :num_children) ? schelem.num_children : 0

function max_repetition_level(sch::Schema, schname::T) where {T <: AbstractVector{String}}
    lev = isrepeated(sch, schname) ? 1 : 0
    istoplevel(schname) ? lev : (lev + max_repetition_level(sch, parentname(schname)))
end

function max_definition_level(sch::Schema, schname::T) where {T <: AbstractVector{String}}
    lev = isrequired(sch, schname) ? 0 : 1
    istoplevel(schname) ? lev : (lev + max_definition_level(sch, parentname(schname)))
end

logical_decimal_unscaled_type(precision::Int32) = (precision < 5) ? UInt16 :
    (precision < 10) ? UInt32 :
    (precision < 19) ? UInt64 : UInt128

function map_logical_decimal(precision::Int32, scale::Int32; use_float::Bool=false)
    T = logical_decimal_unscaled_type(precision)
    if scale == 0
        # integral values
        return (signed(T), (bytes)->logical_decimal_integer(bytes, T))
    elseif use_float
        # use Float64
        return (Float64, (bytes)->logical_decimal_float64(bytes, T, scale))
    else
        # use Decimal
        return (Decimal, (bytes)->logical_decimal_scaled(bytes, T, scale))
    end
end


text:
export Day1, Day1_1, Day1_2
abstract type Day1 <: Day end
abstract type Day1_1 <: Day1 end
abstract type Day1_2 <: Day1 end

line_parser(::Type{<:Day1}, line) = parse(Int, line)

solve(::Type{Day1_1}, lines::Vector) = solve(Day1_1, lines, 2020, 2)
solve(::Type{Day1_2}, lines::Vector) = solve(Day1_2, lines, 2020, 3)

solve(::Type{T}, lines::Vector, sum_num, count) where T <: Day1 = lines |> i -> find_sum_numbers(i, sum_num, count) |> prod

function find_sum_numbers(items, sum_to_find, how_many)
    sort!(items)
    for nums in Iterators.product([items for i in 1:how_many]...)
        val = sum(nums)
        if val == sum_to_find
            return nums
        end
    end
end


text:
using UnPack
function staticstep!(F, x, η, dη, q_old, GS, VVlp, region)
    if region == 1
        q   = x[1]
        ψ   = x[2]
        ssq = x[3]

        F[1] = log(q) / γ + log(GS) - log(aₑ * ψ + aₕ * (1. - ψ) - (χ₁ * q - 1) / χ₂)
        F[2] = ssq * (q - (q - q_old) * (ψ - η) / dη) - σ * q
        F[3] = aₑ - aₕ - q * (ψ - η) * ssq^2 * VVlp
    elseif region == 2
        q    = x[1]
        F[1] = log(q) / γ + log(GS) - log(aₑ - (χ₁ * q - 1.) / χ₂)
    end
    nothing
end

function fixedpoint_timestep(F, vₑvₕ, η, dη, Q, Qp, Psi, SSQ, MU, S, Svₑ, Svₕ, Mvₑ, Mvₕ, vₑ₁, vₕ₁, θ;
                             use_nlsolve::Bool = false, method::Symbol = :newton, uniform::Bool = true,
                             my_upwind::Bool = false)
    @unpack ρₑ, ρₕ, τ, aₑ, aₕ, σ, γ, χ₁, χ₂, δ, Δ = θ
    N = length(η)

    vₑ = vₑvₕ[1:N]
    vₕ = vₑvₕ[(N + 1):end]
#= # For diagnosing problems. Typically boundary values becoming negative are the problem
    if any(vₑ .< 0.) || any(vₕ .< 0.)
        @show vₑ[1:5]
        @show vₑ[end - 4:end]
        @show vₕ[1:5]
        @show vₕ[end - 4:end]
        @show η[vₑ .< 0.]
        @show η[vₕ .< 0.]
    end
=#

    ## Initialization
    GS1 = (η        ./ vₑ ) .^ (1. / γ)
    GS2 = ((1 .- η) ./ vₕ ) .^ (1. / γ)
    GS  =  GS1 + GS2

    vₑp  = diff(vₑ)  ./ dη
    vₕp = diff(vₕ) ./ dη

    vₑpl  = vₑp  ./ vₑ[2:N]
    vₕpl = vₕp ./ vₕ[2:N]
    VVlp = vₕpl - vₑpl + 1. ./ (η[2:N] .* (1. .- η[2:N]))

    ## Initial condition at η[1]
    Psi[1] = 0.
    SSQ[1] = σ
    qL     = [0.]
    qR     = [aₑ * χ₂ + 1.]
    q₀     = [0.]
    for k = 1:30,
        q₀[1] = (qL[1] + qR[1]) / 2.
        ι     = (χ₁ * q₀[1] - 1.) / χ₂
        A₀   = aₕ - ι

        if log(q₀[1]) / γ + log(GS[1]) > log(A₀)
            qR[1] = q₀[1]
        else
            qL[1] = q₀[1]
        end
    end
    Q[1]  = q₀[1]

    ## Newton method
    #  find q, psi and ssq = σ + σ^q from value functions
    break_pt = Int[0]
    for n = 2:N
        q     = Q[n - 1]
        q_old = Q[n - 1]
        ψ     = Psi[n - 1]
        ssq   = SSQ[n - 1]

        xₙ₋₁  = [q, ψ, ssq]

        # errors given guess
        if use_nlsolve
            EN = nlsolve((F, x) -> staticstep!(F, x, η[n], dη[n - 1], q_old, GS[n], VVlp[n - 1], 1), xₙ₋₁,
                         method = method).zero
        elseif default
            ER = [log(q) / γ + log(GS[n]) - log(aₑ * ψ + aₕ * (1. - ψ) - (χ₁ * q - 1) / χ₂),
                  ssq * (q - (q - q_old) * (ψ - η[n]) / dη[n - 1]) - σ * q,
                  aₑ - aₕ - q * (ψ - η[n]) * ssq^2 * VVlp[n - 1]]

            # matrix of derivatives of errors
            # could shorten it since q_old = q

            QN = zeros(3,3)

            QN[1, :] = [1 / (q * γ) + 1/((aₑ - aₕ) * ψ + aₕ - (q - 1.) / χ₂) / χ₂,
                        -(aₑ - aₕ) / ((aₑ - aₕ) * ψ + aₕ - (q - 1.) / χ₂), 0.]

            QN[2, :] = [ssq * (1. - (ψ - η[n]) / dη[n - 1]) - σ,
                        -ssq * (q - q_old) / dη[n-1],  q - (q - q_old) * (ψ - η[n]) / dη[n-1]]

            QN[3, :] = [- (ψ - η[n]) * ssq^2 * VVlp[n-1],
                        -q * ssq^2 * VVlp[n - 1], -2 * q * (ψ - η[n]) * ssq * VVlp[n - 1]]

            EN = [q, ψ, ssq] - QN \ ER
        end

        # if the boundary of the crisis regime has been reached, we have
        # ψ = 1 from now on
        if EN[2] > 1.
            if EN[1] < 0
                out = nlsolve((F, x) -> staticstep!(F, x, η[n], dη[n - 1], q_old, GS[n], VVlp[n - 1], 1), xₙ₋₁,
                         method = method)
                @show out
            end
            break_pt[1] = n
            break
        end

        # save results
        Q[n]   = EN[1]
        Psi[n] = EN[2]
        SSQ[n] = EN[3]
        Qp[n]  = (Q[n] - q_old) / dη[n-1]
    end

    ## Newton method for ψ = 1, for remaining eta
    for n in break_pt[1]:N
        q = max(Q[n - 1], 1e-6)

        #=ER = log(q) / γ + log(GS[n]) - log(aₑ - (χ₁ * q - 1.) / χ₂)
        QN = 1. / (q * γ)  + 1. / (aₑ - (q - 1.) / χ₂) / χ₂
        EN = q - ER / QN
        EN = find_zero(x -> x^(1 / γ) * GS[n] - (aₑ - (χ₁ * x - 1.) / χ₂), (0., q * 2))=#
        EN = nlsolve(x -> x.^(1 / γ) .* GS[n] - (aₑ .- (χ₁ .* x .- 1.) ./ χ₂), [q], method = :newton).zero[1]

        Qp[n]  = (EN - q) / dη[n-1]
        Q[n]   =  EN
        Psi[n] = 1.
        SSQ[n] = 1. / (1. - (1. - η[n]) * Qp[n] / Q[n]) * σ
    end

    ## Computing the PDE
    # Volatility of η
    S    .= (Psi - η) .* SSQ # This is σ_η * η
    S[1]  = 0.
    S[N]  = 0.
    Iota .= (χ₁ .* Q .- 1.) ./ χ₂
    Φ    .= log.(χ₁ .* Q) / χ₂
    A    .= aₑ .* Psi + aₕ .* (1. .- Psi) - Iota

    # for the law of motion of η, C/N = (η Q)^(1/γ - 1)/V^(1/γ)
    # recall that GS1 = (η./V).^(1/γ)
    CN  = GS1 .* Q .^ (1. / γ - 1.) ./ η # C/N*η
    C_N = GS2 .* Q.^(1. / γ - 1.) ./ (1. .- η)

    Svₑ[2:N] .= vₑpl .* S[2:N]
    Svₕ[2:N] .= vₕpl .* S[2:N]

    VarSig  = -Svₑ + S ./ η        + SSQ .+ (γ - 1) * σ
    VarSig_ = -Svₕ + S ./ (1. .- η) + SSQ .+ (γ - 1) * σ

    ind = 2:(N - 1)
    MU[ind] .= ((aₑ .- Iota[ind]) ./ Q[ind] - CN[ind] .- τ) .* η[ind] + S[ind] .* (VarSig[ind] - SSQ[ind]) # This is μ_η * η

    Mvₑ .= ρₑ .- CN  - (1 - γ) .* (Φ .- δ + Svₑ .* σ .- γ * σ^2 / 2.)
    Mvₕ .= ρₕ .- C_N - (1 - γ) .* (Φ .- δ + Svₕ .* σ .- γ * σ^2 / 2.)

    ## Updating vₑ and vₕ
    # last parameter is dt*ρₑ if dt is small, can be at most 1 (1 corresponds to policy iteration)
    # it is more agressive to set the last parameter closer to 1, but code may not converge

    if my_upwind
        pseudo_transient_relaxation!(stategrid, (vₑ₁, vₕ₁), (vₑ, vₕ), (Mvₑ, Mvₕ),
                                     MU ./ η, (S ./ η).^2, Δ; uniform = uniform)
    else
        vₑ₁ .= upwind_parabolic_pde(η, Mvₑ, MU, S.^2, G, vₑ, Δ)
        vₕ₁ .= upwind_parabolic_pde(η, Mvₕ, MU, S.^2, G, vₕ, Δ)
    end

    F[1:N]         .= vₑ₁ - vₑ
    F[(N + 1):end] .= vₕ₁ - vₕ
end


text:
"""
    SpatialAPI

Module responsible for defining core methods for accessing and manipulating spatial
components data structures in Julia.
"""
module SpatialAPI

using AxisIndices
using AxisIndices.ObservationDims
using AxisIndices.Interface: NamedDimsArray

using MappedArrays
using MappedArrays: AbstractMultiMappedArray, AbstractMappedArray

using GeometryBasics
using TimeAxes

using NeuroCore.ColorChannels
using NeuroCore.ColorChannels: ZeroArray, is_channel

using Base: tail

export
    height,
    pixel_spacing,
    sdims,
    spatialdims,
    spatial_axes,
    spatial_indices,
    spatial_keys,
    spatial_size,
    spatial_directions,
    spatial_order,
    width,
    width_height

widths_height(img::AbstractArray) = (widths(img), height(img))

GeometryBasics.width(img::AbstractArray) = size(img, 2)

GeometryBasics.height(img::AbstractArray) = size(img, 1)

# yes this abuses @pure but each of these is a pure method and it's the only way I
# could get names to propagate through other operations.
Base.@pure function is_spatial(x::Symbol)
    return !(is_time(x) || is_observation(x) || is_channel(x))
end

"""
    spatial_order(x) -> Tuple{Vararg{Symbol}}

Returns the `dimnames` of `x` that correspond to spatial dimensions.
"""
@inline function spatial_order(::T) where {T}
    if has_dimnames(T)
        return _spatial_order(Val(dimnames(T)))
    else
        return AxisIndices.NamedAxes.default_names(Val(min(3, ndims(T))))
    end
end

@generated function _spatial_order(::Val{L}) where {L}
    keep_names = Symbol[]
    i = 1
    itr = 0
    while (itr < 3) && (i <= length(L))
        n = getfield(L, i)
        if is_spatial(n)
            push!(keep_names, n)
            itr += 1
        end
        i += 1
    end
    quote
        return $(keep_names...,)
    end
end

"""
    spatialdims(x) -> Tuple{Vararg{Int}}

Return a tuple listing the spatial dimensions of `img`.
Note that a better strategy may be to use ImagesAxes and take slices along the time axis.
"""
@inline function spatialdims(x::T) where {T}
    if has_dimnames(x)
        return _spatialdims(Val(dimnames(x)))
    else
        return ntuple(+, Val(min(ndims(T), 3)))
    end
end

spatialdims(img::AbstractMappedArray) = spatialdims(parent(img))
function spatialdims(img::AbstractMultiMappedArray)
    ps = traititer(spatialdims, parent(img)...)
    checksame(ps)
end
# TODO remove spatialdims(img::OffsetArray) = coords_spatial(parent(img))
@inline spatialdims(img::SubArray) =
    _subarray_offset(0, spatialdims(parent(img)), img.indices...)
@inline function spatialdims(img::Base.PermutedDimsArrays.PermutedDimsArray{T,N,perm,iperm}) where {T,N,perm,iperm}
    return _getindex_tuple(spatialdims(parent(img)), iperm)
end

@generated function _spatialdims(::Val{L}) where {L}
    keep_names = Int[]
    i = 1
    itr = 0
    while (itr < 3) && (i <= length(L))
        if is_spatial(getfield(L, i))
            push!(keep_names, i)
            itr += 1
        end
        i += 1
    end
    quote
        return $(keep_names...,)
    end
end

"""
    spatial_axes(x) -> Tuple

Returns a tuple of each axis corresponding to a spatial dimensions.
"""
@inline spatial_axes(x) = _filter_axes(named_axes(x), spatial_order(x))

@inline function _filter_axes(naxs::NamedTuple, d::Tuple{Vararg{Any,N}}) where {N}
    return ntuple(i -> getfield(naxs, getfield(d, i)), Val(N))
end

"""
    spatial_axes(x, sz; kwargs...)

Returns an `AxesIterator` along the spatial axes.
"""
@inline function spatial_axes(x, sz; kwargs...)
    return AxesIterator(CartesianIndices(spatial_axes(x)), sz; kwargs...)
end

"""
    spatial_size(x) -> Tuple{Vararg{Int}}

Return a tuple listing the sizes of the spatial dimensions of the image.
"""
@inline spatial_size(x) =  map(i -> size(x, i), spatialdims(x))
spatial_size(img::AbstractMappedArray) = size_spatial(parent(img))
function spatial_size(img::AbstractMultiMappedArray)
    ps = traititer(spatial_size, parent(img)...)
    checksame(ps)
end
# TODO remove spatial_size(img::OffsetArray) = size_spatial(parent(img))
@inline spatial_size(img::SubArray) =
    _subarray_filter(spatial_size(parent(img)), img.indices...)
@inline function spatial_size(img::Base.PermutedDimsArrays.PermutedDimsArray{T,N,perm,iperm}) where {T,N,perm,iperm}
    return _getindex_tuple(spatial_size(parent(img)), iperm)
end

"""
    spatial_indices(x)

Return a tuple with the indices of the spatial dimensions of the
image. Defaults to the same as `indices`, but using `NamedDimsArray` you can
mark some axes as being non-spatial.
"""
@inline spatial_indices(x) = map(values, spatial_axes(x))

"""
    spatial_keys(x)
"""
@inline spatial_keys(x) = map(keys, spatial_axes(x))

"""
    pixel_spacing(x)

Return a tuple representing the separation between adjacent pixels along each axis
of the image. Derived from the step size of each element of `spatial_keys`.
"""
@inline function pixel_spacing(x)
    map(spatial_keys(x)) do ks_i
        if AxisIndices.StaticRanges.has_step(ks_i)
            return step(ks_i)
        else
            return 1
        end
    end
end

"""
    spatial_offset(x)

The offset of each dimension (i.e., where each spatial axis starts).
"""
spatial_offset(x) = map(first, spatial_keys(x))

"""
    sdims(x)

Return the number of spatial dimensions in the image. Defaults to the same as
`ndims`, but with `NamedDimsArray` you can specify that some dimensions correspond
to other quantities (e.g., time) and thus not included by `sdims`.
"""
@inline function sdims(x)
    cnt = 0
    for name in dimnames(x)
        if is_spatial(name)
            cnt += 1
        end
    end
    return cnt
end

"""
    affine_map(x) -> AffineMap

Returns and affine map. By default using `spatial_directions` and `pixel_spacing`
are used to constuct the mapping.
"""
function affine_map(x)
    return AffineMap(_spatial_directions_to_rotation(RotMatrix, spatial_directions(x)),
                     _pixelspacing_to_linearmap(pixel_spacing(x)))
end

function _pixelspacing_to_linearmap(ps::NTuple{2,T}) where {T}
    return @inbounds LinearMap(SVector(Float64(ps[1]), Float64(ps[2]), 0.0))
end

function _pixelspacing_to_linearmap(ps::NTuple{3,T}) where {T}
    return @inbounds LinearMap(SVector(Float64(ps[1]), Float64(ps[2]), Float64(ps[3])))
end

function _spatial_directions_to_rotation(::Type{R}, sd::NTuple{2,NTuple{2,T}}) where {R,T}
    return @inbounds R(SMatrix{3,3,Float64,9}(
        sd[1][1], sd[2][1], 0,
        sd[1][2], sd[2][2], 0,
               0,        0, 1)
    )
end

function _spatial_directions_to_rotation(::Type{R}, sd::NTuple{3,NTuple{3,T}}) where {R,T}
    return @inbounds R(SMatrix{3,3,Float64,9}(
        sd[1][1], sd[2][1], sd[3][1],
        sd[1][2], sd[2][2], sd[3][2],
        sd[1][3], sd[2][3], sd[3][3])
    )
end

"""
    spatial_directions(img) -> (axis1, axis2, ...)

Return a tuple-of-tuples, each `axis[i]` representing the displacement
vector between adjacent pixels along spatial axis `i` of the image
array, relative to some external coordinate system ("physical
coordinates").

By default this is computed from `pixel_spacing`, but you can set this
manually using ImagesMeta.
"""
spatial_directions(img::AbstractArray) = _spatial_directions(pixel_spacing(img))

# FIXME
spatial_directions(img::AbstractMappedArray) = spatial_directions(parent(img))
function spatial_directions(img::AbstractMultiMappedArray)
    ps = traititer(spatial_directions, parent(img)...)
    checksame(ps)
end
@inline function spatial_directions(img::SubArray)
    return _subarray_filter(spatial_directions(parent(img)), getfield(img, :indices)...)
end

@inline function spatial_directions(img::Base.PermutedDimsArray{T,N,perm}) where {T,N,perm}
    return _getindex_tuple(spatial_directions(parent(img)), perm)
end

function _spatial_directions(ps::NTuple{N,Any}) where N
    return ntuple(i->ntuple(d->d==i ? ps[d] : zero(ps[d]), Val(N)), Val(N))
end

@inline _subarray_filter(x, i::Real, inds...) =
    _subarray_filter(tail(x), inds...)
@inline _subarray_filter(x, i, inds...) =
    (x[1], _subarray_filter(tail(x), inds...)...)
_subarray_filter(x::Tuple{}) = ()

@inline _subarray_offset(off, x, i::Real, inds...) =
    _subarray_offset(off-1, tail(x), inds...)
@inline _subarray_offset(off, x, i, inds...) =
    (x[1]+off, _subarray_offset(off, tail(x), inds...)...)
_subarray_offset(off, x::Tuple{}) = ()

@inline _getindex_tuple(t::Tuple, inds::Tuple) =
    (t[first(inds)], _getindex_tuple(t, tail(inds))...)
_getindex_tuple(t::Tuple, ::Tuple{}) = ()

@inline traititer(f, A, rest...) = (f(A), traititer(f, rest...)...)
@inline traititer(f, A::ZeroArray, rest...) = traititer(f, rest...)
traititer(f) = ()

function checksame(t::Tuple)
    val1 = t[1]
    @assert all(p -> p == val1, t)
    return val1
end

###
### Deprecations
###

# deprecations since "The Great Refactoring"

@deprecate pixelspacing pixel_spacing

@deprecate spacedirections spatial_directions

@deprecate coords_spatial spatialdims

@deprecate indices_spatial spatial_axes

@deprecate size_spatial spatial_size

@deprecate spatialorder spatial_order

@deprecate namedaxes named_axes

@deprecate nimages ntime

end


text:
################################################################################
#  Copyright 2020, Marta Vanin, Tom Van Acker                                  #
################################################################################
# PowerModelsDistributionStateEstimation.jl                                    #
# An extention package of PowerModels(Distribution).jl for Static Power System #
# State Estimation.                                                            #
################################################################################
"""
    variable_mc_residual

This is the residual variable, which is later associated to the residual (in)equality constraint(s), depending on the chosen
state estimation criterion.
If bounded, the lower bound is set to zero, while the upper bound defaults to Inf, unless the user provides
a different value in the measurement dictionary.
"""
function variable_mc_residual(  pm::_PMD.AbstractUnbalancedPowerModel;
                                nw::Int=_IM.nw_id_default, bounded::Bool=true,
                                report::Bool=true)

    connections = Dict(i => length(meas["dst"]) for (i,meas) in _PMD.ref(pm, nw, :meas) )

    res = _PMD.var(pm, nw)[:res] = Dict(i => JuMP.@variable(pm.model,
        [c in 1:connections[i]], base_name = "$(nw)_res_$(i)",
        start = _PMD.comp_start_value(_PMD.ref(pm, nw, :meas, i), "res_start", c, 0.0)
        ) for i in _PMD.ids(pm, nw, :meas)
    )

    if bounded
        for i in _PMD.ids(pm, nw, :meas), c in 1:connections[i]
            JuMP.set_lower_bound(res[i][c], 0.0)
            res_max = haskey(_PMD.ref(pm, nw, :meas, i), "res_max") ? meas["res_max"] : Inf
            JuMP.set_upper_bound(res[i][c], res_max)
        end
    end

    report && _IM.sol_component_value(pm,:pmd, nw, :meas, :res, _PMD.ids(pm, nw, :meas), res)
end
"""
    variable_mc_load in terms of power, for ACR and ACP
"""
function variable_mc_load(pm::_PMD.AbstractUnbalancedPowerModel; kwargs...)
    variable_mc_load_active(pm; kwargs...)
    variable_mc_load_reactive(pm; kwargs...)
end

function variable_mc_load_active(pm::_PMD.AbstractUnbalancedPowerModel;
                                 nw::Int=_IM.nw_id_default, bounded::Bool=true, report::Bool=true)

    connections = Dict(i => load["connections"] for (i,load) in _PMD.ref(pm, nw, :load))

    pd = _PMD.var(pm, nw)[:pd] = Dict(i => JuMP.@variable(pm.model,
            [c in connections[i]], base_name="$(nw)_pd_$(i)"
            #,start = _PMD.comp_start_value(_PMD.ref(pm, nw, :load, i), "pd_start", c, 0.0) #findall(idx -> idx == c, connections[i])[1]
        ) for i in _PMD.ids(pm, nw, :load)
    )

    if bounded
        for (i,load) in _PMD.ref(pm, nw, :load)
            if haskey(load, "pmin")
                for (idx, c) in enumerate(connections[i])
                    JuMP.set_lower_bound(pd[i][c], load["pmin"][idx])
                end
            end
            if haskey(load, "pmax")
                for (idx, c) in enumerate(connections[i])
                    JuMP.set_upper_bound(pd[i][c], load["pmax"][idx])
                end
            end
        end
    end

    report && _IM.sol_component_value(pm, :pmd, nw, :load, :pd, _PMD.ids(pm, nw, :load), pd)
end

function variable_mc_load_reactive(pm::_PMD.AbstractUnbalancedPowerModel;
                                   nw::Int=_IM.nw_id_default, bounded::Bool=true, report::Bool=true)

    connections = Dict(i => load["connections"] for (i,load) in _PMD.ref(pm, nw, :load))

    qd = _PMD.var(pm, nw)[:qd] = Dict(i => JuMP.@variable(pm.model,
            [c in connections[i]], base_name="$(nw)_qd_$(i)"
            #,start = _PMD.comp_start_value(_PMD.ref(pm, nw, :load, i), "qd_start", c, 0.0) #findall(idx -> idx == c, connections[i])[1]
        ) for i in _PMD.ids(pm, nw, :load)
    )

    if bounded
        for (i,load) in _PMD.ref(pm, nw, :load)
            if haskey(load, "qmin")
                for (idx, c) in enumerate(connections[i])
                    JuMP.set_lower_bound(qd[i][c], load["qmin"][idx])
                end
            end
            if haskey(load, "qmax")
                for (idx, c) in enumerate(connections[i])
                    JuMP.set_upper_bound(qd[i][c], load["qmax"][idx])
                end
            end
        end
    end

    report && _IM.sol_component_value(pm, :pmd, nw, :load, :qd, _PMD.ids(pm, nw, :load), qd)

end
"""
    variable_mc_load_current, IVR current equivalent of variable_mc_load
"""
function variable_mc_load_current(pm::_PMD.AbstractUnbalancedIVRModel; kwargs...)
    variable_mc_load_current_real(pm; kwargs...)
    variable_mc_load_current_imag(pm; kwargs...)
end

function variable_mc_load_current_real(pm::_PMD.AbstractUnbalancedIVRModel;
                                 nw::Int=_IM.nw_id_default, bounded::Bool=true, report::Bool=true)

    connections = Dict(i => load["connections"] for (i,load) in _PMD.ref(pm, nw, :load))

    crd = _PMD.var(pm, nw)[:crd] = Dict(i => JuMP.@variable(pm.model,
            [c in connections[i]], base_name="$(nw)_crd_$(i)"
            #,start = _PMD.comp_start_value(_PMD.ref(pm, nw, :load, i), "crd_start", c, 0.0)
        ) for i in _PMD.ids(pm, nw, :load)
    )

    report && _IM.sol_component_value(pm, :pmd, nw, :load, :crd, _PMD.ids(pm, nw, :load), crd)
end

function variable_mc_load_current_imag(pm::_PMD.AbstractUnbalancedIVRModel; nw::Int=_IM.nw_id_default, bounded::Bool=true, report::Bool=true, meas_start::Bool=false)

    connections = Dict(i => load["connections"] for (i,load) in _PMD.ref(pm, nw, :load))

    cid = _PMD.var(pm, nw)[:cid] = Dict(i => JuMP.@variable(pm.model,
            [c in connections[i]], base_name="$(nw)_cid_$(i)"
            #,start = _PMD.comp_start_value(_PMD.ref(pm, nw, :load, i), "cid_start",c, 0.0)
        ) for i in _PMD.ids(pm, nw, :load)
    )

    report && _IM.sol_component_value(pm, :pmd, nw, :load, :cid, _PMD.ids(pm, nw, :load), cid)

end
"""
    variable_mc_measurement
checks for every measurement if the measured
quantity belongs to the formulation's variable space. If not, the function
`create_conversion_constraint' is called, that adds a constraint that
associates the measured quantity to the formulation's variable space.
"""
function variable_mc_measurement(pm::_PMD.AbstractUnbalancedPowerModel; nw::Int=_IM.nw_id_default, bounded::Bool=false)
    for i in _PMD.ids(pm, nw, :meas)
        msr_var = _PMD.ref(pm, nw, :meas, i, "var")
        cmp_id = _PMD.ref(pm, nw, :meas, i, "cmp_id")
        cmp_type = _PMD.ref(pm, nw, :meas, i, "cmp")
        connections = get_active_connections(pm, nw, cmp_type, cmp_id)
        if no_conversion_needed(pm, msr_var)
            #no additional variable is created, it is already by default in the formulation
        else
            cmp_type == :branch ? id = (cmp_id, _PMD.ref(pm,nw,:branch, cmp_id)["f_bus"], _PMD.ref(pm,nw,:branch, cmp_id)["t_bus"]) : id = cmp_id
            if haskey(_PMD.var(pm, nw), msr_var)
                push!(_PMD.var(pm, nw)[msr_var], id => JuMP.@variable(pm.model,
                    [c in connections], base_name="$(nw)_$(String(msr_var))_$id"))
            else
                _PMD.var(pm, nw)[msr_var] = Dict(id => JuMP.@variable(pm.model,
                    [c in connections], base_name="$(nw)_$(String(msr_var))_$id"))
            end
            msr_type = assign_conversion_type_to_msr(pm, i, msr_var; nw=nw)
            create_conversion_constraint(pm, _PMD.var(pm, nw)[msr_var], msr_type; nw=nw)
        end
    end
end

function variable_mc_generator_current_se(pm::_PMD.AbstractUnbalancedIVRModel; nw::Int=_IM.nw_id_default, bounded::Bool=true, report::Bool=true, kwargs...)
    #NB: the difference with PowerModelsDistributions is that pg and qg expressions are not created
    _PMD.variable_mc_generator_current_real(pm, nw=nw, bounded=bounded, report=report; kwargs...)
    _PMD.variable_mc_generator_current_imaginary(pm, nw=nw, bounded=bounded, report=report; kwargs...)
end


text:
using Bioview
using Base.Test

TESTS = [
    "errors.jl", "options.jl", "MSA.jl"
]

function main()
    root = dirname(@__FILE__)
    for test in TESTS
        include(test)
    end
end

main()


text:
function count_solns(lim)
    counts = Dict()

    for a in 3:Int(lim * 5/4)
        for b in Int(floor(a/5) + 1):Int(floor((a - 1)/2))

            n = (a - b) * (5 * b - a)

            if n >= lim
                break
            end

            counts[n] = get(counts, n, 0) + 1

        end
    end

    println(count(i -> i == 1, values(counts)))
end

count_solns(50000000)

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
# NOTE: these can be standalone windows, or subplots within a Subplots window

abstract PlotWidget <: Widget

type Plot <: PlotWidget
  widget::PyObject  # BasicPlot
  lines::Vector{PlotItem}
  numLeft::Int
  numRight::Int
  autoscale_x::Bool
  autoscale_y::Bool
end

Plot() = Plot(PLOT.BasicPlot(), PlotItem[], 0, 0, true, true)

Base.string(plt::Plot) = "Plot{lines=$(plt.lines)}"
Base.print(io::IO, plt::Plot) = print(io, string(plt))
Base.show(io::IO, plt::Plot) = print(io, string(plt))

function Base.push!(plt::Plot, idx::Int, x::Number, y::Number)
  push!(plt.lines[idx], x, y)
end

# add one data point for each line, all with the same x value
function Base.push!(plt::Plot, x::Number, ys::Vector)
  @assert length(ys) == length(plt.lines)
  for (i,y) in enumerate(ys)
    push!(plt.lines[i], x, y)
  end
end

# refreshes the plot object
function refresh(plt::Plot)

  plt.widget[:startUpdate](plt.autoscale_x, plt.autoscale_y)

  for l in plt.lines
    ca = get(plt.widget[isleft(l) ? "curvesAxis1" : "curvesAxis2"], l.idx-1)
    updateWidgetData(l, ca)
    # ca[:setData](l.x, l.y)
  end

  plt.widget[:finishUpdate]()
  # plt.widget[:replot]()

  nothing
end

function Base.empty!(plt::Plot)
  for l in plt.lines
    empty!(l)
  end
end

getplot(plt::Plot, c::Int) = plt
getline(plt::Plot, c::Int) = plt.lines[c]

title(plt::Plot, title::AbstractString) = plt.widget[:setPlotTitle](title)
xlabel(plt::Plot, label::AbstractString) = plt.widget[:setXAxisTitle](label)
ylabel(plt::Plot, label::AbstractString) = plt.widget[:setYAxisTitle](label)
yrightlabel(plt::Plot, label::AbstractString) = plt.widget[:setYAxisTitleRight](label)
hidelegend(plt::Plot) = plt.widget[:hideLegend]()
showlegend(plt::Plot) = plt.widget[:showLegend]()

# QwtScaleWidget *qwtsw = myqwtplot.axisWidget(QwtPlot::xBottom);
#   QPalette palette = qwtsw->palette();
#   palette.setColor( QPalette::WindowText, Qt::gray);  // for ticks
#   palette.setColor( QPalette::Text, Qt::gray);                  // for ticks' labels
#   qwtsw->setPalette( palette );

function updatePalette(plt::Plot, color, isbackground::Bool)
  qcolor = convertRGBToQColor(color)
  palette = plt.widget[:palette]()

  # note: 0 is for QPalette::WindowText, which is the foreground (axis borders, ticks),
  #       6 is for QPalette::Text
  #       10 is for QPalette::Window, which is background color (outside of axis canvas)
  enumvals = isbackground ? [10] : [0,6]
  for enumval in enumvals
    palette[:setColor](enumval, qcolor)
  end
  plt.widget[:setPalette](palette)
  plt.widget[:setAutoFillBackground](true)

  if isbackground
    plt.widget[:setCanvasBackground](qcolor)
  end
end

function foreground!(plt::Plot, color)
  updatePalette(plt, color, false)
end

function background!(plt::Plot, color)
  updatePalette(plt, color, true)
end

# ----------------------------------------------------------------

type CurrentPlot
  nullableplot::Nullable{Plot}
end
const CURRENT_PLOT = CurrentPlot(Nullable{Plot}())

function currentPlot()
  # create a new plot if it doesn't exist yet
  isnull(CURRENT_PLOT.nullableplot) && currentPlot!(Plot())
  get(CURRENT_PLOT.nullableplot)
end
currentPlot!(plot::Plot) = (CURRENT_PLOT.nullableplot = Nullable(plot))

# ----------------------------------------------------------------

# kvs is a list of (key,value) tuples, where key is a Symbol.
# valid keys:
#   :x (can be vector, range, or matrix... if matrix, y must be matrix too, and nc must match)
#   :y (can be vector, range, or matrix... if matrix, series must go down the columns)

# you can specify any of the following optional arguments
# Note they are all Symbol's (except label, which is a string)
#   :axis, :color, :label, :linetype, :linestyle, :marker, :markercolor

# you can also specify lists of these values, 1 per series, by adding an s to the end of the symbol:
#   :axiss, :colors, etc

const COLORS = [:black, :blue, :green, :red, :darkGray, :darkCyan, :darkYellow, :darkMagenta,
                :darkBlue, :darkGreen, :darkRed, :gray, :cyan, :yellow, :magenta]
const NUMCOLORS = length(COLORS)

# these are valid choices... first one is default value if unset
const LINE_AXES = (:left, :right)
const LINE_TYPES = (:line, :step, :stepinverted, :sticks, :dots, :none, :heatmap)
const LINE_STYLES = (:solid, :dash, :dot, :dashdot, :dashdotdot)
const LINE_MARKERS = (:none, :ellipse, :rect, :diamond, :utriangle, :dtriangle, :cross, :xcross, :star1, :star2, :hexagon)

const DEFAULT_axis = LINE_AXES[1]
const DEFAULT_color = :auto
const DEFAULT_label = "AUTO"
const DEFAULT_width = 2
const DEFAULT_linetype = LINE_TYPES[1]
const DEFAULT_linestyle = LINE_STYLES[1]
const DEFAULT_marker = LINE_MARKERS[1]
const DEFAULT_markercolor = :auto
const DEFAULT_markersize = 10
const DEFAULT_heatmap_n = 100
const DEFAULT_heatmap_c = (0.15, 0.5)
const DEFAULT_fillto = nothing
const DEFAULT_fillcolor = :auto

const DEFAULT_title = ""
const DEFAULT_xlabel = ""
const DEFAULT_ylabel = ""
const DEFAULT_yrightlabel = ""
const DEFAULT_legend = true

makematrix(i::Int) = zeros(Float64, 0, i)
makematrix{T<:Number}(z::Matrix{T}) = convert(Matrix{Float64}, z)
makematrix{T<:Number}(z::Union{StepRange{T,T},UnitRange{T},FloatRange{T},Vector{T}}) = convert(Matrix{Float64}, reshape(z, length(z), 1))
buildX(Y::Matrix{Float64}) = makematrix(1:size(Y,1))

makeplural(s::Symbol) = Symbol(string(s,"s"))
makedefault(s::Symbol) = Symbol(string("DEFAULT_",s))

convertRGBToQColor(rgb::RGB) = QT.QColor(Float64[f(rgb)*255 for f in (red,green,blue)]...)
convertRGBToQColor(color::Colorant) = QT.QColor(Float64[f(color)*255 for f in (red,green,blue,alpha)]...)

"duplicate a single value, or pass the 2-tuple through"
maketuple(x::Real) = (x,x)
maketuple{T,S}(x::Tuple{T,S}) = x

# get the corresponding plural arg, or the regular arg, or the default
# example: if :colors is set, then grab the i_th color, otherwise if :color is set, return that color, otherwise return the default :auto
function getarg(s::Symbol, d::Dict, c::Int)
  plural = makeplural(s)
  if haskey(d, plural)
    return d[plural][c]
  end
  get(d, s, eval(makedefault(s)))
end

autocolor(idx::Integer) = COLORS[mod1(idx,NUMCOLORS)]

# add one line to plot
function addline(plt::Plot, x, y, color, markercolor, fillcolor,
                                   axis::Symbol, label::AbstractString, width::Int, linetype::Symbol,
                                   linestyle::Symbol, marker::Symbol, markersize::Int,
                                   heatmap_n, heatmap_c::Tuple{Float64,Float64},
                                   tit::AbstractString, xlab::AbstractString, ylab::AbstractString, yrightlab::AbstractString, fillto)

  leftaxis = axis == :left
  isheatmap = linetype == :heatmap

  if leftaxis
    plt.numLeft += 1
    idx = plt.numLeft
  else
    plt.numRight += 1
    idx = plt.numRight
  end

  color = (isa(color, Symbol) ? (color == :auto ? autocolor(idx) : string(color)) : color)
  markercolor = (isa(markercolor, Symbol) ? (markercolor == :auto ? autocolor(idx) : string(markercolor)) : markercolor)
  fillcolor = (isa(fillcolor, Symbol) ? (fillcolor == :auto ? autocolor(idx) : string(fillcolor)) : fillcolor)
  label = string(label == "AUTO" ? "y_$idx" : label, leftaxis ? "" : " (R)")

  # println(color)

  # check our inputs
  # @assert color in COLORS
  # @assert width > 0
  # @assert markersize > 0
  @assert linetype in LINE_TYPES
  @assert linestyle in LINE_STYLES
  @assert marker in LINE_MARKERS
  # @assert markercolor in COLORS
  # @assert heatmap_n > 0
  @assert heatmap_c[1] >= 0.0 && heatmap_c[2] >= heatmap_c[1]

  heatmap_nx, heatmap_ny = maketuple(heatmap_n)

  # create a new plotitem
  plotitem = (isheatmap ? HeatMap(axis, label, idx, heatmap_nx, heatmap_ny, plt) : Series(axis, label, idx, color, plt))

  setdata(plotitem, x, y)
  push!(plt.lines, plotitem)

  tit != "" && title(plt, tit)
  xlab != "" && xlabel(plt, xlab)
  ylab != "" && ylabel(plt, ylab)
  yrightlab != "" && yrightlabel(plt, yrightlab)

  # add it to the figure
  if isheatmap
    plt.widget[:addHeatMap](leftaxis, string(label), heatmap_c...)
  else
    args = map(string, (label, linetype, linestyle, marker))
    plt.widget[:addLine](leftaxis, width, markersize, color, args..., markercolor, fillto, fillcolor)
  end

  plotitem
end

function addRegressionLine(line)
  x = [minimum(line.x), maximum(line.x)]
  reg = [line.x ones(length(line.x))] \ line.y
  y = reg[1] * x + reg[2]
  oplot(line.plt, x, y, label = split(line.label, " (")[1] * " REG", color = getreglinecolor(line), width = 3)
  nothing
end

oplot(plt::PlotWidget, y::AbstractArray; kvs...) = oplot(plt; y = y, kvs...)
oplot(plt::PlotWidget, x::AbstractArray, y::AbstractArray; kvs...) = oplot(plt; x = x, y = y, kvs...)
oplot(plt::PlotWidget, f::Function, x::AbstractArray; kvs...) = oplot(plt; x = x, y = map(f, x), kvs...)
oplot(args...; kwargs...) = oplot(currentPlot(), args...; kwargs...)

# generic way to add to plot
function oplot(plotwidget::PlotWidget; kvs...)
  d = Dict(kvs)

  # create a Matrix{Float64} with series y's going down columns
  # NOTE: you can leave y out of the parameter list and instead set n=numberOfBlankSeries
  Y = makematrix(get(d, :y, get(d, :n, 0)))

  # create a Matrix{Float64} with series x's going down columns
  # X can have 1 column, in which case it is reused for each series
  X = makematrix(get(d, :x, buildX(Y)))

  # duplicate the column of Y to match the columns of X
  xcols = size(X,2)
  ycols = size(Y,2)
  if xcols > 1 && ycols == 1
    Y = repmat(Y, 1, xcols)
  end

  # add the lines
  for c in 1:ycols
    x = X[:,(xcols==1 ? 1 : c)]
    plt = getplot(plotwidget, c)  # get the correct Plot... nop when plotwidget is Plot, gets sp.plots[c] if MySubplot
    # if isempty(x)
    #   continue
    # end

    # convert Colorant to QColor
    color = getarg(:color, d, c)
    if isa(color, Colors.Colorant)
      color = convertRGBToQColor(color)
    end
    markercolor = getarg(:markercolor, d, c)
    if isa(markercolor, Colors.Colorant)
      markercolor = convertRGBToQColor(markercolor)
    end
    fillcolor = getarg(:fillcolor, d, c)
    if isa(fillcolor, Colors.Colorant)
      fillcolor = convertRGBToQColor(fillcolor)
    end

    line = addline(plt, x, Y[:,c], color, markercolor, fillcolor, [getarg(s,d,c) for s in (:axis, :label, :width, :linetype, :linestyle, :marker, :markersize, :heatmap_n, :heatmap_c, :title, :xlabel, :ylabel, :yrightlabel, :fillto)]...)

    if haskey(d, :reg) && d[:reg]
      addRegressionLine(line)
    end
  end

  updateWindow(plotwidget, d)
  plotwidget
end

function updateWindow(plotwidget::PlotWidget, d::Dict)
  if haskey(d, :size)
    resizewidget(plotwidget, d[:size]...)
  end
  if haskey(d, :windowtitle)
    windowtitle(plotwidget, d[:windowtitle])
  end

  if haskey(d, :screen)
    moveToScreen(plotwidget, d[:screen])
  elseif haskey(d, :pos)
    movewidget(plotwidget, d[:pos]...)
  end

  if haskey(d, :legend)
    d[:legend] ? showlegend(plotwidget) : hidelegend(plotwidget)
  end

  if haskey(d, :background_color)
    background!(plotwidget, d[:background_color])
  end
  if haskey(d, :foreground_color)
    foreground!(plotwidget, d[:foreground_color])
  end

  refresh(plotwidget)

  if haskey(d, :show) && d[:show]
    showwidget(plotwidget)
  end
end

plot(y::AbstractArray; kvs...) = plot(; y = y, kvs...)
plot(x::AbstractArray, y::AbstractArray; kvs...) = plot(; x = x, y = y, kvs...)
plot(f::Function, x::AbstractArray; kvs...) = plot(; x = x, y = map(f, x), kvs...)

# plot!(plt::PlotWidget, )

function plot(; kvs...)

  plt = Plot()
  currentPlot!(plt)
  resizewidget(plt, 800, 600)
  moveToLastScreen(plt)  # partial hack so it goes to my center monitor... sorry

  # show the plot (unless show=false)
  if !((:show, false) in kvs)
    push!(kvs, (:show, true))
  end

  oplot(plt; kvs...)
  plt
end

scatter(args...; kwargs...) = plot(args...; kwargs..., linetype=:dots)
heatmap(args...; kwargs...) = plot(args...; kwargs..., linetype=:heatmap)


text:
# Author: Lukasz Bratos
include("methods.jl")
using Test

f(x) = x^2 - 16
pf(x) = 2*x
g(x) = x * exp(-x)
pg(x) = -exp(-x) * (x - 1)
delta = 10^-5
epsilon = 10^-5
maxit = 100

@testset "Metoda bisekcji" begin
    @test abs(Methods.mbisekcji(f, 2.0, 5.5, delta, epsilon)[1]) <= 4 + epsilon
    @test abs(Methods.mbisekcji(f, 2.0, 3.0, delta, epsilon)[4]) == 1
end

@testset "Metoda Newtona" begin
    @test abs(Methods.mstycznych(g, pg, -1.0, delta, epsilon, maxit)[1]) <= 0 + epsilon
    @test abs(Methods.mstycznych(g, pg, 1.0, delta, epsilon, maxit)[4]) == 2
    @test abs(Methods.mstycznych(g, pg, -1.0, delta, epsilon, 1)[4]) == 1
end

@testset "Metoda siecznych" begin
    @test abs(Methods.msiecznych(g, -1.0, 1.0, delta, epsilon, maxit)[1]) <= 0 + epsilon
    @test abs(Methods.msiecznych(g, -1.0, 1.0, delta, epsilon, 2)[4]) == 1
end


text:
struct RDataFrame <: AbstractColumns
   ellist::ElementList
end

istable(::Type{RDataFrame}) = true
columnaccess(::Type{RDataFrame}) = true
columns(df::RDataFrame) = df
columnnames(df::RDataFrame) = getfield(getfield(df, :ellist), :names)

function getcolumn(df::RDataFrame, nm::Symbol)
   getfield(getfield(df, :ellist), :namedelements)[nm]
end

function getcolumn(df::RDataFrame, i::Int)
   getfield(getfield(df, :ellist), :namedelements)[getfield(df, :names)[i]]
end

function getcolumn(df::RDataFrame, ::Type{T}, col::Int, nm::Symbol) where {T}
   getcolumn(df, nm)
end

function r_compatible_type(t::Type{Char})
   Int32
end

function r_compatible_type(t::Type{T}) where T <: AbstractString
   String
end

function r_compatible_type(t::Type{Int32})
   Int32
end

function r_compatible_type(t::Type{Bool})
   Bool
end

function r_compatible_type(t::Type{T}) where {T2, T <: Complex{T2}}
   Complex{T2}
end

function r_compatible_type(t::Type{T}) where T <: Number
   Float64
end

function r_compatible_type(t::Type{Missing})
   Missing
end

function r_compatible_type(t::Type{Union{Missing, T}}) where T
   Union{Missing, r_compatible_type(T)}
end

function convert_to_r_compatible_type(x::AbstractArray{T}) where T
   convert(Array{r_compatible_type(T)}, x)
end

function get_df(obj)
   coldict = Dict{Symbol, Any}()
   attributes = Dict{String, Any}()
   cols = columns(obj)
   colnames = columnnames(cols)
   for col in columnnames(cols)
      coldict[col] = convert_to_r_compatible_type(getcolumn(cols, col))
   end
   colnamesarr = collect(colnames)

   if length(colnamesarr) > 0
      allsame(x) = all(y -> y == first(x), x)
      if !allsame(map(length, values(coldict)))
         error("Lengths of columns not equal")
      end
   end

   RDataFrame(ElementList([], colnamesarr, coldict, attributes))
end