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

#--------  random generation

function GenSymArrow(n::Integer,i::Integer)
    # generates symmetric n x n arrowhad matrix with arrow top at (i,i)
    SymArrow(rand(n-1),rand(n-1),rand(),i)
end

function GenSymDPR1(n::Integer)
    # generates symmetric n x n DPR1 matrix
    SymDPR1(rand(n),rand(n),rand())
end

#-------- Inverses

function inv{T}(A::SymArrow{T},i::Integer,tols::Vector{Float64})

    # COMPUTES: inverse of a SymArrow matrix A, inv(A-A.D[i]*I) which is again SymArrow
    # uses higher precision to compute top of the arrow element accurately, if
    # needed.
    # tols=[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,n)
    z=Array(T,n)
    wz=1/A.z[i]
    shift=A.D[i]

    for k=1:i-1
        D[k]=one(T)/(A.D[k]-shift)
        z[k]=-A.z[k]*D[k]*wz
    end
    for k=i+1:n
        D[k-1]=one(T)/(A.D[k]-shift)
        z[k-1]=-A.z[k]*D[k-1]*wz
    end

    # D=[1./(A.D[1:i-1]-shift),1./(A.D[i+1:end]-shift)]
    # wz=1/A.z[i]
    # z=-[A.z[1:i-1], A.z[i+1:end]].*D1*wz1
    a=A.a-A.D[i]

    # 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+=A.z[k]^2*D[k] : Q+=A.z[k]^2*D[k]
    end
    for k=i:nn
        D[k]>0.0 ? P+=A.z[k+1]^2*D[k] : Q+=A.z[k+1]^2*D[k]
    end
    a<0 ? P=P-a : Q=Q-a

    Kb=(P-Q)/abs(P+Q)
    Kz=maxabs(A.z)*abs(wz)

    if Kb<tols[1] ||  Kz<tols[2]
        b=(P+Q)*wz*wz
    else  # recompute in Double
        Qout=1
        shiftd=map(Double,A.D[i])
        Dd=[Double{Float64}[Double(A.D[k])-shiftd for k=1:i-1];
            Double{Float64}[Double(A.D[k])-shiftd for
                            k=i+1:length(A.D)]]
        wzd=Double(A.z[i])
        ad=Double(A.a)-shiftd

        Pd,Qd=map(Double,(0.0,0.0))

        for k=1:i-1
            Dd[k].hi>0.0 ? Pd+=Double(A.z[k])^2/Dd[k] :
            Qd+=Double(A.z[k])^2/Dd[k]
        end

        for k=i:nn
            Dd[k].hi>0.0 ? Pd+=Double(A.z[k+1])^2/Dd[k] :
            Qd+=Double(A.z[k+1])^2/Dd[k]
            # @show P,Q
        end

        ad.hi<0 ?   Pd=Pd-ad : Qd=Qd-ad

        bd=(Pd+Qd)/(wzd*wzd)
        b=bd.hi+bd.lo
    end

    if i<A.i
        # 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
        for k=n-1:-1:A.i-1
            D[k+1]=D[k]
            z[k+1]=z[k]
        end
        D[A.i-1]=zero(T)
        z[A.i-1]=wz
        SymArrow(D,z,b,i),Kb,Kz,Qout
    else
        # SymArrow([D[1:A.i-1],zero(T),D[A.i:end]],[z[1:A.i-1],wz,z[A.i:end]],b,i+1),Kb,Kz,Qout
        for k=n-1:-1:A.i
            D[k+1]=D[k]
            z[k+1]=z[k]
        end
        D[A.i]=zero(T)
        z[A.i]=wz
        SymArrow(D,z,b,i+1),Kb,Kz,Qout
    end

end # inv

function inv{T}(A::SymArrow{T}, shift::Float64, tolr::Float64)

    # COMPUTES: inverse of the shifted SymArrow A, inv(A-shift*I) which is SymDPR1
    # uses DoubleDouble to compute rho accurately, if needed.
    # tolr is tolerance, usually 1e3,  0.0 forces Double, 1e50 would never use it
    # RETURNS: SymDPR1(D,u,rho), Krho, Qout
    # Krho - condition Krho, Qout = 1 / 0 - Double was / was not used

    n=length(A.D)
    D=Array(T,n+1)
    u=Array(T,n+1)

    # Ds=A.D-shift

    # D=[1./Ds[1:A.i-1];0.0;1./Ds[A.i:end]];
    # u=[A.z[1:A.i-1]./Ds[1:A.i-1];-1.0;A.z[A.i:end]./Ds[A.i:end]];

    for k=1:A.i-1
        D[k]=one(T)/(A.D[k]-shift)
        u[k]=A.z[k]*D[k]
    end
    D[A.i]=zero(T)
    u[A.i]=-one(T)
    for k=A.i:n
        D[k+1]=one(T)/(A.D[k]-shift)
        u[k+1]=A.z[k]*D[k+1]
    end

    # compute rho and Krho
    #--- compute the sum in a plain loop

    a=A.a-shift
    P=zero(T)
    Q=zero(T)
    Qout=0

    for k=1:A.i-1
        D[k]>0.0 ? P=P+A.z[k]^2*D[k] : Q=Q+A.z[k]^2*D[k]
    end
    for k=A.i:n
        D[k+1]>0.0 ? P=P+A.z[k]^2*D[k+1] : Q=Q+A.z[k]^2*D[k+1]
    end

    a<0 ? P=P-a : Q=Q-a

    # Condition of rho
    Krho=(P-Q)/abs(P+Q)

    if Krho<tolr
        rho=-1.0/(P+Q)

    else  # recompute in Double
        Qout=1
        Pd,Qd=map(Double,(0.0,0.0))
        shiftd=Double(shift)
        ad=Double(A.a)-shiftd
        #   Dd=[Double(A.D[k])-shiftd for k=1:length(A.D) ]
        #   Dd=map(Double,A.D)-shiftd
        for k=1:A.i-1
            D[k]>0.0 ? Pd=Pd+Double(A.z[k])^2/(Double(A.D[k])-shiftd) : Qd=Qd+Double(A.z[k])^2/(Double(A.D[k])-shiftd)
        end
        for k=A.i:n
            D[k+1]>0.0 ? Pd=Pd+Double(A.z[k])^2/(Double(A.D[k])-shiftd) : Qd=Qd+Double(A.z[k])^2/(Double(A.D[k])-shiftd)
        end

        # for k=1:length(A.D)
        #    Dd[k].hi>0.0 ? Pd=Pd+Double(A.z[k])^2/Dd[k] : Qd=Qd+Double(A.z[k])^2/Dd[k]
        # end
        ad.hi+ad.lo<0 ? Pd=Pd-ad : Qd=Qd-ad
        r=Double(1.0)/(Pd+Qd)
        rho=-(r.hi+r.lo)
    end

    # returns the following
    SymDPR1(D,u,rho), Krho, Qout

end # inv

function inv{T}(A::SymArrow{T}, shift::Double)

    # COMPUTES: inverse of the shifted SymArrow A, inv(A-shift*I), which is a SymDPR1
    # here shift is Double so it uses Double to compute everything
    # RETURNS: SymDPR1(D1,u1,rho), Qout
    # Qout = 1 on exit meaning Double was used

    n=length(A.D)
    D=Array(Double,n+1)
    u=Array(Double,n+1)

    # D[1:n]=map(Double,A.D)-shift
    for k=1:n
        D[k]=Double(A.D[k])-shift
    end

    u[1:n]=map(Double,A.z)
    a=map(Double,A.a)-shift
    i=A.i

    oned=Double(1.0)
    zerod=Double(0.0)

    # Ds=A.D-shift
    # D=[1./Ds[1:A.i-1];0.0;1./Ds[A.i:end]];
    # u=[A.z[1:A.i-1]./Ds[1:A.i-1];-1.0;A.z[A.i:end]./Ds[A.i:end]];

    for k=1:i-1
        #   tmp=Double(A.D[k])-shift
        #    D[k]=oned/(Double(A.D[k])-shift)
        u[k]=u[k]/D[k]
        D[k]=oned/D[k]
        #   u[k]=Double(A.z[k])/tmp # *D[k]
    end

    for k=i:n
        u[k+1]=u[k]/D[k]
        D[k+1]=oned/D[k+1]
        # D[k+1]=oned/(Double(A.D[k])-shift)
        # u[k+1]=Double(A.z[k])*D[k+1]
    end

    D[i]=zerod
    u[i]=Double(-1.0)

    # D=[Double(1.0)./D1[1:A.i-1];Double(0.0);Double(1.0)./D1[A.i:end]];
    # u=[z1[1:A.i-1]./D1[1:A.i-1];Double(-1.0);z1[A.i:end]./D1[A.i:end]];

    # compute rho and Krho
    #--- compute the sum in a plain loop

    P,Q=zerod,zerod
    Qout=1

    for k=1:i-1
        # D[k].hi+D[k].lo > 0.0  # .lo cannot change sign, I think
        # D[k].hi > 0.0 ? P=P+Double(A.z[k])^2*D[k] : Q=Q+Double(A.z[k])^2*D[k]
        D[k].hi > 0.0 ? P=P+Double(A.z[k])*u[k] : Q=Q+Double(A.z[k])*u[k]
    end
    for k=i:n
        D[k+1].hi >0.0 ? P=P+Double(A.z[k])*u[k+1] : Q=Q+Double(A.z[k])*u[k+1]
    end

    a.hi+a.lo<0.0  ? P=P-a : Q=Q-a
    r=oned/(P+Q)

    # for k=1:length(A.D)
    #    D1[k].hi+D1[k].lo>0.0 ? P=P+z1[k]^2/D1[k] : Q=Q+z1[k]^2/D1[k]
    # end

    rho=-(r.hi+r.lo)

    # returns the following
    D1=Array(T,n+1)
    u1=Array(T,n+1)
    for k=1:n+1
        D1[k]=D[k].hi+D[k].lo
    end
    for k=1:n+1
        u1[k]=u[k].hi+u[k].lo
    end

    # SymDPR1(T[x.hi+x.lo for x=D],T[x.hi+x.lo for x=u],rho), Qout

    SymDPR1(D1,u1,rho), Qout

end # inv

function  eig{T}( A::SymArrow{T},k::Integer,tols::Vector{Float64})

    # COMPUTES: k-th eigenpair of an ordered irreducible SymArrow
    # A = [diag (D) z; z' alpha]
    # tols=[tolb,tolz,tolnu,tolrho,tollambda] = [1e3,10.0*n,1e3,1e3,1e3]
    # RETURNS: lambda, v, Sind, Kb, Kz, Knu, Krho, Qout
    # lambda - k-th eigenvalue in descending order
    # v - lambda's normalized eigenvector
    # Kb, Kz, Knu, Krho - condition numbers
    # Qout = 1 / 0 - Double was / was not used

    # If Qout>0, quad precision was used
    # i was the shift index

    # Set the dimension
    n = length(A.D) + 1

    # Set all conditions initially to zero
    Kb,Kz,Knu,Krho=0.0,0.0,0.0,0.0
    Qout=0
    v=zeros(T,n)
    # Kz is former kappa_nu

    # Determine the shift sigma, the shift index i, and whether lambda
    # is on the left or the right side of the nearest pole
    # Exterior eigenvalues (k = 1 or k = n):

    if k == 1
        sigma,i,side = A.D[1],1,'R'
    elseif k==n
        sigma,i,side = A.D[n-1],n-1,'L'
    else
        # Interior eigenvalues (k in (2,...n-1) ):
        Dtemp = A.D-A.D[k]
        atemp = A.a-A.D[k]
        middle = Dtemp[k-1]/2.0
        Fmiddle = (atemp-middle)-sum(A.z.^2./(Dtemp-middle))
        # middle=(A.D[k-1]-A.D[k])/2.0
        # Fmiddle=0.0
        # for l=1:n-1
        #     Fmiddle+=A.z[l]^2/((A.D[l]-A.D[k])-middle)
        # end
        # Fmiddle=((A.a-A.D[k])-middle)-Fmiddle
        # @show Fmiddle
        sigma,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_i^(-1), Kb and Kz

    Ainv,Kb,Kz,Qout = inv(A,i,tols[1:2])

    # Compute the eigenvalue of the inverse shifted matrix

    nu = bisect( Ainv,side )

    #  nu=fastzero([invD1; 0; invD2], [w1;wz;w2], b, side); # implement later
    #  [nu-nu1]/nu, [nu-nueva]/nu, pause  # and compare

    if abs(nu)==Inf
        # this is nonstandard
        # Deflation in aheig (nu=Inf)
        v[i]=1.0
        lambda=sigma
    else
        # standard case, full computation
        # nu1 is the F- or 1-norm of the inverse of the shifted matrix
        # nu10=maximum([sum(abs(Ainv.z))+abs(Ainv.a), maximum(abs(Ainv.D)+abs(Ainv.z))])
        nu1=0.0
        for l=1:n-1
            nu1=max(nu1,abs(Ainv.D[l])+abs(Ainv.z[l]))
        end
        nu1=max(sumabs(Ainv.z)+abs(Ainv.a), nu1)

        Knu=nu1/abs(nu)
        if Knu<tols[3]
            # Accuracy is fine, compute the eigenvector
            mu = 1.0/nu
            # v=[ A.z./((A.D-sigma)-mu);-1.0]
            for l=1:n-1
            	v[l] = A.z[l]/((A.D[l]-sigma)-mu)
            end
            v[n]=-1.0
            lambda, v = mu+sigma, v/norm(v)
        else
            # Remedies according to Remark 3 - we shift between original
            # eigenvalues and compute DPR1 matrix
            # 1/nu1+sigma, 1/nu+sigma
            # println("Remedy 3 ")
            nu = side=='R' ? abs(nu) : -abs(nu)
            nu1=-sign(nu)*nu1
            sigma1=(nu1+nu)/(2.0*nu*nu1)+sigma

            if findfirst(A.D-sigma1,0.0)>0 # we came back with a pole
                # recompute sigmav more accurately according with dekker
                sigmav=(Double(nu1)+Double(nu))/(Double(2.0)*Double(nu)*Double(nu1))+Double(sigma)
                # Compute the inverse of the shifted arrowhead (DPR1)
                Ainv,Qout1=inv(A,sigmav) # Ainv is Float64
                nu1=bisect(Ainv,side)
                mu1 = 1.0/nu1
                D0=map(Double,A.D)-sigmav
                D1=[D0[l].hi+D0[l].lo for l=1:length(D0)]
                if findfirst(D1-mu1,0.0)>0
                    ind=find(D1-mu1==0.0)
                    v=zeros(n)
                    v[ind]=1.0
                else
                    v=[ A.z./(D1-mu1);-1.0]
                end
                # Shift the eigenvalue back in Double
                lam = Double(1.0)/Double(nu1)+sigmav
                # Return this
                lambda,v = lam.hi+lam.lo, v/norm(v)
            else
                # Compute the inverse of the shifted arrowhead (DPR1)
                Ainv, Krho,Qout1=inv(A,sigma1,tols[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
                nu1= bisect(Ainv,side)
                mu1=1.0/nu1
                # standard v
                v=[ A.z./((A.D-sigma1)-mu1);-1.0]
                # Return this - shift the eigenvalue back and normalize the vector
                lambda, v = mu1+sigma1, v/norm(v)
            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).
        if (abs(A.D[i])+abs(1.0/nu))/abs(lambda)>tols[5]

            if k==1 && A.D[1]<0.0 || k==n && A.D[n-1]>0.0 || i<n-1 && side=='L' && sign(A.D[i])+sign(A.D[i+1])==0 ||
                side=='R' && sign(A.D[i])+sign(A.D[i-1])==0
                # println("Remedy 1 ")
                # Compute the inverse of the original arrowhead (DPR1)
                Ainv,Krho,Qout1 = inv(A,0.0,tols[4]) # Ainv is Float64
                Qout==1 && (Qout=Qout+4)
                if abs(Ainv.r)==Inf
                    lambda=0.0
                else
                    # Here we do not need bisection. We compute the Rayleigh
                    # quotient by using already computed vectors which is
                    # componentwise accurate
                    nu1=sum(v.^2.*Ainv.D)+Ainv.r*sum(v.*Ainv.u)^2;
                    lambda=1.0/nu1
                end
            end

        end
    end

    # Return this
    lambda,v,i,Kb,Kz,Knu,Krho,Qout

end # eig (k)

function eig{T}(A::SymArrow{T}, tols::Vector{Float64})

    # COMPUTES: all eigenvalues and eigenvectors of a real symmetric SymArrow
    # A = [diag(D) z;z' alpha] (notice, here we assume A.i==n)
    # tols = [tolb,tolz,tolnu,tolrho,tollambda] = [1e3,10.0*n,1e3,1e3,1e3] or similar
    # RETURNS: U, E, Sind, Kb, Kz, Knu, Krho, Qout
    # U = eigenvectors, E = eigenvalues in decreasing order
    # Sind[k] - shift index i for the k-th eigenvalue
    # Kb, Kz, Knu, Krho [k] - respective conditions for the k-th eigenvalue
    # Qout[k] = 1 / 0 - Double was / was not used when computing k-th eigenvalue

    n=length(A.D)+1
    n0=n

    # Ordering the matrix
    is=sortperm(A.D,rev=true)
    D=A.D[is]
    z=A.z[is]
    U=eye(T,n,n)
    E=zeros(T,n)
    Kb=zeros(T,n); Kz=zeros(T,n); Knu=zeros(T,n); Krho=zeros(T,n)
    Qout=zeros(Int,n); Sind=zeros(Int,n)

    # Quick return for 1x1
    if n==1
        return U,[A.a],Sind,Kb,Kz,Knu,Krho,Qout
    end

    #  test for deflation in z
    z0=find(z.==0)
    zx=find(z.!=0)

    if isempty(zx)  # nothing to do
        E=[A.D;A.a]
        isE=sortperm(E,rev=true)
        E=E[isE]
        U=U[:,isE]
        return U,E,Sind,Kb,Kz,Knu,Krho,Qout
    end

    if !isempty(z0)
        E[z0]=D[z0]
        D=D[zx]
        z=z[zx]
        if isempty(z)
            E[n]=A.a
            # return U,E,Sind,Kb,Kz,Knu,Krho,Qout
        else
            n=length(z)+1
        end
    end
    zx=[zx;n0]

    #  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=find(g.==0.0)
    gx=find(g.!=0.0)
    if !isempty(g0)
        # Deflation
        Dgx=D[gx]; zgx=z[gx]
        lg0=length(g0)
        R=Array(Tuple{Givens{Float64},Float64},lg0)
        for l=lg0:-1:1
            # This will be changed in v0.4
            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
            E[zx[g0[l]+1]]=D[g0[l]+1]
        end
        # remains
        gx=[0;gx]+1
        nn=length(gx)

        zxx=zx[[gx;n]]
        for k=1:nn+1
            E[zxx[k]],U[zxx,zxx[k]],Sind[zxx[k]],Kb[zxx[k]],Kz[zxx[k]],Knu[zxx[k]],Krho[zxx[k]],Qout[zxx[k]]=
            eig(SymArrow(D[gx],z[gx],A.a,nn+1),k,tols)
        end

        for l=1:lg0
            U=R[l][1]'*U
        end

    else

        # No deflation in D
        for k=1:n
            E[zx[k]],U[zx,zx[k]],Sind[zx[k]],Kb[zx[k]],Kz[zx[k]],Knu[zx[k]],Krho[zx[k]],Qout[zx[k]]=
            eig(SymArrow(D,z,A.a,n),k,tols)
        end
    end
    # end
    # back premutation of vectors
    isi=sortperm(is)

    # must sort E once more
    es=sortperm(E,rev=true)
    E=E[es]

    U=U[[isi[1:A.i-1];n0;isi[A.i:n0-1]],es]

    # Return this
    U,E,Sind[es],Kb[es],Kz[es],Knu[es],Krho[es],Qout[es]

end # eig (all)

function bisect{T}(A::SymArrow{T}, side::Char)
    # COMPUTES: the leftmost (for side='L') or the rightmost (for side='R') eigenvalue
    # of a SymArrow A = [diag (D) z; z'] by bisection.
    # RETURNS: the eigenvalue

    # Determine the starting interval for bisection, [left; right]
    # left, right = side == 'L' ? {minimum([A.D-abs(A.z),A.a-sum(abs(A.z))]), minimum(A.D)} :
    #   {maximum(A.D),maximum([A.D+abs(A.z),A.a+sum(abs(A.z))])}

    absAz = abs(A.z)
    if side == 'L'
        left  = minimum(A.D - absAz)::T
        left  = min(left, A.a - sum(absAz))::T
        right = minimum(A.D)::T
    else
        left  = maximum(A.D)::T
        right = maximum(A.D + absAz)::T
        right = max(right, A.a + sum(absAz))::T
    end

    # Bisection
    middle = (left + right) / convert(T,2)
    z2 = A.z .^ 2
    count, n = 0, length(A.D)

    while (right-left) > 2.0 * eps() * max(abs(left), abs(right))
        # in @time 50% of the time was garbage collection. The following line
        # assigns new vector every time it is called, so it is much better in the
        # loop?? Original timing were 30 secs for n=4000, 2.17 sec for n=1000

        # Fmiddle = A.a-middle-sum(z2./(A.D-middle))

        Fmiddle = zero(T)
        for k=1:n
            Fmiddle += z2[k] / (A.D[k] - middle)
        end
        Fmiddle = A.a - middle - Fmiddle

        if Fmiddle > zero(T)
            left = middle
        else
            right = middle
        end
        middle = (left + right) / convert(T,2)
    end
    # Return the eigenvalue
    right

end # bisect

    function bisect( A::SymDPR1, side::Char )
        # COMPUTES: the leftmost (for side='L') or the rightmost (for side='R') eigenvalue
        # of a SymDPR1 matrix A = diagm(A.D) + A.r*A.u*(A.u)' by bisection.
        # RETURNS: the eigenvalue

        n=length(A.D)
        # Determine the starting interval for bisection, [left; right]
        indD=sortperm(A.D,rev=true)

        if A.r>0.0
            left, right = side == 'L' ? (A.D[indD[n]], A.D[indD[n-1]]) :  (A.D[indD[1]], A.D[indD[1]]+A.r*dot(A.u,A.u))
        else # rho<=0
            left, right = side == 'L' ? (A.D[indD[n]]+A.r*dot(A.u,A.u), A.D[indD[n]]) : (A.D[indD[2]], A.D[indD[1]])
        end

        # Bisection
        middle = (left + right)/2.0
        u2=A.u.^2
        n = length(A.D)
        while (right-left) > 2.0*eps()*maximum([abs(left),abs(right)])

            # Fmiddle = 1.0+A.r*sum(u2./(A.D-middle))

            Fmiddle=0.0
            for k=1:n
                Fmiddle=Fmiddle+u2[k]/(A.D[k]-middle)
            end
            Fmiddle=1.0+A.r*Fmiddle

            sign(A.r)*Fmiddle < 0.0 ? left = middle : right = middle
            middle = (left + right)/2.0
        end
        # Return the eigenvalue
        right

    end # bisect

    # this is just for convenience
    # tols=[1e3,1000.0,1e3,1e3,1e3]


text:
module CUDNN

using ..CUDA

if is_windows()
    const libcudnn = Libdl.find_library(["cudnn64_5","cudnn64_4"])
else
    const libcudnn = Libdl.find_library(["libcudnn"])
end
isempty(libcudnn) && throw("CUDNN library cannot be found.")

const version = ccall((:cudnnGetVersion,libcudnn),Cint,())
const major = div(version, 1000)
const minor = div(version - major*1000, 100)

info("CUDNN version: $(version)")
include("../../lib/$(CUDA.major).$(CUDA.minor)/libcudnn$(major)$(minor).jl")
include("../../lib/$(CUDA.major).$(CUDA.minor)/libcudnn$(major)$(minor)_types.jl")

function checkstatus(status)
    status == CUDNN_STATUS_SUCCESS && return
    Base.show_backtrace(STDOUT, backtrace())
    throw(bytestring(cudnnGetErrorString(status)))
end

datatype(::Type{Float32}) = CUDNN_DATA_FLOAT
datatype(::Type{Float64}) = CUDNN_DATA_DOUBLE
datatype(::Type{Float16}) = CUDNN_DATA_HALF

const handles = Ptr{Void}[]
function handle(x::CuArray)
    dev = device(x) + 1
    while dev > length(handles)
        p = Ptr{Void}[0]
        cudnnCreate(p)
        push!(handles, p[1])
    end
    handles[dev]
end
atexit(() -> foreach(cudnnDestroy, handles))

include("tensor.jl")
include("activation.jl")
#include("batchnorm.jl")
include("convolution.jl")
include("dropout.jl")
include("filter.jl")
##include("lrn.jl")
include("pooling.jl")
include("softmax.jl")
##include("rnn.jl")

end

ext: jl


text:
import NumericExtensions; const ne = NumericExtensions

include("ocl_utils.jl")
include("potential.jl")
include("diffusion.jl")
include("laplacian.jl")
include("align.jl")
include("calcRow.jl")
include("delta.jl")
include("deltaStep2.jl")
include("area.jl")

function createComputeContext{T <: FloatingPoint}(:: Type{T}, ctx, queue)
end

function destroyComputeContext()
end

diffusion!(buff_Xfield, buff_Xpot, buff_Xlap) = diffusionJl!(buff_Xfield, buff_Xpot, buff_Xlap)
potential!(buff_Xfield, buff_Yfield, buff_Zfield, buff_Xpot, buff_Ypot, repulsion) = potentialJl!(buff_Xfield, buff_Yfield, buff_Zfield, buff_Xpot, buff_Ypot, repulsion)
area!(buff_in, buff_out) = areaJl!(buff_in, buff_out, long_direction)
align!(buff_Xfield, buff_Yfield, buff_OUTfield) = alignJl!(buff_Xfield, buff_Yfield, buff_OUTfield, attractionRate, stepIntegration)
laplacian!(buff_in, buff_out) = LaPlacianJl!(buff_in, buff_out)
calcRow!(buff_Xfield, buff_Yfield, buff_Zfield, buff_Wfield, buff_OUT, x, y, z, w) = calcRowJl!(buff_Xfield, buff_Yfield, buff_Zfield, buff_Wfield, buff_OUT, x, y, z, w)
delta!(buff1, buff2, buff3, buff4, buff5, buff6, buff_out, x1, x2, x3, x4, x5, x6) = deltaJl!(buff1, buff2, buff3, buff4, buff5, buff6, buff_out, x1, x2, x3, x4, x5, x6)
delta2!(buff_D, buff_Xfield, buff_Yfield, buff_L,buff_out, decay) = delta2Jl!(buff_D, buff_Xfield, buff_Yfield, buff_L, buff_out, decay)

smul!(X, buff_in, buff_out) = ne.multiply!(copy!(bouff_out, buff_in), X)
c_add!(in1, in2, out) = ne.add!(copy!(out, in1), in2)

#Julia style smul and add
smul!(out, x) = ne.multiply!(out, x)
c_add!(out, in2) = ne.add!(out, in2)

c_read(source) = copy(source)
c_copy!(target, source) = Base.copy!(target, source)

create{T <: FloatingPoint}(:: Type{T}) = Array(T, fieldResY, fieldResX)
create_n4{T <: FloatingPoint}(:: Type{T}) = Array(Number4{T}, fieldResY, fieldResX)
create_const_n4{T <: FloatingPoint}(x :: T) = fill(Number4(x, x, x, x), fieldResY, fieldResX)


text:
using DiffEqJump, StochasticDiffEq, Test

# Requires threads to be effective
# https://github.com/SciML/DifferentialEquations.jl/issues/854
using Base.Threads
@test Threads.nthreads() > 1

function testdrift!(du, u, p, t)
    du[1] = u[1]
end

function testdiffusion!(du, u, p, t)
    du[1] = u[1]
end

function testrate(u, p, t)
    return 0.32
end

function testaffect!(integrator)
    integrator.u[1] += integrator.u[1] * randn()
end

testjump = ConstantRateJump(testrate, testaffect!)

test_sdeprob = SDEProblem(testdrift!,
    testdiffusion!,
    [1.0],
    (0.0, 1.0)
)

test_jumpprob = JumpProblem(test_sdeprob,
    Direct(),
    testjump
)

test_ensemprob = EnsembleProblem(test_jumpprob)

test_ensemsim = solve(test_ensemprob,
    EM(),
    dt=0.01,
    EnsembleThreads();
    trajectories=1_000_000,
    alias_jump = false,
)

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq

max_issues_repo_name: JuliaDiffEq/StochasticDiffEq

max_forks_repo_name: JuliaDiffEq/StochasticDiffEq


text:
"""
Code declaring SpectralTimeSeriesCommonWavelengths and their abstract versions.

Author: Eric Ford
Created: August 2020
"""

""" Abstract type for a time series of spectra that share a common wavelength grid. """
abstract type AbstractSpectralTimeSeriesCommonWavelengths <: AbstractSpectra1D   end

""" Time series of spectra that share a common wavelength grid. """
struct SpectralTimeSeriesCommonWavelengths{T1<:Real,T2<:Real,T3<:Real,AA1<:AbstractArray{T1,1},AA2<:AbstractArray{T2,2},AA3<:AbstractArray{T3,2},
            AA4<:AbstractArray{UnitRange{Int64},1}, InstT<:AbstractInstrument } <: AbstractSpectralTimeSeriesCommonWavelengths
    λ::AA1
    flux::AA2
    var::AA3
    chunk_map::AA4
    inst::InstT
    metadata::MetadataT
end

function SpectralTimeSeriesCommonWavelengths(λ::A1, flux::A2, var::A3, chunk_map::A4, inst::InstT;
        metadata::MetadataT = MetadataT() ) where {
          T1<:Real, T2<:Real, T3<:Real, A1<:AbstractArray{T1,1}, A2<:AbstractArray{T2,2}, A3<:AbstractArray{T3,2},
          A4<:AbstractArray{UnitRange{Int64},1}, InstT<:AbstractInstrument }
    #println("len(λ) = ", length(λ))
    #println("size(flux) = ", size(flux))
    @assert length(λ) == size(flux,1)
    @assert length(λ) == size(var,1)
    @assert 1 <= length(λ)
    SpectralTimeSeriesCommonWavelengths{eltype(λ),eltype(flux),eltype(var),typeof(λ),typeof(flux),typeof(var),typeof(chunk_map),typeof(inst)}(λ,flux,var,chunk_map,inst,metadata)
end

"""   make_spectral_time_series_common_wavelengths_with_selected_times( input, time_idx )
Returns a SpectralTimeSeriesCommonWavelengths, retaining only those times and spectra specified by time_idx.
"""
function make_spectral_time_series_common_wavelengths_with_selected_times(input::STSCWT, time_idx::AbstractVector{T1} ) where { STSCWT<:AbstractSpectralTimeSeriesCommonWavelengths, T1<:Integer }
    @assert minimum(time_idx) >= 1
    @assert maximum(time_idx) <= size(input.flux,2)
    metadata = length(input.metadata) == size(input.flux,2) ? input.metadata[time_idx] : MetadataT()
    SpectralTimeSeriesCommonWavelengths(input.λ, input.flux[:,time_idx], input.var[:,time_idx], input.chunk_map, input.inst, metadata=metadata )
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
# This diagram shows what is meant by "lower", "upper", and "large":
#      +1   +1
#       |    |
# upper |    |
#       |    |
#      -1    |
#            | large
#      +1    |
#       |    |
# lower |    |
#       |    |
#      -1   -1
#
# That is, we are only concerned with 2:1 subdivision of a surface/element.

# Calculate forward projection matrix for discrete L2 projection from large to upper
#
# Note: This is actually an interpolation.
function calc_forward_upper(n_nodes)
  # Calculate nodes, weights, and barycentric weights
  nodes, weights = gauss_lobatto_nodes_weights(n_nodes)
  wbary = barycentric_weights(nodes)

  # Calculate projection matrix (actually: interpolation)
  operator = zeros(n_nodes, n_nodes)
  for j in 1:n_nodes
    poly = lagrange_interpolating_polynomials(1/2 * (nodes[j] + 1), nodes, wbary)
    for i in 1:n_nodes
      operator[j, i] = poly[i]
    end
  end

  return operator
end

# Calculate forward projection matrix for discrete L2 projection from large to lower
#
# Note: This is actually an interpolation.
function calc_forward_lower(n_nodes)
  # Calculate nodes, weights, and barycentric weights
  nodes, weights = gauss_lobatto_nodes_weights(n_nodes)
  wbary = barycentric_weights(nodes)

  # Calculate projection matrix (actually: interpolation)
  operator = zeros(n_nodes, n_nodes)
  for j in 1:n_nodes
    poly = lagrange_interpolating_polynomials(1/2 * (nodes[j] - 1), nodes, wbary)
    for i in 1:n_nodes
      operator[j, i] = poly[i]
    end
  end

  return operator
end

# Calculate reverse projection matrix for discrete L2 projection from upper to large (Gauss version)
#
# Note: To not make the L2 projection exact, first convert to Gauss nodes,
# perform projection, and convert back to Gauss-Lobatto.
function calc_reverse_upper(n_nodes, ::Val{:gauss})
  # Calculate nodes, weights, and barycentric weights for Legendre-Gauss
  gauss_nodes, gauss_weights = gauss_nodes_weights(n_nodes)
  gauss_wbary = barycentric_weights(gauss_nodes)

  # Calculate projection matrix (actually: discrete L2 projection with errors)
  operator = zeros(n_nodes, n_nodes)
  for j in 1:n_nodes
    poly = lagrange_interpolating_polynomials(1/2 * (gauss_nodes[j] + 1), gauss_nodes, gauss_wbary)
    for i in 1:n_nodes
      operator[i, j] = 1/2 * poly[i] * gauss_weights[j]/gauss_weights[i]
    end
  end

  # Calculate Vandermondes
  lobatto_nodes, lobatto_weights = gauss_lobatto_nodes_weights(n_nodes)
  gauss2lobatto = polynomial_interpolation_matrix(gauss_nodes, lobatto_nodes)
  lobatto2gauss = polynomial_interpolation_matrix(lobatto_nodes, gauss_nodes)

  return gauss2lobatto * operator * lobatto2gauss
end

# Calculate reverse projection matrix for discrete L2 projection from lower to large (Gauss version)
#
# Note: To not make the L2 projection exact, first convert to Gauss nodes,
# perform projection, and convert back to Gauss-Lobatto.
function calc_reverse_lower(n_nodes, ::Val{:gauss})
  # Calculate nodes, weights, and barycentric weights for Legendre-Gauss
  gauss_nodes, gauss_weights = gauss_nodes_weights(n_nodes)
  gauss_wbary = barycentric_weights(gauss_nodes)

  # Calculate projection matrix (actually: discrete L2 projection with errors)
  operator = zeros(n_nodes, n_nodes)
  for j in 1:n_nodes
    poly = lagrange_interpolating_polynomials(1/2 * (gauss_nodes[j] - 1), gauss_nodes, gauss_wbary)
    for i in 1:n_nodes
      operator[i, j] = 1/2 * poly[i] * gauss_weights[j]/gauss_weights[i]
    end
  end

  # Calculate Vandermondes
  lobatto_nodes, lobatto_weights = gauss_lobatto_nodes_weights(n_nodes)
  gauss2lobatto = polynomial_interpolation_matrix(gauss_nodes, lobatto_nodes)
  lobatto2gauss = polynomial_interpolation_matrix(lobatto_nodes, gauss_nodes)

  return gauss2lobatto * operator * lobatto2gauss
end

# Calculate reverse projection matrix for discrete L2 projection from upper to large (Gauss-Lobatto
# version)
function calc_reverse_upper(n_nodes, ::Val{:gauss_lobatto})
  # Calculate nodes, weights, and barycentric weights
  nodes, weights = gauss_lobatto_nodes_weights(n_nodes)
  wbary = barycentric_weights(nodes)

  # Calculate projection matrix (actually: discrete L2 projection with errors)
  operator = zeros(n_nodes, n_nodes)
  for j in 1:n_nodes
    poly = lagrange_interpolating_polynomials(1/2 * (nodes[j] + 1), nodes, wbary)
    for i in 1:n_nodes
      operator[i, j] = 1/2 * poly[i] * weights[j]/weights[i]
    end
  end

  return operator
end

# Calculate reverse projection matrix for discrete L2 projection from lower to large (Gauss-Lobatto
# version)
function calc_reverse_lower(n_nodes, ::Val{:gauss_lobatto})
  # Calculate nodes, weights, and barycentric weights
  nodes, weights = gauss_lobatto_nodes_weights(n_nodes)
  wbary = barycentric_weights(nodes)

  # Calculate projection matrix (actually: discrete L2 projection with errors)
  operator = zeros(n_nodes, n_nodes)
  for j in 1:n_nodes
    poly = lagrange_interpolating_polynomials(1/2 * (nodes[j] - 1), nodes, wbary)
    for i in 1:n_nodes
      operator[i, j] = 1/2 * poly[i] * weights[j]/weights[i]
    end
  end

  return operator
end


text:
# Autogenerated wrapper script for Ipopt_jll for x86_64-apple-darwin-libgfortran3
export amplexe, libipopt, libipoptamplinterface

using ASL_jll
using MUMPS_seq_jll
using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Ipopt")
JLLWrappers.@declare_executable_product(amplexe)
JLLWrappers.@declare_library_product(libipopt, "@rpath/libipopt.3.dylib")
JLLWrappers.@declare_library_product(libipoptamplinterface, "@rpath/libipoptamplinterface.3.dylib")
function __init__()
    JLLWrappers.@generate_init_header(ASL_jll, MUMPS_seq_jll, OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_executable_product(
        amplexe,
        "bin/ipopt",
    )

    JLLWrappers.@init_library_product(
        libipopt,
        "lib/libipopt.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libipoptamplinterface,
        "lib/libipoptamplinterface.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Summation
export sum_naive, sum_kbn, sum_oro
export dot_naive, dot_oro

import VectorizationBase

import ..SIMDops
using  ..SIMDops: Vec, vload, vsum, vzero, fptype

using ..EFT: two_sum, fast_two_sum, two_prod

include("accumulators/sum.jl")
include("accumulators/dot.jl")
include("accumulators/compSum.jl")
include("accumulators/compDot.jl")

# T. Ogita, S. Rump and S. Oishi, "Accurate sum and dot product",
# SIAM Journal on Scientific Computing, 6(26), 2005.
# DOI: 10.1137/030601818
@generated function accumulate(x::NTuple{A, AbstractArray{T}},
                                 accType::F,
                                 rem_handling = Val(:scalar),
                                 ::Val{Ushift} = Val(2)
                                 )  where {F, A, T <: Union{Float32,Float64}, Ushift}
    @assert 0 ≤ Ushift < 6
    U = 1 << Ushift

    W, shift = VectorizationBase.pick_vector_width_shift(T)
    sizeT = sizeof(T)
    WT = W * sizeT
    WU = W * U
    V = Vec{W,T}

    quote
        $(Expr(:meta,:inline))
        px = pointer.(x)
        N = length(first(x))
        Base.Cartesian.@nexprs $U u -> begin
            acc_u = accType($V)
        end

        Nshift = N >> $(shift + Ushift)
        offset = 0
        for n in 1:Nshift
            Base.Cartesian.@nexprs $U u -> begin
                xi = vload.($V, px.+offset)
                add!(acc_u, xi...)
                offset += $WT
            end
        end

        rem = N & $(WU-1)
        for n in 1:(rem >> $shift)
            xi = vload.($V, px.+offset)
            add!(acc_1, xi...)
            offset += $WT
        end

        if $rem_handling <: Val{:mask}
            rem &= $(W-1)
            if rem > 0
                mask = VectorizationBase.mask(Val{$W}(), rem)
                xi = vload.($V, px.+offset, mask)
                add!(acc_1, xi...)
            end
        end

        Base.Cartesian.@nexprs $(U-1) u -> begin
            add!(acc_1, acc_{u+1})
        end

        acc = sum(acc_1)

        if $rem_handling <: Val{:scalar}
            offset = div(offset, $sizeT) + 1
            while offset <= N
                @inbounds xi = getindex.(x, offset)
                add!(acc, xi...)
                offset += 1
            end
        end

        sum(acc)
    end
end

sum_naive(x) = accumulate((x,), sumAcc,                   Val(:scalar), Val(3))
sum_kbn(x)   = accumulate((x,), compSumAcc(fast_two_sum), Val(:scalar), Val(2))
sum_oro(x)   = accumulate((x,), compSumAcc(two_sum),      Val(:scalar), Val(2))

dot_naive(x, y) = accumulate((x,y), dotAcc,     Val(:scalar), Val(3))
dot_oro(x, y)   = accumulate((x,y), compDotAcc, Val(:scalar), Val(2))

end


text:
CPLEX_studio_binaries_key="CPLEX_STUDIO_BINARIES"

if haskey(ENV, CPLEX_studio_binaries_key)
    if isdir(ENV[CPLEX_studio_binaries_key])
        @info "Checked 1 location, found 1 location for CPLEX, using $(ENV[CPLEX_studio_binaries_key])."
        println(ENV[CPLEX_studio_binaries_key])
        exit()
    end
end
locations = Dict(
    "~/.local/bin/ibm/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "~/.local/bin/ibm/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "~/.local/bin/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "~/.local/bin/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "~/.local/bin/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "~/.local/bin/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    # ---
    "~/.local/share/ibm/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "~/.local/share/ibm/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "~/.local/share/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "~/.local/share/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "~/.local/share/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "~/.local/share/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    # ---
    "~/.local/lib/ibm/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "~/.local/lib/ibm/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "~/.local/lib/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "~/.local/lib/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "~/.local/lib/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "~/.local/lib/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    # ---
    "/usr/lib/ibm/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/usr/lib/ibm/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "/usr/lib/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/usr/lib/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "/usr/lib/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/usr/lib/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    # ---
    "/usr/share/ibm/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/usr/share/ibm/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "/usr/share/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/usr/share/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "/usr/share/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/usr/share/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    # ---
    "/usr/local/lib/ibm/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/usr/local/lib/ibm/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "/usr/local/lib/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/usr/local/lib/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "/usr/local/lib/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/usr/local/lib/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    # ---
    "/usr/local/share/ibm/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/usr/local/share/ibm/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "/usr/local/share/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/usr/local/share/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "/usr/local/share/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/usr/local/share/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    # ---
    "/opt/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/opt/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "/opt/ibm/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/opt/ibm/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "/opt/ILOG/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/opt/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
    "/opt/ibm/CPLEX_Studio201/cplex/bin/x86-64_linux/" => false,
    "/opt/ibm/CPLEX_Studio1210/cplex/bin/x86-64_linux/" => false,
)
@assert !isempty(locations) "Locations should not be empty before the search."
prel = length(locations)

for (location, cplexpresent) in locations
    locations[location] = cplexpresent || isdir(location)
end

filter!(l -> l.second, locations)

if isempty(locations)
    # try to find other locations
    error("No CPLEX installation found.")
end
@assert !isempty(locations) "Locations should not be empty."

l = length(locations)
loc = first(locations).first

if l > 1
    # Do some special handling for multiple files.
end

@info "Checked $prel location$(prel == 1 ? "" : "s"), found $l location$(l == 1 ? "" : "s") for CPLEX, using \"$loc\"."

println(loc)

max_stars_repo_name: Entze/jkMST

max_issues_repo_name: Entze/jkMST

max_forks_repo_name: Entze/jkMST


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "onnxruntime"
version = v"1.10.0"

# Collection of sources required to complete build
source_sources = [
    GitSource("https://github.com/microsoft/onnxruntime.git", "0d9030e79888d1d5828730b254fedc53c7b640c1")
]

# Bash recipe for building across all platforms
source_script = raw"""
cd $WORKSPACE/srcdir

wget https://github.com/Kitware/CMake/releases/download/v3.22.1/cmake-3.22.1-linux-x86_64.sh
chmod u+x cmake-3.22.1-linux-x86_64.sh
./cmake-3.22.1-linux-x86_64.sh --skip-license --include-subdir
CMAKE=`pwd`/cmake-3.22.1-linux-x86_64/bin/cmake

wget https://github.com/protocolbuffers/protobuf/releases/download/v3.16.1/protoc-3.16.1-linux-x86_64.zip
unzip -d protoc protoc-3.16.1-linux-x86_64.zip

cd onnxruntime/

git submodule update --init --recursive

mkdir -p build
cd build
$CMAKE $WORKSPACE/srcdir/onnxruntime/cmake \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DONNX_CUSTOM_PROTOC_EXECUTABLE=$WORKSPACE/srcdir/protoc/bin/protoc \
    -Donnxruntime_BUILD_SHARED_LIB=ON \
    -Donnxruntime_BUILD_UNIT_TESTS=OFF \
    -Donnxruntime_CROSS_COMPILING=ON
make -j $nproc
make install
install_license $WORKSPACE/srcdir/onnxruntime/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
binaries = Dict(
    Platform("aarch64", "macOS") => "onnxruntime-osx-arm64-$version.tgz",
    Platform("x86_64", "Windows") => "onnxruntime-win-x64-$version.zip",
    Platform("i686", "Windows") => "onnxruntime-win-x86-$version.zip",
)
function source_platform_exclude_filter(p::Platform)
    libc(p) == "musl" ||
    p == Platform("i686", "Linux") ||
    p == Platform("x86_64", "FreeBSD") ||
    p in keys(binaries)
end
source_platforms = supported_platforms(; exclude=source_platform_exclude_filter)
source_platforms = expand_cxxstring_abis(source_platforms)
@info "Re-packaging binaries for:\n$(join(keys(binaries), "\n"))"
@info "Building binaries for:\n$(join(source_platforms, "\n"))"

# The products that we will ensure are always built
products = Product[
    LibraryProduct(["libonnxruntime", "onnxruntime"], :libonnxruntime)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

for (platform, dist_name) in binaries
    binary_sources = []
    binary_script = """
    cd \$WORKSPACE/srcdir

    wget https://github.com/microsoft/onnxruntime/releases/download/v$version/$dist_name
    if [[ $dist_name == *.zip ]]; then
        unzip $dist_name
    else
        tar xfz $dist_name
    fi
    mkdir -p \$includedir \$libdir
    cp -a onnxruntime*/include/* \$includedir
    if [[ \${target} == *w64* ]]; then
        chmod 775 onnxruntime*/lib/*
    fi
    cp -a onnxruntime*/lib/* \$libdir
    install_license onnxruntime*/LICENSE
    """
    binary_platforms = [platform]
    build_tarballs(ARGS, name, version, binary_sources, binary_script, binary_platforms, products, dependencies;
        preferred_gcc_version = v"8",
        julia_compat = "1.6")
end

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, source_sources, source_script, source_platforms, products, dependencies;
    preferred_gcc_version = v"8",
    julia_compat = "1.6")


text:
## Link

@adddevice(Link, ActiveBranch, LinearExpansionForm{LinearDispatchForm},
           :links, (:L, :T=>:snapshots), joinpath(@__DIR__, "attrs", "links.csv"))

function addtoexpr!(injection::AxisArray, ::NodalActivePower,
                    m::EnergyModel, d::Link)
    B, T = AxisArrays.axes(injection)
    L = axis(m, d)

    bus0 = get(d, :bus0, L)
    bus1 = get(d, :bus1, L)
    p = get(d, :p, L, T)
    eff = get(d, :efficiency, L, T)

    for (i, idx0, idx1) in zip(1:length(L), indexin(bus0, B), indexin(bus1, B)), t in T
        !isnothing(idx0) && add_to_expression!(injection[idx0,t], -1, p[i,t])
        !isnothing(idx1) && add_to_expression!(injection[idx1,t], eff[i,t], p[i,t])
    end
end

cost(d::Link) = sum(d[:marginal_cost] .* d[:p]) + sum(d[:capital_cost] .* (d[:p_nom] - getparam(d, :p_nom)))

function addto!(jm::ModelView, m::EnergyModel, d::Link{DF}) where {DDF, DF <: LinearExpansionForm{DDF}}
    addto!(jm, m, with_formulation(d, LinearExpansionInvestmentForm))
    addto!(jm, m, with_formulation(d, LinearExpansionDispatchForm{DDF}))
end

function addto!(jm::ModelView, m::EnergyModel, d::Link{DF}) where DF <: LinearExpansionInvestmentForm
    L = axis(d)

    p_nom_min = get(d, :p_nom_min, L)
    p_nom_max = get(d, :p_nom_max, L)

    @variable(jm, p_nom_min[l] <= p_nom[l=L] <= p_nom_max[l])
end

function addto!(jm::ModelView, m::EnergyModel, d::Link{DF}) where
      {DDF, DF <: LinearExpansionDispatchForm{DDF}}

    L = axis(m, d)
    T = axis(m, :snapshots)

    p_nom = get(d, :p_nom)
    p_min_pu = get(d, :p_min_pu, L, T)
    p_max_pu = get(d, :p_max_pu, L, T)

    @variable(jm, p[l=L,t=T])

    @constraints jm begin
        p_lower[l=L,t=T], p[l,t] >= p_min_pu[l,t] * p_nom[l]
        p_upper[l=L,t=T], p[l,t] <= p_max_pu[l,t] * p_nom[l]
    end
end

function addto!(jm::ModelView, m::EnergyModel, d::Link{DF}) where
      {DF <: DispatchForm}

    L = axis(m, d)
    T = axis(m, :snapshots)

    p_nom = get(d, :p_nom)
    p_min_pu = get(d, :p_min_pu, L, T)
    p_max_pu = get(d, :p_max_pu, L, T)

    @variable(jm, p_min_pu[l,t] * p_nom[l] <= p[l=L,t=T] <= p_max_pu[l,t] * p_nom[l])
end

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
# CheckpointCondition
abstract type CheckpointCondition end

struct CheckpointAny <: CheckpointCondition end
mutable struct CheckpointLowest <: CheckpointCondition
    lowest::Real
end
CheckpointLowest() = CheckpointLowest(Inf)

(::CheckpointAny)(loss) = true

function (checklowest::CheckpointLowest)(loss)
    cond = loss < checklowest.lowest
    if cond
        checklowest.lowest = loss
    end
    return cond
end

# Checkpointer
struct Checkpointer <: SafeCallback
    condition::CheckpointCondition
    deleteprevious::Bool
end
Checkpointer(condition = CheckpointLowest(); deleteprevious = false) = Checkpointer(
    condition, deleteprevious)

# TODO: refactor to use `Metrics`
function on(::EpochEnd, ::ValidationPhase, checkpointer::Checkpointer, learner)
    loss = epochvalue(getloss(learner.callbacks))
    if checkpointer.condition(loss)
        if checkpointer.deleteprevious
            previousfiles = glob("model-chckpnt-E*", artifactpath(learner))
            foreach(rm, previousfiles)
        end
        epochs = learner.cbstate.history.epochs
        filename = "model-chckpnt-E$epochs-L$loss.bson"
        path = joinpath(artifactpath(learner), filename)
        savemodel(learner.model, path)
    end
end

stateaccess(::Checkpointer) = (callbacks = Read(), cbstate = (history = Read()))


text:
__precompile__()
module Video4Linux

function __init__()
    @show depfile = joinpath(dirname(@__FILE__),"../deps/loadpath.jl")
    isfile(depfile) ? include(depfile) : error("Video4Linux not properly installed. Please run: Pkg.build(\"Video4Linux\")")
end

export
#wrappers
set_io_method,
open_device,
init_device,
start_capturing,
mainloop,
stop_capturing,
uninit_device,
close_device,
copy_buffer_bytes,

#convertions TODO: extend convert
convertY10BtoU16,
convertUYVYtoYCbCr,
convertYUYVtoYCbCr,
#decoders

#videotasks
videoproducer,
stopVideoProducer

include("v4l2header.jl")

include("wrappers.jl")
include("encodings.jl")
include("videotasks.jl")

end # module


text:
"""
GeometricProblems.jl is a collection of ODEs and DAEs with interesting geometric structures.
"""
module GeometricProblems

    include("diagnostics.jl")
    include("plot_recipes.jl")

    include("double_pendulum.jl")
    include("duffing_oscillator.jl")
    include("exponential_growth.jl")
    include("harmonic_oscillator.jl")
    include("kubo_oscillator.jl")
    include("lennard_jones_oscillator.jl")
    include("lorenz_attractor.jl")
    include("lotka_volterra_2d.jl")
    include("lotka_volterra_2d_gauge.jl")
    include("lotka_volterra_2d_singular.jl")
    include("lotka_volterra_2d_symmetric.jl")
    include("lotka_volterra_2d_plots.jl")
    include("lotka_volterra_3d.jl")
    include("lotka_volterra_3d_plots.jl")
    include("lotka_volterra_4d.jl")
    include("lotka_volterra_4d_lagrangian.jl")
    include("lotka_volterra_4d_plots.jl")
    include("massless_charged_particle.jl")
    include("massless_charged_particle_plots.jl")
    include("mathews_lakshmanan_oscillator.jl")
    include("morse_oscillator.jl")
    include("pendulum.jl")
    include("point_vortices.jl")
    include("point_vortices_linear.jl")

end


text:
using QuNet
using LightGraphs, SimpleWeightedGraphs, GraphPlot, MetaGraphs, Plots, Colors, Statistics
using LaTeXStrings
using JLD
using Parameters

num_trials = 1::Int64
grid_size = 5::Int64
time_depth = 5::Int64
num_pairs = 12::Int64
asynchronus_weight = 100::Int64
####### USERPAIR ########

# Generate ixi graph and extend it in time
G = GridNetwork(grid_size, grid_size)

raw_data = []
for j in 1:num_trials
    # Extend in time without memory
    T = QuNet.TemporalGraph(G, time_depth, memory_prob=0.0)
    # Extend in time with memory links:
    T_mem = QuNet.TemporalGraph(G, time_depth, memory_prob=1.0, memory_costs = unit_costvector())

    # Get i random userpairs. Ensure src nodes are fixed on T=1, dst nodes are asynchronus.
    mem_user_pairs = make_user_pairs(T, num_pairs, src_layer=1, dst_layer=-1)
    user_pairs = make_user_pairs(T, num_pairs, src_layer=-1, dst_layer=-1)

    # Add async nodes
    QuNet.add_async_nodes!(T_mem, mem_user_pairs, ϵ=asynchronus_weight)
    QuNet.add_async_nodes!(T, user_pairs, ϵ=asynchronus_weight)

    # Get pathset data
    pathset_mem, dum1, dum2 = QuNet.greedy_multi_path!(T_mem, QuNet.purify, mem_user_pairs, 4)
    pathset, dum1, dum2 = QuNet.greedy_multi_path!(T, QuNet.purify, user_pairs, 4)

    # Pathset is an array of vectors containing edges describing paths between end-user pairs
    # Objective: find the largest timedepth used in the pathsets

    max_depth_mem = QuNet.max_timedepth(pathset_mem, T)
    max_depth = QuNet.max_timedepth(pathset, T)

    # Get the ratio of these two quantities. Add it to data array
    push!(raw_data, max_depth / max_depth_mem )
end

println("")
println(mean(raw_data))

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
"""Struct for the ideal gas thermo.

nasa_low: Array with low temperature nasa coeff. for each species

nasa_high: Array with high temperature nasa coeff. for each species

Trange: Array with temperature ranges for each species

isTcommon: bool which indicates if both polynoms share same T at intersection

"""
struct IdealGasThermo <: Thermo
    nasa_low::Array{Float64,2}
    nasa_high::Array{Float64,2}
    Trange::Array{Float64,2}
    isTcommon::Bool

end
"""
Constructor for the idealGasThermo:

yaml:: Dict of the input yaml file
"""
function IdealGasThermo(yaml::AbstractDict)
    n_species, n_reactions, species_names,
    elements, n_elements, ele_matrix = read_species_basics(yaml)

    nasa_low = zeros(n_species, 7)
    nasa_high = zeros(n_species, 7)
    Trange = zeros(n_species, 3)

    _species_names =
        [yaml["species"][i]["name"] for i = 1:length(yaml["species"])]

    for (i, species) in enumerate(species_names)
        spec = yaml["species"][findfirst(x -> x == species, _species_names)]
        nasa_low[i, :] = spec["thermo"]["data"][1]
        nasa_high[i, :] = spec["thermo"]["data"][2]
        Trange[i, :] .= spec["thermo"]["temperature-ranges"]

        for j = 1:n_elements
            if haskey(spec["composition"], elements[j])
                ele_matrix[j, i] = spec["composition"][elements[j]]
            end
        end
    end
    isTcommon = (maximum(Trange[:, 2]) - minimum(Trange[:, 2])) < 0.01
    return  IdealGasThermo(nasa_low, nasa_high, Trange, isTcommon)
end
"""
    cal_h_RT(gas, T, p, X)

calculates the dimensionless mole based enthalpy (h) for each species
"""
function cal_h_RT(gas::Solution,thermo::IdealGasThermo, T::Real, p::Real, X::AbstractArray)
    H_T = [1.0, T / 2.0, T^2 / 3.0, T^3 / 4.0, T^4 / 5.0, 1.0 / T]
    if T <= 1000.0
        h_mole = @view(thermo.nasa_low[:, 1:6]) * H_T
    else
        h_mole = @view(thermo.nasa_high[:, 1:6]) * H_T
    end
    if !thermo.isTcommon
        ind_correction = @. (T > 1000.0) & (T < thermo.Trange[:, 2])
        h_mole[ind_correction] .=
            @view(thermo.nasa_low[ind_correction, 1:6]) * H_T
    end
    # H_mole = dot(h_mole, X)
    return h_mole
end

"""
    cal_s0_R(gas, T, p, X)

calculates the dimensionless mole based reference state entropy (s0) for each species
"""
function cal_s0_R(gas::Solution,thermo::IdealGasThermo, T::Real, p::Real, X::AbstractArray)
    S_T = [log(T), T, T^2 / 2.0, T^3 / 3.0, T^4 / 4.0, 1.0]
    if T <= 1000.0
        S0 = @view(thermo.nasa_low[:, [1, 2, 3, 4, 5, 7]]) * S_T
    else
        S0 = @view(thermo.nasa_high[:, [1, 2, 3, 4, 5, 7]]) * S_T
    end
    if !thermo.isTcommon
        ind_correction = @. (T > 1000.0) & (T < thermo.Trange[:, 2])
        S0[ind_correction] .=
            @view(thermo.nasa_low[ind_correction, [1, 2, 3, 4, 5, 7]]) *
            S_T
    end
    return S0
end
export cal_s0_R

"""
    cal_s_R(gas, T, p, X)

calculates the dimensionless mole based entropy (s) for each species
"""
function cal_s_R(gas::Solution,thermo::IdealGasThermo, T::Real, p::Real, X::AbstractArray)
   return cal_s0_R(gas,thermo, T,p,X) - log.(max.(X,1e-30)) .- log(p/one_atm)
end

"""
    cal_g_RT(gas, T, p, X)

calculates the dimensionless mole based free gibbs energy (g) for each species
"""
function cal_g_RT(gas::Solution,thermo::IdealGasThermo, T::Real, p::Real, X::AbstractArray)
    return cal_h_RT(gas,thermo, T, p, X) - cal_s_R(gas,thermo, T, p, X)
end

"""
    cal_u_RT(gas, T, p, X)

calculates the dimensionless mole based internal energy (u) for each species
"""
function cal_u_RT(gas::Solution,thermo::IdealGasThermo, T::Real, p::Real, X::AbstractArray)
    return cal_h_RT(gas,thermo, T, p, X) .- 1
end

"""
    cal_a_RT(gas, T, p, X)

calculates the dimensionless mole based helmholz free energy (a) for each species
"""
function cal_a_RT(gas::Solution,thermo::IdealGasThermo, T::Real, p::Real, X::AbstractArray)
    return cal_u_RT(gas,thermo, T, p, X) - cal_s_R(gas,thermo, T, p, X)
end

"""
    cal_cp_R(gas, T, p, X)

calculates the dimensionless mole based heat capacity
at constant pressure (cp) for each species
"""
function cal_cp_R(gas::Solution,thermo::IdealGasThermo, T::Real, p::Real, X::AbstractArray)
    cp_T = [1.0, T, T^2, T^3, T^4]
    if T <= 1000.0
        cp = @view(thermo.nasa_low[:, 1:5]) * cp_T
    else
        cp = @view(thermo.nasa_high[:, 1:5]) * cp_T
    end
    # TODO: not sure if inplace operation will be an issue for AD
    if !thermo.isTcommon
        ind_correction = @. (T > 1000.0) & (T < thermo.Trange[:, 2])
        cp[ind_correction] .=
            @view(thermo.nasa_low[ind_correction, 1:5]) * cp_T
    end
    return cp
end

"""
    cal_cv_R(gas, T, p, X)

calculates the dimensionless mole based heat capacity
at constant volume (cv) for each species
"""
function cal_cv_R(gas::Solution,thermo::IdealGasThermo, T::Real, p::Real, X::AbstractArray)
   return cal_cp_R(gas,thermo, T, p, X) .- 1
end


text:
export sites, edges, surround, faces, neighbors

# This file defines the interface of a lattice iterator
# All lattice iterators are duck-typed

"""
    sites(lattice) -> iterator

Returns an iterator for all sites on the lattice.
"""
function sites end

"""
    edges(lattice; [order=1]) -> iterator

Returns an iterator for all edges of the lattice of a given `order`.
For example, `order = 1` means nearest neighbor edges and `order = 2`
means next-nearest neighbors edges.
"""
function edges end

"""
    surround(lattice) -> iterator

Returns an iterator of surroundings of all sites.

    surround(lattice, site) -> iterator

Returns an iterator of surrounding sites of given `site`.
"""
function surround end

"""
    faces(lattice) -> iterator

Returns an iterator of all faces.
"""
function faces end

"""
    neighbors(lattice, s; order=1) -> iterator

Returns an iterator of the surrounding sites of given site `s`.
"""
function neighbors end


text:
using HTrg
using Test
using Geodesy

@testset "HTrg.jl" begin
    # Write your own tests here.

    p1 = LLA(51.5214588, -0.1729636)
    e1 = ECEFfromLLA(wgs84)(p1)

    code, city, dis, ci = reverse_geocode(e1)
    @test code == "GB"
    @test city == "Bayswater"
    @test isapprox(dis, 1388.453678654349)

    p2 = LLA(22.2710, 113.5767)
    e2 = ECEFfromLLA(wgs84)(p2)
    @show code, city, dis, ci = reverse_geocode(e2)
    @test code == "CN"
    @test city == "Zhuhai"
    @test isapprox(dis, 1130.3958719970858)

end


text:
"Generate a file name to be given to the parser (will show up in stack traces)."
pluto_filename(notebook::Notebook, cell::Cell)::String = notebook.path * "#==#" * string(cell.uuid)

"Is Julia new enough to support filenames in parsing?"
const can_insert_filename = (Base.parse_input_line("1;2") != Base.parse_input_line("1\n2"))

"Parse the code from `cell.code` into a Julia expression (`Expr`). Equivalent to `Meta.parse_input_line` in Julia v1.3, no matter the actual Julia version.

1. Turn multiple expressions into a error expression.
2. Will always produce an expression of the form: `Expr(:toplevel, LineNumberNode(..), root)`. It gets transformed (i.e. wrapped) into this form if needed. A `LineNumberNode` contains a line number and a file name. We use the cell UUID as a 'file name', which makes the stack traces easier to interpret. (Otherwise it would be impossible to tell from which cell a stack frame originates.) Not all Julia versions insert these `LineNumberNode`s, so we insert it ourselves if Julia doesn't.
3. Apply `preprocess_expr` (below) to `root` (from rule 2)."
function parse_custom(notebook::Notebook, cell::Cell)::Expr
    # 1.
    raw = if can_insert_filename
        filename = pluto_filename(notebook, cell)
        ex = Base.parse_input_line(cell.code, filename=filename)
        if (ex isa Expr) && (ex.head == :toplevel)
            # if there is more than one expression:
            if count(a -> !(a isa LineNumberNode), ex.args) > 1
                Expr(:error, "extra token after end of expression")
            else
                ex
            end
        else
            ex
        end
    else
        # Meta.parse returns the "extra token..." like we want, but also in cases like "\n\nx = 1\n# comment", so we need to do the multiple expressions check ourselves after all
        parsed1, next_ind1 = Meta.parse(cell.code, 1, raise=false)
        parsed2, next_ind2 = Meta.parse(cell.code, next_ind1, raise=false)

        if parsed2 === nothing
            # only whitespace or comments after the first expression
            parsed1
        else
            Expr(:error, "extra token after end of expression")
        end
    end

    # 2.
    topleveled = if ExploreExpression.is_toplevel_expr(raw)
        raw
    else
        filename = pluto_filename(notebook, cell)
        Expr(:toplevel, LineNumberNode(1, Symbol(filename)), raw)
    end

    # 3.
    Expr(topleveled.head, topleveled.args[1], preprocess_expr(topleveled.args[2]))
end

"Make some small adjustments to the `expr` to make it work nicely inside a timed, wrapped expression:

1. If `expr` is a `:toplevel` expression (this is the case iff the expression is a combination of expressions using semicolons, like `a = 1; b` or `123;`), then it gets turned into a `:block` expression. The reason for this transformation is that `:toplevel` does not return/relay the output of its last argument, unlike `begin`, `let`, `if`, etc. (But we want it!)
2. If `expr` is a `:module` expression, wrap it in a `:toplevel` block - module creation needs to be at toplevel. Rule 1. is not applied."
function preprocess_expr(expr::Expr)
    if expr.head == :toplevel
		Expr(:block, expr.args...)
    elseif expr.head == :module
        Expr(:toplevel, expr)
    else
        expr
    end
end

# for expressions that are just values, like :(1) or :(x)
preprocess_expr(val::Any) = val

"Wrap `expr` inside a timing block."
function timed_expr(expr::Expr)::Expr
    # @assert ExploreExpression.is_toplevel_expr(expr)

    linenumbernode = expr.args[1]
    root = expr.args[2] # pretty much equal to what `Meta.parse(cell.code)` would give

    # we don't use `quote ... end` here to avoid the LineNumberNodes that it adds (these would taint the stack trace).
    Expr(:block,
        :(local elapsed_ns = time_ns()),
        linenumbernode,
        :(local result = $root),
        :(elapsed_ns = time_ns() - elapsed_ns),
        :((result, elapsed_ns))
    )
end

"Wrap `expr` inside a timing block, and then inside a try ... catch block."
function trycatch_expr(expr::Expr, module_name::Symbol, cell_id::UUID)
    quote
        ans, runtime = try
            # We eval `expr` in the global scope of the workspace module:
            Core.eval($(module_name), $(timed_expr(expr) |> QuoteNode))
        catch ex
            bt = stacktrace(catch_backtrace())
            (CapturedException(ex, bt), missing)
        end
        setindex!(Main.PlutoRunner.cell_results, WeakRef(ans), $(cell_id))
    end
end


text:
struct UnableToParseSSHKey <: Exception
    message::AbstractString
end
Base.show(io::IO, e::UnableToParseSSHKey) = println(io, e.message)

struct UnableToDetectCIService <: Exception
    message::AbstractString
end
Base.show(io::IO, e::UnableToDetectCIService) = println(io, e.message)


text:
using OpenQuantumBase, Test, Random
import LinearAlgebra: Diagonal, diag

#= =============== Define  test functions ================= =#
# calculate the AME linblad term using the formula in reference paper
function ame_update_term(op, ρ, w, v, γ, S)
    L_ij = Array{Array{Complex{Float64},2},2}(undef, (4, 4))
    A_ij = Array{Complex{Float64},2}(undef, (4, 4))
    for i = 1:4
        for j = 1:4
            A_ij[i, j] = v[:, i]' * op * v[:, j]
            L_ij[i, j] = v[:, i]' * op * v[:, j] * v[:, i] * v[:, j]'
        end
    end
    dρ = []
    hls = []
    w_ab = Array{Float64,2}(undef, (4, 4))
    for i = 1:4
        for j = 1:4
            w_ab[i, j] = w[j] - w[i]
            if i != j
                T = L_ij[i, j]' * L_ij[i, j]
                push!(
                    dρ,
                    γ(w[j] - w[i]) * (
                        L_ij[i, j] * ρ * L_ij[i, j]' -
                        0.5 * (T * ρ + ρ * T)
                    ),
                )
                push!(hls, T * S(w[j] - w[i]))
            end
            T = L_ij[i, i]' * L_ij[j, j]
            push!(
                dρ,
                γ(0) *
                (L_ij[i, i] * ρ * L_ij[j, j]' - 0.5 * (T * ρ + ρ * T)),
            )
            push!(hls, T * S(0))
        end
    end
    dρ = sum(dρ)
    hls = sum(hls)
    dρ = dρ - 1.0im * (hls * ρ - ρ * hls)
    dρ, A_ij
end

# calculate the ame trajectory term
function ame_trajectory_update_term(op, w, v, γ, S)
    A_ij = Array{Complex{Float64},2}(undef, (4, 4))
    for i = 1:4
        for j = 1:4
            A_ij[i, j] = v[:, i]' * op * v[:, j]
        end
    end
    cache = zeros(ComplexF64, size(v, 1), size(v, 1))
    for i = 1:4
        for j = 1:4
            cache -=
                abs2(A_ij[i, j]) *
                (0.5 * γ(w[j] - w[i]) + 1.0im * S(w[j] - w[i])) *
                v[:, j] *
                v[:, j]'
        end
    end
    cache
end

# calculate the one sided AME term
function onesided_ame_udpate(op, ρ, w, v, γ, S)
    ω_ba = transpose(w) .- w
    Γ =  0.5 * γ.(ω_ba) + 1.0im * S.(ω_ba)
    L = v * (Γ .* (v' * op * v)) * v'
    K = L * ρ * op - op * L * ρ
    K + K'
end

γ(x) = x >= 0 ? x + 1 : (1 - x) * exp(x)
S(x) = x + 0.1
#= ================ Define test functions ================= =#

H = DenseHamiltonian(
    [(s) -> 1 - s, (s) -> s],
    [-standard_driver(2), (0.1 * σz ⊗ σi + 0.5 * σz ⊗ σz)],
)
coupling = ConstantCouplings(["ZI+IZ"])
davies = OpenQuantumBase.DaviesGenerator(coupling, γ, S)
onesided = OpenQuantumBase.OneSidedAMELiouvillian(coupling, OpenQuantumBase.SingleFunctionMatrix(γ), OpenQuantumBase.SingleFunctionMatrix(S), [(1, 1)])
op = 2π * (σz ⊗ σi + σi ⊗ σz)

w, v = eigen_decomp(H, 0.5, lvl=4)
w = 2π * w
ψ = (v[:, 1] + v[:, 2] + v[:, 3]) / sqrt(3)
ρ = ψ * ψ'
u = v' * ρ * v
w_ab = transpose(w) .- w
gm = γ.(w_ab)
sm = S.(w_ab)

# expected result
dρ, A_ij = ame_update_term(op, ρ, w, v, γ, S)

du = zeros(ComplexF64, (4, 4))
OpenQuantumBase.davies_update!(du, u, A_ij, gm, sm)
@test v * du * v' ≈ dρ atol = 1e-6 rtol = 1e-6

du = zeros(ComplexF64, (4, 4))
davies(du, u, w_ab, v, 0.5)
@test v * du * v' ≈ dρ atol = 1e-6 rtol = 1e-6

onesided_dρ = onesided_ame_udpate(op, ρ, w, v, γ, S)
du = zeros(ComplexF64, (4, 4))
onesided(du, u, w_ab, v, 0.5)
@test v * du * v' ≈ onesided_dρ atol = 1e-6 rtol = 1e-6

onesided_dρ = onesided_ame_udpate(v*op*v', ρ, w, v, γ, S)
du = zeros(ComplexF64, (4, 4))
onesided(du, u, w_ab, 0.5)
@test  du ≈ v' * onesided_dρ * v atol = 1e-6 rtol = 1e-6

cache = zeros(ComplexF64, (4, 4))
exp_effective_H = ame_trajectory_update_term(op, w, v, γ, S)
update_cache!(cache, davies, w_ab, v, 0.5)
@test v * cache * v' ≈ exp_effective_H atol = 1e-6 rtol = 1e-6

# test for DiffEqLiouvillian
ame_op = DiffEqLiouvillian(H, [davies], [], 4)
p = ODEParams(H, 2.0, (tf, t) -> t / tf)
exp_effective_H = ame_trajectory_update_term(op, w, v, γ, S) - 1.0im * H(0.5)
cache = zeros(ComplexF64, 4, 4)
update_cache!(cache, ame_op, p, 1)
@test cache ≈ exp_effective_H atol = 1e-6 rtol = 1e-6

hmat = H(0.5)
expected_drho = dρ - 1.0im * (hmat * ρ - ρ * hmat)
du = zeros(ComplexF64, (4, 4))
ame_op(du, ρ, p, 1)
@test du ≈ expected_drho atol = 1e-6 rtol = 1e-6

Random.seed!(1234)
jump_op = OpenQuantumBase.lindblad_jump(ame_op, ψ, p, 1)
exp_res = Complex{Float64}[6.672340269678421 + 0.0im 0.37611761184098264 + 0.0im 0.30757886113480604 + 0.0im -7.009918137704409 + 0.0im; 8.329733374366892 + 0.0im 0.46954431240210126 + 0.0im 0.38398070261603034 + 0.0im -8.751164764267996 + 0.0im; 9.465353222476072 + 0.0im 0.5335588272449766 + 0.0im 0.4363300501381911 + 0.0im -9.944239734825716 + 0.0im; 7.213269209357397 + 0.0im 0.4066095970732551 + 0.0im 0.33251438607759143 + 0.0im -7.578214632218711 + 0.0im]
@test size(jump_op) == (4, 4)
@test jump_op ≈ exp_res

# Sparse Hamiltonian test
Hd = standard_driver(4; sp=true)
Hp = q_translate("-0.9ZZII+IZZI-0.9IIZZ"; sp=true)
H = SparseHamiltonian([(s) -> 1 - s, (s) -> s], [Hd, Hp])
op = 2π * q_translate("ZIII+IZII+IIZI+IIIZ")
coupling = ConstantCouplings(["ZIII+IZII+IIZI+IIIZ"])
davies = OpenQuantumBase.DaviesGenerator(coupling, γ, S)
w, v = eigen_decomp(H, 0.5, lvl=4)
w = 2π * real(w)

ψ = (v[:, 1] + v[:, 2] + v[:, 3]) / sqrt(3)
ρ = ψ * ψ'
u = v' * ρ * v
dρ, = ame_update_term(op, ρ, w, v, γ, S)
exp_effective_H =
    ame_trajectory_update_term(op, w, v, γ, S) - 1.0im * v * Diagonal(w) * v'

ame_op = DiffEqLiouvillian(H, [davies], [], 4)
du = zeros(ComplexF64, (16, 16))
ame_op(du, ρ, p, 1)
hmat = H(0.5)
@test isapprox(
    du,
    dρ - 1.0im * (hmat * ρ - ρ * hmat),
    atol=1e-6,
    rtol=1e-6,
)

cache = zeros(ComplexF64, 16, 16)
update_cache!(cache, ame_op, p, 1)
@test cache ≈ exp_effective_H atol = 1e-6 rtol = 1e-6

# test for adiabatc frame Hamiltonian
H = AdiabaticFrameHamiltonian((s) -> [0, s, 1 - s, 1], nothing)
hmat =  H(2.0, 0.4)
w = diag(hmat)
v = collect(Diagonal(ones(4)))
coupling = CustomCouplings([(s) -> s * (σx ⊗ σi + σi ⊗ σx) + (1 - s) * (σz ⊗ σi + σi ⊗ σz)])
davies = OpenQuantumBase.DaviesGenerator(coupling, γ, S)

ψ = (v[:, 1] + v[:, 2] + v[:, 3]) / sqrt(3)
ρ = ψ * ψ'

dρ, = ame_update_term(coupling(0.4)[1], ρ, w, v, γ, S)
p = ODEParams(H, 2.0, (tf, t) -> t / tf)
ame_op = DiffEqLiouvillian(H, [davies], [], 4)
du = zeros(ComplexF64, (4, 4))
ame_op(du, ρ, p, 0.4 * 2)

@test isapprox(
    du,
    dρ - 1.0im * (hmat * ρ - ρ * hmat),
    atol=1e-6,
    rtol=1e-6,
)

# test suite for CorrelatedDaviesGenerator
coupling = ConstantCouplings([σ₊, σ₋], unit=:ħ)
γfun(w) = w >= 0 ? 1.0 : exp(-0.5)
cbath = CorrelatedBath(((1, 2), (2, 1)), spectrum=[(w) -> 0 γfun; γfun (w) -> 0])
D = OpenQuantumBase.davies_from_interactions(InteractionSet(Interaction(coupling, cbath)), 1:10, false, nothing)[1]
@test typeof(D) <: OpenQuantumBase.CorrelatedDaviesGenerator
du = zeros(ComplexF64, 2, 2)
ρ = [0.5 0;0 0.5]
ω = [1, 2]
ω =  ω' .- ω
D(du, ρ, ω, 0.5)
@test du ≈ [(1 - exp(-0.5)) * 0.5 0; 0 -(1 - exp(-0.5)) * 0.5]


text:
function get_chi_square(obj::OptimumImpl)
    return jcall(obj, "getChiSquare", jdouble, ())
end

function get_cost(obj::OptimumImpl)
    return jcall(obj, "getCost", jdouble, ())
end

function get_covariances(obj::OptimumImpl, arg0::jdouble)
    return jcall(obj, "getCovariances", RealMatrix, (jdouble,), arg0)
end

function get_evaluations(obj::OptimumImpl)
    return jcall(obj, "getEvaluations", jint, ())
end

function get_iterations(obj::OptimumImpl)
    return jcall(obj, "getIterations", jint, ())
end

function get_jacobian(obj::OptimumImpl)
    return jcall(obj, "getJacobian", RealMatrix, ())
end

function get_point(obj::OptimumImpl)
    return jcall(obj, "getPoint", RealVector, ())
end

function get_reduced_chi_square(obj::OptimumImpl, arg0::jint)
    return jcall(obj, "getReducedChiSquare", jdouble, (jint,), arg0)
end

function get_residuals(obj::OptimumImpl)
    return jcall(obj, "getResiduals", RealVector, ())
end

function get_rms(obj::OptimumImpl)
    return jcall(obj, "getRMS", jdouble, ())
end

function get_sigma(obj::OptimumImpl, arg0::jdouble)
    return jcall(obj, "getSigma", RealVector, (jdouble,), arg0)
end


text:
struct MultipatchBSplineBasis{T} <: Basis where {T}
    B::BSplineBasis #Defined on interval [0,1]
    interval::AbstractInterval #Domain
    el::Int #Number of elements in the domain
end

function DirichletMultipatchBSpline(interval::AbstractInterval, k::Int, N::Int, el::Int)
    B = DirichletBSpline(0..1, k, N)

    return MultipatchBSplineBasis{dirichlet}(B, interval, el)
end

function derivative(basis::MultipatchBSplineBasis{T}) where {T}

    return MultipatchBSplineBasis{T}(basis.B', basis.interval, basis.el)
end

function adjoint(basis::MultipatchBSplineBasis{T}) where {T}

    return derivative(basis)
end

function L2_prod(f::Function, basis::MultipatchBSplineBasis{dirichlet})
    a = basis.interval.left
    b = basis.interval.right
    Δx = (b-a)/basis.el
    n = length(basis.B.B)
    uh = zeros(n, basis.el)

    for i in 1:basis.el
        g = x -> f( a + Δx * x  )
        uh[:, i] = 1/Δx * L2_prod(g, basis.B)
        a += Δx
    end

    # Which representation is better?
    uh2 = zeros((n-1)*basis.el+1)
    for i in 1:basis.el
        uh2[(i-1)*(n-1)+1:i*(n-1)+1] += uh[1:end, i]
    end

    return uh2
end

function L2_prod_d(f::Function, basis::MultipatchBSplineBasis{dirichlet})
    a = basis.interval.left
    b = basis.interval.right
    Δx = (b-a)/basis.el
    n = length(basis.B.B)
    uh = zeros(n, basis.el)

    for i in 1:basis.el
        g = x -> f( a + Δx * x  )
        uh[:, i] = 1/Δx^2 * L2_prod(g, basis.B)
        a += Δx
    end

    # Which representation is better?
    uh2 = zeros(n*basis.el)
    for i in 1:basis.el
        uh2[(i-1)*n+1:i*n] += uh[1:end, i]
    end

    return uh2
end

function mass_matrix(basis::MultipatchBSplineBasis{dirichlet})

    a = basis.interval.left
    b = basis.interval.right
    Δx = (b-a)/basis.el
    n = length(basis.B.B)

    M = mass_matrix(basis.B)

    MM = zeros((n-1)*basis.el+1, (n-1)*basis.el+1)
    for i in 1:basis.el
        MM[(i-1)*(n-1)+1:i*(n-1)+1, (i-1)*(n-1)+1:i*(n-1)+1] += M
    end

    return 1/Δx .* MM
end

function mass_matrix_d(basis::MultipatchBSplineBasis{dirichlet})

    a = basis.interval.left
    b = basis.interval.right
    Δx = (b-a)/basis.el
    n = length(basis.B.B)

    M = mass_matrix(basis.B)

    MM = zeros(n*basis.el, n*basis.el)
    for i in 1:basis.el
        MM[(i-1)*n+1:i*n, (i-1)*n+1:i*n] += M
    end

    return 1/Δx^3 .* MM
end

function stiffness_matrix(basis::MultipatchBSplineBasis{dirichlet})
    a = basis.interval.left
    b = basis.interval.right
    Δx = (b-a)/basis.el
    n = length(basis.B.B)

    M = stiffness_matrix(basis.B)
    MM = zeros(n*basis.el+1, n*basis.el+1)

    for i in 1:basis.el
        MM[(i-1)*n+1:i*(n)+1, (i-1)*n+1:i*n+1] += M
    end

    return 1/Δx^3 .* MM
end

function L2_proj(f::Function, basis::MultipatchBSplineBasis{T}) where {T}
    M = mass_matrix(basis)

    return M\L2_prod(f, basis)
end

function L2_proj_d(f::Function, basis::MultipatchBSplineBasis{T}) where {T}
    M = mass_matrix_d(basis)

    return M\L2_prod_d(f, basis)
end

function der_mat(basis::MultipatchBSplineBasis{dirichlet})
    el = basis.el
    n = length(basis'.B.B)
    n2 = length(basis.B.B)
    DD = zeros(n*basis.el, (n2-1)*basis.el+1)
    D = der_mat(n2)[2:end, :]

    # n-1 x n matrix
    for i in 1:el
        DD[1 + (i-1)*n:i*n,1 + (i-1)*(n2-1):i*(n2-1)+1] += D
    end

    DD[end, end] = 1
    return DD
end

function der_mat_d(basis::MultipatchBSplineBasis{dirichlet})
    el = basis.el
    n = length(basis'.B.B)
    n2 = length(basis.B.B)
    DD = zeros(n*basis.el, n2*basis.el)
    D = der_mat(n2)[2:end, :]

    for i in 1:el
        DD[1 + (i-1)*n:i*n,1 + (i-1)*n2:i*n2] += D
    end

    return DD
end

function disc_grad(basis::MultipatchBSplineBasis{T}) where {T}
    return der_mat(basis)
end

function col_mat(basis::MultipatchBSplineBasis{dirichlet}, m::Int)
    a = basis.interval.left
    b = basis.interval.right
    Δx = (b-a)/basis.el
    n = length(basis.B.B)
    k = basis.B.k
    N = basis.B.N
    el = basis.el

    CM = zeros((m-1) * el+1, (n-1)*el+1)
    cm = col_mat(basis.B, m)[2]

    for i in 1:el
        CM[(i-1)*(m-1)+1:i*(m-1)+1, (i-1)*(n-1)+1:i*(n-1)+1] = cm
    end

    xx = range(a, b, length = (m-1)*el+1)
    return xx, CM
end

function col_mat_d(basis::MultipatchBSplineBasis{dirichlet}, m::Int)
    a = basis.interval.left
    b = basis.interval.right
    Δx = (b-a)/basis.el
    n = length(basis.B.B)
    k = basis.B.k
    N = basis.B.N
    el = basis.el

    CM = zeros(m*el, n*el)
    cm = col_mat(basis.B, m)[2]

    for i in 1:el
        CM[(i-1)*m+1:i*m, (i-1)*(n)+1:i*n]  = 1/Δx * cm
    end

    xx = range(0, Δx, length = m)
    x = zeros(m*el)

    for i in 1:el
        x[(i-1)*m + 1:i*m] += a .+ Δx * (i-1) .+ xx
    end

    return x, CM
end

function conforming_averaging_projection(c::Vector, basis::MultipatchBSplineBasis{dirichlet})
    n = length(basis.B.B)
    el = basis.el

    new_c = zeros(el *(n-1)+1)
    for i in 1:el
        new_c[1+(i-1)*(n-1):i*(n-1)+1] += c[1+(i-1)*(n-1)+(i-1):i*(n-1)+(i-1)+1]
    end

    new_c[(n-1)+1:(n-1):(el-1)*(n-1)+1] *= 1/2

    return new_c
end

function conforming_averaging_projection_matrix(basis::MultipatchBSplineBasis{dirichlet})
    n = length(basis.B.B)
    el = basis.el

    new_c = zeros(el *(n-1)+1)

    M = zeros(el*(n-1)+1, el*n)
    for i in 1:el
        M[1+(i-1)*(n-1):i*(n-1)+1, 1+(i-1)*(n-1)+(i-1):i*(n-1)+1+(i-1)] += I
    end

    M[(n-1)+1:(n-1):(el-1)*(n-1)+1, :] *= 1/2

    return M
end

function weighted_mass_matrix(basis::MultipatchBSplineBasis{dirichlet}, ω::Function)
    a = basis.interval.left
    b = basis.interval.right
    Δx = (b-a)/basis.el
    n = length(basis.B.B)

    M = zeros((n-1)*basis.el+1, (n-1)*basis.el+1)

    for i in 1:basis.el
        g = x -> ω(a + Δx * x)
        M[(i-1)*(n-1)+1:i*(n-1)+1, (i-1)*(n-1)+1:i*(n-1)+1] = 1/Δx * weighted_mass_matrix(basis.B, g)
        a += Δx
    end

    return M
end

function weighted_mass_matrix_d(basis::MultipatchBSplineBasis{dirichlet}, ω::Function)
    a = basis.interval.left
    b = basis.interval.right
    Δx = (b-a)/basis.el
    n = length(basis.B.B)

    M = zeros(n*basis.el, n*basis.el)

    for i in 1:basis.el
        g = x -> ω(a + Δx * x)
        M[(i-1)*n+1:i*n, (i-1)*n+1:i*n] += 1/Δx^2 * weighted_mass_matrix(basis.B, g)
        a += Δx
    end

    return M
end

#-----

function PeriodicMultipatchBSpline(interval::AbstractInterval, k::Int, N::Int, el::Int)
    B = DirichletBSpline(0..1, k, N)

    return MultipatchBSplineBasis{periodic}(B, interval, el)
end

function ToDirichlet(basis::MultipatchBSplineBasis{periodic})

    return MultipatchBSplineBasis{dirichlet}(basis.B, basis.interval, basis.el)
end

function L2_prod(f::Function, basis::MultipatchBSplineBasis{periodic})

    fh = L2_prod(f, ToDirichlet(basis))
    fh_p = fh[1:end-1]
    fh_p[1] += fh[end]

    return fh_p
end

function L2_prod_d(f::Function, basis::MultipatchBSplineBasis{periodic})
    fh = L2_prod_d(f, ToDirichlet(basis))

    return fh
end

function mass_matrix(basis::MultipatchBSplineBasis{periodic})

    M = mass_matrix(ToDirichlet(basis))
    M_p = zero(M[1:end-1, 1:end-1])
    M_p += M[1:end-1, 1:end-1]
    M_p[1, 1:end] += M[end, 1:end-1]
    M_p[1:end, 1] += M[1:end-1, end]
    M_p[1,1] += M[end,end]

    return M_p
end

function mass_matrix_d(basis::MultipatchBSplineBasis{periodic})

    M = mass_matrix_d(ToDirichlet(basis))

    return M
end

function stiffness_matrix(basis::MultipatchBSplineBasis{periodic})

    M = stiffness_matrix(ToDirichlet(basis))
    M_p = zero(M[1:end-1, 1:end-1])
    M_p += M[1:end-1, 1:end-1]
    M_p[1, 1:end] += M[end, 1:end-1]
    M_p[1:end, 1] += M[1:end-1, end]
    M_p[1,1] += M[end,end]

    return M_p

end

function col_mat(basis::MultipatchBSplineBasis{periodic}, m::Int)
    xx, cm = col_mat(ToDirichlet(basis), m)

    cm_p = zero(cm[:, 1:end-1])
    cm_p += cm[:, 1:end-1]
    cm_p[:, 1] += cm[:, end]

    return xx, cm_p
end

function col_mat_d(basis::MultipatchBSplineBasis{periodic}, m::Int)
    xx, cm = col_mat_d(ToDirichlet(basis), m)

    return xx, cm
end

#this is not fully correct I think, as I don't know how the derivative of the periodic function should behave...
function der_mat(basis::MultipatchBSplineBasis{periodic})
    DD = SplineFEEC.der_mat(SplineFEEC.ToDirichlet(basis))
    D_p = zero(DD[ :, 1:end-1])
    D_p[ 1:end-1, 1:end] += DD[ 1:end-1, 1:end-1]
    D_p[end, end] = -1
    D_p[end, 1] = 1

    return D_p
end

function der_mat_d(basis::MultipatchBSplineBasis{periodic})
    D = der_mat_d(ToDirichlet(basis))

    return D
end

function conforming_averaging_projection_matrix(basis::MultipatchBSplineBasis{periodic})
    P = conforming_averaging_projection_matrix(ToDirichlet(basis))
    P = P[1:end-1, :]
    P[1, end] = 1
    P[1, :] *= 1/2

    return P
end

function weighted_mass_matrix(basis::MultipatchBSplineBasis{periodic}, ω::Function)
    M = weighted_mass_matrix(ToDirichlet(basis))
    M_p = zero(M[1:end-1, 1:end-1])
    M_p += M[1:end-1, 1:end-1]
    M_p[1, 1:end] += M[end, 1:end-1]
    M_p[1:end, 1] += M[1:end-1, end]
    M_p[1,1] += M[end,end]

    return M_p
end

function weighted_mass_matrix_d(basis::MultipatchBSplineBasis{periodic}, ω::Function)
     M = weighted_mass_matrix_d(ToDirichlet(basis))

    return M
end

#-----

struct MultipatchBSplineTensorProductBasis{N, T} <: TensorProductBasis{N}
    B::SArray{Tuple{N}, MultipatchBSplineBasis{T}}
end

function DirichletMultipatchTensorProductBSpline(domain::ProductDomain, k::NTuple{d,Int}, N::NTuple{d,Int}, el::NTuple{d,Int}) where {d}
    basis =[]
    for i in 1:d
        push!(basis, DirichletMultipatchBSpline(domain.domains[i], k[i], N[i], el[i]))
    end

    return MultipatchBSplineTensorProductBasis{d, dirichlet}(SVector(basis...))
end

function PeriodicMultipatchTensorProductBSpline(domain::ProductDomain, k::NTuple{d,Int}, N::NTuple{d,Int}, el::NTuple{d,Int}) where {d}
    basis =[]
    for i in 1:d
        push!(basis, PeriodicMultipatchBSpline(domain.domains[i], k[i], N[i], el[i]))
    end

    return MultipatchBSplineTensorProductBasis{d, periodic}(SVector(basis...))
end

function ToDirichlet(basis::MultipatchBSplineTensorProductBasis{d, periodic}) where {d}

    return MultipatchBSplineTensorProductBasis{d, dirichlet}(ToDirichlet.(basis.B))
end

getindex(B::MultipatchBSplineTensorProductBasis{N, T}, j::Integer) where {N, T} = B.B[j]
getindex(B::MultipatchBSplineTensorProductBasis{N, T}, j::Colon) where {N, T} = B.B

#tensorproduct
function tensorproduct(a::MultipatchBSplineBasis{T}, b::MultipatchBSplineBasis{T}) where T
    Bs = SVector(a, b)
    MultipatchBSplineTensorProductBasis{2, T}(SVector(a, b))
end

function tensorproduct(a::MultipatchBSplineTensorProductBasis{d, T}, b::MultipatchBSplineBasis{T}) where {d, T}
    B = SVector(a.B..., b)
    MultipatchBSplineTensorProductBasis{d+1, T}(B)
end

function tensorproduct(a::MultipatchBSplineBasis{T}, b::MultipatchBSplineTensorProductBasis{d, T}) where {d, T}
    B = SVector(a, b.B...)
    MultipatchBSplineTensorProductBasis{d+1, T}(B)
end

function tensorproduct(a::MultipatchBSplineTensorProductBasis{d1, T}, b::MultipatchBSplineTensorProductBasis{d2, T}) where {d1, d2, T}
    B = SVector(a.B..., b.B...)
    MultipatchBSplineTensorProductBasis{d1+d2, T}(B)
end

⊗(a::Union{MultipatchBSplineTensorProductBasis, MultipatchBSplineBasis}, b::Union{MultipatchBSplineTensorProductBasis, MultipatchBSplineBasis}) = tensorproduct(a, b)

function derivative(basis::MultipatchBSplineTensorProductBasis{d, T}) where {d, T}
    b = derivative(basis.B[1])
    for i in 2:d
        bn = derivative(basis.B[i])
        b = b ⊗ bn
    end

    return b
end

function adjoint(basis::MultipatchBSplineTensorProductBasis)
    return derivative(basis)
end

#Methods
function L2_prod(f::Function, basis::MultipatchBSplineTensorProductBasis{2, dirichlet})
    a1 = basis[1].interval.left
    b1 = basis[1].interval.right
    el1 = basis[1].el
    Δx1 = (b1-a1)/basis[1].el
    n1 = length(basis[1].B.B)

    a2 = basis[2].interval.left
    b2 = basis[2].interval.right
    el2 = basis[2].el
    Δx2 = (b2-a2)/el2
    n2 = length(basis[2].B.B)

    uh = zeros(((n1-1)*el1+1)*((n2-1)*el2+1))

    for i in 1:el1
        for j in 1:el2
            g = (x, y) -> f(a1 + Δx1 * x, a2 + Δx2*y)
            val = 1/(Δx1 * Δx2) * L2_prod(g, basis[1].B ⊗ basis[2].B)

            for k in 1:n1
                uh[1+(k-1)*n2+(k-1)*(el2-1)*(n2-1)+(j-1)*el2*(n2-1)+(i-1)*(el2-1)*(n2-1)*(n1-1)+(i-1)*n2*(n1-1)-(j-1)*(n2-1)*(el2-1):k*n2+(k-1)*(el2-1)*(n2-1)+(j-1)*el2*(n2-1)+(i-1)*(el2-1)*(n1-1)*(n2-1)-(j-1)*(n2-1)*(el2-1)+(i-1)*n2*(n1-1)] += val[1+(k-1)*n2:k*n2]
            end

            a2 += Δx2
        end

        a1 += Δx1
    end
    return uh
end

function L2_prod_mixed_x(f::Function, basis::MultipatchBSplineTensorProductBasis{2, dirichlet})
    a1 = basis[1].interval.left
    b1 = basis[1].interval.right
    el1 = basis[1].el
    Δx1 = (b1-a1)/basis[1].el
    n1 = length(basis[1].B.B)

    a2 = basis[2].interval.left
    b2 = basis[2].interval.right
    el2 = basis[2].el
    Δx2 = (b2-a2)/el2
    n2 = length(basis[2].B.B)

    uh = zeros((n1*el1)*((n2-1)*el2+1))

    for i in 1:el1
        for j in 1:el2
            g = (x, y) -> f(a1 + Δx1 * x, a2 + Δx2*y)
            val = 1/(Δx1^2 * Δx2) * L2_prod(g, basis[1].B ⊗ basis[2].B)

            for k in 1:n1
                uh[1+(k-1)*n2+(k-1)*(el2-1)*(n2-1)+(j-1)*el2*(n2-1)+(i-1)*(el2-1)*(n2-1)*n1+(i-1)*n2*n1-(j-1)*(n2-1)*(el2-1):k*n2+(k-1)*(el2-1)*(n2-1)+(j-1)*el2*(n2-1)+(i-1)*(el2-1)*n1*(n2-1)-(j-1)*(n2-1)*(el2-1)+(i-1)*n2*n1] += val[1+(k-1)*n2:k*n2]
            end

            a2 += Δx2
        end

        a1 += Δx1
    end
    return uh
end

function L2_prod_mixed_y(f::Function, basis::MultipatchBSplineTensorProductBasis{2, dirichlet})
    a1 = basis[1].interval.left
    b1 = basis[1].interval.right
    el1 = basis[1].el
    Δx1 = (b1-a1)/basis[1].el
    n1 = length(basis[1].B.B)

    a2 = basis[2].interval.left
    b2 = basis[2].interval.right
    el2 = basis[2].el
    Δx2 = (b2-a2)/el2
    n2 = length(basis[2].B.B)

    uh = zeros(((n1-1)*el1+1)*n2*el2)

    for i in 1:el1
        for j in 1:el2
            g = (x, y) -> f(a1 + Δx1 * x, a2 + Δx2*y)
            val = 1/(Δx1 * Δx2^2) * L2_prod(g, basis[1].B ⊗ basis[2].B)

            for k in 1:n1
                uh[1+(k-1)*n2+(k-1)*(el2-1)*n2+(j-1)*el2*n2+(i-1)*(el2-1)*n2*(n1-1)+(i-1)*n2*(n1-1)-(j-1)*n2*(el2-1):k*n2+(k-1)*(el2-1)*n2+(j-1)*el2*n2+(i-1)*(el2-1)*(n1-1)*n2-(j-1)*n2*(el2-1)+(i-1)*n2*(n1-1)] += val[1+(k-1)*n2:k*n2]
            end

            a2 += Δx2
        end

        a1 += Δx1
    end
    return uh
end

function L2_prod_d(f::Function, basis::MultipatchBSplineTensorProductBasis{2, dirichlet})
    a1 = basis[1].interval.left
    b1 = basis[1].interval.right
    el1 = basis[1].el
    Δx1 = (b1-a1)/basis[1].el
    n1 = length(basis[1].B.B)

    a2 = basis[2].interval.left
    b2 = basis[2].interval.right
    el2 = basis[2].el
    Δx2 = (b2-a2)/el2
    n2 = length(basis[2].B.B)

    uh = zeros(n1*el1*n2*el2)

    for i in 1:el1
        for j in 1:el2
            g = (x, y) -> f(a1 + Δx1 * x, a2 + Δx2*y)
            val = 1/(Δx1 * Δx2)^2 * L2_prod(g, basis[1].B ⊗ basis[2].B)

            for k in 1:n1
                uh[1+(k-1)*n2+(k-1)*(el2-1)*n2+(j-1)*el2*n2+(i-1)*(el2-1)*n2*n1+(i-1)*n2*n1-(j-1)*n2*(el2-1):k*n2+(k-1)*(el2-1)*n2+(j-1)*el2*n2+(i-1)*(el2-1)*n1*n2-(j-1)*n2*(el2-1)+(i-1)*n2*n1] += val[1+(k-1)*n2:k*n2]
            end

            a2 += Δx2
        end

        a1 += Δx1
    end

    return uh
end

function col_mat(basis::MultipatchBSplineTensorProductBasis{2}, mm::Int)
    X,Y = basis.B

    x, M1 = col_mat(X, mm)
    y, M2 = col_mat(Y, mm)

    xy = [(a,b) for a in x for b in y]

    return xy, M1 ⊗ M2
end

function col_mat_mixed_x(basis::MultipatchBSplineTensorProductBasis{2}, mm::Int)
    X,Y = basis.B

    x, M1 = col_mat_d(X, mm)
    y, M2 = col_mat(Y, mm)

    xy = [(a,b) for a in x for b in y]

    return xy, M1 ⊗ M2
end

function col_mat_mixed_y(basis::MultipatchBSplineTensorProductBasis{2}, mm::Int)
    X,Y = basis.B

    x, M1 = col_mat(X, mm)
    y, M2 = col_mat_d(Y, mm)

    xy = [(a,b) for a in x for b in y]

    return xy, M1 ⊗ M2
end
function col_mat_d(basis::MultipatchBSplineTensorProductBasis{2}, mm::Int)
    X,Y = basis.B

    x, M1 = col_mat_d(X, mm)
    y, M2 = col_mat_d(Y, mm)

    xy = [(a,b) for a in x for b in y]

    return xy, M1 ⊗ M2
end

#n-D Calculations
function mass_matrix(basis::MultipatchBSplineTensorProductBasis{d, T}) where {d, T}
    matrices = []
    for i in 1:d
        push!(matrices, mass_matrix(basis[i]))
    end

    return ⊗(matrices...)
end

function mass_matrix_d(basis::MultipatchBSplineTensorProductBasis{d, T}) where {d, T}
    matrices = []
    for i in 1:d
        push!(matrices, mass_matrix_d(basis[i]))
    end

    return ⊗(matrices...)
end

function stiffness_matrix(basis::MultipatchBSplineTensorProductBasis{d, T}) where {d, T}
    matrices = []
    for i in 1:d
        push!(matrices,  stiffness_matrix(basis[i]))
    end

    return ⊗(matrices...)
end

function disc_diff(T, basis::MultipatchBSplineTensorProductBasis{2, dirichlet})
    if T == grad
        N1 = (length(basis[1].B.B)-1)*basis[1].el+1
        N2 = (length(basis[2].B.B)-1)*basis[2].el+1

        D1 = der_mat(basis[1]) ⊗ circ_id(N2)
        D2 = circ_id(N1) ⊗ der_mat(basis[2])
    elseif T == curl
        N1 = length(basis[1]'.B.B)*basis[1].el
        N2 = length(basis[2]'.B.B)*basis[2].el

        D1 = der_mat(basis[1]) ⊗ circ_id(N2)
        D2 = circ_id(N1) ⊗ der_mat(basis[2])
    end

    return discrete_differential{T, 2}(SVector(D1, D2))
end

function disc_diff_d(T, basis::MultipatchBSplineTensorProductBasis{2, ty}) where ty
    if T == grad
        N1 = length(basis[1].B.B)*basis[1].el
        N2 = length(basis[2].B.B)*basis[2].el

        D1 = der_mat_d(basis[1]) ⊗ circ_id(N2)
        D2 = circ_id(N1) ⊗ der_mat_d(basis[2])
    elseif T == curl
        N1 = length(basis[1]'.B.B)*basis[1].el
        N2 = length(basis[2]'.B.B)*basis[2].el

        D1 = der_mat_d(basis[1]) ⊗ circ_id(N2)
        D2 = circ_id(N1) ⊗ der_mat_d(basis[2])
    end

    return discrete_differential{T, 2}(SVector(D1, D2))
end

function L2_proj(f::Function, basis::MultipatchBSplineTensorProductBasis{N, T}) where {N, T}
    M = mass_matrix(basis)
    fh = M\L2_prod(f, basis)
    return SVector{length(fh)}(fh)
end

function L2_proj_d(f::Function, basis::MultipatchBSplineTensorProductBasis{N, T}) where {N, T}
    M = mass_matrix_d(basis)
    fh = M\L2_prod_d(f, basis)
    return SVector{length(fh)}(fh)
end

function L2_proj(f::Function, basis::SArray{Tuple{2}, MultipatchBSplineTensorProductBasis{2, T}}) where {T}
    fh = []
    M = mass_matrix.(basis)

    for i in 1:2
        push!(fh, M[i]\L2_prod((x,y) -> f(x,y)[i], basis[i]))
    end

    return SVector(fh...)
end

function L2_proj_d(f::Function, basis::SArray{Tuple{2}, MultipatchBSplineTensorProductBasis{2, T}}) where {T}
    fh = []
    M = mass_matrix_d.(basis)

    for i in 1:2
        push!(fh, M[i]\L2_prod_d((x,y) -> f(x,y)[i], basis[i]))
    end

    return SVector(fh...)
end

function col_mat_proj(basis::MultipatchBSplineTensorProductBasis{2, T}) where {T}
    X,Y = basis.B
    N = (length(X.B.B), length(Y.B.B))

    x, M1 = col_mat(X, N[1])
    y, M2 = col_mat(Y, N[2])

    xy = [(a,b) for a in x for b in y]

    return xy, M1 ⊗ M2
end

function col_mat_proj_mixed_x(basis::MultipatchBSplineTensorProductBasis{2, T}) where {T}
    X,Y = basis.B
    N = (length(X.B.B), length(Y.B.B))

    x, M1 = col_mat_d(X, N[1])
    y, M2 = col_mat(Y, N[2])

    xy = [(a,b) for a in x for b in y]

    return xy, M1 ⊗ M2
end

function col_mat_proj_mixed_y(basis::MultipatchBSplineTensorProductBasis{2, T}) where {T}
    X,Y = basis.B
    N = (length(X.B.B), length(Y.B.B))

    x, M1 = col_mat(X, N[1])
    y, M2 = col_mat_d(Y, N[2])

    xy = [(a,b) for a in x for b in y]

    return xy, M1 ⊗ M2
end

function col_mat_proj_d(basis::MultipatchBSplineTensorProductBasis{2, T}) where {T}
    X,Y = basis.B
    N = (length(X.B.B), length(Y.B.B))

    x, M1 = col_mat_d(X, N[1])
    y, M2 = col_mat_d(Y, N[2])

    xy = [(a,b) for a in x for b in y]

    return xy, M1 ⊗ M2
end

function conforming_averaging_projection_matrix(basis::MultipatchBSplineTensorProductBasis{2, T}) where {T}
    X,Y = basis.B
    Mx = conforming_averaging_projection_matrix(X)
    My = conforming_averaging_projection_matrix(Y)

    return Mx ⊗ My
end

function conforming_averaging_projection_matrix_x(basis::MultipatchBSplineTensorProductBasis{2, T}) where {T}
    X,Y = basis.B

    Mx = conforming_averaging_projection_matrix(X)
    My = circ_id(Y.el*(length(Y.B.B)-1))

    return Mx ⊗ My
end

function conforming_averaging_projection_matrix_y(basis::MultipatchBSplineTensorProductBasis{2, T}) where {T}
    X,Y = basis.B

    Mx = circ_id(X.el*(length(X.B.B)-1))
    My = conforming_averaging_projection_matrix(Y)

    return Mx ⊗ My
end

#todo weighted mass matrices

function L2_prod(f::Function, basis::MultipatchBSplineTensorProductBasis{2, periodic})
    a1 = basis[1].interval.left
    b1 = basis[1].interval.right
    el1 = basis[1].el
    Δx1 = (b1-a1)/basis[1].el
    n1 = length(basis[1].B.B)

    a2 = basis[2].interval.left
    b2 = basis[2].interval.right
    el2 = basis[2].el
    Δx2 = (b2-a2)/el2
    n2 = length(basis[2].B.B)

    uh = zeros(((n1-1)*el1+1)*((n2-1)*el2+1))

    for i in 1:el1
        for j in 1:el2
            g = (x, y) -> f(a1 + Δx1 * x, a2 + Δx2*y)
            val = 1/(Δx1 * Δx2) * L2_prod(g, basis[1].B ⊗ basis[2].B)

            for k in 1:n1
                uh[1+(k-1)*n2+(k-1)*(el2-1)*(n2-1)+(j-1)*el2*(n2-1)+(i-1)*(el2-1)*(n2-1)*(n1-1)+(i-1)*n2*(n1-1)-(j-1)*(n2-1)*(el2-1):k*n2+(k-1)*(el2-1)*(n2-1)+(j-1)*el2*(n2-1)+(i-1)*(el2-1)*(n1-1)*(n2-1)-(j-1)*(n2-1)*(el2-1)+(i-1)*n2*(n1-1)] += val[1+(k-1)*n2:k*n2]
            end

            a2 += Δx2
        end

        a1 += Δx1
    end

    uh_new = zeros(((n1-1)*el1)*((n2-1)*el2))
    for i in 1:(n1-1)*el1
            uh_new[1+(n2-1)*el2*(i-1):i*el2*(n2-1)] += uh[1+(n2-1)*el2*(i-1)+(i-1):i*el2*(n2-1)+(i-1)]
            uh_new[1+(n2-1)*el2*(i-1)] += uh[1+(n2-1)*el2*i+(i-1)]
    end

    uh_new[1:el2*(n2-1)] += uh[end-el2*(n2-1):end-1]
    uh_new[1] += uh[end]

    return uh_new
end

function L2_prod_mixed_x(f::Function, basis::MultipatchBSplineTensorProductBasis{2, periodic})
    a1 = basis[1].interval.left
    b1 = basis[1].interval.right
    el1 = basis[1].el
    Δx1 = (b1-a1)/basis[1].el
    n1 = length(basis[1].B.B)

    a2 = basis[2].interval.left
    b2 = basis[2].interval.right
    el2 = basis[2].el
    Δx2 = (b2-a2)/el2
    n2 = length(basis[2].B.B)

    uh = zeros((n1*el1)*((n2-1)*el2+1))

    for i in 1:el1
        for j in 1:el2
            g = (x, y) -> f(a1 + Δx1 * x, a2 + Δx2*y)
            val = 1/(Δx1^2 * Δx2) * L2_prod(g, basis[1].B ⊗ basis[2].B)

            for k in 1:n1
                uh[1+(k-1)*n2+(k-1)*(el2-1)*(n2-1)+(j-1)*el2*(n2-1)+(i-1)*(el2-1)*(n2-1)*n1+(i-1)*n2*n1-(j-1)*(n2-1)*(el2-1):k*n2+(k-1)*(el2-1)*(n2-1)+(j-1)*el2*(n2-1)+(i-1)*(el2-1)*n1*(n2-1)-(j-1)*(n2-1)*(el2-1)+(i-1)*n2*n1] += val[1+(k-1)*n2:k*n2]
            end

            a2 += Δx2
        end

        a1 += Δx1
    end

    uh_new = zeros((n1*el1)*((n2-1)*el2))
    for i in 1:n1*el1
            uh_new[1+(n2-1)*el2*(i-1):i*el2*(n2-1)] += uh[1+(n2-1)*el2*(i-1)+(i-1):i*el2*(n2-1)+(i-1)]
            uh_new[1+(n2-1)*el2*(i-1)] += uh[1+(n2-1)*el2*i+(i-1)]
    end

    return uh_new

end

function L2_prod_mixed_y(f::Function, basis::MultipatchBSplineTensorProductBasis{2, periodic})
    a1 = basis[1].interval.left
    b1 = basis[1].interval.right
    el1 = basis[1].el
    Δx1 = (b1-a1)/basis[1].el
    n1 = length(basis[1].B.B)

    a2 = basis[2].interval.left
    b2 = basis[2].interval.right
    el2 = basis[2].el
    Δx2 = (b2-a2)/el2
    n2 = length(basis[2].B.B)

    uh = zeros(((n1-1)*el1+1)*n2*el2)

    for i in 1:el1
        for j in 1:el2
            g = (x, y) -> f(a1 + Δx1 * x, a2 + Δx2*y)
            val = 1/(Δx1 * Δx2^2) * L2_prod(g, basis[1].B ⊗ basis[2].B)

            for k in 1:n1
                uh[1+(k-1)*n2+(k-1)*(el2-1)*n2+(j-1)*el2*n2+(i-1)*(el2-1)*n2*(n1-1)+(i-1)*n2*(n1-1)-(j-1)*n2*(el2-1):k*n2+(k-1)*(el2-1)*n2+(j-1)*el2*n2+(i-1)*(el2-1)*(n1-1)*n2-(j-1)*n2*(el2-1)+(i-1)*n2*(n1-1)] += val[1+(k-1)*n2:k*n2]
            end

            a2 += Δx2
        end

        a1 += Δx1
    end

    uh_new = zeros(((n1-1)*el1)*n2*el2)
    for i in 1:(n1-1)*el1
            uh_new[1+(n2)*el2*(i-1):i*el2*(n2)] += uh[1+(n2)*el2*(i-1):i*el2*(n2)]
    end

    uh_new[1:(n2)*el2] += uh[end-(n2)*el2+1:end]

    return uh_new
end

function L2_prod_d(f::Function, basis::MultipatchBSplineTensorProductBasis{2, periodic})
    a1 = basis[1].interval.left
    b1 = basis[1].interval.right
    el1 = basis[1].el
    Δx1 = (b1-a1)/basis[1].el
    n1 = length(basis[1].B.B)

    a2 = basis[2].interval.left
    b2 = basis[2].interval.right
    el2 = basis[2].el
    Δx2 = (b2-a2)/el2
    n2 = length(basis[2].B.B)

    uh = zeros(n1*el1*n2*el2)

    for i in 1:el1
        for j in 1:el2
            g = (x, y) -> f(a1 + Δx1 * x, a2 + Δx2*y)
            val = 1/(Δx1 * Δx2)^2 * L2_prod(g, basis[1].B ⊗ basis[2].B)

            for k in 1:n1
                uh[1+(k-1)*n2+(k-1)*(el2-1)*n2+(j-1)*el2*n2+(i-1)*(el2-1)*n2*n1+(i-1)*n2*n1-(j-1)*n2*(el2-1):k*n2+(k-1)*(el2-1)*n2+(j-1)*el2*n2+(i-1)*(el2-1)*n1*n2-(j-1)*n2*(el2-1)+(i-1)*n2*n1] += val[1+(k-1)*n2:k*n2]
            end

            a2 += Δx2
        end

        a1 += Δx1
    end

    return uh
end

function disc_diff(T, basis::MultipatchBSplineTensorProductBasis{2, periodic})
    if T == grad
        N1 = (length(basis[1].B.B)-1)*basis[1].el
        N2 = (length(basis[2].B.B)-1)*basis[2].el

        D1 = der_mat(basis[1]) ⊗ circ_id(N2)
        D2 = circ_id(N1) ⊗ der_mat(basis[2])
    elseif T == curl
        N1 = length(basis[1]'.B.B)*basis[1].el
        N2 = length(basis[2]'.B.B)*basis[2].el

        D1 = der_mat(basis[1]) ⊗ circ_id(N2)
        D2 = circ_id(N1) ⊗ der_mat(basis[2])
    end

    return discrete_differential{T, 2}(SVector(D1, D2))
end

#(X', Y) and (X, Y')
function mixed_mass_matrix_x(basis1::MultipatchBSplineTensorProductBasis{2, periodic}, basis2::MultipatchBSplineTensorProductBasis{2, periodic}) #where T
    X1, Y1 = basis1.B
    X2, Y2 = basis2.B

    # x basis functions
    a1 = X1.interval.left
    b1 = X1.interval.right
    n1 = length(X1.B.B)
    el1 = X1.el
    Δx1 = (b1-a1)/el1

    a2 = X2.interval.left
    b2 = X2.interval.right
    n2 = length(X2.B.B)
    el2 = X2.el
    Δx2 = (b2-a2)/el2

    Mx = mass_matrix(X1.B, X2.B)
    MMx = zeros(n1*el1, (n2-1)*el2+1)

    for i in 1:el1
        for j in 1:el2
            MMx[1+(i-1)*n1:i*n1,1+(j-1)*(n2-1):j*(n2-1)+1] += Mx
        end
    end

    #periodic in x in 2nd basis function
    M_px = zero(MMx[1:end, 1:end-1])
    M_px += MMx[1:end, 1:end-1]
    M_px[1:end, 1] += MMx[:, end]

    # y basis functions
    a1 = Y1.interval.left
    b1 = Y1.interval.right
    n1 = length(Y1.B.B)
    el1 = Y1.el
    Δy1 = (b1-a1)/el1

    a2 = Y2.interval.left
    b2 = Y2.interval.right
    n2 = length(Y2.B.B)
    el2 = Y2.el
    Δy2 = (b2-a2)/el2

    My = mass_matrix(Y1.B, Y2.B)
    MMy = zeros((n1-1)*el1+1, n2*el2)

    for i in 1:el1
        for j in 1:el2
            MMy[1+(i-1)*(n1-1):i*(n1-1)+1,1+(j-1)*n2:j*n2] += My
        end
    end

    #periodic in x in 2nd basis function
    M_py = zero(MMy[1:end-1, 1:end])
    M_py += MMy[1:end-1, 1:end]
    M_py[1, 1:end] += MMy[end, :]

    return 1/Δx1 * 1/Δx2 * M_px ⊗ M_py
end

#(X', Y) and (X, Y')
function mixed_mass_matrix_x_d(basis1::MultipatchBSplineTensorProductBasis{2, periodic}, basis2::MultipatchBSplineTensorProductBasis{2, periodic}) #where T
    X1, Y1 = basis1.B
    X2, Y2 = basis2.B

    # x basis functions
    a1 = X1.interval.left
    b1 = X1.interval.right
    n1 = length(X1.B.B)
    el1 = X1.el
    Δx1 = (b1-a1)/el1

    a2 = X2.interval.left
    b2 = X2.interval.right
    n2 = length(X2.B.B)
    el2 = X2.el
    Δx2 = (b2-a2)/el2

    Mx = mass_matrix(X1.B, X2.B)
    MMx = zeros(n1*el1, n2*el2)

    for i in 1:el1
        for j in 1:el2
            MMx[1+(i-1)*n1:i*n1,1+(j-1)*n2:j*n2] += Mx
        end
    end

    # y basis functions
    a1 = Y1.interval.left
    b1 = Y1.interval.right
    n1 = length(Y1.B.B)
    el1 = Y1.el
    Δy1 = (b1-a1)/el1

    a2 = Y2.interval.left
    b2 = Y2.interval.right
    n2 = length(Y2.B.B)
    el2 = Y2.el
    Δy2 = (b2-a2)/el2

    My = mass_matrix(Y1.B, Y2.B)
    MMy = zeros(n1*el1, n2*el2)

    for i in 1:el1
        for j in 1:el2
            MMy[1+(i-1)*n1:i*n1,1+(j-1)*n2:j*n2] += My
        end
    end

    return 1/Δx1^2 * 1/Δx2^2 * MMx ⊗ MMy
end

function L2_proj_mixed_V1(f::Function, basis::SArray{Tuple{2}, MultipatchBSplineTensorProductBasis{2, T}}) where {T}
    X, Y = basis[1].B
    M1x = mass_matrix_d(X)
    M1y = mass_matrix(Y)
    M1_1 = M1x ⊗ M1y

    fh_1 = M1_1\L2_prod_mixed_x((x,y) -> f(x,y)[1], basis[1])

    X, Y = basis[2].B
    M1x = mass_matrix(X)
    M1y = mass_matrix_d(Y)
    M1_2 = M1x ⊗ M1y

    fh_2 = M1_2\L2_prod_mixed_y((x,y) -> f(x,y)[2], basis[2])

    return SVector(fh_1, fh_2)
end