text:
#=
module Gaugefields_2D_nowing_module
    using LinearAlgebra
    import ..AbstractGaugefields_module:AbstractGaugefields,Shifted_Gaugefields,shift_U,
                        Adjoint_Gaugefields,set_wing_U!,Abstractfields,construct_staple!,clear_U!,
                        calculate_Plaquette,substitute_U!,mul_skiplastindex!,partial_tr,add_U!,
                        Traceless_antihermitian!
    import Base
    import ..Gaugefields_2D_module:Gaugefields_2D
=#

    struct Gaugefields_2D_nowing{NC} <: Gaugefields_2D{NC}
        U::Array{ComplexF64,4}
        NX::Int64
        NT::Int64
        NDW::Int64
        NV::Int64
        NC::Int64
        mpi::Bool
        verbose_print::Verbose_print
        Ushifted::Array{ComplexF64,4}

        function Gaugefields_2D_nowing(NC::T,NX::T,NT::T;verbose_level = 2) where T<: Integer
            NV = NX*NT
            NDW = 0
            U = zeros(ComplexF64,NC,NC,NX+2NDW,NT+2NDW)
            mpi = false
            verbose_print = Verbose_print(verbose_level )
            Ushifted = zero(U)
            #U = Array{Array{ComplexF64,4}}(undef,2)
            #for μ=1:2
            #    U[μ] = zeros(ComplexF64,NC,NC,NX+2NDW,NY+2NDW,NZ+2NDW,NT+2NDW)
            #end
            return new{NC}(U,NX,NT,NDW,NV,NC,mpi,verbose_print,Ushifted)
        end
    end

    function write_to_numpyarray(U::T,filename) where T <: Gaugefields_2D_nowing
        data = Dict{String,Any}()
        data["U"] = U.U
        data["NX"] = U.NX
        data["NT"] = U.NT
        data["NV"] = U.NV
        data["NDW"] = U.NDW
        data["NC"] = U.NC

        npzwrite(filename, data)
    end

    function Base.setindex!(x::Gaugefields_2D_nowing,v,i1,i2,i3,i6)
        @inbounds x.U[i1,i2,i3,i6] = v
    end

    @inline function Base.getindex(x::Gaugefields_2D_nowing,i1,i2,i3,i6)
        @inbounds return x.U[i1,i2,i3  ,i6  ]
    end

    function Base.setindex!(x::Gaugefields_2D_nowing,v,i1,i2,ii)
        ix,it =  get_latticeindex(ii,x.NX,x.NT)
        @inbounds x.U[i1,i2,ix ,it ] = v
    end

    @inline function Base.getindex(x::Gaugefields_2D_nowing,i1,i2,ii)
        ix,it =  get_latticeindex(ii,x.NX,x.NT)
        @inbounds return x.U[i1,i2,ix  ,it  ]
    end

    function substitute_U!(a::Array{T1,1},b::Array{T2,1}) where {T1 <: Gaugefields_2D_nowing,T2 <: Gaugefields_2D_nowing}
        for μ=1:2
            substitute_U!(a[μ],b[μ])
        end
    end

    function substitute_U!(a::Array{T1,1},b::Array{T2,1},iseven) where {T1 <: Gaugefields_2D_nowing,T2 <: Gaugefields_2D_nowing}
        for μ=1:2
            substitute_U!(a[μ],b[μ],iseven)
        end
    end

    function Base.similar(U::T) where T <: Gaugefields_2D_nowing
        Uout = Gaugefields_2D_nowing(U.NC,U.NX,U.NT;verbose_level = U.verbose_print.level)
        #identityGaugefields_2D_nowing(U.NC,U.NX,U.NY,U.NZ,U.NT,U.NDW)
        return Uout
    end

    function Base.similar(U::Array{T,1}) where T <: Gaugefields_2D_nowing
        Uout = Array{T,1}(undef,2)
        for μ=1:2
            Uout[μ] = similar(U[μ])
        end
        return Uout
    end

    function substitute_U!(a::T,b::T) where T <: Gaugefields_2D_nowing
        for i=1:length(a.U)
            a.U[i] = b.U[i]
        end
        return
    end

    function substitute_U!(a::Gaugefields_2D_nowing{NC},b::T2) where {NC, T2 <: Abstractfields}
        NT = a.NT
        #NZ = a.NZ
        #NY = a.NY
        NX = a.NX
        for it=1:NT
            ##for iz=1:NZ
                ##for iy=1:NY
                    for ix=1:NX
                        for k2=1:NC
                            for k1=1:NC
                                @inbounds a[k1,k2,ix,it] = b[k1,k2,ix,it]
                            end
                        end
                    end
                #end
            #end
        end
        set_wing_U!(a)

    end

    function substitute_U!(a::Gaugefields_2D_nowing{NC},b::T2,iseven::Bool) where {NC, T2 <: Abstractfields}
        NT = a.NT
        #NZ = a.NZ
        #NY = a.NY
        NX = a.NX
        for it=1:NT
            ##for iz=1:NZ
            #    #for iy=1:NY
                    for ix=1:NX
                        evenodd = ifelse( (ix+it) % 2 ==0, true,false)
                        if evenodd == iseven
                            for k2=1:NC
                                for k1=1:NC
                                    @inbounds a[k1,k2,ix,it] = b[k1,k2,ix,it]
                                end
                            end
                        end
                    end
            #    end
            #end
        end
        set_wing_U!(a)

    end

    function m3complv!(a::T) where T <: Gaugefields_2D_nowing
        aa = zeros(Float64,18)
        NX = a.NX
        #NY = a.NY
        #NZ = a.NZ
        NT = a.NT

        @inbounds for it=1:NT
            ##for iz=1:NZ
            #    #for iy=1:NY
                    @simd for ix=1:NX

                        aa[ 1] = real( a[1,1,ix,it])
                        aa[ 2] = imag(a[1,1,ix,it])
                        aa[ 3] = real( a[1,2,ix,it])
                        aa[ 4] = imag(a[1,2,ix,it])
                        aa[ 5] = real( a[1,3,ix,it])
                        aa[ 6] = imag(a[1,3,ix,it])
                        aa[ 7] = real( a[2,1,ix,it])
                        aa[ 8] = imag(a[2,1,ix,it])
                        aa[ 9] = real( a[2,2,ix,it])
                        aa[10] = imag(a[2,2,ix,it])
                        aa[11] = real( a[2,3,ix,it])
                        aa[12] = imag(a[2,3,ix,it])

                        aa[13] = aa[ 3]*aa[11] - aa[ 4]*aa[12] -
                                    aa[ 5]*aa[ 9] + aa[ 6]*aa[10]
                        aa[14] = aa[ 5]*aa[10] + aa[ 6]*aa[ 9] -
                                    aa[ 3]*aa[12] - aa[ 4]*aa[11]
                        aa[15] = aa[ 5]*aa[ 7] - aa[ 6]*aa[ 8] -
                                    aa[ 1]*aa[11] + aa[ 2]*aa[12]
                        aa[16] = aa[ 1]*aa[12] + aa[ 2]*aa[11] -
                                    aa[ 5]*aa[ 8] - aa[ 6]*aa[ 7]
                        aa[17] = aa[ 1]*aa[ 9] - aa[ 2]*aa[10] -
                                    aa[ 3]*aa[ 7] + aa[ 4]*aa[ 8]
                        aa[18] = aa[ 3]*aa[ 8] + aa[ 4]*aa[ 7] -
                                    aa[ 1]*aa[10] - aa[ 2]*aa[ 9]

                        a[3,1,ix,it] = aa[13]+im*aa[14]
                        a[3,2,ix,it] = aa[15]+im*aa[16]
                        a[3,3,ix,it] = aa[17]+im*aa[18]

                        #println(a[:,:,ix,it]'*a[:,:,ix,it] )
                    end
            #    end
            #end
        end
    end

    function randomGaugefields_2D_nowing(NC,NX,NT;verbose_level = 2,randomnumber="Random")
        U = Gaugefields_2D_nowing(NC,NX,NT,verbose_level = verbose_level)
        if randomnumber== "Random"
            rng = MersenneTwister()
            #randomfunc() = rand()
        elseif randomnumber== "Reproducible"
            rng = StableRNG(123)
            #randomfunc() = rand(rng,Float64)
        else
            error("randomnumber should be \"Random\" or \"Reproducible\". Now randomnumber = $randomnumber")
        end

        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        for j=1:NC
                            @simd for i=1:NC
                                U[i,j,ix,it] = rand(rng)-0.5 + im*(rand(rng)-0.5)
                            end
                        end
                    end
                #end
            #end
        end
        normalize_U!(U)
        set_wing_U!(U)

        return U
    end

    function identityGaugefields_2D_nowing(NC,NX,NT;verbose_level = 2)
        U = Gaugefields_2D_nowing(NC,NX,NT,verbose_level = verbose_level)

        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        @simd for ic=1:NC
                            U[ic,ic,ix,it] = 1
                        end
                    end
                #end
            #end
        end
        set_wing_U!(U)
        return U
    end

    function Oneinstanton_2D_nowing(NC,NX,NT;verbose_level = 2)
        @assert NC == 2 "NC should be 2"
        u = Gaugefields_2D_nowing(NC,NX,NT,verbose_level = verbose_level)
        U = Array{typeof(u),1}(undef,2)
        U[1] = u
        for μ = 2:2
            U[μ] = similar(u)
        end
        L = (NX,NT)
        NV = prod(L)

        R = div(NX,2) # instanton radius

        println("# Starting from a instanton backgorund with radius R=$R ")
        inst_cent = [L[1]/2+0.5, L[2]/2+0.5]
        #inst_cent = [L[1]/2+0.5, L[2]/2+0.5, L[3]/2+0.5, L[4]/2+0.5]
        #eps = 1/10000000
        s1 = [
        0.0 1.0
        1.0 0.0]
        s2 = [
        0.0 -im*1.0
        im*1.0 0.0]
        s3 = [
        1.0 0.0
        0.0 -1.0]
        En = [
        1.0 0.0
        0.0 1.0]
        ss=[ im*s1,  im*s2,  im*s3, En]
        sd=[-im*s1, -im*s2, -im*s3, En]
        nn = 0

        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        nn += 0
                        nv=[ix-1, it-1]-inst_cent
                        n2=nv⋅nv
                        for mu = 1:2
                            tau=[
                                0 0
                                0 0]
                            for nu = 1:2
                                smunu = sd[mu]*ss[nu]-sd[nu]*ss[mu]
                                tau+=smunu*nv[nu]
                            end
                            sq = sqrt(n2+R^2)
                            tau = exp(im*tau*(1/2)*(1/(n2))*(im*R^2/(n2+R^2)) ) #1b
                            for j=1:2
                                for i=1:2
                                    U[mu][i,j,ix,it] = tau[i,j]
                                end
                            end

                        end
                    end
                #end
            #end
        end
        set_wing_U!(U)

        return U

    end

    struct Shifted_Gaugefields_2D_nowing{NC} <: Shifted_Gaugefields{NC,2}
        parent::Gaugefields_2D_nowing{NC}
        #parent::T
        shift::NTuple{2,Int8}
        NX::Int64
        #NY::Int64
        #NZ::Int64
        NT::Int64
        NDW::Int64

        #function Shifted_Gaugefields(U::T,shift,Dim) where {T <: AbstractGaugefields}
        function Shifted_Gaugefields_2D_nowing(U::Gaugefields_2D_nowing{NC},shift) where NC
            shifted_U!(U,shift)
            return new{NC}(U,shift,U.NX,U.NT,U.NDW)
        end
    end

    function shifted_U!(U::Gaugefields_2D_nowing{NC},shift) where NC
        NT = U.NT
        NX = U.NX
        for it=1:NT
            it_shifted  = it + shift[2]
            it_shifted += ifelse(it_shifted > NT,-NT,0)
            it_shifted += ifelse(it_shifted < 1,NT,0)

                    for ix=1:NX
                        ix_shifted  = ix + shift[1]
                        ix_shifted += ifelse(ix_shifted > NX,-NX,0)
                        ix_shifted += ifelse(ix_shifted < 1,NX,0)
                        for k2=1:NC
                            for k1=1:NC
                                U.Ushifted[k1,k2,ix,it] = U[k1,k2,ix_shifted,it_shifted]
                            end
                        end
                    end

        end
    end

        #lattice shift
    function shift_U(U::Gaugefields_2D_nowing{NC},ν::T) where {T <: Integer,NC}
            if ν == 1
                shift = (1,0)
            elseif ν == 2
                shift = (0,1)
            #elseif ν == 3
            #    shift = (0,0,1,0)
            #elseif ν == 4
            #    shift = (0,0,0,1)
            elseif ν == -1
                    shift = (-1,0)
            elseif ν == -2
                    shift = (0,-1)
            #elseif ν == -3
            #        shift = (0,0,-1,0)
            #elseif ν == -4
            #        shift = (0,0,0,-1)
            end

            return Shifted_Gaugefields_2D_nowing(U,shift)
    end

    function shift_U(U::TU,shift::NTuple{Dim,T}) where {Dim,T <: Integer,TU <: Gaugefields_2D_nowing}
        return Shifted_Gaugefields_2D_nowing(U,shift)
    end

    @inline function Base.getindex(U::Shifted_Gaugefields_2D_nowing{NC},i1,i2,i3,i6) where NC
        #function Base.getindex(U::Shifted_Gaugefields{T,4},i1,i2,i3,i4,i5,i6) where T <: Gaugefields_2D_nowing
        @inbounds return U.parent.Ushifted[i1,i2,i3 ,i6]
    end

    function Base.getindex(u::Staggered_Gaugefields{T,μ},i1,i2,i3,i6) where {T <: Gaugefields_2D_nowing,μ}
        NT = u.parent.NT
        #NZ = u.parent.NZ
        #NY = u.parent.NY
        NX = u.parent.NX

        t = i6-1
        t += ifelse(t<0,NT,0)
        t += ifelse(t ≥ NT,-NT,0)
        #boundary_factor_t = ifelse(t == NT -1,BoundaryCondition[4],1)

        x = i3-1
        x += ifelse(x<0,NX,0)
        x += ifelse(x ≥ NX,-NX,0)
        #boundary_factor_x = ifelse(x == NX -1,BoundaryCondition[1],1)
        if μ ==1
            η = 1
        elseif μ ==2
            #η = (-1.0)^(x)
            η = ifelse(x%2 == 0,1,-1)
        else
            error("η should be positive but η = $η")
        end

        @inbounds return η*u.parent[i1,i2,i3,i6]
    end

    function map_U!(U::Gaugefields_2D_nowing{NC},f!::Function,V::Gaugefields_2D_nowing{NC},iseven::Bool) where {NC}
        NT = U.NT
        #NZ = U.NZ
        #NY = U.NY
        NX = U.NX
        A = zeros(ComplexF64,NC,NC)
        B = zeros(ComplexF64,NC,NC)
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        evenodd = ifelse( (ix+it) % 2 ==0, true,false)
                        if evenodd == iseven
                            for k2=1:NC
                                for k1=1:NC
                                    A[k1,k2] = V[k1,k2,ix,it]
                                    B[k1,k2] = U[k1,k2,ix,it]
                                end
                            end
                            f!(B,A)
                            for k2=1:NC
                                for k1=1:NC
                                    U[k1,k2,ix,it] = B[k1,k2]
                                end
                            end
                        end
                    #end
                #end
            end
        end
        set_wing_U!(U)
    end

    function map_U_sequential!(U::Gaugefields_2D_nowing{NC},f!::Function,Uin) where {NC}
        NT = U.NT
        #NZ = U.NZ
        #NY = U.NY
        NX = U.NX
        #A = zeros(ComplexF64,NC,NC)
        B = zeros(ComplexF64,NC,NC)
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX

                        for k2=1:NC
                            for k1=1:NC
                                B[k1,k2] = U[k1,k2,ix,it]
                            end
                        end
                        f!(B,Uin,ix,it)
                        for k2=1:NC
                            for k1=1:NC
                                U[k1,k2,ix,it] = B[k1,k2]
                            end
                        end

                    end
                #end
            #end
        end
        set_wing_U!(U)
    end

    function clear_U!(Uμ::Gaugefields_2D_nowing{NC}) where NC
        NT = Uμ.NT
        #NZ = Uμ.NZ
        #NY = Uμ.NY
        NX = Uμ.NX
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        for k2=1:NC
                            for k1=1:NC
                                @inbounds Uμ[k1,k2,ix,it] = 0
                            end
                        end
                    end
                #end
            #end
        end
        set_wing_U!(Uμ)

    end

    function clear_U!(Uμ::Gaugefields_2D_nowing{NC},iseven::Bool) where NC
        NT = Uμ.NT
        #NZ = Uμ.NZ
        #NY = Uμ.NY
        NX = Uμ.NX
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        evenodd = ifelse( (ix+it) % 2 ==0, true,false)
                        if evenodd == iseven
                            for k2=1:NC
                                for k1=1:NC
                                    @inbounds Uμ[k1,k2,ix,it] = 0
                                end
                            end
                        end
                    end
                #end
            #end
        end
        set_wing_U!(Uμ)

    end

    function unit_U!(Uμ::Gaugefields_2D_nowing{NC}) where NC
        NT = Uμ.NT
        #NZ = Uμ.NZ
        #NY = Uμ.NY
        NX = Uμ.NX
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX

                        for k2=1:NC
                            for k1=1:NC
                                @inbounds Uμ[k1,k2,ix,it] = ifelse(k1 == k2,1,0)
                            end
                        end
                    end
                #end
            #end
        end
        set_wing_U!(Uμ)

    end

    """
    M = (U*δ_prev) star (dexp(Q)/dQ)
    Λ = TA(M)
    """
    function construct_Λmatrix_forSTOUT!(Λ,δ_current::Gaugefields_2D_nowing{NC},Q,u::Gaugefields_2D_nowing{NC}) where NC
        NT = u.NT
        #NY = u.NY
        #NZ = u.NZ
        NX = u.NX
        Qn = zeros(ComplexF64,NC,NC)
        Un = zero(Qn)
        Mn = zero(Qn)
        Λn = zero(Qn)
        δn_current = zero(Qn)
        temp1 = similar(Qn)
        temp2 = similar(Qn)
        temp3 = similar(Qn)

        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX

                        for jc=1:NC
                            for ic=1:NC
                                Un[ic,jc] = u[ic,jc,ix,it]
                                Qn[ic,jc] = Q[ic,jc,ix,it]#*im
                                #if (ix,it) == (1,1,1,1)
                                #    println("Qij $ic $jc ",Qn[ic,jc],"\t ",Q[ic,jc,ix,it])
                                #end
                                δn_current[ic,jc] = δ_current[ic,jc,ix,it]
                            end
                        end

                        #=
                        if (ix,it) == (1,1,1,1)
                            println("Qn = ",Qn[1:3,1:3])
                        end
                        =#

                        calc_Mmatrix!(Mn,δn_current,Qn,Un,u,[temp1,temp2,temp3])
                        #=
                        if (ix,it) == (1,1,1,1)
                            println(" Un ",  Un[1,1])
                            println("δn_current ", δn_current[1,1])
                            #println("Qn[1,1] = ",Qn[1,1])
                            println("M[1,1] = ",Mn[1,1])

                            #println("Qn = ",Qn[1:3,1:3])
                        end
                        =#
                        calc_Λmatrix!(Λn,Mn,NC)

                        for jc=1:NC
                            for ic=1:NC
                                Λ[ic,jc,ix,it] = Λn[ic,jc]
                            end
                        end

                    end
                #end
            #end
        end
        set_wing_U!(Λ)
    end

    function set_wing_U!(u::Array{Gaugefields_2D_nowing{NC},1}) where NC
        return
    end

    function set_wing_U!(u::Gaugefields_2D_nowing{NC}) where NC
        return
    end

    function exptU!(uout::T,t::N,u::Gaugefields_2D_nowing{NC},temps::Array{T,1}) where {N <: Number, T <: Gaugefields_2D_nowing, NC} #uout = exp(t*u)
        @assert NC != 3 && NC != 2 "This function is for NC != 2,3"

        NT = u.NT
        #NZ = u.NZ
        #NY = u.NY
        NX = u.NX
        V0 = zeros(ComplexF64,NC,NC)
        V1 = zeros(ComplexF64,NC,NC)
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        for k2=1:NC
                            for k1=1:NC
                                @inbounds V0[k1,k2] = im*t*u[k1,k2,ix,it]
                            end
                        end
                        V1[:,:] = exp(V0)
                        for k2=1:NC
                            for k1=1:NC
                                @inbounds uout[k1,k2,ix,it] = V1[k1,k2]
                            end
                        end

                    end
                #end
            #end
        end
        #error("exptU! is not implemented in type $(typeof(u)) ")
    end

    const tinyvalue =1e-100
    const pi23 = 2pi/3
    const fac13 = 1/3

    # #=
    function exptU!(uout::T,t::N,v::Gaugefields_2D_nowing{3},temps::Array{T,1}) where {N <: Number, T <: Gaugefields_2D_nowing} #uout = exp(t*u)
        ww = temps[1]
        w = temps[2]

        NT = v.NT
        #NZ = v.NZ
        #NY = v.NY
        NX = v.NX
        #t = 1

        @inbounds for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        v11 = v[1,1,ix,it]
                        v22 = v[2,2,ix,it]
                        v33 = v[3,3,ix,it]

                        tri = fac13*(imag(v11)+imag(v22)+imag(v33))

                        #=
                        vout[1,1,ix,it] = (imag(v11)-tri)*im
                        vout[2,2,ix,it] = (imag(v22)-tri)*im
                        vout[3,3,ix,it] = (imag(v33)-tri)*im
                        =#
                        y11 = (imag(v11)-tri)*im
                        y22 = (imag(v22)-tri)*im
                        y33 = (imag(v33)-tri)*im

                        v12 = v[1,2,ix,it]
                        v13 = v[1,3,ix,it]
                        v21 = v[2,1,ix,it]
                        v23 = v[2,3,ix,it]
                        v31 = v[3,1,ix,it]
                        v32 = v[3,2,ix,it]

                        x12 = v12 - conj(v21)
                        x13 = v13 - conj(v31)
                        x23 = v23 - conj(v32)

                        x21 = - conj(x12)
                        x31 = - conj(x13)
                        x32 = - conj(x23)

                        y12 = 0.5  * x12
                        y13 = 0.5  * x13
                        y21 = 0.5  * x21
                        y23 = 0.5  * x23
                        y31 = 0.5  * x31
                        y32 =  0.5  * x32

                        c1_0 = ( imag(y12) + imag(y21) )
                        c2_0 = ( real(y12) - real(y21) )
                        c3_0 = ( imag(y11) - imag(y22) )
                        c4_0= ( imag(y13) + imag(y31) )
                        c5_0 = ( real(y13) - real(y31) )

                        c6_0 = ( imag(y23) + imag(y32) )
                        c7_0 = ( real(y23) - real(y32) )
                        c8_0 = sr3i *
                                ( imag(y11) + imag(y22) -
                                        2*imag(y33) )

                        c1 = t*c1_0 * 0.5
                        c2 = t*c2_0 * 0.5
                        c3 = t*c3_0 * 0.5
                        c4 = t*c4_0 * 0.5
                        c5 = t*c5_0 * 0.5
                        c6 = t*c6_0 * 0.5
                        c7 = t*c7_0 * 0.5
                        c8 = t*c8_0 * 0.5
                        csum = c1+c2+c3+c4+c5+c6+c7+c8
                        if csum == 0
                            w[1,1,ix,it]=   1
                            w[1,2,ix,it]=   0
                            w[1,3,ix,it]=   0
                            w[2,1,ix,it]=  0
                            w[2,2,ix,it]=   1
                            w[2,3,ix,it]=   0
                            w[3,1,ix,it]=   0
                            w[3,2,ix,it]=   0
                            w[3,3,ix,it]=   1

                            ww[1,1,ix,it]=   1
                            ww[1,2,ix,it]=   0
                            ww[1,3,ix,it]=   0
                            ww[2,1,ix,it]=   0
                            ww[2,2,ix,it]=   1
                            ww[2,3,ix,it]=   0
                            ww[3,1,ix,it]=   0
                            ww[3,2,ix,it]=   0
                            ww[3,3,ix,it]=   1
                            continue
                        end

                        #x[1,1,icum] =  c3+sr3i*c8 +im*(  0.0 )
                        v1 = c3+sr3i*c8
                        v2 = 0.0
                        #x[1,2,icum] =  c1         +im*( -c2   )
                        v3 = c1
                        v4 = -c2
                        #x[1,3,icum] =  c4         +im*(-c5   )
                        v5 = c4
                        v6 = -c5

                        #x[2,1,icum] =  c1         +im*(  c2   )
                        v7 = c1
                        v8 = c2

                        #x[2,2,icum] =  -c3+sr3i*c8+im*(  0.0 )
                        v9 =-c3+sr3i*c8
                        v10 = 0.0

                        #x[2,3,icum] =  c6         +im*( -c7   )
                        v11 = c6
                        v12 = -c7

                        #x[3,1,icum] =  c4         +im*(  c5   )
                        v13 = c4
                        v14 = c5

                        #x[3,2,icum] =  c6         +im*(  c7   )
                        v15 = c6
                        v16 = c7
                        #x[3,3,icum] =  -sr3i2*c8  +im*(  0.0 )
                        v17 = -sr3i2*c8
                        v18 = 0.0

                #c find eigenvalues of v
                        trv3 = (v1 + v9 + v17) / 3.0
                        cofac = v1 *  v9  - v3^2  - v4^2 +
                                v1 * v17 -  v5^2  - v6^2 +
                                v9 * v17 - v11^2 - v12^2
                        det = v1 * v9 * v17 -
                                v1 * (v11^2 + v12^2) -
                                v9 * (v5^2 + v6^2) -
                                v17 * (v3^2 + v4^2) +
                                (v5 * (v3 * v11 - v4 * v12) +
                                v6 * (v3 * v12 + v4 * v11)) * 2.0
                        p3 = cofac / 3.0 - trv3^2
                        q = trv3 * cofac - det - 2.0 * trv3^3
                        x = sqrt(-4.0 * p3) + tinyvalue
                        arg = q / (x * p3)

                        arg = min(1, max(-1, arg))
                        theta = acos(arg) / 3.0
                        e1 = x * cos(theta) + trv3
                        theta = theta + pi23
                        e2 = x * cos(theta) + trv3
                #       theta = theta + pi23
                #       e3 = x * cos(theta) + trv3
                        e3 = 3.0 * trv3 - e1 - e2

                # solve for eigenvectors

                        w1 =   v5 * (v9 - e1) - v3 * v11 + v4 * v12
                        w2 = - v6 * (v9 - e1) + v4 * v11 + v3 * v12
                        w3 =   (v1 - e1) * v11 - v3 * v5 - v4 * v6
                        w4 = - (v1 - e1) * v12 - v4 * v5 + v3 * v6
                        w5 = - (v1 - e1) * (v9 - e1) +  v3^2 + v4^2
                        w6 = 0.0

                        coeff = 1.0 / sqrt(w1^2 + w2^2 + w3^2 + w4^2 + w5^2)

                        w1 = w1 * coeff
                        w2 = w2 * coeff
                        w3 = w3 * coeff
                        w4 = w4 * coeff
                        w5 = w5 * coeff

                        w7 =   v5 * (v9 - e2) - v3 * v11 + v4 * v12
                        w8 = - v6 * (v9 - e2) + v4 * v11 + v3 * v12
                        w9 = (v1 - e2) * v11 - v3 * v5 - v4 * v6
                        w10 = - (v1 - e2) * v12 - v4 * v5 + v3 * v6
                        w11 = - (v1 - e2) * (v9 - e2) +  v3^2 + v4^2
                        w12 = 0.0

                        coeff = 1.0 / sqrt(w7^2  + w8^2 + w9^2+ w10^2 + w11^2)

                        w7  = w7  * coeff
                        w8  = w8  * coeff
                        w9  = w9  * coeff
                        w10 = w10 * coeff
                        w11 = w11 * coeff

                        w13 =   v5 * (v9 - e3) - v3 * v11 + v4 * v12
                        w14 = - v6 * (v9 - e3) + v4 * v11 + v3 * v12
                        w15 =   (v1 - e3) * v11 - v3 * v5 - v4 * v6
                        w16 = - (v1 - e3) * v12 - v4 * v5 + v3 * v6
                        w17 = - (v1 - e3) * (v9 - e3) +  v3^2 + v4^2
                        w18 = 0.0

                        coeff = 1.0 / sqrt(w13^2 + w14^2 + w15^2+ w16^2 + w17^2)
                        w13 = w13 * coeff
                        w14 = w14 * coeff
                        w15 = w15 * coeff
                        w16 = w16 * coeff
                        w17 = w17 * coeff

                # construct the projection v
                        c1 = cos(e1)
                        s1 = sin(e1)
                        ww1  = w1  * c1 - w2 * s1
                        ww2  = w2  * c1 + w1 * s1
                        ww3  = w3  * c1 - w4 * s1
                        ww4  = w4  * c1 + w3 * s1
                        ww5  = w5  * c1 - w6 * s1
                        ww6  = w6  * c1 + w5 * s1

                        c2 = cos(e2)
                        s2 = sin(e2)
                        ww7  = w7  * c2 - w8 * s2
                        ww8  = w8  * c2 + w7 * s2
                        ww9  = w9  * c2 - w10 * s2
                        ww10 = w10 * c2 + w9 * s2
                        ww11 = w11 * c2 - w12 * s2
                        ww12 = w12 * c2 + w11 * s2

                        c3 = cos(e3)
                        s3 = sin(e3)
                        ww13 = w13 * c3 - w14 * s3
                        ww14 = w14 * c3 + w13 * s3
                        ww15 = w15 * c3 - w16 * s3
                        ww16 = w16 * c3 + w15 * s3
                        ww17 = w17 * c3 - w18 * s3
                        ww18 = w18 * c3 + w17 * s3

                        w[1,1,ix,it]=   w1+im*w2
                        w[1,2,ix,it]=   w3+im*w4
                        w[1,3,ix,it]=   w5+im*w6
                        w[2,1,ix,it]=   w7+im*w8
                        w[2,2,ix,it]=   w9+im*w10
                        w[2,3,ix,it]=   w11+im*w12
                        w[3,1,ix,it]=   w13+im*w14
                        w[3,2,ix,it]=   w15+im*w16
                        w[3,3,ix,it]=   w17+im*w18

                        ww[1,1,ix,it]=   ww1+im*ww2
                        ww[1,2,ix,it]=   ww3+im*ww4
                        ww[1,3,ix,it]=   ww5+im*ww6
                        ww[2,1,ix,it]=   ww7+im*ww8
                        ww[2,2,ix,it]=   ww9+im*ww10
                        ww[2,3,ix,it]=   ww11+im*ww12
                        ww[3,1,ix,it]=   ww13+im*ww14
                        ww[3,2,ix,it]=   ww15+im*ww16
                        ww[3,3,ix,it]=   ww17+im*ww18

                        #a = ww[:,:,ix,it]
                        #b = w[:,:,ix,it]
                        #println(b'*a)
                        #println(exp(im*t*v[:,:,ix,it]))
                        #error("d")
                    end
                #end
            #end
        end

        #mul!(v,w',ww)
        mul!(uout,w',ww)

        #error("exptU! is not implemented in type $(typeof(u)) ")
    end
    # =#

    """
-----------------------------------------------------c
     !!!!!   vin and vout should be different vectors

     Projectin of the etraceless antiermite part
     vout = x/2 - Tr(x)/6
     wher   x = vin - Conjg(vin)
-----------------------------------------------------c
    """

    #Q = -(1/2)*(Ω' - Ω) + (1/(2NC))*tr(Ω' - Ω)*I0_2
    #Omega' - Omega = -2i imag(Omega)
    function Traceless_antihermitian!(vout::Gaugefields_2D_nowing{3},vin::Gaugefields_2D_nowing{3})
        #error("Traceless_antihermitian! is not implemented in type $(typeof(vout)) ")
        fac13 = 1/3
        NX = vin.NX
        #NY = vin.NY
        #NZ = vin.NZ
        NT = vin.NT

        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    @simd for ix=1:NX
                        v11 = vin[1,1,ix,it]
                        v21 = vin[2,1,ix,it]
                        v31 = vin[3,1,ix,it]

                        v12 = vin[1,2,ix,it]
                        v22 = vin[2,2,ix,it]
                        v32 = vin[3,2,ix,it]

                        v13 = vin[1,3,ix,it]
                        v23 = vin[2,3,ix,it]
                        v33 = vin[3,3,ix,it]

                        tri = fac13*(imag(v11)+imag(v22)+imag(v33))

                        #=
                        vout[1,1,ix,it] = (imag(v11)-tri)*im
                        vout[2,2,ix,it] = (imag(v22)-tri)*im
                        vout[3,3,ix,it] = (imag(v33)-tri)*im
                        =#
                        y11 = (imag(v11)-tri)*im
                        y22 = (imag(v22)-tri)*im
                        y33 = (imag(v33)-tri)*im

                        x12 = v12 - conj(v21)
                        x13 = v13 - conj(v31)
                        x23 = v23 - conj(v32)

                        x21 = - conj(x12)
                        x31 = - conj(x13)
                        x32 = - conj(x23)

                        #=
                        vout[1,2,ix,it] = 0.5  * x12
                        vout[1,3,ix,it] = 0.5  * x13
                        vout[2,1,ix,it] = 0.5  * x21
                        vout[2,3,ix,it] = 0.5  * x23
                        vout[3,1,ix,it] = 0.5  * x31
                        vout[3,2,ix,it] = 0.5  * x32
                        =#
                        y12 = 0.5  * x12
                        y13 = 0.5  * x13
                        y21 = 0.5  * x21
                        y23 = 0.5  * x23
                        y31 = 0.5  * x31
                        y32 =  0.5  * x32

                        vout[1,1,ix,it] = y11
                        vout[2,1,ix,it] = y21
                        vout[3,1,ix,it] = y31

                        vout[1,2,ix,it] = y12
                        vout[2,2,ix,it] = y22
                        vout[3,2,ix,it] = y32

                        vout[1,3,ix,it] = y13
                        vout[2,3,ix,it] = y23
                        vout[3,3,ix,it] = y33

                    end
                #end
            #end
        end

    end

    function Traceless_antihermitian!(vout::Gaugefields_2D_nowing{2},vin::Gaugefields_2D_nowing{2})
        #error("Traceless_antihermitian! is not implemented in type $(typeof(vout)) ")
        fac13 = 1/3
        NX = vin.NX
        #NY = vin.NY
        #NZ = vin.NZ
        NT = vin.NT

        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    @simd for ix=1:NX

                        v11 = vin[1,1,ix,it]
                        v22 = vin[2,2,ix,it]

                        tri = (imag(v11)+imag(v22))*0.5

                        v12 = vin[1,2,ix,it]
                        #v13 = vin[1,3,ix,it]
                        v21 = vin[2,1,ix,it]

                        x12 = v12 - conj(v21)

                        x21 = - conj(x12)

                        vout[1,1,ix,it] = (imag(v11)-tri)*im
                        vout[1,2,ix,it] = 0.5  * x12
                        vout[2,1,ix,it] = 0.5  * x21
                        vout[2,2,ix,it] = (imag(v22)-tri)*im
                    end
                #end
            #end
        end

    end

    function Traceless_antihermitian!(vout::Gaugefields_2D_nowing{1},vin::Gaugefields_2D_nowing{1})
        #error("Traceless_antihermitian! is not implemented in type $(typeof(vout)) ")
        fac13 = 1/3
        NX = vin.NX
        #NY = vin.NY
        #NZ = vin.NZ
        NT = vin.NT

        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    @simd for ix=1:NX
                        v11 = vin[1,1,ix,it]
                        vout[1,1,ix,it] = imag(v11)*im

                    end
                #end
            #end
        end

    end

    function Traceless_antihermitian!(vout::Gaugefields_2D_nowing{NC},vin::Gaugefields_2D_nowing{NC}) where NC
        #NC = vout.NC
        fac1N = 1/NC
        nv = vin.NV

        NX = vin.NX
        #NY = vin.NY
        #NZ = vin.NZ
        NT = vin.NT

        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    @simd for ix=1:NX
                        tri = 0.0
                        @simd for k=1:NC
                            tri += imag(vin[k,k,ix,it])
                        end
                        tri *= fac1N
                        @simd for k=1:NC
                            vout[k,k,ix,it] = (imag(vin[k,k,ix,it])-tri)*im
                        end
                    end
                #end
            #end
        end

        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    @simd for ix=1:NX
                        for k1=1:NC
                            @simd for k2=k1+1:NC
                                vv = 0.5*(vin[k1,k2,ix,it] - conj(vin[k2,k1,ix,it]))
                                vout[k1,k2,ix,it] = vv
                                vout[k2,k1,ix,it] = -conj(vv)
                            end
                        end
                    end
                #end
            #end
        end

    end

    function LinearAlgebra.tr(a::Gaugefields_2D_nowing{NC},b::Gaugefields_2D_nowing{NC}) where NC
        NX=a.NX
        #NY=a.NY
        #NZ=a.NZ
        NT=a.NT

        s = 0
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        for k=1:NC
                            for k2=1:NC
                                s += a[k,k2,ix,it]*b[k2,k,ix,it]
                            end
                        end
                    end
                #end
            #end
        end
        #println(3*NT*NZ*NY*NX*NC)
        return s
    end

    function LinearAlgebra.tr(a::Gaugefields_2D_nowing{1},b::Gaugefields_2D_nowing{1})
        NX=a.NX
        #NY=a.NY
        #NZ=a.NZ
        NT=a.NT

        s = 0
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        a11 = a[1,1,ix,it]

                        b11 = b[1,1,ix,it]

                        s += a11*b11

                    end
                #end
            #end
        end
        #println(3*NT*NZ*NY*NX*NC)
        return s
    end

    function LinearAlgebra.tr(a::Gaugefields_2D_nowing{2},b::Gaugefields_2D_nowing{2})
        NX=a.NX
        #NY=a.NY
        #NZ=a.NZ
        NT=a.NT

        s = 0
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        a11 = a[1,1,ix,it]
                        a21 = a[2,1,ix,it]
                        a12 = a[1,2,ix,it]
                        a22 = a[2,2,ix,it]

                        b11 = b[1,1,ix,it]
                        b21 = b[2,1,ix,it]
                        b12 = b[1,2,ix,it]
                        b22 = b[2,2,ix,it]

                        s += a11*b11+a12*b21+
                            a21*b12+a22*b22

                    end
                #end
            #end
        end
        #println(3*NT*NZ*NY*NX*NC)
        return s
    end

    function LinearAlgebra.tr(a::Gaugefields_2D_nowing{3},b::Gaugefields_2D_nowing{3})
        NX=a.NX
        #NY=a.NY
        #NZ=a.NZ
        NT=a.NT

        s = 0
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        a11 = a[1,1,ix,it]
                        a21 = a[2,1,ix,it]
                        a31 = a[3,1,ix,it]
                        a12 = a[1,2,ix,it]
                        a22 = a[2,2,ix,it]
                        a32 = a[3,2,ix,it]
                        a13 = a[1,3,ix,it]
                        a23 = a[2,3,ix,it]
                        a33 = a[3,3,ix,it]
                        b11 = b[1,1,ix,it]
                        b21 = b[2,1,ix,it]
                        b31 = b[3,1,ix,it]
                        b12 = b[1,2,ix,it]
                        b22 = b[2,2,ix,it]
                        b32 = b[3,2,ix,it]
                        b13 = b[1,3,ix,it]
                        b23 = b[2,3,ix,it]
                        b33 = b[3,3,ix,it]
                        s += a11*b11+a12*b21+a13*b31 +
                            a21*b12+a22*b22+a23*b32 +
                            a31*b13+a32*b23+a33*b33
                    end
                #end
            #end
        end
        #println(3*NT*NZ*NY*NX*NC)
        return s
    end

    function LinearAlgebra.tr(a::Gaugefields_2D_nowing{NC}) where NC
        NX=a.NX
        #NY=a.NY
        #NZ=a.NZ
        NT=a.NT

        s = 0
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        @simd for k=1:NC
                            s += a[k,k,ix,it]
                            #println(a[k,k,ix,it])
                        end
                    end
                #end
            #end
        end
        #println(3*NT*NZ*NY*NX*NC)
        return s
    end

    function LinearAlgebra.tr(a::Gaugefields_2D_nowing{1})
        NX=a.NX
        #NY=a.NY
        #NZ=a.NZ
        NT=a.NT

        s = 0
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        s += a[1,1,ix,it]
                    end
                #end
            #end
        end
        #println(3*NT*NZ*NY*NX*NC)
        return s
    end

    function partial_tr(a::Gaugefields_2D_nowing{NC},μ) where NC
        NX=a.NX
        #NY=a.NY
        #NZ=a.NZ
        NT=a.NT

        if μ == 1
            s = 0
            ix = 1
            for it=1:NT
                #for iz=1:NZ
                    #for iy=1:NY
                        #for ix=1:NX
                            @simd for k=1:NC
                                s += a[k,k,ix,it]
                                #println(a[k,k,ix,it])
                            end

                        #end
                    #end
                #end
            end
        elseif μ == 2
            s = 0
            iy =1
            for it=1:NT
                #for iz=1:NZ
                    ##for iy=1:NY
                        for ix=1:NX
                            @simd for k=1:NC
                                s += a[k,k,ix,it]
                                #println(a[k,k,ix,it])
                            end
                        end
                    #end
                #end
            end
        elseif μ == 3
            s = 0
            iz = 1
            for it=1:NT
                ##for iz=1:NZ
                    #for iy=1:NY
                        for ix=1:NX
                            @simd for k=1:NC
                                s += a[k,k,ix,it]
                                #println(a[k,k,ix,it])
                            end
                        end
                    #end
                #end
            end
        else
            s = 0
            it = 1
                #for iz=1:NZ
                    #for iy=1:NY
                        for ix=1:NX
                            @simd for k=1:NC
                                s += a[k,k,ix,it]
                                #println(a[k,k,ix,it])
                            end
                        end
                    #end
                #end

        end

        #println(3*NT*NZ*NY*NX*NC)
        return s
    end

    function add_U!(c::Gaugefields_2D_nowing{NC},a::T1) where {NC,T1 <: Abstractfields}
        @inbounds for i=1:length(c.U)
            c.U[i] += a.U[i]
        end
        return

        NT = c.NT
        #NZ = c.NZ
        #NY = c.NY
        NX = c.NX
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        for k2=1:NC
                            @simd for k1=1:NC
                                c[k1,k2,ix,it] += a[k1,k2,ix,it]
                            end
                        end
                    end
                #end
            #end
        end
    end

    function add_U!(c::Gaugefields_2D_nowing{NC},a::T1,iseven::Bool) where {NC,T1 <: Abstractfields}
        NT = c.NT
        #NZ = c.NZ
        #NY = c.NY
        NX = c.NX
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        evenodd = ifelse( (ix+it) % 2 ==0, true,false)
                        if evenodd == iseven
                            for k2=1:NC
                                @simd for k1=1:NC
                                    c[k1,k2,ix,it] += a[k1,k2,ix,it]
                                end
                            end
                        end
                    end
                #end
            #end
        end
    end

    function add_U!(c::Gaugefields_2D_nowing{NC},α::N,a::T1) where {NC,T1 <: Abstractfields, N<:Number}
        #@inbounds for i=1:length(c.U)
        #    c.U[i] += α*a.U[i]
        #end
        #return

        NT = c.NT
        #NZ = c.NZ
        #NY = c.NY
        NX = c.NX
        @inbounds for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        for k2=1:NC
                            @simd for k1=1:NC
                                c[k1,k2,ix,it] += α*a[k1,k2,ix,it]
                            end
                        end
                    end
                #end
            #end
        end
        set_wing_U!(c)
    end

    function LinearAlgebra.mul!(c::Gaugefields_2D_nowing{1},a::T1,b::T2) where {T1 <: Abstractfields,T2 <: Abstractfields}
        #@assert NC != 2 && NC != 3 "This function is for NC != 2,3"
        NT = c.NT
        #NZ = c.NZ
        #NY = c.NY
        NX = c.NX
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        @inbounds c[1,1,ix,it] = a[1,1,ix,it]*b[1,1,ix,it]
                    end
                #end
            #end
        end
        set_wing_U!(c)

    end

    function LinearAlgebra.mul!(c::Gaugefields_2D_nowing{NC},a::T1,b::T2) where {NC,T1 <: Abstractfields,T2 <: Abstractfields}
        @assert NC != 2 && NC != 3 "This function is for NC != 2,3"
        NT = c.NT
        #NZ = c.NZ
        #NY = c.NY
        NX = c.NX
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        for k2=1:NC
                            for k1=1:NC
                                c[k1,k2,ix,it] = 0

                                @simd for k3=1:NC
                                    @inbounds c[k1,k2,ix,it] += a[k1,k3,ix,it]*b[k3,k2,ix,it]
                                end
                            end
                        end
                    end
                #end
            #end
        end
        set_wing_U!(c)
    end

    function LinearAlgebra.mul!(c::Gaugefields_2D_nowing{NC},a::T1,b::T2,iseven::Bool) where {NC,T1 <: Abstractfields,T2 <: Abstractfields}
        @assert NC != 2 && NC != 3 && NC != 1 "This function is for NC != 1,2,3"
        NT = c.NT
        #NZ = c.NZ
        #NY = c.NY
        NX = c.NX
        @inbounds for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        evenodd = ifelse( (ix+it) % 2 ==0, true,false)
                        if evenodd == iseven
                            for k2=1:NC
                                for k1=1:NC
                                    c[k1,k2,ix,it] = 0
                                    @simd for k3=1:NC
                                        c[k1,k2,ix,it] += a[k1,k3,ix,it]*b[k3,k2,ix,it]
                                    end
                                end
                            end
                        end
                    end
                #end
            #end
        end
        set_wing_U!(c)
    end

    function LinearAlgebra.mul!(c::Gaugefields_2D_nowing{NC},a::T1,b::T2) where {NC,T1 <: Number,T2 <: Abstractfields}
        @assert NC != 2 && NC != 3 "This function is for NC != 2,3"
        @inbounds for i=1:length(c)
            c.U[i] = a*b.U[i]
        end
        return

    end

    function LinearAlgebra.mul!(c::Gaugefields_2D_nowing{NC},a::T1,b::T2,α::Ta,β::Tb) where {NC,T1 <: Abstractfields,T2 <: Abstractfields,Ta <: Number, Tb <: Number}
        @assert NC != 2 && NC != 3 && NC != 1 "This function is for NC != 1,2,3"
        NT = c.NT
        #NZ = c.NZ
        #NY = c.NY
        NX = c.NX
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        for k2=1:NC
                            for k1=1:NC
                                c[k1,k2,ix,it] = β*c[k1,k2,ix,it]
                                @simd for k3=1:NC
                                    @inbounds c[k1,k2,ix,it] += α*a[k1,k3,ix,it]*b[k3,k2,ix,it]
                                end
                            end
                        end
                    end
                #end
            #end
        end
        set_wing_U!(c)
    end

    function LinearAlgebra.mul!(c::Gaugefields_2D_nowing{1},a::T1,b::T2,α::Ta,β::Tb) where {T1 <: Abstractfields,T2 <: Abstractfields,Ta <: Number, Tb <: Number}
        #@assert NC != 2 && NC != 3 "This function is for NC != 2,3"
        NT = c.NT
        #NZ = c.NZ
        #NY = c.NY
        NX = c.NX
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        @inbounds c[1,1,ix,it] = β*c[1,1,ix,it] + α*a[1,1,ix,it]*b[1,1,ix,it]
                    end
                #end
            #end
        end
        set_wing_U!(c)
    end

#=
evenodd = ifelse( (ix+iy+iz+it) % 2 ==0, true,false)
                        if evenodd == iseven
=#

    function LinearAlgebra.mul!(c::Gaugefields_2D_nowing{2},a::T1,b::T2,iseven::Bool) where {T1 <: Abstractfields,T2 <: Abstractfields}
        NT = c.NT
        #NZ = c.NZ
        #NY = c.NY
        NX = c.NX
        #println("threads = ", Threads.nthreads())
        @inbounds  for it=1:NT#,iz=1:NZ,iy=1:NY
            #for iz=1:NZ
                #for iy=1:NY
                    @simd for ix=1:NX
                        evenodd = ifelse( (ix+it) % 2 ==0, true,false)
                        if evenodd == iseven
                            a11 = a[1,1,ix,it]
                            a21 = a[2,1,ix,it]
                            a12 = a[1,2,ix,it]
                            a22 = a[2,2,ix,it]

                            b11 = b[1,1,ix,it]
                            b21 = b[2,1,ix,it]
                            b12 = b[1,2,ix,it]
                            b22 = b[2,2,ix,it]

                            c[1,1,ix,it] = a11*b11+a12*b21
                            c[2,1,ix,it] = a21*b11+a22*b21
                            c[1,2,ix,it] = a11*b12+a12*b22
                            c[2,2,ix,it] = a21*b12+a22*b22
                        end

                    end
                #end
            #end
        end
        set_wing_U!(c)
    end

    function LinearAlgebra.mul!(c::Gaugefields_2D_nowing{2},a::T1,b::T2) where {T1 <: Abstractfields,T2 <: Abstractfields}
        NT = c.NT
        #NZ = c.NZ
        #NY = c.NY
        NX = c.NX
        #println("threads = ", Threads.nthreads())
        @inbounds  for it=1:NT#,iz=1:NZ,iy=1:NY
            #for iz=1:NZ
                #for iy=1:NY
                    @simd for ix=1:NX
                        a11 = a[1,1,ix,it]
                        a21 = a[2,1,ix,it]
                        a12 = a[1,2,ix,it]
                        a22 = a[2,2,ix,it]

                        b11 = b[1,1,ix,it]
                        b21 = b[2,1,ix,it]
                        b12 = b[1,2,ix,it]
                        b22 = b[2,2,ix,it]

                        c[1,1,ix,it] = a11*b11+a12*b21
                        c[2,1,ix,it] = a21*b11+a22*b21
                        c[1,2,ix,it] = a11*b12+a12*b22
                        c[2,2,ix,it] = a21*b12+a22*b22

                    end
                #end
            #end
        end
        set_wing_U!(c)
    end

    function LinearAlgebra.mul!(c::Gaugefields_2D_nowing{2},a::T1,b::T2,α::Ta,β::Tb) where {NC,T1 <: Abstractfields,T2 <: Abstractfields,Ta <: Number, Tb <: Number}
        NT = c.NT
        #NZ = c.NZ
        #NY = c.NY
        NX = c.NX
        if β == zero(β)
            if α == one(α)
                mul!(c,a,b)
                return
            end
        end

        if α != 0
            @inbounds for it=1:NT
                #for iz=1:NZ
                    #for iy=1:NY
                        @simd for ix=1:NX
                            a11 = a[1,1,ix,it]
                            a21 = a[2,1,ix,it]
                            a12 = a[1,2,ix,it]
                            a22 = a[2,2,ix,it]

                            b11 = b[1,1,ix,it]
                            b21 = b[2,1,ix,it]
                            b12 = b[1,2,ix,it]
                            b22 = b[2,2,ix,it]

                            c[1,1,ix,it] = (a11*b11+a12*b21)*α + β*c[1,1,ix,it]
                            c[2,1,ix,it] = (a21*b11+a22*b21)*α + β*c[2,1,ix,it]
                            c[1,2,ix,it] = (a11*b12+a12*b22)*α + β*c[1,2,ix,it]
                            c[2,2,ix,it] = (a21*b12+a22*b22)*α + β*c[2,2,ix,it]

                        end
                    #end
                #end
            end
        end
        set_wing_U!(c)
    end

    function LinearAlgebra.mul!(c::Gaugefields_2D_nowing{3},a::T1,b::T2) where {T1 <: Abstractfields,T2 <: Abstractfields}
        NT = c.NT
        #NZ = c.NZ
        #NY = c.NY
        NX = c.NX
        #println("threads = ", Threads.nthreads())
        @inbounds  for it=1:NT#,iz=1:NZ,iy=1:NY
            #for iz=1:NZ
                #for iy=1:NY
                    @simd for ix=1:NX
                        a11 = a[1,1,ix,it]
                        a21 = a[2,1,ix,it]
                        a31 = a[3,1,ix,it]
                        a12 = a[1,2,ix,it]
                        a22 = a[2,2,ix,it]
                        a32 = a[3,2,ix,it]
                        a13 = a[1,3,ix,it]
                        a23 = a[2,3,ix,it]
                        a33 = a[3,3,ix,it]
                        b11 = b[1,1,ix,it]
                        b21 = b[2,1,ix,it]
                        b31 = b[3,1,ix,it]
                        b12 = b[1,2,ix,it]
                        b22 = b[2,2,ix,it]
                        b32 = b[3,2,ix,it]
                        b13 = b[1,3,ix,it]
                        b23 = b[2,3,ix,it]
                        b33 = b[3,3,ix,it]
                        c[1,1,ix,it] = a11*b11+a12*b21+a13*b31
                        c[2,1,ix,it] = a21*b11+a22*b21+a23*b31
                        c[3,1,ix,it] = a31*b11+a32*b21+a33*b31
                        c[1,2,ix,it] = a11*b12+a12*b22+a13*b32
                        c[2,2,ix,it] = a21*b12+a22*b22+a23*b32
                        c[3,2,ix,it] = a31*b12+a32*b22+a33*b32
                        c[1,3,ix,it] = a11*b13+a12*b23+a13*b33
                        c[2,3,ix,it] = a21*b13+a22*b23+a23*b33
                        c[3,3,ix,it] = a31*b13+a32*b23+a33*b33

                    end
                #end
            #end
        end
        set_wing_U!(c)
        #error("mul")
    end

    function LinearAlgebra.mul!(c::Gaugefields_2D_nowing{3},a::T1,b::T2,iseven::Bool) where {T1 <: Abstractfields,T2 <: Abstractfields}
        NT = c.NT
        #NZ = c.NZ
        #NY = c.NY
        NX = c.NX
        #println("threads = ", Threads.nthreads())
        @inbounds  for it=1:NT#,iz=1:NZ,iy=1:NY
            #for iz=1:NZ
                #for iy=1:NY
                    @simd for ix=1:NX
                        evenodd = ifelse( (ix+it) % 2 ==0, true,false)
                        if evenodd == iseven

                            a11 = a[1,1,ix,it]
                            a21 = a[2,1,ix,it]
                            a31 = a[3,1,ix,it]
                            a12 = a[1,2,ix,it]
                            a22 = a[2,2,ix,it]
                            a32 = a[3,2,ix,it]
                            a13 = a[1,3,ix,it]
                            a23 = a[2,3,ix,it]
                            a33 = a[3,3,ix,it]
                            b11 = b[1,1,ix,it]
                            b21 = b[2,1,ix,it]
                            b31 = b[3,1,ix,it]
                            b12 = b[1,2,ix,it]
                            b22 = b[2,2,ix,it]
                            b32 = b[3,2,ix,it]
                            b13 = b[1,3,ix,it]
                            b23 = b[2,3,ix,it]
                            b33 = b[3,3,ix,it]
                            c[1,1,ix,it] = a11*b11+a12*b21+a13*b31
                            c[2,1,ix,it] = a21*b11+a22*b21+a23*b31
                            c[3,1,ix,it] = a31*b11+a32*b21+a33*b31
                            c[1,2,ix,it] = a11*b12+a12*b22+a13*b32
                            c[2,2,ix,it] = a21*b12+a22*b22+a23*b32
                            c[3,2,ix,it] = a31*b12+a32*b22+a33*b32
                            c[1,3,ix,it] = a11*b13+a12*b23+a13*b33
                            c[2,3,ix,it] = a21*b13+a22*b23+a23*b33
                            c[3,3,ix,it] = a31*b13+a32*b23+a33*b33
                        end

                    end
                #end
            #end
        end
        set_wing_U!(c)
    end

    function LinearAlgebra.mul!(c::Gaugefields_2D_nowing{3},a::T1,b::T2,α::Ta,β::Tb) where {NC,T1 <: Abstractfields,T2 <: Abstractfields,Ta <: Number, Tb <: Number}
        NT = c.NT
        #NZ = c.NZ
        #NY = c.NY
        NX = c.NX
        if β == zero(β)
            if α == one(α)
                mul!(c,a,b)
                return
            end
        end

        if α != 0
            @inbounds for it=1:NT
                #for iz=1:NZ
                    #for iy=1:NY
                        @simd for ix=1:NX
                            a11 = a[1,1,ix,it]
                            a21 = a[2,1,ix,it]
                            a31 = a[3,1,ix,it]
                            a12 = a[1,2,ix,it]
                            a22 = a[2,2,ix,it]
                            a32 = a[3,2,ix,it]
                            a13 = a[1,3,ix,it]
                            a23 = a[2,3,ix,it]
                            a33 = a[3,3,ix,it]
                            b11 = b[1,1,ix,it]
                            b21 = b[2,1,ix,it]
                            b31 = b[3,1,ix,it]
                            b12 = b[1,2,ix,it]
                            b22 = b[2,2,ix,it]
                            b32 = b[3,2,ix,it]
                            b13 = b[1,3,ix,it]
                            b23 = b[2,3,ix,it]
                            b33 = b[3,3,ix,it]
                            c[1,1,ix,it] = (a11*b11+a12*b21+a13*b31)*α + β*c[1,1,ix,it]
                            c[2,1,ix,it] = (a21*b11+a22*b21+a23*b31)*α + β*c[2,1,ix,it]
                            c[3,1,ix,it] = (a31*b11+a32*b21+a33*b31)*α + β*c[3,1,ix,it]
                            c[1,2,ix,it] = (a11*b12+a12*b22+a13*b32)*α + β*c[1,2,ix,it]
                            c[2,2,ix,it] = (a21*b12+a22*b22+a23*b32)*α + β*c[2,2,ix,it]
                            c[3,2,ix,it] = (a31*b12+a32*b22+a33*b32)*α + β*c[3,2,ix,it]
                            c[1,3,ix,it] = (a11*b13+a12*b23+a13*b33)*α + β*c[1,3,ix,it]
                            c[2,3,ix,it] = (a21*b13+a22*b23+a23*b33)*α + β*c[2,3,ix,it]
                            c[3,3,ix,it] = (a31*b13+a32*b23+a33*b33)*α + β*c[3,3,ix,it]

                        end
                    #end
                #end
            end
        end
        set_wing_U!(c)
    end

    function mul_skiplastindex!(c::Gaugefields_2D_nowing{NC},a::T1,b::T2) where {NC,T1 <: Abstractfields,T2 <: Abstractfields}
        #@assert NC != 2 && NC != 3 "This function is for NC != 2,3"
        NT = c.NT
        #NZ = c.NZ
        #NY = c.NY
        NX = c.NX
        #for it=1:NT
        it = 1
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        for k2=1:NC
                            for k1=1:NC
                                c[k1,k2,ix,it] = 0

                                @simd for k3=1:NC
                                    c[k1,k2,ix,it] += a[k1,k3,ix,it]*b[k3,k2,ix,it]
                                end
                            end
                        end
                    end
                #end
            #end
        #end
    end

    function normalize_U!(u::Gaugefields_2D_nowing{NC}) where NC
        NX = u.NX
        #NY = u.NY
        #NZ = u.NZ
        NT = u.NT

        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    @simd for ix=1:NX
                        A = u[:,:,ix,it]
                        gramschmidt!(A)
                        u[:,:,ix,it] = A[:,:]
                    end
                #end
            #end
        end

    end

    function normalize_U!(u::Gaugefields_2D_nowing{2})
        NX = u.NX
        #NY = u.NY
        #NZ = u.NZ
        NT = u.NT

        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX

                        α = u[1,1,ix,it]
                        β = u[2,1,ix,it]
                        detU = sqrt(abs(α)^2 + abs(β)^2)

                        u[1,1,ix,it] = α/detU
                        u[2,1,ix,it] = β/detU
                        u[1,2,ix,it] = -conj(β)/detU
                        u[2,2,ix,it] = conj(α)/detU

                    end
                #end
            #end
        end
    end

    function normalize_U!(u::Gaugefields_2D_nowing{3})
        NX = u.NX
        #NY = u.NY
        #NZ = u.NZ
        NT = u.NT

        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX

                        w1 = 0
                        w2 = 0
                        @simd for ic=1:3
                            w1 += u[2,ic,ix,it]*conj(u[1,ic,ix,it])
                            w2 += u[1,ic,ix,it]*conj(u[1,ic,ix,it])
                        end
                        zerock2 = w2
                        if zerock2 == 0
                            println("w2 is zero  !!  (in normlz)")
                            println("u[1,1),u[1,2),u[1,3) : ",u[1,1,ix,it], "\t",u[1,2,ix,it],"\t", u[1,3,ix,it])
                        end

                        w1 = -w1/w2

                        x4 = (u[2,1,ix,it]) + w1*u[1,1,ix,it]
                        x5 = (u[2,2,ix,it]) + w1*u[1,2,ix,it]
                        x6 = (u[2,3,ix,it]) + w1*u[1,3,ix,it]

                        w3 = x4*conj(x4) + x5*conj(x5) + x6*conj(x6)

                        zerock3 = w3
                        if zerock3 == 0
                            println("w3 is zero  !!  (in normlz)")
                            println("x4, x5, x6 : $x4, $x5, $x6")
                            exit()
                        end

                        u[2,1,ix,it] = x4
                        u[2,2,ix,it] = x5
                        u[2,3,ix,it] = x6

                        w3 = 1/sqrt(w3)
                        w2 = 1/sqrt(w2)

                        u[1,1,ix,it] = u[1,1,ix,it]*w2
                        u[1,2,ix,it] = u[1,2,ix,it]*w2
                        u[1,3,ix,it] = u[1,3,ix,it]*w2
                        u[2,1,ix,it] = u[2,1,ix,it]*w3
                        u[2,2,ix,it] = u[2,2,ix,it]*w3
                        u[2,3,ix,it] = u[2,3,ix,it]*w3

                        if zerock2*zerock3 == 0
                            println("!! devided by zero !! (in normalize)")
                            println("w2 or w3 in normlz is zero !!")
                            println("w2, w3 : $w2, $w3   ")
                            exit()
                        end
                        #println(u[:,:,ix,it]'*u[:,:,ix,it])
                    end
                #end
            #end
        end
        m3complv!(u)
    end

    function m3complv!(a::Gaugefields_2D_nowing{3})
        aa = zeros(Float64,18)
        NX = a.NX
        #NY = a.NY
        #NZ = a.NZ
        NT = a.NT

        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    @simd for ix=1:NX

                        aa[ 1] = real( a[1,1,ix,it])
                        aa[ 2] = imag(a[1,1,ix,it])
                        aa[ 3] = real( a[1,2,ix,it])
                        aa[ 4] = imag(a[1,2,ix,it])
                        aa[ 5] = real( a[1,3,ix,it])
                        aa[ 6] = imag(a[1,3,ix,it])
                        aa[ 7] = real( a[2,1,ix,it])
                        aa[ 8] = imag(a[2,1,ix,it])
                        aa[ 9] = real( a[2,2,ix,it])
                        aa[10] = imag(a[2,2,ix,it])
                        aa[11] = real( a[2,3,ix,it])
                        aa[12] = imag(a[2,3,ix,it])

                        aa[13] = aa[ 3]*aa[11] - aa[ 4]*aa[12] -
                                    aa[ 5]*aa[ 9] + aa[ 6]*aa[10]
                        aa[14] = aa[ 5]*aa[10] + aa[ 6]*aa[ 9] -
                                    aa[ 3]*aa[12] - aa[ 4]*aa[11]
                        aa[15] = aa[ 5]*aa[ 7] - aa[ 6]*aa[ 8] -
                                    aa[ 1]*aa[11] + aa[ 2]*aa[12]
                        aa[16] = aa[ 1]*aa[12] + aa[ 2]*aa[11] -
                                    aa[ 5]*aa[ 8] - aa[ 6]*aa[ 7]
                        aa[17] = aa[ 1]*aa[ 9] - aa[ 2]*aa[10] -
                                    aa[ 3]*aa[ 7] + aa[ 4]*aa[ 8]
                        aa[18] = aa[ 3]*aa[ 8] + aa[ 4]*aa[ 7] -
                                    aa[ 1]*aa[10] - aa[ 2]*aa[ 9]

                        a[3,1,ix,it] = aa[13]+im*aa[14]
                        a[3,2,ix,it] = aa[15]+im*aa[16]
                        a[3,3,ix,it] = aa[17]+im*aa[18]

                        #println(a[:,:,ix,it]'*a[:,:,ix,it] )
                    end
                #end
            #end
        end
    end

# end


text:
# Support for ZeroMQ, a network and interprocess communication library

VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module ZMQ
using Compat
import Compat: String, unsafe_string
if VERSION >= v"0.4.0-dev+3710"
    import Base.unsafe_convert
else
    const unsafe_convert = Base.convert
end
if VERSION >= v"0.4.0-dev+3844"
    using Base.Libdl, Base.Libc
    using Base.Libdl: dlopen_e
    using Base.Libc: EAGAIN
else
    using Base: EAGAIN
end
if VERSION >= v"0.5.0-dev+1229"
    import Base.Filesystem: UV_READABLE, uv_pollcb
else
    import Base: UV_READABLE
    if isdefined(Base, :uv_pollcb)
        import Base: uv_pollcb
    end
end

const depfile = joinpath(dirname(@__FILE__),"..","deps","deps.jl")
if isfile(depfile)
    include(depfile)
else
    error("ZMQ not properly installed. Please run Pkg.build(\"ZMQ\")")
end

import Base:
    convert, get, bytestring,
    length, size, stride, similar, getindex, setindex!,
    fd, wait, notify, close, connect,
    bind, send, recv

export
    #Types
    StateError,Context,Socket,Message,
    #functions
    set, subscribe, unsubscribe,
    #Constants
    IO_THREADS,MAX_SOCKETS,PAIR,PUB,SUB,REQ,REP,ROUTER,DEALER,PULL,PUSH,XPUB,XSUB,XREQ,XREP,UPSTREAM,DOWNSTREAM,MORE,POLLIN,POLLOUT,POLLERR,STREAMER,FORWARDER,QUEUE,SNDMORE
const SNDMORE = true

# A server will report most errors to the client over a Socket, but
# errors in ZMQ state can't be reported because the socket may be
# corrupted. Therefore, we need an exception type for errors that
# should be reported locally.
type StateError <: Exception
    msg::AbstractString
end
show(io, thiserr::StateError) = print(io, "ZMQ: ", thiserr.msg)

# Basic functions
zmq_errno() = ccall((:zmq_errno, zmq), Cint, ())
function jl_zmq_error_str()
    errno = zmq_errno()
    c_strerror = ccall ((:zmq_strerror, zmq), Ptr{UInt8}, (Cint,), errno)
    if c_strerror != C_NULL
        strerror = unsafe_string(c_strerror)
        return strerror
    else
        return "Unknown error"
    end
end

if VERSION >= v"0.5-" && isdefined(Base, :Filesystem)
    if is_windows()
        using Base.Libc: WindowsRawSocket
    end
    const _FDWatcher = Base.Filesystem._FDWatcher
    const _have_good_fdwatcher = true
elseif VERSION >= v"0.4-" && isdefined(Base, :_FDWatcher)
    if is_windows()
        using Base.Libc: WindowsRawSocket
    end
    const _FDWatcher = Base._FDWatcher
    const _have_good_fdwatcher = true
else
    if is_windows()
        using Base: WindowsRawSocket
    end
    const _FDWatcher = Base.FDWatcher
    const _have_good_fdwatcher = false
end

## Sockets ##
type Socket
    data::Ptr{Void}
    pollfd::_FDWatcher

    # ctx should be ::Context, but forward type references are not allowed
    function Socket(ctx, typ::Integer)
        p = ccall((:zmq_socket, zmq), Ptr{Void}, (Ptr{Void}, Cint), ctx.data, typ)
        if p == C_NULL
            throw(StateError(jl_zmq_error_str()))
        end
        socket = new(p)
        if _have_good_fdwatcher
            socket.pollfd = _FDWatcher(fd(socket), #=readable=#true, #=writable=#false)
        else
            socket.pollfd = _FDWatcher(fd(socket))
        end
        finalizer(socket, close)
        push!(ctx.sockets, socket)
        return socket
    end
end

function close(socket::Socket)
    if socket.data != C_NULL
        data = socket.data
        socket.data = C_NULL
        if _have_good_fdwatcher
            close(socket.pollfd, #=readable=#true, #=writable=#false)
        else
            notify(socket.pollfd.notify)
            Base.stop_watching(socket.pollfd)
        end
        rc = ccall((:zmq_close, zmq), Cint,  (Ptr{Void},), data)
        if rc != 0
            throw(StateError(jl_zmq_error_str()))
        end
    end
end

## Contexts ##
# Provide the same constructor API for version 2 and version 3, even
# though the underlying functions are changing
type Context
    data::Ptr{Void}

    # need to keep a list of sockets for this Context in order to
    # close them before finalizing (otherwise zmq_term will hang)
    sockets::Vector{Socket}

    function Context()
        p = ccall((:zmq_ctx_new, zmq), Ptr{Void},  ())
        if p == C_NULL
            throw(StateError(jl_zmq_error_str()))
        end
        zctx = new(p, Array(Socket,0))
        finalizer(zctx, close)
        return zctx
    end
end

@deprecate Context(n::Integer) Context()

function close(ctx::Context)
    if ctx.data != C_NULL # don't close twice!
        data = ctx.data
        ctx.data = C_NULL
        for s in ctx.sockets
            close(s)
        end
        rc = ccall((:zmq_ctx_destroy, zmq), Cint,  (Ptr{Void},), data)
        if rc != 0
            throw(StateError(jl_zmq_error_str()))
        end
    end
end
term(ctx::Context) = close(ctx)

function get(ctx::Context, option::Integer)
    val = ccall((:zmq_ctx_get, zmq), Cint, (Ptr{Void}, Cint), ctx.data, option)
    if val < 0
        throw(StateError(jl_zmq_error_str()))
    end
    return val
end

function set(ctx::Context, option::Integer, value::Integer)
    rc = ccall((:zmq_ctx_set, zmq), Cint, (Ptr{Void}, Cint, Cint), ctx.data, option, value)
    if rc != 0
        throw(StateError(jl_zmq_error_str()))
    end
end

# Getting and setting socket options
# Socket options of integer type
const u64p = zeros(UInt64, 1)
const i64p = zeros(Int64, 1)
const ip = zeros(Cint, 1)
const u32p = zeros(UInt32, 1)
const sz = zeros(UInt, 1)
const pp = fill(C_NULL, 1)

for (fset, fget, k, p) in [
    (:set_affinity,                :get_affinity,                 4, u64p)
    (:set_type,                    :get_type,                    16,   ip)
    (:set_linger,                  :get_linger,                  17,   ip)
    (:set_reconnect_ivl,           :get_reconnect_ivl,           18,   ip)
    (:set_backlog,                 :get_backlog,                 19,   ip)
    (:set_reconnect_ivl_max,       :get_reconnect_ivl_max,       21,   ip)
    (:set_rate,                    :get_rate,                     8,   ip)
    (:set_recovery_ivl,            :get_recovery_ivl,             9,   ip)
    (:set_sndbuf,                  :get_sndbuf,                  11,   ip)
    (:set_rcvbuf,                  :get_rcvbuf,                  12,   ip)
    (nothing,                      :_zmq_getsockopt_rcvmore,     13,   ip)
    (nothing,                      :get_events,                  15,   ip)
    (:set_maxmsgsize,              :get_maxmsgsize,              22,   ip)
    (:set_sndhwm,                  :get_sndhwm,                  23,   ip)
    (:set_rcvhwm,                  :get_rcvhwm,                  24,   ip)
    (:set_multicast_hops,          :get_multicast_hops,          25,   ip)
    (:set_ipv4only,                :get_ipv4only,                31,   ip)
    (:set_tcp_keepalive,           :get_tcp_keepalive,           34,   ip)
    (:set_tcp_keepalive_idle,      :get_tcp_keepalive_idle,      35,   ip)
    (:set_tcp_keepalive_cnt,       :get_tcp_keepalive_cnt,       36,   ip)
    (:set_tcp_keepalive_intvl,     :get_tcp_keepalive_intvl,     37,   ip)
    (:set_rcvtimeo,                :get_rcvtimeo,                27,   ip)
    (:set_sndtimeo,                :get_sndtimeo,                28,   ip)
    (nothing,                      :get_fd,                      14, is_windows() ? pp : ip)
    ]
    if fset != nothing
        @eval function ($fset)(socket::Socket, option_val::Integer)
            ($p)[1] = option_val
            rc = ccall((:zmq_setsockopt, zmq), Cint,
                       (Ptr{Void}, Cint, Ptr{Void}, UInt),
                       socket.data, $k, $p, sizeof(eltype($p)))
            if rc != 0
                throw(StateError(jl_zmq_error_str()))
            end
        end
    end
    if fget != nothing
        @eval function ($fget)(socket::Socket)
            ($sz)[1] = sizeof(eltype($p))
            rc = ccall((:zmq_getsockopt, zmq), Cint,
                       (Ptr{Void}, Cint, Ptr{Void}, Ptr{UInt}),
                       socket.data, $k, $p, $sz)
            if rc != 0
                throw(StateError(jl_zmq_error_str()))
            end
            return @compat Int(($p)[1])
        end
    end
end

# For some functions, the publicly-visible versions should require &
# return boolean:
get_rcvmore(socket::Socket) = @compat Bool(_zmq_getsockopt_rcvmore(socket))
# And a convenience function
ismore(socket::Socket) = get_rcvmore(socket)

# subscribe/unsubscribe options take an arbitrary byte array
for (f,k) in ((:subscribe,6), (:unsubscribe,7))
    f_ = @compat Symbol(string(f, "_"))
    @eval begin
        function $f_{T}(socket::Socket, filter::Ptr{T}, len::Integer)
            rc = ccall((:zmq_setsockopt, zmq), Cint,
                       (Ptr{Void}, Cint, Ptr{T}, UInt),
                       socket.data, $k, filter, len)
            if rc != 0
                throw(StateError(jl_zmq_error_str()))
            end
        end
        @compat $f(socket::Socket, filter::Union{Array,AbstractString}) =
            $f_(socket, pointer(filter), sizeof(filter))
        $f(socket::Socket) = $f_(socket, C_NULL, 0)
    end
end

# Raw FD access
if is_unix()
    fd(socket::Socket) = RawFD(get_fd(socket))
end
if is_windows()
    fd(socket::Socket) = WindowsRawSocket(convert(Ptr{Void}, get_fd(socket)))
end

if _have_good_fdwatcher
    wait(socket::Socket) = wait(socket.pollfd, readable=true, writable=false)
    notify(socket::Socket) = uv_pollcb(socket.pollfd.handle, Int32(0),
                                       Int32(UV_READABLE))
else
    wait(socket::Socket) = Base._wait(socket.pollfd, #=readable=#true, #=writable=#false)
    notify(socket::Socket) = Base._uv_hook_pollcb(socket.pollfd, int32(0),
                                                  int32(UV_READABLE))
end

# Socket options of string type
const u8ap = zeros(UInt8, 255)
for (fset, fget, k) in [
    (:set_identity,                :get_identity,                5)
    (:set_subscribe,               nothing,                      6)
    (:set_unsubscribe,             nothing,                      7)
    (nothing,                      :get_last_endpoint,          32)
    (:set_tcp_accept_filter,       nothing,                     38)
    ]
    if fset != nothing
        @eval function ($fset)(socket::Socket, option_val::String)
            if length(option_val) > 255
                throw(StateError("option value too large"))
            end
            rc = ccall((:zmq_setsockopt, zmq), Cint,
                       (Ptr{Void}, Cint, Ptr{UInt8}, UInt),
                       socket.data, $k, option_val, length(option_val))
            if rc != 0
                throw(StateError(jl_zmq_error_str()))
            end
        end
    end
    if fget != nothing
        @eval function ($fget)(socket::Socket)
            ($sz)[1] = length($u8ap)
            rc = ccall((:zmq_getsockopt, zmq), Cint,
                       (Ptr{Void}, Cint, Ptr{UInt8}, Ptr{UInt}),
                       socket.data, $k, $u8ap, $sz)
            if rc != 0
                throw(StateError(jl_zmq_error_str()))
            end
            return unsafe_string(unsafe_convert(Ptr{UInt8}, $u8ap), @compat Int(($sz)[1]))
        end
    end
end

function bind(socket::Socket, endpoint::AbstractString)
    rc = ccall((:zmq_bind, zmq), Cint, (Ptr{Void}, Ptr{UInt8}), socket.data, endpoint)
    if rc != 0
        throw(StateError(jl_zmq_error_str()))
    end
end

function connect(socket::Socket, endpoint::AbstractString)
    rc=ccall((:zmq_connect, zmq), Cint, (Ptr{Void}, Ptr{UInt8}), socket.data, endpoint)
    if rc != 0
        throw(StateError(jl_zmq_error_str()))
    end
end

# in order to support zero-copy messages that share data with Julia
# arrays, we need to hold a reference to the Julia object in a dictionary
# until zeromq is done with the data, to prevent it from being garbage
# collected.  The gc_protect dictionary is keyed by a uv_async_t* pointer,
# used in uv_async_send to tell Julia to when zeromq is done with the data.
const gc_protect = Dict{Ptr{Void},Any}()
# 0.2 compatibility
gc_protect_cb(work, status) = gc_protect_cb(work)
if VERSION < v"0.4.0-dev+3970"
    function close_handle(work)
        Base.disassociate_julia_struct(work.handle)
        ccall(:jl_close_uv,Void,(Ptr{Void},),work.handle)
        Base.unpreserve_handle(work)
    end
else
    close_handle(work) = Base.close(work)
end
gc_protect_cb(work) = (pop!(gc_protect, work.handle, nothing); close_handle(work))

function gc_protect_handle(obj::Any)
    work = Compat.AsyncCondition(gc_protect_cb)
    gc_protect[work.handle] = (work,obj)
    work.handle
end

# Thread-safe zeromq callback when data is freed, passed to zmq_msg_init_data.
# The hint parameter will be a uv_async_t* pointer.
function gc_free_fn(data::Ptr{Void}, hint::Ptr{Void})
    ccall(:uv_async_send,Cint,(Ptr{Void},),hint)
end

## Messages ##
bitstype 64 * 8 MsgPadding

type Message <: AbstractArray{UInt8,1}
    # Matching the declaration in the header: char _[64];
    w_padding::MsgPadding
    handle::Ptr{Void} # index into gc_protect, if any

    # Create an empty message (for receive)
    function Message()
        zmsg = new()
        zmsg.handle = C_NULL
        rc = ccall((:zmq_msg_init, zmq), Cint, (Ptr{Message},), &zmsg)
        if rc != 0
            throw(StateError(jl_zmq_error_str()))
        end
        finalizer(zmsg, close)
        return zmsg
    end
    # Create a message with a given buffer size (for send)
    function Message(len::Integer)
        zmsg = new()
        zmsg.handle = C_NULL
        rc = ccall((:zmq_msg_init_size, zmq), Cint, (Ptr{Message}, Csize_t), &zmsg, len)
        if rc != 0
            throw(StateError(jl_zmq_error_str()))
        end
        finalizer(zmsg, close)
        return zmsg
    end

    # low-level function to create a message (for send) with an existing
    # data buffer, without making a copy.  The origin parameter should
    # be the Julia object that is the origin of the data, so that
    # we can hold a reference to it until zeromq is done with the buffer.
    function Message{T}(origin::Any, m::Ptr{T}, len::Integer)
        zmsg = new()
        zmsg.handle = gc_protect_handle(origin)
        rc = ccall((:zmq_msg_init_data, zmq), Cint, (Ptr{Message}, Ptr{T}, Csize_t, Ptr{Void}, Ptr{Void}), &zmsg, m, len, gc_free_fn_c::Ptr{Void}, zmsg.handle)
        if rc != 0
            throw(StateError(jl_zmq_error_str()))
        end
        finalizer(zmsg, close)
        return zmsg
    end

    # Create a message with a given AbstractString or Array as a buffer (for send)
    # (note: now "owns" the buffer ... the Array must not be resized,
    #        or even written to after the message is sent!)
    Message(m::String) = Message(m, unsafe_convert(Ptr{UInt8}, pointer(m)), sizeof(m))
    Message{T<:String}(p::SubString{T}) =
        Message(p, pointer(p.string.data)+p.offset, sizeof(p))
    Message(a::Array) = Message(a, pointer(a), sizeof(a))
    function Message(io::IOBuffer)
        if !io.readable || !io.seekable
            error("byte read failed")
        end
        Message(io.data)
    end
end

# check whether zeromq has called our free-function, i.e. whether
# we are save to reclaim ownership of any buffer object
isfreed(m::Message) = haskey(gc_protect, m.handle)

# AbstractArray behaviors:
similar(a::Message, T, dims::Dims) = Array(T, dims) # ?
length(zmsg::Message) = @compat Int(ccall((:zmq_msg_size, zmq), Csize_t, (Ptr{Message},), &zmsg))
size(zmsg::Message) = (length(zmsg),)
unsafe_convert(::Type{Ptr{UInt8}}, zmsg::Message) = ccall((:zmq_msg_data, zmq), Ptr{UInt8}, (Ptr{Message},), &zmsg)
function getindex(a::Message, i::Integer)
    if i < 1 || i > length(a)
        throw(BoundsError())
    end
    unsafe_load(pointer(a), i)
end
function setindex!(a::Message, v, i::Integer)
    if i < 1 || i > length(a)
        throw(BoundsError())
    end
    unsafe_store(pointer(a), v, i)
end

# Convert message to string (copies data)
unsafe_string(zmsg::Message) = Compat.unsafe_string(pointer(zmsg), length(zmsg))
if VERSION < v"0.5-dev+4341"
    bytestring(zmsg::Message) = unsafe_string(zmsg)
else
    @deprecate bytestring(zmsg::Message) unsafe_string(zmsg::Message)
end

# Build an IOStream from a message
# Copies the data
function convert(::Type{IOStream}, zmsg::Message)
    s = IOBuffer()
    write(s, zmsg)
    return s
end
# Close a message. You should not need to call this manually (let the
# finalizer do it).
function close(zmsg::Message)
    rc = ccall((:zmq_msg_close, zmq), Cint, (Ptr{Message},), &zmsg)
    if rc != 0
        throw(StateError(jl_zmq_error_str()))
    end
end

function get(zmsg::Message, property::Integer)
    val = ccall((:zmq_msg_get, zmq), Cint, (Ptr{Message}, Cint), &zmsg, property)
    if val < 0
        throw(StateError(jl_zmq_error_str()))
    end
    val
end
function set(zmsg::Message, property::Integer, value::Integer)
    rc = ccall((:zmq_msg_set, zmq), Cint, (Ptr{Message}, Cint, Cint), &zmsg, property, value)
    if rc < 0
        throw(StateError(jl_zmq_error_str()))
    end
end

## Send/receive messages
#
# Julia defines two types of ZMQ messages: "raw" and "serialized". A "raw"
# message is just a plain ZeroMQ message, used for sending a sequence
# of bytes. You send these with the following:
#   send(socket, zmsg)
#   zmsg = recv(socket)

#Send/Recv Options
const ZMQ_DONTWAIT = 1
const ZMQ_SNDMORE = 2

function send(socket::Socket, zmsg::Message, SNDMORE::Bool=false)
    while true
        rc = ccall((:zmq_msg_send, zmq), Cint, (Ptr{Message}, Ptr{Void}, Cint),
                    &zmsg, socket.data, (ZMQ_SNDMORE*SNDMORE) | ZMQ_DONTWAIT)
        if rc == -1
            zmq_errno() == EAGAIN || throw(StateError(jl_zmq_error_str()))
            while (get_events(socket) & POLLOUT) == 0
                wait(socket)
            end
        else
            get_events(socket) != 0 && notify(socket)
            break
        end
    end
end

# strings are immutable, so we can send them zero-copy by default
send(socket::Socket, msg::AbstractString, SNDMORE::Bool=false) = send(socket, Message(msg), SNDMORE)

# Make a copy of arrays before sending, by default, since it is too
# dangerous to require that the array not change until ZMQ is done with it.
# For zero-copy array messages, construct a Message explicitly.
send(socket::Socket, msg::AbstractArray, SNDMORE::Bool=false) = send(socket, Message(copy(msg)), SNDMORE)

function send(f::Function, socket::Socket, SNDMORE::Bool=false)
    io = IOBuffer()
    f(io)
    send(socket, Message(io), SNDMORE)
end

function recv(socket::Socket)
    zmsg = Message()
    rc = -1
    while true
        rc = ccall((:zmq_msg_recv, zmq), Cint, (Ptr{Message}, Ptr{Void}, Cint),
                    &zmsg, socket.data, ZMQ_DONTWAIT)
        if rc == -1
            zmq_errno() == EAGAIN || throw(StateError(jl_zmq_error_str()))
            while (get_events(socket) & POLLIN) == 0
                wait(socket)
            end
        else
            get_events(socket) != 0 && notify(socket)
            break
        end
    end
    return zmsg
end

## Constants

# Context options
const IO_THREADS = 1
const MAX_SOCKETS = 2

#Socket Types
const PAIR = 0
const PUB = 1
const SUB = 2
const REQ = 3
const REP = 4
const DEALER = 5
const ROUTER = 6
const PULL = 7
const PUSH = 8
const XPUB = 9
const XSUB = 10
const XREQ = DEALER
const XREP = ROUTER
const UPSTREAM = PULL
const DOWNSTREAM = PUSH

#Message options
const MORE = 1

#IO Multiplexing
const POLLIN = 1
const POLLOUT = 2
const POLLERR = 4

#Built in devices
const STREAMER = 1
const FORWARDER = 2
const QUEUE = 3

function __init__()
    major = Array(Cint,1)
    minor = Array(Cint,1)
    patch = Array(Cint,1)
    ccall((:zmq_version, zmq), Void, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), major, minor, patch)
    global const version = VersionNumber(major[1], minor[1], patch[1])
    if version < v"3"
        error("ZMQ version $version < 3 is not supported")
    end
    global const gc_free_fn_c = cfunction(gc_free_fn, Cint, (Ptr{Void}, Ptr{Void}))
end

end

ext: jl


text:
__precompile__()

module Boot

  using Compat
  using StringCases

  include("file_util_types.jl")

  include("check_module_node.jl")
  include("purge_corrupted_methods.jl")
  include("strip_output_color.jl")

  include("do_main_load_loop.jl")
  include("make_initial_load.jl")
  include("attempt_file_load.jl")

  include("print_load_info.jl")
  include("start_load_failure.jl")
  include("load_invalid_file.jl")

  include("shard_list_constants.jl")
  include("attempt_shard_load.jl")
  include("clean_shard.jl")

  include("include_folder.jl")

  include("get_method_count.jl")
  include("get_all_files.jl")
  include("parse_file.jl")

  export include_folder

  export get_method_count
  export get_all_files
  export parse_file

end


text:
# geometry engine for vm2d
import LinearAlgebra
la = LinearAlgebra
include("vm2d_particle.jl")

# circle
"""
`circle(center,radius,points,spacing)`

returns an array of 2-element arrays describing the coordinates of a circle
"""
function circlexy(center::Array{Float64,1},radius::Float64,refinement::Int64,thetastart::Float64,thetaend::Float64)
    thetas = range(thetastart,stop=thetaend,length=refinement)
    points = Array{Float64,1}[]
    if abs(thetaend - 2*pi) < 0.00000001
        thetas = thetas[1:end-1]
    end
    for theta in thetas
        x = center[1] + radius*cos(theta)
        y = center[2] + radius*sin(theta)
        if length(center) == 3
            z = center[3]
            push!(points,[x,y,z])
        elseif length(center) == 2
            push!(points,[x,y,0])
        else
            throw(ArgumentError(center, "argument must be a 2- or 3-element array"))
        end
    end

    return points
end

# ring
"""
`vortexring(radius::Float64,center::Array{Float64,1},refinement::Int64,gamma::Float64; velocity::Array{Float64,1}=[0.0,0.0,0.0])`

* returns an array of vortex particles defining a vortex ring in an xy-plane
* positive values of gamma produce clockwise vorticity from the vantage a traditional x-y view
"""
function vortexring(radius::Float64,center::Array{Float64,1},refinement::Int64,gamma::Float64;velocity::Array{Float64,1}=[0.0,0.0,0.0])
    thetastart = 0.0
    thetaend = 2*pi
    points = circlexy(center,radius,refinement,thetastart,thetaend)
    # get gammas
    gammas = Array{Float64,1}[]
    for point in points
        rhat = (point - center)/la.norm(point-center)# radius unit vector
        # cross product results in a vector tangent to the circle
        # gamma scales to size
        Gamma = la.cross(rhat,[0,0,1.0])*gamma
        push!(gammas,Gamma)
    end
    # get velocities
    vs = Array{Float64,1}[]
    for point in points
        push!(vs,velocity)
    end
    particles = place(points,gammas,vs)

    return particles
end


text:
using sem, Feather, ModelingToolkit, Statistics, LinearAlgebra,
    Optim, SparseArrays, Test, Zygote, LineSearches, ForwardDiff,
    Distributions, FiniteDiff

## Observed Data
three_path_dat = Feather.read("test/comparisons/three_path_dat.feather")
three_path_par = Feather.read("test/comparisons/three_path_par.feather")
three_path_start = Feather.read("test/comparisons/three_path_start.feather")
fitmeasures = Feather.read("test/comparisons/three_path_fitm.feather")

semobserved = SemObsCommon(data = Matrix(three_path_dat))

diff_fin = SemFiniteDiff(BFGS(), Optim.Options())

## Model definition
@ModelingToolkit.variables x[1:31]

S =[x[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
    0     x[2]  0     0     0     0     0     0     0     0     0     0     0     0
    0     0     x[3]  0     0     0     0     0     0     0     0     0     0     0
    0     0     0     x[4]  0     0     0     x[15] 0     0     0     0     0     0
    0     0     0     0     x[5]  0     x[16] 0     x[17] 0     0     0     0     0
    0     0     0     0     0     x[6]  0     0     0     x[18] 0     0     0     0
    0     0     0     0     x[16] 0     x[7]  0     0     0     x[19] 0     0     0
    0     0     0     x[15] 0     0     0     x[8]  0     0     0     0     0     0
    0     0     0     0     x[17] 0     0     0     x[9]  0     x[20] 0     0     0
    0     0     0     0     0     x[18] 0     0     0     x[10] 0     0     0     0
    0     0     0     0     0     0     x[19] 0     x[20] 0     x[11] 0     0     0
    0     0     0     0     0     0     0     0     0     0     0     x[12] 0     0
    0     0     0     0     0     0     0     0     0     0     0     0     x[13] 0
    0     0     0     0     0     0     0     0     0     0     0     0     0     x[14]]

F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0]

A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[21] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[22] 0     0
    0  0  0  0  0  0  0  0  0  0  0     0     1     0
    0  0  0  0  0  0  0  0  0  0  0     0     x[23] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[24] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[25] 0
    0  0  0  0  0  0  0  0  0  0  0     0     0     1
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[26]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[27]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[28]
    0  0  0  0  0  0  0  0  0  0  0     0     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[29] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[30] x[31] 0]

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

start_val = vcat(
    vec(var(Matrix(three_path_dat), dims = 1))./2,
    fill(0.05, 3),
    fill(0.0, 6),
    fill(1.0, 8),
    fill(0, 3)
    )

loss = Loss([SemML(semobserved, [0.0], similar(start_val))])

imply = ImplySymbolic(A, S, F, x, start_val)

model_fin = Sem(semobserved, imply, loss, diff_fin)

@benchmark solution_fin = sem_fit(model_fin)

solution_fin.minimum

par_order = [collect(21:34); collect(15:20); 2;3; 5;6;7; collect(9:14)]

all(
    abs.(solution_fin.minimizer .- three_path_par.est[par_order]
        ) .< 0.001*abs.(three_path_par.est[par_order]))

model_fin.observed.n_man

size(solution_fin.minimizer, 1)

# bollen
function sem_chi2(model, solution)
    F = solution.minimum
    S = model.observed.obs_cov
    N = model.observed.n_obs
    n = model.observed.n_man
    q = Float64(size(solution.minimizer, 1))
    df = 0.5*n*(n+1) - q

    if isnothing(model.imply.imp_mean)
        χ = (N - 1.0)*(F - logdet(S) - n)
    end

    d = Chisq(df)
    p = 1 - cdf(d, χ)

    return χ, p
end

using BenchmarkTools

value, p = sem_chi2(model_fin, solution_fin)

value ≈ fitmeasures.chisq[1]
p ≈ fitmeasures.pvalue[1]

function sem_fi(model, solution)
    F = solution.minimum
    S = model.observed.obs_cov
    N = model.observed.n_obs
    n = model.observed.n_man
    q = Float64(size(solution.minimizer, 1))

    df = 0.5*n*(n+1) - q
    if isnothing(model.imply.imp_mean)
        χ = (N-1.0)*(F - logdet(S) - n)
    end

    df₀ = 0.5*n*(n-1)
    Σ₀= Diagonal(S)
    χ₀ = (N-1.0)*(tr(inv(Σ₀)*S) + logdet(Σ₀) - logdet(S) - n)

    FI = (χ₀ - df₀ - χ + df)/(χ₀ - df₀)

    return FI

end

sem_fi(model_fin, solution_fin) ≈ fitmeasures.cfi[1]

function sem_RMSEA(model, solution)
    F = solution.minimum
    S = model.observed.obs_cov
    N = model.observed.n_obs
    n = model.observed.n_man
    q = Float64(size(solution.minimizer, 1))

    df = 0.5*n*(n+1) - q
    if isnothing(model.imply.imp_mean)
        χ = (N-1.0)*(F - logdet(S) - n)
    end

    R = (χ - df)/((N-1.0)*df)
    R < 0.0 ? R = 0.0 : R = R
    R = √(R)

    return R
end

sem_RMSEA(model_fin, solution_fin) ≈ fitmeasures.rmsea[1]

function sem_logl(model, solution)
    let F = solution.minimum,
        N = model.observed.n_obs,
        n = model_fin.observed.n_man
            l = - N*0.5*(F + n*log(2pi))
    end
end

sem_logl(model_fin, solution_fin) ≈ fitmeasures.logl[1]

function logl(obs_mean, exp_cov, data_matr)
    exp_cov = Matrix(Hermitian(exp_cov))
    likelihood = -loglikelihood(MvNormal(obs_mean, exp_cov), transpose(data_matr))
    return likelihood
end

loglikelihood(MvNormal(vcat(Statistics.mean(Matrix(three_path_dat), dims = 1)...), model_fin.imply.imp_cov),transpose(model_fin.observed.data))

logl(vcat(Statistics.mean(Matrix(three_path_dat), dims = 1)...), model_fin.imply.imp_cov, model_fin.observed.data)

function sem_aic(model, solution)
    q = Float64(size(solution_fin.minimizer, 1))
    logl = sem_logl(model, solution)
    AIC = 2q - 2logl
    return AIC
end

function sem_aic_diff(model, solution)
    n = model.observed.n_man
    q = Float64(size(solution.minimizer, 1))
    df = 0.5*n*(n+1) - q

    χ = sem_chi2(model, solution)[1]
    AIC = χ - 2df
    return AIC
end

function sem_aic_bollen(model, solution)
    q = Float64(size(solution.minimizer, 1))
    χ = sem_chi2(model, solution)[1]
    AIC = χ + 2q
    return AIC
end

function sem_bic_andreas(model, solution)
    q = Float64(size(solution.minimizer, 1))
    χ = sem_chi2(model, solution)[1]
    N = model_fin.observed.n_obs

    BIC = χ + q*log(N)
    return BIC
end

function sem_bic(model, solution)
    q = Float64(size(solution_fin.minimizer, 1))
    logl = sem_logl(model, solution)

    BIC = q*log(N) - 2logl
    return BIC
end

# tests only work if lavaan meanstructure = FALSE

sem_aic(model_fin, solution_fin) ≈ fitmeasures.aic[1]

sem_bic(model_fin, solution_fin)

#

fitmeasures |> names

function f(func)
    getproperty(LinearAlgebra, Symbol(func))
end

f(inv)

getproperty(LinearAlgebra, Symbol("inv"))(rand(5,5))

v = inv

s = Symbol(v)

getproperty(LinearAlgebra, s)

pre = zeros(31,31)

FiniteDiff.finite_difference_hessian!(pre, model_fin, solution_fin.minimizer)

# hessian! from NLSolversBase is faster if the objective already was
# TwiceDifferentiable, otherwise it is unclear
function sem_se(td::T, solution) where {T <: TwiceDifferentiable}
    H = hessian!(td, solution.minimizer)
    return H
end

function sem_se(model, solution)
    H = FiniteDiff.finite_difference_hessian(model, solution.minimizer)
    return H
end

sem_se(func, solution_fin) ≈ pre

using NLSolversBase

func = TwiceDifferentiable(model_fin, start_val)

obj = OnceDifferentiable(model_fin, start_val)

stored_state = Optim.initial_state(model_fin.diff.algorithm, model_fin.diff.options,
    obj, model_fin.imply.start_val)

result = optimize(
    obj,
    model_fin.imply.start_val,
    model_fin.diff.algorithm,
    model_fin.diff.options,
    stored_state)

@btime f()

diff = stored_state.invH - inv(pre)

maximum(diff)

se1 = .√diag(inv(pre))
se2 = .√diag(stored_state.invH)

diff = se1 - se2

maximum(diff)

z = solution_fin.minimizer / se1

all(
    abs.(se1 .- three_path_par.se[par_order]
        ) .< 0.001*abs.(three_path_par.se[par_order]))

all(
    abs.(z .- three_path_par.z[par_order]
        ) .< 0.001*abs.(three_path_par.z[par_order]))

using LinearAlgebra, BenchmarkTools

x = rand(50)

y = rand(50)

A = rand(50,50)
A = A*A'
A = cholesky(A)
x = permutedims(x)

@benchmark dot($x, $A, $y)

x'*A*y

dot(x, A, y)

function mydot(x, A, y, pre)
    mul!(pre, x, A)
    return dot(pre, y)
end

pre = zeros(50,1)
mydot(x, A, y, pre)