text:
module AbstractGaugefields_module
    using LinearAlgebra
    import ..Wilsonloops_module:Wilson_loop_set,calc_coordinate,make_plaq_staple_prime,calc_shift,make_plaq,make_plaq_staple,
                        Tensor_wilson_lines_set,Tensor_wilson_lines,Tensor_derivative_set,
                        get_leftstartposition,get_rightstartposition,Wilson_loop,calc_loopset_μν_name
    import Wilsonloop:loops_staple_prime,Wilsonline,get_position,get_direction,GLink,isdag
    using Requires
    using Distributions
    using StableRNGs
    import ..Verboseprint_mpi:Verbose_print,println_verbose_level1,println_verbose_level2,println_verbose_level3

    #using MPI
    using InteractiveUtils

    abstract type Abstractfields end

    abstract type AbstractGaugefields{NC,Dim} <: Abstractfields
    end

    abstract type Adjoint_fields{T} <: Abstractfields
    end

    struct Adjoint_Gaugefields{T} <: Adjoint_fields{T}
        parent::T
    end

    function Base.adjoint(U::T) where T <: Abstractfields
        Adjoint_Gaugefields{T}(U)
    end

    abstract type Unit_Gaugefield{NC,Dim} <: Abstractfields
    end

    function Base.adjoint(U::Unit_Gaugefield) where T <: Unit_Gaugefield
        return U
    end

    abstract type Shifted_Gaugefields{NC,Dim} <: Abstractfields
    end

    struct Staggered_Gaugefields{T,μ} <: Abstractfields
        parent::T

        function Staggered_Gaugefields(u::T,μ) where T <: Abstractfields
            return new{T,μ}(u)
        end
    end

    function staggered_U(u::T,μ) where T <: Abstractfields
        return Staggered_Gaugefields(u,μ)
    end

    function Base.setindex!(U::T,v...)  where T <: Staggered_Gaugefields
        error("type $(typeof(U)) has no setindex method. This type is read only.")
    end

    function Base.setindex!(U::T,v...)  where T <:  Unit_Gaugefield
        error("type $(typeof(U)) has no setindex method. This type is read only.")
    end

    function Base.getindex(U::T,v...)  where T <:  Staggered_Gaugefields
        error("type $(typeof(U)) has no getindex method")
    end

    struct Gaugefield_latticeindices{Dim,NC,T}
        NN::NTuple{Dim,Int64}
        NC::Int8

        function Gaugefield_latticeindices(u::AbstractGaugefields{Dim,NC}) where {Dim,NC}
            _,_,NN... = size(u)
            return new{Dim,NC,typeof{u}}(NN,NC)
        end
    end

    function Gaugefield_latticeindices(U::Array{<: AbstractGaugefields{Dim,NC},1}) where {Dim,NC}
        return Gaugefield_latticeindices(U[1])
    end

    function Base.iterate(g::Gaugefield_latticeindices{Dim,NC,T}) where {Dim,NC,T}
        N<1 && return nothing
        return (1, 2)
    end

    mutable struct Data_sent{NC} #data format for MPI
        count::Int64
        data::Array{ComplexF64,3}
        positions::Vector{Int64}

        function Data_sent(N,NC)
            data = zeros(ComplexF64,NC,NC,N)
            count = 0
            positions = zeros(Int64,N)

            return new{NC}(count,data,positions)
        end
    end

    include("./2D/gaugefields_2D.jl")
    include("./4D/gaugefields_4D.jl")

    include("TA_Gaugefields.jl")
    include("Adjoint_rep_Gaugefields.jl")

    function LinearAlgebra.mul!(C,A::T,B) where T <: Unit_Gaugefield
        substitute_U!(C,B)
    end

    function LinearAlgebra.mul!(C,A,B::T) where T <: Unit_Gaugefield
        substitute_U!(C,A)
    end

    function Staggered_Gaugefields(u::AbstractGaugefields{NC,Dim}) where {NC,Dim}
        if Dim == 4
            return Staggered_Gaugefields_4D(u)
        else
            error("Dim = $Dim is not supported")
        end
    end

    function println_verbose_level1(u::T,val...)  where T <: AbstractGaugefields
        println_verbose_level1(u.verbose_print,val...)
    end

    function println_verbose_level2(u::T,val...)  where T <: AbstractGaugefields
        println_verbose_level2(u.verbose_print,val...)
    end

    function println_verbose_level3(u::T,val...)  where T <: AbstractGaugefields
        println_verbose_level3(u.verbose_print,val...)
    end

    function get_myrank(U::T) where T <: AbstractGaugefields
        return 0
    end

    function get_myrank(U::Array{T,1}) where T <: AbstractGaugefields
        return 0
    end

    function get_nprocs(U::T) where T <: AbstractGaugefields
        return 1
    end

    function get_nprocs(U::Array{T,1}) where T <: AbstractGaugefields
        return 1
    end

    function barrier(x::T) where T <: AbstractGaugefields
        return
    end

    function getvalue(U::T,i1,i2,i3,i4,i5,i6) where T <: Abstractfields
        return U[i1,i2,i3,i4,i5,i6]
    end

    function getvalue(U::T,i1,i2,i3,i6) where T <: Abstractfields
        return U[i1,i2,i3,i6]
    end

    #include("./gaugefields_4D_wing_mpi.jl")

    using NPZ

    function write_to_numpyarray(U::T,filename) where T <: AbstractGaugefields
        error("write_to_numpyarray! is not implemented in type $(typeof(U)) ")
    end

    function Base.similar(U::T) where T <: AbstractGaugefields
        error("similar! is not implemented in type $(typeof(U)) ")
    end

    function substitute_U!(a::Array{<: AbstractGaugefields{NC,Dim},1},b::Array{<: AbstractGaugefields{NC,Dim},1}) where {NC,Dim}
        error("substitute_U! is not implemented in type $(typeof(a)) and $(typeof(b))")
        for i=1:Dim
            substitute_U!(a[i],b[i])
        end
    end

    function substitute_U!(a::T1,b::T2) where {T1 <: Abstractfields,T2 <: Abstractfields}
        error("substitute_U! is not implemented in type $(typeof(a)) and $(typeof(b))")
        return
    end

    function substitute_U!(a::Array{T1,1},b::Array{T2,1},iseven::Bool) where {T1 <: AbstractGaugefields,T2 <: AbstractGaugefields}
        error("substitute_U! is not implemented in type $(typeof(a)) and $(typeof(b))")
    end

    function substitute_U!(a::T1,b::T2,iseven::Bool) where {T1 <: Abstractfields,T2 <: Abstractfields}
        error("substitute_U!(a,b,iseven) is not implemented in type $(typeof(a)) and $(typeof(b))")
        return
    end

    function Initialize_4DGaugefields(NC,NDW,NN...;condition = "cold",verbose_level =2,randomnumber="Random")
        if condition == "cold"
            if NDW == 0

                u1 = IdentityGauges_4D(NC,NN...,verbose_level = verbose_level)
            else
                u1 = IdentityGauges_4D(NC,NDW,NN...,verbose_level = verbose_level)
            end
        elseif condition == "hot"
            if NDW == 0
                u1 = RandomGauges_4D(NC,NN...,verbose_level = verbose_level,randomnumber=randomnumber)
            else
                u1 = RandomGauges_4D(NC,NDW,NN...,verbose_level = verbose_level,randomnumber=randomnumber)
            end
        else
            error("not supported")
        end
        @assert length(NN) == 4 "Dimension should be 4. "
        Dim = 4

        U = Array{typeof(u1),1}(undef,Dim)
        U[1] = u1

        for μ=2:Dim
            if condition == "cold"
                if NDW == 0
                    U[μ] = IdentityGauges_4D(NC,NN...,verbose_level = verbose_level)
                else
                    U[μ] = IdentityGauges_4D(NC,NDW,NN...,verbose_level = verbose_level)
                end
            elseif condition == "hot"
                if NDW == 0
                    U[μ] = RandomGauges_4D(NC,NN...,verbose_level = verbose_level,randomnumber=randomnumber)
                else
                    U[μ] = RandomGauges_4D(NC,NDW,NN...,verbose_level = verbose_level,randomnumber=randomnumber)
                end
            else
                error("not supported")
            end
        end
        return U
    end
"""
```
Initialize_Gaugefields(NC,NDW,NN...;
    condition = "cold",mpi = false,PEs=nothing,mpiinit = nothing,verbose_level = 2,randomnumber="Random")
```
Initialize gaugefields. 2D or 4D is supported.

### Examples
If you want to have randomely distributed gauge fields (so-called "hot start") in four dimension, just do:

```julia
U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "hot")
```

If you want to have uniform gauge fields (so-called "cold start") in four dimension, just do:

```julia
U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")
```

### Return values:
*   `U`: Dim dimensional vector. Here, Dim is a dimension (Dim = 2 or 4).

### Keyword arguments:
* `condition`: ```cold```(cold start) or ```hot```(hot start). The default is ```cold```
* `mpi`: Using MPI or not. The default is ```false```. If you want to use MPI, you should call ```using MPI```.
* `PEs`: If ```mpi = true```, we have to set ```PEs = [2,2,2,2]```, which is numbers of regions for MPI process. ```prod(PEs)``` should be the number of total MPI processes.

### MPI
Gaugefields with using MPI is not well tested.
"""
    function Initialize_Gaugefields(NC,NDW,NN...;condition = "cold",mpi = false,PEs=nothing,mpiinit = nothing,verbose_level = 2,randomnumber="Random")
        Dim = length(NN)
        if condition == "cold"
            u1 = IdentityGauges(NC,NDW,NN...,mpi =mpi, PEs=PEs, mpiinit = mpiinit,verbose_level =verbose_level)
        elseif condition == "hot"
            u1 = RandomGauges(NC,NDW,NN...,mpi =mpi, PEs=PEs, mpiinit = mpiinit,verbose_level =verbose_level,randomnumber=randomnumber)
        else
            error("not supported")
        end

        U = Array{typeof(u1),1}(undef,Dim)
        U[1] = u1

        for μ=2:Dim
            if condition == "cold"
                U[μ] = IdentityGauges(NC,NDW,NN...,mpi =mpi, PEs=PEs, mpiinit = false,verbose_level =verbose_level)
            elseif condition == "hot"
                U[μ] = RandomGauges(NC,NDW,NN...,mpi =mpi, PEs=PEs, mpiinit = false,verbose_level =verbose_level,randomnumber=randomnumber)
            else
                error("not supported")
            end
        end
        return U
    end

    function RandomGauges(NC,NDW,NN...;mpi = false,PEs=nothing,mpiinit = nothing,verbose_level =2,randomnumber="Random")
        dim = length(NN)
        if mpi
            if PEs == nothing || mpiinit == nothing
                error("not implemented yet!")
            else
                if dim==4
                    if NDW == 0
                        U = randomGaugefields_4D_nowing_mpi(NC,NN[1],NN[2],NN[3],NN[4],PEs,mpiinit = mpiinit,verbose_level = verbose_level,randomnumber=randomnumber)
                    else
                        U = randomGaugefields_4D_wing_mpi(NC,NN[1],NN[2],NN[3],NN[4],NDW,PEs,mpiinit = mpiinit,verbose_level = verbose_level,randomnumber=randomnumber)
                    end
                else
                    error("not implemented yet!")
                end

            end
        else
            if dim == 4
                if NDW == 0
                    U = randomGaugefields_4D_nowing(NC,NN[1],NN[2],NN[3],NN[4],verbose_level= verbose_level,randomnumber=randomnumber)

                else
                    U = randomGaugefields_4D_wing(NC,NN[1],NN[2],NN[3],NN[4],NDW,verbose_level= verbose_level,randomnumber=randomnumber)
                end
            elseif dim == 2
                if NDW == 0
                    U = randomGaugefields_2D_nowing(NC,NN[1],NN[2],verbose_level  = verbose_level,randomnumber=randomnumber)
                else
                    U = randomGaugefields_2D_wing(NC,NN[1],NN[2],NDW,verbose_level  = verbose_level,randomnumber=randomnumber)
                end
            else
                error("not implemented yet!")
            end
        end
        return U
    end

    function IdentityGauges(NC,NDW,NN...;mpi = false,PEs=nothing,mpiinit = nothing,verbose_level =2)
        dim = length(NN)

        if mpi
            if PEs == nothing || mpiinit == nothing
                error("not implemented yet!")
            else
                if dim==4
                    if NDW == 0
                        U = identityGaugefields_4D_nowing_mpi(NC,NN[1],NN[2],NN[3],NN[4],PEs,mpiinit = mpiinit,verbose_level=verbose_level)
                    else
                        U = identityGaugefields_4D_wing_mpi(NC,NN[1],NN[2],NN[3],NN[4],NDW,PEs,mpiinit = mpiinit,verbose_level=verbose_level)
                    end
                elseif dim == 2
                    if NDW == 0
                        U = identityGaugefields_2D_nowing_mpi(NC,NN[1],NN[2],PEs,mpiinit = mpiinit,verbose_level=verbose_level)
                    end
                else
                    error("$dim dimension with $NDW  is not implemented yet! set NDW = 0")
                end

            end
        else
            if dim == 4
                if NDW == 0
                    U = identityGaugefields_4D_nowing(NC,NN[1],NN[2],NN[3],NN[4],verbose_level =verbose_level)
                else
                    U = identityGaugefields_4D_wing(NC,NN[1],NN[2],NN[3],NN[4],NDW,verbose_level =verbose_level)
                end
            elseif dim == 2
                if NDW == 0

                    U = identityGaugefields_2D_nowing(NC,NN[1],NN[2],verbose_level = verbose_level)
                else
                    U = identityGaugefields_2D_wing(NC,NN[1],NN[2],NDW,verbose_level = verbose_level)
                end
            else
                error("$dim dimension system is not implemented yet!")
            end
        end
        set_wing_U!(U)
        return U
    end

    function Oneinstanton(NC,NDW,NN...;mpi = false,PEs=nothing,mpiinit = nothing)
        dim = length(NN)
        if mpi
            if PEs == nothing || mpiinit == nothing
                error("not implemented yet!")
            else
                if dim==4

                    U = Oneinstanton_4D_wing_mpi(NC,NN[1],NN[2],NN[3],NN[4],NDW,PEs,mpiinit = mpiinit)
                else
                    error("not implemented yet!")
                end

            end
        else
            if dim == 4
                if NDW == 0
                    U = Oneinstanton_4D_nowing(NC,NN[1],NN[2],NN[3],NN[4])
                else
                    U = Oneinstanton_4D_wing(NC,NN[1],NN[2],NN[3],NN[4],NDW)
                end
            elseif dim == 2
                if NDW == 0
                    U = Oneinstanton_2D_nowing(NC,NN[1],NN[2])
                else
                    U = Oneinstanton_2D_wing(NC,NN[1],NN[2],NDW)

                end
            else
                error("not implemented yet!")
            end
        end
        return U
    end

    function construct_gauges(NC,NDW,NN...;mpi = false,PEs=nothing,mpiinit = nothing)
        dim = length(NN)
        if mpi
            if PEs == nothing || mpiinit == nothing
                error("not implemented yet!")
            else
                if dim==4
                    U = identityGaugefields_4D_wing_mpi(NC,NN[1],NN[2],NN[3],NN[4],NDW,PEs,mpiinit = mpiinit)
                else
                    error("not implemented yet!")
                end

            end
        else
            if dim == 4
                if NDW == 0
                    U = identityGaugefields_4D_nowing(NC,NN[1],NN[2],NN[3],NN[4])
                else
                    U = identityGaugefields_4D_wing(NC,NN[1],NN[2],NN[3],NN[4],NDW)
                end
            elseif dim == 2
                if NDW == 0
                    U = identityGaugefields_2D_nowing(NC,NN[1],NN[4])
                else
                    U = identityGaugefields_2D_wing(NC,NN[1],NN[4],NDW)
                end
            else
                error("not implemented yet!")
            end
        end
        return U
    end

    function clear_U!(U::T) where T <: AbstractGaugefields
        error("clear_U! is not implemented in type $(typeof(U)) ")
    end

    function clear_U!(U::Array{<: AbstractGaugefields{NC,Dim},1}) where {NC,Dim}
        for μ=1:Dim
            clear_U!(U[μ])
        end
    end

    function clear_U!(U::T,iseven::Bool) where T <: AbstractGaugefields
        error("clear_U! is not implemented in type $(typeof(U)) ")
    end

    function clear_U!(U::Array{<: AbstractGaugefields{NC,Dim},1},iseven::Bool) where {NC,Dim}
        for μ=1:Dim
            clear_U!(U[μ],iseven)
        end
    end

    function unit_U!(U::T) where T <: AbstractGaugefields
        error("unit_U! is not implemented in type $(typeof(U)) ")
    end

    function unit_U!(U::Array{<: AbstractGaugefields{NC,Dim},1}) where {NC,Dim}
        for μ=1:Dim
            unit_U!(U[μ])
        end
    end

    function shift_U(U::AbstractGaugefields{NC,Dim},ν) where {NC,Dim}
        error("shift_U is not implemented in type $(typeof(U)) ")
        return nothing
    end

    function map_U!(U::AbstractGaugefields{NC,Dim},f::Function,V::AbstractGaugefields{NC,Dim},iseven::Bool) where {NC,Dim, T<: Abstractfields}
        error("map_U! is not implemented in type $(typeof(U)) ")
        return nothing
    end

    function map_U!(U::AbstractGaugefields{NC,Dim},f::Function,V::AbstractGaugefields{NC,Dim}) where {NC,Dim, T<: Abstractfields}
        error("map_U! is not implemented in type $(typeof(U)) ")
        return nothing
    end

    function map_U_sequential!(U::AbstractGaugefields{NC,Dim},f::Function,Uin) where {NC,Dim, T<: Abstractfields}
        error("map_U_sequential! is not implemented in type $(typeof(U)) ")
        return nothing
    end

    function identitymatrix(U::T) where T <: AbstractGaugefields
        error("identitymatrix is not implemented in type $(typeof(U)) ")
    end

    function set_wing_U!(U::Array{<: AbstractGaugefields{NC,Dim},1}) where {NC,Dim}
        for μ=1:Dim
            set_wing_U!(U[μ])
        end
    end

    function set_wing_U!(U::T) where T <: AbstractGaugefields
        error("set_wing_U! is not implemented in type $(typeof(U)) ")
    end

    function evaluate_gaugelinks_evenodd!(uout::T,w::Wilsonline{Dim},U::Array{T,1},temps::Array{T,1},iseven) where {T<: AbstractGaugefields,Dim}

        #Uold = temps[1]
        Unew = temps[1]
        #Utemp2 = temps[2]
        #clear_U!(uout)
        origin = Tuple(zeros(Int64,Dim))

        glinks = w
        numlinks = length(glinks)
        if numlinks == 0
            unit_U!(uout)
            return
        end

        j = 1
        U1link = glinks[1]
        direction = get_direction(U1link)
        position = get_position(U1link)
        isU1dag = isdag(U1link)
        #isU1dag = ifelse(typeof(U1link) <: Adjoint_GLink,true,false)

        if numlinks == 1
            substitute_U!(Unew,U[direction])
            Ushift1 = shift_U(Unew,position)
            if isU1dag
                #println("Ushift1 ",Ushift1'[1,1,1,1,1,1])
                substitute_U!(uout,Ushift1',iseven)
            else
                substitute_U!(uout,Ushift1,iseven)
            end
            return
        end

        substitute_U!(Unew,U[direction])
        #println(position)
        Ushift1 = shift_U(Unew,position)

        for j=2:numlinks
            Ujlink = glinks[j]
            isUkdag = isdag(Ujlink)
            #isUkdag = ifelse(typeof(Ujlink) <: Adjoint_GLink,true,false)
            position = get_position(Ujlink)
            direction = get_direction(Ujlink)

            Ushift2 = shift_U(U[direction],position)
            multiply_12!(uout,Ushift1,Ushift2,j,isUkdag,isU1dag,iseven)

            substitute_U!(Unew,uout)

            Ushift1 = shift_U(Unew,origin)
        end

    end

    function evaluate_gaugelinks!(uout::T,w::Wilsonline{Dim},U::Array{T,1},temps::Array{T,1}) where {T<: AbstractGaugefields,Dim}
        #println_verbose_level3(uout,"evaluating Wilson loops")
        #Uold = temps[1]
        #set_wing_U!(U)
        Unew = temps[1]
        #Utemp2 = temps[2]
        #clear_U!(uout)
        origin = Tuple(zeros(Int64,Dim))

        glinks = w
        numlinks = length(glinks)
        if numlinks == 0
            unit_U!(uout)
            return
        end

        j = 1
        U1link = glinks[1]
        direction = get_direction(U1link)
        position = get_position(U1link)
        isU1dag = isdag(U1link)
        #isU1dag = ifelse(typeof(U1link) <: Adjoint_GLink,true,false)

        #show(glinks)
        #println("in evaluate_gaugelinks!")
        #show(w)
        #println("numlinks = $numlinks")
        if numlinks == 1
            substitute_U!(Unew,U[direction])
            Ushift1 = shift_U(Unew,position)
            if isU1dag
                #println("Ushift1 ",Ushift1'[1,1,1,1,1,1])
                substitute_U!(uout,Ushift1')
            else
                substitute_U!(uout,Ushift1)
            end

            return
        end

        #j = 1
        #U1link = glinks[1]
        #direction = get_direction(U1link)
        #position = get_position(U1link)
        #println("i = $i j = $j position = $position")
        substitute_U!(Unew,U[direction])
        Ushift1 = shift_U(Unew,position)

        #ix,iy,iz,it=(2,2,2,2)
        #println("position = $position")
        #pos = Tuple([ix,iy,iz,it] .+ collect(position))
        #U1 = Unew[:,:,pos...]
        #println("U1, ",Unew[:,:,pos...])
        #isU1dag = U1link.isdag
        #isU1dag = ifelse(typeof(U1link) <: Adjoint_GLink,true,false)

        for j=2:numlinks
            Ujlink = glinks[j]
            isUkdag = isdag(Ujlink)
            #isUkdag = ifelse(typeof(Ujlink) <: Adjoint_GLink,true,false)
            position = get_position(Ujlink)
            direction = get_direction(Ujlink)
            #println("j = $j position = $position")
            #println("a,b, $isUkdag , $isU1dag")
            Ushift2 = shift_U(U[direction],position)

            #zerocheck(U,Ushift1.parent.Ushifted,"Ushift1")
            #zerocheck(U,Ushift2.parent.Ushifted,"Ushift2 position $position")

            multiply_12!(uout,Ushift1,Ushift2,j,isUkdag,isU1dag)

            #zerocheck(U,uout.U,"uout")

            #pos = Tuple([ix,iy,iz,it] .+ collect(position))
            #U2 = U[direction][:,:,pos...]
            #println("U1U2dag ", U1*U2')
            substitute_U!(Unew,uout)

            #println("Unew ", Unew[:,:,ix,iy,iz,it])

            Ushift1 = shift_U(Unew,origin)
            #println("uout ", uout[:,:,ix,iy,iz,it])
        end

        #zerocheck(U,uout.U,"uoutfinal")

        #println("uout2 ", uout[:,:,ix,iy,iz,it])

    end

    function zerocheck(U,U2bare,Uname)
        for myrank=0:(get_nprocs(U)-1)
            #println(get_nprocs(U))
            if get_myrank(U) == myrank

                #println("rank = $myrank Ushift1",Ushift1.parent.Ushifted)
                #println("rank = $myrank Ushift2",Ushift2.parent.Ushifted)

                nc,_,nx,ny,nz,nt = size(U2bare)
                for it=1:nt
                    for iz=1:nz
                        for iy=1:ny
                            for ix = 1:nx
                                ic = 1
                                jc = 1
                                #for ic=1:nc
                                #    for jc = 1:nc
                                        v = U2bare[jc,ic,ix,iy,iz,it]
                                        if abs(v) <= 1e-15
                                            println("rank = $myrank ",Uname)
                                            println("$position type: $(typeof(U))\n ",v," at $((jc,ic,ix,iy,iz,it))")
                                        end

                                #    end
                                #end
                            end
                        end
                    end
                end

                #println("rank $myrank j = $j $position Ushift1 $(typeof(Ushift1)) ",getvalue(Ushift1,1,1,1,1,1,1))
                #println("rank $myrank j = $j $position Ushift2 $(typeof(Ushift2)) ",getvalue(Ushift2,1,1,1,1,1,1))
            end
            barrier(U[1])
        end
    end
    function evaluate_gaugelinks_evenodd!(xout::T,w::Array{<:Wilsonline{Dim},1},U::Array{T,1},temps::Array{T,1},iseven) where {T<: AbstractGaugefields,Dim}
        num = length(w)
        temp1 = temps[1]
        temp2 = temps[2]

        #ix,iy,iz,it=(2,2,2,2)

        clear_U!(xout,iseven)
        for i=1:num
            glinks = w[i]
            evaluate_gaugelinks_evenodd!(temp2,glinks,U,[temp1],iseven)
            #println("uout2 ", temp2[:,:,ix,iy,iz,it])
            add_U!(xout,temp2,iseven)
            #println("xout ", xout[:,:,ix,iy,iz,it])
        end

        #println("xout2 ", xout[:,:,ix,iy,iz,it])
        return

    end

    function evaluate_gaugelinks!(xout::T,w::Array{<:Wilsonline{Dim},1},U::Array{T,1},temps::Array{T,1}) where {T<: AbstractGaugefields,Dim}
        num = length(w)
        temp1 = temps[1]
        temp2 = temps[2]

        #ix,iy,iz,it=(2,2,2,2)
        #ix,iy,iz,it=(1,1,1,1)

        clear_U!(xout)
        for i=1:num
            glinks = w[i]
            evaluate_gaugelinks!(temp2,glinks,U,[temp1])
            #println("uout2 ", temp2[:,:,ix,iy,iz,it])
            add_U!(xout,temp2)
            #println("xout ", xout[:,:,ix,iy,iz,it])
        end

        #println("xout2 ", xout[:,:,ix,iy,iz,it])
        return

    end

    #=
    function evaluate_gaugelinks_inside!(U,glinks,Ushift1,Uold,Unew,numlinks)
        for k=2:numlinks
            position =get_position(glinks[k])
            direction = get_direction(glinks[k])
            Ushift2 = shift_U(U[direction],position)
            multiply_12!(Unew,Ushift1,Ushift2,k,loopk,loopk1_2)

            Unew,Uold = Uold,Unew
            Ushift1 = shift_U(Uold,(0,0,0,0))
        end
    end
    =#

    function multiply_12!(temp3,temp1,temp2,k,isUkdag::Bool,isU1dag::Bool)
        if k==2
            if isUkdag
                if isU1dag
                    mul!(temp3,temp1',temp2')
                else
                    mul!(temp3,temp1,temp2')
                end
            else
                if isU1dag
                    mul!(temp3,temp1',temp2)
                else
                    mul!(temp3,temp1,temp2)
                end
            end
        else
            if isUkdag
                mul!(temp3,temp1,temp2')
            else
                mul!(temp3,temp1,temp2)
            end
        end
        return
    end

    function multiply_12!(temp3,temp1,temp2,k,isUkdag::Bool,isU1dag::Bool,iseven)
        if k==2
            if isUkdag
                if isU1dag
                    mul!(temp3,temp1',temp2',iseven)
                else
                    mul!(temp3,temp1,temp2',iseven)
                end
            else
                if isU1dag
                    mul!(temp3,temp1',temp2,iseven)
                else
                    mul!(temp3,temp1,temp2,iseven)
                end
            end
        else
            if isUkdag
                mul!(temp3,temp1,temp2',iseven)
            else
                mul!(temp3,temp1,temp2,iseven)
            end
        end
        return
    end

    function evaluate_gaugelinks_eachsite!(uout_mat::T1,w::Array{<:Wilsonline{Dim},1},U::Array{<:AbstractGaugefields{NC,Dim},1},
        temps,indices...) where {Dim,T1 <: Matrix{ComplexF64},NC}

        #xout::T,w::Array{<:Wilsonline{Dim},1},U::Array{T,1},temps::Array{T,1},iseven) where {T<: AbstractGaugefields,Dim}
        num = length(w)
        temp = temps[4]

        #ix,iy,iz,it=(2,2,2,2)

        uout_mat .= 0
        for i=1:num
            glinks = w[i]
            evaluate_gaugelinks_eachsite!(temp,glinks,U,view(temps,1:3),indices...)

            #println("uout2 ", temp2[:,:,ix,iy,iz,it])
            uout_mat .+= temp
            #println("xout ", xout[:,:,ix,iy,iz,it])
        end

        #println("xout2 ", xout[:,:,ix,iy,iz,it])
        return

    end

    function evaluate_gaugelinks_eachsite!(uout_mat::T1,w::Wilsonline{Dim},U::Array{<:AbstractGaugefields{NC,Dim},1},
                            temps,indices...) where {Dim,T1 <: Matrix{ComplexF64},NC}
        _,_,NN... = size(U[1])
        Unew = temps[1]
        Ushift1 = temps[2]
        Ushift2 = temps[3]

        #origin = Tuple(zeros(Int64,Dim))
        shifted = zeros(Int64,Dim)

        glinks = w
        numlinks = length(glinks)
        if numlinks == 0
            for i=1:NC
                for j=1:NC
                    uout_mat[j,i] = ifelse(i==j,1,0)
                end
            end
            return
        end

        j = 1
        U1link = glinks[1]
        direction = get_direction(U1link)
        position = get_position(U1link)
        isU1dag = isdag(U1link)#.isdag
        #isU1dag = ifelse(typeof(U1link) <: Adjoint_GLink,true,false)

        #show(glinks)
        #println("in evaluate_gaugelinks!")
        #show(w)
        #println("numlinks = $numlinks")
        if numlinks == 1
            for i=1:NC
                for j=1:NC
                    Unew[j,i] = U[direction][j,i,indices...]
                end
            end

            #substitute_U!(Unew,U[direction])
            for μ=1:Dim
                shifted[μ] = indices[μ] + position[μ]# collect(position)
            end
            #shifted = collect(indices) .+ collect(position)
            for μ=1:Dim
                shifted[μ] += ifelse(shifted[μ] < 1,NN[μ],0)
                shifted[μ] += ifelse(shifted[μ] > NN[μ],-NN[μ],0)
            end
            #shiftedposition = Tuple(shifted)
            #shiftedposition = Tuple(collect(indices) .+ collect(position))
            for i=1:NC
                for j=1:NC
                    #Ushift1[j,i] = Unew[j,i,shiftedposition...]
                    Ushift1[j,i] = Unew[j,i,shifted...]
                end
            end

            #Ushift1 = shift_U(Unew,position)
            if isU1dag
                #println("Ushift1 ",Ushift1'[1,1,1,1,1,1])
                for i=1:NC
                    for j=1:NC
                        Unew[j,i] = Ushift1'[j,i]#,indices]
                    end
                end
                #substitute_U!(uout,Ushift1')
            else
                for i=1:NC
                    for j=1:NC
                        Unew[j,i] = Ushift1[j,i]#,indices]
                    end
                end
                #substitute_U!(uout,Ushift1)
            end
            return
        end

        #j = 1
        #U1link = glinks[1]
        #direction = get_direction(U1link)
        #position = get_position(U1link)
        #println("i = $i j = $j position = $position")
        #substitute_U!(Unew,U[direction])
        #for b=1:NC
        #    for a=1:NC
        #        Unew[a,b] = U[direction][a,b,indices...]
        #    end
        #end

        for μ=1:Dim
            shifted[μ] = indices[μ] + position[μ]# collect(position)
        end
        #shifted = collect(indices) .+ collect(position)
        for μ=1:Dim
            shifted[μ] += ifelse(shifted[μ] < 1,NN[μ],0)
            shifted[μ] += ifelse(shifted[μ] > NN[μ],-NN[μ],0)
        end
        #shiftedposition = Tuple(shifted)
        for i=1:NC
            for j=1:NC
                #Ushift1[j,i] = U[direction][j,i,shiftedposition...]
                Ushift1[j,i] = U[direction][j,i,shifted...]
            end
        end

        #Ushift1 = shift_U(Unew,position)

        #ix,iy,iz,it=(2,2,2,2)
        #println("posotion = $position")
        #pos = Tuple([ix,iy,iz,it] .+ collect(position))
        #U1 = Unew[:,:,pos...]
        #println("U1, ",Unew[:,:,pos...])
        #isU1dag = ifelse(typeof(U1link) <: Adjoint_GLink,true,false)

        for j=2:numlinks
            Ujlink = glinks[j]
            isUkdag = isdag(Ujlink)
            #isUkdag = ifelse(typeof(Ujlink) <: Adjoint_GLink,true,false)
            position = get_position(Ujlink)
            direction = get_direction(Ujlink)
            #println("j = $j position = $position")
            #println("a,b, $isUkdag , $isU1dag")
            for μ=1:Dim
                shifted[μ] = indices[μ] + position[μ]# collect(position)
            end
            #shifted = collect(indices) .+ collect(position)
            for μ=1:Dim
                shifted[μ] += ifelse(shifted[μ] < 1,NN[μ],0)
                shifted[μ] += ifelse(shifted[μ] > NN[μ],-NN[μ],0)
            end
            #shiftedposition2 = Tuple(shifted)
            ##  shiftedposition2 = Tuple(collect(indices) .+ collect(position))
            for b=1:NC
                for a=1:NC
                    #Ushift2[a,b] = U[direction][a,b,shiftedposition2...]
                    Ushift2[a,b] = U[direction][a,b,shifted...]
                end
            end

            #Ushift2 = shift_U(U[direction],position)
            multiply_12!(uout_mat,Ushift1,Ushift2,j,isUkdag,isU1dag)

            #pos = Tuple([ix,iy,iz,it] .+ collect(position))
            #U2 = U[direction][:,:,pos...]
            #println("U1U2dag ", U1*U2')
            #substitute_U!(Unew,uout)
            for b=1:NC
                for a=1:NC
                    Unew[a,b] =uout_mat[a,b]
                end
            end

            #println("Unew ", Unew[:,:,ix,iy,iz,it])
            for b=1:NC
                for a=1:NC
                    Ushift1[a,b] = Unew[a,b]
                end
            end

            #Ushift1 = shift_U(Unew,origin)
            #println("uout ", uout[:,:,ix,iy,iz,it])
        end

        #println("uout2 ", uout[:,:,ix,iy,iz,it])

    end

    const shifted_4_temp = zeros(Int64,4)
    const NN_4_temp = zeros(Int64,4)

    function evaluate_gaugelinks_eachsite!(uout_mat::T1,w::Wilsonline{4},U::Array{<:AbstractGaugefields{NC,4},1},
                            temps,indices...) where {T1 <: Matrix{ComplexF64},NC}
        #NN = NN_4_temp
        #NN[1] = U[1].NX
        #NN[2] = U[1].NY
        #NN[3] = U[1].NZ
        #NN[4] = U[1].NT
        _,_,NN... = size(U[1])

        Unew = temps[1]
        Ushift1 = temps[2]
        Ushift2 = temps[3]

        #origin = Tuple(zeros(Int64,Dim))
        shifted = shifted_4_temp #zeros(Int64,Dim)

        glinks = w
        numlinks = length(glinks)
        if numlinks == 0
            for i=1:NC
                for j=1:NC
                    uout_mat[j,i] = ifelse(i==j,1,0)
                end
            end
            return
        end

        j = 1
        U1link = glinks[1]
        direction = get_direction(U1link)
        position = get_position(U1link)
        isU1dag = isdag(U1link)#.isdag
        #isU1dag = ifelse(typeof(U1link) <: Adjoint_GLink,true,false)

        #show(glinks)
        #println("in evaluate_gaugelinks!")
        #show(w)
        #println("numlinks = $numlinks")
        if numlinks == 1
            for i=1:NC
                for j=1:NC
                    Unew[j,i] = U[direction][j,i,indices...]
                end
            end

            #substitute_U!(Unew,U[direction])
            for μ=1:4
                shifted[μ] = indices[μ] + position[μ]# collect(position)
            end
            #shifted = collect(indices) .+ collect(position)
            for μ=1:4
                shifted[μ] += ifelse(shifted[μ] < 1,NN[μ],0)
                shifted[μ] += ifelse(shifted[μ] > NN[μ],-NN[μ],0)
            end
            #shiftedposition = Tuple(shifted)
            #shiftedposition = Tuple(collect(indices) .+ collect(position))
            for i=1:NC
                for j=1:NC
                    #Ushift1[j,i] = Unew[j,i,shiftedposition...]
                    #Ushift1[j,i] = Unew[j,i,shifted...]
                    Ushift1[j,i] = Unew[j,i,shifted[1],shifted[2],shifted[3],shifted[4]]
                end
            end

            #Ushift1 = shift_U(Unew,position)
            if isU1dag
                #println("Ushift1 ",Ushift1'[1,1,1,1,1,1])
                for i=1:NC
                    for j=1:NC
                        Unew[j,i] = Ushift1'[j,i]#,indices]
                    end
                end
                #substitute_U!(uout,Ushift1')
            else
                for i=1:NC
                    for j=1:NC
                        Unew[j,i] = Ushift1[j,i]#,indices]
                    end
                end
                #substitute_U!(uout,Ushift1)
            end
            return
        end

        #j = 1
        #U1link = glinks[1]
        #direction = get_direction(U1link)
        #position = get_position(U1link)
        #println("i = $i j = $j position = $position")
        #substitute_U!(Unew,U[direction])
        #for b=1:NC
        #    for a=1:NC
        #        Unew[a,b] = U[direction][a,b,indices...]
        #    end
        #end

        for μ=1:4
            shifted[μ] = indices[μ] + position[μ]# collect(position)
        end
        #shifted = collect(indices) .+ collect(position)
        for μ=1:4
            shifted[μ] += ifelse(shifted[μ] < 1,NN[μ],0)
            shifted[μ] += ifelse(shifted[μ] > NN[μ],-NN[μ],0)
        end
        #shiftedposition = Tuple(shifted)
        for i=1:NC
            for j=1:NC
                #Ushift1[j,i] = U[direction][j,i,shiftedposition...]
                #Ushift1[j,i] = U[direction][j,i,shifted...]
                Ushift1[j,i] = U[direction][j,i,shifted[1],shifted[2],shifted[3],shifted[4]]

            end
        end

        #Ushift1 = shift_U(Unew,position)

        #ix,iy,iz,it=(2,2,2,2)
        #println("posotion = $position")
        #pos = Tuple([ix,iy,iz,it] .+ collect(position))
        #U1 = Unew[:,:,pos...]
        #println("U1, ",Unew[:,:,pos...])
        #isU1dag = ifelse(typeof(U1link) <: Adjoint_GLink,true,false)

        for j=2:numlinks
            Ujlink = glinks[j]
            isUkdag = isdag(Ujlink)
            #isUkdag = ifelse(typeof(Ujlink) <: Adjoint_GLink,true,false)
            position = get_position(Ujlink)
            direction = get_direction(Ujlink)
            #println("j = $j position = $position")
            #println("a,b, $isUkdag , $isU1dag")
            for μ=1:4
                shifted[μ] = indices[μ] + position[μ]# collect(position)
            end
            #shifted = collect(indices) .+ collect(position)
            for μ=1:4
                shifted[μ] += ifelse(shifted[μ] < 1,NN[μ],0)
                shifted[μ] += ifelse(shifted[μ] > NN[μ],-NN[μ],0)
            end
            #shiftedposition2 = Tuple(shifted)
            ##  shiftedposition2 = Tuple(collect(indices) .+ collect(position))
            for b=1:NC
                for a=1:NC
                    #Ushift2[a,b] = U[direction][a,b,shiftedposition2...]
                    #Ushift2[a,b] = U[direction][a,b,shifted...]
                    Ushift2[a,b] = U[direction][a,b,shifted[1],shifted[2],shifted[3],shifted[4]]
                end
            end

            #Ushift2 = shift_U(U[direction],position)
            multiply_12!(uout_mat,Ushift1,Ushift2,j,isUkdag,isU1dag)

            #pos = Tuple([ix,iy,iz,it] .+ collect(position))
            #U2 = U[direction][:,:,pos...]
            #println("U1U2dag ", U1*U2')
            #substitute_U!(Unew,uout)
            for b=1:NC
                for a=1:NC
                    Unew[a,b] =uout_mat[a,b]
                end
            end

            #println("Unew ", Unew[:,:,ix,iy,iz,it])
            for b=1:NC
                for a=1:NC
                    Ushift1[a,b] = Unew[a,b]
                end
            end

            #Ushift1 = shift_U(Unew,origin)
            #println("uout ", uout[:,:,ix,iy,iz,it])
        end

        #println("uout2 ", uout[:,:,ix,iy,iz,it])

    end

    function evaluate_wilson_loops!(xout::T,w::Wilson_loop_set,U::Array{T,1},temps::Array{T,1}) where T<: AbstractGaugefields
        num = length(w)
        clear_U!(xout)
        Uold = temps[1]
        Unew = temps[2]

        for i=1:num
            wi = w[i]
            numloops = length(wi)

            shifts = calc_shift(wi)
            #println("shift ",shifts)

            loopk = wi[1]
            k = 1
            #println("k = $k shift: ",shifts[k])
            substitute_U!(Uold,U[loopk[1]])
            Ushift1 = shift_U(Uold,shifts[1])

            #gauge_shift_all!(temp1,shifts[1],U[loopk[1]])

            loopk1_2 = loopk[2]
            evaluate_wilson_loops_inside!(U,shifts,wi,Ushift1,Uold,Unew,numloops,loopk,loopk1_2)

            #=
            for k=2:numloops

                loopk = wi[k]
                #println("k = $k shift: ",shifts[k])
                #println("gauge_shift!(temp2,$(shifts[k]),$(loopk[1]) )")
                #clear!(temp2)
                Ushift2 = shift_U(U[loopk[1]],shifts[k])
                #gauge_shift_all!(temp2,shifts[k],U[loopk[1]])

                #multiply_12!(temp3,temp1,temp2,k,loopk,loopk1_2)
                multiply_12!(Unew,Ushift1,Ushift2,k,loopk,loopk1_2)
                Unew,Uold = Uold,Unew
                #Ushift1 = shift_U(Uold,(0,0,0,0))
                Ushift1 = Uold
                #temp1,temp3 = temp3,temp1
            end
            =#
            add_U!(xout,Uold)
            #println("i = $i ",Uold[1,1,1,1,1,1])
            #add_U!(xout,Ushift1)
            #add!(xout,temp1)

        end
    end

    function evaluate_wilson_loops_inside!(U,shifts,wi,Ushift1,Uold,Unew,numloops,loopk,loopk1_2)
        for k=2:numloops

            loopk = wi[k]
            #println("k = $k shift: ",shifts[k])
            #println("gauge_shift!(temp2,$(shifts[k]),$(loopk[1]) )")
            #clear!(temp2)
            Ushift2 = shift_U(U[loopk[1]],shifts[k])
            #gauge_shift_all!(temp2,shifts[k],U[loopk[1]])

            #multiply_12!(temp3,temp1,temp2,k,loopk,loopk1_2)
            multiply_12!(Unew,Ushift1,Ushift2,k,loopk,loopk1_2)

            Unew,Uold = Uold,Unew
            Ushift1 = shift_U(Uold,(0,0,0,0))

            #Ushift1 = Uold
            #temp1,temp3 = temp3,temp1
        end
    end

    function multiply_12!(temp3,temp1,temp2,k,loopk,loopk1_2)
        if loopk[2] == 1
            if k==2
                if loopk1_2 == 1
                    mul!(temp3,temp1,temp2)
                else
                    mul!(temp3,temp1',temp2)
                end
            else
                mul!(temp3,temp1,temp2)
            end
        elseif loopk[2] == -1
            if k==2
                if loopk1_2 == 1
                    mul!(temp3,temp1,temp2')
                else
                    mul!(temp3,temp1',temp2')
                end
            else
                mul!(temp3,temp1,temp2')
            end
        else
            error("Second element should be 1 or -1 but now $(loopk)")
        end
        return
    end

    function calculate_Plaquette(U::Array{T,1}) where T <: AbstractGaugefields
        error("calculate_Plaquette is not implemented in type $(typeof(U)) ")
    end

    function calculate_Plaquette(U::Array{T,1},temps::Array{T1,1}) where {T <: AbstractGaugefields,T1 <: AbstractGaugefields}
        return calculate_Plaquette(U,temps[1],temps[2])
    end

    function calculate_Plaquette(U::Array{T,1},temp::AbstractGaugefields{NC,Dim},staple::AbstractGaugefields{NC,Dim}) where {NC,Dim,T <: AbstractGaugefields}
        plaq = 0
        V = staple
        for μ=1:Dim
            construct_staple!(V,U,μ,temp)
            mul!(temp,U[μ],V')
            plaq += tr(temp)

        end

        return real(plaq*0.5)
    end

    function construct_staple!(staple::AbstractGaugefields,U,μ) where T <: AbstractGaugefields
        error("construct_staple! is not implemented in type $(typeof(U)) ")
    end

    function add_force!(F::Array{T1,1},U::Array{T2,1},temps::Array{<: AbstractGaugefields{NC,Dim},1};
        plaqonly = false,staplefactors::Union{Array{<: Number,1},Nothing} = nothing,factor = 1) where {NC,Dim,T1 <: AbstractGaugefields,T2 <: AbstractGaugefields}
        error("add_force! is not implemented in type $(typeof(F)) ")
    end

    function add_force!(F::Array{T1,1},U::Array{T2,1},temps::Array{<: AbstractGaugefields{NC,Dim},1};
                plaqonly = false,staplefactors::Union{Array{<: Number,1},Nothing} = nothing,factor = 1) where {NC,Dim,T1 <: TA_Gaugefields,T2 <: AbstractGaugefields}
        @assert length(temps) >= 3 "length(temps) should be >= 3. But $(length(temps))"
        #println("add force, plaqonly = $plaqonly")

        V = temps[3]
        temp1 = temps[1]
        temp2 = temps[2]

        for μ=1:Dim
            if plaqonly

                construct_double_staple!(V,U,μ,temps[1:2])

                mul!(temp1,U[μ],V') #U U*V
            else
                clear_U!(V)
                for i=1:gparam.numactions
                    loops = gparam.staples[i][μ]
                    evaluate_wilson_loops!(temp3,loops,U,[temp1,temp2])
                    add_U!(V,staplefactors[i],temp3)
                    #add_U!(V,gparam.βs[i]/gparam.β,temp3)
                end
                mul!(temp1,U[μ],V) #U U*V
            end

            Traceless_antihermitian_add!(F[μ],factor,temp1)
            #add_U!(F[μ],factor,temp2)
        end

    end

    #=
    function add_force!(F::Array{T,1},U::Array{T,1},temps::Array{<: AbstractGaugefields{NC,Dim},1},factor = 1) where {NC,Dim,T <: AbstractGaugefields,GP}
        @assert length(temps) >= 3 "length(temps) should be >= 3. But $(length(temps))"
        clear_U!(F)
        V = temps[3]
        temp1 = temps[1]
        temp2 = temps[2]

        for μ=1:Dim
            construct_double_staple!(V,U,μ,temps[1:2])

            mul!(temp1,U[μ],V') #U U*V

            a = temp1[:,:,1,1,1,1]
            println(a'*a)

            Traceless_antihermitian!(temp2,temp1)
            #println(temp2[1,1,1,1,1,1])
            a = temp2[:,:,1,1,1,1]
            println(a'*a)
            error("a")
            add_U!(F[μ],factor,temp2)
        end
    end
    =#

    function exptU!(uout::T,t::N,f::T1,temps::Array{T,1}) where {N <: Number, T <: AbstractGaugefields,T1 <: AbstractGaugefields} #uout = exp(t*u)
        error("expUt! is not implemented in type $(typeof(f)) uout: $(typeof(uout))")
    end

    function exptU!(uout::T,f::T1,temps::Array{T,1}) where {T <: AbstractGaugefields,T1 <: AbstractGaugefields}
        expU!(uout,1,f,temps)
    end

    function exp_aF_U!(W::Array{<: AbstractGaugefields{NC,Dim},1},a::N,F::Array{T1,1},U::Array{T,1},temps::Array{T,1}) where {NC,Dim,N <: Number, T <: AbstractGaugefields,T1 <: AbstractGaugefields} #exp(a*F)*U
        @assert a != 0 "Δτ should not be zero in expF_U! function!"
        expU = temps[1]
        temp1 = temps[2]
        temp2 = temps[3]
        #clear_U!(temp1)
        #clear_U!(temp2)
        #clear_U!(expU)
        #clear_U!(W)

        for μ=1:Dim
            exptU!(expU,a,F[μ],[temp1,temp2])
            mul!(W[μ],expU,U[μ])
        end

        set_wing_U!(W)
    end

    function staple_prime()
        loops_staple_prime = Array{Wilson_loop_set,2}(undef,4,4)
        for Dim=1:4
            for μ=1:Dim
                loops_staple_prime[Dim,μ] = make_plaq_staple_prime(μ,Dim)
            end
        end
        return loops_staple_prime
    end
    const loops_staple_prime_old = staple_prime()

    function construct_double_staple!(staple::AbstractGaugefields{NC,Dim},U::Array{T,1},μ,temps::Array{<: AbstractGaugefields{NC,Dim},1}) where {NC,Dim,T <: AbstractGaugefields}
        #println("mu = ",μ)
        #loops = loops_staple_prime_old[Dim,μ] #make_plaq_staple_prime(μ,Dim)
        #println("staple")
        #@time evaluate_wilson_loops!(staple,loops,U,temps)
        #println(staple[1,1,1,1,1,1])
        loops = loops_staple_prime[(Dim,μ)]
        evaluate_gaugelinks!(staple,loops,U,temps)
        #println(staple[1,1,1,1,1,1])
        #error("construct!!")
    end

    function construct_staple!(staple::AbstractGaugefields{NC,Dim},U::Array{T,1},μ,temp::AbstractGaugefields{NC,Dim}) where {NC,Dim,T <: AbstractGaugefields}
        U1U2 = temp
        firstterm = true

        for ν=1:Dim
            if ν == μ
                continue
            end

            #=
                    x+nu temp2
                    .---------.
                    I         I
            temp1 I         I
                    I         I
                    .         .
                    x        x+mu
            =#
            U1 = U[ν]
            U2 = shift_U(U[μ],ν)
            #println(typeof(U1))
            mul!(U1U2,U1,U2)
            #error("test")

            U3 = shift_U(U[ν],μ)
            #mul!(staple,temp,Uμ')
            #  mul!(C, A, B, α, β) -> C, A B α + C β
            if firstterm
                β = 0
                firstterm = false
            else
                β = 1
            end
            mul!(staple,U1U2,U3',1,β) #C = alpha*A*B + beta*C

            #println("staple ",staple[1,1,1,1,1,1])

            #mul!(staple,U0,Uν,Uμ')
        end
        set_wing_U!(staple)
    end

    function Base.size(U::T) where T <: AbstractGaugefields
        error("Base.size is not implemented in type $(typeof(U)) ")
    end

    function calculate_Polyakov_loop(U::Array{T,1},temp1::AbstractGaugefields{NC,Dim},temp2::AbstractGaugefields{NC,Dim}) where {NC,Dim,T <: AbstractGaugefields}
        Uold = temp1
        Unew = temp2
        shift = zeros(Int64,Dim)

        μ = Dim
        _,_,NN... = size(U[1]) #NC,NC,NX,NY,NZ,NT 4D case
        lastaxis = NN[end]
        #println(lastaxis)

        substitute_U!(Uold,U[μ])
        for i=2:lastaxis
            shift[μ] = i-1
            U1 = shift_U(U[μ],Tuple(shift))
            mul_skiplastindex!(Unew,Uold,U1)
            Uold,Unew = Unew,Uold
        end

        set_wing_U!(Uold)
        poly = partial_tr(Uold,μ)/prod(NN[1:Dim-1])
        return poly

    end

    function mul_skiplastindex!(c::T,a::T1,b::T2) where {T <: AbstractGaugefields, T1 <: Abstractfields,T2 <: Abstractfields}
        error("mul_skiplastindex! is not implemented in type $(typeof(c)) ")
    end

    function Traceless_antihermitian(vin::T) where T <: AbstractGaugefields
        #error("Traceless_antihermitian is not implemented in type $(typeof(vin)) ")
        vout = deepcopy(vin)
        Traceless_antihermitian!(vout,vin)
        return vout
    end

    function Traceless_antihermitian_add!(U::T,factor,temp1) where T <: AbstractGaugefields
        error("Traceless_antihermitian_add! is not implemented in type $(typeof(U)) ")
    end

    function Traceless_antihermitian!(vout::T,vin::T) where T <: AbstractGaugefields
        error("Traceless_antihermitian! is not implemented in type $(typeof(vout)) ")
    end

    function add_U!(c::T,a::T1) where {T<: AbstractGaugefields,T1 <: Abstractfields}
        error("add_U! is not implemented in type $(typeof(c)) ")
    end

    function add_U!(c::Array{<: AbstractGaugefields{NC,Dim},1},α::N,a::Array{T1,1}) where {NC,Dim,T1 <: Abstractfields, N<:Number}
        for μ=1:Dim
            add_U!(c[μ],α,a[μ])
        end
    end

    function add_U!(c::T,a::T1,iseven) where {T<: AbstractGaugefields,T1 <: Abstractfields}
        error("add_U! is not implemented in type $(typeof(c)) ")
    end

    function add_U!(c::Array{<: AbstractGaugefields{NC,Dim},1},α::N,a::Array{T1,1},iseven) where {NC,Dim,T1 <: Abstractfields, N<:Number}
        for μ=1:Dim
            add_U!(c[μ],α,a[μ],iseven)
        end
    end

    function add_U!(c::T,α::N,a::T1) where {T<: AbstractGaugefields,T1 <: Abstractfields, N<:Number}
        error("add_U! is not implemented in type $(typeof(c)) ")
    end

    function LinearAlgebra.mul!(c::T,a::T1,b::T2,α::Ta,β::Tb) where {T<: AbstractGaugefields,T1 <: Abstractfields,T2 <: Abstractfields,Ta <: Number, Tb <: Number}
        error("LinearAlgebra.mul! is not implemented in type $(typeof(c)) ")
    end

    function LinearAlgebra.mul!(c::T,a::N,b::T2) where {T<: AbstractGaugefields,N <: Number ,T2 <: Abstractfields}
        error("LinearAlgebra.mul! is not implemented in type $(typeof(c)) ")
    end

    function partial_tr(a::T,μ) where T<: Abstractfields
        error("partial_tr is not implemented in type $(typeof(a)) ")
    end

    function LinearAlgebra.tr(a::T) where T<: Abstractfields
        error("LinearAlgebra.tr! is not implemented in type $(typeof(a)) ")
    end

    """
        Tr(A*B)
    """
    function LinearAlgebra.tr(a::T,b::T) where T<: Abstractfields
        error("LinearAlgebra.tr! is not implemented in type $(typeof(a)) ")
    end

    """
    M = (U*δ_prev) star (dexp(Q)/dQ)
    Λ = TA(M)
    """
    function construct_Λmatrix_forSTOUT!(Λ,δ_prev::T,Q,u::T) where T <: AbstractGaugefields
        error("construct_Λmatrix_forSTOUT! is not implemented in type $(typeof(u)) ")
    end

    const eps_Q = 1e-18

    function calc_Λmatrix!(Λ,M,NC)
        #println("M= ", M)
        if NC == 1
            #Λ = -(M - M')
            @. Λ[:,:] = (M - M')
        elseif NC == 2
            #Λ = (1/2)*(M - M') - (1/(2NC))*tr(M - M')*I0_2
            @. Λ[:,:] = (1/2)*(M - M')
            trM = (1/(2NC))*(M[1,1]-conj(M[1,1]) + M[2,2] - conj(M[2,2]))#  tr(M - M')
            #trM = (1/(2NC))*tr(M - M')
            for i=1:NC
                Λ[i,i] += - trM
            end
            #Λ = 2*Λ
        else
            @. Λ[:,:] = (1/2)*(M - M')
            trM = (1/(2NC))*tr(M - M')
            for i=1:NC
                Λ[i,i] += - trM
            end
        end
        #display(Λ)
        #println("\t")
        #exit()
        return
    #        return Λ
    end

    function calc_Mmatrix!(Mn,δn_prev,Qn,Un,u::AbstractGaugefields{2,Dim},tempmatrices) where {Dim}
        Unδn = tempmatrices[1]
        B = tempmatrices[2]
        tmp_matrix1 = tempmatrices[3]
        NC = 2

        trQ2 = 0.0
        for i=1:2
            for j=1:2
                trQ2 += Qn[i,j]*Qn[j,i]
            end
        end

        if abs(trQ2) > eps_Q
            q = sqrt((-1/2)*trQ2)
            calc_Bmatrix!(B,q,Qn,NC)
            for i=1:2
                for j=1:2
                    tmp_matrix1[j,i] = Un[j,i]
                end
            end

            mul!(Unδn,tmp_matrix1,δn_prev)
            trsum = 0.0im
            for i=1:2
                for j=1:2
                    trsum += Unδn[i,j]*B[j,i]
                end
            end
            for i=1:2
                for j=1:2
                    Mn[j,i] = (sin(q)/q)*Unδn[j,i] + trsum*Qn[j,i]
                end
            end
        end
    end

    function calc_Mmatrix!(Mn,δn_prev,Qn,Un,u::AbstractGaugefields{3,Dim},tempmatrices) where {Dim}
        Unδn = tempmatrices[1]
        tmp_matrix1 = tempmatrices[2]
        tmp_matrix2 = tempmatrices[3]
        #println("Qn ", Qn)
        trQ2 = 0.0
        for i=1:3
            for j=1:3
                trQ2 += Qn[i,j]*Qn[j,i]
            end
        end
        #println("tr", trQ2)

        if abs(trQ2) > eps_Q
            Qn ./= im
            #println("Qn b ",Qn)
            f0,f1,f2,b10,b11,b12,b20,b21,b22 = calc_coefficients_Q(Qn)
            for i=1:3
                for j=1:3
                    tmp_matrix1[j,i] = Un[j,i]
                end
            end
            mul!(Unδn,tmp_matrix1,δn_prev)

            B1 = tmp_matrix1
            B1 .= 0
            B2 = tmp_matrix2
            B2 .= 0
            for i=1:3
                B1[i,i] = b10
                B2[i,i] = b20
            end
            for j=1:3
                for i=1:3
                    B1[i,j] += b11*Qn[i,j]
                    B2[i,j] += b21*Qn[i,j]
                    for k=1:3
                        B1[i,j] += b12*Qn[i,k]*Qn[k,j]
                        B2[i,j] += b22*Qn[i,k]*Qn[k,j]
                    end
                end
            end
            #println("coeff, ",(f0,f1,f2,b10,b11,b12,b20,b21,b22))
            #println("B1 ",B1)
            #println("B2 ",B2)

            trB1 = 0.0
            trB2 = 0.0
            for i=1:3
                for j=1:3
                    trB1 += Unδn[i,j]*B1[j,i]
                    trB2 += Unδn[i,j]*B2[j,i]
                end
            end

            for j=1:3
                for i=1:3
                    Mn[i,j] = trB1*Qn[i,j] + f1*Unδn[i,j]
                    for k=1:3
                        Mn[i,j] += trB2*Qn[i,k]*Qn[k,j]+f2*(Qn[i,k]*Unδn[k,j]+Unδn[i,k]*Qn[k,j])
                    end
                end
            end
            Mn ./= im
        end
    end

    function calc_Mmatrix!(Mn,δn_prev,Qn,Un,u::AbstractGaugefields{NC,Dim},tempmatrices) where {NC,Dim}
        error("not supported yet")

        @assert NC > 3 "NC > 3 not NC = $NC"
        Unδn = tempmatrices[1]
        B = tempmatrices[2]
        tempmatrix = tempmatrices[3]

        trQ2 = 0.0
        for i=1:NC
            for j=1:NC
                trQ2 += Qn[i,j]*Qn[j,i]
            end
        end

        if abs(trQ2) > eps_Q
            e,v = eigen(Qn)
            mul!(Unδn,Un,δn_prev)
            #=
                    A star dexp(Q)/dQ = \sum_{n=0}^{infty} \frac{1}{n!}
                                            \sum_{k=0}^{n-1} i^{n-1-k}P^+ D^{n-1-k} P A P^+ D^k P i^k
                                    = P^+ (\sum_{n=0}^{infty} (i^{n-1}/n!) sum_{k=0}^{n-1} D^{n-1-k} B D^k)  P
                                    B = P A P+
                =#

            mul!(tempmatrix,Unδn,v)
            mul!(B,v',tempmatrix)

        end

    end

    function calc_Bmatrix!(B,q,Q,NC)
        @assert NC == 2 "NC should be 2! now $NC"
        mul!(B,cos(q)/q -sin(q)/q^2,Q)
        for i=1:NC
            B[i,i] += -sin(q)
        end
        B .*= -1/2q
        #B[:,:] .= (cos(q)/q -sin(q)/q^2 )*Q

        #q = sqrt((-1/2)*tr(Q^2))
        #B = -(-sin(q)*I0_2 +(cos(q)/q -sin(q)/q^2 )*Q)*(1/2q)
    end

    function calc_coefficients_Q(Q)
        @assert size(Q) == (3,3)
        c0 = Q[1,1]*Q[2,2]*Q[3,3]+Q[1,2]*Q[2,3]*Q[3,1]+Q[1,3]*Q[2,1]*Q[3,2]-Q[1,3]*Q[2,2]*Q[3,1]-Q[1,2]*Q[2,1]*Q[3,3]-Q[1,1]*Q[2,3]*Q[3,2]
        #@time cdet = det(Q)
        ##println(c0,"\t",cdet)
        #exit()

        c1 = 0.0
        for i=1:3
            for j=1:3
                c1 += Q[i,j]*Q[j,i]
            end
        end
        c1 /= 2
        c0max = 2*(c1/3)^(3/2)
        θ = acos(c0/c0max)
        u = sqrt(c1/3)*cos(θ/3)
        w = sqrt(c1)*sin(θ/3)
        ξ0 = sin(w)/w
        ξ1 = cos(w)/w^2 - sin(w)/w^3

        emiu = exp(-im*u)
        e2iu = exp(2*im*u)

        h0 = (u^2-w^2)*e2iu + emiu*(
            8u^2*cos(w)+2*im*u*(3u^2+w^2)* ξ0
        )
        h1 = 2u*e2iu-emiu*(
            2u*cos(w)-im*(3u^2-w^2)* ξ0
        )
        h2 = e2iu - emiu*(cos(w)+3*im*u*ξ0)

        denom = 9u^2-w^2

        f0 = h0/denom
        f1 = h1/denom
        f2 = h2/denom

        r10 = 2*(u+im*(u^2-w^2))*e2iu +
                2*emiu*(
                    4u*(2-im*u)*cos(w) +
                    im*(9u^2+w^2-im*u*(3u^2+w^2))*ξ0
                )
        r11 = 2*(1+2*im*u)*e2iu+
                emiu*(
                    -2*(1-im*u)*cos(w)+
                    im*(6u+im*(w^2-3u^2))*ξ0
                )
        r12 = 2*im*e2iu + im*emiu*(
            cos(w) -3*(1-im*u)*ξ0
        )
        r20 = -2*e2iu+2*im*u*emiu*(
            cos(w)+(1+4*im*u)*ξ0+3u^2*ξ1
        )
        r21 = -im*emiu*(
            cos(w)+(1+2*im*u)*ξ0 -
            3*u^2*ξ1
        )
        r22 = emiu*(
            ξ0-3*im*u*ξ1
        )
        b10 = (
            2*u*r10+(3u^2-w^2)*r20 - 2*(15u^2+w^2)*f0
            )/(
                2*denom^2
            )

        b11 = (
                2*u*r11+(3u^2-w^2)*r21 - 2*(15u^2+w^2)*f1
                )/(
                    2*denom^2
                )
        b12 = (
            2*u*r12+(3u^2-w^2)*r22 - 2*(15u^2+w^2)*f2
            )/(
                2*denom^2
            )
        b20 = (
            r10 - 3*u*r20 - 24*u*f0
            )/(2*denom^2)
        b21 = (
                r11 - 3*u*r21 - 24*u*f1
                )/(2*denom^2)
        b22 = (
            r12 - 3*u*r22 - 24*u*f2
            )/(2*denom^2)

        return f0,f1,f2,b10,b11,b12,b20,b21,b22
    end

    function staggered_phase(μ,iii...)
        error("staggered_phase is not implemented")
    end

    function normalize!(u::Array{ComplexF64,2})
        NC,_ = size(u)
        if NC == 2
            normalize2!(u)
        elseif NC == 3
            normalize3!(u)
        else
            normalizeN!(u)
        end
        return
    end

    function normalize2!(u)
        α = u[1,1]
        β = u[2,1]
        detU = abs(α)^2 + abs(β)^2
        u[1,1] = α/detU
        u[2,1] = β/detU
        u[1,2] = -conj(β)/detU
        u[2,2] = conj(α)/detU
    end

    function normalizeN!(u)
        gramschmidt!(u)
    end

    function normalize3!(u)
        w1 = 0
        w2 = 0
        for ic=1:3
            w1 += u[2,ic]*conj(u[1,ic])
            w2 += u[1,ic]*conj(u[1,ic])
        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], "\t",u[1,2],"\t", u[1,3])
        end

        w1 = -w1/w2

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

        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] = x4
        u[2,2] = x5
        u[2,3] = x6

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

        u[1,1] = u[1,1]*w2
        u[1,2] = u[1,2]*w2
        u[1,3] = u[1,3]*w2
        u[2,1] = u[2,1]*w3
        u[2,2] = u[2,2]*w3
        u[2,3] = u[2,3]*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

        m3complv3!(u)
    end

    function m3complv3!(a)
        aa = zeros(Float64,18)
        aa[ 1] = real( a[1,1])
        aa[ 2] = imag(a[1,1])
        aa[ 3] = real( a[1,2])
        aa[ 4] = imag(a[1,2])
        aa[ 5] = real( a[1,3])
        aa[ 6] = imag(a[1,3])
        aa[ 7] = real( a[2,1])
        aa[ 8] = imag(a[2,1])
        aa[ 9] = real( a[2,2])
        aa[10] = imag(a[2,2])
        aa[11] = real( a[2,3])
        aa[12] = imag(a[2,3])

        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] = aa[13]+im*aa[14]
        a[3,2] = aa[15]+im*aa[16]
        a[3,3] = aa[17]+im*aa[18]
        return
    end

    function gramschmidt!(v)
        n = size(v)[1]
        for i=1:n
            for j=1:i-1
                v[:,i] = v[:,i] - v[:,j]'*v[:,i]*v[:,j]
            end
            v[:,i] = v[:,i]/norm(v[:,i])
        end
    end

    function gramschmidt_special!(v)
        n = size(v)[1]
        #vdet = det(v)

        vnorm1 = norm(v[:,1])
        for i=1:n
            #vnorm[i] = norm(v[:,i])
            for j=1:i-1
                v[:,i] = v[:,i] - v[:,j]'*v[:,i]*v[:,j]
            end
            v[:,i] = v[:,i]/norm(v[:,i])
        end
        for i=1:n
            #v[:,i] = v[:,i]*vnorm[i]
            v[:,i] = v[:,i]*vnorm1
        end
    end

end


text:
@testset "CHSF Descriptor" begin

@info("Testing CHSF Descriptor for DC 4x4x4 Si with cutoff=6.5, n=2, l=2")
using CHSF, JuLIP, Test

at = bulk(:Si, cubic=true)
desc = chsf(at, 6.5, n=2, l=2)
#round.(desc[1], digit=7) == chsfpy
chsf_ref = [10.3698237,1.4503467,-8.2118063,51.6882200,-53.0113716,69.8233316] #n=2,l=2 case
chsf_now = vcat(desc[1,:]...)
println("CHSF.jl returns:",chsf_now)
println("Reference:",chsf_ref)
println(@test chsf_now  ≈  chsf_ref)

end

ext: jl


text:
include("$(@__DIR__)/../src/FiniteGroups.jl")
using .FiniteGroups, LinearAlgebra, Test
import .FiniteGroups: findprime, primitive_root, vec_split, normalize_modp!, maptocomplex, eigvecs_modp
import .FiniteGroups: hermiteform!, column_reduce!, kill_column!, mod!, set_prime_column!
g = pointgroup(32)
h = class_multab(g)
n, p = findprime(g)
Z = primitive_root(n, p)
ξ = exp(1im * 2π/n)
NC = size(h, 1)
v = [I(NC)]

for i = 1:NC
    hi = h[i, :, :]
    v = vcat([vec_split(hi, vi, p) for vi in v]...)
    all(isone(size(vi, 2)) for vi in v) && break
end
"""
mat = vcat(transpose.(v)...)
for i = 1:size(mat, 1)
    normalize_modp!(g, view(mat, i, :), p)
end
maptocomplex.(mat, Z, n, ξ)
"""


text:
using FinancialToolbox, Test, FinancialFFT, FinancialMonteCarlo, ReverseDiff

A = 600.0;
N = 18;

S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01
sigma = 0.2;
lam = 15.0;
mu1 = 0.03;
sigma1 = 0.02;
zero_rate = ZeroRate(r);

Model = BlackScholesProcess(sigma, Underlying(S0, d));

EUData = EuropeanOption(T, K)
Nsim = 100000;
Nstep = 30;
mc = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());
analytic_result = blsprice(S0, K, r, T, sigma, d);
function test_reverse_diff(x, y, toll)
    @test abs(sum(x .- y)) < toll
    nothing
end
x = Float64[sigma, S0, r]
method_cm = CarrMadanMethod(A, N);
method_lewis = LewisMethod(700.0, 200000);
method_cm_lewis = CarrMadanLewisMethod(A, N);
f(x) = blsprice(x[2], K, x[3], T, x[1], d);
f_lewis(x) = pricer(BlackScholesProcess(x[1], Underlying(x[2], d)), ZeroRate(x[3]), method_lewis, EUData);
analytic_result = ReverseDiff.gradient(f, x);
@show result_lewis = ReverseDiff.gradient(f_lewis, x);
toll = 1e-2;
@testset "ReverseDiff test" begin
    test_reverse_diff(analytic_result, result_lewis, toll)
end


text:
# Unit tests for unit_convention.jl

@testset "unit_convention.jl" begin
    @test Atomistic.MASS_UNIT == aunit(u"kg")
    @test Atomistic.LENGTH_UNIT == aunit(u"m")
    @test Atomistic.ENERGY_UNIT == aunit(u"J")
    @test Atomistic.TIME_UNIT == aunit(u"s")
    @test Atomistic.VELOCITY_UNIT == aunit(u"m/s")
    @test Atomistic.TEMPERATURE_UNIT == aunit(u"K")
end


text:
using Test
using DataFrames
using CSV
using ChemistryFeaturization.FeatureDescriptor
using ChemistryFeaturization.Featurization

@testset "GraphNodeFeaturization" begin

    # custom lookup table...
    test_df = CSV.read(abspath(@__DIR__, "..", "test_data", "lookup_table.csv"), DataFrame)

    local fnames = ["X", "Block", "Atomic mass"]

    @testset "Encode" begin
        # make sure both constructors give the same results
        features = ElementFeatureDescriptor.(fnames)
        fzn1 = GraphNodeFeaturization(features)
        fzn2 = GraphNodeFeaturization(fnames)

        @test output_shape(fzn1) == 24 # 10 (default length for X) + 4 (default for block) + 10 (default for atomic mass)

        triangle_C_1 = AtomGraph(Float32.([0 1 1; 1 0 1; 1 1 0]), ["C", "C", "C"])
        triangle_C_2 = AtomGraph(Float32.([0 1 1; 1 0 1; 1 1 0]), ["C", "C", "C"])
        # @test fzn1 == fzn2

        featurized_1, featurized_2 = encode.([fzn1], [triangle_C_1, triangle_C_2]) # encode can be broadcasted
        @test featurized_1 == featurized_2

        featurized_1, featurized_2 = encode.([fzn1, fzn2], [triangle_C_1, triangle_C_2]) # encode can be broadcasted
        @test featurized_1 == featurized_2
    end

    @testset "Decode" begin
        # test that other options work properly

        # TODO - Figure out a way to robustly test this. Would creating a custom FD type whose encoded/decoded value you know beforehand and directly compare with be the best thing to do here?
        # fzn3 = GraphNodeFeaturization(fnames, nbins = 2)
        # F2 = AtomGraph(Float32.([0 1; 1 0]), ["F", "F"])

        # encode(F2, fzn3)
        # decoded_matrix = decode(fzn3, F2.encoded_features)
        # decoded_ag = decode(F2)
        # enc1 = F2.encoded_features
        # @test all(
        #     map(d -> d[1]["Block"] == d[2]["Block"] == "p", [decoded_matrix, decoded_ag]),
        # )
        # fzn4 = GraphNodeFeaturization(fnames, nbins = [2, 4, 2])
        # F2 = AtomGraph(Float32.([0 1; 1 0]), ["F", "F"])
        # encode(F2, fzn4)
        # @test all(F2.encoded_features .== enc1)
    end

    # encodable_elements
    @testset "Encodable Elements" begin
        fzn3 = GraphNodeFeaturization(ElementFeatureDescriptor.(["Boiling point", "6d"]))
        fzn4 = GraphNodeFeaturization(ElementFeatureDescriptor.(["7s", "Valence"]))
        @test encodable_elements(fzn3) == ["Ac", "Th", "U"]
        @test encodable_elements(fzn4) == ["Fr", "Ra", "Ac", "Th"]

        # Custom lookup_table
        feature_1 = ElementFeatureDescriptor("MeaningOfLife", test_df)   # zero-value case - `As` has a value = 0
        feature_2 = ElementFeatureDescriptor("neg_nums", test_df)
        feature_3 = ElementFeatureDescriptor("first_letter", test_df)
        feature_4 = ElementFeatureDescriptor("noarsenic", test_df) # missing value case - `As` has a missing value

        @test encodable_elements(GraphNodeFeaturization([feature_1, feature_2])) ==
              ["C", "As", "Tc"]
        @test encodable_elements(GraphNodeFeaturization([feature_1, feature_4])) ==
              ["C", "Tc"]
        @test encodable_elements(GraphNodeFeaturization([feature_2, feature_3])) ==
              ["C", "As", "Tc"]
        @test encodable_elements(GraphNodeFeaturization([feature_2, feature_4])) ==
              ["C", "Tc"]
    end

    # chunk_vec helper fcn
    @testset "chunk_vec" begin
        vec = [1, 1, 0, 1, 0, 1, 0]
        @test_throws AssertionError chunk_vec(vec, [3, 3])
        @test chunk_vec(vec, [4, 1, 2]) == [[1, 1, 0, 1], [0], [1, 0]]
    end

end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__,"../.."))
include(joinpath(@__DIR__,"model.jl"))
include(joinpath(@__DIR__,"plotting.jl"))
include(joinpath(@__DIR__,"utilities.jl"))
using ITensors
using ITensors.HDF5

# Load the results in
## NOTE: NEED TO CHANGE UTILS BACK TO TAKE IN PHONON ARGUMENTS TO PARAMS ##
#loadpath = "/Users/nicole/sherlock/ITensors.jl/projects/1dhubbardholstein/outputs/07:28:01_80N_1t_8U_0ω_0g0_0g1_0.1doping_0phonons.h5"
loadpath = "/Users/nicole/Dropbox/Grad school/Devereaux lab/ITensors.jl/projects/1dhubbardholstein/outputs/2-5-2_16:36:24_8N_1t_8U_0ω_0g0_0g1_0λ_0doping_1phonons.h5"
params, hubbholst, dmrg_results, equilibrium_corr, tebd_results = load_structs(loadpath)

# Visualize things we want
plot_equilibrium_correlations(dmrg_results,"spin",hubbholst)
plot_equilibrium_correlations(dmrg_results,"sSC",hubbholst)
plot_entropy(tebd_results)
plot_correlation_function(tebd_results)
plot_spectral_function(tebd_results,params)


text:
using SymbolServer
using Test

@testset "SymbolServer" begin
    s = SymbolServerProcess()

    @test isa(s, SymbolServerProcess)
end


text:
"""
`optimalPLD_solveSAO_LP_Benders( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, x_initial, A_feas, b_feas)`

Solve the LP relaxation of the Bertsimas and Mišić first-choice PLD formulation using Benders decomposition.

`n` = number of candidate products.

`revenues` = array of n numbers representing marginal revenue/profit of each of the n products. (These are the pi_i's in the Bertsimas and Mišić paper.)

`numPermutations` = number of rankings / customer types.

`lambda` = array specifying the probability of each ranking / customer type; should be nonnegative and sum to 1.

`orderings` = 2D array of integers encoding each ranking. orderings[k,j] contains the index of ranking k's jth top product; (n+1) encodes the no-purchase option. For example, if n = 5 and

	orderings[k,:] = [3,5,2,6,1,4]

then the customer prefers 3 < 5 < 2 < 0 (no-purchase option).

`MIP_timeLimit` = maximum time limit on each solve of the Benders LP master problem. If set to 0, no time limit is imposed. (This argument is not used in any experiment scripts, since master LP generally solves quickly.)

`x_initial` = array with with n values (between 0 and 1) indicating the initial values of the x_i variables to use to warm-start the Benders decomposition method. If empty, no warm-start solution will be used.

`A_feas`, `b_feas` = 2D array and 1D array specifying the system of linear inequalities Ax <= b for the PLD problem. This can be used to impose side constraints (this is the Cx <= d constraint system in the formulation in Bertsimas and Mišić). If b_feas is empty, then no constraints will be added (i.e., it will solve the unconstrained problem).

Returns:

`x_val` = array of optimal x_i variables for the LP relaxation.

`expRevenue` = expected per-customer revenue/profit of the LP solution.

`solutiontime` = time required to solve the problem.

`rankingid_list`, `alpha_list`, `beta_list`, `gamma_list` = 1D, 2D, 2D and 1D arrays containing information on cuts generated in Benders LP solution process. Each row of these arrays corresponds to a cut for a particular customer type. rankingid_list[i] is the customer type of the ith cut; alpha_list[i,:], beta_list[i,:] and gamma_list[i] are the dual variable values defining the ith cut (see Bertsimas and Mišić for more details.)

"""
function optimalPLD_solveSAO_LP_Benders( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, x_initial, A_feas, b_feas)

	tic()
	maxRev = maximum(revenues);

	if (MIP_timeLimit > 0)
		m = Model( solver = GurobiSolver(TimeLimit = MIP_timeLimit, Method = 2, Crossover = 0))
	else
		m = Model( solver = GurobiSolver(Method = 2, Crossover = 0))
	end

	@variable(m, 0 <= x[1:n] <= 1)
	@variable(m, 0 <= subobj[1:numPermutations] <= maxRev)

	for t=1:length(b_feas)
		@constraint(m, sum( A_feas[t,i]*x[i] for i=1:n) <= b_feas[t])
	end

	@objective(m, Max, sum( lambda[k] * subobj[k] for k=1:numPermutations))

	solve(m)

	# BEGIN preprocessing for inline
	temp_rev = [revenues;0];
	revenues_ordered_mat = zeros(numPermutations, n+1);
	sorted_rev_inds_mat = zeros(Int64, numPermutations, n+1);
	p2r_mat = zeros(Int64, numPermutations, n+1);
	for k in 1:numPermutations
		single_ordering = convert(Array{Int64}, orderings[k,:][:]);
		revenues_ordered =  temp_rev[single_ordering[:]]
		revenues_ordered_mat[k,:] = revenues_ordered
		sorted_rev_inds_mat[k,:] = sortperm(revenues_ordered[:], rev= true);
		p2r_mat[k,:] = sortperm(single_ordering[:]); # single_ordering is r2p
	end

	# END preprocessing for inline test

	bestUB = maxRev;
	bestLB = 0;
	x_val = getvalue(x)[:];
	if (!isempty(x_initial))
		oracle_subobj_vec = zeros(numPermutations)
		for k=1:numPermutations
			temp_ordering = convert(Array{Int64}, orderings[k,:][:]);
			temp_subobj_val, alpha_val, beta_val, gamma_val = oracle_inline_fn( n, x_initial, temp_ordering,
																				revenues_ordered_mat[k,:][:], temp_rev,
																				sorted_rev_inds_mat[k,:][:], p2r_mat[k,:][:] )

			oracle_subobj_vec[k] = temp_subobj_val;
		end
		x_val = x_initial
		bestLB = sum( oracle_subobj_vec .* lambda);
	end
	best_x_val = x_val;

	rankingid_list = zeros(Int64,0);
	alpha_list = zeros(Float64,0,n+1);
	beta_list = zeros(Float64,0,n+1);
	gamma_list = zeros(Float64,0);

	tol_UB_LB = 0.01;
	tol_subobj = 0.01;
	foundAcceptableSolution = false;
	itercounter = 0;

	while ( !foundAcceptableSolution )
		itercounter = itercounter + 1;
		bestUB = getobjectivevalue(m);
		solution_subobj_vec = getvalue(subobj)[:];
		oracle_subobj_vec = zeros(numPermutations)
		if (itercounter != 1)
			x_val = getvalue(x)[:];
		end

		loop_start = time();
		alpha_mat = zeros(Float64, numPermutations,n+1);
		beta_mat = zeros(Float64, numPermutations,n+1);
		gamma_mat = zeros(Float64, numPermutations);
		oracle_subobj_vec = zeros(Float64, numPermutations)

		for k=1:numPermutations
			temp_ordering = convert(Array{Int64}, orderings[k,:][:]);
			temp_subobj_val, alpha_val, beta_val, gamma_val = oracle_inline_fn( n, x_val, temp_ordering,
																				revenues_ordered_mat[k,:][:], temp_rev,
																				sorted_rev_inds_mat[k,:][:], p2r_mat[k,:][:] )
			alpha_mat[k,:] = alpha_val
			beta_mat[k,:] = beta_val
			gamma_mat[k] = gamma_val
			oracle_subobj_vec[k] = temp_subobj_val;
		end

		for k=1:numPermutations
			if (oracle_subobj_vec[k] < solution_subobj_vec[k] - tol_subobj)
				alpha_val = alpha_mat[k,:][:]
				beta_val = beta_mat[k,:][:]
				gamma_val = gamma_mat[k];
				@constraint(m,  subobj[k] <= gamma_val + sum(alpha_val[i] * x[i] for i=1:n) + sum(beta_val[i] * (1-x[i]) for i=1:n) )

				push!(rankingid_list, k)
				alpha_list = vcat(alpha_list, alpha_val' )
				beta_list = vcat(beta_list, beta_val')
				push!(gamma_list, gamma_val)
			end
		end
		loop_etime = time() - loop_start

		@show loop_etime
		@show oracle_subobj_vec
		@show solution_subobj_vec

		currentLB = sum( oracle_subobj_vec .* lambda);
		if (currentLB > bestLB)
			bestLB = currentLB
			best_x_val = x_val;
		end
		gap = round( (bestUB - bestLB) / bestUB * 100, 2);
		println("Iter $itercounter : bestLB = $bestLB, bestUB = $bestUB, gap = $gap")

		if (gap < tol_UB_LB)
			foundAcceptableSolution = true;
		else
			solve(m);
		end
	end

	solutiontime = toc()

	println(best_x_val)

	expRevenue = getobjectivevalue(m)

	return best_x_val, expRevenue, solutiontime, rankingid_list, alpha_list, beta_list, gamma_list, itercounter
end

function oracle_inline_fn( n::Int64, x_val::Array{Float64}, single_ordering::Array{Int64}, revenues_ordered::Array{Float64}, temp_rev::Array{Float64}, sorted_rev_inds::Array{Int64}, p2r::Array{Int64} )
#function oracle_inline_fn( n, x_val, single_ordering, revenues_ordered, temp_rev, sorted_rev_inds, p2r )
	temp_x = [x_val; 1.0];
	x_ordered = temp_x[single_ordering];

	remainingbudget = 1.0;
	ub = 1 - x_ordered;

	#minimizing_inds = zeros(Int32, size(y))
	minimizing_ind = 0;

	is_alpha_nonneg = zeros(Int64, n+1)
	is_beta_nonneg = zeros(Int64, n+1)
	is_beta_link = zeros(Int64, n+1);
	#is_gamma = zeros(Int64, n+1);
	gamma_ind = -1;
	# prev_remainingbudget = remainingbudget;

	# tic()
	# Primal phase
	guess_revenue = 0.0;
	y_rev_ind = 0.0;
	@inbounds for cur_prod=1:(n+1)
		# println("Loop body")
		rev_ind = sorted_rev_inds[cur_prod];

		if (rev_ind == 1)
			min_ub = Inf;
		else
			min_ub = minimum(ub[1:rev_ind-1]);
		end

		y_rev_ind = 0.0;
		if ( remainingbudget <= min_ub )
			if (remainingbudget <= x_ordered[rev_ind])
				minimizing_ind = 1
				y_rev_ind = remainingbudget
			else
				minimizing_ind = 2
				y_rev_ind = x_ordered[rev_ind]
			end
		else
			if (x_ordered[rev_ind] <= min_ub) ###
				minimizing_ind = 2
				y_rev_ind = x_ordered[rev_ind]
			else
				minimizing_ind = 2 + indmin(ub[1:rev_ind-1]);
				y_rev_ind = min_ub;
			end
		end

		guess_revenue += y_rev_ind * revenues_ordered[rev_ind];

		# update
		# prev_remainingbudget = remainingbudget
		remainingbudget = remainingbudget - y_rev_ind; #y[rev_ind];
		#ub[1:rev_ind-1] = ub[1:rev_ind-1] - y[rev_ind];
		ub[1:rev_ind-1] -=  y_rev_ind; #y[rev_ind];

		# println("Before dual conditionals")
		if (minimizing_ind == 2)
			is_alpha_nonneg[rev_ind] = 1
		elseif (minimizing_ind > 2)
			temp = minimizing_ind - 2; #minimizing_inds[rev_ind] - 2;

			if (is_beta_nonneg[temp] == 0)
				is_beta_nonneg[temp] = 1
				is_beta_link[rev_ind] = 1;
				# @show temp
				# @show rev_ind
			end
		else
			#is_gamma[rev_ind] = 1;
			gamma_ind = rev_ind;
			break;
		end
		#cur_prod = cur_prod+1;
	end

	#guess_revenue = sum( y .* revenues_ordered);
	#guess_revenue = dot(y, revenues_ordered);
	# println("Guess revenue = $guess_revenue")

	# Dual phase
	# println("inline: Start of dual phase")
	#gamma_ind =  find(is_gamma)
	if (gamma_ind == -1) #(isempty(gamma_ind))
		gamma_val = 0.0;
	else
		gamma_val = revenues_ordered[gamma_ind];
	end

	beta_inds = find(is_beta_nonneg .> 0)
	beta_link_inds = find(is_beta_link .> 0)
	beta_val = zeros(Float64, n+1);

	if (!isempty(beta_inds))
		beta_val[beta_inds[1]] = max(revenues_ordered[beta_link_inds[1]] - gamma_val, 0)
		for j=2:length(beta_inds)
			beta_val[beta_inds[j]] = max(revenues_ordered[beta_link_inds[j]] - gamma_val - sum(beta_val[1:(beta_link_inds[j]-1)]) , 0);
		end
		# println("here")
	end

	alpha_inds = find(is_alpha_nonneg)
	alpha_val = zeros(Float64,n+1);
	for i=1:length(alpha_inds)
		# alpha_val[alpha_inds[i]] = maximum( [ revenues_ordered[alpha_inds[i]] - gamma_val - sum( beta_val[1:(alpha_inds[i]-1)] ), 0])
		alpha_val[alpha_inds[i]] = max( revenues_ordered[alpha_inds[i]] - gamma_val - sum( beta_val[1:(alpha_inds[i]-1)] ), 0)
	end

	# guess_dual_rev = gamma_val + sum(alpha_val .* x_ordered[:]) + sum(beta_val .* (1 - x_ordered)[:])

	alpha_val = alpha_val[p2r[:]];
	beta_val = beta_val[p2r[:]];

	#guess_dual_rev_unsorted = gamma_val + sum(alpha_val .* temp_x[:]) + sum(beta_val .* (1 - temp_x)[:])
	subobj_val = guess_revenue;

	return subobj_val, alpha_val, beta_val, gamma_val  #, y
end

max_stars_repo_name: vvmisic/EFCPLO

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: vvmisic/EFCPLO

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: vvmisic/EFCPLO

max_forks_repo_licenses: MIT Unlicense


text:
# syntax: proto2
using Compat
using ProtoBuf
import ProtoBuf.meta

struct __enum_RpcKindProto <: ProtoEnum
    RPC_BUILTIN::Int32
    RPC_WRITABLE::Int32
    RPC_PROTOCOL_BUFFER::Int32
    __enum_RpcKindProto() = new(0,1,2)
end #struct __enum_RpcKindProto
const RpcKindProto = __enum_RpcKindProto()

mutable struct RPCTraceInfoProto <: ProtoType
    traceId::Int64
    parentId::Int64
    RPCTraceInfoProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RPCTraceInfoProto

struct __enum_RpcRequestHeaderProto_OperationProto <: ProtoEnum
    RPC_FINAL_PACKET::Int32
    RPC_CONTINUATION_PACKET::Int32
    RPC_CLOSE_CONNECTION::Int32
    __enum_RpcRequestHeaderProto_OperationProto() = new(0,1,2)
end #struct __enum_RpcRequestHeaderProto_OperationProto
const RpcRequestHeaderProto_OperationProto = __enum_RpcRequestHeaderProto_OperationProto()

mutable struct RpcRequestHeaderProto <: ProtoType
    rpcKind::Int32
    rpcOp::Int32
    callId::Int32
    clientId::Array{UInt8,1}
    retryCount::Int32
    traceInfo::RPCTraceInfoProto
    RpcRequestHeaderProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RpcRequestHeaderProto
const __req_RpcRequestHeaderProto = Symbol[:callId,:clientId]
const __val_RpcRequestHeaderProto = Dict(:retryCount => -1)
const __wtype_RpcRequestHeaderProto = Dict(:callId => :sint32, :retryCount => :sint32)
meta(t::Type{RpcRequestHeaderProto}) = meta(t, __req_RpcRequestHeaderProto, ProtoBuf.DEF_FNUM, __val_RpcRequestHeaderProto, true, ProtoBuf.DEF_PACK, __wtype_RpcRequestHeaderProto, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

struct __enum_RpcResponseHeaderProto_RpcStatusProto <: ProtoEnum
    SUCCESS::Int32
    ERROR::Int32
    FATAL::Int32
    __enum_RpcResponseHeaderProto_RpcStatusProto() = new(0,1,2)
end #struct __enum_RpcResponseHeaderProto_RpcStatusProto
const RpcResponseHeaderProto_RpcStatusProto = __enum_RpcResponseHeaderProto_RpcStatusProto()

struct __enum_RpcResponseHeaderProto_RpcErrorCodeProto <: ProtoEnum
    ERROR_APPLICATION::Int32
    ERROR_NO_SUCH_METHOD::Int32
    ERROR_NO_SUCH_PROTOCOL::Int32
    ERROR_RPC_SERVER::Int32
    ERROR_SERIALIZING_RESPONSE::Int32
    ERROR_RPC_VERSION_MISMATCH::Int32
    FATAL_UNKNOWN::Int32
    FATAL_UNSUPPORTED_SERIALIZATION::Int32
    FATAL_INVALID_RPC_HEADER::Int32
    FATAL_DESERIALIZING_REQUEST::Int32
    FATAL_VERSION_MISMATCH::Int32
    FATAL_UNAUTHORIZED::Int32
    __enum_RpcResponseHeaderProto_RpcErrorCodeProto() = new(1,2,3,4,5,6,10,11,12,13,14,15)
end #struct __enum_RpcResponseHeaderProto_RpcErrorCodeProto
const RpcResponseHeaderProto_RpcErrorCodeProto = __enum_RpcResponseHeaderProto_RpcErrorCodeProto()

mutable struct RpcResponseHeaderProto <: ProtoType
    callId::UInt32
    status::Int32
    serverIpcVersionNum::UInt32
    exceptionClassName::AbstractString
    errorMsg::AbstractString
    errorDetail::Int32
    clientId::Array{UInt8,1}
    retryCount::Int32
    RpcResponseHeaderProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RpcResponseHeaderProto
const __req_RpcResponseHeaderProto = Symbol[:callId,:status]
const __val_RpcResponseHeaderProto = Dict(:retryCount => -1)
const __wtype_RpcResponseHeaderProto = Dict(:retryCount => :sint32)
meta(t::Type{RpcResponseHeaderProto}) = meta(t, __req_RpcResponseHeaderProto, ProtoBuf.DEF_FNUM, __val_RpcResponseHeaderProto, true, ProtoBuf.DEF_PACK, __wtype_RpcResponseHeaderProto, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

struct __enum_RpcSaslProto_SaslState <: ProtoEnum
    SUCCESS::Int32
    NEGOTIATE::Int32
    INITIATE::Int32
    CHALLENGE::Int32
    RESPONSE::Int32
    WRAP::Int32
    __enum_RpcSaslProto_SaslState() = new(0,1,2,3,4,5)
end #struct __enum_RpcSaslProto_SaslState
const RpcSaslProto_SaslState = __enum_RpcSaslProto_SaslState()

mutable struct RpcSaslProto_SaslAuth <: ProtoType
    method::AbstractString
    mechanism::AbstractString
    protocol::AbstractString
    serverId::AbstractString
    challenge::Array{UInt8,1}
    RpcSaslProto_SaslAuth(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RpcSaslProto_SaslAuth
const __req_RpcSaslProto_SaslAuth = Symbol[:method,:mechanism]
meta(t::Type{RpcSaslProto_SaslAuth}) = meta(t, __req_RpcSaslProto_SaslAuth, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct RpcSaslProto <: ProtoType
    version::UInt32
    state::Int32
    token::Array{UInt8,1}
    auths::Base.Vector{RpcSaslProto_SaslAuth}
    RpcSaslProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RpcSaslProto
const __req_RpcSaslProto = Symbol[:state]
meta(t::Type{RpcSaslProto}) = meta(t, __req_RpcSaslProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

export RpcKindProto, RPCTraceInfoProto, RpcRequestHeaderProto_OperationProto, RpcRequestHeaderProto, RpcResponseHeaderProto_RpcStatusProto, RpcResponseHeaderProto_RpcErrorCodeProto, RpcResponseHeaderProto, RpcSaslProto_SaslState, RpcSaslProto_SaslAuth, RpcSaslProto


text:
module KernelRecursiveLeastSquares

using KernelLibrary

export krls

include("krls.jl")

end # module

max_stars_repo_name: the-moliver/KernelRecursiveLeastSquares

max_issues_repo_name: the-moliver/KernelRecursiveLeastSquares

max_forks_repo_name: the-moliver/KernelRecursiveLeastSquares


text:
using BinaryBuilder

name = "MKL"
version = v"2020.2.254"

sources = [
    ArchiveSource("https://anaconda.org/intel/mkl/2020.2/download/linux-64/mkl-2020.2-intel_254.tar.bz2",
                  "930d67bb4298c6da2dabb8ea068170e69e9304fc0ce1b4fc094af01851102a35"; unpack_target = "mkl-x86_64-linux-gnu"),
    ArchiveSource("https://anaconda.org/intel/mkl/2020.2/download/osx-64/mkl-2020.2-intel_258.tar.bz2",
                  "628c54329ab3c088b4c55d047ccb58feeb2ade7d6e50bf982aa51ac088cccd45"; unpack_target = "mkl-x86_64-apple-darwin14"),
    ArchiveSource("https://anaconda.org/intel/mkl/2020.2/download/win-32/mkl-2020.2-intel_254.tar.bz2",
                  "19dfae9402e764e507ee143b2c22cbf091fd6b778b03620504addc2b8135f987"; unpack_target = "mkl-i686-w64-mingw32"),
    ArchiveSource("https://anaconda.org/intel/mkl/2020.2/download/win-64/mkl-2020.2-intel_254.tar.bz2",
                  "99efbdd8014668f1683aec61ea190aa7182a90e14969439a7dfb7b3dfef55693"; unpack_target = "mkl-x86_64-w64-mingw32"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/mkl-${target}
if [[ ${target} == *mingw* ]]; then
    cp -r Library/bin/* ${libdir}
    install_license info/*.txt
else
    cp -r lib/* ${libdir}
    install_license info/licenses/*.txt
fi
"""

platforms = [
    Linux(:x86_64, libc=:glibc),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64),
]

# The products that we will ensure are always built
products = [
    LibraryProduct(["libmkl_core", "mkl_core"], :libmkl_core),
    LibraryProduct(["libmkl_rt", "mkl_rt"], :libmkl_rt),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("IntelOpenMP_jll"),
]

non_reg_ARGS = filter(arg -> arg != "--register", ARGS)
include("../../fancy_toys.jl")
no_autofix_platforms = [Windows(:i686), Windows(:x86_64), MacOS(:x86_64)]
autofix_platforms = [Linux(:x86_64)]
if any(should_build_platform.(triplet.(no_autofix_platforms)))
    # Need to disable autofix: updating linkage of libmkl_intel_thread.dylib on
    # macOS causes runtime issues:
    # https://github.com/JuliaPackaging/Yggdrasil/issues/915.
    build_tarballs(non_reg_ARGS, name, version, sources, script, no_autofix_platforms, products, dependencies; lazy_artifacts = true, autofix = false)
end
if any(should_build_platform.(triplet.(autofix_platforms)))
    # ... but we need to run autofix on Linux, because here libmkl_rt doesn't
    # have a soname, so we can't ccall it without specifying the path:
    # https://github.com/JuliaSparse/Pardiso.jl/issues/69
    build_tarballs(ARGS, name, version, sources, script, autofix_platforms, products, dependencies; lazy_artifacts = true)
end

max_stars_repo_name: jpsamaroo/Yggdrasil

max_issues_repo_name: jpsamaroo/Yggdrasil

max_forks_repo_name: jpsamaroo/Yggdrasil


text:
using POMDPs
using POMDPModelTools
using RLInterface

# Define a test environment
# it has 2 states, it ends up after taking 5 action
# Visiting the second state multiply the next reward by 10
# Optimal value 2.1
# Optimal Policy [2,1,2,1,3]

mutable struct TestMDP <: MDP{Tuple{Vector{Int32}, Int32}, Int64}
    shape
    stack::Int32
    o_stack::Int32
    max_time::Int32
    bad_state::Array{Int32}
    normal_state::Array{Int32}
    good_state::Array{Int32}
    _observation_space::Array{Array{Float32}}
    _rewards::Array{Float32}
    discount_factor::Float32
end

RLInterface.obsvector_type(mdp::TestMDP) = Array{Float32, length(mdp.shape) + 1}

function TestMDP(shape=(6,), stack=4, max_time=6, discount_factor=0.99)
    bad_state =  convert.(Int32, rand(1:50, shape...))
    normal_state =  convert.(Int32, rand(100:150, shape...))
    good_state =  convert.(Int32, rand(150:200, shape...))
    _observation_space = [bad_state, normal_state, good_state]
    _rewards = [-0.1f0, 0.0f0, 0.1f0]
    return TestMDP(shape, 4, stack, max_time, bad_state, normal_state, good_state, _observation_space, _rewards, discount_factor)
end

POMDPs.discount(mdp::TestMDP) = mdp.discount_factor

function POMDPs.observations(mdp::TestMDP)
    return mdp._observation_space
end

function POMDPs.actions(mdp::TestMDP)
    return 1:4
end

POMDPs.actionindex(mdp::TestMDP, a::T) where T<:Integer = a

function POMDPs.initialstate(mdp::TestMDP)
    ImplicitDistribution() do rng
        init_t = Int32(1)
        init_s = fill(Int32(1), mdp.stack)
        return (init_s, init_t)
    end
end

function POMDPs.convert_s(t::Type{V}, s::Tuple{Vector{T}, T}, mdp::TestMDP) where {T,V<:AbstractArray}
    obs = zeros(mdp.shape..., mdp.o_stack)
    for i=1:mdp.o_stack
        obs[Base.setindex(axes(obs), i, ndims(obs))...] = observations(mdp)[s[1][end-i+1]]
    end
    return convert(t, obs./255.0f0)
end

function was_in_second(s::Tuple{Vector{T}, T}) where T<:Integer
    s[1][end] == 2
end

function POMDPs.gen(mdp::TestMDP, s::Tuple{Vector{T}, T}, a::N, rng::AbstractRNG) where {T <: Integer, N <: Integer}
    t_new = s[2] + convert(Int32, 1) # increment time
    s_new = circshift(s[1], -1)
    if a < 4
        s_new[end] = a
    else
        s_new[end] = s_new[end-1]
    end
    return (sp=(s_new, t_new),)
end

function POMDPs.reward(mdp::TestMDP, s::Tuple{Vector{T}, T}, a::N, sp::Tuple{Vector{T}, T}) where {T <: Integer, N <: Integer}
    r = mdp._rewards[sp[1][end]]
    if was_in_second(s)
        r *= -10
    end
    return r
end

function POMDPs.isterminal(mdp::TestMDP, s)
    return s[2] >= mdp.max_time
end


text:
export Sprite, draw!, parent, position, position, position!, visible!, visible

"""
A visualization of a game character, static object or tile
in the game world. While you can use this directly, ideally
you use it with higher level objects like `Actor` to manage
its location in the world.
"""
mutable struct Sprite
    parent::Union{Sprite, Nothing}
    pos::Point{Float64}
    frames::Vector{Image}
    frame::Int
    visible::Bool
end

"""
    Sprite(frames)
Create a sprite with given frames and an optional parent sprite,
for which sprite position will be relative to.
"""
function Sprite(frames::AbstractVector{Image}, pos = Point(0.0, 0.0))
    Sprite(nothing, pos, frames, 1, true)
end

"""
    draw!(canvas::Image, sprite)
Draws sprite on canvas using position `world_position(sprite)` and
current frame.
"""
function draw!(canvas::Image, sprite::Sprite)
    if sprite.visible
        frame = sprite.frame
        draw!(canvas, sprite.frames[frame], sprite.pos)
    end
end

"""
    visible!(sprite::Sprite, visible = true)
Set a sprite visible or hidden. If it is hidden it will not be drawn when
the `draw!` call is made.
"""
function visible!(sprite::Sprite, visible = true)
    sprite.visible = visible
end

"""
    visible(sprite::Sprite)
Check if sprite is visible. If it is true, then a `draw!(canvas, sprite)` call
will make it appear on the canvas.
"""
function visible(sprite::Sprite)
    sprite.visible
end